diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl index 7344e9a46..67479e481 100644 --- a/src/Adaptivity/EdgeBasedRefinement.jl +++ b/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 diff --git a/src/Geometry/GridTopologies.jl b/src/Geometry/GridTopologies.jl index 7d00b7a15..72a92b95b 100644 --- a/src/Geometry/GridTopologies.jl +++ b/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 """ diff --git a/test/AdaptivityTests/EdgeBasedRefinementTests.jl b/test/AdaptivityTests/EdgeBasedRefinementTests.jl index 45997bfa5..a16b6c47a 100644 --- a/test/AdaptivityTests/EdgeBasedRefinementTests.jl +++ b/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 \ No newline at end of file