Support for FESpaces and Facet integration on Triangulations covering a subdomain of that of the Octree model

Project.toml

@@ -18,7 +18,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.19.2"
+Gridap = "0.19.6"
 GridapDistributed = "0.4"
 Logging = "1"
 MPI = "0.20"

src/PXestTypeMethods.jl

@@ -1914,36 +1914,15 @@ function generate_cell_faces_and_non_conforming_glue(pXest_type::PXestType,
     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,
-        n_cell_vertices,
-        n_cell_edges,
-        n_cell_faces,
-        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(
-            n_cell_vertices,
-            n_cell_edges,
-            n_cell_faces,
-            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
+    _fill_owner_faces_pindex_and_lids!(Dc,
+                                       pXest_refinement_rule,
+                                       owner_faces_pindex,
+                                       owner_faces_lids,
+                                       num_hanging_faces,
+                                       hanging_faces_glue,
+                                       hanging_faces_to_cell,
+                                       hanging_faces_to_lface,
+                                       gridap_cell_faces)
 
     return num_regular_faces, 
            num_hanging_faces,

test/DarcyNonConformingOctreeModelsTests.jl

@@ -10,6 +10,10 @@ module DarcyNonConformingOctreeModelsTests
   using Logging
   using Test
 
+  include("CoarseDiscreteModelsTools.jl")
+  include("GenerateTriangulationPortionTools.jl")
+
+
   function test_transfer_ops_and_redistribute(ranks,
                                               dmodel::GridapDistributed.DistributedDiscreteModel{Dc},
                                               order) where Dc
@@ -166,6 +170,87 @@ module DarcyNonConformingOctreeModelsTests
     end
   end
 
+  function test_2d_fe_space_on_triangulation(ranks,order;num_amr_steps=5,num_ghost_layers=1)
+    coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1))
+    dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2;num_ghost_layers=num_ghost_layers)
+    dtrian=generate_triangulation_portion(ranks,dmodel)
+    for i=1:num_amr_steps
+     dmodel,dtrian=test_fe_space_on_triangulation(ranks,dmodel,dtrian,order,i)
+    end
+  end
+
+  function test_fe_space_on_triangulation(ranks,cmodel,ctrian,order,amr_step)
+    ref_coarse_flags=map(ranks,partition(get_cell_gids(cmodel.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(cmodel,ref_coarse_flags)
+    ftrian = generate_triangulation_portion(ranks,fmodel,ctrian=ctrian,glue=glue)
+
+    u_ex, p_ex, f_ex=get_analytical_functions(2)
+
+    V = FESpace(ftrian,
+                ReferenceFE(raviart_thomas,Float64,order),
+                conformity=:Hdiv,
+                dirichlet_tags="boundary")
+
+    Q = FESpace(ftrian,
+                ReferenceFE(lagrangian,Float64,order); 
+                conformity=:L2)
+
+    U = TrialFESpace(V,u_ex)
+    P = TrialFESpace(Q)
+
+    Y = MultiFieldFESpace([V, Q])
+    X = MultiFieldFESpace([U, P])
+
+    degree = 2*(order+1)
+    dΩ = Measure(ftrian,degree)
+
+    btrian = Boundary(ftrian,tags="interior_boundary")
+    degree = 2*(order+1)
+    dΓ = Measure(btrian,degree)
+    nb = get_normal_vector(btrian)
+
+    a((u, p),(v, q)) = ∫(u⋅v)dΩ +∫(q*(∇⋅u))dΩ-∫((∇⋅v)*p)dΩ
+    b(( v, q)) = ∫( v⋅f_ex + q*(∇⋅u_ex))dΩ - ∫((v⋅nb)*p_ex )dΓ
+
+    op = AffineFEOperator(a,b,X,Y)
+    xh = solve(op)
+
+    uh, ph = xh
+    eu = u_ex - uh
+    ep = p_ex - ph
+
+    l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
+
+    eu_l2 = l2(eu)
+    ep_l2 = l2(ep)
+
+    # writevtk(fmodel, "fmodel_amr_level_$(amr_step)")
+    # writevtk(ftrian, "ftrian_amr_level_$(amr_step)")
+
+    tol=1e-5
+    println("[SOLVE FINE] eu_l2 < tol: $(eu_l2) < $(tol)")
+    println("[SOLVE FINE] ep_l2 < tol: $(ep_l2) < $(tol)")
+
+    tol = 1.0e-6
+    @assert eu_l2 < tol
+    @assert ep_l2 < tol
+
+    fmodel, ftrian
+  end
+
+
   u_ex_2D(x) = VectorValue(2*x[1],x[1]+x[2])
   p_ex_2D(x) = x[1]-x[2]
   f_ex_2D(x) = u_ex_2D(x) + ∇(p_ex_2D)(x)
@@ -182,7 +267,6 @@ module DarcyNonConformingOctreeModelsTests
     end
   end
 
-  include("CoarseDiscreteModelsTools.jl")
 
   function solve_darcy(model::GridapDistributed.DistributedDiscreteModel{Dc},order) where {Dc}
     if (Dc==2)
@@ -241,16 +325,13 @@ module DarcyNonConformingOctreeModelsTests
     ep = p_ex - ph
 
     l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
-    h1(v) = sqrt(sum(∫(v*v + ∇(v)⋅∇(v))*dΩ))
 
     eu_l2 = l2(eu)
     ep_l2 = l2(ep)
-    ep_h1 = h1(ep)
 
     tol = 1.0e-6
     @test eu_l2 < tol
     @test ep_l2 < tol
-    @test ep_h1 < tol
   end
 
   function run(distribute)
@@ -264,5 +345,9 @@ module DarcyNonConformingOctreeModelsTests
 
     order=2
     test_3d(ranks,order)
+
+    order=1
+    test_2d_fe_space_on_triangulation(ranks,order;num_amr_steps=5)
+
   end
 end

test/GenerateTriangulationPortionTools.jl

@@ -0,0 +1,71 @@
+function _generate_triangulation_portion(ranks, fmodel)
+        trians = map(ranks, 
+             local_views(fmodel.dmodel), 
+                         partition(get_cell_gids(fmodel))) do rank, lmodel, indices
+            mask = Vector{Bool}(undef,num_cells(lmodel))
+            mask .= false
+            cell_to_part = local_to_owner(indices)
+            graph = GridapDistributed.compute_cell_graph(lmodel)
+            ncells = num_cells(lmodel)
+            icell_to_jcells_ptrs = graph.colptr
+            icell_to_jcells_data = graph.rowval
+            for icell in 1:ncells
+               if cell_to_part[icell] == rank
+                  pini = icell_to_jcells_ptrs[icell]
+                  pend = icell_to_jcells_ptrs[icell+1]-1
+                  for p in pini:pend
+                     jcell = icell_to_jcells_data[p]
+                     if cell_to_part[jcell] != rank
+                        mask[icell] = true
+                     end                 
+                  end
+               end
+            end
+            Triangulation(lmodel, mask)
+        end
+        GridapDistributed.DistributedTriangulation(trians,fmodel)
+  end
+
+  function _generate_triangulation_portion(ranks, fmodel, ctrian, glue)
+        trians = map(ranks,
+                 local_views(fmodel.dmodel),
+                 local_views(ctrian),
+                 glue) do rank, fmodel, ctrian, glue
+            mask = Vector{Bool}(undef,num_cells(fmodel))
+            mask .= false
+
+            cglue = get_glue(ctrian, Val{num_cell_dims(ctrian)}())
+            for ccell in cglue.tface_to_mface
+                fcells = glue.o2n_faces_map[ccell]
+                for fcell in fcells
+                    mask[fcell] = true
+                end
+            end
+            Triangulation(fmodel, mask)
+        end
+        GridapDistributed.DistributedTriangulation(trians,fmodel)
+  end
+
+  function generate_triangulation_portion(ranks,fmodel; ctrian=nothing, glue=nothing)
+    if (length(ranks)==1 && ctrian==nothing)
+       trians = map(ranks, local_views(fmodel.dmodel)) do rank, fmodel
+            mask = Vector{Bool}(undef,num_cells(fmodel))
+            mask .= false
+            for cell=1:Int(round(num_cells(fmodel)*0.25))
+               mask[cell] = true
+            end
+            for cell=Int(round(num_cells(fmodel)*0.75)):num_cells(fmodel)
+               mask[cell] = true
+            end
+            Triangulation(fmodel, mask)
+        end
+        return GridapDistributed.DistributedTriangulation(trians,fmodel)
+    else 
+        if ctrian==nothing
+            @assert glue==nothing
+            return _generate_triangulation_portion(ranks, fmodel)
+        else
+            return _generate_triangulation_portion(ranks, fmodel, ctrian, glue)
+        end
+    end
+  end

test/MaxwellNonConformingOctreeModelsTests.jl

@@ -175,7 +175,7 @@ module MaxwellNonConformingOctreeModelsTests
     u_ex, f_ex=get_analytical_functions(Dc)
 
     V = FESpace(model,
-                ReferenceFE(nedelec,order),
+                ReferenceFE(nedelec,Float64, order),
                 conformity=:Hcurl,
                 dirichlet_tags="boundary")