Introduce a function SBP_determine_onesided_stencil which determines

which faces should use one-sided stencils, depending on which
boundaries are inter-processor boundaries, symmetry boundaries, and
multi-patch boundaries.

Use this function everywhere.  Remove the previous mechinisms; some
were not in all cases correct.


git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@75 f69c4107-0314-4c4f-9ad4-17e986b73f4a

11 files changed, 196 insertions, 193 deletions

diff --git a/src/GetScalProdDiag.F90 b/src/GetScalProdDiag.F90
index e2a7e99..00e3dff 100644
--- a/src/GetScalProdDiag.F90
+++ b/src/GetScalProdDiag.F90
@@ -19,12 +19,7 @@ subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad )
 
   CCTK_REAL, parameter :: zero = 0.0
   integer, parameter :: wp = kind(zero)
-  CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np
-  CCTK_INT, dimension(6) :: symbnd
-  CCTK_POINTER :: psym_name
-  character(len=256) :: symmetry_name
-  integer :: status
-  integer, dimension(6) :: bbox
+  integer, dimension(6) :: onesided
 
   CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 0.5_wp, 1.0_wp /) 
   CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 17.0_wp/48.0_wp, &
@@ -50,21 +45,11 @@ subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad )
     call CCTK_WARN(0, 'Error: dir is outside the legal range')
   end if
 
-  call CCTK_GroupbboxGN ( status, cctkGH, 6, bbox, &
-                                  'SummationByParts::normmask' )
-  if ( status < 0 ) then
-    call CCTK_WARN(0,'Error: unable to get bounding box information')
-  end if
-
-  symtable = SymmetryTableHandleForGrid ( cctkGH )
-  call Util_TableGetIntArray ( n_elements, symtable, 6, &
-                                           symbnd, 'symmetry_handle' )
+  call SBP_determine_onesided_stencil (cctkGH, onesided)
   
-  pen_sym_handle = SymmetryHandleOfName ( 'multipatch' )
-
   sigmad = 1.0_wp
 
-  if ( symbnd(dir+1) == pen_sym_handle .and. bbox(dir*2+1) == 1 ) then
+  if ( onesided(dir*2+1) == 1 ) then
     select case (order)
     case (2)
       sigmad(1:2) = bmask_2
@@ -76,7 +61,7 @@ subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad )
       sigmad(1:8) = bmask_8
     end select
   end if
-  if ( symbnd(dir+2) == pen_sym_handle .and. bbox(dir*2+2) == 1 ) then
+  if ( onesided(dir*2+2) == 1 ) then
     select case (order)
     case (2)
       sigmad(nsize:nsize-1:-1) = bmask_2
diff --git a/src/Get_Coeff.F90 b/src/Get_Coeff.F90
index 3c7a09f..60dab5b 100644
--- a/src/Get_Coeff.F90
+++ b/src/Get_Coeff.F90
@@ -22,6 +22,8 @@ CCTK_REAL function GetCoeff ()
       GetCoeff = 13649.0_wp/43200.0_wp
     case (8)
       GetCoeff =  1498139.0_wp/5080320.0_wp
+   case default
+      call CCTK_WARN (0, "operator not implemented")
     end select
   else
     if ( CCTK_EQUALS(operator_type,'Minimal Bandwidth') ) then
@@ -30,6 +32,8 @@ CCTK_REAL function GetCoeff ()
         GetCoeff = 3.0_wp/11.0_wp
       case(6)
         GetCoeff = 30.0_wp/137.0_wp
+     case default
+        call CCTK_WARN (0, "operator not implemented")
       end select
     else
       select case (order)
@@ -37,6 +41,8 @@ CCTK_REAL function GetCoeff ()
         GetCoeff = 0.2388575707774486064323157210922003533466_wp
       case(6)
         GetCoeff = 0.2028105550720356346665604029847379994496_wp
+     case default
+        call CCTK_WARN (0, "operator not implemented")
       end select
     end if
   end if
diff --git a/src/call_derivs.c b/src/call_derivs.c
index 5dd4530..22d6527 100644
--- a/src/call_derivs.c
+++ b/src/call_derivs.c
@@ -5,28 +5,7 @@
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
 
