From 5ce2a2f3f7402ea8dc4f702a09645d9f4a477af5 Mon Sep 17 00:00:00 2001
From: Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>
Date: Thu, 21 Oct 2021 17:06:26 +0200
Subject: [PATCH] bfacemask! for 3D, introduce allow_new flag

---
 Project.toml          |  2 +-
 src/extendablegrid.jl |  8 ++++
 src/regionedit.jl     | 88 ++++++++++++++++++++-----------------------
 3 files changed, 49 insertions(+), 49 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2d537969..885e27ee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ExtendableGrids"
 uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>, Christian Merdon  <christian.merdon@wias-berlin.de>"]
-version = "0.8.5"
+version = "0.8.6"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
diff --git a/src/extendablegrid.jl b/src/extendablegrid.jl
index 8ed3939f..1b21bf21 100644
--- a/src/extendablegrid.jl
+++ b/src/extendablegrid.jl
@@ -207,6 +207,14 @@ Set new grid component
 """
 Base.setindex!(grid::ExtendableGrid,v,T::Type{<:AbstractGridComponent})= grid.components[T]=veryform(grid,v,T)
 
+"""
+$(TYPEDSIGNATURES)
+
+Remove grid component
+"""
+Base.delete!(grid::ExtendableGrid,T::Type{<:AbstractGridComponent})= delete!(grid.components,T)
+
+
 
 """
 $(TYPEDSIGNATURES)
diff --git a/src/regionedit.jl b/src/regionedit.jl
index ac6b3984..c75bb350 100644
--- a/src/regionedit.jl
+++ b/src/regionedit.jl
@@ -40,85 +40,77 @@ end
 $(TYPEDSIGNATURES)
 
 Edit region numbers of grid  boundary facets  via rectangular mask.
-For 1D grids, inner boundaries can be added by this method.
+If `allow_new` is true (default), new facets are added
 """
 function bfacemask!(grid::ExtendableGrid,
                     maskmin::AbstractArray,
                     maskmax::AbstractArray,
                     ireg::Int;
+                    allow_new=true,
                     tol=1.0e-10)
 
     xmaskmin=maskmin.-tol
     xmaskmax=maskmax.+tol
 
 
-    nbfaces=num_bfaces(grid)
     bfacenodes=grid[BFaceNodes]
+    nbfaces=size(bfacenodes,2)
+
+
+
+    
     Ti=eltype(bfacenodes)
     dim=dim_space(grid)
     bfaceregions=grid[BFaceRegions]
-    new_bfacenodes=ElasticArray{Ti,2}(bfacenodes)
     coord=grid[Coordinates]
-    
-    function isbface(ix)
-        for ibface=1:num_bfaces(grid)
-            if bfacenodes[1,ibface]==ix
-                return ibface
-            end
-        end
-        return 0
-    end
 
-    function isbface(ix,iy)
-        for ibface=1:num_bfaces(grid)
-            if (bfacenodes[1,ibface] == ix && bfacenodes[2,ibface] == iy) ||
-                (bfacenodes[1,ibface] == iy && bfacenodes[2,ibface] == ix)
-                return ibface
-            end
+    xfacenodes=bfacenodes
+    if allow_new
+        new_bfacenodes=ElasticArray{Ti,2}(bfacenodes)
+        facenodes=grid[FaceNodes]
+        bfacefaces=grid[BFaceFaces]
+        nfaces=size(facenodes,2)
+        bmark=zeros(Int,nfaces)
+        for ibface=1:nbfaces
+            bmark[bfacefaces[ibface]]=ibface
         end
-        return 0
+        xfacenodes=facenodes
     end
-
-    if dim_space(grid)==1
-        for inode=1:num_nodes(grid)
-            x=coord[1,inode]
-            if x>xmaskmin[1] && x<xmaskmax[1]
-                ibface=isbface(inode)
-                if ibface>0
-                    bfaceregions[ibface]=ireg
-                else
-                    push!(bfaceregions,ireg)
-                    append!(new_bfacenodes,[inode])
+    
+    for ixface=1:size(xfacenodes,2)
+        in_region=true
+        for inode=1:num_targets(xfacenodes,ixface)
+            ignode=xfacenodes[inode,ixface]
+            for idim=1:dim_space(grid)
+                if coord[idim,ignode]<xmaskmin[idim]
+                    in_region=false
+                elseif coord[idim,ignode]>xmaskmax[idim]
+                    in_region=false
                 end
             end
         end
-    else
-        facenodes=grid[FaceNodes]
-        for iface=1:size(facenodes,2)
-            in_region=true
-            for inode=1:num_targets(facenodes,iface)
-                ignode=facenodes[inode,iface]
-                for idim=1:dim_space(grid)
-                    if coord[idim,ignode]<xmaskmin[idim]
-                        in_region=false
-                    elseif coord[idim,ignode]>xmaskmax[idim]
-                        in_region=false
-                    end
-                end
-            end
-            if in_region
-                ibface=isbface(facenodes[1,iface],facenodes[2,iface])
+
+        if in_region
+            if allow_new
+                ibface=bmark[ixface]
                 if ibface>0
                     bfaceregions[ibface]=ireg
                 else
                     push!(bfaceregions,ireg)
-                    append!(new_bfacenodes,[facenodes[1,iface],facenodes[2,iface]])
+                    @views append!(new_bfacenodes,xfacenodes[:,ixface])
                 end
+            else
+                bfaceregions[ixface]=ireg
             end
         end
     end
+
+
+    delete!(grid,BFaceFaces)
     btype=grid[BFaceGeometries][1]
-    grid[BFaceNodes]=Array{Ti,2}(new_bfacenodes)
+    if allow_new
+        grid[BFaceNodes]=Array{Ti,2}(new_bfacenodes)
+    end
     grid[BFaceGeometries]=VectorOfConstants(btype,length(bfaceregions))
     grid[NumBFaceRegions]=max(num_bfaceregions(grid),ireg)
     return grid