Provide infrastructure to take symmetries into account for global

interpolations.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/SymBase/trunk@11 906471b6-c639-44d1-9ea0-3e3d6879f074

10 files changed, 1068 insertions, 12 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index c37694e..8e1808b 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -45,11 +45,11 @@
 %
 % 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}
+%       \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:
@@ -296,6 +296,179 @@ end do
 
 
 
+\section{Symmetry Interpolation}
+
+When the computational domain has symmetries, then it is often very
+convenient to be able to interpolate at points that are not present on
+the actual computational grid, but can be mapped into the grid through
+the symmetries.  Thorn SymBase provides a mechanism by which symmetry
+conditions can register routines that handle this mapping when a
+global interpolator is called.  Additionally, the driver has to be
+aware that is calls thorn SymBase's mapping routine before it actually
+interpolates.  The whole mechanism is transparent for the user.
+
+The mechanism by which the grid points are mapped into the domain
+works as follows:
+\begin{enumerate}
+\item The user calls \texttt{CCTK\_Interp} with a list of coordinates.
+\item The flesh forwards this call to the driver.
+\item The driver calls SymBase's aliased function.
+  \texttt{SymmetryInterpolate}, passing along all arguments.
+\item SymBase sets a flag for all boundaries that have a symmetry
+  condition associated with it, and then calls
+  \texttt{SymmetryInterpolateFaces}, passing along all arguments.
+  This is the beginning of a chain of recursive calls.
+\item \texttt{SymmetryInterpolateFaces} checks whether any faces are
+  flagged.
+\item If no faces are flagged, SymBase calls the driver's aliased
+  function \texttt{DriverInterpolate}, which performs the actual
+  interpolation.  This ends the chain of recursive calls.
+\item If there are faces with symmetry conditions flagged, SymBase
+  chooses one such face, and then calls this face's symmetry
+  condition's ``symmetry interpolation'' routine, passing along all
+  arguments.
+\item The ``symmetry interpolation'' routine maps the coordinates into
+  the domain by applying the symmetry condition for this face.  It
+  then removes the flag for the corresponding face, and calls
+  \texttt{SymmetryInterpolateFaces}, passing along the arguments with
+  the changed interpolation locations.
+\item After the actual interpolation has happened in the driver, the
+  recursive call will return.  The ``symmetry interpolation'' routine
+  then examines the tensor types of the interpolated quantities and
+  un-maps the values back onto their original locations.  That is,
+  e.g., after a reflection on the lower $x$-boundary, $x$-components
+  of vectors need their sign changed.
+\item The chain of recursive calls unravels until the call to
+  \texttt{CCTK\_Interp} returns.
+\end{enumerate}
+
+This mechanism has thus four players:
+\begin{itemize}
+\item The driver forwards the call to SymBase so that the list of
+  interpolation points can be mapped into the domain.
+\item Thorn SymBase controls which symmetry conditions perform this
+  mapping on which faces.
+\item Each symmetry boundary condition has to register a ``symmetry
+  interpolation'' routine that first maps the points into the domain,
+  then calls SymBase recursively, and finally corrects the tensor
+  types of the interpolated quantities.
+\item Finally, the user calls \texttt{CCTK\_Interp} as before, and
+  everything happens transparently for him.
+\end{itemize}
+
+
+
+\subsection{Driver Interaction}
+
+The driver has to call SymBase's aliased function
+\texttt{SymmetryInterpolate}, and has to provide an aliased function
+\texttt{DriverInterpolate}.  Both functions have essentially the same
+prototype as \texttt{CCTK\_Interp}, which are:
+
+\begin{verbatim}
+CCTK_INT FUNCTION
+    SymmetryInterpolate
+        (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)
+\end{verbatim}
+
+\begin{verbatim}
+CCTK_INT FUNCTION
+    DriverInterpolate
+        (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)
+\end{verbatim}
+
+
+
+\subsection{Interaction With Symmetry Conditions}
+
+The symmetry conditions have to register their ``symmetry
+interpolation'' routines by calling SymBase's aliased function
+\texttt{SymmetryRegisterGridInterpolator}.  The ``symmetry
+interpolation'' routine must use C linkage and must have the prototype
+
+\begin{verbatim}
+CCTK_INT 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)
+\end{verbatim}
+
+These arguments are essentially the same as those for
+\texttt{CCTK\_Interp}, except that the bit field \texttt{faces} flags
+those faces that still their symmetry boundary condition applied to
+the interpolation points.  The aliased function
+\texttt{SymmetryRegisterGridInterpolator} has the prototype
+
+\begin{verbatim}
+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))
+\end{verbatim}
+
+which takes a function pointer to the aforementioned ``symmetry
+interpolation'' routine, while \texttt{sym\_handle} specifies which
+symmetry condition this routine is for.  This handle must have been
+obtained from \texttt{SymmetryRegister}.
+
+\begin{quote}
+  The routine \texttt{SymmetryRegisterGridInterpolator} must be called
+  \emph{after} the symmetry faces have been selected by the call to
+  \texttt{SymmetryRegisterGrid}.
+\end{quote}
+
+
+
 % Do not delete next line
 % END CACTUS THORNGUIDE
 
