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

6 files changed, 326 insertions, 1 deletions

diff --git a/interface.ccl b/interface.ccl
index 6631138..f4fe515 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,7 +2,7 @@
 # $Header$
 
 implements: ehfinder
-inherits: grid admbase coordgauge staticconformal boundary
+inherits: grid admbase coordgauge staticconformal spacemask boundary
 
 USES INCLUDE: Boundary.h
 USES INCLUDE: MoL.h
@@ -44,6 +44,16 @@ CCTK_INT mask_bak TYPE=GF TIMELEVELS=1
   eh_mask_bak
 } "Temporary placeholder for the mask during re-parametrization"
 
+CCTK_INT surface_index TYPE=GF TIMELEVELS=1
+{
+  sc
+}
+
+CCTK_INT surface_counter TYPE=SCALAR
+CCTK_INT points_counter TYPE=SCALAR
+CCTK_INT more_surfaces TYPE=SCALAR
+CCTK_INT more_points TYPE=SCALAR
+ 
 CCTK_REAL surface_arrays TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=ntheta,nphi GHOSTSIZE=n_array_ghosts,n_array_ghosts DISTRIB=DEFAULT
 {
   ctheta, cphi, rsurf, sintheta, costheta, sinphi, cosphi,
diff --git a/param.ccl b/param.ccl
index 66cc38d..ca7335f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -172,11 +172,16 @@ INT nphi "Number of points in the phi direction when finding points on the surfa
 {
 1:*:2            :: "Positive and odd please"
 } 51
+
 INT n_array_ghosts "Number of ghost points in the surface grid array"
 {
 1:               :: "Positive please"
 } 1
 
+INT maximum_surface_number "The maximum number of surfaces expected in the data"{
+  1:*            :: "Positive please"
+} 1
+
 BOOLEAN use_user_center "Should the user prescribed center be used"
 {
 } "no"
@@ -221,3 +226,6 @@ EXTENDS KEYWORD initial_shift
   "read from file" :: "Read in shift from a file"
 }
 
+shares: spacemask
+
+USES KEYWORD use_mask
diff --git a/schedule.ccl b/schedule.ccl
index 47bafca..572ea93 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -18,6 +18,7 @@ else
   STORAGE: surface_arrays
   STORAGE: center_arrays
   STORAGE: eh_area2
+#  STORAGE: surface_index
 }
 
 # Check for metric_state
@@ -63,6 +64,14 @@ if (CCTK_Equals(metric_type,"static conformal") && CCTK_Equals(eh_metric_type,"n
   } "Read in conformal factor from file"
 }
 
+if (use_mask)
+{
+  schedule EHFinder_Read_Mask at CCTK_INITIAL
+  {
+    LANG: Fortran
+  } "Read in excision mask from file"
+}
+
 # Set up the initial level set function
 
 schedule EHFinder_Init at CCTK_POSTINITIAL
@@ -76,6 +85,45 @@ schedule EHFinder_InitWeights at CCTK_POSTINITIAL
   LANG: FORTRAN
 } "Setup weights for Simpson integration"
 
