extended subgrid functionality (#49)

* extended subgrid method (new argument 'support' can be ON_CELLS, ON_BFACES, ON_FACES) * SubGrid type has now a parameter that knows the support relative to parent grid, removed BoundarySubGrid type * removed get_facegrid, get_bfacegrid, get_edgegrid in derived.jl (that offer no new functionality compared to subgrid) * fixed default behaviour for coordinatesystem and restricted 1D coordinate sorting to project = true
WIAS-PDELib · Jun 20, 2024 · f3ff8f0 · f3ff8f0
1 parent 4729632
commit f3ff8f0
Show file tree

Hide file tree

Showing 4 changed files with 46 additions and 68 deletions.
diff --git a/src/ExtendableGrids.jl b/src/ExtendableGrids.jl
@@ -85,11 +85,22 @@ export num_pcolors, num_partitions, pcolors, pcolor_partitions, partition_cells,
 export partition, num_partitions_per_color,check_partitioning
 export AbstractPartitioningAlgorithm, TrivialPartitioning, PlainMetisPartitioning
 
+include("assemblytypes.jl")
+export AssemblyType
+export AT_NODES, ON_CELLS, ON_FACES, ON_IFACES, ON_BFACES, ON_EDGES, ON_BEDGES
+export ItemType4AssemblyType
+export GridComponentNodes4AssemblyType
+export GridComponentVolumes4AssemblyType
+export GridComponentGeometries4AssemblyType
+export GridComponentRegions4AssemblyType
+export GridComponentUniqueGeometries4AssemblyType
+export GridComponentAssemblyGroups4AssemblyType
+
 include("subgrid.jl")
 export subgrid
 export ParentGrid
 export NodeParents, CellParents, FaceParents, BFaceParents
-export ParentGridRelation, SubGrid, BoundarySubGrid, RefinedGrid
+export ParentGridRelation, SubGrid, RefinedGrid
 
 include("shape_specs.jl")
 export refcoords_for_geometry
@@ -118,7 +129,6 @@ export CellAssemblyGroups, FaceAssemblyGroups, BFaceAssemblyGroups, EdgeAssembly
 export GridComponent4TypeProperty
 export ITEMTYPE_CELL, ITEMTYPE_FACE, ITEMTYPE_BFACE, ITEMTYPE_EDGE, ITEMTYPE_BEDGE
 export PROPERTY_NODES, PROPERTY_REGION, PROPERTY_VOLUME, PROPERTY_UNIQUEGEOMETRY, PROPERTY_GEOMETRY, PROPERTY_ASSEMBLYGROUP
-export get_facegrid, get_bfacegrid, get_edgegrid
 export GridEGTypes
 export GridRegionTypes
 
@@ -138,17 +148,6 @@ include("adaptive_meshrefinements.jl")
 export bulk_mark
 export RGB_refine
 
-include("assemblytypes.jl")
-export AssemblyType
-export AT_NODES, ON_CELLS, ON_FACES, ON_IFACES, ON_BFACES, ON_EDGES, ON_BEDGES
-export ItemType4AssemblyType
-export GridComponentNodes4AssemblyType
-export GridComponentVolumes4AssemblyType
-export GridComponentGeometries4AssemblyType
-export GridComponentRegions4AssemblyType
-export GridComponentUniqueGeometries4AssemblyType
-export GridComponentAssemblyGroups4AssemblyType
-
 include("l2gtransformations.jl")
 export L2GTransformer, update_trafo!, eval_trafo!, mapderiv!
 

diff --git a/src/derived.jl b/src/derived.jl
@@ -125,42 +125,6 @@ GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_GEOMETRY}) = B
 GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_UNIQUEGEOMETRY}) = UniqueBEdgeGeometries
 GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_ASSEMBLYGROUP}) = BEdgeAssemblyGroups
 
