From 78ba61b65fb29a4b0e0217fc51fbdc59c6fb9487 Mon Sep 17 00:00:00 2001 From: miguel Date: Tue, 3 Apr 2001 15:15:15 +0000 Subject: 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 --- src/findnormals.F90 | 56 +++++++++++++++++++++++++++++++++++------------------ 1 file 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 -- cgit v1.2.3