diff --git a/interface.ccl b/interface.ccl
index 2c07bd8..c8ddacc 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -28,10 +28,10 @@ PROVIDES FUNCTION SymmetryNameOfHandle \
 
 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]
+        (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]
 PROVIDES FUNCTION SymmetryRegisterGrid \
      WITH SymBase_SymmetryRegisterGrid \
      LANGUAGE C
@@ -60,6 +60,33 @@ PROVIDES FUNCTION SymmetryRegisterGN \
 
 
 
+# Register a symmetry interpolator:
+
+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))
+PROVIDES FUNCTION SymmetryRegisterGridInterpolator \
+     WITH SymBase_SymmetryRegisterGridInterpolator \
+     LANGUAGE C
+
+
+
 # Get the symmetry table handle for a grid or grid array:
 
 CCTK_INT FUNCTION \
@@ -83,3 +110,63 @@ CCTK_INT FUNCTION                          \
 PROVIDES FUNCTION SymmetryTableHandleForGN \
      WITH SymBase_SymmetryTableHandleForGN \
      LANGUAGE C
+
+
+
+# Interpolation
+
+CCTK_INT FUNCTION                                           \
+    SymmetryInterpolate                                     \
+        (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)
+PROVIDES FUNCTION SymmetryInterpolate \
+    WITH SymBase_SymmetryInterpolate  \
+    LANGUAGE C
+
+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)
+PROVIDES FUNCTION SymmetryInterpolateFaces \
+    WITH SymBase_SymmetryInterpolateFaces  \
+    LANGUAGE C
+
+CCTK_INT FUNCTION                                           \
+    DriverInterpolate                                       \
+        (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)
+USES FUNCTION DriverInterpolate
diff --git a/param.ccl b/param.ccl
index e82ab1e..7803536 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,2 +1,6 @@
 # Parameter definitions for thorn SymBase
 # $Header$
+
+BOOLEAN verbose "Output symmetry boundary face descriptions after registration"
+{
+} "yes"
diff --git a/schedule.ccl b/schedule.ccl
index f3d7f54..350d1f9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -6,6 +6,15 @@ SCHEDULE SymBase_Startup AT CCTK_STARTUP
   LANG: C
 } "Register GH Extension for SymBase"
 
-SCHEDULE GROUP SymmetryRegister AT CCTK_WRAGH
+SCHEDULE GROUP SymBase_Wrapper AT CCTK_WRAGH
+{
+} "Wrapper group for SymBase"
+
+SCHEDULE GROUP SymmetryRegister IN SymBase_Wrapper
 {
 } "Register your symmetries here"
+
+SCHEDULE SymBase_Statistics IN SymBase_Wrapper AFTER SymmetryRegister
+{
+  LANG: C
+} "Print symmetry boundary face descriptions"
diff --git a/src/Faces.c b/src/Faces.c
index c1a668f..844dfc2 100644
--- a/src/Faces.c
+++ b/src/Faces.c
@@ -360,3 +360,168 @@ SymBase_SymmetryRegisterGN (CCTK_POINTER const cctkGH_,
   return SymBase_SymmetryRegisterGI
     (cctkGH_, sym_handle, which_faces, new_symmetry_zone_width, group_index);
 }
