Adding Erik's excision routine.

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

2 files changed, 294 insertions, 1 deletions

diff --git a/src/excision.F90 b/src/excision.F90
new file mode 100644
index 0000000..2b7a670
--- /dev/null
+++ b/src/excision.F90
@@ -0,0 +1,293 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+subroutine excision_boundary (ierr, cctkgh, var, oldvar, mask)
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  
+  ! Arguments
+  integer,      intent(out)   :: ierr
+  CCTK_POINTER, intent(in)    :: cctkgh
+  CCTK_REAL,    intent(inout) :: var(:,:,:)
+  CCTK_REAL,    intent(in)    :: oldvar(:,:,:)
+  CCTK_REAL,    intent(in)    :: mask(:,:,:)
+  
+  ! 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
+  
+  
+  
+  !
+  ! Check the mask
+  !
+  do k=1,size(mask,3)
+     do j=1,size(mask,2)
+        do i=1,size(mask,1)
+           if (mask(i,j,k)/=0 .and. mask(i,j,k)/=1) then
+              write (message, '("Mask contains illegal value at (",i4,",",i4,",",i4,")")') i,j,k
+              call CCTK_WARN (0, trim(message))
+           end if
+        end do
+     end do
+  end do
+  
+  
+  
+  !
+  ! 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
+  ! Extrapolate
+  !
+  do k=2,size(mask,3)-1
+     do j=2,size(mask,2)-1
+        do i=2,size(mask,1)-1
+           
+           if (mask(i,j,k)==0 .and. any(mask(i-1:i+1,j-1:j+1,k-1:k+1)/=0)) then
+
+              !
+              ! Calculate normal direction
+              !
+              
+              ! 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(:)
+              
+              
+              
+              !
+              ! Extrapolate
+              !
+              var(i,j,k) = oldvar(i,j,k) + var(i+dir(1), j+dir(2), k+dir(3)) &
+                                         - oldvar(i+dir(1), j+dir(2), k+dir(3))
+              
+           end if
+           
+        end do
+     end do
+  end do
+  
+  ! Success
+  ierr = 0
+  
+end subroutine excision_boundary
diff --git a/src/make.code.defn b/src/make.code.defn
index 573e79f..af22b60 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = 
+SRCS = excision.F90
 
 # Subdirectories containing source files
 SUBDIRS =