diff --git a/NEWS.md b/NEWS.md index 94a5eeed9..f4e6ec555 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,11 +10,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - 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 + `src/Adaptivity/EdgeBasedRefinement.jl`. Since PR [#901](https://github.com/gridap/Gridap.jl/pull/901). +- When refining `DiscreteModels`, the `FaceLabeling` of the resulting `AdaptedDiscreteModel` will now correctly inhering the tags of the parent model. This has been made possible by the addition of the method `get_d_to_face_to_parent_face`. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886). +- Added support for mixed adaptivity (i.e coarsening and refining), as well as non-conforming adaptivity. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886). ### Changed @@ -26,9 +27,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fixed the method `get_normal_vector` for `AdaptedTriangulation`. The method `get_facet_normal` 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) - Fixed generation of Modal C0 bases for Julia 1.9. Since PR [#918](https://github.com/gridap/Gridap.jl/pull/918). +- Fixed some edge cases for `change_domain` between `AdaptedTriangulations` where inneficient coordinate transformations would be applied between physical and reference domains. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886). +- Fixed: Domain limits can now be of any type (notably, floats) when refining `CartesianDiscreteModels`. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886). ## [0.17.17] - 2023-02-24 diff --git a/src/Adaptivity/AdaptedDiscreteModels.jl b/src/Adaptivity/AdaptedDiscreteModels.jl index 25e6cea96..e00e8ddbc 100644 --- a/src/Adaptivity/AdaptedDiscreteModels.jl +++ b/src/Adaptivity/AdaptedDiscreteModels.jl @@ -103,24 +103,31 @@ function refine(model::CartesianDiscreteModel{Dc}, cell_partition::Tuple) where desc = Geometry.get_cartesian_descriptor(model) nC = desc.partition - # Refined model - domain = _get_cartesian_domain(desc) - model_ref = CartesianDiscreteModel(domain,cell_partition.*nC) - - # Glue + # Refinement Glue f2c_cell_map, fcell_to_child_id = _create_cartesian_f2c_maps(nC,cell_partition) - faces_map = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc] - reffe = LagrangianRefFE(Float64,first(get_polytopes(model)),1) - rrules = RefinementRule(reffe,cell_partition) + faces_map = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc] + reffe = LagrangianRefFE(Float64,first(get_polytopes(model)),1) + rrules = RefinementRule(reffe,cell_partition) glue = AdaptivityGlue(faces_map,fcell_to_child_id,rrules) + # Refined model + domain = _get_cartesian_domain(desc) + _model_ref = CartesianDiscreteModel(domain,cell_partition.*nC) + + # Propagate face labels + coarse_labels = get_face_labeling(model) + coarse_topo = get_grid_topology(model) + fine_topo = get_grid_topology(_model_ref) + fine_labels = _refine_face_labeling(coarse_labels,glue,coarse_topo,fine_topo) + + model_ref = CartesianDiscreteModel(get_grid(_model_ref),fine_topo,fine_labels) return AdaptedDiscreteModel(model_ref,model,glue) end function _get_cartesian_domain(desc::CartesianDescriptor{D}) where D origin = desc.origin corner = origin + VectorValue(desc.sizes .* desc.partition) - domain = Vector{Int}(undef,2*D) + domain = Vector{eltype(origin)}(undef,2*D) for d in 1:D domain[d*2-1] = origin[d] domain[d*2] = corner[d] @@ -159,3 +166,10 @@ end return f2c_map, child_map end) end + +function get_d_to_fface_to_cface(model::AdaptedDiscreteModel) + ftopo = get_grid_topology(get_model(model)) + ctopo = get_grid_topology(get_parent(model)) + glue = get_adaptivity_glue(model) + return get_d_to_fface_to_cface(glue,ctopo,ftopo) +end diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index 6f803ab6e..d092b7a64 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -183,8 +183,14 @@ end for sdomain in [:ReferenceDomain,:PhysicalDomain] for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTriangulation,:Triangulation),(:Triangulation,:AdaptedTriangulation)] + sstrian = (stype==:AdaptedTriangulation) ? :(strian.trian) : :(strian) + tttrian = (ttype==:AdaptedTriangulation) ? :(ttrian.trian) : :(ttrian) @eval begin - function CellData.change_domain(a::CellField,strian::$stype,::$sdomain,ttrian::$ttype,::PhysicalDomain) + function CellData.change_domain(a::CellField,strian::$stype,sd::$sdomain,ttrian::$ttype,::PhysicalDomain) + if (get_background_model(strian) === get_background_model(ttrian)) + return change_domain(a,$sstrian,sd,$tttrian,PhysicalDomain()) + end + a_ref = change_domain(a,ReferenceDomain()) atrian = change_domain(a_ref,strian,ReferenceDomain(),ttrian,ReferenceDomain()) return change_domain(atrian,PhysicalDomain()) @@ -281,11 +287,36 @@ function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTri return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain()) else - f_fine = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian)) + f_fine = Fill(Fields.ConstantField(0.0),num_cells(ftrian)) return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain()) end end +function change_domain_o2n( + f_old, + old_trian::Triangulation{Dc}, + new_trian::AdaptedTriangulation, + glue::AdaptivityGlue{<:MixedGlue}) where {Dc} + + oglue = get_glue(old_trian,Val(Dc)) + nglue = get_glue(new_trian,Val(Dc)) + + @notimplementedif num_point_dims(old_trian) != Dc + @notimplementedif isa(nglue,Nothing) + + if (num_cells(old_trian) != 0) + # If mixed refinement/coarsening, then f_c2f is a Table + f_old_data = CellData.get_data(f_old) + f_c2f = c2f_reindex(f_old_data,glue) + new_rrules = get_new_cell_refinement_rules(glue) + field_array = lazy_map(OldToNewField, f_c2f, new_rrules, glue.n2o_cell_to_child_id) + return CellData.similar_cell_field(f_old,field_array,new_trian,ReferenceDomain()) + else + f_new = Fill(Fields.ConstantField(0.0),num_cells(new_trian)) + return CellData.similar_cell_field(f_old,f_new,new_trian,ReferenceDomain()) + end +end + """ Given a AdaptivityGlue and a CellField defined on the child(new) mesh, returns an equivalent CellField on the parent(old) mesh. @@ -310,8 +341,9 @@ function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Trian # f_c2f[i_coarse] = [f_fine[i_fine_1], ..., f_fine[i_fine_nChildren]] f_c2f = f2c_reindex(fine_mface_to_field,glue) - rrules = get_old_cell_refinement_rules(glue) - coarse_mface_to_field = lazy_map(FineToCoarseField,f_c2f,rrules) + child_ids = f2c_reindex(glue.n2o_cell_to_child_id,glue) + rrules = get_old_cell_refinement_rules(glue) + coarse_mface_to_field = lazy_map(FineToCoarseField,f_c2f,rrules,child_ids) ### Old Model -> Old Triangulation coarse_tface_to_field = lazy_map(Reindex(coarse_mface_to_field),cglue.tface_to_mface) @@ -319,7 +351,7 @@ function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Trian return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) else - f_coarse = Fill(Gridap.Fields.ConstantField(0.0),num_cells(fcoarse)) + f_coarse = Fill(Fields.ConstantField(0.0),num_cells(fcoarse)) return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) end end diff --git a/src/Adaptivity/Adaptivity.jl b/src/Adaptivity/Adaptivity.jl index 9cc5971a4..f16d80d8e 100644 --- a/src/Adaptivity/Adaptivity.jl +++ b/src/Adaptivity/Adaptivity.jl @@ -42,6 +42,7 @@ export change_domain, move_contributions include("RefinementRules.jl") include("FineToCoarseFields.jl") +include("OldToNewFields.jl") include("FineToCoarseReferenceFEs.jl") include("AdaptivityGlues.jl") include("AdaptedDiscreteModels.jl") diff --git a/src/Adaptivity/AdaptivityGlues.jl b/src/Adaptivity/AdaptivityGlues.jl index b65f0ca77..2b08b90a1 100644 --- a/src/Adaptivity/AdaptivityGlues.jl +++ b/src/Adaptivity/AdaptivityGlues.jl @@ -22,13 +22,20 @@ struct AdaptivityGlue{GT,Dc,A,B,C,D,E} <: GridapType is_refined :: E function AdaptivityGlue(n2o_faces_map::Vector{<:Union{AbstractVector{<:Integer},Table{<:Integer}}}, - n2o_cell_to_child_id::AbstractVector{<:Integer}, + n2o_cell_to_child_id::Union{AbstractVector{<:Integer},Table{<:Integer}}, refinement_rules::AbstractVector{<:RefinementRule}) Dc = length(n2o_faces_map)-1 is_refined = select_refined_cells(n2o_faces_map[Dc+1]) o2n_faces_map = get_o2n_faces_map(n2o_faces_map[Dc+1]) GT = all(is_refined) ? RefinementGlue : MixedGlue + if isa(GT(),RefinementGlue) + @assert isa(n2o_faces_map,Vector{<:AbstractVector{<:Integer}}) + @assert isa(n2o_cell_to_child_id,AbstractVector{<:Integer}) + else + @assert isa(n2o_faces_map,Vector{<:Table{<:Integer}}) + @assert isa(n2o_cell_to_child_id,Table{<:Integer}) + end A = typeof(n2o_faces_map) B = typeof(n2o_cell_to_child_id) @@ -67,7 +74,26 @@ end but the algorithm is optimized for Vectors (refinement only). """ function get_o2n_faces_map(ncell_to_ocell::Table{T}) where {T<:Integer} - @notimplemented + nC = maximum(ncell_to_ocell.data) + + ptrs = fill(0,nC+1) + for ccell in ncell_to_ocell.data + ptrs[ccell+1] += 1 + end + Arrays.length_to_ptrs!(ptrs) + + data = Vector{Int}(undef,ptrs[end]-1) + for fcell = 1:length(ncell_to_ocell.ptrs)-1 + for j = ncell_to_ocell.ptrs[fcell]:ncell_to_ocell.ptrs[fcell+1]-1 + ccell = ncell_to_ocell.data[j] + data[ptrs[ccell]] = fcell + ptrs[ccell] += 1 + end + end + Arrays.rewind_ptrs!(ptrs) + + ocell_to_ncell = Table(data,ptrs) + return ocell_to_ncell end function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer} @@ -95,12 +121,19 @@ function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer} return ocell_to_ncell end -function get_new_cell_refinement_rules(g::AdaptivityGlue) +function get_new_cell_refinement_rules(g::AdaptivityGlue{<:RefinementGlue}) old_rrules = g.refinement_rules n2o_faces_map = g.n2o_faces_map[end] return lazy_map(Reindex(old_rrules),n2o_faces_map) end +function get_new_cell_refinement_rules(g::AdaptivityGlue{<:MixedGlue}) + old_rrules = g.refinement_rules + n2o_faces_map = g.n2o_faces_map[end] + new_idx = lazy_map(Reindex(n2o_faces_map.data),n2o_faces_map.ptrs[1:end-1]) + return lazy_map(Reindex(old_rrules), new_idx) +end + function get_old_cell_refinement_rules(g::AdaptivityGlue) return g.refinement_rules end @@ -129,3 +162,104 @@ function _reindex(data,idx::Vector) m = Reindex(data) return lazy_map(m,idx) end + + +# New to old face glues + +function get_d_to_fface_to_cface(::AdaptivityGlue,::GridTopology,::GridTopology) + @notimplemented +end + +" + For each child/fine face, returns the parent/coarse face containing it. The parent + face might have higher dimension. + + Returns two arrays: + - [dimension][fine face gid] -> coarse parent face gid + - [dimension][fine face gid] -> coarse parent face dimension +" +function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue}, + ctopo::GridTopology{Dc}, + ftopo::GridTopology{Dc}) where Dc + + # Local data for each coarse cell, at the RefinementRule level + rrules = Adaptivity.get_old_cell_refinement_rules(glue) + ccell_to_d_to_faces = lazy_map(rr->map(d->Geometry.get_faces(get_grid_topology(rr.ref_grid),Dc,d),0:Dc),rrules) + ccell_to_d_to_fface_to_parent_face = lazy_map(get_d_to_face_to_parent_face,rrules) + + # Global data, concerning the complete meshes + ccell_to_fcell = glue.o2n_faces_map + d_to_ccell_to_cface = map(d->Geometry.get_faces(ctopo,Dc,d),0:Dc) + d_to_fcell_to_fface = map(d->Geometry.get_faces(ftopo,Dc,d),0:Dc) + + d_to_fface_to_cface = [fill(Int32(0),num_faces(ftopo,d)) for d in 0:Dc] + d_to_fface_to_cface_dim = [fill(Int32(0),num_faces(ftopo,d)) for d in 0:Dc] + + # For each coarse cell + for ccell in 1:num_cells(ctopo) + local_d_to_fface_to_parent_face, + local_d_to_fface_to_parent_dim = ccell_to_d_to_fface_to_parent_face[ccell] + # For each fine subcell: + # child_id -> Local Id of the fine cell within the refinement rule (ccell) + for (child_id,fcell) in enumerate(ccell_to_fcell[ccell]) + # For each fine face on the fine subcell: + # d -> Dimension of the fine face + # iF -> Local Id of the fine face within the fine cell + # fface -> Global Id of the fine face + for d in 0:Dc + for (iF,fface) in enumerate(d_to_fcell_to_fface[d+1][fcell]) + # Local Id of the fine face within the refinement rule + fface_child_id = ccell_to_d_to_faces[ccell][d+1][child_id][iF] + # Local Id of the coarse parent face within the coarse cell + parent = local_d_to_fface_to_parent_face[d+1][fface_child_id] + + # Global Id of the coarse parent face, and it's dimension + cface_dim = local_d_to_fface_to_parent_dim[d+1][fface_child_id] + cface = d_to_ccell_to_cface[cface_dim+1][ccell][parent] + d_to_fface_to_cface[d+1][fface] = cface + d_to_fface_to_cface_dim[d+1][fface] = cface_dim + end + end + end + end + + return (d_to_fface_to_cface,d_to_fface_to_cface_dim) +end + +# FaceLabeling refinement + +function _refine_face_labeling(coarse_labeling::FaceLabeling, + glue :: AdaptivityGlue, + ctopo :: GridTopology, + ftopo :: GridTopology) + d_to_fface_to_cface, + d_to_fface_to_cface_dim = get_d_to_fface_to_cface(glue,ctopo,ftopo) + + return _refine_face_labeling(coarse_labeling,d_to_fface_to_cface,d_to_fface_to_cface_dim) +end + +function _refine_face_labeling(coarse_labeling::FaceLabeling, + d_to_fface_to_cface, + d_to_fface_to_cface_dim) + tag_to_name = copy(coarse_labeling.tag_to_name) + tag_to_entities = copy(coarse_labeling.tag_to_entities) + + Dc = num_dims(coarse_labeling) + d_to_dface_to_entity = Vector{Vector{Int32}}(undef,Dc+1) + for d in 0:Dc + nF = length(d_to_fface_to_cface[d+1]) + dface_to_entity = Vector{Int32}(undef,nF) + + for fface in 1:nF + cface = d_to_fface_to_cface[d+1][fface] + cface_dim = d_to_fface_to_cface_dim[d+1][fface] + + cface_entity = coarse_labeling.d_to_dface_to_entity[cface_dim+1][cface] + dface_to_entity[fface] = cface_entity + end + + d_to_dface_to_entity[d+1] = dface_to_entity + end + + return Geometry.FaceLabeling(d_to_dface_to_entity,tag_to_entities,tag_to_name) +end \ No newline at end of file diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl index 7347a1356..53452b995 100644 --- a/src/Adaptivity/EdgeBasedRefinement.jl +++ b/src/Adaptivity/EdgeBasedRefinement.jl @@ -86,19 +86,21 @@ function refine(method::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc, # a) nothing -> All cells get refined # b) AbstractArray{<:Bool} of size num_cells(model) # -> Only cells such that cells_to_refine[iC] == true get refined - # c) AbstractArray{<:Integer} + # 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(method, model.grid_topology,cells_to_refine) - topo = _refine_unstructured_topology(model.grid_topology,rrules,faces_list) + 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)) - labels = FaceLabeling(topo) - ref_model = UnstructuredDiscreteModel(grid,topo,labels) - ## Create ref glue + 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 @@ -415,14 +417,10 @@ end Edge-based RefinementRule where a new vertex is added to each edge of the original Polytope. """ -function RedRefinementRule(p::Polytope{2}) - @notimplementedif (p ∉ [TRI,QUAD]) +function RedRefinementRule(p::Polytope) + @notimplementedif (p ∉ [TRI,QUAD,HEX]) - if p == TRI - faces_list = ([1,2,3],[1,2,3],[]) - elseif p == QUAD - faces_list = ([1,2,3,4],[1,2,3,4],[1]) - end + faces_list = _get_red_refined_faces_list(p) coords = get_new_coordinates_from_faces(p,faces_list) polys, cell_types, conn = _get_red_refined_connectivity(p) @@ -432,24 +430,48 @@ function RedRefinementRule(p::Polytope{2}) return RefinementRule(RedRefinement(),p,ref_grid) end -function _get_red_refined_connectivity(p::Polytope{2}) +function _get_red_refined_faces_list(p::Polytope) + if p == TRI + return (Int32[1,2,3],Int32[1,2,3],Int32[]) + elseif p == QUAD + return (Int32[1,2,3,4],Int32[1,2,3,4],Int32[1]) + elseif p == HEX + return (Int32[1,2,3,4,5,6,7,8],Int32[1,2,3,4,5,6,7,8,9,10,11,12],Int32[1,2,3,4,5,6],Int32[1]) + end + @notimplemented +end + +function _get_red_refined_connectivity(p::Polytope) if p == TRI polys = [TRI] - cell_type = [1,1,1,1] - conn_data = [1,4,5, - 2,4,6, - 3,5,6, - 4,5,6] - conn_ptrs = [1,4,7,10,13] + cell_type = Int32[1,1,1,1] + conn_data = Int32[1,4,5, + 2,4,6, + 3,5,6, + 4,5,6] + conn_ptrs = Int32[1,4,7,10,13] return polys, cell_type, Table(conn_data,conn_ptrs) elseif p == QUAD polys = [QUAD] - cell_type = [1,1,1,1] - conn_data = [1,5,7,9, - 5,2,9,8, - 7,9,3,6, - 9,8,6,4] - conn_ptrs = [1,5,9,13,17] + cell_type = Int32[1,1,1,1] + conn_data = Int32[1,5,7,9, + 5,2,9,8, + 7,9,3,6, + 9,8,6,4] + conn_ptrs = Int32[1,5,9,13,17] + return polys, cell_type, Table(conn_data,conn_ptrs) + elseif p == HEX + polys = [HEX] + cell_type = Int32[1,1,1,1,1,1,1,1] + conn_data = Int32[ 1, 9,13,21,17,23,25,27, + 9, 2,21,14,23,18,27,26, + 13,21, 3,10,25,27,19,24, + 21,14,10, 4,27,26,24,20, + 17,23,25,27, 5,11,15,22, + 23,18,27,26,11, 6,22,16, + 25,27,19,24,15,22, 7,12, + 27,26,24,20,22,16,12, 8] + conn_ptrs = Int32[1,9,17,25,33,41,49,57,65] return polys, cell_type, Table(conn_data,conn_ptrs) end @notimplemented @@ -457,12 +479,117 @@ end function _has_interior_point(rr::RefinementRule,::RedRefinement) p = get_polytope(rr) - if p == QUAD + if p ∈ [QUAD,HEX] return true end return false end +# [Face dimension][Coarse Face id] -> [Fine faces] +function get_d_to_face_to_child_faces(rr::RefinementRule,::RedRefinement) + p = get_polytope(rr) + if p == QUAD + return [ + [Int32[1],Int32[2],Int32[3],Int32[4]], # [Coarse Node] -> [Fine Node] + [Int32[1,5],Int32[8,11],Int32[3,9],Int32[7,12]], # [Coarse Edge] -> [Fine Edge] + [Int32[1,2,3,4]] # [Coarse Cell] -> [Fine Cells] + ] + elseif p == TRI + return [ + [Int32[1],Int32[2],Int32[3]], # [Coarse Node] -> [Fine Node] + [Int32[1,4],Int32[2,7],Int32[5,8]], # [Coarse Edge] -> [Fine Edge] + [Int32[1,2,3]] # [Coarse Cell] -> [Fine Cells] + ] + elseif p == HEX + return [ + [Int32[1],Int32[2],Int32[3],Int32[4], + Int32[5],Int32[6],Int32[7],Int32[8]], # [Coarse Node] -> [Fine Node] + [Int32[1,13],Int32[21,29],Int32[34,42],Int32[47,52], + Int32[5,23],Int32[17,31],Int32[36,48],Int32[44,53], + Int32[9,38],Int32[19,45],Int32[27,50],Int32[33,54]], # [Coarse Edge] -> [Fine Edge] + [Int32[1,7,12,17] ,Int32[21,26,30,34], + Int32[3,9,22,27] ,Int32[14,19,31,35], + Int32[5,15,24,32],Int32[11,20,29,36]], # [Coarse Face] -> [Fine Face] + [Int32[1,2,3,4,5,6,7,8]], # [Coarse Cell] -> [Fine Cells] + ] + else + @notimplemented + end +end + +# 1 - [Face dimension][Fine Face id] -> [Parent Face] +# 2 - [Face dimension][Fine Face id] -> [Parent Face Dimension] +function get_d_to_face_to_parent_face(rr::RefinementRule,::RedRefinement) + p = get_polytope(rr) + if p == QUAD + parent_faces = [ + Int32[1,2,3,4,1,2,3,4,1], # [Fine node] -> [Coarse face] + Int32[1,1,3,1,1,1,4,2,3,1,2,4], # [Fine edge] -> [Coarse face] + Int32[1,1,1,1] # [Fine cell] -> [Coarse face] + ] + parent_dims = [ + Int32[0,0,0,0,1,1,1,1,2], # [Fine node] -> [Coarse face dim] + Int32[1,2,1,2,1,2,1,1,1,2,1,1], # [Fine edge] -> [Coarse face dim] + Int32[2,2,2,2] # [Fine cell] -> [Coarse face dim] + ] + elseif p == TRI + parent_faces = [ + Int32[1,2,3,1,2,3], # [Fine node] -> [Coarse face] + Int32[1,2,1,1,3,1,2,3,1], # [Fine edge] -> [Coarse face] + Int32[1,1,1,1] # [Fine cell] -> [Coarse face] + ] + parent_dims = [ + Int32[0,0,0,1,1,1], # [Fine node] -> [Coarse face dim] + Int32[1,1,2,1,1,2,1,1,2], # [Fine edge] -> [Coarse face dim] + Int32[2,2,2,2] # [Fine cell] -> [Coarse face dim] + ] + elseif p == HEX + parent_faces = [ + Int32[1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,1], # [Fine node] -> [Coarse face] + Int32[1,1,3,1,5,1,5,1,9,3,5,1, + 1,1,3,1,6,6,10,6, + 2,4,5,1,5,1,11,4, + 2,4,6,6,12, + 3,2,7,2,9,3,5,1, + 3,2,8,10,6, + 4,7,2,11,4, + 4,8,12], # [Fine edge] -> [Coarse face] + Int32[1,1,3,1,5,1, + 1,1,3,1,6, + 1,1,4,5,1, + 1,1,4,6, + 2,3,1,5,1, + 2,3,1,6, + 2,4,5,1, + 2,4,6], # [Fine face] -> [Coarse face] + Int32[1,1,1,1,1,1,1,1] # [Fine cell] -> [Coarse face] + ] + parent_dims = [ + Int32[0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3], # [Fine node] -> [Coarse face dim] + Int32[1,2,2,3,1,2,2,3,1,2,2,3, + 1,2,2,3,1,2,1,2, + 1,2,1,2,2,3,1,2, + 1,2,1,2,1, + 1,2,1,2,1,2,2,3, + 1,2,1,1,2, + 1,1,2,1,2, + 1,1,1], # [Fine edge] -> [Coarse face dim] + Int32[2,3,2,3,2,3, + 2,3,2,3,2, + 2,3,2,2,3, + 2,3,2,2, + 2,2,3,2,3, + 2,2,3,2, + 2,2,2,3, + 2,2,2], # [Fine face] -> [Coarse face dim] + Int32[3,3,3,3,3,3,3,3] # [Fine cell] -> [Coarse face dim] + ] + else + parent_faces, parent_dims = get_d_to_face_to_parent_face(rr,GenericRefinement()) + end + return parent_faces, parent_dims +end + """ Edge-based RefinementRule where a new vertex is added to a single edge of the original Polytope. diff --git a/src/Adaptivity/OldToNewFields.jl b/src/Adaptivity/OldToNewFields.jl new file mode 100644 index 000000000..02b468742 --- /dev/null +++ b/src/Adaptivity/OldToNewFields.jl @@ -0,0 +1,55 @@ + +abstract type NewFieldType end; +struct CoarsenedNewFieldType <: NewFieldType end; +struct RefinedOrUntouchedNewFieldType <: NewFieldType end; + +# Unfortunately, I cannot use traits in OldToNewField as its +# usage results in conversion errors when leveraging it with +# the lazy_map infrastructure +struct OldToNewField <: Fields.Field + new_field_type::NewFieldType + fine_to_coarse_field + refined_or_untouched_field + function OldToNewField(a::NewFieldType,b::Fields.Field,c::Fields.Field) + new(a,b,c) + end +end + +function OldToNewField(old_fields::AbstractArray{<:Fields.Field}, + rrule::RefinementRule, + child_ids::AbstractVector{<:Integer}) + @assert length(old_fields)==length(child_ids) + if length(old_fields) == 1 + cell_map = get_cell_map(rrule)[child_ids[1]] + old_field = old_fields[1] + fine_to_coarse_field = FineToCoarseField( + [old_field for i = 1:num_subcells(rrule)], + rrule, + [i for i=1:num_subcells(rrule)]) + refined_or_untouched_field = old_field∘cell_map + return OldToNewField(RefinedOrUntouchedNewFieldType(),fine_to_coarse_field,refined_or_untouched_field) + else + @assert length(old_fields) <= num_subcells(rrule) + fine_to_coarse_field = FineToCoarseField(old_fields,rrule,child_ids) + cell_map = get_cell_map(rrule)[1] + refined_or_untouched_field = old_fields[1]∘cell_map + return OldToNewField(CoarsenedNewFieldType(),fine_to_coarse_field,refined_or_untouched_field) + end +end + +function Fields.return_cache(a::OldToNewField,x::AbstractArray{<:Point}) + f2c_cache = Fields.return_cache(a.fine_to_coarse_field,x) + rou_cache = Fields.return_cache(a.refined_or_untouched_field,x) + return (f2c_cache,rou_cache) +end + +function Fields.evaluate!(cache,a::OldToNewField,x::AbstractArray{<:Point}) + if isa(a.new_field_type,CoarsenedNewFieldType) + f2c_cache,rou_cache = cache + Fields.evaluate!(f2c_cache,a.fine_to_coarse_field,x) + else + @assert isa(a.new_field_type,RefinedOrUntouchedNewFieldType) + f2c_cache,rou_cache = cache + Fields.evaluate!(rou_cache,a.refined_or_untouched_field,x) + end +end diff --git a/src/Adaptivity/RefinementRules.jl b/src/Adaptivity/RefinementRules.jl index 5acc1df62..183584eae 100644 --- a/src/Adaptivity/RefinementRules.jl +++ b/src/Adaptivity/RefinementRules.jl @@ -104,6 +104,224 @@ function bundle_points_by_subcell(rr::RefinementRule,x::AbstractArray{<:Point}) return Table(data,ptrs) end + +# Faces to child faces, dof maps + +" +Given a `RefinementRule`, returns for each parent/coarse face the child/fine faces of the +same dimension that it contains. Therefore, only fine faces at the coarse cell boundary are +listed in the returned structure. + +Returns: [Face dimension][Coarse Face id] -> [Fine faces] +" +function get_d_to_face_to_child_faces(rr::RefinementRule) + get_d_to_face_to_child_faces(rr,RefinementRuleType(rr)) +end + +# Generic version of the function. Spetializations may exist for some other ref rule types. +# This generic method relies on get_d_to_face_to_parent_face, and simply inverts the map. +function get_d_to_face_to_child_faces(rr::RefinementRule,::RefinementRuleType) + d_to_face_to_parent_face, d_to_face_to_parent_face_dim = get_d_to_face_to_parent_face(rr) + + poly = get_polytope(rr) + Dc = num_cell_dims(poly) + d_to_face_to_child_faces = Vector{Vector{Vector{Int32}}}(undef,Dc+1) + + for cface_dim in 0:Dc + num_cfaces = num_faces(poly,cface_dim) + cface_to_child_faces = Vector{Vector{Int32}}(undef,num_cfaces) + + parent_faces = d_to_face_to_parent_face[cface_dim+1] + parent_faces_dim = d_to_face_to_parent_face_dim[cface_dim+1] + parent_pairs = collect(zip(parent_faces,parent_faces_dim)) + + for cface in 1:num_cfaces + cface_to_child_faces[cface] = findall(p -> (p[1] == cface) && (p[2] == cface_dim),parent_pairs) + end + + d_to_face_to_child_faces[cface_dim+1] = cface_to_child_faces + end + + return d_to_face_to_child_faces +end + +""" +Given a `RefinementRule`, returns for each fine/child face the parent/coarse face +containing it. The parent face can have higher dimension. + +Returns the tuple (A,B) with + - A = [Face dimension][Fine Face id] -> [Parent Face] + - B = [Face dimension][Fine Face id] -> [Parent Face Dimension] +""" +function get_d_to_face_to_parent_face(rr::RefinementRule) + get_d_to_face_to_parent_face(rr,RefinementRuleType(rr)) +end + +# Generic version of the function. Spetializations may exist for some other ref rule types. +function get_d_to_face_to_parent_face(rr::RefinementRule,::RefinementRuleType) + # WARNING: The functions below are NOT general for any point and any polytope. + # They are only valid for the specific case of a refinement rule. + function belongs_to_face(::Val{0},::Val{0},fface_coords,cface_coords) + return fface_coords[1] == cface_coords[1] + end + function belongs_to_face(::Val{0},::Val{1},fface_coords,cface_coords) + norm(cross(fface_coords[1] - cface_coords[1], fface_coords[1] - cface_coords[2])) ≈ 0.0 + end + function belongs_to_face(::Val{0},::Val{2},fface_coords,cface_coords) + n = cross(cface_coords[2] - cface_coords[1], cface_coords[3] - cface_coords[1]) + return sum(map(ccoords -> dot(n,ccoords - fface_coords[1]), cface_coords)) ≈ 0.0 + end + function belongs_to_face(::Val{fface_dim},::Val{cface_dim},fface_coords,cface_coords) where {fface_dim,cface_dim} + return all(map(p -> belongs_to_face(Val(0),Val(cface_dim),[p],cface_coords),fface_coords)) + end + + ref_grid = get_ref_grid(rr) + topo = get_grid_topology(ref_grid) + poly = get_polytope(rr) + + fnode_coords = get_node_coordinates(ref_grid) + cnode_coords = get_vertex_coordinates(poly) + + Dc = num_cell_dims(ref_grid) + d_to_face_to_parent_face = Vector{Vector{Int32}}(undef,Dc+1) + d_to_face_to_parent_face_dim = Vector{Vector{Int32}}(undef,Dc+1) + + # For each fface dimension + for fface_dim in 0:Dc + fface_nodes = Geometry.get_faces(topo,fface_dim,0) + fface_node_coords = lazy_map(nodes -> lazy_map(Reindex(fnode_coords),nodes),fface_nodes) + + num_ffaces = Geometry.num_faces(topo,fface_dim) + fface_to_parent_face = fill(Int32(-1),num_ffaces) + fface_to_parent_face_dim = fill(Int32(-1),num_ffaces) + + # For each fface find the parent face containing it + for (fface,fcoords) in enumerate(fface_node_coords) + found = false + cface_dim = fface_dim + # Start with cfaces of the same dimension as the fface, and go up until reaching Dc-1 + while (!found) && (cface_dim < Dc) + cface_nodes = get_faces(poly,cface_dim,0) + cface_node_coords = lazy_map(nodes -> lazy_map(Reindex(cnode_coords),nodes),cface_nodes) + + for (cface,ccoords) in enumerate(cface_node_coords) + if !found && belongs_to_face(Val(fface_dim),Val(cface_dim),fcoords,ccoords) + found = true + fface_to_parent_face[fface] = cface + fface_to_parent_face_dim[fface] = cface_dim + end + end + cface_dim += 1 + end + if !found # Belongs to the cell itself (dimension Dc) + fface_to_parent_face[fface] = 1 + fface_to_parent_face_dim[fface] = Dc + end + end + + d_to_face_to_parent_face[fface_dim+1] = fface_to_parent_face + d_to_face_to_parent_face_dim[fface_dim+1] = fface_to_parent_face_dim + end + + return d_to_face_to_parent_face, d_to_face_to_parent_face_dim +end + +function _get_terms(poly::Polytope,orders) + _nodes, facenodes = ReferenceFEs._compute_nodes(poly,orders) + terms = ReferenceFEs._coords_to_terms(_nodes,orders) + return terms +end + +function _get_face_orders(p::Polytope{Dc},D::Int,orders::Tuple) where Dc + @check length(orders) == Dc + @check 1 <= D < Dc + @check is_n_cube(p) + + if D == 1 # Edges (2D, 3D) + tangents = get_edge_tangent(p) + face_orders = map(tangents) do t + axis = findfirst(i -> abs(t[i]) > 0.5 ,1:Dc) + return [orders[axis]] + end + elseif D == Dc-1 # Faces (3D) + normals = get_facet_normal(p) + face_orders = map(normals) do n + mask = map(i -> abs(n[i]) < 1.e-3,1:Dc) + return [orders[mask]...] + end + else + @notimplemented + end + + return face_orders +end + +function _get_local_dof_ranges(p::Polytope{Dc},orders) where Dc + @check length(orders) == Dc + @check is_n_cube(p) + idx = CartesianIndices(Tuple(fill(2,Dc))) + ranges = map(idx) do ii + map(Tuple(ii),orders) do i,o + (i-1)*o+1:i*o+1 + end + end + return ranges +end + +""" +Given a `RefinementRule` of dimension Dc and a Dc-Tuple `fine_orders` of approximation orders, +returns a map between the fine nodal dofs of order `fine_orders` in the reference grid and the +coarse nodal dofs of order `2⋅fine_orders` in the coarse parent cell. + +The result is given for each coarse/parent face of dimension `D` as a list of the corresponding +fine dof lids, i.e +- [coarse face][coarse dof lid] -> fine dof lid +""" +function get_face_subface_ldof_to_cell_ldof(rr::RefinementRule{<:ExtrusionPolytope{Dc}}, + fine_orders::NTuple{Dc,<:Integer}, + D::Int) where Dc + poly = get_polytope(rr) + coarse_orders = 2 .* fine_orders + coarse_reffe = ReferenceFE(poly,lagrangian,Float64,coarse_orders) + coarse_face_polys = CompressedArray(ReferenceFEs._compute_reffaces_and_face_types(poly,Val(D))...) + c_edge_to_coarse_dof = coarse_reffe.face_nodes[get_dimranges(poly)[D+1]] + + model = get_ref_grid(rr) + fine_face_grid = Grid(ReferenceFE{D},model) + fine_face_polys = CompressedArray(map(get_polytope,get_reffes(fine_face_grid)),get_cell_type(fine_face_grid)) + + d_to_face_to_child_faces = get_d_to_face_to_child_faces(rr) + face_to_child_faces = d_to_face_to_child_faces[D+1] + + coarse_face_orders = _get_face_orders(poly,D,coarse_orders) + fine_face_orders = _get_face_orders(poly,D,fine_orders) + + num_coarse_faces = num_faces(coarse_reffe,D) + coarse_dofs_above_fine_dofs = Vector{Vector{Vector{Int32}}}(undef,num_coarse_faces) + for cF in 1:num_coarse_faces + coarse_face_poly = coarse_face_polys[cF] + coarse_terms = _get_terms(coarse_face_poly,coarse_face_orders[cF]) + coarse_dofs = zeros(Int32,Tuple(maximum(coarse_terms))) + coarse_dofs[coarse_terms] .= c_edge_to_coarse_dof[cF] + fine_face_to_dof_range = _get_local_dof_ranges(coarse_face_poly,fine_face_orders[cF]) + + child_faces = face_to_child_faces[cF] + fine_dofs = Vector{Vector{Int32}}(undef,length(child_faces)) + for (i,fF) in enumerate(child_faces) + fine_face_poly = fine_face_polys[fF] + fine_terms = _get_terms(fine_face_poly,fine_face_orders[cF]) + + local_dof_range = fine_face_to_dof_range[i] + local_coarse_dofs = view(coarse_dofs,local_dof_range...) + fine_dofs[i] = map(Reindex(local_coarse_dofs),fine_terms) + end + coarse_dofs_above_fine_dofs[cF] = fine_dofs + end + + return coarse_dofs_above_fine_dofs +end + + # GenericRefinement Rule function RefinementRule(reffe::LagrangianRefFE{D},nrefs::Integer;kwargs...) where D diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl index 68c50ab38..ed9e7dcd5 100644 --- a/src/Geometry/CartesianDiscreteModels.jl +++ b/src/Geometry/CartesianDiscreteModels.jl @@ -8,65 +8,62 @@ struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D} grid::CartesianGrid{D,T,F} grid_topology::UnstructuredGridTopology{D,D,T,Oriented} face_labeling::FaceLabeling - @doc """ - CartesianDiscreteModel(desc::CartesianDescriptor) - - Inner constructor - """ - function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}) where {D,T,F} - grid = CartesianGrid(desc) - _grid = UnstructuredGrid(grid) - if any(desc.isperiodic) - topo = _cartesian_grid_topology_with_periodic_bcs(_grid, desc.isperiodic, desc.partition) - else - topo = UnstructuredGridTopology(_grid) - end - nfaces = [num_faces(topo,d) for d in 0:num_cell_dims(topo)] - labels = FaceLabeling(nfaces) - _fill_cartesian_face_labeling!(labels,topo) - new{D,T,F}(grid,topo,labels) +end + +""" + CartesianDiscreteModel(desc::CartesianDescriptor) +""" +function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}) where {D,T,F} + grid = CartesianGrid(desc) + _grid = UnstructuredGrid(grid) + if any(desc.isperiodic) + topo = _cartesian_grid_topology_with_periodic_bcs(_grid, desc.isperiodic, desc.partition) + else + topo = UnstructuredGridTopology(_grid) end + nfaces = [num_faces(topo,d) for d in 0:num_cell_dims(topo)] + labels = FaceLabeling(nfaces) + _fill_cartesian_face_labeling!(labels,topo) + return CartesianDiscreteModel(grid,topo,labels) +end + +""" + CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}, + cmin::CartesianIndex, + cmax::CartesianIndex) + + Builds a CartesianDiscreteModel object which represents a subgrid of + a (larger) grid represented by desc. This subgrid is described by its + D-dimensional minimum (cmin) and maximum (cmax) CartesianIndex + identifiers. +""" +function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}, + cmin::CartesianIndex, + cmax::CartesianIndex, + remove_boundary=map(i->false,desc.sizes)) where {D,T,F} + + suborigin = Tuple(desc.origin) .+ (Tuple(cmin) .- 1) .* desc.sizes + subpartition = Tuple(cmax) .- Tuple(cmin) .+ 1 + subsizes = desc.sizes + subisperiodic = map(subpartition,desc.partition,desc.isperiodic) do subn,n,p + # Periodic only make sense in dims in which we take the full span. + # Otherwise, we would artificially modify the periodicity period. + subn==n ? p : false + end + subdesc = + CartesianDescriptor(Point(suborigin), subsizes, subpartition; map=desc.map, isperiodic=subisperiodic) - @doc """ - CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}, - cmin::CartesianIndex, - cmax::CartesianIndex) - - Builds a CartesianDiscreteModel object which represents a subgrid of - a (larger) grid represented by desc. This subgrid is described by its - D-dimensional minimum (cmin) and maximum (cmax) CartesianIndex - identifiers. - - Inner constructor - """ - function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}, - cmin::CartesianIndex, - cmax::CartesianIndex, - remove_boundary=map(i->false,desc.sizes)) where {D,T,F} - - suborigin = Tuple(desc.origin) .+ (Tuple(cmin) .- 1) .* desc.sizes - subpartition = Tuple(cmax) .- Tuple(cmin) .+ 1 - subsizes = desc.sizes - subisperiodic = map(subpartition,desc.partition,desc.isperiodic) do subn,n,p - # Periodic only make sense in dims in which we take the full span. - # Otherwise, we would artificially modify the periodicity period. - subn==n ? p : false - end - subdesc = - CartesianDescriptor(Point(suborigin), subsizes, subpartition; map=desc.map, isperiodic=subisperiodic) - - grid = CartesianGrid(subdesc) - _grid = UnstructuredGrid(grid) - if any(subdesc.isperiodic) - topo = _cartesian_grid_topology_with_periodic_bcs(_grid, subdesc.isperiodic, subdesc.partition) - else - topo = UnstructuredGridTopology(_grid) - end - nfaces = [num_faces(topo, d) for d = 0:num_cell_dims(topo)] - labels = FaceLabeling(nfaces) - _fill_subgrid_cartesian_face_labeling!(labels,topo,subdesc,desc,cmin,remove_boundary) - new{D,T,F}(grid, topo, labels) + grid = CartesianGrid(subdesc) + _grid = UnstructuredGrid(grid) + if any(subdesc.isperiodic) + topo = _cartesian_grid_topology_with_periodic_bcs(_grid, subdesc.isperiodic, subdesc.partition) + else + topo = UnstructuredGridTopology(_grid) end + nfaces = [num_faces(topo, d) for d = 0:num_cell_dims(topo)] + labels = FaceLabeling(nfaces) + _fill_subgrid_cartesian_face_labeling!(labels,topo,subdesc,desc,cmin,remove_boundary) + return CartesianDiscreteModel(grid, topo, labels) end """ diff --git a/src/ReferenceFEs/ExtrusionPolytopes.jl b/src/ReferenceFEs/ExtrusionPolytopes.jl index 37025d5b0..fed341ac3 100644 --- a/src/ReferenceFEs/ExtrusionPolytopes.jl +++ b/src/ReferenceFEs/ExtrusionPolytopes.jl @@ -723,6 +723,9 @@ function _nullspace(v) end end +# Returns the first vertex not belonging to the facet i_f, or -1 if all vertices +# belong to the facet. +# nf_vs is the array of arrays of vertices of the facets of the polytope. function _vertex_not_in_facet(p::DFace, i_f, nf_vs) for i in p.nf_dimranges[end][1] is_in_f = false diff --git a/test/AdaptivityTests/FaceLabelingTests.jl b/test/AdaptivityTests/FaceLabelingTests.jl new file mode 100644 index 000000000..bf5ef97b6 --- /dev/null +++ b/test/AdaptivityTests/FaceLabelingTests.jl @@ -0,0 +1,46 @@ +module FaceLabelingTests + +using Test +using Gridap +using Gridap.Geometry +using Gridap.Adaptivity + +Dc = 2 +model1 = CartesianDiscreteModel((0,1,0,1),(2,2)) +model2 = UnstructuredDiscreteModel(model1) + +for model in [model1]#,model2] + coarse_labeling = get_face_labeling(model) + add_tag_from_tags!(coarse_labeling,"left_boundary",[7]) + add_tag_from_tags!(coarse_labeling,"bottom_corners",[1,2]) +end + +fmodel1 = refine(model1,3) +fmodel2 = refine(model2) + +for fmodel in [fmodel1,fmodel2] + ftopo = get_grid_topology(fmodel) + fine_labeling = get_face_labeling(fmodel) + + # Compute face centroids + vertex_coordinates = Geometry.get_node_coordinates(fmodel) + d_to_face_to_vertex = map(d->Geometry.get_faces(ftopo,d,0),0:Dc) + + d_to_centroid = map(d_to_face_to_vertex) do face_to_vertex + map(face_to_vertex) do vertices + sum(vertex_coordinates[vertices])/length(vertices) + end + end + + # Check tags + left_boundary_by_tags = findall(get_face_mask(fine_labeling,"left_boundary",1)) + left_boundary_by_coords = findall(p->abs(p[1]) < 1.e-3,d_to_centroid[2]) + @test sort(left_boundary_by_tags) == sort(left_boundary_by_coords) + + bottom_corners_by_tags = findall(get_face_mask(fine_labeling,"bottom_corners",0)) + bottom_corners_by_coords = findall(p->(abs(p[2]) < 1.e-3) && ((abs(p[1]) < 1.e-3) || (abs(p[1]-1.0) < 1.e-3)),d_to_centroid[1]) + @test sort(bottom_corners_by_tags) == sort(bottom_corners_by_coords) +end + + +end \ No newline at end of file diff --git a/test/AdaptivityTests/RefinementRuleBoundaryTests.jl b/test/AdaptivityTests/RefinementRuleBoundaryTests.jl new file mode 100644 index 000000000..a30d6dcc2 --- /dev/null +++ b/test/AdaptivityTests/RefinementRuleBoundaryTests.jl @@ -0,0 +1,20 @@ +module RefinementRuleBoundaryTests + +using Test +using Gridap +using Gridap.Helpers +using Gridap.Adaptivity +using Gridap.Geometry +using Gridap.ReferenceFEs +using Gridap.Arrays + +# Test the final function in 2D +rr = Gridap.Adaptivity.RedRefinementRule(QUAD) +res = Adaptivity.get_face_subface_ldof_to_cell_ldof(rr,(2,2),1) + +# Test the final function in 3D +rr = Gridap.Adaptivity.RedRefinementRule(HEX) +res = Adaptivity.get_face_subface_ldof_to_cell_ldof(rr,(2,2,2),1) +res = Adaptivity.get_face_subface_ldof_to_cell_ldof(rr,(2,2,2),2) + +end \ No newline at end of file