diff --git a/NEWS.md b/NEWS.md index 1fb74152c..dbbd7eb46 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Jacobi polynomial bases. Since PR [#896](https://github.com/gridap/Gridap.jl/pull/896). +- Replaced newest vertex bisection mesh adaptation in + `src/Geometry/NewestVertexBisection.jl` with appropriate changes to + `src/Adaptivity/EdgeBasedRefinement.jl`. see PR + [#901](https://github.com/gridap/Gridap.jl/pull/901). + ### Fixed - ODE operators cache linear system at initial time or the time stored by the operator. Before, the linear system was cached at time `t = 0.0`, which cannot be done if the operator is not well-defined at `t = 0.0`. Since PR [#891](https://github.com/gridap/Gridap.jl/pull/891). @@ -18,6 +23,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 was using default, it's now using the spetialized implementation for the underlying triangulation type. Since PR [#884](https://github.com/gridap/Gridap.jl/pull/884). +- Fixed `cell_dof_ids` for the case of vectorial `ConstantFESpace`. Since PR [#888](https://github.com/gridap/Gridap.jl/pull/888) + ## [0.17.17] - 2023-02-24 ### Added diff --git a/Project.toml b/Project.toml index 0d6673399..59d16a9d5 100644 --- a/Project.toml +++ b/Project.toml @@ -38,7 +38,7 @@ DataStructures = "0.18.13" DocStringExtensions = "0.8.1, 0.9" FastGaussQuadrature = "0.4.2" FileIO = "1.2.2, 1.3, 1.4" -FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 0.13" +FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 0.13, 1" ForwardDiff = "0.10.10" JLD2 = "0.1.11, 0.3, 0.4" JSON = "0.21.0" diff --git a/src/Adaptivity/AdaptedDiscreteModels.jl b/src/Adaptivity/AdaptedDiscreteModels.jl index 01456c311..cd3c67ea0 100644 --- a/src/Adaptivity/AdaptedDiscreteModels.jl +++ b/src/Adaptivity/AdaptedDiscreteModels.jl @@ -81,8 +81,15 @@ function refine(model::UnstructuredDiscreteModel,::AdaptivityMethod,args...;kwar @abstractmethod end -function refine(model::UnstructuredDiscreteModel,args...;refinement_method=EdgeBasedRefinement(),kwargs...) - return refine(refinement_method,model,args...;kwargs...) +# Handle the user's requested choice for refinement +function string_to_refinement(refinement_method::String, model) + refinement_method == "red_green" && return RedGreenRefinement() + refinement_method == "nvb" && return NVBRefinement(model) + error("refinement_method $refinement_method not recognized") +end + +function refine(model::UnstructuredDiscreteModel,args...;refinement_method="red_green",kwargs...) + return refine(string_to_refinement(refinement_method, model),model,args...;kwargs...) end # CartesianDiscreteModel refining diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl index 22804d44b..d283fb204 100644 --- a/src/Adaptivity/EdgeBasedRefinement.jl +++ b/src/Adaptivity/EdgeBasedRefinement.jl @@ -1,57 +1,106 @@ """ -Note on RefinementRules and Orientation of the refined grids: +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 + 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 + 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). + 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 + 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: + 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 + then we would have X_coarse = Φ∘β(X_fine) """ -struct EdgeBasedRefinement <: AdaptivityMethod end +abstract type EdgeBasedRefinement <: AdaptivityMethod end -function refine(::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,Dp};cells_to_refine=nothing) where {Dc,Dp} - # cells_to_refine can be + +struct NVBRefinement{T} <: EdgeBasedRefinement + cell_to_longest_edge_lid::Vector{T} + cell_to_longest_edge_gid::Vector{T} +end + +_calculate_edge_length(node_coords, nodes) = norm(node_coords[nodes[1]] - node_coords[nodes[2]]) + +# Create Vectors containing the longest edge gids and lids +function _get_longest_edge_ids(c2e_map, e2n_map, node_coords) + nC = length(c2e_map) + e2n_map_cache = array_cache(e2n_map) + c2e_map_cache = array_cache(c2e_map) + longest_edge_gids = Vector{eltype(eltype(c2e_map))}(undef, nC) + longest_edge_lids = Vector{eltype(eltype(c2e_map))}(undef, nC) + # Loop over cells + for c in 1:length(c2e_map) + es = getindex!(c2e_map_cache, c2e_map, c) + e_length_max = 0.0 + e_gid_max = -1 # Longest edge index must be found + e_lid_max = -1 # Longest edge index must be found + for (e_lid, e_gid) in enumerate(es) + ns = getindex!(e2n_map_cache, e2n_map, e_gid) + e_length = _calculate_edge_length(node_coords, ns) + if e_length > e_length_max + e_length_max = e_length + e_gid_max = e_gid + e_lid_max = e_lid + end + end + longest_edge_gids[c] = e_gid_max + longest_edge_lids[c] = e_lid_max + end + return longest_edge_lids, longest_edge_gids +end + +function NVBRefinement(model::DiscreteModel{Dc,Dp}) where {Dc, Dp} + topo = model.grid_topology + @check all(p->p==TRI,get_polytopes(topo)) + c2e_map = get_faces(topo,Dc,1) + e2n_map = get_faces(topo,1 ,0) + node_coords = get_node_coordinates(model) + longest_edge_lids, longest_edge_gids = _get_longest_edge_ids(c2e_map, e2n_map, node_coords) + NVBRefinement(longest_edge_lids, longest_edge_gids) +end + +NVBRefinement(model::AdaptedDiscreteModel) = NVBRefinement(model.model) + +struct RedGreenRefinement <: EdgeBasedRefinement end + +function refine(method::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,Dp};cells_to_refine=nothing) where {Dc,Dp} + # cells_to_refine can be # a) nothing -> All cells get refined - # b) AbstractArray{<:Bool} of size num_cells(model) + # b) AbstractArray{<:Bool} of size num_cells(model) # -> Only cells such that cells_to_refine[iC] == true get refined # c) AbstractArray{<:Integer} # -> Cells for which gid ∈ cells_to_refine get refined ctopo = get_grid_topology(model) coarse_labels = get_face_labeling(model) - # Create new model - rrules, faces_list = setup_edge_based_rrules(ctopo,cells_to_refine) + rrules, faces_list = setup_edge_based_rrules(method, model.grid_topology,cells_to_refine) topo = _refine_unstructured_topology(ctopo,rrules,faces_list) reffes = map(p->LagrangianRefFE(Float64,p,1),get_polytopes(topo)) - grid = UnstructuredGrid(get_vertex_coordinates(topo),get_faces(topo,Dc,0),reffes,get_cell_type(topo),OrientationStyle(topo)) - + grid = UnstructuredGrid(get_vertex_coordinates(topo), + get_faces(topo,Dc,0),reffes, + get_cell_type(topo), + OrientationStyle(topo)) glue = _get_refinement_glue(topo,model.grid_topology,rrules) - labels = _refine_face_labeling(coarse_labels,glue,model.grid_topology,topo) ref_model = UnstructuredDiscreteModel(grid,topo,labels) - return AdaptedDiscreteModel(ref_model,model,glue) end @@ -65,6 +114,7 @@ function _refine_unstructured_topology(topo::UnstructuredGridTopology{Dc}, # We can guarantee the new topology is oriented if # 1 - the old topology was oriented # 2 - we have a single type of polytope (i.e new topo is not mixed) + #orientation = NonOriented() orientation = NonOriented() return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,orientation) @@ -93,27 +143,116 @@ function _get_refinement_glue(ftopo ::UnstructuredGridTopology{Dc}, end """ -Given an UnstructuredTopology and a list of cells_to_refine, provides an edge-based coloring -for the topology: +Given an UnstructuredTopology and a list of cells_to_refine, provides an edge-based coloring +for the topology: - If a cell is completely refined, it is colored RED - If a cell touches more than one refined (RED) cell, it becomes RED. - If a cell touches a single refined (RED) cell, it is colored GREEN. -The method returns a vector of Red/Green refinement rules, as well a the list of +The method returns a vector of Red/Green refinement rules, as well the list of vertices, edges and cells which correspond to new vertices in the refined topology. """ -function setup_edge_based_rrules(topo::UnstructuredGridTopology{Dc},::Nothing) where Dc +function setup_edge_based_rrules(::EdgeBasedRefinement, topo::UnstructuredGridTopology{Dc},::Nothing) where Dc rrules = lazy_map(RedRefinementRule,CompressedArray(topo.polytopes,topo.cell_type)) faces_list = _redgreen_refined_faces_list(topo,rrules,[true]) return rrules, faces_list end -function setup_edge_based_rrules(topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Bool}) where Dc +function setup_edge_based_rrules(::EdgeBasedRefinement,topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Bool}) where Dc return setup_edge_based_rrules(topo,findall(cells_to_refine)) end -function setup_edge_based_rrules(topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc +_shift_to_first(v::AbstractVector{T}, i::T) where {T<:Integer} = circshift(v, -(i - 1)) + +function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc + nE = num_faces(topo,1) + c2e_map = get_faces(topo,Dc,1) + c2e_map_cache = array_cache(c2e_map) + e2c_map = get_faces(topo,1,Dc) + polys = topo.polytopes + cell_types = topo.cell_type + cell_color = copy(cell_types) # WHITE + # Hardcoded for TRI + WHITE, GREEN, BLUE, BLUE_DOUBLE = Int8(1), Int8(2), Int8(5), Int(11) + # Create the inverse mapping from long/short edge pairs to unique indices + BLUE_dict = Dict((1,2) => 1, (1, 3) => 2, (2, 1) => 3, (2, 3) => 4, (3, 1) => 5, (3, 2) => 6) + green_offset = 3 + blue_offset = 3*2 + blue_double_offset = 3 + c_to_longest_edge_gid = method.cell_to_longest_edge_gid + c_to_longest_edge_lid = method.cell_to_longest_edge_lid + is_refined = falses(nE) + # Loop over cells and mark edges to refine i.e. is_refined + # The reason to not loop directly on c is that we need to change c within + # a given iteration of the for loop + for i in 1:length(cells_to_refine) + c = cells_to_refine[i] + e_longest = c_to_longest_edge_gid[c] + # Has to terminate because an edge is marked each iteration or we skip an + # iteration due to a boundary cell + while !is_refined[e_longest] + is_refined[e_longest] = true + c_nbor_lid = findfirst(c′ -> c′ != c, e2c_map[e_longest]) + if isnothing(c_nbor_lid) # We've reach the boundary + continue + else + # Get the longest edge of the neighbor + c_nbor_gid = e2c_map[e_longest][c_nbor_lid] + e_longest = c_to_longest_edge_gid[c_nbor_gid] + # Set the current cell gid to that of the neighbor + c = c_nbor_gid + end + end + end + # Loop over cells and refine based on marked edges + for c in 1:length(c2e_map) + c_edges = getindex!(c2e_map_cache, c2e_map, c) + refined_edge_lids = findall(is_refined[c_edges]) + # GREEN refinement because only one edge should be bisected + if length(refined_edge_lids) == 1 + ref_edge = refined_edge_lids[1] + cell_color[c] = GREEN + Int8(ref_edge-1) + # BLUE refinement: two bisected + elseif length(refined_edge_lids) == 2 + long_ref_edge_lid = c_to_longest_edge_lid[c] + short_ref_edge_lid = setdiff(refined_edge_lids, long_ref_edge_lid)[1] + blue_idx = BLUE_dict[(long_ref_edge_lid, short_ref_edge_lid)] + cell_color[c] = BLUE + Int8(blue_idx - 1) + # DOUBLE BLUE refinement: three bisected edges (somewhat rare) + elseif length(refined_edge_lids) == 3 + long_ref_edge_lid = c_to_longest_edge_lid[c] + cell_color[c] = BLUE_DOUBLE + Int(long_ref_edge_lid - 1) + end + end + # Setup color_rrules for the CompressedArray + num_rr = 1 + green_offset + blue_offset + blue_double_offset # GREEN+BLUE+DOUBLE_BLUE + T = typeof(WhiteRefinementRule(first(polys))) # Needed to make eltype(color_rrules) concrete + color_rrules = Vector{T}(undef,num_rr) + p = polys[1] # Hardcoded for TRI + color_rrules[WHITE] = WhiteRefinementRule(p) + for e in 1:num_faces(p,1) + color_rrules[GREEN+e-1] = GreenRefinementRule(p,e) + end + # blue_idx corresponds to offset into BLUE_dict + blue_idx = 0 + for e1 in 1:num_faces(p, 1) + for e2 in 1:num_faces(p, 1) + if e1 != e2 + color_rrules[BLUE + blue_idx] = BlueRefinementRule(p, e1, e2) + blue_idx += 1 + end + end + end + for e in 1:num_faces(p, 1) + color_rrules[BLUE_DOUBLE + e - 1] = BlueDoubleRefinementRule(p, e) + end + rrules = CompressedArray(color_rrules,cell_color) + face_list = _redgreen_refined_faces_list(topo, rrules, is_refined) + return rrules, face_list +end + +function setup_edge_based_rrules(::RedGreenRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc nC = num_cells(topo) nE = num_faces(topo,1) c2e_map = get_faces(topo,Dc,1) @@ -259,6 +398,8 @@ end abstract type EdgeBasedRefinementRule <: RefinementRuleType end struct RedRefinement <: EdgeBasedRefinementRule end +struct BlueRefinement{P, Q} <: EdgeBasedRefinementRule end +struct BlueDoubleRefinement{P} <: EdgeBasedRefinementRule end struct GreenRefinement{P} <: EdgeBasedRefinementRule end _has_interior_point(rr::RefinementRule) = _has_interior_point(rr,RefinementRuleType(rr)) @@ -273,8 +414,8 @@ function WhiteRefinementRule(p::Polytope) end """ -Edge-based RefinementRule where a new vertex is added to -each edge of the original Polytope. +Edge-based RefinementRule where a new vertex is added to +each edge of the original Polytope. """ function RedRefinementRule(p::Polytope) @notimplementedif (p ∉ [TRI,QUAD,HEX]) @@ -450,7 +591,7 @@ function get_d_to_face_to_parent_face(rr::RefinementRule,::RedRefinement) end """ -Edge-based RefinementRule where a new vertex is added to +Edge-based RefinementRule where a new vertex is added to a single edge of the original Polytope. """ function GreenRefinementRule(p::Polytope{2},ref_edge::Integer) @@ -468,8 +609,8 @@ 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. + # 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] @@ -494,7 +635,7 @@ function _get_green_vertex_permutation(p::Polytope{2},ref_edge::Integer) if p == TRI perm = circshift([1,2,3],ref_edge-3) elseif p == QUAD - # Vertices and edges are not circularly labeled in QUAD, so + # Vertices and edges are not circularly labeled in QUAD, so # we have to be inventive: nodes = [1,2,4,3] # Vertices in circular order nperm = [2,0,-1,1] # Number of perms depending on ref edge @@ -506,9 +647,94 @@ function _get_green_vertex_permutation(p::Polytope{2},ref_edge::Integer) end """ -Provided the gids of the coarse cell faces as a Tuple(vertices,edges,cell) and the -edge-based RefinementRule we want to apply, returns the connectivity of all the refined -children cells. +Edge-based RefinementRule where two vertex edges are added, +where the long edge is bisected first +""" +function BlueRefinementRule(p::Polytope{2}, long_ref_edge::Integer, short_ref_edge::Integer) + @notimplementedif (p ∉ [TRI]) + + faces_list = (collect(1:num_faces(p,0)),[long_ref_edge, short_ref_edge],[]) + coords = get_new_coordinates_from_faces(p,faces_list) + polys, cell_types, conn = _get_blue_refined_connectivity(p,long_ref_edge, short_ref_edge) + reffes = map(x->LagrangianRefFE(Float64,x,1),polys) + + ref_grid = UnstructuredDiscreteModel(UnstructuredGrid(coords,conn,reffes,cell_types)) + return RefinementRule(BlueRefinement{long_ref_edge, short_ref_edge}(),p,ref_grid) +end + +function _get_blue_refined_connectivity(p::Polytope{2}, long_ref_edge, short_ref_edge) + # Only implemented for triangles + if p == TRI + polys = [TRI] + cell_type = [1, 1, 1] + unmarked_edge = setdiff([1,2,3], [long_ref_edge, short_ref_edge])[1] + # Correspondence between an edge and the vertex opposite it + edge_to_opposite_vertex = Dict(1 => 3, 2 => 2, 3 => 1) + e2on = edge_to_opposite_vertex + conn_data = [e2on[long_ref_edge], 4, 5, + e2on[unmarked_edge], 4, 5, + sort([e2on[long_ref_edge], e2on[short_ref_edge], 4])...] + conn_ptrs = [1, 4, 7, 10] + return polys, cell_type, Table(conn_data,conn_ptrs) + end + @notimplemented +end + +function _get_blue_vertex_permutation(p::Polytope{2},ref_edge::Integer) + if p == TRI + perm = circshift([1,2,3],ref_edge-3) + else + @notimplemented + end + return perm +end + +""" +Edge-based RefinementRule where three vertex edges are added, +where the long edge is bisected first +""" +function BlueDoubleRefinementRule(p::Polytope{2}, long_ref_edge::Integer) + @notimplementedif (p ∉ [TRI]) + otherfaces = setdiff(collect(1:num_faces(p,0)), long_ref_edge) |> sort + faces_list = (collect(1:num_faces(p,0)),[otherfaces..., long_ref_edge],[]) + coords = get_new_coordinates_from_faces(p,faces_list) + polys, cell_types, conn = _get_blue_double_refined_connectivity(p,long_ref_edge) + reffes = map(x->LagrangianRefFE(Float64,x,1),polys) + ref_grid = UnstructuredDiscreteModel(UnstructuredGrid(coords,conn,reffes,cell_types)) + return RefinementRule(BlueDoubleRefinement{long_ref_edge}(),p,ref_grid) +end + +function _get_blue_double_refined_connectivity(p::Polytope{2}, long_ref_edge) + if p == TRI + polys = [TRI] + cell_type = [1, 1, 1, 1] + # Hardcoded since there are only three cases + if long_ref_edge == 3 + conn_data = [3, 5, 6, + 1, 5, 6, + 1, 4, 6, + 2, 4, 6] + elseif long_ref_edge == 2 + conn_data = [3, 5, 6, + 2, 5, 6, + 2, 4, 6, + 1, 4, 6] + elseif long_ref_edge == 1 + conn_data = [1, 4, 6, + 3, 4, 6, + 3, 5, 6, + 2, 5, 6] + end + conn_ptrs = [1, 4, 7, 10, 13] + return polys, cell_type, Table(conn_data,conn_ptrs) + end + @notimplemented +end + +""" +Provided the gids of the coarse cell faces as a Tuple(vertices,edges,cell) and the +edge-based RefinementRule we want to apply, returns the connectivity of all the refined +children cells. """ function get_relabeled_connectivity(rr::RefinementRule,faces_gids) return get_relabeled_connectivity(RefinementRuleType(rr),rr,faces_gids) @@ -537,6 +763,22 @@ function get_relabeled_connectivity(::GreenRefinement{N},rr::RefinementRule{P},f return Table(new_data,conn.ptrs) end +function get_relabeled_connectivity(::BlueRefinement{N, M},rr::RefinementRule{P},faces_gids) where {N, M, P<:Polytope{2}} + conn = rr.ref_grid.grid.cell_node_ids + gids = [faces_gids[1]...,faces_gids[2][N], faces_gids[2][M]] + new_data = lazy_map(Reindex(gids),conn.data) + return Table(new_data,conn.ptrs) +end + +function get_relabeled_connectivity(::BlueDoubleRefinement{N},rr::RefinementRule{P},faces_gids) where {N, P<:Polytope{2}} + conn = rr.ref_grid.grid.cell_node_ids + # Sort needed for non-longest edge gids + other_ids = setdiff(collect(1:3), N) |> sort + gids = [faces_gids[1]..., faces_gids[2][other_ids]..., faces_gids[2][N]] + new_data = lazy_map(Reindex(gids),conn.data) + return Table(new_data,conn.ptrs) +end + function counts_to_ptrs(counts::Vector{<:Integer}) n = length(counts) ptrs = Vector{Int32}(undef,n+1) diff --git a/src/FESpaces/ConstantFESpaces.jl b/src/FESpaces/ConstantFESpaces.jl index a491e47bf..d8251f341 100644 --- a/src/FESpaces/ConstantFESpaces.jl +++ b/src/FESpaces/ConstantFESpaces.jl @@ -30,7 +30,7 @@ cell_basis=SingleFieldFEBasis( cell_dof_basis_array=lazy_map(get_dof_basis,cell_reffe) cell_dof_basis=CellDof(cell_dof_basis_array,Triangulation(model),ReferenceDomain()) -cell_dof_ids=Fill(Int32[1],num_cells(model)) +cell_dof_ids=Fill(Int32(1):Int32(num_components(field_type)),num_cells(model)) A=typeof(cell_basis) B=typeof(cell_dof_basis) C=typeof(cell_dof_ids) diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl index bf6f5f16a..8b3ad50c4 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -157,7 +157,6 @@ export compute_vertex_node export get_node_face_owner export compute_node_face_owner export get_triangulation -export newest_vertex_bisection export UnstructuredDiscreteModel export CartesianDiscreteModel @@ -222,8 +221,6 @@ include("GridPortions.jl") include("DiscreteModelPortions.jl") -include("NewestVertexBisection.jl") - include("Triangulations.jl") include("BoundaryTriangulations.jl") diff --git a/src/Geometry/NewestVertexBisection.jl b/src/Geometry/NewestVertexBisection.jl deleted file mode 100644 index 23cfe2c5b..000000000 --- a/src/Geometry/NewestVertexBisection.jl +++ /dev/null @@ -1,638 +0,0 @@ -# This implementation is based on the following article: -# -# LONG CHEN (2008). SHORT IMPLEMENTATION OF BISECTION IN MATLAB. -# In Recent Advances in Computational Sciences. WORLD SCIENTIFIC. -# DOI: 10.1142/9789812792389_0020 - - -using SparseArrays -using Random -using AbstractTrees - -# The following was taken directly from Tim Holy's example at -# https://github.com/JuliaCollections/AbstractTrees.jl/blob/master/examples/binarytree_core.jl -##################################################################################### -mutable struct BinaryNode{T} - data::T - parent::BinaryNode{T} - left::BinaryNode{T} - right::BinaryNode{T} - - # Root constructor - BinaryNode{T}(data) where {T} = new{T}(data) - # Child node constructor - BinaryNode{T}(data, parent::BinaryNode{T}) where {T} = new{T}(data, parent) -end -BinaryNode(data) = BinaryNode{typeof(data)}(data) - -function leftchild(data, parent::BinaryNode) - !isdefined(parent, :left) || error("left child is already assigned") - node = typeof(parent)(data, parent) - parent.left = node -end -function rightchild(data, parent::BinaryNode) - !isdefined(parent, :right) || error("right child is already assigned") - node = typeof(parent)(data, parent) - parent.right = node -end - -function AbstractTrees.children(node::BinaryNode) - if isdefined(node, :left) - if isdefined(node, :right) - return (node.left, node.right) - end - return (node.left,) - end - isdefined(node, :right) && return (node.right,) - return () -end - -AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data) - -Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where {T} = BinaryNode{T} -Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where {T} = Base.HasEltype() -##################################################################################### - -_are_parallel(v, w) = v[1] * w[2] == v[2] * w[1] - -_get_midpoint(ngon::AbstractArray{<:VectorValue}) = sum(ngon) / length(ngon) - -_shift_to_first(v::AbstractVector{T}, i::T) where {T<:Integer} = circshift(v, -(i - 1)) - -function _print_forest(forest::AbstractArray{<:BinaryNode}) - num_leaves = 0 - println("_____________________") - for root_cell in forest - print_tree(root_cell) - num_leaves += length(collect(Leaves(root_cell))) - end - println("_____________________") - @show num_leaves -end - -function _set_d_to_dface_to_old_node!( - d_to_dface_to_oldid::AbstractArray{<:AbstractArray}, - d_to_dface_to_olddim::AbstractArray{<:AbstractArray}, - topo::GridTopology, - node_to_bis_edge::Dict, - num_new_nodes::Integer, -) - - edge_to_node = get_faces(topo, 1, 0) - push!(d_to_dface_to_oldid, Vector{}(undef, num_new_nodes)) - push!(d_to_dface_to_olddim, Vector{}(undef, num_new_nodes)) - # Loop over all new node IDs - for node_id = 1:num_new_nodes - # New node bisected an edge - if node_id in keys(node_to_bis_edge) - bis_edge = node_to_bis_edge[node_id] - bis_edge_index = findfirst(edge -> edge == bis_edge, edge_to_node) - d_to_dface_to_oldid[1][node_id] = bis_edge_index - d_to_dface_to_olddim[1][node_id] = 1 - # Node corresponds to previous node with same ID - # since node IDs are appended to the end - else - d_to_dface_to_oldid[1][node_id] = node_id - d_to_dface_to_olddim[1][node_id] = 0 - end - end -end - -function _set_d_to_dface_to_old_edge!( - d_to_dface_to_oldid::AbstractArray{<:AbstractArray}, - d_to_dface_to_olddim::AbstractArray{<:AbstractArray}, - forest::AbstractArray{<:BinaryNode}, - topo::GridTopology, - topo_ref::GridTopology, - vertices::AbstractArray{<:VectorValue}, -) - # Topological mappings for the refined mesh - edge_to_node_ref = get_faces(topo_ref, 1, 0) - cell_to_edge_ref = get_faces(topo_ref, 2, 1) - num_edges = length(edge_to_node_ref) - # Topological mappings for the old mesh - cell_to_edge = get_faces(topo, 2, 1) - edge_to_node = get_faces(topo, 1, 0) - cell_to_node = get_faces(topo, 2, 0) - push!(d_to_dface_to_oldid, Vector{}(undef, num_edges)) - push!(d_to_dface_to_olddim, Vector{}(undef, num_edges)) - for root_cell in forest - root_edges = cell_to_edge[root_cell.data] - root_nodes = cell_to_node[root_cell.data] - for leaf_cell in Leaves(root_cell) - # Set because of possible repeats for new cells that are neighbors - leaf_edges = Set(cell_to_edge_ref[leaf_cell.data]) - for leaf_edge in leaf_edges - leaf_edge_nodes = edge_to_node_ref[leaf_edge] - is_on_top_of_former_edge = false - # We need to check if this edge is parallel to any of the - # other edges in the old cell to see if it is on top of a former - # edge - if !isempty(intersect(leaf_edge_nodes, root_nodes)) - for root_edge in root_edges - root_edge_nodes = edge_to_node[root_edge] - root_edge_vector = vertices[root_edge_nodes[2]] - vertices[root_edge_nodes[1]] - leaf_edge_vector = vertices[leaf_edge_nodes[2]] - vertices[leaf_edge_nodes[1]] - if _are_parallel(root_edge_vector, leaf_edge_vector) - d_to_dface_to_oldid[2][leaf_edge] = root_edge - d_to_dface_to_olddim[2][leaf_edge] = 1 - is_on_top_of_former_edge = true - end - end - end - # Default case: this edge bisects a cell, i.e. not on top of - # former edge - if !is_on_top_of_former_edge - d_to_dface_to_oldid[2][leaf_edge] = root_cell.data - d_to_dface_to_olddim[2][leaf_edge] = 2 - end - end - end - end -end - -function _set_d_to_dface_to_old_cell!( - d_to_dface_to_oldid::AbstractArray{<:AbstractArray}, - d_to_dface_to_olddim::AbstractArray{<:AbstractArray}, - forest::AbstractArray{<:BinaryNode}, - topo_ref::GridTopology, -) - cell_to_node_ref = get_faces(topo_ref, 2, 0) - num_cells = length(cell_to_node_ref) - push!(d_to_dface_to_oldid, Vector{}(undef, num_cells)) - push!(d_to_dface_to_olddim, Vector{}(undef, num_cells)) - # Can basically just use the forest. All new cells should point - # to the root of their tree - leaf_ids = Int32[] - for root_cell in forest - root_id = root_cell.data - for leaf_cell in Leaves(root_cell) - leaf_id = leaf_cell.data - push!(leaf_ids, leaf_id) - d_to_dface_to_oldid[3][leaf_id] = root_id - d_to_dface_to_olddim[3][leaf_id] = 2 - end - end -end - -function _create_d_to_dface_to_old( - forest::AbstractArray{<:BinaryNode}, - topo::GridTopology, - topo_ref::GridTopology, - vertices::AbstractArray{<:VectorValue}, - node_to_bis_edge::Dict, -) - d_to_dface_to_oldid = Vector{Vector}() - d_to_dface_to_olddim = Vector{Vector}() - num_new_nodes = length(vertices) - # The order is important here: - # Node - _set_d_to_dface_to_old_node!( - d_to_dface_to_oldid, - d_to_dface_to_olddim, - topo, - node_to_bis_edge, - num_new_nodes, - ) - # Edge - _set_d_to_dface_to_old_edge!( - d_to_dface_to_oldid, - d_to_dface_to_olddim, - forest, - topo, - topo_ref, - vertices, - ) - # Cell - _set_d_to_dface_to_old_cell!( - d_to_dface_to_oldid, - d_to_dface_to_olddim, - forest, topo_ref) - # Checks to make sure all new dims and ids are filled - for d = 0:2 - @test undef ∉ d_to_dface_to_oldid[d+1] - @test undef ∉ d_to_dface_to_olddim[d+1] - @test length(d_to_dface_to_olddim[d+1]) == length(d_to_dface_to_olddim[d+1]) - end - d_to_dface_to_olddim, d_to_dface_to_oldid -end - -function _propagate_labeling(model, d_to_dface_to_olddim, d_to_dface_to_oldid) - labels = get_face_labeling(model) - labels_ref = FaceLabeling(length.(d_to_dface_to_oldid)) - @test isempty(labels_ref.tag_to_entities) - # Copy entities - for entities in labels.tag_to_entities - push!(labels_ref.tag_to_entities, copy(entities)) - end - # Copy refs - @test isempty(labels_ref.tag_to_name) - for name in labels.tag_to_name - push!(labels_ref.tag_to_name, name) - end - # Populate new entity IDs using d_to_dface_to_oldid and - # d_to_dface_to_olddim - for d = 0:num_dims(model) - oldid_d = d_to_dface_to_oldid[d+1] - olddim_d = d_to_dface_to_olddim[d+1] - @test length(oldid_d) == length(olddim_d) - for i = 1:length(oldid_d) - old_id = oldid_d[i] - old_dim = olddim_d[i] - old_entity_id = labels.d_to_dface_to_entity[old_dim+1][old_id] - labels_ref.d_to_dface_to_entity[d+1][i] = old_entity_id - end - end - labels_ref -end - -function _sort_longest_edge!( - elem::AbstractVector{<:AbstractVector{T}}, - node::AbstractVector{<:VectorValue}, -) where {T<:Integer} - NT = length(elem) - edgelength = zeros(NT, 3) - for t = 1:NT - elem_t = elem[t][:] - for (j, e) in enumerate(elem_t) - arr = filter(x -> x != e, elem_t) - diff = norm(node[arr[1]] - node[arr[2]]) .^ 2 - edgelength[t, j] = diff - end - end - max_indices = findmax(edgelength, dims = 2)[2] - for t = 1:NT - shifted = _shift_to_first(elem[t][:], T(max_indices[t][2])) - elem[t] = shifted - end -end - -function _setup_markers_and_nodes!( - node::AbstractVector{<:VectorValue}, - elem::AbstractVector{<:AbstractVector{T}}, - d2p::SparseMatrixCSC{T,T}, - dualedge::SparseMatrixCSC{T}, - NE::T, - η_arr::AbstractArray, - θ::AbstractFloat, -) where {T<:Integer} - total_η = sum(η_arr) - partial_η = 0 - sorted_η_idxs = sortperm(-η_arr) - markers = zeros(T, NE) - # Loop over global triangle indices - for t in sorted_η_idxs - if (partial_η >= θ * total_η) - break - end - need_to_mark = true - # Get triangle index with next largest error - while (need_to_mark) - # Base point - base = d2p[elem[t][2], elem[t][3]] - # Already marked - if markers[base] > 0 - need_to_mark = false - else - # Get the estimator contribution for the current triangle - partial_η = partial_η + η_arr[t] - # Increase the number of nodes to add new midpoint - N = size(node, 1) + 1 - # The markers of the current elements is this node - markers[d2p[elem[t][2], elem[t][3]]] = N - # Coordinates of new node - elem2 = elem[t][2] - elem3 = elem[t][3] - midpoint = _get_midpoint(node[[elem2, elem3], :]) - node = push!(node, midpoint) - t = dualedge[elem[t][3], elem[t][2]] - # There is no dual edge here, go to next triangle index - if t == 0 - need_to_mark = false - end - end - end - end - node, markers -end - -function _divide!(elem::AbstractVector, t::T, p::AbstractVector{T}) where {T<:Integer} - new_row = [p[4], p[3], p[1]] - update_row = [p[4], p[1], p[2]] - push!(elem, new_row) - elem[t] = update_row - elem -end - -function _bisect( - d2p::SparseMatrixCSC{T,T}, - elem::AbstractVector, - markers::AbstractVector{T}, - NT::T, -) where {T<:Integer} - forest = Vector{BinaryNode}() - for t in UnitRange{T}(1:NT) - base = d2p[elem[t][2], elem[t][3]] - if (markers[base] > 0) - newnode = BinaryNode(t) - p = vcat(elem[t][:], markers[base]) - elem = _divide!(elem, t, p) - cur_size::T = size(elem, 1) - # l and r are BinaryNodes for building forest - l = leftchild(t, newnode) - r = rightchild(cur_size, newnode) - # left and right are point indices - left = d2p[p[1], p[2]] - right = d2p[p[3], p[1]] - if (markers[right] > 0) - leftchild(cur_size, r) - # Need to handle edge case where both markers are positive. - # In this case, we have to offset by 2 because the other branch - # already took the next index. - if markers[left] > 0 - rightchild(cur_size + 2, r) - else - rightchild(cur_size + 1, r) - end - elem = _divide!(elem, cur_size, [p[4], p[3], p[1], markers[right]]) - end - if (markers[left] > 0) - leftchild(t, l) - rightchild(cur_size + 1, l) - elem = _divide!(elem, t, [p[4], p[1], p[2], markers[left]]) - end - push!(forest, newnode) - else - push!(forest, BinaryNode(t)) - end - end - elem, forest -end - -function _build_edges(elem::AbstractVector{<:AbstractVector{T}}) where {T<:Integer} - edge = zeros(T, 3 * length(elem), 2) - for t = 1:length(elem) - off = 3 * (t - 1) - edge[off+1, :] = [elem[t][1] elem[t][2]] - edge[off+2, :] = [elem[t][1] elem[t][3]] - edge[off+3, :] = [elem[t][2] elem[t][3]] - end - edge = unique(sort!(edge, dims = 2), dims = 1) -end - -function _build_directed_dualedge( - elem::AbstractVector{<:AbstractVector{T}}, - N::T, - NT::T, -) where {T<:Integer} - dualedge = spzeros(T, N, N) - for t = 1:NT - dualedge[elem[t][1], elem[t][2]] = t - dualedge[elem[t][2], elem[t][3]] = t - dualedge[elem[t][3], elem[t][1]] = t - end - dualedge -end - -function _dual_to_primal(edge::Matrix{T}, NE::T, N::T) where {T<:Integer} - d2p = spzeros(T, T, N, N) - for k = 1:NE - i = edge[k, 1] - j = edge[k, 2] - d2p[i, j] = k - d2p[j, i] = k - end - d2p -end - -function _is_against_top(face::Table{<:Integer}, top::GridTopology, d::Integer) - face_vec = sort.([face[i][:] for i = 1:size(face, 1)]) - face_top = sort.(get_faces(top, d, 0)) - issetequal_bitvec = issetequal(face_vec, face_top) - all(issetequal_bitvec) -end - -function _is_against_top(face::Matrix{<:Integer}, top::GridTopology, d::Integer) - face_vec = sort.([face[i, :] for i = 1:size(face, 1)]) - face_top = sort.(get_faces(top, d, 0)) - issetequal_bitvec = issetequal(face_vec, face_top) - all(issetequal_bitvec) -end - -# For initial labeling -function _sort_ccw(cell_coords::AbstractVector{<:VectorValue}) - midpoint = _get_midpoint(cell_coords) - offset_coords = cell_coords .- midpoint - sortperm(offset_coords, by = v -> atan(v[2], v[1])) -end - -# For initial labeling -function _sort_cell_node_ids_ccw!( - cell_node_ids::AbstractVector{<:AbstractVector{T}}, - node_coords::AbstractVector{<:VectorValue}, -) where {T<:Integer} - for (i, cell) in enumerate(cell_node_ids) - cell_coords = node_coords[cell] - perm = _sort_ccw(cell_coords) - cell_node_ids[i] = cell[perm] - end -end - -function _markers_to_node_to_bis_edge(markers, edge) - node_to_bis_edge = Dict() - for (i, marker) in enumerate(markers) - # Edge was bisected - if marker != 0 - node_to_bis_edge[marker] = edge[i, :] - end - end - node_to_bis_edge -end - -""" -Lowest level interface to the newest vertex bisection algorithm. This step -takes places after unpacking the Grid. It also uses the `should_sort` to -determine if an initial sorting by the longest edge (for nonuniform meshes) is -necessary. - -`node_coords, cell_node_ids -> node_coords, cell_node_ids` - -# Arguments - - -`node_coords::AbstractVector{<:VectorValue}`: The vector of d-dimensional - nodal coordinates stored as `VectorValue`s - - -`cell_node_ids::AbstractVector{<:AbstractVector{T}}`: Contains the elements - as ordered tuples of nodal indices. - - -`η_arr::AbstractArray`: The values of the estimator on each cell of the domain - - -`θ::AbstractArray`: Dörfler marking parameter: 0 = no refinement, - 1=uniform refinement. - -""" -function newest_vertex_bisection( - node_coords::AbstractVector{<:VectorValue}, - cell_node_ids::AbstractVector{<:AbstractVector{T}}, - η_arr::AbstractArray, - θ::AbstractFloat, -) where {T<:Integer} - # Number of nodes - N::T = size(node_coords, 1) - # Long Chen terminology - elem = cell_node_ids - # Number of cells (triangles in 2D) - NT::T = size(elem, 1) - @assert length(η_arr) == NT - # Long Chen terminology - edge = _build_edges(elem) - NE::T = size(edge, 1) - # Long Chen terminology - dualedge = _build_directed_dualedge(elem, N, NT) - d2p = _dual_to_primal(edge, NE, N) - node_coords, markers = - _setup_markers_and_nodes!(node_coords, elem, d2p, dualedge, NE, η_arr, θ) - # For performance since we need to push! to elem - elem = Vector{Vector}(elem) - cell_node_ids, forest = _bisect(d2p, elem, markers, NT) - node_to_bis_edge = _markers_to_node_to_bis_edge(markers, edge) - node_coords, cell_node_ids, forest, node_to_bis_edge -end - -""" -Middle level interface to the newest vertex bisection algorithm. This step -takes places after unpacking the DiscreteModel and performs sorting if -necessary. It maps - -`Grid -> Grid, buffer` - -# Arguments - - -`grid::Grid`: The current Grid. - - -`η_arr::AbstractArray`: The values of the estimator on each cell of the domain - - -`θ::AbstractArray`: Dörfler marking parameter: 0 = no refinement, - 1=uniform refinement. - -""" -function newest_vertex_bisection(grid::Grid, η_arr::AbstractArray, θ::AbstractFloat) - node_coords = get_node_coordinates(grid) - # Need "un lazy" version for resize! - node_coords = [v for v in node_coords] - #@show cell_node_ids = get_cell_node_ids(grid) - cell_node_ids = get_cell_node_ids(grid) - # Convert from Table to Vector{Vector} - cell_node_ids = [v for v in cell_node_ids] - # Should always sort on the first iteration - _sort_cell_node_ids_ccw!(cell_node_ids, node_coords) - _sort_longest_edge!(cell_node_ids, node_coords) - node_coords_ref, cell_node_ids_unsort, forest, node_to_bis_edge = - newest_vertex_bisection(node_coords, cell_node_ids, η_arr, θ) - reffes = get_reffes(grid) - cell_types = get_cell_type(grid) - cell_types = [c for c in cell_types] - # TODO : Is this ok since we only handle triangular mesh for now? - new_cell_types = fill(1, length(cell_node_ids_unsort) - length(cell_node_ids)) - append!(cell_types, new_cell_types) - cell_node_ids_ref = Table([c for c in cell_node_ids_unsort]) - grid_ref_unsort = UnstructuredGrid(node_coords_ref, cell_node_ids_ref, reffes, cell_types) - # Other buffer information can be added later - buffer = (; grid_ref_unsort) - sort!.(cell_node_ids_unsort) - cell_node_ids_ref = Table([c for c in cell_node_ids_unsort]) - grid_ref = UnstructuredGrid(node_coords_ref, cell_node_ids_ref, reffes, cell_types) - grid_ref, buffer, forest, node_to_bis_edge -end - -""" -The newest vertex bisection algorithm provides a method of local refinement -without creating hanging nodes. For now, only 2D simplicial meshes are -supported. - -This is the highest level version of the function, it maps - -`DiscreteModel -> DiscreteModel, buffer` - - -# Arguments - - -`model::DiscreteModel`: The current DiscreteModel to be refined. - - -`η_arr::AbstractArray`: The values of the estimator on each cell of the domain - i.e. one should have `length(η_arr) == num_cells(model)` - - -`θ::AbstractArray=1.0`: Dörfler marking parameter: 0 = no refinement, - 1=uniform refinement. - -""" -function newest_vertex_bisection( - model::DiscreteModel, - η_arr::AbstractArray; - θ = 1.0, # corresponds to uniform refinement -) - @assert length(η_arr) == num_cells(model) - grid = get_grid(model) - topo = GridTopology(grid) - grid_ref, buffer, forest, node_to_bis_edge = newest_vertex_bisection(grid, η_arr, θ) - topo_ref = GridTopology(grid_ref) - d_to_dface_to_olddim, d_to_dface_to_oldid = _create_d_to_dface_to_old( - forest, - topo, - topo_ref, - get_node_coordinates(grid_ref), - node_to_bis_edge, - ) - labels_ref = _propagate_labeling(model, d_to_dface_to_olddim, d_to_dface_to_oldid) - DiscreteModel(grid_ref, topo_ref, labels_ref), buffer -end - -""" -The newest vertex bisection algorithm provides a method of local refinement -without creating hanging nodes. For now, only 2D simplicial meshes are -supported. - -This is the highest level version of the function, it maps - -`DiscreteModel, buffer -> DiscreteModel, buffer` - - -# Arguments - - -`model::DiscreteModel`: The current DiscreteModel to be refined. - - -`buffer::NamedTuple`: The buffer providing all the itermediate information - between refinements. For now just includes a Grid with orderings necessary - for the algorithm. - - -`η_arr::AbstractArray`: The values of the estimator on each cell of the domain - i.e. one should have `length(η_arr) == num_cells(model)` - - -`θ::AbstractArray=1.0`: Dörfler marking parameter: 0 = no refinement, - 1=uniform refinement. - -""" -function newest_vertex_bisection( - model::DiscreteModel, - buffer::NamedTuple, - η_arr::AbstractArray; - θ = 1.0, # corresponds to uniform refinement -) - @assert length(η_arr) == num_cells(model) - grid = get_grid(model) - grid_ref_unsort = buffer.grid_ref_unsort - topo = GridTopology(grid) - grid_ref, buffer, forest, node_to_bis_edge = - newest_vertex_bisection(grid_ref_unsort, η_arr, θ) - topo_ref = GridTopology(grid_ref) - d_to_dface_to_olddim, d_to_dface_to_oldid = _create_d_to_dface_to_old( - forest, - topo, - topo_ref, - get_node_coordinates(grid_ref), - node_to_bis_edge, - ) - labels_ref = _propagate_labeling(model, d_to_dface_to_olddim, d_to_dface_to_oldid) - model_ref = DiscreteModel(grid_ref, topo_ref, labels_ref) - model_ref, buffer -end diff --git a/test/AdaptivityTests/EdgeBasedRefinementTests.jl b/test/AdaptivityTests/EdgeBasedRefinementTests.jl index c77977107..9e525318e 100644 --- a/test/AdaptivityTests/EdgeBasedRefinementTests.jl +++ b/test/AdaptivityTests/EdgeBasedRefinementTests.jl @@ -145,5 +145,21 @@ trian6 = Triangulation(ref_model6.model) visualize && writevtk(trian6,"test/AdaptivityTests/ref_model6") test_grid_transfers(2,model2,ref_model6,1) +# Test newest vertex bisection i.e. longest edge bisection +# Refine all edges using NVB +ref_model7 = refine(model2, refinement_method = "nvb") +trian7 = Triangulation(ref_model7.model) +visualize && writevtk(trian7, "test/AdaptivityTests/ref_model7") +test_grid_transfers(2, model2, ref_model7, 1) + +# Refine all edges using NVB +# Mark edges such that blue, and double_blue refinement are triggered +ref_model8 = refine(model2, refinement_method = "nvb", cells_to_refine = [4, 9]) +trian8 = Triangulation(ref_model8) +visualize && writevtk(trian8, "test/AdaptivityTests/ref_model8") +ref_model9 = refine(ref_model8, refinement_method = "nvb", cells_to_refine = [1, 3, 4, 11]) +trian9 = Triangulation(ref_model9.model) +visualize && writevtk(trian9, "test/AdaptivityTests/ref_model9") +test_grid_transfers(2, ref_model8, ref_model9, 1) -end \ No newline at end of file +end diff --git a/test/FESpacesTests/ConstantFESpaceTests.jl b/test/FESpacesTests/ConstantFESpaceTests.jl index ecb84c340..980ee6124 100644 --- a/test/FESpacesTests/ConstantFESpaceTests.jl +++ b/test/FESpacesTests/ConstantFESpaceTests.jl @@ -1,4 +1,3 @@ - module ConstantFESpacesTests using Gridap @@ -30,4 +29,13 @@ uh = solve(op) @assert sum(∫((uh[1]-u)*(uh[1]-u))dΩ) < 1.0e-14 abs(sum(∫(uh[2])dΩ)) < 1.0e-12 +Λ2=ConstantFESpace(model,field_type=VectorValue{2,Float64}) +Gridap.FESpaces.test_fe_space(Λ2) +M2=TrialFESpace(Λ2) +a2(μ,λ) = ∫(λ⋅μ)dΩ +l2(λ) = ∫(VectorValue(0.0,0.0)⋅λ)dΩ +op2 = AffineFEOperator(a2,l2,M2,Λ2) +μ2h = solve(op2) +@assert sum(∫(μ2h⋅μ2h)dΩ) < 1.0e-12 + end # module diff --git a/test/GeometryTests/NewestVertexBisectionTests.jl b/test/GeometryTests/NewestVertexBisectionTests.jl deleted file mode 100644 index 85467c91c..000000000 --- a/test/GeometryTests/NewestVertexBisectionTests.jl +++ /dev/null @@ -1,125 +0,0 @@ -module NewestVertexBisectionTests - -using Test -using Random -using Gridap.Arrays -using Gridap.Fields -using Gridap.ReferenceFEs -using Gridap.Geometry -using Gridap.Visualization - - -# For testing only -abstract type Estimator end - -struct ConstantEst <: Estimator - val::Float64 -end - -struct RandomEst <: Estimator - function RandomEst(seed) - Random.seed!(seed) - new() - end -end - -function make_nvb_levels( - model::DiscreteModel, - Nsteps::Integer, - θ::AbstractFloat, - est::Estimator, -) - model_refs = Vector{DiscreteModel}(undef, Nsteps) - cell_map = get_cell_map(get_triangulation(model)) - ncells = length(cell_map) - η_arr = compute_estimator(est, ncells) - model_refs[1], buffer = newest_vertex_bisection(model, η_arr; θ = θ) - buffer = deepcopy(buffer) - for i = 1:(Nsteps - 1) - cell_map = get_cell_map(get_triangulation(model_refs[i])) - ncells = length(cell_map) - η_arr = compute_estimator(est, ncells) - model_refs[i + 1], buffer = - newest_vertex_bisection(model_refs[i], buffer, η_arr; θ = θ) - # Necessary to have each level stored separately - buffer = deepcopy(buffer) - end - model_refs -end - -function d_get_num_boundary_labels(labels, d) - entities = labels.d_to_dface_to_entity[d + 1] - names = labels.tag_to_name[entities] - bdry_names = filter(name -> name != "interior", names) - length(bdry_names) -end - - -# For testing only -compute_estimator(est::RandomEst, ncells) = rand(ncells) -compute_estimator(est::ConstantEst, ncells) = fill(est.val, ncells) - -domain = (0, 1, 0, 1) -partition = (1, 1) # Initial partition -Nsteps = 10 -est = ConstantEst(1.0) -θ = 1.0 -uniform_write_to_vtk = false -# Uniform refinement -let model = simplexify(CartesianDiscreteModel(domain, partition)) - model_refs = make_nvb_levels(model, Nsteps, θ, est) - for (n, model_ref) in enumerate(model_refs) - trian_ref = get_triangulation(model_ref) - if uniform_write_to_vtk - writevtk(trian_ref, "uniform$(string(n, pad=2))") - end - # Combinatorial checks for cells - let cell_map = get_cell_map(trian_ref) - ncells = length(cell_map) - @test ncells == 2^(n + 1) - end - # Combinatorial checks for nodes - let node_coords = get_node_coordinates(trian_ref) - ncoords = length(node_coords) - if isodd(n) - ncoords_true = Int(2 * (4^((n - 1) / 2) + 2^((n - 1) / 2)) + 1) - else - ncoords_true = Int(2^(n / 2) + 1)^2 - end - @test ncoords_true == ncoords - end - # Test labels are propagating properly - let labels = get_face_labeling(model_ref) - let d = 0 - @test Int(4*2^(floor(n / 2))) == d_get_num_boundary_labels(labels, d) - end - # Combinatorial check for face_labels (edge) - let d = 1 - @test Int(4*2^(floor(n / 2))) == d_get_num_boundary_labels(labels, d) - end - # Combinatorial check for cell labels - let d = 2 - @test 0 == d_get_num_boundary_labels(labels, d) - end - end - end -end - -# Nonuniform refinement. For now only visually checking conformity -#domain = (0, 1, 0, 1) -#partition = (1, 1) # Initial partition -#Nsteps = 15 -#seed = 5 -#est = RandomEst(seed) -#θ = 0.5 -#nonuniform_write_to_vtk = true -#model = simplexify(CartesianDiscreteModel(domain, partition)) -#@time model_refs = make_nvb_levels(model, Nsteps, θ, est) -#if nonuniform_write_to_vtk -# for (n, model_ref) in enumerate(model_refs) -# trian_ref = get_triangulation(model_ref) -# writevtk(trian_ref, "nonuniform$(string(n, pad=2))") -# end -#end - -end diff --git a/test/GeometryTests/runtests.jl b/test/GeometryTests/runtests.jl index f3499927b..ecb449ddf 100644 --- a/test/GeometryTests/runtests.jl +++ b/test/GeometryTests/runtests.jl @@ -38,6 +38,4 @@ using Test @testset "CompressedCellArrays" begin include("CompressedCellArraysTests.jl") end -@time @testset "Refinement" begin include("NewestVertexBisectionTests.jl") end - end # module diff --git a/test/ReferenceFEsTests/NedelecRefFEsTests.jl b/test/ReferenceFEsTests/NedelecRefFEsTests.jl index ec931f171..b94eeb134 100644 --- a/test/ReferenceFEsTests/NedelecRefFEsTests.jl +++ b/test/ReferenceFEsTests/NedelecRefFEsTests.jl @@ -156,6 +156,9 @@ ux = evaluate(shapes,x) gux = evaluate(gshapes,x) ndat = ["s$i"=>ux[:,i] for i in 1:num_dofs(reffe)] gndat = ["g$i"=>gux[:,i] for i in 1:num_dofs(reffe)] -writevtk(grid,"nede_tet_1",nodaldata=vcat(ndat,gndat)) +d = mktempdir() +f = joinpath(d, "nede_tet_1") +writevtk(grid,f,nodaldata=vcat(ndat,gndat)) +rm(d,recursive=true) end # module