Changed the grid function used in the flooding alogrithm (to find the number

of separate surfaces) from integer to real type in preparation to get the
area calculation to work with Carpet. Introduced a module with fuzzy
comparison routines. Testsuite still passes. Still needs to be tested in
more complicated circumstances. Please report any problems.


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

6 files changed, 118 insertions, 25 deletions

diff --git a/interface.ccl b/interface.ccl
index ae4d0e7..e63187f 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -75,7 +75,7 @@ CCTK_INT niter_reinit TYPE=SCALAR
 #CCTK_INT carpet_re_init_control TYPE=SCALAR
 
 # Grid function used in counting surfaces.
-CCTK_INT surface_index TYPE=GF TIMELEVELS=1 tags='tensortypealias="Scalar" Prolongation="None"'
+CCTK_REAL surface_index TYPE=GF TIMELEVELS=1 tags='tensortypealias="Scalar"
 {
   sc
 }
diff --git a/src/EHFinder_Constants.F90 b/src/EHFinder_Constants.F90
index 2162631..c2c316e 100644
--- a/src/EHFinder_Constants.F90
+++ b/src/EHFinder_Constants.F90
@@ -16,4 +16,6 @@ module EHFinder_Constants
   CCTK_REAL, parameter :: fourthirds = four / three
   CCTK_REAL, parameter :: quarter = 0.25d0
   CCTK_REAL, parameter :: pi = 3.1415926535897932385d0
+  CCTK_REAL, parameter :: small = 1.0d-10
+  integer, parameter :: wp = kind(zero)
 end module EHFinder_Constants
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 63981bf..4662360 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -9,6 +9,7 @@
 subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
   use EHFinder_mod
+  use EHFinder_fuzzy
 
   implicit none
 
@@ -115,13 +116,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                    ( f(i,j,k,l) * f(i,j+1,k,l) .lt. zero ) .or. &
                    ( f(i,j,k,l) * f(i,j,k-1,l) .lt. zero ) .or. &
                    ( f(i,j,k,l) * f(i,j,k+1,l) .lt. zero ) ) .and. &
-                 ( ( sc(i,j,k) .eq. sn ) .or. &
-                   ( sc(i-1,j,k) .eq. sn ) .or. & 
-                   ( sc(i+1,j,k) .eq. sn ) .or. &
-                   ( sc(i,j-1,k) .eq. sn ) .or. &
-                   ( sc(i,j+1,k) .eq. sn ) .or. &
-                   ( sc(i,j,k-1) .eq. sn ) .or. &
-                   ( sc(i,j,k+1) .eq. sn ) ) ) then
+                 ( fuzzy_eq ( sc(i,j,k), real(sn,wp), small ) .or. &
+                   fuzzy_eq ( sc(i-1,j,k), real(sn,wp), small ) .or. & 
+                   fuzzy_eq ( sc(i+1,j,k), real(sn,wp), small ) .or. &
+                   fuzzy_eq ( sc(i,j-1,k), real(sn,wp), small ) .or. &
+                   fuzzy_eq ( sc(i,j+1,k), real(sn,wp), small ) .or. &
+                   fuzzy_eq ( sc(i,j,k-1), real(sn,wp), small ) .or. &
+                   fuzzy_eq ( sc(i,j,k+1), real(sn,wp), small ) ) ) then
               nc_loc = nc_loc + 1
               center_loc(1) = center_loc(1) + x(i,j,k)
               center_loc(2) = center_loc(2) + y(i,j,k)
@@ -144,12 +145,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
 ! Try to figure out the symmetries. First find the range in coordinates
 ! for the points inside the current surface.
