Changed the way in which the normal was found so that when the "real"

normal vector is close to a diagonal, we get in fact the diagnal and
not a vector along the coordinate lines.


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

1 files changed, 37 insertions, 19 deletions

diff --git a/src/findnormals.F90 b/src/findnormals.F90
index 3d6a7d0..85f56c1 100644
--- a/src/findnormals.F90
+++ b/src/findnormals.F90
@@ -30,9 +30,11 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
   ! out: normal directions to use for interpolation
   CCTK_REAL    :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk)
 
-! Extra variables.
+! Internal variables.
+
   integer i,j,k,ii,jj,kk
-  integer sx,sy,sz,smax
+
+  CCTK_REAL sx,sy,sz,smax
 
 ! Initialise direction arrays to zero.
 
@@ -40,56 +42,72 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
   diry = 0
   dirz = 0
 
-! Loop over grid points.
+! Loop over all grid points.
 
   do k=2,nk-1
      do j=2,nj-1
         do i=2,ni-1
 
-!          Check if the point is in the boundary.
+!          Check if the point is in the boundary, if it
+!          is not forget about ir.
 
            if (mask(i,j,k)==MASK_BOUNDARY) then
 
-!             Initialise integer vector (sx,sy,sz) to zero.
+!             Initialise the vector (sx,sy,sz) to zero.
 
               sx = 0; sy = 0; sz = 0
 
-!             Loop around nearest neighbours.
+!             Loop around all nearest neighbours.
 
               do kk=-1,1
                  do jj=-1,1
                     do ii=-1,1
 
-!                      If a neighbour is active, add its relative
-!                      coordinates to (sx,sy,sz).
+!                      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
-                          sx = sx + ii; sy = sy + jj; sz = sz + kk
+                          sx = sx + dble(ii)
+                          sy = sy + dble(jj)
+                          sz = sz + dble(kk)
                        end if
 
                     end do
                  end do
               end do
 
-!             Find maximum component of (sx,sy,sz).
+!             Find the maximum component of (sx,sy,sz).
+
+              smax = max(dabs(sx),dabs(sy),dabs(sz))
+
+!             Divide (sx,sy,sz) with smax.
+
+              if (smax /= 0.0D0) then
+                 sx = sx/smax
+                 sy = sy/smax
+                 sz = sz/smax
+              end if
 
-              smax = max(abs(sx),abs(sy),abs(sz))
+!             Now we have a real vector where the largest component
+!             has absolute value 1.  We want to change this into
+!             the closest vector with components given by 0 or +-1.
+!             To do this, I add 1/2 to all components (with the
+!             correct sign), ang get the integer part.
 
-!             Do an integer division of (sx,sy,sz) with smax.
-!             This ensures that we end up with a vector that
-!             has only 0 and +-1 as components.
+              sx = dble(int(dabs(sx)+0.5D0))*sx/dabs(sx)
+              sy = dble(int(dabs(sy)+0.5D0))*sy/dabs(sy)
+              sz = dble(int(dabs(sz)+0.5D0))*sz/dabs(sz)
 
-              sx = sx/smax
-              sy = sy/smax
-              sz = sz/smax
+!             Copy the vector (sx,sy,sz) into the normal arrays.
 
               dirx(i,j,k) = sx
               diry(i,j,k) = sy
               dirz(i,j,k) = sz
 
-!             Sanity check.
+!             Sanity check.  If the check fails then the mask really
+!             has a crazy shape.
 
-              if (mask(i+sx,j+sy,k+sz) /= MASK_ACTIVE) then
+              if (mask(i+int(sx),j+int(sy),k+int(sz)) /= MASK_ACTIVE) then
                  call CCTK_WARN (0, "Mask boundary layer is too thick")
               end if