diff --git a/Project.toml b/Project.toml index edf8d4a..1c7724c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridapP4est" uuid = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9" authors = ["Alberto F. Martin "] -version = "0.3.6" +version = "0.3.9" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" @@ -17,8 +17,8 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] ArgParse = "1" FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 1" -Gridap = "0.17.22, 0.18" -GridapDistributed = "0.3.1, 0.4" +Gridap = "0.18.6" +GridapDistributed = "0.4" MPI = "0.20" P4est_wrapper = "0.2.2" PartitionedArrays = "0.3.3" diff --git a/README.md b/README.md index dab97de..75124e8 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,8 @@ [![Build Status](https://github.com/gridap/GridapP4est.jl/workflows/CI/badge.svg)](https://github.com/gridap/GridapP4est.jl/actions) [![Coverage](https://codecov.io/gh/gridap/GridapP4est.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/gridap/GridapP4est.jl) -| ![](https://user-images.githubusercontent.com/38347633/134634010-2be9b499-201b-4166-80ac-e161f6adceb0.png) | ![](https://user-images.githubusercontent.com/38347633/134634023-83f37646-f6b9-435c-9f9f-291dea9f86c2.png) -|:-------------:|:-------------:| +![amr_cubed_sphere](https://github.com/gridap/GridapP4est.jl/assets/38347633/596ab00e-58a8-4aeb-bb0f-efbe63cb2b59) + ## Purpose 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 7fda1fe..93b36c2 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -1,3 +1,7 @@ + +const LagragianOrNedelec = + Union{Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any},Tuple{Gridap.ReferenceFEs.Nedelec,Any,Any}} + function _generate_ref_constraints(modelH,modelh,cell_reffe) VH=TestFESpace(modelH,cell_reffe) Vh=TestFESpace(modelh,cell_reffe) @@ -31,7 +35,7 @@ end function _build_constraint_coefficients_matrix_in_ref_space(::PXestUniformRefinementRuleType, Dc, - reffe::Tuple{<:Lagrangian,Any,Any}) + reffe::LagragianOrNedelec) cell_polytope = Dc == 2 ? QUAD : HEX basis, reffe_args, reffe_kwargs = reffe cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) @@ -94,7 +98,7 @@ function _fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof, end function _generate_face_subface_ldof_to_cell_ldof(ref_rule::PXestUniformRefinementRuleType, - Df,Dc,reffe::Tuple{<:Lagrangian,Any,Any}) + Df,Dc,reffe::LagragianOrNedelec) cell_polytope = (Dc == 2) ? QUAD : HEX _coarse_faces_to_child_ids = coarse_faces_to_child_ids(ref_rule,Df,Dc) @@ -276,9 +280,6 @@ function coarse_faces_to_child_ids(::PXestVerticalRefinementRuleType,Df,Dc) end end - - - function coarse_faces_to_child_ids(::PXestHorizontalRefinementRuleType,Df,Dc) @assert Dc==3 if (Df==2) @@ -351,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, @@ -450,7 +420,10 @@ function _restrict_face_dofs_to_face_dim(cell_reffe,Df) first_face = Gridap.ReferenceFEs.get_offset(polytope,Df)+1 face_own_dofs=get_face_own_dofs(cell_reffe) first_face_faces = Gridap.ReferenceFEs.get_faces(polytope)[first_face] - face_dofs_to_face_dim = face_own_dofs[first_face_faces] + face_dofs_to_face_dim = Vector{Vector{Int}}() + for face in first_face_faces + push!(face_dofs_to_face_dim,copy(face_own_dofs[face])) + end touched=Dict{Int,Int}() dofs=Int64[] current=1 @@ -599,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}) @@ -815,7 +616,22 @@ 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 + order = length(get_face_dofs(reffe)[first_face+1])-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 get_face_dofs_permutations( + reffe::Gridap.ReferenceFEs.GenericRefFE{Gridap.ReferenceFEs.Nedelec, Dc},Df::Integer) where Dc + first_face = get_offset(get_polytope(reffe),Df) + order = length(get_face_dofs(reffe)[first_face+1])-1 nfaces=num_faces(reffe,Df) if (Df==Dc-1) facet_polytope = Dc == 2 ? SEGMENT : QUAD @@ -849,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...) @@ -916,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] = @@ -1142,7 +916,7 @@ end # Generates a new DistributedSingleFieldFESpace composed # by local FE spaces with linear multipoint constraints added function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc}, - reffe::Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any}; + reffe::LagragianOrNedelec; kwargs...) where {Dc} spaces_wo_constraints = map(local_views(model)) do m FESpace(m, reffe; kwargs...) 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 144a76a..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, @@ -1599,9 +1838,23 @@ end # Assumptions. Either: # A) model.parts MPI tasks are included in parts_redistributed_model MPI tasks; or # B) model.parts MPI tasks include parts_redistributed_model MPI tasks +const WeightsArrayType=Union{Nothing,MPIArray{<:Vector{<:Integer}}} function GridapDistributed.redistribute(model::OctreeDistributedDiscreteModel{Dc,Dp}, - parts_redistributed_model=model.parts) where {Dc,Dp} + parts_redistributed_model=model.parts; + weights::WeightsArrayType=nothing) where {Dc,Dp} parts = (parts_redistributed_model === model.parts) ? model.parts : parts_redistributed_model + _weights=nothing + if (weights !== nothing) + Gridap.Helpers.@notimplementedif parts!==model.parts + _weights=map(model.dmodel.models,weights) do lmodel,weights + # The length of the local weights array has to match the number of + # cells in the model. This includes both owned and ghost cells. + # Only the flags for owned cells are actually taken into account. + @assert num_cells(lmodel)==length(weights) + convert(Vector{Cint},weights) + end + end + comm = parts.comm if (GridapDistributed.i_am_in(model.parts.comm) || GridapDistributed.i_am_in(parts.comm)) if (parts_redistributed_model !== model.parts) @@ -1610,7 +1863,7 @@ function GridapDistributed.redistribute(model::OctreeDistributedDiscreteModel{Dc @assert A || B end if (parts_redistributed_model===model.parts || A) - _redistribute_parts_subseteq_parts_redistributed(model,parts_redistributed_model) + _redistribute_parts_subseteq_parts_redistributed(model,parts_redistributed_model,_weights) else _redistribute_parts_supset_parts_redistributed(model, parts_redistributed_model) end @@ -1619,7 +1872,9 @@ function GridapDistributed.redistribute(model::OctreeDistributedDiscreteModel{Dc end end -function _redistribute_parts_subseteq_parts_redistributed(model::OctreeDistributedDiscreteModel{Dc,Dp}, parts_redistributed_model) where {Dc,Dp} +function _redistribute_parts_subseteq_parts_redistributed(model::OctreeDistributedDiscreteModel{Dc,Dp}, + parts_redistributed_model, + _weights::WeightsArrayType) where {Dc,Dp} parts = (parts_redistributed_model === model.parts) ? model.parts : parts_redistributed_model if (parts_redistributed_model === model.parts) ptr_pXest_old = model.ptr_pXest @@ -1631,7 +1886,15 @@ function _redistribute_parts_subseteq_parts_redistributed(model::OctreeDistribut parts.comm) end ptr_pXest_new = pXest_copy(model.pXest_type, ptr_pXest_old) - pXest_partition!(model.pXest_type, ptr_pXest_new) + if (_weights !== nothing) + init_fn_callback_c = pXest_reset_callbacks(model.pXest_type) + map(_weights) do _weights + pXest_reset_data!(model.pXest_type, ptr_pXest_new, Cint(sizeof(Cint)), init_fn_callback_c, pointer(_weights)) + end + pXest_partition!(model.pXest_type, ptr_pXest_new; weights_set=true) + else + pXest_partition!(model.pXest_type, ptr_pXest_new; weights_set=false) + end # Compute RedistributeGlue parts_snd, lids_snd, old2new = pXest_compute_migration_control_data(model.pXest_type,ptr_pXest_old,ptr_pXest_new) @@ -1648,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, @@ -1880,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, @@ -1893,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 @@ -1980,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 6397742..4317dfd 100644 --- a/src/PXestTypeMethods.jl +++ b/src/PXestTypeMethods.jl @@ -309,16 +309,31 @@ function pXest_balance!(::P8estType, ptr_pXest; k_2_1_balance=0) end end -function pXest_partition!(::P4estType, ptr_pXest) - p4est_partition(ptr_pXest, 0, C_NULL) +function pXest_partition!(pXest_type::P4estType, ptr_pXest; weights_set=false) + if (!weights_set) + p4est_partition(ptr_pXest, 0, C_NULL) + else + wcallback=pXest_weight_callback(pXest_type) + p4est_partition(ptr_pXest, 0, wcallback) + end end -function pXest_partition!(::P6estType, ptr_pXest) - p6est_partition(ptr_pXest, C_NULL) +function pXest_partition!(pXest_type::P6estType, ptr_pXest; weights_set=false) + if (!weights_set) + p6est_partition(ptr_pXest, C_NULL) + else + wcallback=pXest_weight_callback(pXest_type) + p6est_partition(ptr_pXest, wcallback) + end end -function pXest_partition!(::P8estType, ptr_pXest) - p8est_partition(ptr_pXest, 0, C_NULL) +function pXest_partition!(pXest_type::P8estType, ptr_pXest; weights_set=false) + if (!weights_set) + p8est_partition(ptr_pXest, 0, C_NULL) + else + wcallback=pXest_weight_callback(pXest_type) + p8est_partition(ptr_pXest, 0, wcallback) + end end @@ -805,6 +820,30 @@ function pXest_refine_callbacks(::P8estType) refine_callback_c, refine_replace_callback_c end +function pXest_weight_callback(::P4estType) + function weight_callback(::Ptr{p4est_t}, + which_tree::p4est_topidx_t, + quadrant_ptr::Ptr{p4est_quadrant_t}) + quadrant = quadrant_ptr[] + return unsafe_wrap(Array, Ptr{Cint}(quadrant.p.user_data), 1)[] + end + @cfunction($weight_callback, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t})) +end + +function pXest_weight_callback(::P6estType) + Gridap.Helpers.@notimplemented +end + +function pXest_weight_callback(::P8estType) + function weight_callback(::Ptr{p8est_t}, + which_tree::p4est_topidx_t, + quadrant_ptr::Ptr{p8est_quadrant_t}) + quadrant = quadrant_ptr[] + return unsafe_wrap(Array, Ptr{Cint}(quadrant.p.user_data), 1)[] + end + @cfunction($weight_callback, Cint, (Ptr{p8est_t}, p4est_topidx_t, Ptr{p8est_quadrant_t})) +end + function _unwrap_ghost_quadrants(::P4estType, pXest_ghost) Ptr{p4est_quadrant_t}(pXest_ghost.ghosts.array) end @@ -1022,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 @@ -1216,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) @@ -1278,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, @@ -1835,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 @@ -1850,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/DarcyNonConformingOctreeModelsTests.jl b/test/DarcyNonConformingOctreeModelsTests.jl index 2d513cf..d1e909c 100644 --- a/test/DarcyNonConformingOctreeModelsTests.jl +++ b/test/DarcyNonConformingOctreeModelsTests.jl @@ -142,7 +142,7 @@ module DarcyNonConformingOctreeModelsTests dΩh = Measure(trian,degree) el2 = sqrt(sum( ∫( e⋅e )*dΩh )) - tol=1e-7 + tol=1e-6 @assert el2 < tol end @@ -247,25 +247,11 @@ module DarcyNonConformingOctreeModelsTests ep_l2 = l2(ep) ep_h1 = h1(ep) - tol = 1.0e-9 + tol = 1.0e-6 @test eu_l2 < tol @test ep_l2 < tol @test ep_h1 < tol - end - - function driver(ranks,coarse_model,order) - model=OctreeDistributedDiscreteModel(ranks,coarse_model,1) - ref_coarse_flags=map(ranks,partition(get_cell_gids(model.dmodel))) do rank,indices - flags=zeros(Cint,length(indices)) - flags.=nothing_flag - #flags[1:2^2].=coarsen_flag - flags[own_length(indices)]=refine_flag - flags - end - fmodel,glue=Gridap.Adaptivity.adapt(model,ref_coarse_flags) - xh,Xh = solve_darcy(fmodel,order) - check_error_darcy(fmodel,order,xh) - end + end function run(distribute) # debug_logger = ConsoleLogger(stderr, Logging.Debug) @@ -276,7 +262,7 @@ module DarcyNonConformingOctreeModelsTests order=1 test_2d(ranks,order) - order=1 + order=2 test_3d(ranks,order) end end \ No newline at end of file diff --git a/test/MaxwellNonConformingOctreeModelsTests.jl b/test/MaxwellNonConformingOctreeModelsTests.jl new file mode 100644 index 0000000..3abc6f8 --- /dev/null +++ b/test/MaxwellNonConformingOctreeModelsTests.jl @@ -0,0 +1,241 @@ +module MaxwellNonConformingOctreeModelsTests + using P4est_wrapper + using GridapP4est + using Gridap + using PartitionedArrays + using GridapDistributed + using MPI + using Gridap.FESpaces + using FillArrays + using Logging + using Test + + function test_transfer_ops_and_redistribute(ranks, + dmodel::GridapDistributed.DistributedDiscreteModel{Dc}, + order) where Dc + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices + flags=zeros(Cint,length(indices)) + flags.=nothing_flag + + flags[1]=refine_flag + flags[own_length(indices)]=refine_flag + + # To create some unbalance + if (rank%2==0 && own_length(indices)>1) + flags[div(own_length(indices),2)]=refine_flag + end + flags + end + fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); + + # Solve coarse + uH,UH=solve_maxwell(dmodel,order) + check_error_maxwell(dmodel,order,uH) + + # Solve fine + uh,Uh=solve_maxwell(fmodel,order) + check_error_maxwell(fmodel,order,uh) + + Ωh = Triangulation(fmodel) + degree = 2*(order+1) + dΩh = Measure(Ωh,degree) + + # prolongation via interpolation + uHh=interpolate(uH,Uh) + e = uh - uHh + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + tol=1e-6 + println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + # prolongation via L2-projection + # Coarse FEFunction -> Fine FEFunction, by projection + ahp(u,v) = ∫(v⋅u)*dΩh + lhp(v) = ∫(v⋅uH)*dΩh + oph = AffineFEOperator(ahp,lhp,Uh,Uh) + uHh = solve(oph) + e = uh - uHh + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + println("[L2 PROJECTION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + # restriction via interpolation + uhH=interpolate(uh,UH) + e = uH - uhH + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + # restriction via L2-projection + ΩH = Triangulation(dmodel) + degree = 2*(order+1) + dΩH = Measure(ΩH,degree) + + dΩhH = Measure(ΩH,Ωh,2*order) + aHp(u,v) = ∫(v⋅u)*dΩH + lHp(v) = ∫(v⋅uh)*dΩhH + oph = AffineFEOperator(aHp,lHp,UH,UH) + uhH = solve(oph) + e = uH - uhH + el2 = sqrt(sum( ∫( e⋅e )*dΩH )) + + fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); + uh_red,Uh_red=solve_maxwell(fmodel_red,order) + check_error_maxwell(fmodel_red,order,uh_red) + + trian = Triangulation(fmodel_red) + degree = 2*(order+1) + dΩhred = Measure(trian,degree) + + u_ex, f_ex=get_analytical_functions(Dc) + + uhred2 = GridapDistributed.redistribute_fe_function(uh,Uh_red,fmodel_red,red_glue) + + e = u_ex - uhred2 + el2 = sqrt(sum( ∫( e⋅e )*dΩhred )) + println("[REDISTRIBUTE SOLUTION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + fmodel_red + end + + function test_refine_and_coarsen_at_once(ranks, + dmodel::OctreeDistributedDiscreteModel{Dc}, + order) where Dc + degree = 2*order+1 + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices + flags=zeros(Cint,length(indices)) + flags.=nothing_flag + if (rank==1) + flags[1:min(2^Dc,own_length(indices))].=coarsen_flag + end + flags[own_length(indices)]=refine_flag + flags + end + fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); + + # Solve coarse + uH,UH=solve_maxwell(dmodel,order) + check_error_maxwell(dmodel,order,uH) + + # Solve fine + uh,Uh=solve_maxwell(fmodel,order) + check_error_maxwell(fmodel,order,uh) + + # # prolongation via interpolation + uHh=interpolate(uH,Uh) + e = uh - uHh + + trian = Triangulation(fmodel) + degree = 2*(order+1) + dΩh = Measure(trian,degree) + + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + tol=1e-6 + @assert el2 < tol + end + + function test_2d(ranks,order) + coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) + dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,1) + test_refine_and_coarsen_at_once(ranks,dmodel,order) + rdmodel=dmodel + for i=1:5 + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) + end + end + + function test_3d(ranks,order) + coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(2,2,2)) + dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,0) + test_refine_and_coarsen_at_once(ranks,dmodel,order) + rdmodel=dmodel + for i=1:5 + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) + end + end + + u_ex_2D((x,y)) = 2*VectorValue(-y,x) + f_ex_2D(x) = u_ex_2D(x) + u_ex_3D((x,y,z)) = 2*VectorValue(-y,x,0.) - VectorValue(0.,-z,y) + f_ex_3D(x) = u_ex_3D(x) + + function get_analytical_functions(Dc) + if Dc==2 + return u_ex_2D, f_ex_2D + else + @assert Dc==3 + return u_ex_3D, f_ex_3D + end + end + + include("CoarseDiscreteModelsTools.jl") + + function solve_maxwell(model::GridapDistributed.DistributedDiscreteModel{Dc},order) where {Dc} + u_ex, f_ex=get_analytical_functions(Dc) + + V = FESpace(model, + ReferenceFE(nedelec,order), + conformity=:Hcurl, + dirichlet_tags="boundary") + + U = TrialFESpace(V,u_ex) + + trian = Triangulation(model) + degree = 2*(order+1) + dΩ = Measure(trian,degree) + + a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ + l(v) = ∫(f_ex⋅v)dΩ + + op = AffineFEOperator(a,l,U,V) + if (num_free_dofs(U)==0) + # UMFPACK cannot handle empty linear systems + uh = zero(U) + else + uh = solve(op) + end + + # uh_ex=interpolate(u_ex_3D,U) + # map(local_views(get_free_dof_values(uh_ex)), local_views(op.op.matrix), local_views(op.op.vector)) do U_ex, A, b + # r_ex = A*U_ex - b + # println(r_ex) + # end + uh,U + end + + function check_error_maxwell(model::GridapDistributed.DistributedDiscreteModel{Dc},order,uh) where {Dc} + trian = Triangulation(model) + degree = 2*(order+1) + dΩ = Measure(trian,degree) + + u_ex, f_ex = get_analytical_functions(Dc) + + eu = u_ex - uh + + l2(v) = sqrt(sum(∫(v⋅v)*dΩ)) + hcurl(v) = sqrt(sum(∫(v⋅v + (∇×v)⋅(∇×v))*dΩ)) + + eu_l2 = l2(eu) + eu_hcurl = hcurl(eu) + + tol = 1.0e-6 + @test eu_l2 < tol + @test eu_hcurl < tol + end + + function run(distribute) + # debug_logger = ConsoleLogger(stderr, Logging.Debug) + # global_logger(debug_logger); # Enable the debug logger globally + np = MPI.Comm_size(MPI.COMM_WORLD) + ranks = distribute(LinearIndices((np,))) + + for order=0:2 + test_2d(ranks,order) + end + + for order=0:2 + test_3d(ranks,order) + end + end +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 7ec3cb8..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 @@ -139,19 +206,58 @@ module PoissonNonConformingOctreeModelsTests e = uH - uhH el2 = sqrt(sum( ∫( e⋅e )*dΩH )) + weights=map(ranks,fmodel.dmodel.models) do rank,lmodel + if (rank%2==0) + zeros(Cint,num_cells(lmodel)) + else + ones(Cint,num_cells(lmodel)) + end + end fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); - Vhred=FESpace(fmodel_red,reffe,conformity=:H1;dirichlet_tags="boundary") - Uhred=TrialFESpace(Vhred,u) + 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 @@ -169,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) @@ -187,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 @@ -214,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 @@ -235,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 @@ -274,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 - diff --git a/test/mpi/MaxwellNonConformingOctreeModelsTests.jl b/test/mpi/MaxwellNonConformingOctreeModelsTests.jl new file mode 100644 index 0000000..69d3f81 --- /dev/null +++ b/test/mpi/MaxwellNonConformingOctreeModelsTests.jl @@ -0,0 +1,17 @@ + +using MPI +using PartitionedArrays + +include("../MaxwellNonConformingOctreeModelsTests.jl") +import .MaxwellNonConformingOctreeModelsTests as TestModule + +if !MPI.Initialized() + MPI.Init() +end + +with_mpi() do distribute + TestModule.run(distribute) +end + +MPI.Finalize() + diff --git a/test/runtests.jl b/test/runtests.jl index 28168c1..bb038db 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,9 @@ function run_tests(testdir) elseif f in ["DarcyNonConformingOctreeModelsTests.jl"] np = [1,4] extra_args = "" + elseif f in ["MaxwellNonConformingOctreeModelsTests.jl"] + np = [1,4] + extra_args = "" elseif f in ["PoissonNonConformingOctreeModelsTests.jl"] np = [1,2,4] extra_args = ""