EdgeBasedRefinement now yields Oriented simplex meshes

src/Adaptivity/EdgeBasedRefinement.jl

@@ -1,3 +1,33 @@
+"""
+Note on RefinementRules and Orientation of the refined grids: 
+
+  In order to guarantee that the refined grid is Oriented (which is something
+  we want for div- and curl-conforming discretisations), we need to guarantee 
+  for simplices that each fine cell has it's vertex gids sorted in increasing 
+  order.
+  
+  In the case of refined meshes, this has an additional constraint: the sorting 
+  of the gids CANNOT be done after the refinement. This would make the glue 
+  inconsistent. This is an attempt to explain why: 
+  If we change the order of the gids of a fine cell, we are basically applying a 
+  rotation(+symmetry) to the reference space of that cell. The mesh itself will be 
+  aware of this rotation, but the RefinementRule will not. So the reference space 
+  corresponding to that fine cell within the RefinementRule will NOT be rotated 
+  accordingly. This creates an inconsistency, and the fine-to-coarse/coarse-to-fine
+  mesh transfers will therefore be broken (the glue is no longer valid). 
+
+  How to fix this in the future? The glue should probably keep a record of the 
+  cell-wise rotations. This could be an auxiliary cell_map that maps the reference 
+  space of the fine cell within the fine mesh to the reference space of the fine cell 
+  within the refinement rule.
+  For instance: 
+    - Φ: Cell_map given by the RefinementRule, going from the fine reference space
+         to the coarse reference space.
+    - β: Cell_map encoding the rotation from the reference space of the fine mesh to
+         the reference space of the RefinementRule.
+  then we would have 
+      X_coarse = Φ∘β(X_fine)
+"""
 
 struct EdgeBasedRefinement <: AdaptivityMethod end
 
@@ -30,7 +60,7 @@ function refine_unstructured_topology(topo::UnstructuredGridTopology{Dc},
   c2n_map_new = get_refined_cell_to_vertex_map(topo,rrules,faces_list)
   polys_new, cell_type_new = get_refined_polytopes(rrules)
 
-  return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,topo.orientation_style)
+  return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,OrientationStyle(topo))
 end
 
 function get_refined_polytopes(rrules::AbstractArray{<:RefinementRule})
@@ -228,6 +258,7 @@ function get_refined_cell_to_vertex_map(topo::UnstructuredGridTopology{2},
   return Table(data_new,ptrs_new)
 end
 
+###############################################################################
 # EdgeBasedRefinementRules
 
 abstract type EdgeBasedRefinementRule <: RefinementRuleType end
@@ -271,9 +302,9 @@ function _get_red_refined_connectivity(p::Polytope{2})
     polys     = [TRI]
     cell_type = [1,1,1,1]
     conn_data = [1,4,5,
-                4,2,6,
-                5,6,3,
-                6,5,4]
+                2,4,6,
+                3,5,6,
+                4,5,6]
     conn_ptrs = [1,4,7,10,13]
     return polys, cell_type, Table(conn_data,conn_ptrs)
   elseif p == QUAD
@@ -315,20 +346,23 @@ function GreenRefinementRule(p::Polytope{2},ref_edge::Integer)
 end
 
 function _get_green_refined_connectivity(p::Polytope{2},ref_edge)
+  # Note: Sorting is necessary in order to guarantee that the gids
+  #       of the refined mesh are sorted (and therefore that the fine 
+  #       grid is Oriented). See the note at top of the file. 
   P = _get_green_vertex_permutation(p,ref_edge)
   if p == TRI
     polys     = [TRI]
     cell_type = [1,1]
-    conn_data = [4,P[1],P[2],
-                 4,P[3],P[1]]
+    conn_data = [sort([4,P[1],P[2]])...,
+                 sort([4,P[3],P[1]])...]
     conn_ptrs = [1,4,7]
     return polys, cell_type, Table(conn_data,conn_ptrs)
   elseif p == QUAD
     polys     = [TRI]
     cell_type = [1,1,1]
-    conn_data = [P[1],P[2],5,
-                 P[3],P[1],5,
-                 P[4],5,P[2]]
+    conn_data = [sort([P[1],P[2],5])...,
+                 sort([P[3],P[1],5])...,
+                 sort([P[4],5,P[2]])...]
     conn_ptrs = [1,4,7,10]
     return polys, cell_type, Table(conn_data,conn_ptrs)
   end

