Removed the original file that handled extrapolation directly.

There are now three files for three purposes:

findboundary takes a mask with only 0s and 1s and marks the boundary
with 0.5s.

findnormals takes a mask where the boundary has been marked and
returns the normal directions to the mask for these locations.

extrapolate takes a mask the normals and extrapolates a grid function.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/LegoExcision/trunk@6 f75ba9e5-694f-0410-ac2c-87ea7ce7132b

5 files changed, 418 insertions, 1 deletions

diff --git a/src/extrapolate.F90 b/src/extrapolate.F90
new file mode 100644
index 0000000..b9b0b99
--- /dev/null
+++ b/src/extrapolate.F90
@@ -0,0 +1,62 @@
+! $Header$
+
+! Extrapolate the difference between the grid function var and oldvar to
+! var where indicated by mask.  Use the normal direction given by dirx,
+! diry, and dirz for extrapolation.
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "maskvalues.h"
+
+subroutine excision_extrapolate (ierr, cctkgh, var, oldvar, &
+     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
+  
+  ! inout: grid function that should be interpolated
+  CCTK_REAL    :: var(ni,nj,nk)
+  
+  ! in: other grid function for interpolation
+  CCTK_REAL    :: oldvar(ni,nj,nk)
+  
+  ! in: mask
+  CCTK_REAL    :: mask(ni,nj,nk)
+  
+  ! in: normal directions to use for interpolation
+  CCTK_REAL    :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk)
+  
+  integer i,j,k
+  integer ii,jj,kk
+  
+  do k=2,nk-1
+     do j=2,nj-1
+        do i=2,ni-1
+           if (mask(i,j,k)==MASK_BOUNDARY) then
+              ii = i + dirx(i,j,k)
+              jj = j + diry(i,j,k)
+              kk = k + dirz(i,j,k)
+              var(i,j,k) = oldvar(i,j,k) + var(ii,jj,kk) - oldvar(ii,jj,kk)
+           end if
+        end do
+     end do
+  end do
+  
+  ierr = 0
+  
+end subroutine excision_extrapolate
diff --git a/src/findboundary.F90 b/src/findboundary.F90
new file mode 100644
index 0000000..1e16e97
--- /dev/null
+++ b/src/findboundary.F90
@@ -0,0 +1,57 @@
+! $Header$
+
+! Take a mask that contains only the values 0 and 1.
+! Return a mask where the outermost 0s have been replaced by 0.5s.
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "maskvalues.h"
+
+subroutine excision_findboundary (ierr, cctkgh, mask, 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
+  
+  ! inout: mask
+  CCTK_REAL    :: mask(ni,nj,nk)
+  
+  integer i,j,k
+  integer ii,jj,kk
+  logical bnd
+  
+  do k=2,nk-1
+     do j=2,nj-1
+        do i=2,ni-1
+           if (mask(i,j,k)==MASK_EXCISED) then
+              bnd = .false.
+              do kk=k-1,k+1
+                 do jj=j-1,j+1
+                    do ii=i-1,i+1
+                       bnd = bnd .or. mask(i,j,k)/=MASK_EXCISED
+                    end do
+                 end do
+              end do
+              if (bnd) mask(i,j,k) = MASK_BOUNDARY
+           end if
+        end do
+     end do
+  end do
+  
+  ierr = 0
+  
+end subroutine excision_findboundary
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
diff --git a/src/make.code.defn b/src/make.code.defn
index af22b60..854a5e8 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = excision.F90
+SRCS = extrapolate.F90 findboundary.F90 findnormals.F90
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/maskvalues.h b/src/maskvalues.h
new file mode 100644
index 0000000..c2822f2
--- /dev/null
+++ b/src/maskvalues.h
@@ -0,0 +1,8 @@
+/* -*-f90-*- */
+/* $Header$ */
+
+/* Values of the mask */
+
+#define MASK_EXCISED  0
+#define MASK_BOUNDARY 0.5
+#define MASK_ACTIVE   1