+
+
+
+/*@@
+  @routine    SymBase_SymmetryRegisterInterpolatorFaces
+  @author     Erik Schnetter
+  @date       2004-05-25
+  @desc
+              Register a symmetry interpolator for certain faces
+  @enddesc
+  @var        sym_table
+  @vtype      CCTK_INT
+  @vdesc      Table which describes the grid or grid array
+  @vio        in
+  @endvar
+  @var        group_dim
+  @vtype      CCTK_INT
+  @vdesc      Dimension of the grid or grid array
+  @vio        in
+  @endvar
+  @var        sym_handle
+  @vtype      CCTK_INT
+  @vdesc      Symmetry handle
+  @vio        in
+  @endvar
+  @var        symmetry_interpolate
+  @vtype      CCTK_FPOINTER
+  @vdesc      Routine that applies symmetries to the interpolation points
+              and then calls back
+              (may be NULL)
+  @vio        in
+  @endvar
+  @returntype CCTK_INT
+  @returndesc
+              status code
+              0 for success
+              -1 if sym_table has an illegal value
+              -9 if group_dim has an illegal value
+              -2 if sym_handle has an illegal value
+              -10 symmetry_interpolate is NULL
+  @endreturndesc
+@@*/
+
+CCTK_INT
+SymBase_SymmetryRegisterInterpolatorFaces (CCTK_INT const sym_table,
+                                           CCTK_INT const group_dim,
+                                           CCTK_INT const sym_handle,
+                                           SymmetryInterpolateFPointer const new_symmetry_interpolate)
+{
+  CCTK_INT symmetry_handle[100];
+  CCTK_FPOINTER symmetry_interpolate[100];
+  int face;
+  int ierr;
+
+  /* Check arguments */
+  if (sym_table < 0)
+  {
+    return -1;                  /* illegal argument */
+  }
+  if (group_dim < 0)
+  {
+    return -9;                  /* illegal argument */
+  }
+  if (sym_handle < 0 || sym_handle >= SymBase_num_symmetries)
+  {
+    return -2;                  /* illegal argument */
+  }
+  if (! new_symmetry_interpolate)
+  {
+    return -10;
+  }
+  
+  /* Get table entries */
+  if (2 * group_dim > 100)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  ierr = Util_TableGetIntArray
+    (sym_table, 2 * group_dim, symmetry_handle, "symmetry_handle");
+  if (ierr != 2 * group_dim)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  ierr = Util_TableGetFPointerArray
+    (sym_table, 2 * group_dim, symmetry_interpolate, "symmetry_interpolate");  if (ierr != 2 * group_dim)  {
+    CCTK_WARN (0, "internal error");
+  }
+  
+  /* Update table entries */
+  for (face = 0; face < 2 * group_dim; ++face)
+  {
+    if (symmetry_handle[face] == sym_handle)
+    {
+      symmetry_interpolate[face] = (CCTK_FPOINTER)new_symmetry_interpolate;
+    }
+  }
+
+  /* Set table entries */
+  ierr = Util_TableSetFPointerArray
+    (sym_table, 2 * group_dim, symmetry_interpolate, "symmetry_interpolate");
+  if (ierr != 1)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+
+  return 0;
+}
+
+
+
+/*@@
+  @routine    SymBase_SymmetryRegisterGridInterpolator
+  @author     Erik Schnetter
+  @date       2004-05-25
+  @desc
+              Register a symmetry interpolator
+              for certain faces of the grid hierarchy
+  @enddesc
+  @var        cctkGH
+  @vtype      cGH *
+  @vdesc      The grid hierarchy
+  @vio        in
+  @endvar
+  @var        sym_handle
+  @vtype      CCTK_INT
+  @vdesc      Symmetry handle
+  @vio        in
+  @endvar
+  @var        symmetry_interpolate
+  @vtype      CCTK_FPOINTER
+  @vdesc      Routine that applies symmetries to the interpolation points
+              and then calls back
+  @vio        in
+  @endvar
+  @returntype CCTK_INT
+  @returndesc
+              status code
+              0 for success
+              Error codes of SymBase_SymmetryRegisterInterpolatorFaces
+  @endreturndesc
+@@*/
+
+CCTK_INT
+SymBase_SymmetryRegisterGridInterpolator (CCTK_POINTER const cctkGH_,
+                                          CCTK_INT const sym_handle,
+                                          SymmetryInterpolateFPointer const new_symmetry_interpolate)
+{
+  cGH const *const cctkGH = cctkGH_;
+  struct SymBase const *symdata;
+
+  if (!cctkGH)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+
+  symdata = CCTK_GHExtension (cctkGH, "SymBase");
+  if (!symdata)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  
+  return SymBase_SymmetryRegisterInterpolatorFaces
+    (symdata->sym_table, cctkGH->cctk_dim, sym_handle,
+     new_symmetry_interpolate);
+}
diff --git a/src/Interpolation.c b/src/Interpolation.c
new file mode 100644
index 0000000..6873ebe
--- /dev/null
+++ b/src/Interpolation.c
@@ -0,0 +1,382 @@
+/*@@
+  @file      Interpolation.c
+  @author    Erik Schnetter
+  @date      2004-05-11
+  @desc
+             Apply symmetry conditions during interpolation
+  @version   $Header$
+  @enddesc
+@@*/
+
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "SymBase.h"
+
+
+
+/* the rcs ID and its dummy function to use it */
+static const char *const rcsid = "$Header$";
+CCTK_FILEVERSION (CactusBase_SymBase_Interpolation_c);
+
+
+
+/*@@
+  @routine    SymBase_SymmetryInterpolate
+  @author     Erik Schnetter
+  @date       2004-05-11
+  @desc
+              Adjust the coordinates of the interpolation points,
+              call the driver's interpolation function,
+              and then adjust the interpolated tensor components.
+              All arguments are equivalent to the CCTK_InterpGridArrays
+              function call.
+  @enddesc
+  @var        cctkGH
+  @vtype      CCTK_POINTER_TO_CONST
+  @vdesc      Grid hierarchy
+  @vio        in
+  @endvar
+  @var        N_dims
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        local_interp_handle
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        param_table_handle
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        coord_system_handle
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        N_interp_points
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        interp_coords_type
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        interp_coords
+  @vtype      CCTK_POINTER_TO_CONST[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        N_input_arrays
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        input_array_indices
+  @vtype      CCTK_INT[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        N_output_arrays
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        output_array_types
+  @vtype      CCTK_INT[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        output_arrays
+  @vtype      CCTK_POINTER[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @returntype CCTK_INT
+  @returndesc
+              status code
+              0 for success
+              -91 GH is NULL
+              -92 N_dims is not equal to cctkGH->cctk_dim
+  @endreturndesc
+@@*/
+
+CCTK_INT
+SymBase_SymmetryInterpolate (CCTK_POINTER_TO_CONST 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 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 const output_arrays[])
+{
+  cGH const * restrict const cctkGH = cctkGH_;
+  CCTK_INT sym_table;
+  CCTK_FPOINTER symmetry_interpolate[100];
+  CCTK_INT faces;
+  int f;
+  int ierr;
+  
+  /* Check arguments */
+  if (! cctkGH)
+  {
+    CCTK_WARN (1, "cctkGH is NULL");
+    return -91;
+  }
+  
+  if (N_dims != cctkGH->cctk_dim)
+  {
+    CCTK_WARN (1, "The number of dimensions is not equal to the GH's number of dimensions");
+    return -92;
+  }
+  
+  /* Get table */
+  sym_table = SymmetryTableHandleForGrid (cctkGH);
+  if (sym_table < 0) CCTK_WARN (0, "internal error");
+
+  /* Get table entries */
+  if (2 * cctkGH->cctk_dim > 100)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  ierr = Util_TableGetFPointerArray
+    (sym_table, 2 * cctkGH->cctk_dim,
+     symmetry_interpolate, "symmetry_interpolate");
+  if (ierr != 2 * cctkGH->cctk_dim)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  
+  /* Find all faces that need symmetries applied */
+  faces = 0;
+  for (f=0; f<2*cctkGH->cctk_dim; ++f) {
+    if (symmetry_interpolate[f]) {
+      faces |= (1 << f);
+    }
+  }
+  
+  /* Forward the call */
+  return SymBase_SymmetryInterpolateFaces
+    (cctkGH, N_dims,
+     local_interp_handle, param_table_handle, coord_system_handle,
+     N_interp_points, interp_coords_type, interp_coords,
+     N_input_arrays, input_array_indices,
+     N_output_arrays, output_array_types, output_arrays,
+     faces);
+}
+
+
+
+/*@@
+  @routine    SymBase_SymmetryInterpolateFaces
+  @author     Erik Schnetter
+  @date       2004-05-11
+  @desc
+              Adjust the coordinates of the interpolation points,
+              call the driver's interpolation function,
+              and then adjust the interpolated tensor components,
+              but only for the selected faces.
+              When no faces are selected, forward the call to the
+              driver's interpolator throug the aliased function
+              "DriverInterpolate".
+              All arguments except "faces" are equivalent to the
+              CCTK_InterpGridArrays function call.
+  @enddesc
+  @var        cctkGH
+  @vtype      CCTK_POINTER_TO_CONST
+  @vdesc      Grid hierarchy
+  @vio        in
+  @endvar
+  @var        N_dims
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        local_interp_handle
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        param_table_handle
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        coord_system_handle
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        N_interp_points
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        interp_coords_type
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        interp_coords
+  @vtype      CCTK_POINTER_TO_CONST[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        N_input_arrays
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        input_array_indices
+  @vtype      CCTK_INT[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        N_output_arrays
+  @vtype      CCTK_INT
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        output_array_types
+  @vtype      CCTK_INT[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        output_arrays
+  @vtype      CCTK_POINTER[]
+  @vdesc      same as for the CCTK_InterpGridArrays call
+  @vio        in
+  @endvar
+  @var        faces
+  @vtype      CCTK_INT
+  @vdesc      Bit mask selecting which faces' symmetries should be treated.
+  @vio        in
+  @endvar
+  @returntype CCTK_INT
+  @returndesc
+              status code
+              0 for success
+              -91 GH is NULL
+              -92 N_dims is not equal to cctkGH->cctk_dim
+              -93 faces is NULL
+              or the return value from DriverInterpolate
+  @endreturndesc
+@@*/
+
+CCTK_INT
+SymBase_SymmetryInterpolateFaces (CCTK_POINTER_TO_CONST 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 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 const output_arrays[],
+                                  CCTK_INT const faces)
+{
+  cGH const * restrict const cctkGH = cctkGH_;
+  CCTK_INT sym_table;
+  CCTK_FPOINTER symmetry_interpolate[100];
+  int face;
+  int ierr;
+  
+  /* Check arguments */
+  
+  if (! cctkGH)
+  {
+    CCTK_WARN (1, "cctkGH is NULL");
+    return -91;
+  }
+  
+  if (N_dims != cctkGH->cctk_dim)
+  {
+    CCTK_WARN (1, "The number of dimensions is not equal to the GH's number of dimensions");
+    return -92;
+  }
+  
+  /* Get table */
+  sym_table = SymmetryTableHandleForGrid (cctkGH);
+  if (sym_table < 0) CCTK_WARN (0, "internal error");
+
+  /* Get table entries */
+  if (2 * cctkGH->cctk_dim > 100)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  ierr = Util_TableGetFPointerArray
+    (sym_table, 2 * cctkGH->cctk_dim,
+     symmetry_interpolate, "symmetry_interpolate");
+  if (ierr != 2 * cctkGH->cctk_dim)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  
+  /* Find a face that needs a symmetry applied */
+  for (face=0; face<2*N_dims; ++face)
+  {
+    if (faces & (1 << face))
+    {
+      break;
+    }
+  }
+  
+  if (face == 2*N_dims)
+  {
+    /* All faces are done */
+    
+    /* Call the real interpolator */
+    ierr = CCTK_IsFunctionAliased ("DriverInterpolate");
+    if (! ierr)
+    {
+      CCTK_WARN (0, "The aliased function \"DriverInterpolate\" has not been provided");
+    }
+    ierr = DriverInterpolate
+      (cctkGH, N_dims,
+       local_interp_handle, param_table_handle, coord_system_handle,
+       N_interp_points, interp_coords_type, interp_coords,
+       N_input_arrays, input_array_indices,
+       N_output_arrays, output_array_types, output_arrays);
+    
+  }
+  else
+  {
+    /* Treat the face */
+    
+    /* Recursive call to a symmetry condition */
+    if (! symmetry_interpolate[face])
+    {
+      CCTK_WARN (0, "internal error");
+    }
+    ierr = ((SymmetryInterpolateFPointer)symmetry_interpolate[face])
+      (cctkGH, N_dims,
+       local_interp_handle, param_table_handle, coord_system_handle,
+       N_interp_points, interp_coords_type, interp_coords,
+       N_input_arrays, input_array_indices,
+       N_output_arrays, output_array_types, output_arrays,
+       faces);
+    
+  }
+  
+  return ierr;
+}
diff --git a/src/Startup.c b/src/Startup.c
index 298d0ee..886b4fa 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -95,6 +95,7 @@ SymBase_Setup (tFleshConfig * const config,
   struct SymBase *symdata;
   CCTK_INT symmetry_handle[100];
   CCTK_INT symmetry_zone_width[100];
+  CCTK_FPOINTER symmetry_interpolate[100];
   int group;
   int face;
   int ierr;
@@ -111,6 +112,7 @@ SymBase_Setup (tFleshConfig * const config,
   {
     symmetry_handle[face] = -1;
     symmetry_zone_width[face] = 0;
+    symmetry_interpolate[face] = NULL;
   }
 
   /* Create grid symmetry table */
@@ -139,6 +141,13 @@ SymBase_Setup (tFleshConfig * const config,
   {
     CCTK_WARN (0, "Internal Error: Failed to create symmetry_zone_width array");
   }
+  ierr = Util_TableSetFPointerArray
+    (symdata->sym_table,
+     2 * cctkGH->cctk_dim, symmetry_interpolate, "symmetry_interpolate");
+  if (ierr)
+  {
+    CCTK_WARN (0, "Internal Error: Failed to create symmetry_fold_points array");
+  }
 
   /* Create grid array symmetry tables */
   symdata->array_sym_tables
@@ -190,6 +199,14 @@ SymBase_Setup (tFleshConfig * const config,
         {
           CCTK_WARN (0, "Internal Error: Failed to create symmetry_zone_width array");
         }
+        ierr = Util_TableSetFPointerArray
+            (symdata->array_sym_tables[group],
+             2 * CCTK_GroupDimI (group), symmetry_interpolate,
+             "symmetry_interpolate");
+        if (ierr)
+        {
+          CCTK_WARN (0, "Internal Error: Failed to create symmetry_fold_points array");
+        }
         break;
         
       default:
diff --git a/src/Statistics.c b/src/Statistics.c
new file mode 100644
index 0000000..e31cac6
--- /dev/null
+++ b/src/Statistics.c
@@ -0,0 +1,134 @@
+/*@@
+  @file      Statistics.c
+  @author    Erik Schnetter
+  @date      2004-05-25
+  @desc
+             Provide information about the registered symmetry conditions
+  @version   $Header$
+  @enddesc
+@@*/
+
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "util_Table.h"
+
+#include "SymBase.h"
+
+
+
+/* the rcs ID and its dummy function to use it */
+static const char *const rcsid = "$Header$";
+CCTK_FILEVERSION (CactusBase_SymBase_Statistics_c);
+
+
+
+/*@@
+  @routine    SymBase_Statistics
+  @author     Erik Schnetter
+  @date       2004-05-25
+  @desc
+              If verbose, print the symmetry boundary conditions for each
+              face.
+              Warn if symmetries do not have interpolation operators
+              associated with them.
+  @enddesc
+@@*/
+
+void
+SymBase_Statistics (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT sym_table;
+  CCTK_INT symmetry_handle[100];
+  CCTK_FPOINTER symmetry_interpolate[100];
+  char const * symmetry_name;
+  int face;
+  int ierr;
+  
+  
+  
+  /* Get table */
+  sym_table = SymmetryTableHandleForGrid (cctkGH);
+  if (sym_table < 0) CCTK_WARN (0, "internal error");
+
+  /* Get table entries */
+  if (2 * cctkGH->cctk_dim > 100)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  ierr = Util_TableGetIntArray
+    (sym_table, 2 * cctkGH->cctk_dim, symmetry_handle, "symmetry_handle");
+  if (ierr != 2 * cctkGH->cctk_dim)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  ierr = Util_TableGetFPointerArray
+    (sym_table, 2 * cctkGH->cctk_dim,
+     symmetry_interpolate, "symmetry_interpolate");
+  if (ierr != 2 * cctkGH->cctk_dim)
+  {
+    CCTK_WARN (0, "internal error");
+  }
+  
+  
+  
+  /* Print information about the registered symmetry faces, if desired */
+  
+  if (verbose)
+  {
+    
+    /* Loop over all faces */
+    for (face=0; face<2*cctkGH->cctk_dim; ++face) {
+      if (symmetry_handle[face] >= 0) {
+        symmetry_name = SymBase_SymmetryNameOfHandle (symmetry_handle[face]);
+        if (! symmetry_name)
+        {
+          CCTK_WARN (0, "internal error");
+        }
+        if (face >= 6)
+        {
+          CCTK_WARN (0, "only 3 dimensions are supported so far");
+        }
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Symmetry on %s %c-face: %s",
+                    face % 2 == 0 ? "lower" : "upper",
+                    "xyz"[face/2],
+                    symmetry_name);
+      }
+    }
+    
+  }
+  
+  
+  
+  /* Warn about symmetries without interpolators */
+  
+  /* Loop over all registered symmetries */
+  for (face=0; face<2*cctkGH->cctk_dim; ++face) {
+    if (symmetry_handle[face] >= 0) {
+      if (symmetry_interpolate[face] == NULL)
+      {
+        symmetry_name = SymBase_SymmetryNameOfHandle (symmetry_handle[face]);
+        if (! symmetry_name)
+        {
+          CCTK_WARN (0, "internal error");
+        }
+        if (face >= 6)
+        {
+          CCTK_WARN (0, "only 3 dimensions are supported so far");
+        }
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "The symmetry %s on the %s %c-face has not registered a symmetry interpolator",
+                    symmetry_name,
+                    face % 2 == 0 ? "lower" : "upper",
+                    "xyz"[face/2]);
+      }
+    }
+  }
+  
+}
diff --git a/src/SymBase.h b/src/SymBase.h
index 397c0a1..b4e118e 100644
--- a/src/SymBase.h
+++ b/src/SymBase.h
@@ -14,21 +14,60 @@
 
 
 #include "cctk.h"
+#include "cctk_Arguments.h"
 
 
 
 /* SymBase's GH extension */
 struct SymBase
 {
+  /* The table handles below identify tables that have the following
+     entries:
+     
+     CCTK_INT symmetry_handle[2*dim]
+     CCTK_INT symmetry_zone_width[2*dim]
+     CCTK_FPOINTER symmetry_interpolate[2*dim]
+     
+     dim is here the dimension of the group containing the variable,
+     i.e. there is one entry per face.  The true type of the
+     symmetry_interpolate entries is SymmetryInterpolateFPointer.
+     
+     (Given that user code is not allowed to modify these tables, and
+     that this structure is defined here in this header file, there is
+     no real need to use tables for this.  It would be possible to put
+     these data directly into this structure, thereby simplifying user
+     code and not make anything less flexible.)
+  */
+  
   /* Grid symmetry table handle */
   int sym_table;
 
-  /* Grid array symmetry table handles */
+  /* Grid array symmetry table handles, one handle per grid
+     variable */
   int *array_sym_tables;
 };
 
 
 
+typedef
+CCTK_INT
+(* SymmetryInterpolateFPointer) (CCTK_POINTER_TO_CONST 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 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 const output_arrays[],
+                                 CCTK_INT const faces);
+
+
+
 /* Number of registered symmetries */
 extern size_t SymBase_num_symmetries;
 
@@ -73,6 +112,16 @@ SymBase_SymmetryRegisterGN (CCTK_POINTER const cctkGH_,
                             CCTK_INT const *const new_symmetry_zone_width,
                             CCTK_STRING const group_name);
 
+CCTK_INT
+SymBase_SymmetryRegisterInterpolatorFaces (CCTK_INT const sym_table,
+                                           CCTK_INT const group_dim,
+                                           CCTK_INT const sym_handle,
+                                           SymmetryInterpolateFPointer const new_symmetry_interpolate);
+CCTK_INT
+SymBase_SymmetryRegisterGridInterpolator (CCTK_POINTER const cctkGH_,
+                                          CCTK_INT const sym_handle,
+                                          SymmetryInterpolateFPointer const new_symmetry_interpolate);
+
 /* Table.c */
 CCTK_INT
 SymBase_SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST const cctkGH_);
@@ -83,6 +132,42 @@ CCTK_INT
 SymBase_SymmetryTableHandleForGN (CCTK_POINTER_TO_CONST const cctkGH_,
                                   CCTK_STRING const group_name);
 
+/* Interpolate.c */
+CCTK_INT
+SymBase_SymmetryInterpolate (CCTK_POINTER_TO_CONST 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 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 const output_arrays[]);
+CCTK_INT
+SymBase_SymmetryInterpolateFaces (CCTK_POINTER_TO_CONST 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 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 const output_arrays[],
+                                  CCTK_INT const faces);
+
+
+/* Statistics.c */
+void
+SymBase_Statistics (CCTK_ARGUMENTS);
+
 
 
 #endif /* ! defined SYMBASE_H */
diff --git a/src/make.code.defn b/src/make.code.defn
index 84b3c24..53c838c 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Faces.c Handles.c Startup.c Table.c
+SRCS = Faces.c Handles.c Interpolation.c Startup.c Statistics.c Table.c
 
 # Subdirectories containing source files
 SUBDIRS =