Add parameters to set the origin of the excision region for "fixed

excision".

When checking for (real-valued) mask values, do not check for
equality, but instead allow for a fudge factor.

Check that the thorn is activated when any of its routines are called.

When there is only a single point excised, do not abort with an
internal error, but pretend the normal extends into the x direction.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/LegoExcision/trunk@42 f75ba9e5-694f-0410-ac2c-87ea7ce7132b

5 files changed, 61 insertions, 19 deletions

diff --git a/param.ccl b/param.ccl
index 6769626..20176da 100644
--- a/param.ccl
+++ b/param.ccl
@@ -12,3 +12,18 @@ REAL fixed_size "The size (radius or corner location) of the fixed region"
 {
   0:* :: "Positive"
 } 1.0
+
+REAL fixed_origin_x "The x coordinate of the centre of the fixed region"
+{
+  *:* :: ""
+} 0.0
+
+REAL fixed_origin_y "The y coordinate of the centre of the fixed region"
+{
+  *:* :: ""
+} 0.0
+
+REAL fixed_origin_z "The z coordinate of the centre of the fixed region"
+{
+  *:* :: ""
+} 0.0
diff --git a/src/extrapolate.F90 b/src/extrapolate.F90
index 19ca1b6..c09e35c 100644
--- a/src/extrapolate.F90
+++ b/src/extrapolate.F90
@@ -5,9 +5,8 @@
 ! diry, and dirz for extrapolation.
 
 #include "cctk.h"
-#include "cctk_Parameters.h"
-
 #include "cctk_Functions.h"
+#include "cctk_Parameters.h"
 
 #include "maskvalues.h"
 