-
-function get_facegrid(source_grid::ExtendableGrid{Tc,Ti}) where {Tc,Ti}
-    facegrid=ExtendableGrid{Tc,Ti}()
-    facegrid[Coordinates]=source_grid[Coordinates]
-    facegrid[CellNodes]=source_grid[FaceNodes]
-    facegrid[CoordinateSystem]=source_grid[CoordinateSystem]
-    facegrid[CellGeometries]=source_grid[FaceGeometries]
-    facegrid[UniqueCellGeometries]=source_grid[UniqueFaceGeometries]
-    facegrid[CellRegions] = source_grid[FaceRegions]
-    # todo: facegrid[CellFaces] = source_grid[FaceEdges]
-    return facegrid
-end
-
-function get_bfacegrid(source_grid::ExtendableGrid{Tc,Ti})  where {Tc,Ti}
-    bfacegrid=ExtendableGrid{Tc,Ti}()
-    bfacegrid[Coordinates]=source_grid[Coordinates]
-    bfacegrid[CellNodes]=source_grid[BFaceNodes]
-    bfacegrid[CoordinateSystem]=source_grid[CoordinateSystem]
-    bfacegrid[CellGeometries]=source_grid[BFaceGeometries]
-    bfacegrid[UniqueCellGeometries]=source_grid[UniqueBFaceGeometries]
-    bfacegrid[CellRegions] = source_grid[BFaceRegions]
-    # todo: bfacegrid[CellFaces] = source_grid[BFaceEdges] (or sub-view of FaceEdges ?)
-    return bfacegrid
-end
-
-function get_edgegrid(source_grid::ExtendableGrid{Tc,Ti}) where {Tc,Ti}
-    edgegrid=ExtendableGrid{Tc,Ti}()
-    edgegrid[Coordinates]=source_grid[Coordinates]
-    edgegrid[CellNodes]=source_grid[EdgeNodes]
-    edgegrid[CoordinateSystem]=source_grid[CoordinateSystem]
-    edgegrid[CellGeometries]=source_grid[EdgeGeometries]
-    edgegrid[UniqueCellGeometries]=source_grid[UniqueedgeGeometries]
-    edgegrid[CellRegions] = source_grid[EdgeRegions]
-    return edgegrid
-end
-
 # show function for ExtendableGrids and defined Components in its Dict
 function showmore(io::IO, xgrid::ExtendableGrid{Tc,Ti}) where {Tc,Ti}
 
@@ -198,7 +162,7 @@ end
 function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceNodes}) where {Tc,Ti}
 
     if haskey(xgrid, ParentGrid) && haskey(xgrid, ParentGridRelation)
-        if xgrid[ParentGridRelation] === SubGrid
+        if xgrid[ParentGridRelation] <: SubGrid{ON_CELLS}
             ## get FaceNodes from ParentGrid to keep ordering and orientation
             pgrid = xgrid[ParentGrid]
             pnodes = xgrid[NodeParents]

diff --git a/src/subgrid.jl b/src/subgrid.jl
@@ -46,21 +46,26 @@ $(TYPEDEF)
 
 Grid component key type for indicating that grid is a subgrid of the parentgrid
 """
-abstract type SubGrid <: ParentGridRelation end
+abstract type SubGrid{support} <: ParentGridRelation where {support <: AssemblyType} end
 
 """
 $(TYPEDEF)
 
-Grid component key type for indicating that grid is a boundary subgrid of the parentgrid
+Grid component key type for indicating that grid is a refinement of the parentgrid
 """
-abstract type BoundarySubGrid <: ParentGridRelation end
+abstract type RefinedGrid <: ParentGridRelation end
+
 
 """
-$(TYPEDEF)
+$(TYPEDSIGNATURES)
 
