diff options
Diffstat (limited to 'src/findnormals.F90')
| -rw-r--r-- | src/findnormals.F90 | 26 |
1 files changed, 19 insertions, 7 deletions
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 |