From fe9c96feb0b2281c44b9bf8b39438a4e2744d381 Mon Sep 17 00:00:00 2001 From: miguel Date: Mon, 2 Apr 2001 15:31:14 +0000 Subject: New findboundary routine without templates. It seems to work, but requires more testing yet. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/LegoExcision/trunk@14 f75ba9e5-694f-0410-ac2c-87ea7ce7132b --- src/findnormals.F90 | 291 ++++++++------------------------------------ src/findnormals_erik.F90 | 310 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 363 insertions(+), 238 deletions(-) create mode 100644 src/findnormals_erik.F90 (limited to 'src') diff --git a/src/findnormals.F90 b/src/findnormals.F90 index 2936424..3d6a7d0 100644 --- a/src/findnormals.F90 +++ b/src/findnormals.F90 @@ -23,261 +23,76 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) ! in: array sizes for grid functions ! (you can pass in cctk_lsh(:) for these) integer :: ni,nj,nk - + ! in: mask CCTK_REAL :: mask(ni,nj,nk) ! out: normal directions to use for interpolation CCTK_REAL :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk) - - ! distinguish between faces, edges, and corners. - - ! look at the excised points, and distinguish the following cases, - ! in this order: - - ! [trivial case] - ! 1: if (none) - ! case 0: interior - ! [box shape cases] - ! 2: if (only one corner) - ! case I: corner - ! 3: if (only points on one edge) - ! case II: edge - ! 4: if (only points on one face) - ! case III: face - ! [convex cases] - ! 5: if (only points on two adjacent faces) - ! case IIa: anti-edge - ! 6: if (one points on three adjacent faces) - ! case Ia: anti-corner - ! [illegal cases] - ! error - - logical mask_interior (-1:+1,-1:+1,-1:+1, 1) - logical mask_corner (-1:+1,-1:+1,-1:+1, 8) - logical mask_edge (-1:+1,-1:+1,-1:+1, 12) - logical mask_face (-1:+1,-1:+1,-1:+1, 6) - logical mask_antiedge (-1:+1,-1:+1,-1:+1, 12) - logical mask_anticorner(-1:+1,-1:+1,-1:+1, 8) - integer dir_interior (3, 1), num_interior - integer dir_corner (3, 8), num_corner - integer dir_edge (3, 12), num_edge - integer dir_face (3, 6), num_face - integer dir_antiedge (3, 12), num_antiedge - integer dir_anticorner (3, 8), num_anticorner - - logical msk(-1:+1,-1:+1,-1:+1) - - integer di,dj,dk - integer i,j,k - integer n - integer dir(3) - - character(1000) message - - - - ! - ! Create templates - ! - num_interior = 0 - num_corner = 0 - num_edge = 0 - num_face = 0 - num_antiedge = 0 - num_anticorner = 0 - - do dk=-1,+1 - do dj=-1,+1 - do di=-1,+1 - - select case (count( (/di,dj,dk/) /= 0 )) - - case (0) - ! interior (or error, which is not explicitely checked) - num_interior = num_interior + 1 - dir_interior(:,num_interior) = (/di,dj,dk/) - do k=-1,+1 - do j=-1,+1 - do i=-1,+1 - mask_interior(i,j,k,num_interior) = & - any( (/di,dj,dk/)/=0 & - .and. (/i,j,k/)==(/di,dj,dk/) ) - end do - end do - end do - - case (1) - ! face (or anti-face, which is the same thing) - num_face = num_face + 1 - dir_face(:,num_face) = (/di,dj,dk/) - do k=-1,+1 - do j=-1,+1 - do i=-1,+1 - mask_face(i,j,k,num_face) = & - any( (/di,dj,dk/)/=0 & - .and. (/i,j,k/)==(/di,dj,dk/) ) - end do - end do - end do - - case (2) - ! edge or anti-edge - num_edge = num_edge + 1 - dir_edge(:,num_edge) = (/di,dj,dk/) - do k=-1,+1 - do j=-1,+1 - do i=-1,+1 - mask_edge(i,j,k,num_edge) = & - all( (/di,dj,dk/)==0 & - .or. (/i,j,k/)==(/di,dj,dk/) ) - end do - end do - end do - - num_antiedge = num_antiedge + 1 - dir_antiedge(:,num_antiedge) = (/di,dj,dk/) - do k=-1,+1 - do j=-1,+1 - do i=-1,+1 - mask_antiedge(i,j,k,num_antiedge) = & - any( (/di,dj,dk/)/=0 & - .and. (/i,j,k/)==(/di,dj,dk/) ) - end do - end do - end do - - case (3) - ! corner or anti-corner - num_corner = num_corner + 1 - dir_corner(:,num_corner) = (/di,dj,dk/) - do k=-1,+1 - do j=-1,+1 - do i=-1,+1 - mask_corner(i,j,k,num_corner) = & - all( (/di,dj,dk/)==0 & - .or. (/i,j,k/)==(/di,dj,dk/) ) - end do - end do - end do - - num_anticorner = num_anticorner + 1 - dir_anticorner(:,num_anticorner) = (/di,dj,dk/) - do k=-1,+1 - do j=-1,+1 - do i=-1,+1 - mask_anticorner(i,j,k,num_anticorner) = & - any( (/di,dj,dk/)/=0 & - .and. (/i,j,k/)==(/di,dj,dk/) ) - end do - end do - end do - - case default - call CCTK_WARN (0, "internal error") - end select - - end do - end do - end do - - ! sanity check - if ( num_interior /= 1 & - .or. num_corner /= 8 & - .or. num_edge /= 12 & - .or. num_face /= 6 & - .or. num_antiedge /= 12 & - .or. num_anticorner /= 8 ) then - call CCTK_WARN (0, "internal error") - end if - - - - ! - ! Calculate normal directions - ! - +! Extra variables. + integer i,j,k,ii,jj,kk + integer sx,sy,sz,smax + +! Initialise direction arrays to zero. + dirx = 0 diry = 0 dirz = 0 - + +! Loop over 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 (mask(i,j,k)==MASK_BOUNDARY) then - - ! make a snapshot - msk = mask(i-1:i+1,j-1:j+1,k-1:k+1) /= MASK_EXCISED - - ! for a template to match, all grid points have to be - ! either active, or set in the template. - search: do - - ! 1: case 0: interior - do n=1,num_interior - if (all(msk .or. mask_interior(:,:,:,n))) then - dir = dir_interior(:,n) - exit search - end if - end do - - ! 2: case I: corner - do n=1,num_corner - if (all(msk .or. mask_corner(:,:,:,n))) then - dir = dir_corner(:,n) - exit search - end if - end do - - ! 3: case II: edge - do n=1,num_edge - if (all(msk .or. mask_edge(:,:,:,n))) then - dir = dir_edge(:,n) - exit search - end if - end do - - ! 4: case III: face - do n=1,num_face - if (all(msk .or. mask_face(:,:,:,n))) then - dir = dir_face(:,n) - exit search - end if - end do - - ! 5: case IIa: anti-edge - do n=1,num_antiedge - if (all(msk .or. mask_antiedge(:,:,:,n))) then - dir = dir_antiedge(:,n) - exit search - end if - end do - - ! 6: case Ia: anticorner - do n=1,num_anticorner - if (all(msk .or. mask_anticorner(:,:,:,n))) then - dir = dir_anticorner(:,n) - exit search - end if + +! Initialise integer vector (sx,sy,sz) to zero. + + sx = 0; sy = 0; sz = 0 + +! Loop around 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 (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then + sx = sx + ii; sy = sy + jj; sz = sz + kk + end if + + end do end do - - ! error: unclassified - call CCTK_WARN (0, "Mask is not convex") - - end do search - - dir(:) = -dir(:) - - if (mask(i+dir(1), j+dir(2), k+dir(3)) /= MASK_ACTIVE) then + end do + +! Find maximum component of (sx,sy,sz). + + smax = max(abs(sx),abs(sy),abs(sz)) + +! 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 = sx/smax + sy = sy/smax + sz = sz/smax + + dirx(i,j,k) = sx + diry(i,j,k) = sy + dirz(i,j,k) = sz + +! Sanity check. + + if (mask(i+sx,j+sy,k+sz) /= MASK_ACTIVE) then call CCTK_WARN (0, "Mask boundary layer is too thick") end if - - dirx(i,j,k) = dir(1) - diry(i,j,k) = dir(2) - dirz(i,j,k) = dir(3) - + end if end do diff --git a/src/findnormals_erik.F90 b/src/findnormals_erik.F90 new file mode 100644 index 0000000..8feccf2 --- /dev/null +++ b/src/findnormals_erik.F90 @@ -0,0 +1,310 @@ +! $Header$ + +! Take a mask where the boundary has been been marked. +! Return information about the normals in these locations. + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "maskvalues.h" + +subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + ! Arguments + + ! out: zero for success, nonzero for error + integer :: ierr + + ! in: array sizes for grid functions + ! (you can pass in cctk_lsh(:) for these) + integer :: ni,nj,nk + integer :: ii,jj,kk + + ! in: mask + CCTK_REAL :: mask(ni,nj,nk) + + ! out: normal directions to use for interpolation + CCTK_REAL :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk) + + + ! distinguish between faces, edges, and corners. + + ! look at the excised points, and distinguish the following cases, + ! in this order: + + ! [trivial case] + ! 1: if (none) + ! case 0: interior + ! [box shape cases] + ! 2: if (only one corner) + ! case I: corner + ! 3: if (only points on one edge) + ! case II: edge + ! 4: if (only points on one face) + ! case III: face + ! 5: if (only points on two adjacent faces) + ! case IIa: anti-edge + ! 6: if (one points on three adjacent faces) + ! case Ia: anti-corner + ! [illegal cases] + ! error + + logical mask_interior (-1:+1,-1:+1,-1:+1, 1) + logical mask_corner (-1:+1,-1:+1,-1:+1, 8) + logical mask_edge (-1:+1,-1:+1,-1:+1, 12) + logical mask_face (-1:+1,-1:+1,-1:+1, 6) + logical mask_antiface (-1:+1,-1:+1,-1:+1, 6) + logical mask_antiedge (-1:+1,-1:+1,-1:+1, 12) + logical mask_anticorner(-1:+1,-1:+1,-1:+1, 8) + integer dir_interior (3, 1), num_interior + integer dir_corner (3, 8), num_corner + integer dir_edge (3, 12), num_edge + integer dir_face (3, 6), num_face + integer dir_antiface (3, 6), num_antiface + integer dir_antiedge (3, 12), num_antiedge + integer dir_anticorner (3, 8), num_anticorner + + logical msk(-1:+1,-1:+1,-1:+1) + + integer di,dj,dk + integer i,j,k + integer n + integer dir(3) + + character(1000) message + + + + ! + ! Create templates + ! + num_interior = 0 + num_corner = 0 + num_edge = 0 + num_face = 0 + num_antiedge = 0 + num_anticorner = 0 + + do dk=-1,+1 + do dj=-1,+1 + do di=-1,+1 + + select case (count( (/di,dj,dk/) /= 0 )) + + case (0) + ! interior (or error, which is not explicitely checked) + num_interior = num_interior + 1 + dir_interior(:,num_interior) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_interior(i,j,k,num_interior) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case (1) + ! face (or anti-face, which is the same thing) + num_face = num_face + 1 + dir_face(:,num_face) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_face(i,j,k,num_face) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case (2) + ! edge or anti-edge + num_edge = num_edge + 1 + dir_edge(:,num_edge) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_edge(i,j,k,num_edge) = & + all( (/di,dj,dk/)==0 & + .or. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + num_antiedge = num_antiedge + 1 + dir_antiedge(:,num_antiedge) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_antiedge(i,j,k,num_antiedge) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case (3) + ! corner or anti-corner + num_corner = num_corner + 1 + dir_corner(:,num_corner) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_corner(i,j,k,num_corner) = & + all( (/di,dj,dk/)==0 & + .or. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + num_anticorner = num_anticorner + 1 + dir_anticorner(:,num_anticorner) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_anticorner(i,j,k,num_anticorner) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case default + call CCTK_WARN (0, "internal error") + end select + + end do + end do + end do + + ! sanity check + if ( num_interior /= 1 & + .or. num_corner /= 8 & + .or. num_edge /= 12 & + .or. num_face /= 6 & + .or. num_antiedge /= 12 & + .or. num_anticorner /= 8 ) then + call CCTK_WARN (0, "internal error") + end if + + + + ! + ! Calculate normal directions + ! + + dirx = 0 + diry = 0 + dirz = 0 + + do k=2,nk-1 + do j=2,nj-1 + do i=2,ni-1 + + if (mask(i,j,k)==MASK_BOUNDARY) then + + ! make a snapshot + msk = mask(i-1:i+1,j-1:j+1,k-1:k+1) == MASK_ACTIVE + + ! for a template to match, all grid points have to be + ! either active, or set in the template. + + search: do + + ! 1: case 0: interior + do n=1,num_interior + if (all(msk .eqv. mask_interior(:,:,:,n))) then + dir = dir_interior(:,n) + exit search + end if + end do + + ! 2: case I: corner + do n=1,num_corner + if (all(msk .eqv. mask_corner(:,:,:,n))) then + dir = dir_corner(:,n) + exit search + end if + end do + + ! 3: case II: edge + do n=1,num_edge + if (all(msk .eqv. mask_edge(:,:,:,n))) then + dir = dir_edge(:,n) + exit search + end if + end do + + ! 4: case III: face + do n=1,num_face + if (all(msk .eqv. mask_face(:,:,:,n))) then + dir = dir_face(:,n) + exit search + end if + end do + + ! 5: case IIa: antiedge + do n=1,num_antiedge + if (all(msk .eqv. mask_antiedge(:,:,:,n))) then + dir = dir_antiedge(:,n) + exit search + end if + end do + + ! 6: case Ia: anticorner + do n=1,num_anticorner + if (all(msk .eqv. mask_anticorner(:,:,:,n))) then + dir = dir_anticorner(:,n) + exit search + end if + end do + + ! error: unclassified + print *,i,j,k,dir(1),dir(2),dir(3) + do kk=k-1,k+1 + do jj=j-1,j+1 + do ii=i-1,i+1 + print*,ii-i,jj-j,kk-k,msk(ii-i,jj-j,kk-k) + end do + end do + end do + call CCTK_WARN (0, "Mask is not convex") + + end do search + + ! dir(:) = -dir(:) + + if (mask(i+dir(1), j+dir(2), k+dir(3)) /= MASK_ACTIVE) then + print *,i,j,k,dir(1),dir(2),dir(3) + do kk=k-1,k+1 + do jj=j-1,j+1 + do ii=i-1,i+1 + print*,ii-i,jj-j,kk-k,msk(ii-i,jj-j,kk-k) + end do + end do + end do + call CCTK_WARN (0, "Mask boundary layer is too thick") + end if + + dirx(i,j,k) = dir(1) + diry(i,j,k) = dir(2) + dirz(i,j,k) = dir(3) + + end if + + end do + end do + end do + + ! Success + + ierr = 0 + +end subroutine excision_findnormals -- cgit v1.2.3