diff options
Diffstat (limited to 'src/EHFinder_FindSurface.F90')
| -rw-r--r-- | src/EHFinder_FindSurface.F90 | 50 |
1 files changed, 28 insertions, 22 deletions
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 |