Skip to content

Commit

Permalink
Merge pull request #50 from gridap/adding_vector_laplacian_non_confor…
Browse files Browse the repository at this point in the history 
…ming_meshes

Adding vector laplacian on  non conforming meshes
  • Loading branch information
@amartinhuertas
amartinhuertas authored Oct 28, 2023
2 parents 02c7cb0 + 57a7720 commit a159a22
Show file tree
Hide file tree
Showing 6 changed files with 285 additions and 128 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/ci.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ jobs:
fail-fast: false
matrix:
version:
- '1.9'
- '1.8'
os:
- ubuntu-20.04
arch:
Expand Down
221 changes: 170 additions & 51 deletions src/FESpaces.jl
View file
Original file line number Diff line number Diff line change
@@ -1,28 +1,102 @@

function _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe::Tuple{<:Lagrangian,Any,Any})
function _h_refined_reffe(reffe::Tuple{<:Lagrangian,Any,Any})
(reffe[1], (reffe[2][1], 2 * reffe[2][2]), reffe[3])
end
cell_polytope = Dc == 2 ? QUAD : HEX
basis, reffe_args, reffe_kwargs = reffe
cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...)
h_refined_reffe = _h_refined_reffe(reffe)
basis, reffe_args, reffe_kwargs = h_refined_reffe
cell_reffe_h_refined = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...)
dof_basis_h_refined = Gridap.CellData.get_dof_basis(cell_reffe_h_refined)
coarse_shape_funs = Gridap.ReferenceFEs.get_shapefuns(cell_reffe)
ref_constraints = evaluate(dof_basis_h_refined, coarse_shape_funs)

# TO-DO: How can modelH be created such that it is tailored to cell_polytope?
modelH= _generate_unit_hypercube_model(Dc)
modelh=refine(modelH,2)

VH=TestFESpace(modelH,cell_reffe)
Vh=TestFESpace(modelh,cell_reffe)

uH=get_fe_basis(VH)
uHh=change_domain(uH,get_triangulation(modelh),ReferenceDomain())
σRTh=Gridap.FESpaces.get_fe_dof_basis(Vh)
ref_constraints_contribs=σRTh(uHh)

ref_constraints = Matrix{Float64}(undef,num_free_dofs(Vh),num_free_dofs(VH))
cell_dof_ids = get_cell_dof_ids(Vh)
cache_cell_dof_ids = array_cache(cell_dof_ids)
cache_ref_constraints_contribs = array_cache(ref_constraints_contribs)
for cell=1:length(cell_dof_ids)
current_cell_dof_ids=getindex!(cache_cell_dof_ids,cell_dof_ids,cell)
current_ref_constraints_contribs=getindex!(cache_ref_constraints_contribs,ref_constraints_contribs,cell)
ref_constraints[current_cell_dof_ids,:]=current_ref_constraints_contribs
end
ref_constraints
end

function _fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof,
num_faces,
coarse_faces_to_child_ids,
face_dofs,
cells_dof_ids,
first_face)
for coarse_face_id=1:num_faces
for (subface,child_id) in enumerate(coarse_faces_to_child_ids[coarse_face_id,:])
@debug "coarse_face_id: $(coarse_face_id), subface: $(subface), child_id: $(child_id)"
cell_dof_ids=cells_dof_ids[child_id]
for (i,dof) in enumerate(cell_dof_ids[face_dofs[first_face+coarse_face_id]])
@debug "i: $(i), dof: $(dof)"
face_subface_ldof_to_cell_ldof[coarse_face_id][subface][i]=dof
end
end
end
end


function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:Lagrangian,Any,Any})
cell_polytope = (Dc == 2) ? QUAD : HEX
rr = Gridap.Adaptivity.RedRefinementRule(cell_polytope)
order=reffe[2][2]
Gridap.Adaptivity.get_face_subface_ldof_to_cell_ldof(rr, Tuple(fill(order,Dc)), Df)
coarse_faces_to_child_ids = (Dc == 2) ? _coarse_faces_to_child_ids_2D :
_coarse_faces_to_child_ids_3D

basis, reffe_args, reffe_kwargs = reffe
cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...)

# TO-DO: How can modelH be created such that it is tailored to cell_polytope?
modelH= _generate_unit_hypercube_model(Dc)
modelh=refine(modelH,2)
Vh=TestFESpace(modelh,reffe)
cells_dof_ids=get_cell_dof_ids(Vh)

