Merge pull request #43 from gridap/constraints_rt_fes

Support for RT FEs + non-conforming meshes

Project.toml

@@ -16,9 +16,9 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 ArgParse = "1"
-FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12"
-Gridap = "0.17.18"
-GridapDistributed = "0.3"
+FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 1"
+Gridap = "0.17.20"
+GridapDistributed = "0.3.1"
 MPI = "0.20"
 P4est_wrapper = "0.2.0"
 PartitionedArrays = "0.3.3"

README.md

@@ -10,7 +10,7 @@
 
 ## Purpose 
 
-The purpose of this package is to provide a `DistributedDiscreteModel` implementation (a parallel mesh data structure, see `GridapDistributed.jl` for more details) able to handle forests of quadtrees/octrees of the computational domain. To this end, it leverages the [`p4est` software library](https://p4est.github.io/) meshing engine under the hood (click [here](https://github.com/gridap/GridapP4est.jl/blob/main/test/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl) for an example).
+The purpose of this package is to provide a `DistributedDiscreteModel` implementation (a parallel mesh data structure, see `GridapDistributed.jl` for more details) able to handle forests of quadtrees/octrees of the computational domain. To this end, it leverages the [`p4est` software library](https://p4est.github.io/) meshing engine under the hood (click [here](https://github.com/gridap/GridapP4est.jl/blob/main/test/PoissonUniformlyRefinedOctreeModelsTests.jl) for an example).
 
 ## Build 
 

src/GridapFixes.jl

@@ -4,4 +4,22 @@ function Base.map(::typeof(Gridap.Arrays.testitem),
   a1=Vector{eltype(eltype(a[1]))}(undef,size(a2,1))
   a1.=zero(Gridap.Arrays.testitem(a1))
   (a1,a2)
+end
+
+# Required to transfer fine-grid VECTOR-VALUED fields into coarse-grid 
+function Gridap.Adaptivity.FineToCoarseField(fine_fields::AbstractArray{<:Gridap.Fields.Field},
+                                              rrule::Gridap.Adaptivity.RefinementRule,
+                                              child_ids::AbstractArray{<:Integer})
+
+  grid=Gridap.Adaptivity.get_ref_grid(rrule)
+  D=num_cell_dims(grid)
+  x=zero(Point{D,Float64})
+  ffx=lazy_map(evaluate,fine_fields,Fill([x],length(fine_fields)))
+  ffx=ffx[1]
+  fields = Vector{Gridap.Fields.Field}(undef,Gridap.Adaptivity.num_subcells(rrule))
+  fields = fill!(fields,Gridap.Fields.ConstantField(zero(eltype(ffx))))
+  for (k,id) in enumerate(child_ids)
+    fields[id] = fine_fields[k]
+  end
+  return Gridap.Adaptivity.FineToCoarseField(fields,rrule)
 end

test/CoarseDiscreteModelsTools.jl

@@ -0,0 +1,120 @@
+  function setup_model(::Type{Val{3}}, perm)
+  #               5 +--------+ 7 
+  #                /        /|
+  #               /        / |
+  #            6 +--------+  |
+  #              |        |  |
+  #              |  1     |  + 3 
+  #              |        | /
+  #              |        |/
+  #            2 +--------+ 4
+
+  #     6  +--------+ 8 
+  #       /        /|
+  #      /        / |
+  #  11 +--------+  |
+  #     |        |  |
+  #     |  2     |  + 4
+  #     |        | /
+  #     |        |/
+  #   9 +--------+ 10
+    ptr  = [ 1, 9, 17 ]
+    if (perm==1)
+      data = [ 1,2,3,4,5,6,7,8, 2,9,4,10,6,11,8,12 ]
+    elseif (perm==2)
+      data = [ 1,2,3,4,5,6,7,8, 10,12,4,8,9,11,2,6 ]
+    elseif (perm==3)
+      data = [ 1,2,3,4,5,6,7,8, 12,11,8,6,10,9,4,2 ]
+    elseif (perm==4) 
+      data = [ 1,2,3,4,5,6,7,8, 11,9,6,2,12,10,8,4 ]
+    end  
+    cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+    node_coordinates = Vector{Point{3,Float64}}(undef,12)
+    node_coordinates[1]=Point{3,Float64}(0.0,0.0,0.0)
+    node_coordinates[2]=Point{3,Float64}(1.0,0.0,0.0)
+    node_coordinates[3]=Point{3,Float64}(0.0,1.0,0.0)
+    node_coordinates[4]=Point{3,Float64}(1.0,1.0,0.0)
+    node_coordinates[5]=Point{3,Float64}(0.0,0.0,1.0)
+    node_coordinates[6]=Point{3,Float64}(1.0,0.0,1.0)
+    node_coordinates[7]=Point{3,Float64}(0.0,1.0,1.0)
+    node_coordinates[8]=Point{3,Float64}(1.0,1.0,1.0)
+    node_coordinates[9]=Point{3,Float64}(2.0,0.0,0.0)
+    node_coordinates[10]=Point{3,Float64}(2.0,1.0,0.0)
+    node_coordinates[11]=Point{3,Float64}(2.0,0.0,1.0)
+    node_coordinates[12]=Point{3,Float64}(2.0,1.0,1.0)
+
+    polytope=HEX
+    scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
+    cell_types=collect(Fill(1,length(cell_vertex_lids)))
+    cell_reffes=[scalar_reffe]
+    grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
+                                            cell_vertex_lids,
+                                            cell_reffes,
+                                            cell_types,
+                                            Gridap.Geometry.NonOriented())
+    m=Gridap.Geometry.UnstructuredDiscreteModel(grid)
+    labels = get_face_labeling(m)
+    labels.d_to_dface_to_entity[1].=2
+    if (perm==1 || perm==2)
+      labels.d_to_dface_to_entity[2].=2
+      labels.d_to_dface_to_entity[3].=[2,2,2,2,2,1,2,2,2,2,2]
+    elseif (perm==3 || perm==4)
+      labels.d_to_dface_to_entity[2].=2
+      labels.d_to_dface_to_entity[3].=[2,2,2,2,2,1,2,2,2,2,2]
+    end
+    labels.d_to_dface_to_entity[4].=1  
+    add_tag!(labels,"boundary",[2])
+    add_tag!(labels,"interior",[1])
+    m
+  end
+
+  function setup_model(::Type{Val{2}}, perm)
+    @assert perm ∈ (1,2,3,4)
+    #
+    #  3-------4-------6
+    #  |       |       |
+    #  |       |       |
+    #  |       |       |
+    #  1-------2-------5
+    #
+        ptr  = [ 1, 5, 9 ]
+        if (perm==1)
+          data = [ 1,2,3,4, 2,5,4,6 ]
+        elseif (perm==2)
+          data = [ 1,2,3,4, 6,4,5,2 ]
+        elseif (perm==3)
+          data = [ 4,3,2,1, 2,5,4,6 ]
+        elseif (perm==4) 
+          data = [ 4,3,2,1, 6,4,5,2 ]
+        end  
+        cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+        node_coordinates = Vector{Point{2,Float64}}(undef,6)
+        node_coordinates[1]=Point{2,Float64}(0.0,0.0)
+        node_coordinates[2]=Point{2,Float64}(1.0,0.0)
+        node_coordinates[3]=Point{2,Float64}(0.0,1.0)
+        node_coordinates[4]=Point{2,Float64}(1.0,1.0)
+        node_coordinates[5]=Point{2,Float64}(2.0,0.0)
+        node_coordinates[6]=Point{2,Float64}(2.0,1.0)
+
+        polytope=QUAD
+        scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
+        cell_types=collect(Fill(1,length(cell_vertex_lids)))
+        cell_reffes=[scalar_reffe]
+        grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
+                                                cell_vertex_lids,
+                                                cell_reffes,
+                                                cell_types,
+                                                Gridap.Geometry.NonOriented())
+        m=Gridap.Geometry.UnstructuredDiscreteModel(grid)
+        labels = get_face_labeling(m)
+        labels.d_to_dface_to_entity[1].=2
+        if (perm==1 || perm==2)
+          labels.d_to_dface_to_entity[2].=[2,2,2,1,2,2,2]
+        elseif (perm==3 || perm==4)
+          labels.d_to_dface_to_entity[2].=[2,2,1,2,2,2,2] 
+        end
+        labels.d_to_dface_to_entity[3].=1
+        add_tag!(labels,"boundary",[2])
+        add_tag!(labels,"interior",[1])
+        m
+  end