! $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