From a1ff15bd67bbd22c6fcc986eb623866e545bc629 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 10 Jan 2006 15:37:08 +0000 Subject: 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 --- interface.ccl | 2 +- src/EHFinder_Constants.F90 | 2 ++ src/EHFinder_FindSurface.F90 | 50 ++++++++++++++------------ src/EHFinder_fuzzy.F90 | 83 ++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 3 +- src/make.code.deps | 3 +- 6 files changed, 118 insertions(+), 25 deletions(-) create mode 100644 src/EHFinder_fuzzy.F90 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 -- cgit v1.2.3