-
-
-/* Determine whether a boundary with the symmetry handle symbnd is a
-   "regular" symmetry boundary (where the SBP stencils should not be
-   modified), or an outer boundary (or a multi-block boundary), where
-   the SBP stencils need to be modified.  */
-static int is_regular_symbnd (int const symbnd, int const pen_sym_handle)
-{
-  /* On an outer boundary (which is not a symmetry boundary), symbnd <
-     0.  */
-  if (symbnd < 0) return 0;
-  
-  if (pen_sym_handle >= 0) {
-    /* If the symmetry boundary is a multi-block boundary, then symbnd
-     = pen_sym_handle.  However, we can only check that if thorn
-     MultiPatch is active, i.e., if pen_sym_handle >= 0.  */
-    if (symbnd == pen_sym_handle) return 0;
-  }
-  
-  /* This is a "regular" symmetry boundary.  */
-  return 1;
-}
+#include "stencil.h"
 
 
 
@@ -38,12 +17,10 @@ void DiffGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_ARGUMENTS
 
-  CCTK_INT ni, nj, nk, gsize, ic, loc_order,i;
+  CCTK_INT ni, nj, nk, gsize, loc_order;
   CCTK_REAL delta;
-  CCTK_INT ierr;
-  CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3];
-  CCTK_INT symtable, pen_sym_handle;
-  CCTK_INT symbnd[6];
+  CCTK_INT bb[2];
+  int onesided[6];
 
   int nelements;
   void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(const CCTK_REAL *var,
@@ -151,37 +128,27 @@ void DiffGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
     loc_order = order;
   }
 
-  symtable = SymmetryTableHandleForGrid (cctkGH);
-  if (symtable<0) {
-    CCTK_WARN(0,"symtable is out of bounds");
-  }
-
-  ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle");
-  if (ierr!=6) {
-    CCTK_WARN(0,"Util_TableGetIntArray returned error");
-  }
-
-  pen_sym_handle = SymmetryHandleOfName ( "multipatch" );
+  SBP_determine_onesided_stencil (cctkGH, onesided);
 
   switch(dir) {
   case 0: {
     delta = CCTK_DELTA_SPACE(0);
-    bb[0] = cctk_bbox[0] && ! is_regular_symbnd (symbnd[0], pen_sym_handle);
-    bb[1] = cctk_bbox[1] && ! is_regular_symbnd (symbnd[1], pen_sym_handle);
+    bb[0] = onesided[0];
+    bb[1] = onesided[1];
     gsize = cctk_nghostzones[0];
     break;
   }
   case 1: {
     delta = CCTK_DELTA_SPACE(1);
-    bb[0] = cctk_bbox[2] && ! is_regular_symbnd (symbnd[2], pen_sym_handle);
-    bb[1] = cctk_bbox[3] && ! is_regular_symbnd (symbnd[3], pen_sym_handle);
+    bb[0] = onesided[2];
+    bb[1] = onesided[3];
     gsize = cctk_nghostzones[1];
     break;
   }
   case 2: {
     delta = CCTK_DELTA_SPACE(2);
-    bb[0] = cctk_bbox[4] && ! is_regular_symbnd (symbnd[4], pen_sym_handle);
-    bb[1] = cctk_bbox[5] && ! is_regular_symbnd (symbnd[5], pen_sym_handle);
+    bb[0] = onesided[4];
+    bb[1] = onesided[5];
     gsize = cctk_nghostzones[2];
     break;
   }
diff --git a/src/call_derivs_name.c b/src/call_derivs_name.c
index f7804b0..89e4b39 100644
--- a/src/call_derivs_name.c
+++ b/src/call_derivs_name.c
@@ -1,8 +1,12 @@
+#include <assert.h>
+
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
-#include <assert.h>
+#include "stencil.h"
+
+
 
 void DiffGf ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
               const char *var_name, const char *dvar_name )
@@ -12,10 +16,10 @@ void DiffGf ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   DECLARE_CCTK_ARGUMENTS
 
   CCTK_REAL *var, *dvar;
