Added routine for identifying the boundary (shamelessly copied from

Erik's BoundaryExample). This is probably a bad place to put this, and
function aliasing is probably better than a cross-thorn function call.
But this is easier...

3 files changed, 144 insertions, 2 deletions

diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/interface.ccl b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/interface.ccl
index b1637b4..9bb205c 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/interface.ccl
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/interface.ccl
@@ -8,3 +8,17 @@ implements: GenericFD
 
 INCLUDE HEADER: GenericFD.h in GenericFD.h
 
+#USES INCLUDE: Boundary.h
+#USES INCLUDE: Symmetry.h
+
+CCTK_INT FUNCTION \
+    GetBoundarySpecification (CCTK_INT IN size,                  \
+                              CCTK_INT OUT ARRAY nboundaryzones, \
+                              CCTK_INT OUT ARRAY is_internal,    \
+                              CCTK_INT OUT ARRAY is_staggered,   \
+                              CCTK_INT OUT ARRAY shiftout)
+REQUIRES FUNCTION GetBoundarySpecification
+
+CCTK_INT FUNCTION \
+    SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
+REQUIRES FUNCTION SymmetryTableHandleForGrid
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
index 3750cf1..3c9dc6d 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
@@ -31,9 +31,11 @@
 #include "cctk.h"                             
 #include "cctk_Arguments.h"                   
 #include "cctk_Parameters.h"                  
-                                                
+#include "util_Table.h"
+            
 #include "Symmetry.h"                         
 
+
 #define KRANC_C
                                                 
 #include "GenericFD.h"
@@ -51,3 +53,113 @@ CCTK_INT sgn(CCTK_REAL x)
   else 
     return 0;
 }
+
+/* Return the array indices in imin and imax of each boundary. The
+   indexing is C-style. Also return whether the boundary is a symmetry */
+void GenericFD_GetBoundaryInfo(cGH *cctkGH, CCTK_INT *cctk_lsh, CCTK_INT *cctk_bbox,
+			       CCTK_INT *cctk_nghostzones, CCTK_INT *imin, 
+			       CCTK_INT *imax, CCTK_INT *is_symbnd, 
+			       CCTK_INT *is_physbnd, CCTK_INT *is_ipbnd)
+{
+  CCTK_INT nboundaryzones[6];
+  CCTK_INT is_internal[6];
+  CCTK_INT is_staggered[6];
+  CCTK_INT shiftout[6];
+  CCTK_INT symbnd[6];
+
+  CCTK_INT symtable = 0;
+  CCTK_INT dir = 0;
+  CCTK_INT face = 0;
+  CCTK_INT npoints = 0;
+  CCTK_INT iret = 0;
+  CCTK_INT ierr = 0;
+
+  ierr = GetBoundarySpecification(6, nboundaryzones, is_internal, is_staggered, 
+				  shiftout);
+  if (ierr != 0)
+    CCTK_WARN(0, "Could not obtain boundary specification");
+
+  symtable = SymmetryTableHandleForGrid(cctkGH);
+  if (symtable < 0) 
+  {
+    CCTK_WARN(0, "Could not obtain symmetry table");
+  }  
+  
+  iret = Util_TableGetIntArray(symtable, 6, symbnd, "symmetry_handle");
+  if (iret != 6) CCTK_WARN (0, "Could not obtain symmetry information");
+
+  for (dir = 0; dir < 6; dir++)
+  {
+    is_symbnd[dir] = 0;
+    is_physbnd[dir] = 0;
+    is_ipbnd[dir] = 0;
+
+    if (symbnd[dir] >= 0)
+    {
+      is_symbnd[dir] = 1;
+    }
+    else
+    {
+      is_symbnd[dir] = 0;
+    }
+
+    if (cctk_bbox[dir] == 0)
+    {
+      is_ipbnd[dir] = 1;
+    }
+    else
+    {
+      is_ipbnd[dir] = 0;
+    }
+    
+    if (!is_ipbnd[dir] && !is_symbnd[dir])
+    {
+      is_physbnd[dir] = 1;
+    }
+    else
+    {
+      is_physbnd[dir] = 0;
+    }
+  }
+
+  for (dir = 0; dir < 3; dir++)
+  {
+    for (face = 0; face < 2; face++)
+    {
+      CCTK_INT index = dir*2 + face;
+      if (is_ipbnd[index])
+      {
+	/* Inter-processor boundary */
+	npoints = cctk_nghostzones[dir];
+      }
+      else
+      {
+	/* Symmetry or physical boundary */
+	npoints = nboundaryzones[index];
+             
+	if (is_symbnd[index])
+	{
+	  /* Ensure that the number of symmetry zones is the same
+	     as the number of ghost zones */
+	  if (npoints != cctk_nghostzones[dir])
+	  {
+	    CCTK_WARN (1, "The number of symmetry points is different from the number of ghost points; this is probably an error");
+	  }
+	}
+      }
+
+      switch(face)
+      {
+      case 0: /* Lower boundary */
+	imin[dir] = npoints - 1;
+	break;
+      case 1: /* Upper boundary */
+	imax[dir] = cctk_lsh[dir]  - npoints;
+	break;
+      default:
+	CCTK_WARN(0, "internal error");
+      }
+    }
+  }
+}
+
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
index 2acd0ff..6f10fea 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
@@ -550,7 +550,23 @@
 
 #endif
 
-
 #ifdef KRANC_C
 CCTK_INT sgn(CCTK_REAL x);
+
+#define Dupwind1(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i+dir,j,k)] \
+	 - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
+#define Dupwind2(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i,j+dir,k)] \
+	 - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
+#define Dupwind3(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+dir)] \
+	 - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
+
+void GenericFD_GetBoundaryInfo(cGH *cctkGH, CCTK_INT *cctk_lsh, CCTK_INT *cctk_bbox,
+			       CCTK_INT *cctk_nghostzones, CCTK_INT *imin, 
+			       CCTK_INT *imax, CCTK_INT *is_symbnd, 
+			       CCTK_INT *is_physbnd, CCTK_INT *is_ipbound);
+
 #endif
+
+
+
+