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

2 files changed, 363 insertions, 238 deletions

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