first_face = get_offset(get_polytope(cell_reffe),Df)
face_dofs=get_face_dofs(cell_reffe)
num_dofs_x_face = length(face_dofs[first_face+1])

if (Df==Dc-1) # Facets
num_faces = 2*Dc
num_subfaces = 2^(Dc-1)
face_subface_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_faces,num_subfaces,num_dofs_x_face)
_fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof,
num_faces,
coarse_faces_to_child_ids,
face_dofs,
cells_dof_ids,
first_face)
face_subface_ldof_to_cell_ldof
else
@assert Df==1 && Dc==3 # Edges in 3D
num_edges=12
num_subedges=2
edge_subedge_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_edges,num_subedges,num_dofs_x_face)
_fill_face_subface_ldof_to_cell_ldof!(edge_subedge_ldof_to_cell_ldof,
num_edges,
_coarse_edges_to_child_ids_3D,
face_dofs,
cells_dof_ids,
first_face)
edge_subedge_ldof_to_cell_ldof
end
end



const _coarse_faces_to_child_ids_2D=[1 2; 3 4; 1 3; 2 4]
const _coarse_faces_to_child_ids_3D=[1 2 3 4; 5 6 7 8; 1 2 5 6; 3 4 7 8; 1 3 5 7; 2 4 6 8; ]
const _coarse_edges_to_child_ids_3D=[1 2; 3 4; 5 6; 7 8; 1 3; 2 4; 5 7; 6 8; 1 5; 2 6; 3 7; 4 8; ]


function _generate_unit_hypercube_model(Dc)
@assert Dc==2 || Dc==3
Expand Down Expand Up @@ -58,27 +132,19 @@ function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:RaviartTh
modelH= _generate_unit_hypercube_model(Dc)
modelh=refine(modelH,2)
RTh=TestFESpace(modelh,reffe)

num_faces = 2*Dc
num_subfaces = 2^(Dc-1)
first_face = get_offset(get_polytope(cell_reffe),Df)

face_own_dofs=get_face_own_dofs(cell_reffe)
num_dofs_x_face = length(face_own_dofs[first_face+1])

face_subface_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_faces,num_subfaces,num_dofs_x_face)

RTh_cell_dof_ids=get_cell_dof_ids(RTh)
for coarse_face_id=1:num_faces
for (subface,child_id) in enumerate(coarse_faces_to_child_ids[coarse_face_id,:])
@debug "coarse_face_id: $(coarse_face_id), subface: $(subface), child_id: $(child_id)"
cell_dof_ids=RTh_cell_dof_ids[child_id]
for (i,dof) in enumerate(cell_dof_ids[face_own_dofs[first_face+coarse_face_id]])
@debug "i: $(i), dof: $(dof)"
face_subface_ldof_to_cell_ldof[coarse_face_id][subface][i]=dof
end
end
end
_fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof,
num_faces,
coarse_faces_to_child_ids,
face_own_dofs,
RTh_cell_dof_ids,
first_face)
else
@assert Df==1 && Dc==3 # Edges in 3D
num_edges=12
Expand Down Expand Up @@ -230,16 +296,18 @@ function _restrict_face_dofs_to_face_dim(cell_reffe,Df)
first_face_faces = Gridap.ReferenceFEs.get_faces(polytope)[first_face]
face_dofs_to_face_dim = face_own_dofs[first_face_faces]
touched=Dict{Int,Int}()
dofs=Int64[]
current=1
for (i,face_dofs) in enumerate(face_dofs_to_face_dim)
for (j,dof) in enumerate(face_dofs)
if haskey(touched,dof)
face_dofs[j]=touched[dof]
else
touched[dof]=current
face_dofs[j]=touched[dof]
current=current+1
end
push!(dofs,dof)
end
end
sort!(dofs)
touched = Dict( dofs[i]=>i for i=1:length(dofs))
for (i,face_dofs) in enumerate(face_dofs_to_face_dim)
for (j,dof) in enumerate(face_dofs)
face_dofs[j]=touched[dof]
end
end
face_dofs_to_face_dim
Expand Down Expand Up @@ -508,10 +576,66 @@ function _compute_owner_edges_pindex_and_lids(
owner_edge_lid, _ = owner_edges_lids[owner_edge]
owner_edges_pindex[owner_edge_lid]=pindex
end
end
end
owner_edges_pindex, owner_edges_lids
end

using Gridap.ReferenceFEs

function get_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity})
p = get_polytope(reffe)
face_nodes = get_face_nodes(reffe)
dofs = get_dof_basis(reffe)
interior_nodes = dofs.nodes[face_nodes[end]]
Gridap.ReferenceFEs._compute_node_permutations(p,interior_nodes)
end

