diff options
Diffstat (limited to 'src/findnormals.F90')
| -rw-r--r-- | src/findnormals.F90 | 290 |
1 files changed, 290 insertions, 0 deletions
diff --git a/src/findnormals.F90 b/src/findnormals.F90 new file mode 100644 index 0000000..df0ae69 --- /dev/null +++ b/src/findnormals.F90 @@ -0,0 +1,290 @@ +! $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, cctkgh, 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: pointer to Cactus GH + CCTK_POINTER :: cctkgh + + ! 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 + ! + + 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) /= 0 + + ! 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 + end do + + ! error: unclassified + call CCTK_WARN (0, "Mask is not convex") + + end do search + + dir(:) = -dir(:) + + 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 |