diff --git a/Project.toml b/Project.toml index dbe28dd..1c39130 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] ArgParse = "1" FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 1" -Gridap = "0.18.2" +Gridap = "0.18.6" GridapDistributed = "0.4" MPI = "0.20" P4est_wrapper = "0.2.2" diff --git a/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl b/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl index 8f7d6a1..ec7d689 100644 --- a/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl +++ b/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl @@ -18,6 +18,7 @@ function AnisotropicallyAdapted3DDistributedDiscreteModel( ptr_pXest_lnodes=setup_pXest_lnodes_nonconforming(pXest_type, ptr_pXest, ptr_pXest_ghost) dmodel,non_conforming_glue = setup_non_conforming_distributed_discrete_model(pXest_type, + PXestHorizontalRefinementRuleType(), parts, coarse_model, ptr_pXest_connectivity, @@ -87,6 +88,7 @@ function vertically_adapt(model::OctreeDistributedDiscreteModel{3,3}, # Build fine-grid mesh fmodel,non_conforming_glue = setup_non_conforming_distributed_discrete_model(model.pXest_type, + model.pXest_refinement_rule_type, model.parts, model.coarse_model, model.ptr_pXest_connectivity, @@ -154,6 +156,7 @@ function vertically_uniformly_refine(model::OctreeDistributedDiscreteModel) # Build fine-grid mesh fmodel,non_conforming_glue = setup_non_conforming_distributed_discrete_model(model.pXest_type, + model.pXest_refinement_rule_type, model.parts, model.coarse_model, model.ptr_pXest_connectivity, @@ -201,6 +204,7 @@ end function setup_non_conforming_distributed_discrete_model(pXest_type::P6estType, + pXest_refinement_rule_type::PXestRefinementRuleType, parts, coarse_discrete_model, ptr_pXest_connectivity, @@ -214,7 +218,7 @@ function setup_non_conforming_distributed_discrete_model(pXest_type::P6estType, gridap_cell_faces, non_conforming_glue= - generate_cell_faces_and_non_conforming_glue(pXest_type,ptr_pXest_lnodes, cell_prange) + generate_cell_faces_and_non_conforming_glue(pXest_type,pXest_refinement_rule_type,ptr_pXest_lnodes, cell_prange) nlvertices = map(non_conforming_glue) do ncglue @@ -341,6 +345,7 @@ function generate_face_labeling(pXest_type::P6estType, ptr_p4est_lnodes = setup_pXest_lnodes_nonconforming(p4est_type, ptr_p4est, ptr_p4est_ghost) bottom_boundary_model,_=setup_non_conforming_distributed_discrete_model(p4est_type, + PXestUniformRefinementRuleType(), parts, coarse_discrete_model, ptr_p4est_connectivity, @@ -637,6 +642,7 @@ function horizontally_adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, # Build fine-grid mesh fmodel,non_conforming_glue = setup_non_conforming_distributed_discrete_model(model.pXest_type, + model.pXest_refinement_rule_type, model.parts, model.coarse_model, model.ptr_pXest_connectivity, diff --git a/src/FESpaces.jl b/src/FESpaces.jl index 6a06da1..93b36c2 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -352,37 +352,6 @@ function _build_constraint_coefficients_matrix_in_ref_space(::PXestUniformRefine end -# To-think: might this info go to the glue? -# If it is required in different scenarios, I would say it may make sense -function _generate_hanging_faces_to_cell_and_lface(num_regular_faces, - num_hanging_faces, - gridap_cell_faces) - # Locate for each hanging face the owner cell among the set of cells - # to which it belongs and local position within that cell - # By convention, the owner cell is the cell with minimum identifer among - # all in the set. This convention is used here, and in other parts of the code. - # Breaking this convention here may affect the consistency of other parts of the code. - hanging_faces_to_cell = Vector{Int}(undef, num_hanging_faces) - hanging_faces_to_lface = Vector{Int}(undef, num_hanging_faces) - hanging_faces_to_cell .= -1 - for cell = 1:length(gridap_cell_faces) - s = gridap_cell_faces.ptrs[cell] - e = gridap_cell_faces.ptrs[cell+1] - l = e - s - for j = 1:l - fid = gridap_cell_faces.data[s+j-1] - if fid > num_regular_faces - fid_hanging = fid - num_regular_faces - if (hanging_faces_to_cell[fid_hanging]==-1) - hanging_faces_to_cell[fid_hanging] = cell - hanging_faces_to_lface[fid_hanging] = j - end - end - end - end - hanging_faces_to_cell, hanging_faces_to_lface -end - function _generate_hanging_faces_owner_face_dofs(num_hanging_faces, face_dofs, hanging_faces_glue, @@ -603,178 +572,6 @@ function _generate_constraints!(Df, @debug "sDOF_to_coeffs [$(Df)]= $(sDOF_to_coeffs)" end -function _compute_owner_faces_lids(Df,Dc,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) - if (ocell_dim == Df) - ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) - owner_face = cell_faces[ocell][ocell_lface_within_dim] - if !(haskey(owner_faces_lids, owner_face)) - num_owner_faces += 1 - owner_faces_lids[owner_face] = (num_owner_faces, ocell, ocell_lface) - end - end - end - end - owner_faces_lids -end - -function _subface_to_face_corners(::PXestUniformRefinementRuleType, subface) - (subface,) -end - -function _subface_to_face_corners(::PXestHorizontalRefinementRuleType, subface) - @assert subface==1 || subface==2 - if (subface==1) - (1,2) - else - (3,4) - end -end - -function _subface_to_face_corners(::PXestVerticalRefinementRuleType, subface) - @assert subface==1 || subface==2 - if (subface==1) - (1,3) - else - (2,4) - end -end - - -# count how many different owner faces -# for each owner face -# track the global IDs of its face vertices from the perspective of the subfaces -# for each owner face -# compute permutation id -function _compute_owner_faces_pindex_and_lids(Dc, - pXest_refinement_rule, - num_hanging_faces, - hanging_faces_glue, - hanging_faces_to_cell, - hanging_faces_to_lface, - cell_vertices, - cell_faces, - lface_to_cvertices, - pindex_to_cfvertex_to_fvertex) - - owner_faces_lids =_compute_owner_faces_lids(Dc-1,Dc, - num_hanging_faces, - hanging_faces_glue, - cell_faces) - - @debug "owner_faces_lids [Df=$(Dc-1) Dc=$(Dc)]: $(owner_faces_lids)" - - num_owner_faces = length(keys(owner_faces_lids)) - num_face_vertices = length(first(lface_to_cvertices)) - owner_face_vertex_ids = Vector{Int}(undef, num_face_vertices * num_owner_faces) - owner_face_vertex_ids .= -1 - - for fid_hanging = 1:num_hanging_faces - ocell, ocell_lface, subface = hanging_faces_glue[fid_hanging] - @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) - 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) - 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) - cvertex = lface_to_cvertices[lface][fcorner] - vertex = cell_vertices[cell][cvertex] - @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: cell=$(cell) lface=$(lface) cvertex=$(cvertex) vertex=$(vertex) owner_face=$(owner_face) owner_face_lid=$(owner_face_lid)" - @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: owner_face_vertex_ids[$((owner_face_lid-1)*num_face_vertices+fcorner)] = $(vertex)" - owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fcorner] = vertex - end - end - end - end - @debug "owner_face_vertex_ids [Dc=$(Dc)]: $(owner_face_vertex_ids)" - - 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) - # Compute permutation id by comparing - # 1. cell_vertices[ocell][ocell_lface] - # 2. owner_face_vertex_ids - pindexfound = false - cfvertex_to_cvertex = lface_to_cvertices[ocell_lface_within_dim] - for (pindex, cfvertex_to_fvertex) in enumerate(pindex_to_cfvertex_to_fvertex) - found = true - for (cfvertex, fvertex) in enumerate(cfvertex_to_fvertex) - vertex1 = owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fvertex] - cvertex = cfvertex_to_cvertex[cfvertex] - vertex2 = cell_vertices[ocell][cvertex] - @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: cell_vertices[$(ocell)][$(cvertex)]=$(cell_vertices[ocell][cvertex]) owner_face_vertex_ids[$((owner_face_lid-1)*num_face_vertices+fvertex)]=$(owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fvertex])" - # -1 can only happen in the interface of two - # ghost cells at different refinement levels - if (vertex1 != vertex2) && (vertex1 != -1) - found = false - break - end - end - if found - owner_faces_pindex[owner_face_lid] = pindex - pindexfound = true - break - end - end - @assert pindexfound "Valid pindex not found" - end - @debug "owner_faces_pindex: $(owner_faces_pindex)" - - owner_faces_pindex, owner_faces_lids -end - - -function _compute_owner_edges_pindex_and_lids( - num_hanging_edges, - hanging_edges_glue, - hanging_edges_to_cell, - hanging_edges_to_ledge, - cell_vertices, - cell_edges) - Dc=3 - owner_edges_lids =_compute_owner_faces_lids(1, - Dc, - num_hanging_edges, - hanging_edges_glue, - cell_edges) - - num_owner_edges = length(keys(owner_edges_lids)) - owner_edges_pindex = Vector{Int}(undef, num_owner_edges) - - ledge_to_cvertices = Gridap.ReferenceFEs.get_faces(HEX, 1, 0) - - # Go over hanging edges - # 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) - if (ocell!=-1 && ocell_dim==1) - ocell_ledge_within_dim = face_lid_within_dim(Val{Dc}, 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]] - gvertex2 = cell_vertices[ocell][ledge_to_cvertices[ocell_ledge_within_dim][subedge]] - @debug "fid_hanging=$(fid_hanging) cell=$(cell) ledge=$(ledge) ocell=$(ocell) ocell_ledge=$(ocell_ledge) subedge=$(subedge) gvertex1=$(gvertex1) gvertex2=$(gvertex2)" - pindex = gvertex1==gvertex2 ? 1 : 2 - owner_edge=cell_edges[ocell][ocell_ledge_within_dim] - owner_edge_lid, _ = owner_edges_lids[owner_edge] - owner_edges_pindex[owner_edge_lid]=pindex - end - end - owner_edges_pindex, owner_edges_lids -end - using Gridap.ReferenceFEs function get_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity}) @@ -868,42 +665,37 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, hanging_faces_glue = map(non_conforming_glue) do ncglue Tuple(ncglue.hanging_faces_glue[d] for d = 1:Dc) end + hanging_faces_to_cell = map(non_conforming_glue) do ncglue + Tuple(ncglue.hanging_faces_to_cell[d] for d = 1:Dc) + end + hanging_faces_to_lface = map(non_conforming_glue) do ncglue + Tuple(ncglue.hanging_faces_to_lface[d] for d = 1:Dc) + end + owner_faces_pindex=map(non_conforming_glue) do ncglue + Tuple(ncglue.owner_faces_pindex[d] for d = 1:Dc-1) + end + owner_faces_lids=map(non_conforming_glue) do ncglue + Tuple(ncglue.owner_faces_lids[d] for d = 1:Dc-1) + end sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = map(gridap_cell_faces, num_regular_faces, num_hanging_faces, hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids, dmodel.dmodel.models, spaces_wo_constraints) do gridap_cell_faces, - num_regular_faces, - num_hanging_faces, - hanging_faces_glue, - model, - V - - hanging_faces_to_cell = Vector{Vector{Int}}(undef, Dc) - hanging_faces_to_lface = Vector{Vector{Int}}(undef, Dc) - - # Locate for each hanging vertex a cell to which it belongs - # and local position within that cell - hanging_faces_to_cell[1], - hanging_faces_to_lface[1] = _generate_hanging_faces_to_cell_and_lface(num_regular_faces[1], - num_hanging_faces[1], - gridap_cell_faces[1]) - - if (Dc == 3) - hanging_faces_to_cell[2], - hanging_faces_to_lface[2] = _generate_hanging_faces_to_cell_and_lface(num_regular_faces[2], - num_hanging_faces[2], - gridap_cell_faces[2]) - end - - # Locate for each hanging facet a cell to which it belongs - # and local position within that cell - hanging_faces_to_cell[Dc], - hanging_faces_to_lface[Dc] = - _generate_hanging_faces_to_cell_and_lface(num_regular_faces[Dc], - num_hanging_faces[Dc], - gridap_cell_faces[Dc]) + num_regular_faces, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids, + model, + V basis, reffe_args, reffe_kwargs = reffe cell_reffe = ReferenceFE(Dc == 2 ? QUAD : HEX, basis, reffe_args...; reffe_kwargs...) @@ -935,44 +727,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, sDOF_to_dofs = Vector{Int}[] sDOF_to_coeffs = Vector{Float64}[] - facet_polytope = Dc == 2 ? SEGMENT : QUAD - if (Dc == 3) - edget_polytope = SEGMENT - end - basis, reffe_args, reffe_kwargs = reffe - pindex_to_cfvertex_to_fvertex = Gridap.ReferenceFEs.get_vertex_permutations(facet_polytope) - - if (Dc == 3) - edge_reffe = ReferenceFE(edget_polytope, basis, reffe_args...; reffe_kwargs...) - pindex_to_cevertex_to_evertex = Gridap.ReferenceFEs.get_vertex_permutations(edget_polytope) - end - - owner_faces_pindex = Vector{Vector{Int}}(undef, Dc - 1) - owner_faces_lids = Vector{Dict{Int,Tuple{Int,Int,Int}}}(undef, Dc - 1) - - lface_to_cvertices = Gridap.ReferenceFEs.get_faces(Dc == 2 ? QUAD : HEX, Dc - 1, 0) - owner_faces_pindex[Dc-1], owner_faces_lids[Dc-1] = _compute_owner_faces_pindex_and_lids(Dc, - dmodel.pXest_refinement_rule_type, - num_hanging_faces[Dc], - hanging_faces_glue[Dc], - hanging_faces_to_cell[Dc], - hanging_faces_to_lface[Dc], - gridap_cell_faces[1], - gridap_cell_faces[Dc], - lface_to_cvertices, - pindex_to_cfvertex_to_fvertex) - - if (Dc == 3) - owner_faces_pindex[1], owner_faces_lids[1]= - _compute_owner_edges_pindex_and_lids( - num_hanging_faces[2], - hanging_faces_glue[2], - hanging_faces_to_cell[2], - hanging_faces_to_lface[2], - gridap_cell_faces[1], - gridap_cell_faces[2]) - end face_dofs_permutations = Vector{Vector{Vector{Int}}}(undef, Dc-1) face_dofs_permutations[Dc-1] = diff --git a/src/Geometry.jl b/src/Geometry.jl new file mode 100644 index 0000000..97af0bf --- /dev/null +++ b/src/Geometry.jl @@ -0,0 +1,293 @@ +struct FaceToCellGlueNonConforming{A,B,C,D} <: GridapType + face_is_owner :: Vector{Bool} + face_to_subface :: Vector{Int8} + face_to_cell_glue:: Gridap.Geometry.FaceToCellGlue{A,B,C,D} +end + +function _adjust_cell_to_lface_to_pindex!(f2cg_nc,Df,ncglue,cell_faces) + f2cg=f2cg_nc.face_to_cell_glue + face_to_cell=f2cg.face_to_cell + face_to_lface=f2cg.face_to_lface + cell_to_lface_to_pindex=f2cg.cell_to_lface_to_pindex + for i=1:length(f2cg_nc.face_is_owner) + if (f2cg_nc.face_is_owner[i]) + cell=face_to_cell[i] + lface=face_to_lface[i] + gface=cell_faces[cell][lface] + (oface,_,_)=ncglue.owner_faces_lids[Df][gface] + pindex=ncglue.owner_faces_pindex[Df][oface] + s=f2cg.cell_to_lface_to_pindex.ptrs[cell] + f2cg.cell_to_lface_to_pindex.data[s+lface-1]=pindex + end + end +end + +function FaceToCellGlueNonConforming(topo, + cell_grid, + face_grid, + face_to_bgface, + bgface_to_lcell, + face_is_owner, + face_to_subface, + ncglue) + + face_to_cell_glue=Gridap.Geometry.FaceToCellGlue(topo, + cell_grid, + face_grid, + face_to_bgface, + bgface_to_lcell) + + f2cg_nc=FaceToCellGlueNonConforming(face_is_owner, + face_to_subface, + face_to_cell_glue) + + Dc=num_cell_dims(cell_grid) + Df=num_cell_dims(face_grid) + cell_faces=Gridap.Geometry.get_faces(topo,Dc,Df) + _adjust_cell_to_lface_to_pindex!(f2cg_nc,Df,ncglue,cell_faces) + f2cg_nc +end + + +function _generate_owner_face_to_subfaces(model,ncglue::GridapP4est.NonConformingGlue{D}) where D + num_children=GridapP4est.get_num_children(Val{D-1}) + + # Generate owner_face_to_subfaces_ptrs + num_owner_faces=length(keys(ncglue.owner_faces_lids[D-1])) + owner_face_to_subfaces_ptrs = Vector{Int}(undef, num_owner_faces+1) + owner_face_to_subfaces_ptrs .= num_children + length_to_ptrs!(owner_face_to_subfaces_ptrs) + + topology=Gridap.Geometry.get_grid_topology(model) + cell_faces=Gridap.Geometry.get_faces(topology,D,D-1) + + num_hanging_faces = ncglue.num_hanging_faces[D] + num_regular_faces = ncglue.num_regular_faces[D] + + # Generate owner_face_to_subfaces_data + owner_face_to_subfaces_data=Vector{Int}(undef,owner_face_to_subfaces_ptrs[end]-1) + owner_face_to_subfaces_data.=-1 + for i=1:num_hanging_faces + (ocell,ocell_lface,subface)=ncglue.hanging_faces_glue[D][i] + if (ocell!=-1) + ocell_lface_within_dim =GridapP4est.face_lid_within_dim(Val{D}, ocell_lface) + owner_gface=cell_faces[ocell][ocell_lface_within_dim] + (lowner,_,_)=ncglue.owner_faces_lids[D-1][owner_gface] + spos=owner_face_to_subfaces_ptrs[lowner] + owner_face_to_subfaces_data[spos+subface-1]=num_regular_faces+i + end + end + Gridap.Arrays.Table(owner_face_to_subfaces_data,owner_face_to_subfaces_ptrs) +end + +function Gridap.Geometry.SkeletonTriangulation(model::DiscreteModel{D}, + ncglue::GridapP4est.NonConformingGlue{D}) where {D} + + num_regular_faces=ncglue.num_regular_faces[D] + owner_faces_lids=ncglue.owner_faces_lids + + topo=Gridap.Geometry.get_grid_topology(model) + is_boundary_face=Gridap.Geometry.get_isboundary_face(topo,D-1) + + # Include all regular faces wich are + # either owner or not in the boundary + face_to_regular_bgface=Int[] + for (gface,ibf) in enumerate(is_boundary_face) + if gface <= num_regular_faces + if (!ibf) + push!(face_to_regular_bgface,gface) + else + if (haskey(owner_faces_lids[D-1],gface)) + push!(face_to_regular_bgface,gface) + end + end + end + end + + # IMPORTANT NOTE: Plus side contains the hanging side of all cell interfaces, while + # Minus side the "owner" side of all cell interfaces. This is a MUST. + # For facet integration purposes, the Jacobian of the geometrical mapping + # is extracted from the plus side. + + # Plus side + regular_bgface_to_lcell_minus=Vector{Int8}(undef,num_regular_faces) + regular_bgface_to_lcell_minus.=2 + for i=1:length(is_boundary_face) + if ( i<=num_regular_faces && + is_boundary_face[i] && + !(haskey(owner_faces_lids[D-1],i))) + regular_bgface_to_lcell_minus[i]=1 + end + end + plus=BoundaryTriangulation( + model, + face_to_regular_bgface, + regular_bgface_to_lcell_minus, + ncglue) + + + # Minus side + regular_bgface_to_lcell_plus=Fill(Int8(1),num_regular_faces) + minus=BoundaryTriangulation( + model, + face_to_regular_bgface, + regular_bgface_to_lcell_plus, + ncglue) + + + SkeletonTriangulation(plus,minus) +end + +function Gridap.Geometry.BoundaryTriangulation( + model::DiscreteModel{D}, + face_to_regular_bgface::AbstractVector{<:Integer}, + regular_bgface_to_lcell::AbstractVector{<:Integer}, + ncglue::GridapP4est.NonConformingGlue{D}) where D + + @assert length(regular_bgface_to_lcell)==ncglue.num_regular_faces[D] + + T=eltype(face_to_regular_bgface) + face_to_bgface=T[] + face_is_owner=Bool[] + face_to_subface=Int8[] + bgface_to_lcell=Vector{T}(undef, + ncglue.num_regular_faces[D]+ncglue.num_hanging_faces[D]) + + bgface_to_lcell.=1 + + num_children=GridapP4est.get_num_children(Val{D-1}) + owner_face_to_subfaces=_generate_owner_face_to_subfaces(model,ncglue) + + @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: owner_face_to_subfaces=$(owner_face_to_subfaces)" + + # Generate face_to_bgface AND bgface_to_lcell + for (i,gface) in enumerate(face_to_regular_bgface) + if (haskey(ncglue.owner_faces_lids[D-1],gface)) + # Find all subfaces of current owner face! + (lowner,_,_)=ncglue.owner_faces_lids[D-1][gface] + subfaces=owner_face_to_subfaces[lowner] + + # Regular face is owner of several other hanging faces + if (regular_bgface_to_lcell[gface]==1) + bgface_to_lcell[gface]=1 + for i=1:num_children + if (subfaces[i]!=-1) + push!(face_to_bgface, gface) + push!(face_is_owner,true) + push!(face_to_subface,i) + end + end + else + for i=1:num_children + if (subfaces[i]!=-1) + push!(face_to_bgface, subfaces[i]) + bgface_to_lcell[subfaces[i]]=1 + push!(face_is_owner,false) + push!(face_to_subface,-1) + end + end + end + else + # Regular face which is not owner of any other hanging face + push!(face_to_bgface, gface) + bgface_to_lcell[gface]=regular_bgface_to_lcell[gface] + push!(face_is_owner,false) + push!(face_to_subface,-1) + end + end + + topo = Gridap.Geometry.get_grid_topology(model) + cell_grid = get_grid(model) + bgface_grid = Gridap.Geometry.Grid(ReferenceFE{D-1},model) + bgface_grid = view(bgface_grid,face_to_bgface) + + glue=FaceToCellGlueNonConforming(topo, + cell_grid, + bgface_grid, + face_to_bgface, + bgface_to_lcell, + face_is_owner, + face_to_subface, + ncglue) + + trian = Gridap.Geometry.BodyFittedTriangulation(model,bgface_grid,face_to_bgface) + BoundaryTriangulation(trian,glue) +end + +function Gridap.Geometry.SkeletonTriangulation(model::Gridap.Adaptivity.AdaptedDiscreteModel{D}, + ncglue::GridapP4est.NonConformingGlue{D}) where {D} + trian = SkeletonTriangulation(Gridap.Adaptivity.get_model(model),ncglue) + return Gridap.Adaptivity.AdaptedTriangulation(trian,model) +end + +function Gridap.Geometry.SkeletonTriangulation( + portion,model::OctreeDistributedDiscreteModel{Dc};kwargs...) where Dc + gids = get_face_gids(model.dmodel,Dc) + trians = map(local_views(model.dmodel), + partition(gids), + model.non_conforming_glue) do model, gids, ncglue + trian=SkeletonTriangulation(model,ncglue) + GridapDistributed.filter_cells_when_needed(portion,gids,trian) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + +struct SubfaceRefCoordsMap{A} <: Gridap.Arrays.Map + face_is_owner :: Vector{Bool} + face_to_subface :: Vector{Int8} + ref_rule :: A +end + +function Gridap.Arrays.return_cache(f::SubfaceRefCoordsMap,gface::Integer) + ref_grid=Gridap.Adaptivity.get_ref_grid(f.ref_rule) + poly=Gridap.Adaptivity.get_polytope(f.ref_rule) + poly_vertices=Gridap.Geometry.get_vertex_coordinates(poly) + cell_coordinates=get_cell_coordinates(ref_grid) + cache_cell_coordinates=array_cache(cell_coordinates) + poly_vertices,cell_coordinates,cache_cell_coordinates +end + +function Gridap.Arrays.evaluate!(cache,f::SubfaceRefCoordsMap,gface::Integer) + poly_vertices,cell_coordinates,cache_cell_coordinates=cache + if (f.face_is_owner[gface]) + subface=f.face_to_subface[gface] + getindex!(cache_cell_coordinates,cell_coordinates,subface) + else + poly_vertices + end +end + +function Gridap.Geometry.get_glue(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlueNonConforming}, + ::Val{D},::Val{D}) where {D,Dc,Dp,A} + tface_to_mface = trian.glue.face_to_cell_glue.face_to_cell + + face_to_q_vertex_coords = + Gridap.Geometry._compute_face_to_q_vertex_coords(trian,trian.glue.face_to_cell_glue) + f(p) = + Gridap.ReferenceFEs.get_shapefuns(Gridap.ReferenceFEs.LagrangianRefFE(Float64,Gridap.ReferenceFEs.get_polytope(p),1)) + ftype_to_shapefuns = map( f, Gridap.Geometry.get_reffes(trian) ) + face_to_shapefuns = Gridap.ReferenceFEs.expand_cell_data(ftype_to_shapefuns,trian.glue.face_to_cell_glue.face_to_ftype) + face_s_q = lazy_map(Gridap.Fields.linear_combination,face_to_q_vertex_coords,face_to_shapefuns) + + # Map subface to cell ref space + ref_rule=nothing + if Dc==1 + ref_rule=Gridap.Adaptivity.RefinementRule(SEGMENT,2) + else + @assert Dc==2 + ref_rule=Gridap.Adaptivity.RedRefinementRule(QUAD) + end + + m=SubfaceRefCoordsMap(trian.glue.face_is_owner,trian.glue.face_to_subface,ref_rule) + subface_ref_coords = lazy_map(m, Gridap.Arrays.IdentityVector(length(trian.glue.face_to_subface)) ) + face_to_q_vertex_coords = lazy_map(evaluate,face_s_q,subface_ref_coords) + face_s_q = lazy_map(Gridap.Fields.linear_combination,face_to_q_vertex_coords,face_to_shapefuns) + + tface_to_mface_map = face_s_q + mface_to_tface = nothing + Gridap.Geometry.FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface) +end + +function Gridap.Geometry.get_facet_normal(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlueNonConforming}) where {Dc,Dp,A} + Gridap.Geometry.get_facet_normal(trian, trian.glue.face_to_cell_glue) +end diff --git a/src/GridapP4est.jl b/src/GridapP4est.jl index bcb6c15..ae52258 100644 --- a/src/GridapP4est.jl +++ b/src/GridapP4est.jl @@ -16,10 +16,12 @@ module GridapP4est include("PXestTypeMethods.jl") include("UniformlyRefinedForestOfOctreesDiscreteModels.jl") include("OctreeDistributedDiscreteModels.jl") + include("Geometry.jl") include("AnisotropicallyAdapted3DDistributedDiscreteModels.jl") include("GridapFixes.jl") include("FESpaces.jl") include("AdaptivityFlagsMarkingStrategies.jl") + export UniformlyRefinedForestOfOctreesDiscreteModel export OctreeDistributedDiscreteModel diff --git a/src/OctreeDistributedDiscreteModels.jl b/src/OctreeDistributedDiscreteModels.jl index ee58b56..6a7bd72 100644 --- a/src/OctreeDistributedDiscreteModels.jl +++ b/src/OctreeDistributedDiscreteModels.jl @@ -3,36 +3,272 @@ const nothing_flag = Cint(0) const refine_flag = Cint(1) const coarsen_flag = Cint(2) -struct NonConformingGlue{Dc,A,B,C} - num_regular_faces :: A # <:AbstractVector{<:AbstractVector{<:Integer}} - num_hanging_faces :: B # <:AbstractVector{<:AbstractVector{<:Integer}} - hanging_faces_glue :: C # <:AbstractVector{<:AbstractVector{<:Tuple{<:Integer,:Integer,Integer}}}} - function NonConformingGlue(num_regular_faces,num_hanging_faces,hanging_faces_glue) +struct NonConformingGlue{Dc,A,B,C,D,E,F,G} + num_regular_faces :: A # <:AbstractVector{<:AbstractVector{<:Integer}} + num_hanging_faces :: B # <:AbstractVector{<:AbstractVector{<:Integer}} + hanging_faces_glue :: C # <:AbstractVector{<:AbstractVector{<:Tuple{<:Integer,:Integer,Integer}}}} + hanging_faces_to_cell :: D # <:AbstractVector{<:AbstractVector{<:Integer}} + hanging_faces_to_lface :: E # <:AbstractVector{<:AbstractVector{<:Integer}} + owner_faces_pindex :: F # <:AbstractVector{<:AbstractVector{<:Integer}} + owner_faces_lids :: G # <:AbstractVector{<:Dict{Int,Tuple{Int,Int,Int}}} + function NonConformingGlue(num_regular_faces, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids) Dc = length(num_regular_faces) @assert length(num_hanging_faces)==Dc @assert length(hanging_faces_glue)==Dc A = typeof(num_regular_faces) B = typeof(num_hanging_faces) C = typeof(hanging_faces_glue) - new{Dc,A,B,C}(num_regular_faces,num_hanging_faces,hanging_faces_glue) + D = typeof(hanging_faces_to_cell) + E = typeof(hanging_faces_to_lface) + F = typeof(owner_faces_pindex) + G = typeof(owner_faces_lids) + new{Dc,A,B,C,D,E,F,G}(num_regular_faces, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids) end end +function _compute_owner_faces_lids(Df,Dc,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) + if (ocell_dim == Df) + ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) + owner_face = cell_faces[ocell][ocell_lface_within_dim] + if !(haskey(owner_faces_lids, owner_face)) + num_owner_faces += 1 + owner_faces_lids[owner_face] = (num_owner_faces, ocell, ocell_lface) + end + end + end + end + owner_faces_lids +end + +function _subface_to_face_corners(::PXestUniformRefinementRuleType, subface) + (subface,) +end + +function _subface_to_face_corners(::PXestHorizontalRefinementRuleType, subface) + @assert subface==1 || subface==2 + if (subface==1) + (1,2) + else + (3,4) + end +end + +function _subface_to_face_corners(::PXestVerticalRefinementRuleType, subface) + @assert subface==1 || subface==2 + if (subface==1) + (1,3) + else + (2,4) + end +end + + +# count how many different owner faces +# for each owner face +# track the global IDs of its face vertices from the perspective of the subfaces +# for each owner face +# compute permutation id +function _compute_owner_faces_pindex_and_lids(Dc, + pXest_refinement_rule, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + cell_vertices, + cell_faces, + lface_to_cvertices, + pindex_to_cfvertex_to_fvertex) + + owner_faces_lids =_compute_owner_faces_lids(Dc-1,Dc, + num_hanging_faces, + hanging_faces_glue, + cell_faces) + + @debug "owner_faces_lids [Df=$(Dc-1) Dc=$(Dc)]: $(owner_faces_lids)" + + num_owner_faces = length(keys(owner_faces_lids)) + num_face_vertices = length(first(lface_to_cvertices)) + owner_face_vertex_ids = Vector{Int}(undef, num_face_vertices * num_owner_faces) + owner_face_vertex_ids .= -1 + + for fid_hanging = 1:num_hanging_faces + ocell, ocell_lface, subface = hanging_faces_glue[fid_hanging] + @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) + 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) + 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) + cvertex = lface_to_cvertices[lface][fcorner] + vertex = cell_vertices[cell][cvertex] + @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: cell=$(cell) lface=$(lface) cvertex=$(cvertex) vertex=$(vertex) owner_face=$(owner_face) owner_face_lid=$(owner_face_lid)" + @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: owner_face_vertex_ids[$((owner_face_lid-1)*num_face_vertices+fcorner)] = $(vertex)" + owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fcorner] = vertex + end + end + end + end + @debug "owner_face_vertex_ids [Dc=$(Dc)]: $(owner_face_vertex_ids)" + + 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) + # Compute permutation id by comparing + # 1. cell_vertices[ocell][ocell_lface] + # 2. owner_face_vertex_ids + pindexfound = false + cfvertex_to_cvertex = lface_to_cvertices[ocell_lface_within_dim] + for (pindex, cfvertex_to_fvertex) in enumerate(pindex_to_cfvertex_to_fvertex) + found = true + for (cfvertex, fvertex) in enumerate(cfvertex_to_fvertex) + vertex1 = owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fvertex] + cvertex = cfvertex_to_cvertex[cfvertex] + vertex2 = cell_vertices[ocell][cvertex] + @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: cell_vertices[$(ocell)][$(cvertex)]=$(cell_vertices[ocell][cvertex]) owner_face_vertex_ids[$((owner_face_lid-1)*num_face_vertices+fvertex)]=$(owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fvertex])" + # -1 can only happen in the interface of two + # ghost cells at different refinement levels + if (vertex1 != vertex2) && (vertex1 != -1) + found = false + break + end + end + if found + owner_faces_pindex[owner_face_lid] = pindex + pindexfound = true + break + end + end + @assert pindexfound "Valid pindex not found" + end + @debug "owner_faces_pindex: $(owner_faces_pindex)" + + owner_faces_pindex, owner_faces_lids +end + + + +function _compute_owner_edges_pindex_and_lids( + num_hanging_edges, + hanging_edges_glue, + hanging_edges_to_cell, + hanging_edges_to_ledge, + cell_vertices, + cell_edges) + Dc=3 + owner_edges_lids =_compute_owner_faces_lids(1, + Dc, + num_hanging_edges, + hanging_edges_glue, + cell_edges) + + num_owner_edges = length(keys(owner_edges_lids)) + owner_edges_pindex = Vector{Int}(undef, num_owner_edges) + + ledge_to_cvertices = Gridap.ReferenceFEs.get_faces(HEX, 1, 0) + + # Go over hanging edges + # 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) + if (ocell!=-1 && ocell_dim==1) + ocell_ledge_within_dim = face_lid_within_dim(Val{Dc}, 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]] + gvertex2 = cell_vertices[ocell][ledge_to_cvertices[ocell_ledge_within_dim][subedge]] + @debug "fid_hanging=$(fid_hanging) cell=$(cell) ledge=$(ledge) ocell=$(ocell) ocell_ledge=$(ocell_ledge) subedge=$(subedge) gvertex1=$(gvertex1) gvertex2=$(gvertex2)" + pindex = gvertex1==gvertex2 ? 1 : 2 + owner_edge=cell_edges[ocell][ocell_ledge_within_dim] + owner_edge_lid, _ = owner_edges_lids[owner_edge] + owner_edges_pindex[owner_edge_lid]=pindex + end + end + owner_edges_pindex, owner_edges_lids +end + function _create_conforming_model_non_conforming_glue(model::GridapDistributed.DistributedDiscreteModel{Dc}) where Dc non_conforming_glue = map(local_views(model)) do model num_regular_faces=Vector{Int}(undef,Dc) num_hanging_faces=Vector{Int}(undef,Dc) hanging_faces_glue=Vector{Tuple{Int,Int,Int}}(undef,Dc) + hanging_faces_to_cell=Vector{Vector{Int}}(undef,Dc) + hanging_faces_to_lface=Vector{Vector{Int}}(undef,Dc) + owner_faces_pindex=Vector{Vector{Int}}(undef,Dc) + owner_faces_lids=Vector{Dict{Int,Tuple{Int,Int,Int}}}(undef,Dc) for d=1:Dc num_regular_faces[d]=num_faces(model,d-1) num_hanging_faces[d]=0 + hanging_faces_to_cell[d]=Int[] + hanging_faces_to_lface[d]=Int[] + owner_faces_pindex[d]=Int[] + owner_faces_lids[d]=Dict{Int,Tuple{Int,Int,Int}}() end - NonConformingGlue(num_regular_faces,num_hanging_faces,hanging_faces_glue) + NonConformingGlue(num_regular_faces, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids) end return non_conforming_glue end +function _generate_hanging_faces_to_cell_and_lface(num_regular_faces, + num_hanging_faces, + gridap_cell_faces) + # Locate for each hanging face the owner cell among the set of cells + # to which it belongs and local position within that cell + # By convention, the owner cell is the cell with minimum identifer among + # all in the set. This convention is used here, and in other parts of the code. + # Breaking this convention here may affect the consistency of other parts of the code. + hanging_faces_to_cell = Vector{Int}(undef, num_hanging_faces) + hanging_faces_to_lface = Vector{Int}(undef, num_hanging_faces) + hanging_faces_to_cell .= -1 + for cell = 1:length(gridap_cell_faces) + s = gridap_cell_faces.ptrs[cell] + e = gridap_cell_faces.ptrs[cell+1] + l = e - s + for j = 1:l + fid = gridap_cell_faces.data[s+j-1] + if fid > num_regular_faces + fid_hanging = fid - num_regular_faces + if (hanging_faces_to_cell[fid_hanging]==-1) + hanging_faces_to_cell[fid_hanging] = cell + hanging_faces_to_lface[fid_hanging] = j + end + end + end + end + hanging_faces_to_cell, hanging_faces_to_lface +end + mutable struct OctreeDistributedDiscreteModel{Dc,Dp,A,B,C,D,E,F} <: GridapDistributed.DistributedDiscreteModel{Dc,Dp} parts :: A dmodel :: B @@ -61,7 +297,7 @@ mutable struct OctreeDistributedDiscreteModel{Dc,Dp,A,B,C,D,E,F} <: GridapDistri ptr_pXest_connectivity, ptr_pXest, pXest_type::PXestType, - pXest_refinement_rule_type::Union{Nothing,PXestRefinementRuleType}, + pXest_refinement_rule_type::PXestRefinementRuleType, owns_ptr_pXest_connectivity::Bool, gc_ref) @@ -158,14 +394,14 @@ function OctreeDistributedDiscreteModel(parts::AbstractVector{<:Integer}, ptr_pXest_connectivity, ptr_pXest, pXest_type, - nothing, + PXestUniformRefinementRuleType(), true, nothing) else ## HUGE WARNING: Shouldn't we provide here the complementary of parts ## instead of parts? Otherwise, when calling _free!(...) ## we cannot trust on parts. - return VoidOctreeDistributedDiscreteModel(coarse_model,parts) + return VoidOctreeDistributedDiscreteModel(coarse_model,parts,PXestUniformRefinementRuleType()) end end @@ -179,7 +415,9 @@ end const VoidOctreeDistributedDiscreteModel{Dc,Dp,A,C,D} = OctreeDistributedDiscreteModel{Dc,Dp,A,Nothing,C,D,Nothing} -function VoidOctreeDistributedDiscreteModel(coarse_model::DiscreteModel{Dc,Dp},parts) where {Dc,Dp} +function VoidOctreeDistributedDiscreteModel(coarse_model::DiscreteModel{Dc,Dp}, + parts, + pXest_refinement_rule_type::PXestRefinementRuleType) where {Dc,Dp} ptr_pXest_connectivity = setup_pXest_connectivity(coarse_model) OctreeDistributedDiscreteModel(Dc, Dp, @@ -190,7 +428,7 @@ function VoidOctreeDistributedDiscreteModel(coarse_model::DiscreteModel{Dc,Dp},p ptr_pXest_connectivity, nothing, _dim_to_pXest_type(Dc), - nothing, + pXest_refinement_rule_type, true, nothing) end @@ -205,7 +443,7 @@ function VoidOctreeDistributedDiscreteModel(model::OctreeDistributedDiscreteMode model.ptr_pXest_connectivity, nothing, _dim_to_pXest_type(Dc), - nothing, + model.pXest_refinement_rule_type, false, model) end @@ -1380,6 +1618,7 @@ function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, # Build fine-grid mesh fmodel,non_conforming_glue = setup_non_conforming_distributed_discrete_model(model.pXest_type, + model.pXest_refinement_rule_type, model.parts, model.coarse_model, model.ptr_pXest_connectivity, @@ -1672,6 +1911,7 @@ function _redistribute_parts_subseteq_parts_redistributed(model::OctreeDistribut # Build fine-grid mesh fmodel, non_conforming_glue = setup_non_conforming_distributed_discrete_model(model.pXest_type, + model.pXest_refinement_rule_type, parts, model.coarse_model, model.ptr_pXest_connectivity, @@ -1904,6 +2144,7 @@ function num_cell_faces(::Type{Val{Dc}}) where Dc end function setup_non_conforming_distributed_discrete_model(pXest_type::PXestType, + pXest_refinement_rule_type::PXestRefinementRuleType, parts, coarse_discrete_model, ptr_pXest_connectivity, @@ -1917,7 +2158,7 @@ function setup_non_conforming_distributed_discrete_model(pXest_type::PXestType, gridap_cell_faces, non_conforming_glue= - generate_cell_faces_and_non_conforming_glue(pXest_type,ptr_pXest_lnodes, cell_prange) + generate_cell_faces_and_non_conforming_glue(pXest_type,pXest_refinement_rule_type,ptr_pXest_lnodes, cell_prange) cell_corner_lids = map(gridap_cell_faces[1]) do cell_lids Gridap.Arrays.Table(cell_lids.data,cell_lids.ptrs) # JaggedArray -> Table @@ -2004,5 +2245,3 @@ function _set_hanging_labels!(face_labeling,non_conforming_glue,coarse_face_labe end - - diff --git a/src/PXestTypeMethods.jl b/src/PXestTypeMethods.jl index d20c370..4317dfd 100644 --- a/src/PXestTypeMethods.jl +++ b/src/PXestTypeMethods.jl @@ -1061,11 +1061,11 @@ function pXest_lnodes_decode(::P8estType,face_code, hanging_face, hanging_edge) p8est_lnodes_decode(face_code, hanging_face, hanging_edge) end -function num_cell_dims(::P4estType) +function Gridap.Geometry.num_cell_dims(::P4estType) 2 end -function num_cell_dims(::P6P8estType) +function Gridap.Geometry.num_cell_dims(::P6P8estType) 3 end @@ -1255,7 +1255,8 @@ end const p6est_half_to_regular_vertices = [ 0 1; 2 3; 0 2; 1 3] -function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType, +function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType, + pXest_refinement_rule::PXestRefinementRuleType, ptr_pXest_lnodes, cell_prange) @@ -1317,7 +1318,11 @@ function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType, num_regular_faces, num_hanging_faces, gridap_cell_faces, - hanging_faces_glue = + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids = map(partition(cell_prange), element_nodes_with_ghosts, face_code_with_ghosts) do indices, @@ -1874,11 +1879,54 @@ function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType, end hanging_faces_glue[Dc] = hanging_faces_owner_cell_and_lface + hanging_faces_to_cell = Vector{Vector{Int}}(undef, Dc) + hanging_faces_to_lface = Vector{Vector{Int}}(undef, Dc) + for i=1:Dc + # Locate for each hanging facet a cell to which it belongs + # and local position within that cell + hanging_faces_to_cell[i], + hanging_faces_to_lface[i] = + _generate_hanging_faces_to_cell_and_lface(num_regular_faces[i], + num_hanging_faces[i], + gridap_cell_faces[i]) + end + + owner_faces_pindex = Vector{Vector{Int}}(undef, Dc - 1) + owner_faces_lids = Vector{Dict{Int,Tuple{Int,Int,Int}}}(undef, Dc - 1) + + lface_to_cvertices = Gridap.ReferenceFEs.get_faces(Dc == 2 ? QUAD : HEX, Dc - 1, 0) + pindex_to_cfvertex_to_fvertex = Gridap.ReferenceFEs.get_vertex_permutations(Dc == 2 ? SEGMENT : QUAD) + + owner_faces_pindex[Dc-1], owner_faces_lids[Dc-1] = _compute_owner_faces_pindex_and_lids(Dc, + pXest_refinement_rule, + num_hanging_faces[Dc], + hanging_faces_glue[Dc], + hanging_faces_to_cell[Dc], + hanging_faces_to_lface[Dc], + gridap_cell_faces[1], + gridap_cell_faces[Dc], + lface_to_cvertices, + pindex_to_cfvertex_to_fvertex) + + if (Dc == 3) + owner_faces_pindex[1], owner_faces_lids[1]= + _compute_owner_edges_pindex_and_lids( + num_hanging_faces[2], + hanging_faces_glue[2], + hanging_faces_to_cell[2], + hanging_faces_to_lface[2], + gridap_cell_faces[1], + gridap_cell_faces[2]) + end return num_regular_faces, num_hanging_faces, gridap_cell_faces, - hanging_faces_glue + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids end |> tuple_of_arrays @@ -1889,8 +1937,14 @@ function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType, gridap_cell_faces[i] end end - non_conforming_glue=map(num_regular_faces,num_hanging_faces,hanging_faces_glue) do nrf, nhf, hfg - NonConformingGlue(nrf, nhf, hfg) + non_conforming_glue=map(num_regular_faces, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + owner_faces_pindex, + owner_faces_lids) do nrf, nhf, hfg, hfc, hfl, ofp,ofl + NonConformingGlue(nrf, nhf, hfg, hfc, hfl, ofp, ofl) end gridap_cell_faces_out,non_conforming_glue end diff --git a/test/OctreeDistributedDiscreteModelsTests.jl b/test/OctreeDistributedDiscreteModelsTests.jl index c6446dd..30e0f90 100644 --- a/test/OctreeDistributedDiscreteModelsTests.jl +++ b/test/OctreeDistributedDiscreteModelsTests.jl @@ -64,9 +64,11 @@ module OctreeDistributedDiscreteModelsTests level_parts = generate_level_parts(ranks,num_parts_x_level) coarse_model = CartesianDiscreteModel(domain,nc) model = OctreeDistributedDiscreteModel(level_parts[2],coarse_model,1) - vmodel1 = GridapP4est.VoidOctreeDistributedDiscreteModel(model,ranks) - vmodel2 = GridapP4est.VoidOctreeDistributedDiscreteModel(coarse_model,ranks) - + vmodel1 = GridapP4est.VoidOctreeDistributedDiscreteModel(model, + ranks) + vmodel2 = GridapP4est.VoidOctreeDistributedDiscreteModel(coarse_model, + ranks, + GridapP4est.PXestUniformRefinementRuleType()) # Refining and distributing fmodel , rglue = refine(model,parts=level_parts[1]) dfmodel, dglue = redistribute(fmodel) diff --git a/test/PoissonNonConformingOctreeModelsTests.jl b/test/PoissonNonConformingOctreeModelsTests.jl index 4670669..b473628 100644 --- a/test/PoissonNonConformingOctreeModelsTests.jl +++ b/test/PoissonNonConformingOctreeModelsTests.jl @@ -32,12 +32,23 @@ module PoissonNonConformingOctreeModelsTests u,f end - function test_transfer_ops_and_redistribute(ranks,dmodel,order,T::Type) + function test_transfer_ops_and_redistribute(ranks,dmodel,order,cg_or_dg,T::Type) + @assert cg_or_dg == :cg || cg_or_dg == :dg + if (cg_or_dg==:dg) + @assert T==Float64 + end + + conformity=cg_or_dg==:cg ? :H1 : :L2 + # Define manufactured functions u,f = generate_analytical_problem_functions(T,order) degree = 2*order+1 reffe=ReferenceFE(lagrangian,T,order) - VH=FESpace(dmodel,reffe;dirichlet_tags="boundary") + if (cg_or_dg == :cg) + VH=FESpace(dmodel,reffe,conformity=conformity;dirichlet_tags="boundary") + else + VH=FESpace(dmodel,reffe,conformity=conformity) + end UH=TrialFESpace(VH,u) ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices flags=zeros(Cint,length(indices)) @@ -58,16 +69,46 @@ module PoissonNonConformingOctreeModelsTests # print(glue.n2o_faces_map[end]); print("\n") # end # end - Vh=FESpace(fmodel,reffe,conformity=:H1;dirichlet_tags="boundary") + if (cg_or_dg == :cg) + Vh=FESpace(fmodel,reffe,conformity=conformity;dirichlet_tags="boundary") + else + Vh=FESpace(fmodel,reffe,conformity=conformity) + end Uh=TrialFESpace(Vh,u) ΩH = Triangulation(dmodel) dΩH = Measure(ΩH,degree) - aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH - bH(v) = ∫(v⋅f)*dΩH + if (cg_or_dg==:cg) + aHcg(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH + bHcg(v) = ∫(v⋅f)*dΩH + op = AffineFEOperator(aHcg,bHcg,UH,VH) + else + h = 2 + γ = 10 + + ΛH = Skeleton(dmodel) + ΓH = Boundary(dmodel,tags="boundary") + + dΛH = Measure(ΛH,2*order) + dΓH = Measure(ΓH,2*order) + + n_ΓH = get_normal_vector(ΓH) + n_ΛH = get_normal_vector(ΛH) + + aHdg(u,v) = + ∫( ∇(v)⋅∇(u) )*dΩH + + ∫( (γ/h)*v*u - v*(n_ΓH⋅∇(u)) - (n_ΓH⋅∇(v))*u )*dΓH + + ∫( (γ/h)*jump(v*n_ΛH)⋅jump(u*n_ΛH) - + jump(v*n_ΛH)⋅mean(∇(u)) - + mean(∇(v))⋅jump(u*n_ΛH) )*dΛH + + bHdg(v) = + ∫( v*f )*dΩH + + ∫( (γ/h)*v*u - (n_ΓH⋅∇(v))*u )*dΓH + op = AffineFEOperator(aHdg,bHdg,UH,VH) + end - op = AffineFEOperator(aH,bH,UH,VH) uH = solve(op) e = u - uH @@ -85,10 +126,36 @@ module PoissonNonConformingOctreeModelsTests Ωh = Triangulation(fmodel) dΩh = Measure(Ωh,degree) - ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh - bh(v) = ∫(v⋅f)*dΩh + if (cg_or_dg==:cg) + ahcg(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh + bhcg(v) = ∫(v⋅f)*dΩh + op = AffineFEOperator(ahcg,bhcg,Uh,Vh) + else + h = 2 + γ = 10 + + Λh = Skeleton(fmodel) + Γh = Boundary(fmodel,tags="boundary") + + dΛh = Measure(Λh,2*order) + dΓh = Measure(Γh,2*order) + + n_Γh = get_normal_vector(Γh) + n_Λh = get_normal_vector(Λh) + + ahdg(u,v) = + ∫( ∇(v)⋅∇(u) )*dΩh + + ∫( (γ/h)*v*u - v*(n_Γh⋅∇(u)) - (n_Γh⋅∇(v))*u )*dΓh + + ∫( (γ/h)*jump(v*n_Λh)⋅jump(u*n_Λh) - + jump(v*n_Λh)⋅mean(∇(u)) - + mean(∇(v))⋅jump(u*n_Λh) )*dΛh + + bhdg(v) = + ∫( v*f )*dΩh + + ∫( (γ/h)*v*u - (n_Γh⋅∇(v))*u )*dΓh + op = AffineFEOperator(ahdg,bhdg,Uh,Vh) + end - op = AffineFEOperator(ah,bh,Uh,Vh) uh = solve(op) e = u - uh @@ -146,19 +213,51 @@ module PoissonNonConformingOctreeModelsTests ones(Cint,num_cells(lmodel)) end end - fmodel_red, red_glue=GridapDistributed.redistribute(fmodel,weights=weights); - Vhred=FESpace(fmodel_red,reffe,conformity=:H1;dirichlet_tags="boundary") - Uhred=TrialFESpace(Vhred,u) + fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); + if (cg_or_dg==:cg) + Vhred=FESpace(fmodel_red,reffe,conformity=conformity;dirichlet_tags="boundary") + Uhred=TrialFESpace(Vhred,u) + else + Vhred=FESpace(fmodel_red,reffe,conformity=conformity) + Uhred=TrialFESpace(Vhred,u) + end Ωhred = Triangulation(fmodel_red) dΩhred = Measure(Ωhred,degree) - ahred(u,v) = ∫( ∇(v)⊙∇(u) )*dΩhred - bhred(v) = ∫(v⋅f)*dΩhred + if (cg_or_dg==:cg) + ahcgred(u,v) = ∫( ∇(v)⊙∇(u) )*dΩhred + bhcgred(v) = ∫(v⋅f)*dΩhred + op = AffineFEOperator(ahcgred,bhcgred,Uhred,Vhred) + else + h = 2 + γ = 10 + + Λhred = Skeleton(fmodel_red) + Γhred = Boundary(fmodel_red,tags="boundary") + + dΛhred = Measure(Λhred,2*order) + dΓhred = Measure(Γhred,2*order) + + n_Γhred = get_normal_vector(Γhred) + n_Λhred = get_normal_vector(Λhred) + + ahdgred(u,v) = + ∫( ∇(v)⋅∇(u) )*dΩhred + + ∫( (γ/h)*v*u - v*(n_Γhred⋅∇(u)) - (n_Γhred⋅∇(v))*u )*dΓhred + + ∫( (γ/h)*jump(v*n_Λhred)⋅jump(u*n_Λhred) - + jump(v*n_Λhred)⋅mean(∇(u)) - + mean(∇(v))⋅jump(u*n_Λhred) )*dΛhred + + bhdgred(v) = + ∫( v*f )*dΩhred + + ∫( (γ/h)*v*u - (n_Γhred⋅∇(v))*u )*dΓhred + op = AffineFEOperator(ahdgred,bhdgred,Uhred,Vhred) + end - op = AffineFEOperator(ahred,bhred,Uhred,Vhred) - uhred = solve(op) - e = u - uhred + + uhred = solve(op) + e = u - uhred el2 = sqrt(sum( ∫( e⋅e )*dΩhred )) println("[SOLVE FINE REDISTRIBUTED] el2 < tol: $(el2) < $(tol)") @assert el2 < tol @@ -176,8 +275,14 @@ module PoissonNonConformingOctreeModelsTests function test_refine_and_coarsen_at_once(ranks, dmodel::OctreeDistributedDiscreteModel{Dc}, order, + cg_or_dg, T::Type) where Dc + @assert cg_or_dg == :cg || cg_or_dg == :dg + if (cg_or_dg==:dg) + @assert T==Float64 + end + # Define manufactured functions u,f = generate_analytical_problem_functions(T,order) @@ -194,19 +299,55 @@ module PoissonNonConformingOctreeModelsTests end fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); + conformity=cg_or_dg==:cg ? :H1 : :L2 + reffe=ReferenceFE(lagrangian,T,order) - VH=FESpace(dmodel,reffe,conformity=:H1;dirichlet_tags="boundary") + if (conformity==:L2) + VH=FESpace(dmodel,reffe,conformity=conformity) + else + VH=FESpace(dmodel,reffe,conformity=conformity;dirichlet_tags="boundary") + end UH=TrialFESpace(VH,u) - Vh=FESpace(fmodel,reffe,conformity=:H1;dirichlet_tags="boundary") + if (conformity==:L2) + Vh=FESpace(fmodel,reffe,conformity=conformity) + else + Vh=FESpace(fmodel,reffe,conformity=conformity;dirichlet_tags="boundary") + end Uh=TrialFESpace(Vh,u) ΩH = Triangulation(dmodel) dΩH = Measure(ΩH,degree) - aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH - bH(v) = ∫(v⋅f)*dΩH + if (cg_or_dg==:cg) + aHcg(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH + bHcg(v) = ∫(v⋅f)*dΩH + op = AffineFEOperator(aHcg,bHcg,UH,VH) + else + h = 2 + γ = 10 + + ΛH = Skeleton(dmodel) + ΓH = Boundary(dmodel,tags="boundary") + + dΛH = Measure(ΛH,2*order) + dΓH = Measure(ΓH,2*order) + + n_ΓH = get_normal_vector(ΓH) + n_ΛH = get_normal_vector(ΛH) + + aHdg(u,v) = + ∫( ∇(v)⋅∇(u) )*dΩH + + ∫( (γ/h)*v*u - v*(n_ΓH⋅∇(u)) - (n_ΓH⋅∇(v))*u )*dΓH + + ∫( (γ/h)*jump(v*n_ΛH)⋅jump(u*n_ΛH) - + jump(v*n_ΛH)⋅mean(∇(u)) - + mean(∇(v))⋅jump(u*n_ΛH) )*dΛH + + bHdg(v) = + ∫( v*f )*dΩH + + ∫( (γ/h)*v*u - (n_ΓH⋅∇(v))*u )*dΓH + op = AffineFEOperator(aHdg,bHdg,UH,VH) + end - op = AffineFEOperator(aH,bH,UH,VH) uH = solve(op) e = u - uH @@ -221,10 +362,39 @@ module PoissonNonConformingOctreeModelsTests Ωh = Triangulation(fmodel) dΩh = Measure(Ωh,degree) - ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh - bh(v) = ∫(v⋅f)*dΩh + if (cg_or_dg==:cg) + ahcg(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh + bhcg(v) = ∫(v⋅f)*dΩh + op = AffineFEOperator(ahcg,bhcg,Uh,Vh) + else + h = 2 + γ = 10 + + Λh = Skeleton(fmodel) + Γh = Boundary(fmodel,tags="boundary") + + dΛh = Measure(Λh,2*order) + dΓh = Measure(Γh,2*order) + + n_Γh = get_normal_vector(Γh) + n_Λh = get_normal_vector(Λh) + + ahdg(u,v) = + ∫( ∇(v)⋅∇(u) )*dΩh + + ∫( (γ/h)*v*u - v*(n_Γh⋅∇(u)) - (n_Γh⋅∇(v))*u )*dΓh + + # with p=4 MPI tasks, (γ/h)*jump(v*n_Λh)⋅jump(u*n_Λh) fails! + # did not have time you to understand the cause, workaround + # is the line below + ∫( (γ/h)*(jump(v*n_Λh)⋅jump(u*n_Λh)) - + jump(v*n_Λh)⋅mean(∇(u)) - + mean(∇(v))⋅jump(u*n_Λh) )*dΛh + + bhdg(v) = + ∫( v*f )*dΩh + + ∫( (γ/h)*v*u - (n_Γh⋅∇(v))*u )*dΓh + op = AffineFEOperator(ahdg,bhdg,Uh,Vh) + end - op = AffineFEOperator(ah,bh,Uh,Vh) uh = solve(op) e = u - uh @@ -242,30 +412,30 @@ module PoissonNonConformingOctreeModelsTests @assert el2 < tol end - function test_2d(ranks,order,T::Type;num_amr_steps=5) + function test_2d(ranks,order,cg_or_dg,T::Type;num_amr_steps=5) coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2) - test_refine_and_coarsen_at_once(ranks,dmodel,order,T) + test_refine_and_coarsen_at_once(ranks,dmodel,order,cg_or_dg,T) rdmodel=dmodel for i=1:num_amr_steps - rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,T) + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,cg_or_dg,T) end end - function test_3d(ranks,order,T::Type;num_amr_steps=5) + function test_3d(ranks,order,cg_or_dg,T::Type;num_amr_steps=5) coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(1,1,1)) dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2) - test_refine_and_coarsen_at_once(ranks,dmodel,order,T) + test_refine_and_coarsen_at_once(ranks,dmodel,order,cg_or_dg,T) rdmodel=dmodel for i=1:num_amr_steps - rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,T) + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,cg_or_dg,T) end end - function test(ranks,TVDc::Type{Val{Dc}}, perm, order,T::Type) where Dc + function test(ranks,TVDc::Type{Val{Dc}}, perm, order,cg_or_dg,T::Type) where Dc coarse_model = setup_model(TVDc,perm) model = OctreeDistributedDiscreteModel(ranks, coarse_model, 1) - test_transfer_ops_and_redistribute(ranks,model,order,T) + test_transfer_ops_and_redistribute(ranks,model,order,cg_or_dg,T) end function _field_type(::Val{Dc}, scalar_or_vector::Symbol) where Dc @@ -281,16 +451,19 @@ module PoissonNonConformingOctreeModelsTests #debug_logger = ConsoleLogger(stderr, Logging.Debug) #global_logger(debug_logger); # Enable the debug logger globally ranks = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),))) - for Dc=3:3, perm=1:4, order=1:4, scalar_or_vector in (:scalar,) - test(ranks,Val{Dc},perm,order,_field_type(Val{Dc}(),scalar_or_vector)) + for Dc=2:3, perm in (1,2,4), order=(1,2), scalar_or_vector in (:scalar,) + test(ranks,Val{Dc},perm,order,:dg,_field_type(Val{Dc}(),scalar_or_vector)) end for Dc=2:3, perm in (1,2), order in (1,4), scalar_or_vector in (:vector,) - test(ranks,Val{Dc},perm,order,_field_type(Val{Dc}(),scalar_or_vector)) + test(ranks,Val{Dc},perm,order,:cg,_field_type(Val{Dc}(),scalar_or_vector)) + end + for order=2:2, scalar_or_vector in (:scalar,) + test_2d(ranks,order,:dg,_field_type(Val{2}(),scalar_or_vector), num_amr_steps=5) + test_3d(ranks,order,:dg,_field_type(Val{3}(),scalar_or_vector), num_amr_steps=4) end for order=2:2, scalar_or_vector in (:scalar,:vector) - test_2d(ranks,order,_field_type(Val{2}(),scalar_or_vector), num_amr_steps=5) - test_3d(ranks,order,_field_type(Val{3}(),scalar_or_vector), num_amr_steps=4) + test_2d(ranks,order,:cg,_field_type(Val{2}(),scalar_or_vector), num_amr_steps=5) + test_3d(ranks,order,:cg,_field_type(Val{3}(),scalar_or_vector), num_amr_steps=4) end end end -