function get_nodes_permutations(reffe::ReferenceFE,d::Integer)
p = get_polytope(reffe)
range = get_dimrange(p,d)
get_face_nodes_permutations(reffe)[range]
end

function get_face_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity})
nodes_permutations = get_nodes_permutations(reffe)
reffaces = reffe.reffe.metadata
_reffaces = vcat(reffaces...)
face_nodes_permutations = map(get_nodes_permutations,_reffaces)
push!(face_nodes_permutations,nodes_permutations)
face_nodes_permutations
end

function get_face_dofs_permutations(reffe::LagrangianRefFE)
dofs = get_dof_basis(reffe)
face_nodes_permutations = get_face_nodes_permutations(reffe)
face_nodes = get_face_nodes(reffe)
face_dofs = get_face_dofs(reffe)
face_dofs_permutations = Gridap.ReferenceFEs._generate_face_own_dofs_permutations(
face_nodes_permutations, dofs.node_and_comp_to_dof, face_nodes, face_dofs)
face_dofs_permutations
end

function get_face_dofs_permutations(reffe::ReferenceFE,d::Integer)
p = get_polytope(reffe)
range = get_dimrange(p,d)
get_face_dofs_permutations(reffe)[range]
end

function get_face_dofs_permutations(
reffe::Gridap.ReferenceFEs.GenericRefFE{Gridap.ReferenceFEs.RaviartThomas, Dc},Df::Integer) where Dc
first_face = get_offset(get_polytope(reffe),Df)
order = length(get_face_dofs(reffe)[first_face])-1
nfaces=num_faces(reffe,Df)
if (Df==Dc-1)
facet_polytope = Dc == 2 ? SEGMENT : QUAD
nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope, [order for i = 1:Dc-1])
Fill(Gridap.ReferenceFEs._compute_node_permutations(facet_polytope, nodes),nfaces)
elseif (Dc == 3 && Df==1)
nodes, _ = Gridap.ReferenceFEs.compute_nodes(SEGMENT, [order for i = 1:Dc-2])
Fill(Gridap.ReferenceFEs._compute_node_permutations(SEGMENT, nodes),nfaces)
end
end

function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
spaces_wo_constraints,
reffe,
Expand Down Expand Up @@ -575,11 +699,10 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},

basis, reffe_args, reffe_kwargs = reffe
cell_reffe = ReferenceFE(Dc == 2 ? QUAD : HEX, basis, reffe_args...; reffe_kwargs...)
reffe_cell = cell_reffe

cell_dof_ids = get_cell_dof_ids(V)
face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(reffe_cell)
face_dofs = Gridap.ReferenceFEs.get_face_dofs(reffe_cell)
face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(cell_reffe)
face_dofs = Gridap.ReferenceFEs.get_face_dofs(cell_reffe)

hanging_faces_owner_face_dofs = Vector{Vector{Vector{Int}}}(undef, Dc)

Expand Down Expand Up @@ -610,7 +733,6 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
end

basis, reffe_args, reffe_kwargs = reffe
face_reffe = ReferenceFE(facet_polytope, basis, reffe_args...; reffe_kwargs...)
pindex_to_cfvertex_to_fvertex = Gridap.ReferenceFEs.get_vertex_permutations(facet_polytope)

if (Dc == 3)
Expand Down Expand Up @@ -642,22 +764,19 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
gridap_cell_faces[1],
gridap_cell_faces[2])
end


node_permutations = Vector{Vector{Vector{Int}}}(undef, Dc - 1)
nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope, [reffe_args[2] for i = 1:Dc-1])
node_permutations[Dc-1] = Gridap.ReferenceFEs._compute_node_permutations(facet_polytope, nodes)

