GenericFD: Add function to determine the number of boundary points

2 files changed, 53 insertions, 0 deletions

diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
index 8514d66..5964df0 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
@@ -36,6 +36,7 @@
 #include <assert.h>
 #include <stdlib.h>
 #include <math.h>
+#include <stdio.h>
             
 #include "Symmetry.h"                         
 
@@ -58,6 +59,54 @@ int sgn(CCTK_REAL x)
     return 0;
 }
 
+int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH)
+{
+  int is_internal[6];
+  int is_staggered[6];
+  int nboundaryzones[6];
+  int shiftout[6];
+  int ierr = -1;
+
+  if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) {
+    int const map = MultiPatch_GetMap (cctkGH);
+    if (map < 0)
+      CCTK_WARN(0, "Could not obtain boundary specification");
+    ierr = MultiPatch_GetBoundarySpecification
+      (map, 6, nboundaryzones, is_internal, is_staggered, shiftout);
+    if (ierr != 0)
+      CCTK_WARN(0, "Could not obtain boundary specification");
+  } else if (CCTK_IsFunctionAliased ("GetBoundarySpecification")) {
+    ierr = GetBoundarySpecification
+      (6, nboundaryzones, is_internal, is_staggered, shiftout);
+    if (ierr != 0)
+      CCTK_WARN(0, "Could not obtain boundary specification");
+  } else {
+    CCTK_WARN(0, "Could not obtain boundary specification");
+  }
+
+  int bw = nboundaryzones[0];
+
+  for (int i = 1; i < 6; i++)
+    if (nboundaryzones[i] != bw)
+    CCTK_WARN(0, "Number of boundary points is different on different faces");
+
+  return bw;
+}
+
+/* int GenericFD_BoundaryWidthTable(cGH const * restrict const cctkGH) */
+/* { */
+/*   int nboundaryzones[6]; */
+/*   GenericFD_GetBoundaryWidth(cctkGH, nboundaryzones); */
+
+/*   int table = Util_TableCreate(0); */
+/*   if (table < 0) CCTK_WARN(0, "Could not create table"); */
+
+/*   if (Util_TableSetIntArray(table, 6, nboundaryzones, "BOUNDARY_WIDTH") < 0) */
+/*     CCTK_WARN(0, "Could not set table"); */
+/*   return table; */
+/* } */
+
+
 /* Return the array indices in imin and imax for looping over the
    interior of the grid. imin is the index of the first grid point.
    imax is the index of the last grid point plus 1.  So a loop over
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
index db74360..041347d 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
@@ -673,6 +673,10 @@ int sgn(CCTK_REAL x);
 #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)
 
+int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH);
+
+/* int GenericFD_BoundaryWidthTable(cGH const * restrict const cctkGH); */
+
 void GenericFD_GetBoundaryInfo(cGH const * restrict cctkGH,
                                int const * restrict cctk_lsh,
                                int const * restrict cctk_lssh,