@@ -46,15 +45,19 @@ subroutine excision_extrapolate (ierr, var, oldvar, &
   integer i,j,k
   integer ii,jj,kk
   
+  if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then
+     call CCTK_WARN (0, "The routine excision_extrapolate was called, but thorn " // CCTK_THORNSTRING // " is not active")
+  end if
+  
   do k=2,nk-1
      do j=2,nj-1
         do i=2,ni-1
-           if (mask(i,j,k)==MASK_BOUNDARY) then
+           if (abs(mask(i,j,k)-MASK_BOUNDARY)<0.01) then
               ii = i + int(dirx(i,j,k))
               jj = j + int(diry(i,j,k))
               kk = k + int(dirz(i,j,k))
               var(i,j,k) = oldvar(i,j,k) + var(ii,jj,kk) - oldvar(ii,jj,kk)
-           else if (mask(i,j,k)==MASK_EXCISED) then
+           else if (abs(mask(i,j,k)-MASK_EXCISED)<0.01) then
               var(i,j,k) = var0
            end if
         end do
diff --git a/src/findboundary.F90 b/src/findboundary.F90
index 454c91a..2eaa1e1 100644
--- a/src/findboundary.F90
+++ b/src/findboundary.F90
@@ -4,9 +4,8 @@
 ! Return a mask where the outermost 0s have been replaced by 0.5s.
 
 #include "cctk.h"
-#include "cctk_Parameters.h"
-
 #include "cctk_Functions.h"
+#include "cctk_Parameters.h"
 
 #include "maskvalues.h"
 
@@ -34,6 +33,10 @@ subroutine excision_findboundary (ierr, mask, ni, nj, nk)
   integer i,j,k
   integer ii,jj,kk
   logical bnd
+  
+  if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then
+     call CCTK_WARN (0, "The routine excision_findboundary was called, but thorn " // CCTK_THORNSTRING // " is not active")
+  end if
 
 ! Loop over grid points.
 
@@ -43,7 +46,7 @@ subroutine excision_findboundary (ierr, mask, ni, nj, nk)
 
 !          Check if we are in an excised point.
 
-           if (mask(i,j,k)==MASK_EXCISED) then
+           if (abs(mask(i,j,k)-MASK_EXCISED)<0.01) then
 
               bnd = .false.
 
@@ -53,7 +56,7 @@ subroutine excision_findboundary (ierr, mask, ni, nj, nk)
               do kk=k-1,k+1
                  do jj=j-1,j+1
                     do ii=i-1,i+1
-                       bnd = bnd.or.mask(ii,jj,kk)==MASK_ACTIVE
+                       bnd = bnd.or.abs(mask(ii,jj,kk)-MASK_ACTIVE)<0.01
                     end do
                  end do
               end do
diff --git a/src/findnormals.F90 b/src/findnormals.F90
index ecd4d43..eb8e465 100644
--- a/src/findnormals.F90
+++ b/src/findnormals.F90
@@ -4,9 +4,8 @@
 ! Return information about the normals in these locations.
 
 #include "cctk.h"
-#include "cctk_Parameters.h"
-
 #include "cctk_Functions.h"
+#include "cctk_Parameters.h"
 
 #include "maskvalues.h"
 
@@ -39,6 +38,10 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
   CCTK_REAL sx,sy,sz,smag
   CCTK_REAL vx,vy,vz
   CCTK_REAL dot,dotmax
+  
+  if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then
+     call CCTK_WARN (0, "The routine excision_findnormals was called, but thorn " // CCTK_THORNSTRING // " is not active")
+  end if
 
   ! Initialise direction arrays to zero.
 
@@ -55,7 +58,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
            ! Check if the point is on the boundary, if it
            ! is not forget about it.
 
-           if (mask(i,j,k)==MASK_BOUNDARY) then
+           if (abs(mask(i,j,k)-MASK_BOUNDARY)<0.01) then
 
               ! Initialise the vector (sx,sy,sz) to zero.
 
@@ -72,7 +75,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
                        ! If a given neighbour is active, then add its
                        ! relative coordinates to (sx,sy,sz).
 
-                       if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then
+                       if (abs(mask(i+ii,j+jj,k+kk)-MASK_ACTIVE)<0.01) then
                           sx = sx + dble(ii)
                           sy = sy + dble(jj)
                           sz = sz + dble(kk)
@@ -92,6 +95,11 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
                  sx = sx/smag
                  sy = sy/smag
                  sz = sz/smag
+              else
+                 ! There is a tie.  Break it.
+                 sx = 1
+                 sy = 0
+                 sz = 0
               end if
 
               ! Now we have a unit vector pointing in the average
@@ -118,7 +126,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
 
                        ! If the point is not active, forget about it.
 
-                       if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then
+                       if (abs(mask(i+ii,j+jj,k+kk)-MASK_ACTIVE)<0.01) then
 
                           ! Find normalized dot product.
 
@@ -130,7 +138,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
                               / sqrt(vx**2 + vy**2 + vz**2)
 
 !                         If the new product is larger than the
-!                         largest so far, redifine the normal.
+!                         largest so far, redefine the normal.
 
                           if (dot > dotmax) then
 
@@ -155,7 +163,11 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
               jj = int(diry(i,j,k))
               kk = int(dirz(i,j,k))
 
-              if (mask(i+ii,j+jj,k+kk) /= MASK_ACTIVE) then
+              if (ii==0 .and. jj==0 .and. kk==0) then
+                 call CCTK_WARN (0, "Could not decide for a normal direction")
+              end if
+
+              if (abs(mask(i+ii,j+jj,k+kk)-MASK_ACTIVE)>0.01) then
                  call CCTK_WARN (0, "Mask boundary layer is too thick")
               end if
 
diff --git a/src/fixedsphere.F90 b/src/fixedsphere.F90
index 7d2294d..dafca7d 100644
--- a/src/fixedsphere.F90
+++ b/src/fixedsphere.F90
@@ -2,7 +2,6 @@
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
-
 #include "cctk_Functions.h"
 
 subroutine Lego_FixedSphere(CCTK_ARGUMENTS)
@@ -12,20 +11,30 @@ subroutine Lego_FixedSphere(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
+  
+  integer ierr
+
+  if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then
+     call CCTK_WARN (0, "The routine Lego_FixedSphere was called, but thorn " // CCTK_THORNSTRING // " is not active")
+  end if
+
+  ! Only excise if there would be more than a single grid point excised
+  if (fixed_size < 2*minval(CCTK_DELTA_SPACE(:))) return
 
   if (CCTK_EQUALS(fixed_excision,"sphere")) then
-    where (r < fixed_size)
+    where ((x - fixed_origin_x)**2 + (y - fixed_origin_y)**2 + (z - fixed_origin_z)**2 < fixed_size**2)
       emask = 0.d0
     elsewhere
       emask = 1.d0
     end where
   else if (CCTK_EQUALS(fixed_excision,"cube")) then
-    where ((abs(x) < fixed_size).and.(abs(z) < fixed_size).and.&
-         (abs(y) < fixed_size))
+    where (max(abs(x - fixed_origin_x), abs(y - fixed_origin_y), abs(z - fixed_origin_z)) < fixed_size)
       emask = 0.d0
     elsewhere
       emask = 1.d0
     end where
   end if
   
+  call excision_findboundary (ierr, emask, cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+  
 end subroutine Lego_FixedSphere