face_dofs_permutations = Vector{Vector{Vector{Int}}}(undef, Dc-1)
face_dofs_permutations[Dc-1] =
first(get_face_dofs_permutations(cell_reffe, Dc-1))
if (Dc == 3)
nodes, _ = Gridap.ReferenceFEs.compute_nodes(edget_polytope, [reffe_args[2] for i = 1:Dc-2])
node_permutations[1] = Gridap.ReferenceFEs._compute_node_permutations(edget_polytope, nodes)
end
face_dofs_permutations[1] =
first(get_face_dofs_permutations(cell_reffe, 1))
end

subface_own_dofs = Vector{Vector{Vector{Int}}}(undef, Dc - 1)
subface_own_dofs[Dc-1] = _restrict_face_dofs_to_face_dim(cell_reffe,Dc-1)
#Gridap.ReferenceFEs.get_face_own_dofs(face_reffe)
if (Dc == 3)
subface_own_dofs[1] = _restrict_face_dofs_to_face_dim(cell_reffe,1)
#Gridap.ReferenceFEs.get_face_own_dofs(edge_reffe)
end
if (_num_face_own_dofs(cell_reffe,0)>0)
_generate_constraints!(0,
Expand All @@ -673,7 +792,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
face_own_dofs,
subface_own_dofs,
cell_dof_ids,
node_permutations,
face_dofs_permutations,
owner_faces_pindex,
owner_faces_lids,
ref_constraints,
Expand All @@ -697,7 +816,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
face_own_dofs,
subface_own_dofs,
cell_dof_ids,
node_permutations,
face_dofs_permutations,
owner_faces_pindex,
owner_faces_lids,
ref_constraints,
Expand All @@ -720,7 +839,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
face_own_dofs,
subface_own_dofs,
cell_dof_ids,
node_permutations,
face_dofs_permutations,
owner_faces_pindex,
owner_faces_lids,
ref_constraints,
Expand Down
19 changes: 13 additions & 6 deletions src/OctreeDistributedDiscreteModels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ mutable struct OctreeDistributedDiscreteModel{Dc,Dp,A,B,C,D,E,F} <: GridapDistri

if (isa(dmodel,GridapDistributed.DistributedDiscreteModel))
Gridap.Helpers.@check Dc == Gridap.Geometry.num_cell_dims(dmodel)
Gridap.Helpers.@check Dc == Gridap.Geometry.num_point_dims(dmodel)
Gridap.Helpers.@check Dp == Gridap.Geometry.num_point_dims(dmodel)
end

A = typeof(parts)
Expand Down Expand Up @@ -1230,11 +1230,9 @@ function Gridap.Adaptivity.refine(model::OctreeDistributedDiscreteModel{Dc,Dp};
end
end

function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp},
refinement_and_coarsening_flags::MPIArray{<:Vector};
parts=nothing) where {Dc,Dp}

Gridap.Helpers.@notimplementedif parts!=nothing
function _refine_coarsen_balance!(model::OctreeDistributedDiscreteModel{Dc,Dp},
refinement_and_coarsening_flags::MPIArray{<:Vector}) where {Dc,Dp}

# Variables which are updated accross calls to init_fn_callback_2d
current_quadrant_index_within_tree = Cint(0)
Expand Down Expand Up @@ -1430,6 +1428,16 @@ function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp},
pXest_coarsen!(Val{Dc}, ptr_new_pXest, coarsen_fn_callback_c)
pXest_balance!(Val{Dc}, ptr_new_pXest)
pXest_update_flags!(Val{Dc},model.ptr_pXest,ptr_new_pXest)
ptr_new_pXest
end

function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp},
refinement_and_coarsening_flags::MPIArray{<:Vector};
parts=nothing) where {Dc,Dp}

Gridap.Helpers.@notimplementedif parts!=nothing

ptr_new_pXest = _refine_coarsen_balance!(model, refinement_and_coarsening_flags)

# Extract ghost and lnodes
ptr_pXest_ghost = setup_pXest_ghost(Val{Dc}, ptr_new_pXest)
Expand Down Expand Up @@ -2140,7 +2148,6 @@ function setup_non_conforming_distributed_discrete_model(::Type{Val{Dc}},
face_labeling=generate_face_labeling(parts,
cell_prange,
coarse_discrete_model,
grid,
topology,
ptr_pXest,
ptr_pXest_ghost)
Expand Down
Loading

2 comments on commit a159a22

@amartinhuertas
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Error while trying to register: Version 0.3.2 already exists

Please sign in to comment.