From 66151de369a1cf743574ac03471f259d3a898eaf Mon Sep 17 00:00:00 2001 From: miguel Date: Fri, 27 Apr 2001 11:36:43 +0000 Subject: 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 --- src/findnormals.F90 | 112 ++++++++++++++++++++++++++++++++++++---------------- 1 file 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 -- cgit v1.2.3