A thorn that applies reflection symmetry boundary conditions,

including symmetry interpolation.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@2 082bdb00-0f4f-0410-b49e-b1835e5f2039

10 files changed, 1385 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..0f56f7a
--- /dev/null
+++ b/README
@@ -0,0 +1,10 @@
+CVS info   : $Header$
+
+Cactus Code Thorn ReflectionSymmetry
+Thorn Author(s)     : Erik Schnetter <schnetter@aei.mpg.de>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Provide reflection symmetries, i.e. bitant, quadrant, and octant mode.
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..fbc1f83
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header$
+%
+%  Thorn documentation in the latex file doc/documentation.tex
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include
+%  pages about variables, parameters and scheduling parsed from the
+%  relevant thorn CCL files.
+%
+%  This template contains guidelines which help to assure that your
+%  documentation will be correctly added to ThornGuides. More
+%  information is available in the Cactus UsersGuide.
+%
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations',
+%     'references', and 'image' names should conform to the following
+%     convention:
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package.
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command.
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as
+%     standard sections.
+%   - For the benefit of our Perl scripts, and for future extensions,
+%     please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header$
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+%    \date{$ $Date$ $}
+\date{June 01 2004}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..1a335bc
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,68 @@
+# Interface definition for thorn ReflectionSymmetry
+# $Header$
+
+IMPLEMENTS: ReflectionSymmetry
+
+
+
+CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name)
+REQUIRES FUNCTION SymmetryRegister
+
+CCTK_INT FUNCTION                                           \
+    SymmetryRegisterGrid                                    \
+        (CCTK_POINTER IN cctkGH,                            \
+         CCTK_INT IN sym_handle,                            \
+         CCTK_INT IN ARRAY which_faces,                     \ # array [N_FACES]
+         CCTK_INT IN ARRAY symmetry_zone_width)               # array [N_FACES]
+REQUIRES FUNCTION SymmetryRegisterGrid
+
+CCTK_INT FUNCTION                                           \
+    SymmetryRegisterGridInterpolator                        \
+        (CCTK_POINTER IN cctkGH,                            \
+         CCTK_INT IN sym_handle,                            \
+         CCTK_INT CCTK_FPOINTER IN symmetry_interpolate     \
+             (CCTK_POINTER_TO_CONST IN cctkGH,              \
+              CCTK_INT IN N_dims,                           \
+              CCTK_INT IN local_interp_handle,              \
+              CCTK_INT IN param_table_handle,               \
+              CCTK_INT IN coord_system_handle,              \
+              CCTK_INT IN N_interp_points,                  \
+              CCTK_INT IN interp_coords_type,               \
+              CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \
+              CCTK_INT IN N_input_arrays,                   \
+              CCTK_INT ARRAY IN input_array_indices,        \
+              CCTK_INT IN N_output_arrays,                  \
+              CCTK_INT ARRAY IN output_array_types,         \
+              CCTK_POINTER ARRAY IN output_arrays,          \
+              CCTK_INT IN faces))
+REQUIRES FUNCTION SymmetryRegisterGridInterpolator
+
+CCTK_INT FUNCTION                                           \
+    SymmetryInterpolateFaces                                \
+        (CCTK_POINTER_TO_CONST IN cctkGH,                   \
+         CCTK_INT IN N_dims,                                \
+         CCTK_INT IN local_interp_handle,                   \
+         CCTK_INT IN param_table_handle,                    \
+         CCTK_INT IN coord_system_handle,                   \
+         CCTK_INT IN N_interp_points,                       \
+         CCTK_INT IN interp_coords_type,                    \
+         CCTK_POINTER_TO_CONST ARRAY IN interp_coords,      \
+         CCTK_INT IN N_input_arrays,                        \
+         CCTK_INT ARRAY IN input_array_indices,             \
+         CCTK_INT IN N_output_arrays,                       \
+         CCTK_INT ARRAY IN output_array_types,              \
+         CCTK_POINTER ARRAY IN output_arrays,               \
+         CCTK_INT IN faces)
+REQUIRES FUNCTION SymmetryInterpolateFaces
+
+
+
+CCTK_INT FUNCTION Boundary_SelectedGVs \
+  (CCTK_POINTER_TO_CONST IN  GH, \
+   CCTK_INT IN  array_size, \
+   CCTK_INT ARRAY OUT var_indicies, \
+   CCTK_INT ARRAY OUT faces, \
+   CCTK_INT ARRAY OUT boundary_widths, \
+   CCTK_INT ARRAY OUT table_handles, \
+   CCTK_STRING IN  bc_name)
+REQUIRES FUNCTION Boundary_SelectedGVs
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..f6406ec
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,62 @@
+# Parameter definitions for thorn ReflectionSymmetry
+# $Header$
+
+BOOLEAN verbose "Produce screen output while applying boundary conditions"
+{
+} "no"
+
+
+
+BOOLEAN reflection_x "Reflection symmetry at the lower x boundary"
+{
+} "no"
+
+BOOLEAN reflection_y "Reflection symmetry at the lower y boundary"
+{
+} "no"
+
+BOOLEAN reflection_z "Reflection symmetry at the lower z boundary"
+{
+} "no"
+
+
+
+BOOLEAN reflection_upper_x "Reflection symmetry at the upper x boundary"
+{
+} "no"
+
+BOOLEAN reflection_upper_y "Reflection symmetry at the upper y boundary"
+{
+} "no"
+
+BOOLEAN reflection_upper_z "Reflection symmetry at the upper z boundary"
+{
+} "no"
+
+
+
+BOOLEAN avoid_origin_x "Stagger about the origin on the lower x boundary?"
+{
+} "yes"
+
+BOOLEAN avoid_origin_y "Stagger about the origin on the lower y boundary?"
+{
+} "yes"
+
+BOOLEAN avoid_origin_z "Stagger about the origin on the lower z boundary?"
+{
+} "yes"
+
+
+
+BOOLEAN avoid_origin_upper_x "Stagger about the origin on the upper x boundary?"
+{
+} "yes"
+
+BOOLEAN avoid_origin_upper_y "Stagger about the origin on the upper y boundary?"
+{
+} "yes"
+
+BOOLEAN avoid_origin_upper_z "Stagger about the origin on the upper z boundary?"
+{
+} "yes"
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..8c5db0e
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,12 @@
+# Schedule definitions for thorn ReflectionSymmetry
+# $Header$
+
+SCHEDULE ReflectionSymmetry_Register IN SymmetryRegister
+{
+  LANG: C
+} "Register reflection symmetry boundaries"
+
+SCHEDULE ReflectionSymmetry_Apply IN BoundaryConditions
+{
+  LANG: C
+} "Apply reflection symmetries"
diff --git a/src/apply.c b/src/apply.c
new file mode 100644
index 0000000..114c1f7
--- /dev/null
+++ b/src/apply.c
@@ -0,0 +1,627 @@
+#include <assert.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"
+
+
+
+static const char rcsid[] = "$Header$";
+CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_apply_c);
+
+
+
+#define COPY(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)                                     \
+  {                                                                     \
+    int i, j, k;                                                        \
+                                                                        \
+    for (k=kmin; k<kmax; ++k)                                           \
+    {                                                                   \
+      for (j=jmin; j<jmax; ++j)                                         \
+      {                                                                 \
+        for (i=imin; j<imax; ++i)                                       \
+        {                                                               \
+          int const dstind = i + ni * (j + nj * k);                     \
+          int const ii = ioff + idir * i;                               \
+          int const jj = joff + jdir * j;                               \
+          int const kk = koff + 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 REAL(x) x
+#define IMAG(x) /* nothing */
+#define RE /* nothing */
+#define IM /* nothing */
+
+#ifdef CCTK_INT1
+COPY(CCTK_INT1)
+#endif
+
+#ifdef CCTK_INT2
+     COPY(CCTK_INT2)
+#endif
+
+#ifdef CCTK_INT4
+     COPY(CCTK_INT4)
+#endif
+
+#ifdef CCTK_INT8
+     COPY(CCTK_INT8)
+#endif
+
+#ifdef CCTK_REAL4
+     COPY(CCTK_REAL4)
+#endif
+
+#ifdef CCTK_REAL8
+     COPY(CCTK_REAL8)
+#endif
+
+#ifdef CCTK_REAL16
+     COPY(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 CCTK_REAL4
+     COPY(CCTK_COMPLEX8)
+#endif
+
+#ifdef CCTK_REAL8
+     COPY(CCTK_COMPLEX16)
+#endif
+
+#ifdef CCTK_REAL16
+     COPY(CCTK_COMPLEX32)
+#endif
+
+#undef REAL
+#undef IMAG
+#undef RE
+#undef IM
+
+#undef COPY
+
+
+
+static int
+BndReflectVI (cGH const * restrict const cctkGH,
+              int const * restrict const stencil,
+              int const vi)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  int gi;
+  cGroup group;
+  cGroupDynamicData data;
+  int firstvar, numvars;
+  char * restrict fullname;
+  char * restrict groupname;
+  
+  void * restrict varptr;
+  
+  int table;
+  char tensortypealias[1000];
+  enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR };
+  enum tensortype ttype;
+  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 d;
+  
+  int ierr;
+  
+  
+  
+  /* Check arguments */
+  if (! cctkGH)
+  {
+    CCTK_WARN (0, "Argument cctkGH is NULL");
+  }
+  if (! stencil)
+  {
+    CCTK_WARN (0, "Argument stencil 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) {
+    groupname = CCTK_GroupName(gi);
+    assert (groupname);
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Tensor type alias not declared for group \"%s\" -- assuming a scalar",
+                groupname);
+    free (groupname);
+    strcpy (tensortypealias, "scalar");
+  } else if (ierr<0) {
+    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 */
+    if (numvars != 1) {
+      groupname = CCTK_GroupName(gi);
+      assert (groupname);
+      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element",
+                  groupname);
+      free (groupname);
+    }
+    ttype = SCALAR;
+    tcomponent = 0;
+  }
+  else if (CCTK_EQUALS (tensortypealias, "u"))
+  {
+    /* vector */
+    assert (numvars == 3);
+    ttype = VECTOR;
+    tcomponent = vi - firstvar;
+  } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
+    /* symmetric tensor */
+    assert (numvars == 6);
+    ttype = TENSOR;
+    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 TENSOR:
+    assert (tcomponent>=0 && tcomponent<6);
+    break;
+  default:
+    assert (0);
+  }
+  
+  
+  
+  /* 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 */
+          switch (ttype)
+          {
+          case SCALAR:
+            parity = +1;
+            break;
+          case VECTOR:
+            parity = dir == tcomponent ? -1 : +1;
+            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 = +1; break;
+            case 4: parity = (dir == 1 || dir == 2) ? -1 : +1; break;
+            case 5: parity = +1; break;
+            default: assert (0);
+            }
+            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;
+          }
+          if (face == 0)
+          {
+            imax[d] = cctkGH->cctk_nghostzones[d];
+            ioff[d] = (+ cctkGH->cctk_nghostzones[d] 
+                       + (do_stagger[2*d+face] ? 0 : 1));
+          }
+          else
+          {
+            imin[d] = cctkGH->cctk_lsh[d] - cctkGH->cctk_nghostzones[d];
+            ioff[d] = (- cctkGH->cctk_nghostzones[d]
+                       - (do_stagger[2*d+face] ? 0 : 1));
+          }
+          
+          /* Copy region */
+          switch (group.vartype)
+          {
+            
+#ifdef CCTK_INT1
+          case CCTK_VARIABLE_INT1:
+#ifdef CCTK_INTER_PRECISION_1
+          case CCTK_VARIABLE_INT:
+#endif
+            copy_CCTK_INT1 (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);
+            break;
+#endif
+            
+#ifdef CCTK_INT2
+          case CCTK_VARIABLE_INT2:
+#ifdef CCTK_INTER_PRECISION_2
+          case CCTK_VARIABLE_INT:
+#endif
+            copy_CCTK_INT2 (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);
+            break;
+#endif
+            
+#ifdef CCTK_INT4
+          case CCTK_VARIABLE_INT4:
+#ifdef CCTK_INTER_PRECISION_4
+          case CCTK_VARIABLE_INT:
+#endif
+            copy_CCTK_INT4 (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);
+            break;
+#endif
+            
+#ifdef CCTK_INT8
+          case CCTK_VARIABLE_INT8:
+#ifdef CCTK_INTER_PRECISION_8
+          case CCTK_VARIABLE_INT:
+#endif
+            copy_CCTK_INT8 (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);
+            break;
+#endif
+            
+#ifdef CCTK_REAL4
+          case CCTK_VARIABLE_REAL4:
+#ifdef CCTK_REAL_PRECISION_4
+          case CCTK_VARIABLE_REAL:
+#endif
+            copy_CCTK_REAL4 (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);
+            break;
+#endif
+            
+#ifdef CCTK_REAL8
+          case CCTK_VARIABLE_REAL8:
+#ifdef CCTK_REAL_PRECISION_8
+          case CCTK_VARIABLE_REAL:
+#endif
+            copy_CCTK_REAL8 (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);
+            break;
+#endif
+            
+#ifdef CCTK_REAL16
+          case CCTK_VARIABLE_REAL16:
+#ifdef CCTK_REAL_PRECISION_16
+          case CCTK_VARIABLE_REAL:
+#endif
+            copy_CCTK_REAL16 (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);
+            break;
+#endif
+            
+#ifdef CCTK_REAL4
+          case CCTK_VARIABLE_COMPLEX8:
+#ifdef CCTK_COMPLEX_PRECISION_8
+          case CCTK_VARIABLE_COMPLEX:
+#endif
+            copy_CCTK_COMPLEX8 (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);
+            break;
+#endif
+            
+#ifdef CCTK_REAL8
+          case CCTK_VARIABLE_COMPLEX16:
+#ifdef CCTK_COMPLEX_PRECISION_16
+          case CCTK_VARIABLE_COMPLEX:
+#endif
+            copy_CCTK_COMPLEX16 (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);
+            break;
+#endif
+            
+#ifdef CCTK_REAL16
+          case CCTK_VARIABLE_COMPLEX32:
+#ifdef CCTK_COMPLEX_PRECISION_32
+          case CCTK_VARIABLE_COMPLEX:
+#endif
+            copy_CCTK_COMPLEX32 (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);
+            break;
+#endif
+            
+          default:
+            CCTK_WARN (0, "Unsupported variable type");
+          }
+          
+        } /* if cctk_bbox */
+      } /* if do_reflection */
+    } /* for face */
+  } /* for dir */
+  
+  /* Success */
+  return 0;
+}
+
+
+
+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");
+  }
+  
+  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");
+    }
+    
+    ierr = BndReflectVI (cctkGH, stencil, vi);
+    if (ierr)
+    {
+      CCTK_WARN (0, "Internal error in BndReflectVI");
+    }
+    
+    free (stencil);
+  }
+  
+  free (indices);
+  free (faces);
+  free (widths);
+  free (tables);
+}
diff --git a/src/interpolate.c b/src/interpolate.c
new file mode 100644
index 0000000..8591fd7
--- /dev/null
+++ b/src/interpolate.c
@@ -0,0 +1,359 @@
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "reflection.h"
+
+
+
+static const char rcsid[] = "$Header$";
+CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_interpolate_c);
+
+
+
+CCTK_INT
+ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
+                                CCTK_INT const N_dims,
+                                CCTK_INT const local_interp_handle,
+                                CCTK_INT const param_table_handle,
+                                CCTK_INT const coord_system_handle,
+                                CCTK_INT const N_interp_points,
+                                CCTK_INT const interp_coords_type,
+                                CCTK_POINTER_TO_CONST restrict const interp_coords[],
+                                CCTK_INT const N_input_arrays,
+                                CCTK_INT const input_array_indices[],
+                                CCTK_INT const N_output_arrays,
+                                CCTK_INT const output_array_types[],
+                                CCTK_POINTER restrict const output_arrays[],
+                                CCTK_INT const faces)
+{
+  cGH const * restrict const cctkGH = cctkGH_;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int do_reflection[6];
+  
+  CCTK_POINTER_TO_CONST restrict new_interp_coords[3];
+  CCTK_INT newfaces;
+  
+  CCTK_INT * restrict operand_indices;
+  CCTK_INT * restrict operation_codes;
+  CCTK_INT * restrict output_array_indices;
+  
+  int dir;
+  
+  int iret;
+  
+  int m;
+  int n;
+  int d;
+  
+  int ierr;
+  
+  
+  
+  /* Get 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;
+  
+  newfaces = faces;
+  for (d=0; d<6; ++d)
+  {
+    if (do_reflection[d])
+    {
+      assert (newfaces & (1 << d));
+      newfaces &= ~ (1 << d);
+    }
+  }
+  
+  
+  
+  /* Fold coordinates */
+  assert (interp_coords_type == CCTK_VARIABLE_REAL);
+  for (dir=0; dir<3; ++dir)
+  {
+    assert (! do_reflection[2*dir+1]);
+    
+    if (do_reflection[2*dir])
+    {
+      new_interp_coords[dir]
+        = malloc (N_interp_points * sizeof *new_interp_coords[dir]);
+      assert (new_interp_coords[d]);
+      
+      for (n=0; n<N_interp_points; ++n)
+      {
+        CCTK_REAL const pos = ((CCTK_REAL const *)interp_coords[dir])[n];
+        CCTK_REAL const newpos = fabs(pos);
+        ((CCTK_REAL *)new_interp_coords[dir])[n] = newpos;
+      }
+    }
+    else
+    {
+      new_interp_coords[dir] = interp_coords[dir];
+    }
+  }
+  
+  
+  
+  /* Recursive call */
+  iret = SymmetryInterpolateFaces
+    (cctkGH_, 
+     N_dims, local_interp_handle, param_table_handle, coord_system_handle,
+     N_interp_points, interp_coords_type, new_interp_coords,
+     N_input_arrays, input_array_indices,
+     N_output_arrays, output_array_types, output_arrays,
+     newfaces);
+  
+  
+  
+  /* Free coordinates */
+  for (dir=0; dir<3; ++dir)
+  {
+    if (do_reflection[2*dir])
+    {
+      free (new_interp_coords[dir]);
+    }
+  }
+  
+  
+  /* Find output variable indices */
+  operand_indices = malloc (N_output_arrays * sizeof *operand_indices);
+  assert (operand_indices);
+  ierr = Util_TableGetIntArray
+    (param_table_handle, N_output_arrays, operand_indices, "operand_indices");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    assert (N_output_arrays == N_input_arrays);
+    for (m=0; m<N_output_arrays; ++m) {
+      operand_indices[m] = m;   /* set output index to input index */
+    }
+  } else {
+    assert (ierr == N_output_arrays);
+  }
+  
+  operation_codes = malloc (N_output_arrays * sizeof *operation_codes);
+  assert (operation_codes);
+  ierr = Util_TableGetIntArray
+    (param_table_handle, N_output_arrays, operation_codes, "operation_codes");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    assert (N_output_arrays == N_input_arrays);
+    for (m=0; m<N_output_arrays; ++m) {
+      operation_codes[m] = 0;     /* do not take derivatives */
+    }
+  } else {
+    assert (ierr == N_output_arrays);
+  }
+  
+  output_array_indices
+    = malloc (N_output_arrays * sizeof *output_array_indices);
+  assert (output_array_indices);
+  for (m=0; m<N_output_arrays; ++m) {
+    assert (operand_indices[m]>=0 && operand_indices[m]<N_input_arrays);
+    output_array_indices[m] = input_array_indices[operand_indices[m]];
+     assert (output_array_indices[m]>=0
+             && output_array_indices[m]<CCTK_NumVars());
+  }
+  
+  
+  
+  /* Unfold tensor types */
+  for (m=0; m<N_output_arrays; ++m)
+  {
+    
+    int vi, gi;
+    int numvars, firstvar;
+    char * groupname;
+    
+    int table;
+    char tensortypealias[1000];
+    enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR };
+    enum tensortype ttype;
+    int tcomponent;
+    
+    int parities[3];
+    int check_dir[3];
+    int needs_checking;
+    
+    
+    
+    vi = output_array_indices[m];
+    assert (vi>=0 && vi<CCTK_NumVars());
+    gi = CCTK_GroupIndexFromVarI (vi);
+    assert (gi>=0 && gi<CCTK_NumGroups());
+    numvars = CCTK_NumVarsInGroupI(gi);
+    assert (numvars>0);
+    firstvar = CCTK_FirstVarIndexI(gi);
+    assert (firstvar>=0);
+    
+    
+    
+    /* Get and check tensor type information */
+    ierr = Util_TableGetString
+      (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
+    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+      groupname = CCTK_GroupName(gi);
+      assert (groupname);
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Tensor type alias not declared for group \"%s\" -- assuming a scalar",
+                  groupname);
+      free (groupname);
+      strcpy (tensortypealias, "scalar");
+    } else if (ierr<0) {
+      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 */
+      if (numvars != 1) {
+        groupname = CCTK_GroupName(gi);
+        assert (groupname);
+        CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element",
+                    groupname);
+        free (groupname);
+      }
+      ttype = SCALAR;
+      tcomponent = 0;
+    }
+    else if (CCTK_EQUALS (tensortypealias, "u"))
+    {
+      /* vector */
+      assert (numvars == 3);
+      ttype = VECTOR;
+      tcomponent = vi - firstvar;
+    } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
+      /* symmetric tensor */
+      assert (numvars == 6);
+      ttype = TENSOR;
+      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 TENSOR:
+      assert (tcomponent>=0 && tcomponent<6);
+      break;
+    default:
+      assert (0);
+    }
+    
+    
+    
+    /* Calculate parities */
+    parities[0] = parities[1] = parities[2] = +1;
+    switch (ttype)
+    {
+    case SCALAR:
+      /* do nothing */
+      break;
+    case VECTOR:
+      parities[tcomponent] = -1;
+      break;
+    case TENSOR:
+      switch (tcomponent)
+      {
+      case 0: break;
+      case 1: parities[0] = parities[1] = -1; break;
+      case 2: parities[0] = parities[2] = -1; break;
+      case 3: break;
+      case 4: parities[1] = parities[2] = -1; break;
+      case 5: break;
+      default: assert (0);
+      }
+      break;
+    default:
+      assert (0);
+    }
+    
+    
+    
+    /* Take derivatives into account */
+    {
+      int code = operation_codes[m];
+      while (code) {
+        const int d = code % 10 - 1;
+        code /= 10;
+        assert (d>=0 && d<3);
+        parities[d] *= -1;
+      }
+    }
+    
+    
+    
+    /* Are there negative parities? */
+    needs_checking = 0;
+    for (dir=0; dir<3; ++dir)
+    {
+      check_dir[d] = do_reflection[2*dir] && parities[dir] < 0;
+      needs_checking |= check_dir[d];
+    }
+    
+    
+    
+    /* Loop over all points and unfold */
+    if (needs_checking)
+    {
+      for (n=0; n<N_interp_points; ++n)
+      {
+        int parity = +1;
+        /* Is the point outside the domain? */
+        for (dir=0; dir<3; ++dir)
+        {
+          if (check_dir[dir])
+          {
+            CCTK_REAL const pos = ((CCTK_REAL const *)interp_coords[dir])[n];
+            if (pos < 0)
+            {
+              /* Reflect the tensor component */
+              parity *= -1;
+            }
+          }
+        }
+        ((CCTK_REAL *)output_arrays[m])[n] *= parity;
+      }
+    }
+    
+  } /* for m */
+  
+  
+  
+  /* Free output variable indices */
+  free (operand_indices);
+  free (operation_codes);
+  free (output_array_indices);
+  
+  
+  
+  /* Return */
+  return iret;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..24a04f5
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn ReflectionSymmetry
+# $Header$
+
+# Source files in this directory
+SRCS = apply.c interpolate.c register.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/reflection.h b/src/reflection.h
new file mode 100644
index 0000000..fea4729
--- /dev/null
+++ b/src/reflection.h
@@ -0,0 +1,29 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+#ifndef REFLECTION_H
+#define REFLECTION_H
+
+CCTK_INT
+ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
+                                CCTK_INT const N_dims,
+                                CCTK_INT const local_interp_handle,
+                                CCTK_INT const param_table_handle,
+                                CCTK_INT const coord_system_handle,
+                                CCTK_INT const N_interp_points,
+                                CCTK_INT const interp_coords_type,
+                                CCTK_POINTER_TO_CONST restrict const interp_coords[],
+                                CCTK_INT const N_input_arrays,
+                                CCTK_INT const input_array_indices[],
+                                CCTK_INT const N_output_arrays,
+                                CCTK_INT const output_array_types[],
+                                CCTK_POINTER restrict const output_arrays[],
+                                CCTK_INT const faces);
+
+void
+ReflectionSymmetry_Register (CCTK_ARGUMENTS);
+
+void
+ReflectionSymmetry_Apply (CCTK_ARGUMENTS);
+
+#endif /* REFLECTION_H */
diff --git a/src/register.c b/src/register.c
new file mode 100644
index 0000000..2ae4b58
--- /dev/null
+++ b/src/register.c
@@ -0,0 +1,66 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "reflection.h"
+
+
+
+static const char rcsid[] = "$Header$";
+CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_register_c);
+
+
+
+void
+ReflectionSymmetry_Register (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int do_reflection[6];
+  CCTK_INT handle;
+  CCTK_INT faces[6];
+  CCTK_INT width[6];
+  int f;
+  CCTK_INT ierr;
+  
+  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;
+  
+  for (f=0; f<6; ++f)
+  {
+    if (do_reflection[f])
+    {
+      faces[f] = 1;
+      width[f] = cctk_nghostzones[f/2];
+    }
+    else
+    {
+      faces[f] = 0;
+      width[f] = 0;
+    }
+  }
+  
+  handle = SymmetryRegister ("reflection_symmetry");
+  if (handle < 0)
+  {
+    CCTK_WARN (0, "Could not register symmetry boundary condition");
+  }
+  
+  ierr = SymmetryRegisterGrid (cctkGH, handle, faces, width);
+  if (ierr < 0)
+  {
+    CCTK_WARN (0, "Could not register the symmetry boundaries -- probably some other thorn has already registered the same boundary faces for a different symmetry");
+  }
+  
+  ierr = SymmetryRegisterGridInterpolator
+    (cctkGH, handle, ReflectionSymmetry_Interpolate);
+  if (ierr < 0)
+  {
+    CCTK_WARN (0, "Could not register the symmetry interpolator");
+  }
+}