-  CCTK_INT ni, nj, nk, gsize, ic;
+  CCTK_INT ni, nj, nk, gsize;
   CCTK_REAL delta;
-  CCTK_INT ierr;
-  CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3];
+  CCTK_INT bb[2];
+  int onesided[6];
   void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(const CCTK_REAL *var,
                                        const CCTK_INT *ni,
                                        const CCTK_INT *nj,
@@ -65,22 +69,24 @@ void DiffGf ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   
   ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
   
+  SBP_determine_onesided_stencil (cctkGH, onesided);
+
   switch(dir) {
   case 0: {
     delta = CCTK_DELTA_SPACE(0);
-    bb[0] = cctk_bbox[0]; bb[1] = cctk_bbox[1];
+    bb[0] = onesided[0]; bb[1] = onesided[1];
     gsize = cctk_nghostzones[0];
     break;
   }
   case 1: {
     delta = CCTK_DELTA_SPACE(1);
-    bb[0] = cctk_bbox[2]; bb[1] = cctk_bbox[3];
+    bb[0] = onesided[2]; bb[1] = onesided[3];
     gsize = cctk_nghostzones[1];
     break;
   }
   case 2: {
     delta = CCTK_DELTA_SPACE(2);
-    bb[0] = cctk_bbox[4]; bb[1] = cctk_bbox[5];
+    bb[0] = onesided[4]; bb[1] = onesided[5];
     gsize = cctk_nghostzones[2];
     break;
   }
diff --git a/src/call_up_derivs.c b/src/call_up_derivs.c
index e6f5afa..10607a1 100644
--- a/src/call_up_derivs.c
+++ b/src/call_up_derivs.c
@@ -1,3 +1,5 @@
+#include <assert.h>
+
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
@@ -5,6 +7,10 @@
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
 
+#include "stencil.h"
+
+
+
 void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
               const CCTK_REAL *var, CCTK_REAL *dvar, const CCTK_REAL *up, 
               const CCTK_INT table_handle )
@@ -16,9 +22,8 @@ void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   CCTK_INT ni, nj, nk, gsize, ic, loc_order;
   CCTK_REAL delta;
   CCTK_INT ierr;
-  CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3];
-  CCTK_INT symtable, pen_sym_handle;
-  CCTK_INT symbnd[6];
+  CCTK_INT lsh[3], bb[2], nghostzones[3];
+  int onesided[6];
 
   int nelements;
   void CCTK_FCALL CCTK_FNAME(up_deriv_gf_2_1)(const CCTK_REAL *var,
@@ -76,61 +81,27 @@ void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
     loc_order = order;
   }
 
