Add possibility of a cubical mask.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@58 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 44 insertions, 14 deletions

diff --git a/src/AHFinder_mask.F b/src/AHFinder_mask.F
index 21fac55..0050f2b 100644
--- a/src/AHFinder_mask.F
+++ b/src/AHFinder_mask.F
@@ -14,10 +14,8 @@
       subroutine AHFinder_mask(CCTK_FARGUMENTS,rhor)
 
 !     Set up the horizon mask.  All this routine does is
-!     set up the value of the grid function "ahmask" to
-!     one for points outside the horizon, and to zero
-!     for points inside the horizon.  This is only necessary
-!     if we want to use black hole excision.
+!     set up the value of the grid function "ahmask" to zero
+!     for points inside the horizon.
 !
 !     It also finds out the radius of the largest sphere that
 !     fits inside the horizon and it centered in (xc,yc,zc).
@@ -48,28 +46,31 @@
 !     ***   START ROUTINE   ***
 !     *************************
 
-      call AHFinder_fun(CCTK_FARGUMENTS)
-
-!     If outside of horizon set mask to 1.  For this
+!     If inside of horizon (minus a bufferzone) set the mask to 0.
+!     To find out how big is the buffer zone inside theFor this
 !     I check the parameter ahf_maskshrink.  This parameter
-!     says by which factor of the monopole term to shrink
-!     the mass with respect to the real horizon.  Here I
-!     use the fact that the function ahfgrid is zero at
-!     the horizon, and is linear on the monopole term
-!     (which measures the overall scale).
+!     says by which factor of the monopole term to shrink the
+!     mask with respect to the real horizon.  Here I use the fact
+!     that the function ahfgrid is zero at the horizon, and is linear
+!     on the monopole term (which measures the overall scale).
+!
+!     Notice that if I want a cubical mask, I first find out the radius
+!     of the largest sphere that fits inside the horizon and only after
+!     that I set the mask.
 
       rhor = 1.0D10
 
+      call AHFinder_fun(CCTK_FARGUMENTS)
+
       do k=1,nz
          do j=1,ny
             do i=1,nx
 
                if (ahfgrid(i,j,k).ge.(ahf_maskshrink-one)*c0(0)) then
-                  ahmask(i,j,k) = one
                   aux = dsqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2
      .                + (z(i,j,k)-xc)**2)
                   if (rhor.gt.aux) rhor = aux
-               else
+               else if (ahf_maskcube.eq.0) then
                   ahmask(i,j,k) = zero
                end if
 
@@ -78,6 +79,35 @@
       end do
 
 
+!     ************************
+!     ***   CUBICAL MASK   ***
+!     ************************
+
+!     If we want a cubical mask, set the mask based
+!     on the value of rhor.
+
+      if (ahf_maskcube.ne.0) then
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  if ((x(i,j,k).gt.(xc-0.5D0*rhor)).and.
+     .                (y(i,j,k).gt.(yc-0.5D0*rhor)).and.
+     .                (z(i,j,k).gt.(zc-0.5D0*rhor)).and.
+     .                (x(i,j,k).lt.(xc+0.5D0*rhor)).and.
+     .                (y(i,j,k).lt.(yc+0.5D0*rhor)).and.
+     .                (z(i,j,k).lt.(zc+0.5D0*rhor))) then
+                     ahmask(i,j,k) = zero
+                  end if
+
+               end do
+            end do
+         end do
+
+      end if
+
+
 !     ***************
 !     ***   END   ***
 !     ***************