Added routines to find out how many surfaces are present. At present

this information is not used. Added beginning support for the spacemask.
Only routine to read in the emask. Nothing is done with it.


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

1 files changed, 203 insertions, 0 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 97515cb..fe41fd1 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -432,3 +432,206 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   call Util_TableDestroy ( status, table_handle )
   
 end subroutine EHFinder_FindSurface
+
+
+subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+!  print*,'EHFinder_CountSurfacesInit called'
+  more_surfaces = 1
+  surface_counter = 0
+  sc = 0
+
+end subroutine EHFinder_CountSurfacesInit
+  
+
+subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT, dimension(3) :: minl
+  CCTK_REAL :: minf_loc, minf_glob
+  CCTK_INT :: procpos_loc, procpos_glob
+
+! Find index ranges for points excluding ghost zones and symmetry points for
+! the 3D grid functions.
+#include "include/physical_part.h"
+
+  ! 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), &
+               ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. &
+               ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .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)
+
+  ! Find the minimum over all processors.
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, &
+                              minf_loc, minf_glob, CCTK_VARIABLE_REAL )
+  if ( ierr .ne. 0 ) then
+    call CCTK_WARN(0,'Reduction of minf_glob failed')
+  end if
+
+  ! To figure out on which processor the minimum is located I set procpos_loc
+  ! to the local processor number if the local minimum is equal to the global
+  ! minimum otherwise I set it to the total number of processors. When I then
+  ! find the minimum of procpos_loc over all processors I get the processor
+  ! containing the minimum.
+  if ( minf_loc .eq. minf_glob ) then
+    procpos_loc = CCTK_MyProc ( cctkGH )
+  else
+    procpos_loc = CCTK_nProcs ( cctkGH )
+  end if
+  
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, &
+                              procpos_loc, procpos_glob, CCTK_VARIABLE_INT )
+  if ( ierr .ne. 0 ) then
+    call CCTK_WARN(0,'Reduction of minf_glob failed')
+  end if
+
+  ! If minf_glob<0 there must exist another surface so lets set up the variables
+  ! to start marking and counting these points.
+  if ( minf_glob .lt. zero ) then
+    surface_counter = surface_counter + 1
+    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
+    end if
+  else
+    more_surfaces = 0
+    more_points = 0
+  end if
+
+end subroutine EHFinder_CountSurfaces    
+
+
+! This routine marks all points inside the current surface.
+subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: n_loc, n_glob, i, j, k
+  logical :: marked
+
+! Find index ranges for points excluding ghost zones and symmetry points for
+! the 3D grid functions.
+#include "include/physical_part.h"
+  
+  ! Initialise the local point counter.
+  n_loc = 0
+
+  ! Check if any points have a marked point as a neighbour. Only look at
+  ! active points with f<0 that are not already marked. Take care at the
+  ! boundaries (physical and excision) not to look at inactive or not
+  ! exisiting points.
+  do k = kzl, kzr
+    do j = jyl, jyr
+      do i = ixl, ixr
+        marked = .false.
+        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
+          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
+  end do
+
+  ! Reduce the number of newly marked points over all processors
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
+                              n_loc, n_glob, CCTK_VARIABLE_INT )
+  if ( ierr .ne. 0 ) then
+    call CCTK_WARN(0,'Reduction of n_glob failed')
+  end if
+
+  ! If the total number of newly marked points is 0, set more_points to 0
+  ! to indicate that we are finished. Otherwise continue and add the number
+  ! of points to the points_counter.
+  if ( n_glob .eq. 0 ) then
+    more_points = 0
+  end if
+  points_counter = points_counter + n_glob
+
+  ! This should not be necessary, since I have a SYNC in the schedule.ccl,
+  ! but for some reason that does not work.
+  call CCTK_SyncGroup ( ierr, cctkGH, "ehfinder::surface_index" )
+
+end subroutine EHFinder_MarkSurfaces
+
+
+!This routine prints out information about the found surfaces.
+subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  character(len=80) :: info_message
+
+  ! Write out how many surfaces were found.
+  if ( surface_counter .eq. 1 ) then
+    write(info_message,'(a10,i3,a8)') 'There are ',surface_counter,' surface'
+  else
+    write(info_message,'(a10,i3,a9)') 'There are ',surface_counter,' surfaces'
+  end if
+  call CCTK_INFO(info_message)
+
+end subroutine EHFinder_InfoSurfaces