-Grid component key type for indicating that grid is a refinement of the parentgrid
+Default transform for subgrid creation
 """
-abstract type RefinedGrid <: ParentGridRelation end
+function _copytransform!(a::AbstractArray,b::AbstractArray)
+    for i=1:length(a)
+        a[i]=b[i]
+    end
+end
 
 struct XIPair{Tv, Ti}
     x::Tv
@@ -71,7 +76,6 @@ end
 Base.isless(x::XIPair, y::XIPair) = (x.x < y.x)
 
 
-
 """
     subgrid(parent,                                                             
             subregions::AbstractArray;                                          
@@ -84,7 +88,8 @@ Create subgrid from list of regions.
 
 - `parent`: parent grid 
 - `subregions`:  Array of subregions which define the subgrid
-- `boundary`: if true, create codimension 1 subgrid from boundary regions.
+- 'support': support of subgrid, default is ON_CELLS but can be also ON_FACES or ON_BFACES to create codimension 1 subgrid from face/bfaces region
+- `boundary`: if true, create codimension 1 subgrid from boundary regions (same as support = ON_BFACES)
 - `transform` (kw parameter): transformation function between
    grid and subgrid coordinates acting on one point.
 - `coordinatesystem`: if `boundary==true`, specify coordinate system for the boundary.
@@ -100,8 +105,15 @@ function subgrid(parent,
                  subregions::AbstractArray;
                  transform::T=function(a,b) @views a.=b[1:length(a)] end,                                      
                  boundary=false,
-                 coordinatesystem=codim1_coordinatesystem(parent[CoordinateSystem]),
-                 project=true) where T
+                 support=ON_CELLS,
+                 project=true,
+                 coordinatesystem=project ? codim1_coordinatesystem(parent[CoordinateSystem]) : parent[CoordinateSystem]) where T
+
+    @assert support in [ON_CELLS, ON_FACES, ON_BFACES] "value ($support) for 'support' is not allowed"
+
+    if boundary
+        support = ON_BFACES
+    end
 
     Tc=coord_type(parent)
     Ti=index_type(parent)
@@ -118,14 +130,17 @@ function subgrid(parent,
         return false
     end
 
-
-    if boundary
+    if support == ON_BFACES
         xregions=parent[BFaceRegions]
         xnodes=parent[BFaceNodes]
-        sub_gdim=dim_grid(parent)-1
         xct=parent[BFaceGeometries]
         sub_gdim=dim_grid(parent)-1
-    else
+    elseif support == ON_FACES
+        xregions=parent[FaceRegions]
+        xnodes=parent[FaceNodes]
+        xct=parent[FaceGeometries]
+        sub_gdim=dim_grid(parent)-1
+    elseif support == ON_CELLS
         xregions=parent[CellRegions]
         xnodes=parent[CellNodes]
         xct=parent[CellGeometries]
@@ -193,9 +208,9 @@ function subgrid(parent,
     subgrid[ParentGrid]=parent
     subgrid[NodeParents]=sub_nip
     subgrid[CellParents]=cellparents
-    subgrid[ParentGridRelation]=boundary ? BoundarySubGrid : SubGrid
+    subgrid[ParentGridRelation]=SubGrid{support}
 
-    if boundary
+    if support in [ON_BFACES, ON_FACES]
         subgrid[NumBFaceRegions]=0
         subgrid[BFaceRegions]=Ti[]
         subgrid[BFaceGeometries]=ElementGeometries[]
@@ -260,7 +275,7 @@ function subgrid(parent,
         subgrid[CoordinateSystem]=parent[CoordinateSystem]
     end
 
-    if sub_gdim == 1
+    if sub_gdim == 1 && project
         # Sort nodes of grid for easy plotting
         X=view(subgrid[Coordinates],1,:)
         nx=length(X)

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -243,7 +243,7 @@ end
     grid = grid_unitsquare(Triangle2D)
     grid[CellRegions] = Int32[1,2,2,1]
     sgrid = subgrid(grid, [1])
-    @test sgrid[ParentGridRelation] == SubGrid
+    @test sgrid[ParentGridRelation] == SubGrid{ON_CELLS}
 
     ## check if CellParents are assigned correctly
     @test sgrid[CellParents] == [1,4]