path: root/src/apply.c

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "util_ErrorCodes.h"
#include "util_Table.h"

#include "reflection.h"



#define COPY_PRE(VARTYPE)                                               \
  static void                                                           \
  copy_##VARTYPE (VARTYPE const * restrict const srcvar,                \
                  VARTYPE * restrict const dstvar,                      \
                  int const ni, int const nj, int const nk,             \
                  int const imin, int const jmin, int const kmin,       \
                  int const imax, int const jmax, int const kmax,       \
                  int const ioff, int const joff, int const koff,       \
                  int const idir, int const jdir, int const kdir,       \
                  int const parity)                                     \
  {                                                                     \
    assert (abs(idir)==1);                                              \
    assert (abs(jdir)==1);                                              \
    assert (abs(kdir)==1);                                              \
    int const iioff = ioff + (1 - idir) * imin;                         \
    int const jjoff = joff + (1 - jdir) * jmin;                         \
    int const kkoff = koff + (1 - kdir) * kmin;                         \
    int const iimin = iioff + idir * imin;                              \
    int const jjmin = jjoff + jdir * jmin;                              \
    int const kkmin = kkoff + kdir * kmin;                              \
    int const iimax = iioff + idir * imax;                              \
    int const jjmax = jjoff + jdir * jmax;                              \
    int const kkmax = kkoff + kdir * kmax;                              \
    assert (imin>=0 && imax<=ni);                                       \
    assert (jmin>=0 && jmax<=nj);                                       \
    assert (kmin>=0 && kmax<=nk);                                       \
    assert (iimin>=0 && iimax<=ni);                                     \
    assert (jjmin>=0 && jjmax<=nj);                                     \
    assert (kkmin>=0 && kkmax<=nk);                                     \
    assert (iimax>=-1 && iimin<ni);                                     \
    assert (jjmax>=-1 && jjmin<nj);                                     \
    assert (kkmax>=-1 && kkmin<nk);

#define COPY_LOOP(VARTYPE)                                              \
    for (int k=kmin; k<kmax; ++k) {                                     \
      for (int j=jmin; j<jmax; ++j) {                                   \
        for (int i=imin; i<imax; ++i) {                                 \
          int const dstind = i + ni * (j + nj * k);                     \
          int const ii = iioff + idir * i;                              \
          int const jj = jjoff + jdir * j;                              \
          int const kk = kkoff + kdir * k;                              \
          int const srcind = ii + ni * (jj + nj * kk);                  \
          REAL(dstvar[dstind] RE = parity * srcvar[srcind] RE;)         \
          IMAG(dstvar[dstind] IM = parity * srcvar[srcind] IM;)         \
        }                                                               \
      }                                                                 \
    }                                                                   \

#define COPY_POST(VARTYPE)                      \
  }

#define REAL(x) x
#define IMAG(x) /* nothing */
#define RE /* nothing */
#define IM ERROR ERROR ERROR

#ifdef HAVE_CCTK_INT1
COPY_PRE(CCTK_INT1)
#pragma omp parallel for
COPY_LOOP(CCTK_INT1)
COPY_POST(CCTK_INT1)
#endif

#ifdef HAVE_CCTK_INT2
COPY_PRE(CCTK_INT2)
#pragma omp parallel for
COPY_LOOP(CCTK_INT2)
COPY_POST(CCTK_INT2)
#endif

#ifdef HAVE_CCTK_INT4
COPY_PRE(CCTK_INT4)
#pragma omp parallel for
COPY_LOOP(CCTK_INT4)
COPY_POST(CCTK_INT4)
#endif

#ifdef HAVE_CCTK_INT8
COPY_PRE(CCTK_INT8)
#pragma omp parallel for
COPY_LOOP(CCTK_INT8)
COPY_POST(CCTK_INT8)
#endif

#ifdef HAVE_CCTK_REAL4
COPY_PRE(CCTK_REAL4)
#pragma omp parallel for
COPY_LOOP(CCTK_REAL4)
COPY_POST(CCTK_REAL4)
#endif

