fix for subgrid in case there is empty boundary data (#36)

* fix for subgrid in case there is empty boundary data in the parent grid (or probably also when the subgrid lies completely in the interior)

* modified the subgrid+extrusion test a bit to test if a subgrid for a completely interior region works (and it only does with the previous modifications)

src/subgrid.jl

@@ -194,8 +194,13 @@ function subgrid(parent,
 
         subgrid[BFaceRegions]=sub_bfaceregions
         subgrid[BFaceGeometries]=sub_bfacetypes
-        subgrid[BFaceNodes]=tryfix(sub_bfacenodes)
-        subgrid[NumBFaceRegions]=maximum(sub_bfaceregions)
+        if length(sub_bfaceregions) > 1
+            subgrid[BFaceNodes]=tryfix(sub_bfacenodes)
+            subgrid[NumBFaceRegions]=maximum(sub_bfaceregions)
+        else
+            subgrid[BFaceNodes]=zeros(Ti, 2, 0)
+            subgrid[NumBFaceRegions]=0
+        end
     end
     subgrid[CoordinateSystem]=parent[CoordinateSystem]
 

test/runtests.jl

@@ -146,7 +146,7 @@ end
 end
 
 @testset "subgrid+extrusion" begin
-    function subgen(; h = 0.2)
+    function subgen(; region = 1, h = 0.2)
         X = collect(-1:h:2)
         Y = collect(0:h:1)
         Z = collect(0:h:2)
@@ -158,15 +158,19 @@ end
         # Mark cut-out regions as 2
         cellmask!(g, [-1, 0, 0], [0, 1, 1], 2)
         cellmask!(g, [1, 0, 0], [2, 1, 1], 2)
+        # Mark interior region 3 with region 2
+        cellmask!(g, [-0.75,0.25,0.25], [-0.25,0.75,0.75], 3)
 
         # add new interface elements
         bfacemask!(g, [-1, 0, 1], [2, 1, 1], 3; allow_new = true)
         bfacemask!(g, [0, 0, 0], [0, 1, 1], 3; allow_new = true)
         bfacemask!(g, [1, 0, 0], [1, 1, 1], 3; allow_new = true)
 
         # return subgrid of region 1
-        subgrid(g, [1])
+        subgrid(g, [region])
     end
+    sub2 = subgen(; region = 3)
+
     @test numbers_match(subgen(), 756, 3000, 950)
 
     X = collect(0:0.25:1)