-  symtable = SymmetryTableHandleForGrid (cctkGH);
-  if (symtable<0) {
-    CCTK_WARN(0,"symtable is out of bounds");
-  }
-
-  ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle");
-  if (ierr!=6) {
-    CCTK_WARN(0,"Util_TableGetIntArray returned error");
-  }
-
-  pen_sym_handle = SymmetryHandleOfName ( "multipatch" );
-
+  SBP_determine_onesided_stencil (cctkGH, onesided);
+  
   switch(dir) {
   case 0: {
     delta = CCTK_DELTA_SPACE(0);
-    if ( symbnd[0] == pen_sym_handle && cctk_bbox[0] == 1 ) {
-      bb[0] = 1;
-    } else {
-      bb[0] = 0;
-    }
-    if ( symbnd[1] == pen_sym_handle && cctk_bbox[1] == 1 ) {
-      bb[1] = 1;
-    } else {
-      bb[1] = 0;
-    }
+    bb[0] = onesided[0];
+    bb[1] = onesided[1];
     gsize = cctk_nghostzones[0];
     break;
   }
   case 1: {
     delta = CCTK_DELTA_SPACE(1);
-    if ( symbnd[2] == pen_sym_handle && cctk_bbox[2] == 1) {
-      bb[0] = 1;
-    } else {
-      bb[0] = 0;
-    }
-    if ( symbnd[3] == pen_sym_handle && cctk_bbox[3] == 1) {
-      bb[1] = 1;
-    } else {
-      bb[1] = 0;
-    }
+    bb[0] = onesided[2];
+    bb[1] = onesided[3];
     gsize = cctk_nghostzones[1];
     break;
   }
   case 2: {
     delta = CCTK_DELTA_SPACE(2);
-    if ( symbnd[4] == pen_sym_handle && cctk_bbox[4] == 1) {
-      bb[0] = 1;
-    } else {
-      bb[0] = 0;
-    }
-    if ( symbnd[4] == pen_sym_handle && cctk_bbox[5] == 1) {
-      bb[1] = 1;
-    } else {
-      bb[1] = 0;
-    }
+    bb[0] = onesided[4];
+    bb[1] = onesided[5];
     gsize = cctk_nghostzones[2];
     break;
   }
diff --git a/src/dissipation.c b/src/dissipation.c
index 8f85d32..bf6fc29 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -11,6 +11,8 @@
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
 
+#include "stencil.h"
+
 void CCTK_FCALL CCTK_FNAME(dissipation_4_3_opt) (const CCTK_REAL *var,
                                              const CCTK_INT *lsh,
                                              const CCTK_INT *gsh,
@@ -155,13 +157,12 @@ apply (int const varindex, char const * const optstring, void * const arg)
   cGroup vardata, rhsdata;
   CCTK_REAL const * varptr;
   CCTK_REAL       * rhsptr;
-  CCTK_INT ni, nj, nk, bbox[6], gsize[3];
+  CCTK_INT gsize[3];
   CCTK_REAL dx[3];
-  int n;
   int d;
   int ierr;
-  CCTK_INT symtable, pen_sym_handle;
-  CCTK_INT symbnd[6];
+  CCTK_INT bbox[6];
+  int onesided[6];
   CCTK_INT npatches, patch;
   
   assert (varindex >= 0);
@@ -171,24 +172,9 @@ apply (int const varindex, char const * const optstring, void * const arg)
     gsize[d] = cctk_nghostzones[d];
   }
   
-  symtable = SymmetryTableHandleForGrid (cctkGH);
-  if (symtable<0) {
-    CCTK_WARN(0,"symtable is out of bounds");
-  }
-
-  ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle");
-  if (ierr!=6) {
-    CCTK_WARN(0,"Util_TableGetIntArray returned error");
-  }
-
-  pen_sym_handle = SymmetryHandleOfName ( "multipatch" );
-
+  SBP_determine_onesided_stencil (cctkGH, onesided);
   for (d=0; d<6; ++d) {
-    if ( symbnd[d] == pen_sym_handle && cctk_bbox[d] == 1 ) {
-      bbox[d] = 1;
-    } else {
-      bbox[d] = 0;
-    }
+    bbox[d] = onesided[d];
   }
 
   rhsindex = MoLQueryEvolvedRHS (varindex);
diff --git a/src/get_coeffs.c b/src/get_coeffs.c
index 3d441a3..c2aa224 100644
--- a/src/get_coeffs.c
+++ b/src/get_coeffs.c
@@ -8,6 +8,55 @@
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
 
+#include "stencil.h"
+
+
+
+void CCTK_FCALL CCTK_FNAME(set_coeff_2_1)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
+void CCTK_FCALL CCTK_FNAME(set_coeff_4_2)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
+void CCTK_FCALL CCTK_FNAME(set_coeff_6_3)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
+void CCTK_FCALL CCTK_FNAME(set_coeff_8_4)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
+void CCTK_FCALL CCTK_FNAME(set_coeff_4_3)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
+void CCTK_FCALL CCTK_FNAME(set_coeff_6_5)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
+
+
+
 void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
                  const CCTK_INT nsize, CCTK_INT *imin, CCTK_INT *imax,
                  CCTK_REAL *q, const CCTK_INT table_handle  )
@@ -19,49 +68,7 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   CCTK_INT ni, nj, nk, gsize, loc_order;
   int nelements;
   CCTK_INT bb[2];
