Skip to content

Commit

Permalink
Changes required for generalizing interfaces as per required by
Browse files Browse the repository at this point in the history 
simplexified octrees
  • Loading branch information
@amartinhuertas
amartinhuertas committed Nov 30, 2024
1 parent c85d288 commit 922022e
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 33 deletions.
59 changes: 34 additions & 25 deletions src/FESpaces.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -389,23 +389,17 @@ function _generate_hanging_faces_owner_face_dofs(num_hanging_faces,
Gridap.Arrays.Table(data_owner_face_dofs, ptrs)
end

function face_dim(::Type{Val{Dc}}, face_lid) where {Dc}
num_vertices = GridapP4est.num_cell_vertices(Val{Dc})
num_edges = GridapP4est.num_cell_edges(Val{Dc})
num_faces = GridapP4est.num_cell_faces(Val{Dc})
function face_dim(num_vertices, num_edges, num_faces, face_lid)
if (face_lid <= num_vertices)
return 0
elseif (face_lid <= num_vertices + num_edges)
return 1
elseif (face_lid <= num_vertices + num_edges + num_faces)
return Dc - 1
return (num_edges==0 ? 1 : 2)
end
end

function face_lid_within_dim(::Type{Val{Dc}}, face_lid) where {Dc}
num_vertices = GridapP4est.num_cell_vertices(Val{Dc})
num_edges = GridapP4est.num_cell_edges(Val{Dc})
num_faces = GridapP4est.num_cell_faces(Val{Dc})
function face_lid_within_dim(num_vertices, num_edges, num_faces, face_lid)
if (face_lid <= num_vertices)
return face_lid
elseif (face_lid <= num_vertices + num_edges)
Expand Down Expand Up @@ -450,6 +444,9 @@ end

function _generate_constraints!(Df,
Dc,
num_vertices,
num_edges,
num_faces,
cell_faces,
num_hanging_faces,
hanging_faces_to_cell,
Expand All @@ -472,10 +469,6 @@ function _generate_constraints!(Df,
@assert Dc == 2 || Dc == 3
@assert 0 ≤ Df < Dc

num_vertices = GridapP4est.num_cell_vertices(Val{Dc})
num_edges = GridapP4est.num_cell_edges(Val{Dc})
num_faces = GridapP4est.num_cell_faces(Val{Dc})

offset = 0
if (Df ≥ 1)
offset += num_vertices
Expand All @@ -497,8 +490,8 @@ function _generate_constraints!(Df,
lface = hanging_faces_to_lface[fid_hanging]
ocell, ocell_lface, subface = hanging_faces_glue[fid_hanging]
if (ocell!=-1)
ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface)
oface_dim = face_dim(Val{Dc}, ocell_lface)
ocell_lface_within_dim = face_lid_within_dim(num_vertices, num_edges, num_faces, ocell_lface)
oface_dim = face_dim(num_vertices, num_edges, num_faces, ocell_lface)

if (Df == 0) # Am I a vertex?
hanging_lvertex_within_first_subface = 2^oface_dim
Expand Down Expand Up @@ -544,7 +537,7 @@ function _generate_constraints!(Df,
oface = getindex!(cache_cell_faces[oface_dim+1], cell_faces[oface_dim+1], ocell)[ocell_lface_within_dim]
oface_lid, _ = owner_faces_lids[oface_dim][oface]
pindex = owner_faces_pindex[oface_dim][oface_lid]
@debug "Df=$(Df) Dc=$(Dc) fid_hanging=$(fid_hanging) cell=$(cell) oface=$(oface) oface_lid=$(oface_lid) pindex=$(pindex) ocell=$(ocell) ocell_lface=$(ocell_lface) ocell_lface_within_dim=$(ocell_lface_within_dim) subface=$(subface) oface_dim=$(oface_dim) cur_subface_own_dofs=$(cur_subface_own_dofs) face_own_dofs=$(face_own_dofs[offset+lface])"
@debug "Df=$(Df) Dc=$(Dc) fid_hanging=$(fid_hanging) cell=$(cell) oface=$(oface) oface_lid=$(oface_lid) pindex=$(pindex) ocell=$(ocell) ocell_lface=$(ocell_lface) ocell_lface_within_dim=$(ocell_lface_within_dim) subface=$(subface) oface_dim=$(oface_dim) cur_subface_own_dofs=$(cur_subface_own_dofs) offset=$(offset) lface=$(lface) face_own_dofs[offset+lface]=$(face_own_dofs[offset+lface])"
for ((ldof, dof), ldof_subface) in zip(enumerate(face_own_dofs[offset+lface]), cur_subface_own_dofs)
push!(sDOF_to_dof, current_cell_dof_ids[dof])
push!(sDOF_to_dofs, hanging_faces_owner_face_dofs[fid_hanging])
Expand Down Expand Up @@ -643,9 +636,9 @@ function get_face_dofs_permutations(
end
end

function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
function generate_constraints(dmodel::GridapDistributed.DistributedDiscreteModel{Dc},
spaces_wo_constraints,
reffe,
cell_reffe,
ref_constraints,
face_subface_ldof_to_cell_ldof) where {Dc}

Expand Down Expand Up @@ -697,9 +690,6 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
model,
V

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

cell_dof_ids = get_cell_dof_ids(V)
face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(cell_reffe)
face_dofs = Gridap.ReferenceFEs.get_face_dofs(cell_reffe)
Expand All @@ -726,8 +716,6 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
sDOF_to_dof = Int[]
sDOF_to_dofs = Vector{Int}[]
sDOF_to_coeffs = Vector{Float64}[]

basis, reffe_args, reffe_kwargs = reffe

face_dofs_permutations = Vector{Vector{Vector{Int}}}(undef, Dc-1)
face_dofs_permutations[Dc-1] =
Expand All @@ -742,9 +730,17 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
if (Dc == 3)
subface_own_dofs[1] = _restrict_face_dofs_to_face_dim(cell_reffe,1)
end

num_vertices = Gridap.ReferenceFEs.num_faces(cell_reffe,0)
num_edges = Dc==3 ? Gridap.ReferenceFEs.num_faces(cell_reffe,1) : 0
num_faces = Gridap.ReferenceFEs.Gridap.ReferenceFEs.num_faces(cell_reffe,Dc-1)

if (_num_face_own_dofs(cell_reffe,0)>0)
_generate_constraints!(0,
Dc,
num_vertices,
num_edges,
num_faces,
[gridap_cell_faces[i] for i = 1:Dc],
num_hanging_faces[1],
hanging_faces_to_cell[1],
Expand All @@ -769,6 +765,9 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
if (_num_face_own_dofs(cell_reffe,1)>0)
_generate_constraints!(1,
Dc,
num_vertices,
num_edges,
num_faces,
[gridap_cell_faces[i] for i = 1:Dc],
num_hanging_faces[2],
hanging_faces_to_cell[2],
Expand All @@ -792,6 +791,9 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
if (_num_face_own_dofs(cell_reffe,Dc-1)>0)
_generate_constraints!(Dc - 1,
Dc,
num_vertices,
num_edges,
num_faces,
[gridap_cell_faces[i] for i = 1:Dc],
num_hanging_faces[Dc],
hanging_faces_to_cell[Dc],
Expand Down Expand Up @@ -886,8 +888,15 @@ function _add_constraints(model::OctreeDistributedDiscreteModel{Dc},
face_subface_ldof_to_cell_ldof[1] =
_generate_face_subface_ldof_to_cell_ldof(model.pXest_refinement_rule_type, 1, Dc, reffe)
end
basis, reffe_args, reffe_kwargs = reffe
polytope = Dc==2 ? QUAD : HEX
cell_reffe = ReferenceFE(polytope, basis, reffe_args...; reffe_kwargs...)
sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs =
generate_constraints(model, spaces_wo_constraints, reffe, ref_constraints, face_subface_ldof_to_cell_ldof)
generate_constraints(model,
spaces_wo_constraints,
cell_reffe,
ref_constraints,
face_subface_ldof_to_cell_ldof)
spaces_w_constraints = map(spaces_wo_constraints,
sDOF_to_dof,
sDOF_to_dofs,
Expand Down Expand Up @@ -916,7 +925,7 @@ end
# Generates a new DistributedSingleFieldFESpace composed
# by local FE spaces with linear multipoint constraints added
function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc},
reffe::LagragianOrNedelec;
reffe;
kwargs...) where {Dc}
spaces_wo_constraints = map(local_views(model)) do m
FESpace(m, reffe; kwargs...)
Expand Down
25 changes: 17 additions & 8 deletions src/OctreeDistributedDiscreteModels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,15 +38,17 @@ struct NonConformingGlue{Dc,A,B,C,D,E,F,G}
end
end

function _compute_owner_faces_lids(Df,Dc,num_hanging_faces,hanging_faces_glue,cell_faces)
function _compute_owner_faces_lids(Df,Dc,
num_cell_vertices,num_cell_edges,num_cell_faces,
num_hanging_faces,hanging_faces_glue,cell_faces)
num_owner_faces = 0
owner_faces_lids = Dict{Int,Tuple{Int,Int,Int}}()
for fid_hanging = 1:num_hanging_faces
ocell, ocell_lface, _ = hanging_faces_glue[fid_hanging]
if (ocell!=-1)
ocell_dim = face_dim(Val{Dc}, ocell_lface)
ocell_dim = face_dim(num_cell_vertices,num_cell_edges,num_cell_faces,ocell_lface)
if (ocell_dim == Df)
ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface)
ocell_lface_within_dim = face_lid_within_dim(num_cell_vertices,num_cell_edges,num_cell_faces,ocell_lface)
owner_face = cell_faces[ocell][ocell_lface_within_dim]
if !(haskey(owner_faces_lids, owner_face))
num_owner_faces += 1
Expand Down Expand Up @@ -88,6 +90,9 @@ end
# compute permutation id
function _compute_owner_faces_pindex_and_lids(Dc,
pXest_refinement_rule,
num_cell_vertices,
num_cell_edges,
num_cell_faces,
num_hanging_faces,
hanging_faces_glue,
hanging_faces_to_cell,
Expand All @@ -98,6 +103,9 @@ function _compute_owner_faces_pindex_and_lids(Dc,
pindex_to_cfvertex_to_fvertex)

owner_faces_lids =_compute_owner_faces_lids(Dc-1,Dc,
num_cell_vertices,
num_cell_edges,
num_cell_faces,
num_hanging_faces,
hanging_faces_glue,
cell_faces)
Expand All @@ -114,11 +122,11 @@ function _compute_owner_faces_pindex_and_lids(Dc,
@debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: fid_hanging=$(fid_hanging) ocell=$(ocell) ocell_lface=$(ocell_lface) subface=$(subface)"

if (ocell!=-1)
oface_dim = face_dim(Val{Dc}, ocell_lface)
oface_dim = face_dim(num_cell_vertices, num_cell_edges, num_cell_faces, ocell_lface)
if (oface_dim == Dc-1)
cell = hanging_faces_to_cell[fid_hanging]
lface = hanging_faces_to_lface[fid_hanging]
ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface)
ocell_lface_within_dim = face_lid_within_dim(num_cell_vertices, num_cell_edges, num_cell_faces, ocell_lface)
owner_face = cell_faces[ocell][ocell_lface_within_dim]
owner_face_lid, _ = owner_faces_lids[owner_face]
for fcorner in _subface_to_face_corners(pXest_refinement_rule, subface)
Expand All @@ -136,7 +144,7 @@ function _compute_owner_faces_pindex_and_lids(Dc,
owner_faces_pindex = Vector{Int}(undef, num_owner_faces)
for owner_face in keys(owner_faces_lids)
(owner_face_lid, ocell, ocell_lface) = owner_faces_lids[owner_face]
ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface)
ocell_lface_within_dim = face_lid_within_dim(num_cell_vertices, num_cell_edges, num_cell_faces, ocell_lface)
# Compute permutation id by comparing
# 1. cell_vertices[ocell][ocell_lface]
# 2. owner_face_vertex_ids
Expand Down Expand Up @@ -172,6 +180,7 @@ end


function _compute_owner_edges_pindex_and_lids(
num_cell_vertices, num_cell_edges, num_cell_faces,
num_hanging_edges,
hanging_edges_glue,
hanging_edges_to_cell,
Expand All @@ -194,9 +203,9 @@ function _compute_owner_edges_pindex_and_lids(
# Find the owner hanging edge
for fid_hanging = 1:num_hanging_edges
ocell, ocell_ledge, subedge = hanging_edges_glue[fid_hanging]
ocell_dim = face_dim(Val{Dc}, ocell_ledge)
ocell_dim = face_dim(num_cell_vertices, num_cell_edges, num_cell_faces, ocell_ledge)
if (ocell!=-1 && ocell_dim==1)
ocell_ledge_within_dim = face_lid_within_dim(Val{Dc}, ocell_ledge)
ocell_ledge_within_dim = face_lid_within_dim(num_cell_vertices, num_cell_edges, num_cell_faces, ocell_ledge)
cell = hanging_edges_to_cell[fid_hanging]
ledge = hanging_edges_to_ledge[fid_hanging]
gvertex1 = cell_vertices[cell][ledge_to_cvertices[ledge][subedge]]
Expand Down
6 changes: 6 additions & 0 deletions src/PXestTypeMethods.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1899,6 +1899,9 @@ function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType,

owner_faces_pindex[Dc-1], owner_faces_lids[Dc-1] = _compute_owner_faces_pindex_and_lids(Dc,
pXest_refinement_rule,
n_cell_vertices,
n_cell_edges,
n_cell_faces,
num_hanging_faces[Dc],
hanging_faces_glue[Dc],
hanging_faces_to_cell[Dc],
Expand All @@ -1911,6 +1914,9 @@ function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType,
if (Dc == 3)
owner_faces_pindex[1], owner_faces_lids[1]=
_compute_owner_edges_pindex_and_lids(
n_cell_vertices,
n_cell_edges,
n_cell_faces,
num_hanging_faces[2],
hanging_faces_glue[2],
hanging_faces_to_cell[2],
Expand Down

0 comments on commit 922022e

Please sign in to comment.