Provide symmetry interpolation.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/RotatingSymmetry90/trunk@9 c3c03602-0f4f-0410-b3fa-d2c81c8a7dc5

5 files changed, 810 insertions, 13 deletions

diff --git a/interface.ccl b/interface.ccl
index 46383d9..a011f14 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -10,15 +10,54 @@ USES INCLUDE HEADER: Slab.h
 
 
 CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name)
-USES FUNCTION SymmetryRegister
+REQUIRES FUNCTION SymmetryRegister
 
-CCTK_INT FUNCTION                                  \
-    SymmetryRegisterGrid                           \
-         (CCTK_POINTER IN cctkGH,                  \
-          CCTK_INT IN sym_handle,                  \
-          CCTK_INT IN ARRAY which_faces,           \
-          CCTK_INT IN ARRAY symmetry_zone_width)
-USES FUNCTION SymmetryRegisterGrid
+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
 
 
 
@@ -30,4 +69,4 @@ CCTK_INT FUNCTION Boundary_SelectedGVs \
    CCTK_INT ARRAY OUT boundary_widths, \
    CCTK_INT ARRAY OUT table_handles, \
    CCTK_STRING IN  bc_name)
-USES FUNCTION Boundary_SelectedGVs
+REQUIRES FUNCTION Boundary_SelectedGVs
diff --git a/src/interpolate.c b/src/interpolate.c
new file mode 100644
index 0000000..2d6c369
--- /dev/null
+++ b/src/interpolate.c
@@ -0,0 +1,749 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "rotatingsymmetry90.h"
+
+
+
+/* A tensor description */
+struct tensor
+{
+  const char * name;            /* description */
+  int dim;
+  int rank;
+  int ncomps;                   /* dim^rank */
+  int * restrict vars;          /* map component to variable */
+  int * restrict parity;        /* parity for the above */
+  int nvars;                    /* depends on symmetries */
+  int * restrict comps;         /* map variable to component */
+};
+
+
+
+/* A scalar */
+static int scalar_vars[] = {
+  0
+};
+static int scalar_parity[] = {
+  +1
+};
+static int scalar_comps[] = {
+  0
+};
+static struct tensor scalar = {
+  "scalar",
+  3, 0, 1, scalar_vars, scalar_parity, 1, scalar_comps
+};
+
+
+
+/* A vector */
+static int vector_vars[] = {
+  0,1,2
+};
+static int vector_parity[] = {
+  +1,+1,+1
+};
+static int vector_comps[] = {
+  0,1,2
+};
+static struct tensor vector = {
+  "vector",
+  3, 1, 3, vector_vars, vector_parity, 3, vector_comps
+};
+
+
+
+/* A second rank tensor without symmetries */
+static int tensor_vars[] = {
+  0,1,2,   3,4,5,   6,7,8
+};
+static int tensor_parity[] = {
+  +1,+1,+1,   +1,+1,+1,   +1,+1,+1
+};
+static int tensor_comps[] = {
+  0,1,2,   3,4,5,   6,7,8
+};
+static struct tensor tensor = {
+  "tensor",
+  3, 2, 9, tensor_vars, tensor_parity, 9, tensor_comps
+};
+
+
+
+/* A symmetric second rank tensor */
+static int symmtensor_vars[] = {
+  0,1,2,   1,3,4,   2,4,5
+};
+static int symmtensor_parity[] = {
+  +1,+1,+1,   +1,+1,+1,   +1,+1,+1
+};
+static int symmtensor_comps[] = {
+  0,1,2,   4,5,   8
+};
+static struct tensor symmtensor = {
+  "symmetric tensor T_(ij)",
+  3, 2, 9, symmtensor_vars, symmtensor_parity, 6, symmtensor_comps
+};
+
+
+
+/* A third rank tensor with symmetry in the first two indices */
+static int symmtensor3a_vars[] = {
+   0, 1, 2,
+   3, 4, 5,
+   6, 7, 8,
+  
+   3, 4, 5,
+   9,10,11,
+  12,13,14,
+  
+   6, 7, 8,
+  12,13,14,
+  15,16,17
+};
+static int symmtensor3a_parity[] = {
+  +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,
+};
+static int symmtensor3a_comps[] = {
+   0, 1, 2,
+   3, 4, 5,
+   6, 7, 8,
+  
+  12,13,14,
+  15,16,17,
+  
+  24,25,26
+};
+static struct tensor symmtensor3a = {
+  "symmetric tensor T_(ij)k",
+  3, 3, 27, symmtensor3a_vars, symmtensor3a_parity, 18, symmtensor3a_comps
+};
+
+
+
+/* A third rank tensor with symmetry in the last two indices */
+static int symmtensor3b_vars[] = {
+   0, 1, 2,    1, 3, 4,    2, 4, 5,
+   6, 7, 8,    7, 9,10,    8,10,11,
+  12,13,14,   13,15,16,   14,16,17
+};
+static int symmtensor3b_parity[] = {
+  +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
+};
+static int symmtensor3b_comps[] = {
+   0, 1, 2,    4, 5,    8,
+   9,10,11,   13,14,   17,
+  18,19,20,   22,23,   26
+};
+static struct tensor symmtensor3b = {
+  "symmetric tensor T_i(jk)",
+  3, 3, 27, symmtensor3b_vars, symmtensor3b_parity, 18, symmtensor3b_comps
+};
+
+
+
+/* A fourth rank tensor with symmetries both in its first and last two
+   indices */
+static int symmtensor4_vars[] = {
+   0, 1, 2,    1, 3, 4,    2, 4, 5,
+   6, 7, 8,    7, 9,10,    8,10,11,
+  12,13,14,   13,15,16,   14,16,17,
+  
+   6, 7, 8,    7, 9,10,    8,10,11,
+  18,19,20,   19,21,22,   20,22,23,
+  24,25,26,   25,27,28,   26,28,29,
+  
+  12,13,14,   13,15,16,   14,16,17,
+  24,25,26,   25,27,28,   26,28,29,
+  30,31,32,   31,33,34,   32,34,35
+};
+static int symmtensor4_parity[] = {
+  +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,   +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,   +1,+1,+1,
+  +1,+1,+1,   +1,+1,+1,   +1,+1,+1,
+  +1,+1,+1,   +1,+1,+1,   +1,+1,+1
+};
+static int symmtensor4_comps[] = {
+   0, 1, 2,    4, 5,    8,
+   9,10,11,   13,14,   17,
+  18,19,20,   22,23,   26,
+  
+  36,37,38,   40,41,   44,
+  45,46,47,   49,50,   53,
+  
+  72,73,74,   76,77,   80
+};
+static struct tensor symmtensor4 = {
+  "symmetric tensor T_(ij)(kl)",
+  3, 4, 81, symmtensor4_vars, symmtensor4_parity, 36, symmtensor4_comps
+};
+
+
+
+/* Ensure that all tensor declarations are internally consistent */
+static void
+CheckTensorType (struct tensor const * restrict const atensor)
+{
+  int i, n;
+  assert (atensor->name);
+  assert (atensor->dim>=0);
+  assert (atensor->rank>=0);
+  assert (atensor->ncomps>=0);
+  assert (atensor->ncomps == floor(pow(atensor->dim, atensor->rank) + 0.5));
+  assert (atensor->nvars>=0 && atensor->nvars<=atensor->ncomps);
+  assert (atensor->vars);
+  for (i=0; i<atensor->ncomps; ++i) {
+    assert (atensor->vars[i]>=0 && atensor->vars[i]<atensor->nvars);
+    assert (abs(atensor->parity[i]) <= 1);
+  }
+  assert (atensor->comps);
+  for (n=0; n<atensor->nvars; ++n) {
+    assert (atensor->comps[n]>=0 && atensor->comps[n]<atensor->ncomps);
+    assert (atensor->vars[atensor->comps[n]] == n);
+  }
+}
+
+void
+Rot90_CheckTensorTypes (CCTK_ARGUMENTS)
+{
+  CheckTensorType (&scalar);
+  CheckTensorType (&vector);
+  CheckTensorType (&tensor);
+  CheckTensorType (&symmtensor);
+  CheckTensorType (&symmtensor3a);
+  CheckTensorType (&symmtensor3b);
+  CheckTensorType (&symmtensor4);
+}
+  
+
+
+/* Symmetry interpolation */
+CCTK_INT
+Rot90_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 const faces)
+{
+  cGH const * restrict const cctkGH = cctkGH_;
+  
+  CCTK_POINTER new_interp_coords[3];
+  CCTK_POINTER * restrict new_output_arrays;
+  CCTK_INT new_faces;
+  CCTK_INT * restrict operand_indices;
+  CCTK_INT * restrict operation_codes;
+  CCTK_INT * restrict output_array_indices;
+  
+  struct tensor const * restrict * restrict thetensor[3];
+  int * restrict * restrict thevars[3];
+  int * restrict thebase;
+  int * restrict thevar;
+  
+  int m;                        /* output array */
+  int n;                        /* point */
+  int i;                        /* var */
+  int f;                        /* face */
+  int d;                        /* dimension */
+  int q;                        /* number of derivatives */
+  int r;                        /* rank */
+  
+  int iret;                     /* interpolator return value */
+  int ierr;
+  
+  /* Check arguments */
+  assert (N_dims == 3);
+  assert (N_interp_points >= 0);
+  assert (interp_coords_type >= 0);
+  for (d=0; d<3; ++d) {
+    assert (N_interp_points == 0 || interp_coords[d]);
+  }
+  assert (N_output_arrays >= 0);
+  
+  /* Coordinates must be CCTK_REAL */
+  assert (interp_coords_type == CCTK_VARIABLE_REAL);
+  for (m=0; m<N_output_arrays; ++m) {
+    assert (output_array_types[m] == CCTK_VARIABLE_REAL);
+  }
+  
+  
+  
+  /* Claim faces */
+  assert (faces & (1 << 0));
+  assert (faces & (1 << 2));
+  new_faces = faces;
+  new_faces &= ~ (1 << 0);
+  new_faces &= ~ (1 << 2);
+  
+  
+  
+  /* Copy coordinates */
+  for (d=0; d<3; ++d) {
+    new_interp_coords[d] = malloc (N_interp_points * sizeof(CCTK_REAL));
+    assert (new_interp_coords[d]);
+  }
+  
+  /* Fold coordinates */
+  for (n=0; n<N_interp_points; ++n) {
+    /* Is the point outside the domain? */
+    CCTK_REAL x = ((CCTK_REAL const *)interp_coords[0])[n];
+    CCTK_REAL y = ((CCTK_REAL const *)interp_coords[1])[n];
+    CCTK_REAL z = ((CCTK_REAL const *)interp_coords[2])[n];
+    if (x < 0) {
+      /* Rotate the point by 180 degrees */
+      x = -x;
+      y = -y;
+    }
+    if (y < 0) {
+      /* Rotate the point clockwise by 90 degrees */
+      CCTK_REAL const t = x;
+      x = -y;
+      y =  t;
+    }
+    ((CCTK_REAL *)new_interp_coords[0])[n] = x;
+    ((CCTK_REAL *)new_interp_coords[1])[n] = y;
+    ((CCTK_REAL *)new_interp_coords[2])[n] = z;
+  }
+  
+  
+  
+  /* Allocate new output arrays */
+  new_output_arrays = malloc (N_output_arrays * sizeof *new_output_arrays);
+  assert (new_output_arrays);
+  for (m=0; m<N_output_arrays; ++m) {
+    new_output_arrays[m] = malloc (N_interp_points * sizeof(CCTK_REAL));
+    assert (new_output_arrays[m]);
+  }
+  
+  
+  
+  /* 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, new_output_arrays,
+     new_faces);
+  
+  
+  
+  /* Free coordinates */
+  for (d=0; d<3; ++d) {
+    free (new_interp_coords[d]);
+  }
+  
+  
+  
+  /* 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;   /* output index equals 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());
+  }
+  
+  
+  
+  /* Map Cactus variables to tensor objects */
+  for (q=0; q<3; ++q) {
+    thetensor[q] = malloc (CCTK_NumVars() * sizeof *thetensor[q]);
+    assert (thetensor[q]);
+    thevars[q] = malloc (CCTK_NumVars() * sizeof *thevars[q]);
+    assert (thevars[q]);
+    for (n=0; n<CCTK_NumVars(); ++n) {
+      thetensor[q][n] = NULL;
+      thevars[q][n] = NULL;
+    }
+  }
+  
+  /* Map output arrays to the base Cactus variable (i.e. the first
+     component in the tensor objects) */
+  thebase = malloc (N_output_arrays * sizeof *thebase);
+  assert (thebase);
+  thevar = malloc (N_output_arrays * sizeof *thevar);
+  assert (thevar);
+  
+  for (m=0; m<N_output_arrays; ++m) {
+    
+    int vi, gi;
+    int numvars, firstvar;
+    int table;
+    
+    char tensortypealias[100];
+    
+    struct tensor const * tensortype;
+    int basevar;
+    int var;
+    int indices[10];
+    int oldrank;
+    int num_derivs;
+    
+    /* Get some variable information */
+    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);
+    table = CCTK_GroupTagsTableI(gi);
+    assert (table>=0);
+    
+    /* Get the tensor type alias */
+    ierr = Util_TableGetString
+      (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
+    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+      char * 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) {
+      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);
+    }
+    
+    /* Find the tensor type */
+    tensortype = NULL;
+    basevar = -1;
+    var = -1;
+    if (CCTK_EQUALS (tensortypealias, "scalar")) {
+      /* scalar */
+      if (numvars != 1) {
+        char * 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);
+      }
+      tensortype = &scalar;
+      basevar = vi;
+      var = 0;
+    } else if (CCTK_EQUALS (tensortypealias, "u")) {
+      /* vector */
+      assert (numvars == 3);
+      tensortype = &vector;
+      basevar = firstvar;
+      var = vi - firstvar;
+    } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
+      /* symmetric tensor */
+      assert (numvars == 6);
+      tensortype = &symmtensor;
+      basevar = firstvar;
+      var = 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);
+    }
+    thebase[m] = basevar;
+    
+    /* Find my component */
+    {
+      int component;
+      assert (var>=0 && var<tensortype->nvars);
+      component = tensortype->comps[var];
+      assert (tensortype->rank<10);
+      for (r=tensortype->rank-1; r>=0; --r) {
+        indices[r] = component % tensortype->dim;
+        assert (indices[r]>=0 && indices[r]<tensortype->dim);
+        component /= tensortype->dim;
+      }
+      assert (component == 0);
+      oldrank = tensortype->rank;
+    }
+    
+    /* Take operation code (i.e. derivatives) into account */
+    {
+      int code = operation_codes[m];
+      num_derivs
+        = code == 0 ? 0
+        : code < 10 ? 1
+        : code < 100 ? 2
+        : -1;
+      assert (num_derivs>=0);
+    }
+    
+    if (tensortype == &scalar) {
+      switch (num_derivs) {
+      case 0: break;
+      case 1: tensortype = &vector; break;
+      case 2: tensortype = &symmtensor; break;
+      default: assert (0);
+      }
+    } else if (tensortype == &vector) {
+      switch (num_derivs) {
+      case 0: break;
+      case 1: tensortype = &tensor; break;
+      case 2: tensortype = &symmtensor3b; break;
+      default: assert (0);
+      }
+    } else if (tensortype == &symmtensor) {
+      switch (num_derivs) {
+      case 0: break;
+      case 1: tensortype = &symmtensor3a; break;
+      case 2: tensortype = &symmtensor4; break;
+      default: assert (0);
+      }
+    } else {
+      assert (0);
+    }
+    
+    /* Find the additional indices */
+    {
+      int code;
+      assert (tensortype->rank<10);
+      assert (tensortype->rank >= oldrank);
+      code = operation_codes[m];
+      for (r=oldrank; r<tensortype->rank; ++r) {
+        const int thedir = code % 10 - 1;
+        code /= 10;
+        assert (thedir>=0 && thedir<tensortype->dim);
+        indices[r] = thedir;
+      }
+    }
+    
+    /* Re-calculate component */
+    {
+      int component = 0;
+      for (r=0; r<tensortype->rank; ++r) {
+        component = component * tensortype->dim + indices[r];
+      }
+      assert (component>=0 && component<tensortype->ncomps);
+      var = tensortype->vars[component];
+      assert (var>=0 && var<tensortype->nvars);
+      thevar[m] = var;
+    }
+    
+    /* Create or cross-check the tensor object */
+    if (! thetensor[num_derivs][basevar]) {
+      thetensor[num_derivs][basevar] = tensortype;
+      assert (! thevars[num_derivs][basevar]);
+      thevars[num_derivs][basevar]
+        = malloc (tensortype->nvars * sizeof *thevars[num_derivs][basevar]);
+      assert (thevars[num_derivs][basevar]);
+      for (i=0; i<tensortype->nvars; ++i) {
+        thevars[num_derivs][basevar][i] = -1;
+      }
+    }
+    assert (thetensor[num_derivs][basevar] == tensortype);
+    assert (thevars[num_derivs][basevar]);
+    assert (thevars[num_derivs][basevar][var] == -1);
+    thevars[num_derivs][basevar][var] = m;
+    
+  } /* for m */
+  
+  
+  
+  /* Loop over all output arrays */
+  for (m=0; m<N_output_arrays; ++m) {
+    
+    int num_derivs;
+    int basevar;
+    struct tensor const * restrict tensortype;
+    int const * restrict vars;
+    int var;
+    int indices[10];
+    
+    /* Take operation code (i.e. derivatives) into account */
+    {
+      int code = operation_codes[m];
+      num_derivs
+        = code == 0 ? 0
+        : code < 10 ? 1
+        : code < 100 ? 2
+        : -1;
+      assert (num_derivs>=0);
+    }
+    
+    /* Get the tensor type */
+    basevar = thebase[m];
+    assert (basevar>=0 && basevar<=CCTK_NumVars());
+    tensortype = thetensor[num_derivs][basevar];
+    assert (tensortype);
+    vars = thevars[num_derivs][basevar];
+    assert (vars);
+    var = thevar[m];
+    assert (var>=0 && var<tensortype->nvars);
+    
+    /* Transform into indices */
+    {
+      int component = tensortype->comps[var];
+      assert (tensortype->rank<10);
+      for (r=tensortype->rank-1; r>=0; --r) {
+        indices[r] = component % tensortype->dim;
+        assert (indices[r]>=0 && indices[r]<tensortype->dim);
+        component /= tensortype->dim;
+      }
+    }
+    
+    /* Loop over all grid points */
+    for (n=0; n<N_interp_points; ++n) {
+      
+      CCTK_REAL x = ((CCTK_REAL const *)interp_coords[0])[n];
+      CCTK_REAL y = ((CCTK_REAL const *)interp_coords[1])[n];
+      
+      int myindices[10];
+      int parity;
+      
+      for (r=0; r<tensortype->rank; ++r) {
+        myindices[r] = indices[r];
+      }
+      parity = +1;
+      
+      /* Is the point outside the domain? */
+      if (x<0) {
+        /* Undo 180 degree rotation */
+        for (r=0; r<tensortype->rank; ++r) {
+          switch (myindices[r]) {
+          case 0: parity *= -1; break;
+          case 1: parity *= -1; break;
+          case 2: break;
+          default: assert (0);
+          }
+        }
+        x = -x;
+        y = -y;
+      }
+      if (y<0) {
+        /* Undo clockwise 90 degree rotation (i.e. rotate
+           counterclockwise) */
+        for (r=0; r<tensortype->rank; ++r) {
+          switch (myindices[r]) {
+          case 0: myindices[r]=1; parity *= +1; break;
+          case 1: myindices[r]=0; parity *= -1; break;
+          case 2: break;
+          default: assert (0);
+          }
+        }
+        {
+          CCTK_REAL const t = x;
+          x = -y;
+          y =  t;
+        }
+      }
+      
+      {
+        int mycomponent;
+        int myvar;
+        int mm;
+        
+        mycomponent = 0;
+        for (r=0; r<tensortype->rank; ++r) {
+          mycomponent = mycomponent * tensortype->dim + myindices[r];
+        }
+        assert (mycomponent>=0 && mycomponent<tensortype->ncomps);
+        myvar = tensortype->vars[mycomponent];
+        assert (myvar>=0 && myvar<tensortype->nvars);
+        mm = vars[myvar];
+        assert (mm>=0 && mm<N_output_arrays);
+        
+        ((CCTK_REAL *)output_arrays[m])[n]
+          = parity * ((CCTK_REAL const *)new_output_arrays[mm])[n];
+      }
+      
+    } /* for n */
+    
+  } /* for m */
+  
+  
+  
+  /* Free output variable descriptions */
+  for (m=0; m<N_output_arrays; ++m) {
+    free (new_output_arrays[m]);
+  }
+  free (new_output_arrays);
+  
+  free (operand_indices);
+  free (operation_codes);
+  free (output_array_indices);
+  
+  for (q=0; q<3; ++q) {
+    free ((void *) (thetensor[q]));
+    free ((void *) (thevars[q]));
+  }
+  free (thebase);
+  free (thevar);
+  
+  
+
+  return iret;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index 4b71cce..5907239 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = rotatingsymmetry90.c registersymmetry.c
+SRCS = interpolate.c rotatingsymmetry90.c registersymmetry.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/registersymmetry.c b/src/registersymmetry.c
index da80afa..0d2d7b0 100644
--- a/src/registersymmetry.c
+++ b/src/registersymmetry.c
@@ -4,6 +4,9 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
+#include "rotatingsymmetry90.h"
+
+
 void Rot90_RegisterSymmetry (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
@@ -22,8 +25,8 @@ void Rot90_RegisterSymmetry (CCTK_ARGUMENTS)
   
   faces[0] = 1;
   width[0] = cctk_nghostzones[0];
-  faces[1] = 1;
-  width[1] = cctk_nghostzones[1];
+  faces[2] = 1;
+  width[2] = cctk_nghostzones[1];
   
   handle = SymmetryRegister ("rotating_symmetry_90");
   if (handle < 0) {
@@ -34,4 +37,10 @@ void Rot90_RegisterSymmetry (CCTK_ARGUMENTS)
   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, Rot90_SymmetryInterpolate);
+  if (ierr < 0) {
+    CCTK_WARN (0, "Could not register the symmetry interpolator");
+  }
 }
diff --git a/src/rotatingsymmetry90.h b/src/rotatingsymmetry90.h
index bdc3627..f7b408d 100644
--- a/src/rotatingsymmetry90.h
+++ b/src/rotatingsymmetry90.h
@@ -28,6 +28,6 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
                            CCTK_INT const N_output_arrays,
                            CCTK_INT const output_array_types[],
                            CCTK_POINTER const output_arrays[],
-                           CCTK_INT const faces[]);
+                           CCTK_INT const faces);
 
 #endif /* ! defined ROTATINGSYMMETRY90_H */