-  void CCTK_FCALL CCTK_FNAME(set_coeff_2_1)(const CCTK_INT *nsize,
-                                            const CCTK_INT *loc_order,
-                                            const CCTK_INT *bb,
-                                            const CCTK_INT *gsize,
-                                            CCTK_INT *imin,
-                                            CCTK_INT *imax,
-                                            CCTK_REAL *q);
-  void CCTK_FCALL CCTK_FNAME(set_coeff_4_2)(const CCTK_INT *nsize,
-                                            const CCTK_INT *loc_order,
-                                            const CCTK_INT *bb,
-                                            const CCTK_INT *gsize,
-                                            CCTK_INT *imin,
-                                            CCTK_INT *imax,
-                                            CCTK_REAL *q);
-  void CCTK_FCALL CCTK_FNAME(set_coeff_6_3)(const CCTK_INT *nsize,
-                                            const CCTK_INT *loc_order,
-                                            const CCTK_INT *bb,
-                                            const CCTK_INT *gsize,
-                                            CCTK_INT *imin,
-                                            CCTK_INT *imax,
-                                            CCTK_REAL *q);
-  void CCTK_FCALL CCTK_FNAME(set_coeff_8_4)(const CCTK_INT *nsize,
-                                            const CCTK_INT *loc_order,
-                                            const CCTK_INT *bb,
-                                            const CCTK_INT *gsize,
-                                            CCTK_INT *imin,
-                                            CCTK_INT *imax,
-                                            CCTK_REAL *q);
-  void CCTK_FCALL CCTK_FNAME(set_coeff_4_3)(const CCTK_INT *nsize,
-                                            const CCTK_INT *loc_order,
-                                            const CCTK_INT *bb,
-                                            const CCTK_INT *gsize,
-                                            CCTK_INT *imin,
-                                            CCTK_INT *imax,
-                                            CCTK_REAL *q);
-  void CCTK_FCALL CCTK_FNAME(set_coeff_6_5)(const CCTK_INT *nsize,
-                                            const CCTK_INT *loc_order,
-                                            const CCTK_INT *bb,
-                                            const CCTK_INT *gsize,
-                                            CCTK_INT *imin,
-                                            CCTK_INT *imax,
-                                            CCTK_REAL *q);
-
+  int onesided[6];
   ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
 
   if ( table_handle >=0 ) {
@@ -75,22 +82,24 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
     loc_order = order;
   }
 
+  SBP_determine_onesided_stencil (cctkGH, onesided);
+
   switch(dir) {
   case 0: {
     assert(nsize==ni);
-    bb[0] = cctk_bbox[0]; bb[1] = cctk_bbox[1];
+    bb[0] = onesided[0]; bb[1] = onesided[1];
     gsize = cctk_nghostzones[0];
     break;
   }
   case 1: {
     assert(nsize==nj);
-    bb[0] = cctk_bbox[2]; bb[1] = cctk_bbox[3];
+    bb[0] = onesided[2]; bb[1] = onesided[3];
     gsize = cctk_nghostzones[1];
     break;
   }
   case 2: {
     assert(nsize==nk);
-    bb[0] = cctk_bbox[4]; bb[1] = cctk_bbox[5];
+    bb[0] = onesided[4]; bb[1] = onesided[5];
     gsize = cctk_nghostzones[2];
     break;
   }
diff --git a/src/make.code.defn b/src/make.code.defn
index dd2255b..b9672ed 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -39,7 +39,8 @@ SRCS =  call_derivs.c \
         Coefficients_6_3.F90 \
         Coefficients_8_4.F90 \
         Coefficients_4_3.F90 \
-        Coefficients_6_5.F90
+        Coefficients_6_5.F90 \
+        stencil.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/set_norm_mask.F90 b/src/set_norm_mask.F90
index c9b8ffb..40bae3e 100644
--- a/src/set_norm_mask.F90
+++ b/src/set_norm_mask.F90
@@ -17,8 +17,6 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
   integer, parameter :: wp = kind(zero)
   CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np
   CCTK_INT, dimension(6) :: symbnd
-  CCTK_POINTER :: psym_name
-  character(len=256) :: symmetry_name
   CCTK_REAL, dimension(:), allocatable :: mask_x, mask_y, mask_z
 
 ! Note: The first number is twice the value from the paper, since Carpet
@@ -55,7 +53,7 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
   
   pen_sym_handle = SymmetryHandleOfName ( 'multipatch' )
 
-  if ( any ( symbnd == pen_sym_handle ) ) then
+  if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. any ( symbnd == pen_sym_handle )) ) then
     allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
     mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0
 