#ifdef HAVE_CCTK_REAL8
COPY_PRE(CCTK_REAL8)
#pragma omp parallel for
COPY_LOOP(CCTK_REAL8)
COPY_POST(CCTK_REAL8)
#endif

#ifdef HAVE_CCTK_REAL16
COPY_PRE(CCTK_REAL16)
#pragma omp parallel for
COPY_LOOP(CCTK_REAL16)
COPY_POST(CCTK_REAL16)
#endif

#undef REAL
#undef IMAG
#undef RE
#undef IM

#define REAL(x) x
#define IMAG(x) x
#define RE .Re
#define IM .Im

#ifdef HAVE_CCTK_COMPLEX8
COPY_PRE(CCTK_COMPLEX8)
#pragma omp parallel for
COPY_LOOP(CCTK_COMPLEX8)
COPY_POST(CCTK_COMPLEX8)
#endif

#ifdef HAVE_CCTK_COMPLEX16
COPY_PRE(CCTK_COMPLEX16)
#pragma omp parallel for
COPY_LOOP(CCTK_COMPLEX16)
COPY_POST(CCTK_COMPLEX16)
#endif

#ifdef HAVE_CCTK_COMPLEX32
COPY_PRE(CCTK_COMPLEX32)
#pragma omp parallel for
COPY_LOOP(CCTK_COMPLEX32)
COPY_POST(CCTK_COMPLEX32)
#endif

#undef REAL
#undef IMAG
#undef RE
#undef IM

#undef COPY_PRE
#undef COPY_LOOP
#undef COPY_POST



