Changed the way in which the normals where found. Before, I just chose

that grid point closest to the optimal normal.  But it turned out that
in some cases this grid point was in fact masked.  I now choose the
closest unmasked grid point to the optimal normal.  This means a second
loop over nearest neighbours, but it can't be avoided.


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

1 files changed, 79 insertions, 33 deletions

diff --git a/src/findnormals.F90 b/src/findnormals.F90
index e9aebfe..4c963cf 100644
--- a/src/findnormals.F90
+++ b/src/findnormals.F90
@@ -30,41 +30,45 @@ 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)
 
-! Internal variables.
+  ! Internal variables.
 
   integer i,j,k,ii,jj,kk
 
-  CCTK_REAL sx,sy,sz,smax
+  CCTK_REAL sx,sy,sz,smag
+  CCTK_REAL vx,vy,vz
+  CCTK_REAL dot,dotmax
 
-! Initialise direction arrays to zero.
+  ! Initialise direction arrays to zero.
 
   dirx = 0
   diry = 0
   dirz = 0
 
-! Loop over all 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, if it
-!          is not forget about ir.
+           ! Check if the point is on the boundary, if it
+           ! is not forget about it.
 
            if (mask(i,j,k)==MASK_BOUNDARY) then
 
-!             Initialise the vector (sx,sy,sz) to zero.
+              ! Initialise the vector (sx,sy,sz) to zero.
 
-              sx = 0; sy = 0; sz = 0
+              sx = 0.0D0
+              sy = 0.0D0
+              sz = 0.0D0
 
-!             Loop around all nearest neighbours.
+              ! Loop around nearest neighbours.
 
               do kk=-1,1
                  do jj=-1,1
                     do ii=-1,1
 
-!                      If a given neighbour is active, then 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 + dble(ii)
@@ -76,38 +80,80 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
                  end do
               end do
 
-!             Find the maximum component of (sx,sy,sz).
+              ! Find magnitude of (sx,sy,sz).
 
-              smax = max(dabs(sx),dabs(sy),dabs(sz))
+              smag = dsqrt(sx**2 + sy**2 + sz**2)
 
-!             Divide (sx,sy,sz) with smax.
+              ! Normalize (sx,sy,sz).
 
-              if (smax /= 0.0D0) then
-                 sx = sx/smax
-                 sy = sy/smax
-                 sz = sz/smax
+              if (smag /= 0.0D0) then
+                 sx = sx/smag
+                 sy = sy/smag
+                 sz = sz/smag
               end if
 
-!             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.
+              ! Now we have a unit vector pointing in the average
+              ! direction on the unmasked neighbours.  This is
+              ! as close to the normal direction as we can get.
 
-              if (sx /= 0.0D0) sx = dble(int(dabs(sx)+0.5D0))*sx/dabs(sx)
-              if (sy /= 0.0D0) sy = dble(int(dabs(sy)+0.5D0))*sy/dabs(sy)
-              if (sz /= 0.0D0) sz = dble(int(dabs(sz)+0.5D0))*sz/dabs(sz)
+              ! What we need to do now is find that unmasked
+              ! neighbour that is closest to this direction.
+              ! For this I loop again over unmasked neighbours,
+              ! and find the normalized dot product between
+              ! (sx,sy,sz) and the direction of a given neighbour.
+              ! The largest dot product will give us our best
+              ! normal direction.
 
-!             Copy the vector (sx,sy,sz) into the normal arrays.
+              dotmax = 0.0D0
 
-              dirx(i,j,k) = sx
-              diry(i,j,k) = sy
-              dirz(i,j,k) = sz
+              dirx(i,j,k) = 0.0D0
+              diry(i,j,k) = 0.0D0
+              dirz(i,j,k) = 0.0D0
 
-!             Sanity check.  If the check fails then the mask really
-!             has a crazy shape.
+              do kk=-1,1
+                 do jj=-1,1
+                    do ii=-1,1
+
+                       ! If the point is not active, forget about it.
+
+                       if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then
+
+                          ! Find normalized dot product.
+
+                          vx = dble(ii)
+                          vy = dble(jj)
+                          vz = dble(kk)
+
+                          dot = (sx*vx + sy*vy + sz*vz) &
+                              / dsqrt(vx**2 + vy**2 + vz**2)
+
+!                         If the new product is larger than the
+!                         largest so far, redifine the normal.
+
+                          if (dot > dotmax) then
+
+                             dotmax = dot
+
+                             dirx(i,j,k) = vx
+                             diry(i,j,k) = vy
+                             dirz(i,j,k) = vz
+
+                          end if
+
+                       end if
+
+                    end do
+                 end do
+              end do
+
+              ! Sanity check.  If the check fails then
+              ! something is very wrong.
+
+              ii = int(dirx(i,j,k))
+              jj = int(diry(i,j,k))
+              kk = int(dirz(i,j,k))
 
-              if (mask(i+int(sx),j+int(sy),k+int(sz)) /= MASK_ACTIVE) then
+              if (mask(i+ii,j+jj,k+kk) /= MASK_ACTIVE) then
                  call CCTK_WARN (0, "Mask boundary layer is too thick")
               end if