@@ -79,22 +77,22 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
       end select
     end if
 
-    if ( symbnd(1) == pen_sym_handle ) then
-        mask_x(1:np) = bmask(1:np)
+    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(1) == pen_sym_handle) ) then
+      mask_x(1:np) = bmask(1:np)
     end if
-    if ( symbnd(2) == pen_sym_handle ) then
+    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(2) == pen_sym_handle) ) then
       mask_x(ni:ni-np+1:-1) = bmask(1:np)
     end if
-    if ( symbnd(3) == pen_sym_handle ) then
+    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(3) == pen_sym_handle) ) then
       mask_y(1:np) = bmask(1:np)
     end if
-    if ( symbnd(4) == pen_sym_handle ) then
+    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(4) == pen_sym_handle) ) then
       mask_y(nj:nj-np+1:-1) = bmask(1:np)
     end if
-    if ( symbnd(5) == pen_sym_handle ) then
+    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(5) == pen_sym_handle) ) then
       mask_z(1:np) = bmask(1:np)
     end if
-    if ( symbnd(6) == pen_sym_handle ) then
+    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(6) == pen_sym_handle) ) then
       mask_z(nk:nk-np+1:-1) = bmask(1:np)
     end if
 
diff --git a/src/stencil.c b/src/stencil.c
new file mode 100644
index 0000000..007a7c6
--- /dev/null
+++ b/src/stencil.c
@@ -0,0 +1,66 @@
+#include "cctk.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "stencil.h"
+
+
+
+/* Determine whether a boundary with the symmetry handle symbnd is a
+   "regular" symmetry boundary (where the SBP stencils should not be
+   modified), or an outer boundary (or a multi-block boundary), where
+   the SBP stencils need to be modified.  */
+
+void SBP_determine_onesided_stencil (const cGH * cctkGH, int * onesided)
+{
+  int symtable;
+  int pen_sym_handle;
+  CCTK_INT symbnd[6];
+  int ierr;
+  int d;
+  
+  symtable = SymmetryTableHandleForGrid (cctkGH);
+  if (symtable<0) {
+    CCTK_WARN(0,"Cannot get symmetry table handle -- maybe thorn SymBase is not active?");
+  }
+  
+  ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle");
+  if (ierr!=6) {
+    CCTK_WARN(0,"Cannot get symmetry handles");
+  }
+  
+  pen_sym_handle = SymmetryHandleOfName ( "multipatch" );
+  
+  for (d=0; d<6; ++d) {
+    if (! cctkGH->cctk_bbox[d]) {
+      /* This is an inter-processor boundary */
+      onesided[d] = 0;
+    } else {
+      /* On an outer boundary (which is not a symmetry boundary),
+         symbnd < 0.  */
+      if (symbnd[d] < 0) {
+        onesided[d] = 1;
+      } else {
+        if (pen_sym_handle >= 0) {
+          /* If the symmetry boundary is a multi-block boundary, then
+             symbnd = pen_sym_handle.  However, we can only check that
+             thorn MultiPatch is active, i.e., whether pen_sym_handle
+             >= 0.  */
+          if (symbnd[d] == pen_sym_handle) {
+            onesided[d] = 1;
+          } else {
+            onesided[d] = 0;
+          }
+        } else {
+          onesided[d] = 0;
+        }
+      }
+    }
+  }
+}
+
+CCTK_FCALL void CCTK_FNAME (SBP_determine_onesided_stencil) (CCTK_POINTER_TO_CONST cctkGH, int * onesided)
+{
+  SBP_determine_onesided_stencil (cctkGH, onesided);
+}
diff --git a/src/stencil.h b/src/stencil.h
new file mode 100644
index 0000000..7937e1e
--- /dev/null
+++ b/src/stencil.h
@@ -0,0 +1,8 @@
+#ifndef STENCIL_H
+#define STENCIL_H
+
+#include "cctk.h"
+
+void SBP_determine_onesided_stencil (const cGH * cctkGH, int * onesided);
+
+#endif /* #ifndef STENCIL_H */