Changed EHFinder_MarkSurfaces to check not only the nearest neighbours

when marking points, but also the nearest neighbours of the nearest
neighbours. The reason being that connected surfaces can contain points
where none of the nearest neighbours are inside the surface, but at least
one of the next nearest neighbours are.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@118 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 92 insertions, 24 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 221d969..95e0e63 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -792,6 +792,7 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
 
   CCTK_INT :: n_loc, n_glob, i, j, k
+  CCTK_INT :: i1, i2, j1, j2, k1, k2
   logical :: marked
 
 ! Find index ranges for points excluding ghost zones and symmetry points for
@@ -812,40 +813,107 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
         if ( ( eh_mask(i,j,k) .ge. 0 ) .and. &
              ( f(i,j,k) .lt. 0 ) .and. &
              ( sc(i,j,k) .ne. surface_counter ) ) then
-          if ( .not. btest(eh_mask(i,j,k),0) ) then
-            if ( sc(i-1,j,k) .eq. surface_counter ) then
-              marked = .true.
-            end if
+
+          ! Find out if we are on the boundary
+          if ( i + cctk_lbnd(1) .eq. 1 ) then
+            i1 = 0
+          else
+            i1 = -1
           end if
-          if ( .not. btest(eh_mask(i,j,k),1) ) then
-            if ( sc(i+1,j,k) .eq. surface_counter ) then
-              marked = .true.
-            end if
+          if ( i + cctk_lbnd(1) .eq. cctk_gsh(1) ) then
+            i2 = 0
+          else
+            i2 = 1
           end if
-          if ( .not. btest(eh_mask(i,j,k),2) ) then
-            if ( sc(i,j-1,k) .eq. surface_counter ) then
-              marked = .true.
-            end if
+          if ( j + cctk_lbnd(2) .eq. 1 ) then
+            j1 = 0
+          else
+            j1 = -1
           end if
-          if ( .not. btest(eh_mask(i,j,k),3) ) then
-            if ( sc(i,j+1,k) .eq. surface_counter ) then
-              marked = .true.
-            end if
+          if ( j + cctk_lbnd(2) .eq. cctk_gsh(2) ) then
+            j2 = 0
+          else
+            j2 = 1
           end if
-          if ( .not. btest(eh_mask(i,j,k),4) ) then
-            if ( sc(i,j,k-1) .eq. surface_counter ) then
-              marked = .true.
-            end if
+          if ( k + cctk_lbnd(3) .eq. 1 ) then
+            k1 = 0
+          else
+            k1 = -1
           end if
-          if ( .not. btest(eh_mask(i,j,k),5) ) then
-            if ( sc(i,j,k+1) .eq. surface_counter ) then
-              marked = .true.
-            end if
+          if ( k + cctk_lbnd(3) .eq. cctk_gsh(3) ) then
+            k2 = 0
+          else
+            k2 = 1
           end if
+
+!         Then check if any valid neighbours are marked
+          if ( any ( sc(i+i1:i+i2,j+j1:j+j2,k) .eq. surface_counter ) ) then
+            marked = .true.
+          end if
+          if ( any ( sc(i+i1:i+i2,j,k+k1:k+k2) .eq. surface_counter ) ) then
+            marked = .true.
+          end if
+          if ( any ( sc(i,j+j1:j+j2,k+k1:k+k2) .eq. surface_counter ) ) then
+            marked = .true.
+          end if
+
+!          if ( any ( sc(i+i1:i+i2,j+j1:j+j2,k+k1:k+k2) &
+!                                       .eq. surface_counter ) ) then
+!            marked = .true.
+!          end if
+
+!          ! First look in the -x direction.
+!          if ( .not. btest(eh_mask(i,j,k),0) ) then
+!            if ( sc(i-1,j,k) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!            if ( .not. btest(eh_mask(i-1,j,k),2) ) then
+!              if ( sc(i-1,j-1,k) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!            if ( .not. btest(eh_mask(i-1,j,k),3) ) then
+!              if ( sc(i-1,j+1,k) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!            if ( .not. btest(eh_mask(i-1,j,k),4) ) then
+!              if ( sc(i-1,j,k-1) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!            if ( .not. btest(eh_mask(i-1,j,k),5) ) then
+!              if ( sc(i-1,j,k+1) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!          end if
+!          if ( .not. btest(eh_mask(i,j,k),1) ) then
+!            if ( sc(i+1,j,k) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!          end if
+!          if ( .not. btest(eh_mask(i,j,k),2) ) then
+!            if ( sc(i,j-1,k) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!          end if
+!          if ( .not. btest(eh_mask(i,j,k),3) ) then
+!            if ( sc(i,j+1,k) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!          end if
+!          if ( .not. btest(eh_mask(i,j,k),4) ) then
+!            if ( sc(i,j,k-1) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!          end if
+!          if ( .not. btest(eh_mask(i,j,k),5) ) then
+!            if ( sc(i,j,k+1) .eq. surface_counter ) then
+!              marked = .true.
+!            end if
+!          end if
           if ( marked ) then
             sc(i,j,k) = surface_counter
             n_loc = n_loc + 1
           end if
+
         end if    
       end do
     end do