-  crange_min_loc(1) = minval ( x, mask = sc .eq. sn )
-  crange_min_loc(2) = minval ( y, mask = sc .eq. sn )
-  crange_min_loc(3) = minval ( z, mask = sc .eq. sn )
-  crange_max_loc(1) = maxval ( x, mask = sc .eq. sn )
-  crange_max_loc(2) = maxval ( y, mask = sc .eq. sn )
-  crange_max_loc(3) = maxval ( z, mask = sc .eq. sn )
+  crange_min_loc(1) = minval ( x, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+  crange_min_loc(2) = minval ( y, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+  crange_min_loc(3) = minval ( z, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+  crange_max_loc(1) = maxval ( x, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+  crange_max_loc(2) = maxval ( y, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+  crange_max_loc(3) = maxval ( z, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
 
 ! Reduce over all processors.
   call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
@@ -741,7 +742,7 @@ subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS)
 
   more_surfaces = 1
   surface_counter = 0
-  sc = 0
+  sc = zero
 
 end subroutine EHFinder_CountSurfacesInit
   
@@ -749,6 +750,7 @@ end subroutine EHFinder_CountSurfacesInit
 subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
 
   use EHFinder_mod
+  use EHFinder_fuzzy
 
   implicit none
 
@@ -767,8 +769,8 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
 ! Find the location of the minimum among active points not already found
 ! to be inside a surface.
   minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter), &
-               ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. &
-               ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) )
+                  fuzzy_eq ( sc(ixl:ixr,jyl:jyr,kzl:kzr), zero, small ) .and. &
+                  ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) )
 
 ! Find the value of f at that location.
   minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1,levelset_counter)
@@ -804,7 +806,7 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
     points_counter = 0
     more_points = 1
     if ( procpos_loc .eq. procpos_glob ) then
-      sc(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1) = surface_counter
+      sc(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1) = real(surface_counter,wp)
     end if
   else
     more_surfaces = 0
@@ -818,6 +820,7 @@ end subroutine EHFinder_CountSurfaces
 subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
 
   use EHFinder_mod
+  use EHFinder_fuzzy
 
   implicit none
 
@@ -846,7 +849,7 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
         marked = .false.
         if ( ( eh_mask(i,j,k,levelset_counter) .ge. 0 ) .and. &
              ( f(i,j,k,levelset_counter) .lt. 0 ) .and. &
-             ( sc(i,j,k) .ne. surface_counter ) ) then
+             fuzzy_ne ( sc(i,j,k), real(surface_counter,wp), small ) ) then
 
 !         Find out if we are on the boundary
           if ( i + cctk_lbnd(1) .eq. 1 ) then
@@ -881,18 +884,21 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
           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
+          if ( any ( fuzzy_eq ( sc(i+i1:i+i2,j+j1:j+j2,k), &
+                     real(surface_counter,wp), small ) ) ) then
             marked = .true.
           end if
-          if ( any ( sc(i+i1:i+i2,j,k+k1:k+k2) .eq. surface_counter ) ) then
+          if ( any ( fuzzy_eq ( sc(i+i1:i+i2,j,k+k1:k+k2), &
+                     real(surface_counter,wp), small ) ) ) then
             marked = .true.
           end if
-          if ( any ( sc(i,j+j1:j+j2,k+k1:k+k2) .eq. surface_counter ) ) then
+          if ( any ( fuzzy_eq ( sc(i,j+j1:j+j2,k+k1:k+k2), &
+                     real(surface_counter,wp), small ) ) ) then
             marked = .true.
           end if
 
           if ( marked ) then
-            sc(i,j,k) = surface_counter
+            sc(i,j,k) = real(surface_counter,wp)
             n_loc = n_loc + 1
           end if
 