static int
BndReflectVI (cGH const * restrict const cctkGH,
              int const vi)
{
  DECLARE_CCTK_PARAMETERS;
  
  int gi;
  cGroup group;
  cGroupDynamicData data;
  int firstvar, numvars;
  char * restrict fullname;
  
  void * restrict varptr;
  
  int table;
  char tensortypealias[1000];
  enum tensortype { UNKNOWN,
                    SCALAR, VECTOR, SYMTENSOR, SYMTENSOR3, TENSOR,
                    WEYLSCALARS_REAL, MANUALCARTESIAN };
  enum tensortype ttype;
  CCTK_INT tensorparity;
  int tcomponent;
  
  int do_reflection[6];
  int do_stagger[6];
  
  int dir, face;
  
  int lsh[3], imin[3], imax[3], ioff[3], idir[3];
  
  int parity;
  int manual_parities[3];

  int d;
  
  int ierr;
  
  
  
  /* Check arguments */
  if (! cctkGH) {
    CCTK_WARN (0, "Argument cctkGH is NULL");
  }
  if (vi < 0 || vi >= CCTK_NumVars()) {
    CCTK_WARN (0, "Illegal variable index");
  }
  
  if (verbose) {
    fullname = CCTK_FullName (vi);
    if (! fullname)   {
      CCTK_WARN (0, "Internal error in CCTK_FullName");
    }
    CCTK_VInfo (CCTK_THORNSTRING,
                "Applying reflection boundary conditions to \"%s\"",
                fullname);
    free (fullname);
  }
  
  
  
  /* Get and check group information */
  gi = CCTK_GroupIndexFromVarI (vi);
  if (gi < 0 || gi > CCTK_NumGroups()) {
    CCTK_WARN (0, "Internal error in CCTK_GroupIndexFromVarI");
  }
  
  ierr = CCTK_GroupData (gi, &group);
  assert (!ierr);
  assert (group.grouptype == CCTK_GF);
  assert (group.disttype == CCTK_DISTRIB_DEFAULT);
  assert (group.stagtype == 0);
  
  firstvar = CCTK_FirstVarIndexI (gi);
  assert (firstvar>=0 && firstvar<CCTK_NumVars());
  numvars = CCTK_NumVarsInGroupI (gi);
  assert (numvars>=0);
  
  ierr = CCTK_GroupDynamicData (cctkGH, gi, &data);
  assert (!ierr);
  
  table = CCTK_GroupTagsTableI(gi);
  assert (table>=0);
  
  varptr = CCTK_VarDataPtrI (cctkGH, 0, vi);
  assert (varptr);
  
  
  
  /* Get and check tensor type information */
  ierr = Util_TableGetString
    (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
    /* assume a scalar */
    if (numvars != 1) {
      static int * restrict didwarn = 0;
      if (! didwarn) {
        didwarn = calloc (CCTK_NumGroups(), sizeof *didwarn);
      }
      if (! didwarn[gi]) {
        didwarn[gi] = 1;
        {
          char * groupname = CCTK_GroupName(gi);
          assert (groupname);
          CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                      "Group \"%s\" has no tensor type and contains more than one element -- treating these as \"scalar\"",
                      groupname);
          free (groupname);
        }
      }
    }
    strcpy (tensortypealias, "scalar");
  } else if (ierr<0) {
    char * groupname = CCTK_GroupName(gi);
    assert (groupname);
    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Error in tensor type alias declaration for group \"%s\"",
                groupname);
    free (groupname);
  }
  
  ttype = UNKNOWN;
  tcomponent = 0;
  if (CCTK_EQUALS (tensortypealias, "scalar")) {
    /* scalar */
    ttype = SCALAR;
    tcomponent = 0;
  } else if (CCTK_EQUALS (tensortypealias, "4scalar")) {
    /* 4-scalar */
    ttype = SCALAR;
    tcomponent = 0;
  } else if (CCTK_EQUALS (tensortypealias, "u")
             || CCTK_EQUALS (tensortypealias, "d"))
  {
    /* vector */
    assert (numvars == 3);
    ttype = VECTOR;
    tcomponent = vi - firstvar;
  } else if (CCTK_EQUALS (tensortypealias, "4u")
             || CCTK_EQUALS (tensortypealias, "4d"))
  {
    /* 4-vector */
    assert (numvars == 4);
    if (vi == firstvar) {
      ttype = SCALAR;
      tcomponent = 0;
    } else {
      ttype = VECTOR;
      tcomponent = vi - firstvar - 1;
    }
  } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
             || CCTK_EQUALS (tensortypealias, "dd_sym"))
  {
    /* symmetric tensor */
    assert (numvars == 6);
    ttype = SYMTENSOR;
    tcomponent = vi - firstvar;
  } else if (CCTK_EQUALS (tensortypealias, "uu")
             || CCTK_EQUALS (tensortypealias, "ud")
             || CCTK_EQUALS (tensortypealias, "du")
             || CCTK_EQUALS (tensortypealias, "dd"))
  {
    /* non-symmetric tensor */
    assert (numvars == 9);
    ttype = TENSOR;
    tcomponent = vi - firstvar;
  } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
             || CCTK_EQUALS (tensortypealias, "4dd_sym"))
  {
    /* symmetric 4-tensor */
    assert (numvars == 10);
    if (vi == firstvar) {
      ttype = SCALAR;
      tcomponent = 0;
    } else if (vi <= firstvar+3) {
      ttype = VECTOR;
      tcomponent = vi - firstvar - 1;
    } else {
      ttype = SYMTENSOR;
      tcomponent = vi - firstvar - 4;
    }
  } else if (CCTK_EQUALS (tensortypealias, "ddd_sym")) {
    /* 3rd rank tensor, symmetric in last 2 indices */
    assert (numvars == 18);
    ttype = SYMTENSOR3;
    tcomponent = vi - firstvar;
  } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
    /* Weyl scalars, stored as 10 real values.  NOTE: This assumes
       that Psi_0 comes first, which is NOT the default with
       PsiKadelia.  */
    assert (numvars == 10);
    ttype = WEYLSCALARS_REAL;
    tcomponent = vi - firstvar;
  } else if (CCTK_EQUALS (tensortypealias, "ManualCartesian")) {
      /* Reflection symmetries specified by hand */
      ttype = MANUALCARTESIAN;
      tcomponent = vi - firstvar;
  } else {
    char * groupname = CCTK_GroupName(gi);
    assert (groupname);
    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Illegal tensor type alias for group \"%s\"",
                groupname);
    free (groupname);
  }
  
  switch (ttype) {
  case SCALAR:
    assert (tcomponent>=0 && tcomponent<1);
    break;
  case VECTOR:
    assert (tcomponent>=0 && tcomponent<3);
    break;
  case SYMTENSOR:
    assert (tcomponent>=0 && tcomponent<6);
    break;
  case SYMTENSOR3:
    assert (tcomponent>=0 && tcomponent<18);
    break;
  case TENSOR:
    assert (tcomponent>=0 && tcomponent<9);
    break;
  case WEYLSCALARS_REAL:
    assert (tcomponent>=0 && tcomponent<10);
    break;
  case MANUALCARTESIAN:
    /* No restriction on number of components */
    ReflectionSymmetry_GetManualParities(table, gi, manual_parities);
    break;

  default:
    assert (0);
  }
  
  ierr = Util_TableGetInt (table, & tensorparity, "tensorparity");
  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
    tensorparity = +1;
  } else if (ierr<0) {
    char * groupname = CCTK_GroupName(gi);
    assert (groupname);
    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Error in tensor parity declaration for group \"%s\"",
                groupname);
    free (groupname);
  }
  
  
  
  /* Reflection symmetry information */
  do_reflection[0] = reflection_x;
  do_reflection[1] = reflection_upper_x;
  do_reflection[2] = reflection_y;
  do_reflection[3] = reflection_upper_y;
  do_reflection[4] = reflection_z;
  do_reflection[5] = reflection_upper_z;
  
  do_stagger[0] = avoid_origin_x;
  do_stagger[1] = avoid_origin_upper_x;
  do_stagger[2] = avoid_origin_y;
  do_stagger[3] = avoid_origin_upper_y;
  do_stagger[4] = avoid_origin_z;
  do_stagger[5] = avoid_origin_upper_z;
  
  
  
  /* Loop over all directions and faces */
  for (dir=0; dir<3; ++dir) {
    for (face=0; face<2; ++face) {
      /* If there is a reflection symmetry on that face */
      if (do_reflection[2*dir+face]) {
        /* If we have the outer boundary of that face */
        if (cctkGH->cctk_bbox[2*dir+face]) {
          
          /* Find parity */
          parity = tensorparity;
          switch (ttype) {
          case SCALAR:
            parity *= +1;
            break;
          case VECTOR:
            parity *= dir == tcomponent ? -1 : +1;
            break;
          case SYMTENSOR:
            switch (tcomponent) {
            case 0: parity *= +1; break;
            case 1: parity *= (dir == 0 || dir == 1) ? -1 : +1; break;
            case 2: parity *= (dir == 0 || dir == 2) ? -1 : +1; break;
            case 3: parity *= +1; break;
            case 4: parity *= (dir == 1 || dir == 2) ? -1 : +1; break;
            case 5: parity *= +1; break;
            default: assert (0);
            }
            break;
          case SYMTENSOR3:
            switch (tcomponent % 6) {
            case 0: parity *= +1; break;
            case 1: parity *= (dir == 0 || dir == 1) ? -1 : +1; break;
            case 2: parity *= (dir == 0 || dir == 2) ? -1 : +1; break;
            case 3: parity *= +1; break;
            case 4: parity *= (dir == 1 || dir == 2) ? -1 : +1; break;
            case 5: parity *= +1; break;
            default: assert (0);
            }
            switch (tcomponent / 6) {
            case 0: parity *= dir == 0 ? -1 : +1; break;
            case 1: parity *= dir == 1 ? -1 : +1; break;
            case 2: parity *= dir == 2 ? -1 : +1; break;
            default: assert (0);
            }
            break;
          case TENSOR:
            switch (tcomponent) {
            case 0: parity *= +1; break;
            case 1: parity *= (dir == 0 || dir == 1) ? -1 : +1; break;
            case 2: parity *= (dir == 0 || dir == 2) ? -1 : +1; break;
            case 3: parity *= (dir == 1 || dir == 0) ? -1 : +1; break;
            case 4: parity *= +1; break;
            case 5: parity *= (dir == 1 || dir == 2) ? -1 : +1; break;
            case 6: parity *= (dir == 2 || dir == 0) ? -1 : +1; break;
            case 7: parity *= (dir == 2 || dir == 1) ? -1 : +1; break;
            case 8: parity *= +1; break;
            default: assert (0);
            }
            break;
          case WEYLSCALARS_REAL: {
            static int const weylparities[10][3] =
              {{+1,+1,+1},
               {-1,-1,-1},
               {+1,+1,+1},
               {-1,-1,-1},
               {+1,+1,+1},
               {-1,-1,-1},
               {+1,+1,+1},
               {-1,-1,-1},
               {+1,+1,+1},
               {-1,-1,-1}};
            parity *= weylparities[tcomponent][dir];
            break;
          }
          case MANUALCARTESIAN:
            parity = manual_parities[dir];
            break;
          default:
            assert (0);
          }
          
          /* Find region extent */
          for (d=0; d<3; ++d) {
            lsh[d] = cctkGH->cctk_lsh[d];
            imin[d] = 0;
            imax[d] = cctkGH->cctk_lsh[d];
            ioff[d] = 0;
            idir[d] = 1;
          }
          if (face == 0) {
            imax[dir] = cctkGH->cctk_nghostzones[dir];
            ioff[dir] = (+ 2*cctkGH->cctk_nghostzones[dir] - 1 
                         + (do_stagger[2*dir+face] ? 0 : 1));
            idir[dir] = -1;
          } else {
            imin[dir] = cctkGH->cctk_lsh[dir] - cctkGH->cctk_nghostzones[dir];
            ioff[dir] = (- 2*cctkGH->cctk_nghostzones[dir] + 1
                         - (do_stagger[2*dir+face] ? 0 : 1));
            idir[dir] = -1;
          }
          
          /* Ensure that there are sufficient interior zones, since
             this thorn does not support filling symmetry zones from
             other symmetry zones */
          {
            int const have_points = cctkGH->cctk_gsh[dir];
            int const need_points =
              3 * cctkGH->cctk_nghostzones[dir]
              + !do_stagger[2*dir] + !do_stagger[2*dir+1];
            if (need_points > have_points) {
              CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                          "Cannot apply symmetry boundary zones in the %s %c direction, since there seem to be more symmetry zones than interior zones",
                          (face==0 ? "lower" : "upper"),
                          "xyz"[dir]);
            }
          }
          
          /* Copy region */
          switch (group.vartype) {
            
#define ARGS  varptr, varptr,                   \
              lsh[0], lsh[1], lsh[2],           \
              imin[0], imin[1], imin[2],        \
              imax[0], imax[1], imax[2],        \
              ioff[0], ioff[1], ioff[2],        \
              idir[0], idir[1], idir[2],        \
              parity

#ifdef HAVE_CCTK_INT1
          case CCTK_VARIABLE_INT1:
#ifdef CCTK_INTEGER_PRECISION_1
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT1 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_INT2
          case CCTK_VARIABLE_INT2:
#ifdef CCTK_INTEGER_PRECISION_2
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT2 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_INT4
          case CCTK_VARIABLE_INT4:
#ifdef CCTK_INTEGER_PRECISION_4
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT4 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_INT8
          case CCTK_VARIABLE_INT8:
#ifdef CCTK_INTEGER_PRECISION_8
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT8 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_REAL4
          case CCTK_VARIABLE_REAL4:
#ifdef CCTK_REAL_PRECISION_4
          case CCTK_VARIABLE_REAL:
#endif
            copy_CCTK_REAL4 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_REAL8
          case CCTK_VARIABLE_REAL8:
#ifdef CCTK_REAL_PRECISION_8
          case CCTK_VARIABLE_REAL:
#endif
            copy_CCTK_REAL8 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_REAL16
          case CCTK_VARIABLE_REAL16:
#ifdef CCTK_REAL_PRECISION_16
          case CCTK_VARIABLE_REAL:
#endif
            copy_CCTK_REAL16 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_COMPLEX8
          case CCTK_VARIABLE_COMPLEX8:
#ifdef CCTK_COMPLEX_PRECISION_8
          case CCTK_VARIABLE_COMPLEX:
#endif
            copy_CCTK_COMPLEX8 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_COMPLEX16
          case CCTK_VARIABLE_COMPLEX16:
#ifdef CCTK_COMPLEX_PRECISION_16
          case CCTK_VARIABLE_COMPLEX:
#endif
            copy_CCTK_COMPLEX16 (ARGS);
            break;
#endif
            
#ifdef HAVE_CCTK_COMPLEX32
          case CCTK_VARIABLE_COMPLEX32:
#ifdef CCTK_COMPLEX_PRECISION_32
          case CCTK_VARIABLE_COMPLEX:
#endif
            copy_CCTK_COMPLEX32 (ARGS);
            break;
#endif
            
#undef ARGS
            
          default:
            CCTK_WARN (0, "Unsupported variable type");
          }
          
        } /* if cctk_bbox */
      } /* if do_reflection */
    } /* for face */
  } /* for dir */
  
  /* Success */
  return 0;
}