src/Geometry/GridTopologies.jl

@@ -107,6 +107,25 @@ function test_grid_topology(top::GridTopology{Dc,Dp}) where {Dc,Dp}
   end
 end
 
+"""
+  If OrientationStyle(topo)==true, checks if the topology is indeed oriented.
+  If OrientationStyle(topo)==false or the topology is not oriented, returns false.
+"""
+function test_grid_topology_orientation(topo::GridTopology{Dc,Dp}) where {Dc,Dp}
+  orientation = OrientationStyle(topo)
+  isa(orientation,NonOriented) && (return false)
+
+  nC = num_faces(topo,Dc)
+  c2n_map     = get_faces(topo,Dc,0)
+  is_oriented = true; iC = 1
+  while(is_oriented && iC <= nC)
+    nodes = c2n_map[iC]
+    is_oriented = (is_oriented && issorted(nodes))
+    iC += 1
+  end
+  return is_oriented
+end
+
 # Default API
 
 """

test/AdaptivityTests/EdgeBasedRefinementTests.jl

@@ -12,29 +12,29 @@ visualize = false
 
 # Refining meshes of QUADs
 cart_model = CartesianDiscreteModel((0,1,0,1),(4,4))
-model      = UnstructuredDiscreteModel(cart_model)
+model1     = UnstructuredDiscreteModel(cart_model)
 
 ## Homogeneous refinement
-ref_model1 = refine(model)
+ref_model1 = refine(model1)
 trian1 = Triangulation(ref_model1.model)
 visualize && writevtk(trian1,"test/AdaptivityTests/ref_model1")
 
 ## Propagate to all-red
-ref_model2 = refine(model;cells_to_refine=[1,6,11,16])
+ref_model2 = refine(model1;cells_to_refine=[1,6,11,16])
 trian2 = Triangulation(ref_model2.model)
 visualize && writevtk(trian2,"test/AdaptivityTests/ref_model2")
 
 ## Red-Green refinement
-ref_model3 = refine(model;cells_to_refine=[1,6,16])
+ref_model3 = refine(model1;cells_to_refine=[1,6,16])
 trian3 = Triangulation(ref_model3.model)
 visualize && writevtk(trian3,"test/AdaptivityTests/ref_model3")
 
-ref_model4 = refine(model;cells_to_refine=[6,7,10,11])
+ref_model4 = refine(model1;cells_to_refine=[6,7,10,11])
 trian4 = Triangulation(ref_model4.model)
 visualize && writevtk(trian4,"test/AdaptivityTests/ref_model4")
 
 # Refining meshes of TRIans
-model2 = simplexify(model)
+model2 = simplexify(model1)
 visualize && writevtk(Triangulation(model2),"test/AdaptivityTests/base_model2")
 
 ref_model5 = refine(model2)
@@ -45,10 +45,29 @@ ref_model6 = refine(model2;cells_to_refine=[1,6,16])
 trian6 = Triangulation(ref_model6.model)
 visualize && writevtk(trian6,"test/AdaptivityTests/ref_model6")
 
-# Testing FESpaces
-sol(x) = x[1] + x[2]
-reffe = ReferenceFE(lagrangian,Float64,1)
-V = TestFESpace(ref_model6,reffe;conformity=:H1,dirichlet_tags="boundary")
-U = TrialFESpace(V,sol)
+# Testing FE related functionality
+parent = model2
+model  = ref_model6
+
+topo = get_grid_topology(model)
+@test isa(OrientationStyle(topo),Oriented)
+@test Gridap.Geometry.test_grid_topology_orientation(topo) == true
+
+u((x,y)) = 2*VectorValue(-y,x)
+reffe = ReferenceFE(nedelec,Float64,1)
+Vh = TestFESpace(model,reffe;dirichlet_tags="boundary")
+Uh = TrialFESpace(Vh,u)
+VH = TestFESpace(parent,reffe;dirichlet_tags="boundary")
+UH = TrialFESpace(VH,u)
+
+uh = interpolate(u,Uh)
+uH = interpolate(u,UH)
+
+Ωh = Triangulation(model)
+dΩh = Measure(Ωh,3)
+
+e = u - uH
+el2 = sqrt(sum( ∫( e⋅e )*dΩh ))
+@test el2 < 1.0e-10
 
 end