diff --git a/src/EHFinder_fuzzy.F90 b/src/EHFinder_fuzzy.F90
new file mode 100644
index 0000000..a856fef
--- /dev/null
+++ b/src/EHFinder_fuzzy.F90
@@ -0,0 +1,83 @@
+! Module with routines for fuzzy comparisons.
+! $Header$
+#include "cctk.h"
+
+module EHFinder_fuzzy
+
+  use EHFinder_Constants
+
+  interface fuzzy_eq
+    module procedure fuzzy_eq_0d, fuzzy_eq_1d, fuzzy_eq_2d, fuzzy_eq_3d
+  end interface
+
+  interface fuzzy_ne
+    module procedure fuzzy_ne_0d, fuzzy_ne_1d, fuzzy_ne_2d, fuzzy_ne_3d
+  end interface
+
+contains
+
+  function fuzzy_eq_0d ( x, y, eps )
+    CCTK_REAL, intent(IN) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical :: fuzzy_eq_0d
+
+    fuzzy_eq_0d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+  end function fuzzy_eq_0d 
+
+  function fuzzy_eq_1d ( x, y, eps )
+    CCTK_REAL, intent(IN), dimension(:) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical, dimension(size(x,1)) :: fuzzy_eq_1d
+
+    fuzzy_eq_1d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+  end function fuzzy_eq_1d 
+
+  function fuzzy_eq_2d ( x, y, eps )
+    CCTK_REAL, intent(IN), dimension(:,:) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical, dimension(size(x,1),size(x,2)) :: fuzzy_eq_2d
+
+    fuzzy_eq_2d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+  end function fuzzy_eq_2d 
+
+  function fuzzy_eq_3d ( x, y, eps )
+    CCTK_REAL, intent(IN), dimension(:,:,:) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical, dimension(size(x,1),size(x,2),size(x,3)) :: fuzzy_eq_3d
+
+    fuzzy_eq_3d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+  end function fuzzy_eq_3d 
+
+  function fuzzy_ne_0d ( x, y, eps )
+    CCTK_REAL, intent(IN) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical :: fuzzy_ne_0d
+
+    fuzzy_ne_0d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+  end function fuzzy_ne_0d 
+
+  function fuzzy_ne_1d ( x, y, eps )
+    CCTK_REAL, intent(IN), dimension(:) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical, dimension(size(x,1)) :: fuzzy_ne_1d
+
+    fuzzy_ne_1d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+  end function fuzzy_ne_1d 
+
+  function fuzzy_ne_2d ( x, y, eps )
+    CCTK_REAL, intent(IN), dimension(:,:) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical, dimension(size(x,1),size(x,2)) :: fuzzy_ne_2d
+
+    fuzzy_ne_2d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+  end function fuzzy_ne_2d 
+
+  function fuzzy_ne_3d ( x, y, eps )
+    CCTK_REAL, intent(IN), dimension(:,:,:) :: x
+    CCTK_REAL, intent(IN) :: y, eps
+    logical, dimension(size(x,1),size(x,2),size(x,3)) :: fuzzy_ne_3d
+
+    fuzzy_ne_3d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+  end function fuzzy_ne_3d 
+
+end module EHFinder_fuzzy
diff --git a/src/make.code.defn b/src/make.code.defn
index 701a619..d24f363 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -18,7 +18,8 @@ SRCS = EHFinder_mod.F90 \
        EHFinder_Generator_Sources.F90 \
        EHFinder_Generator_Sources2.F90 \
        EHFinder_ReInitialize.F90 \
-       EHFinder_CarpetSupport.cc
+       EHFinder_CarpetSupport.cc \
+       EHFinder_fuzzy.F90
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/make.code.deps b/src/make.code.deps
index 0b34034..21d5cc2 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -3,6 +3,7 @@
 
 EHFinder_mod.F90.o:			EHFinder_Constants.F90.o
 EHFinder_IsoSurface_mod.F90.o:		EHFinder_Constants.F90.o
+EHFinder_fuzzy.F90.o:			EHFinder_Constants.F90.o
 EHFinder_Init.F90.o:			EHFinder_mod.F90.o
 EHFinder_Sources.F90.o:			EHFinder_mod.F90.o
 EHFinder_Generator_Sources.F90.o:	EHFinder_mod.F90.o
@@ -14,5 +15,5 @@ EHFinder_SetSym.F90.o:			EHFinder_mod.F90.o
 EHFinder_ReadData.F90.o:		EHFinder_mod.F90.o
 EHFinder_Check.F90.o:			EHFinder_mod.F90.o
 EHFinder_Integrate.F90.o:		EHFinder_mod.F90.o
-EHFinder_FindSurface.F90.o:		EHFinder_mod.F90.o
+EHFinder_FindSurface.F90.o:		EHFinder_mod.F90.o EHFinder_fuzzy.F90.o
 EHFinder_IsoSurface.F90.o:		EHFinder_mod.F90.o IsoSurface_mod.F90.o