/* When CoordBase is used to specify the location of the boundary
   points, then ensure that the CoordBase parameters and this thorn's
   parameters are consistent.  */
static
void
CheckBoundaryParameters (cGH const * restrict const cctkGH,
                         int const vi,
                         int const * restrict const stencil)
{
  DECLARE_CCTK_PARAMETERS;
  
  static int did_check = 0;
  
  int type;                     /* Parameter type */
  void const * ptr;             /* Pointer to parameter value */
  char const * coordtype;       /* CartGrid3D::type */
  
  int dim;                      /* Number of dimensions of vi */
  
  CCTK_INT * restrict nboundaryzones; /* CoordBase boundary location */
  CCTK_INT * restrict is_internal;
  CCTK_INT * restrict is_staggered;
  CCTK_INT * restrict shiftout;
  
  int do_reflection[6];         /* This thorn's parameters */
  int do_stagger[6];
  
  int d;
  int ierr;
  
  
  
  /* Check only once to save time */
  if (did_check) return;
  
  /* Check only for grid functions */
  if (CCTK_GroupTypeFromVarI (vi) != CCTK_GF) return;
  
  did_check = 1;
  
  /* Check whether CartGrid3D is active */
  if (! CCTK_IsThornActive ("CartGrid3D")) return;
  
  /* Check whether CoordBase is used */
  ptr = CCTK_ParameterGet ("type", "CartGrid3D", & type);
  assert (ptr != 0);
  assert (type == PARAMETER_KEYWORD);
  coordtype = * (char const * const *) ptr;
  if (! CCTK_EQUALS (coordtype, "coordbase")) return;
  
  /* Get the boundary specification */
  dim = CCTK_GroupDimFromVarI (vi);
  assert (dim >= 0);
  nboundaryzones = malloc (2*dim * sizeof *nboundaryzones);
  is_internal = malloc (2*dim * sizeof *is_internal);
  is_staggered = malloc (2*dim * sizeof *is_staggered);
  shiftout = malloc (2*dim * sizeof *shiftout);
  ierr = GetBoundarySpecification
    (2*dim, nboundaryzones, is_internal, is_staggered, shiftout);
  assert (! ierr);
  
  /* Reflection symmetry information */
  assert (dim == 3);
  do_reflection[0] = reflection_x;
  do_reflection[1] = reflection_upper_x;
  do_reflection[2] = reflection_y;
  do_reflection[3] = reflection_upper_y;
  do_reflection[4] = reflection_z;
  do_reflection[5] = reflection_upper_z;
  
  do_stagger[0] = avoid_origin_x;
  do_stagger[1] = avoid_origin_upper_x;
  do_stagger[2] = avoid_origin_y;
  do_stagger[3] = avoid_origin_upper_y;
  do_stagger[4] = avoid_origin_z;
  do_stagger[5] = avoid_origin_upper_z;
  
  /* Check the boundary sizes */
  for (d=0; d<6; ++d) {
    if (do_reflection[d]) {
      char const dir = "xyz"[d/2];
      char const * const face = (d%2==0) ? "lower" : "upper";
      if (stencil[d/2] != nboundaryzones[d]) {
        CCTK_VWarn (CCTK_WARN_ABORT,
                    __LINE__, __FILE__, CCTK_THORNSTRING,
                    "The %s %c face is a symmetry boundary.  Since there are %d ghost zones in the %c direction, the corresponding CoordBase boundary width must also be %d.  The boundary width is currently %d.",
                    face, dir,
                    stencil[d/2], dir, stencil[d/2],
                    (int) nboundaryzones[d]);
      }
      if (is_internal[d]) {
        CCTK_VWarn (CCTK_WARN_ABORT,
                    __LINE__, __FILE__, CCTK_THORNSTRING,
                    "The %s %c face is a symmetry boundary.  The corresponding CoordBase boundary must not be internal.",
                    face, dir);
      }
      if (do_stagger[d] != is_staggered[d]) {
        CCTK_VWarn (CCTK_WARN_ABORT,
                    __LINE__, __FILE__, CCTK_THORNSTRING,
                    "The %s %c face is a symmetry boundary.  The symmetry condition and the corresponding CoordBase boundary must either be both staggered or both not staggered.",
                    face, dir);
      }
      if ((do_stagger[d] ? 0 : 1) != shiftout[d]) {
        CCTK_VWarn (CCTK_WARN_ABORT,
                    __LINE__, __FILE__, CCTK_THORNSTRING,
                    "The %s %c face is a symmetry boundary.  If the symmetry condition is staggered, then the corresponding CoordBase shiftout must be 0; otherwise it must be 1.",
                    face, dir);
      }
    }
  }
  
  /* Free memory */
  free (nboundaryzones);
  free (is_internal);
  free (is_staggered);
  free (shiftout);
}