+schedule GROUP EHFinder_Surfaces at CCTK_ANALYSIS
+{
+  STORAGE: more_surfaces, surface_counter, more_points, points_counter
+  STORAGE: surface_index
+  TRIGGER: eh_area2, eh_centroid
+} "Schedule group for counting the number of surfaces and integrating over them"
+
+
+schedule EHFinder_CountSurfacesInit in EHFinder_Surfaces
+{
+  LANG: Fortran
+  TRIGGER: eh_area2, eh_centroid
+} "Initialising while loop control for the group EHFinder_CountMarkSurfaces"
+
+schedule GROUP EHFinder_CountMarkSurfaces in EHFinder_Surfaces after EHFinder_CountSurfacesInit WHILE ehfinder::more_surfaces
+{
+  TRIGGER: eh_area2, eh_centroid
+} "Scedule group for counting and marking surfaces"
+
+schedule EHFinder_CountSurfaces in EHFinder_CountMarkSurfaces
+{
+  LANG: Fortran
+  TRIGGER: eh_area2, eh_centroid
+  SYNC: surface_index
+} "Check if there are more surfaces and initialise while loop control for EHFinder_MarkSurfaces"
+
+schedule EHFinder_MarkSurfaces in EHFinder_CountMarkSurfaces after EHFinder_CountSurfaces WHILE ehfinder::more_points
+{
+  LANG: Fortran
+  TRIGGER: eh_area2, eh_centroid
+  SYNC: surface_index
+} "Mark points inside the current surface until they are all marked"
+
+schedule EHFinder_InfoSurfaces in EHFinder_Surfaces after EHFinder_CountMarkSurfaces
+{
+  LANG: Fortran
+  TRIGGER: eh_area2, eh_centroid
+} "Output input about found surfaces"
+
 schedule EHFinder_FindSurface at CCTK_ANALYSIS
 {
   LANG: Fortran
@@ -145,6 +193,14 @@ if (CCTK_Equals(metric_type,"static conformal") && CCTK_Equals(eh_metric_type,"n
   } "Read in conformal factor from file"
 }
 
+if (use_mask)
+{
+  schedule EHFinder_Read_Mask at CCTK_PRESTEP
+  {
+    LANG: Fortran
+  } "Read in excision mask from file"
+}
+
 #schedule EHFinder_ReadData at CCTK_PRESTEP
 #{
 #  LANG: Fortran
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
diff --git a/src/EHFinder_ReadData.F90 b/src/EHFinder_ReadData.F90
index 68516b6..75d633f 100644
--- a/src/EHFinder_ReadData.F90
+++ b/src/EHFinder_ReadData.F90
@@ -157,6 +157,34 @@ subroutine EHFinder_Read_Conformal(CCTK_ARGUMENTS)
 end subroutine EHFinder_Read_Conformal
 
 
+! This routine reads in the excision mask (from SpaceMask).
+subroutine EHFinder_Read_Mask(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  character(len=128) :: in_files, in_vars
+  character(len=10) :: iteration_string
+  CCTK_INT :: nc, res
+
+  write(iteration_string,'(i10)') last_iteration_number - &
+                                  saved_iteration_every *cctk_iteration
+  iteration_string = adjustl(iteration_string)
+  nc = len_trim(iteration_string)
+
+  in_files = 'emask'
+  in_vars = 'spacemask::emask[cctk_iteration='//iteration_string(1:nc)//']'
+
+  call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
+
+end subroutine EHFinder_Read_Mask
+
+
 ! This routine reads in everything.
 subroutine EHFinder_ReadData(CCTK_ARGUMENTS)
 
diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90
index 04ed760..ea8fa57 100644
--- a/src/EHFinder_SetMask.F90
+++ b/src/EHFinder_SetMask.F90
@@ -1,10 +1,24 @@
 ! Mask modification functions.
+! The grid function eh_mask is used to encode excision and boundary
+! information. An excised point has the value -1, while an active point
+! is 0 or positive. If it is zero the point has neighbours in all 
+! directions while if it is positive the value encodes which neighbours are
+! missing. The mask uses 6 bits of an integer to encode this information.
+! If the neighbour in the -x-direction is missing the mask is set 
+! to 1 (b000001). If the neighbour in the +x-direction is missing the mask
+! is set to 2 (b000010). For the y-directions the values are 4 (b000100) and 
+! 8 (b001000). For the z-direction they are 16 (b010000) and 32 (b100000).
+! For example a point with no active neighbours in the +x, -y and +z direction
+! will have a mask value of 38 (b100110).
 ! $Header$
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
+! This routine is called only once to initialise the mask at the physical
+! outer boundaries. The value of the mask in these points should never
+! be changed again.
 subroutine EHFinder_MaskInit(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -17,10 +31,16 @@ subroutine EHFinder_MaskInit(CCTK_ARGUMENTS)
 
   CCTK_INT :: i, j, k
 
+  ! Initialise the mask to 0.
   eh_mask = 0
 
+! Figure out the local dimensions of the grid excluding ghostzones and
+! symmetry zones.
 # include "include/physical_part.h"
 
+  ! Check if the point is located at a physical outer boundary and set the
+  ! eh_mask appropriately. The array ll is defined in EHFinder_mod.F90
+  ! and contains the values ( 1, 2, 4, 8, 16, 32 ).
   if ( ( cctk_bbox(1) .eq. 1 ) .and. ( ixl .eq. 1 ) ) then
     eh_mask(1,:,:) = eh_mask(1,:,:) + ll(0)
   end if