void
ReflectionSymmetry_Apply (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  
  int nvars;
  CCTK_INT * restrict indices;
  CCTK_INT * restrict faces;
  CCTK_INT * restrict widths;
  CCTK_INT * restrict tables;
  int vi;
  int dim;
  int * restrict stencil;
  int i;
  int istat;
  int ierr;
  
  if (! cctkGH) {
    CCTK_WARN (0, "Argument cctkGH is NULL");
  }
  
  nvars = Boundary_SelectedGVs (cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
  if (nvars < 0) {
    CCTK_WARN (0, "Internal error in Boundary_SelectedGVs");
  }
  
  if (nvars == 0) {
    /* Nothing to do */
    return;
  }
  
  indices = malloc (nvars * sizeof *indices);
  if (! indices) {
    CCTK_WARN (0, "Out of memory");
  }
  faces = malloc (nvars * sizeof *faces);
  if (! faces) {
    CCTK_WARN (0, "Out of memory");
  }
  widths = malloc (nvars * sizeof *widths);
  if (! widths) {
    CCTK_WARN (0, "Out of memory");
  }
  tables = malloc (nvars * sizeof *tables);
  if (! tables) {
    CCTK_WARN (0, "Out of memory");
  }
  
  istat =  Boundary_SelectedGVs
    (cctkGH, nvars, indices, faces, widths, tables, 0);
  if (istat != nvars) {
    CCTK_WARN (0, "Internal error in Boundary_SelectedGVs");
  }
  
  for (i=0; i<nvars; ++i) {
    vi = indices[i];
    if (vi < 0 || vi >= CCTK_NumVars()) {
      CCTK_WARN (0, "Illegal variable index");
    }
    
    if (widths[i] < 0) {
      CCTK_WARN (0, "Illegal boundary width");
    }
    
    dim = CCTK_GroupDimFromVarI (vi);
    if (dim < 0) {
      CCTK_WARN (0, "Illegal dimension");
    }
    
    stencil = malloc (dim * sizeof *stencil);
    if (! stencil) {
      CCTK_WARN (0, "Out of memory");
    }
    ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi);
    if (ierr) {
      CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI");
    }
    
    CheckBoundaryParameters (cctkGH, vi, stencil);
    
    ierr = BndReflectVI (cctkGH, vi);
    if (ierr) {
      CCTK_WARN (0, "Internal error in BndReflectVI");
    }
    
    free (stencil);
  }
  
  free (indices);
  free (faces);
  free (widths);
  free (tables);
}