Add support for 4-tensors.

Patch mostly thanks to Bela Szilagyi.


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

2 files changed, 304 insertions, 297 deletions

diff --git a/src/apply.c b/src/apply.c
index 1f71e1b..918957e 100644
--- a/src/apply.c
+++ b/src/apply.c
@@ -47,12 +47,9 @@ CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_apply_c);
     assert(jjmax>=-1&&jjmin<nj);                                        \
     assert(kkmax>=-1&&kkmin<nk);                                        \
                                                                         \
-    for (k=kmin; k<kmax; ++k)                                           \
-    {                                                                   \
-      for (j=jmin; j<jmax; ++j)                                         \
-      {                                                                 \
-        for (i=imin; i<imax; ++i)                                       \
-        {                                                               \
+    for (k=kmin; k<kmax; ++k) {                                         \
+      for (j=jmin; j<jmax; ++j) {                                       \
+        for (i=imin; i<imax; ++i) {                                     \
           int const dstind = i + ni * (j + nj * k);                     \
           int const ii = ioff + idir * i;                               \
           int const jj = joff + jdir * j;                               \
@@ -167,24 +164,19 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   
   /* Check arguments */
-  if (! cctkGH)
-  {
+  if (! cctkGH) {
     CCTK_WARN (0, "Argument cctkGH is NULL");
   }
-  if (! stencil)
-  {
+  if (! stencil) {
     CCTK_WARN (0, "Argument stencil is NULL");
   }
-  if (vi < 0 || vi >= CCTK_NumVars())
-  {
+  if (vi < 0 || vi >= CCTK_NumVars()) {
     CCTK_WARN (0, "Illegal variable index");
   }
   
-  if (verbose)
-  {
+  if (verbose) {
     fullname = CCTK_FullName (vi);
-    if (! fullname)
-    {
+    if (! fullname)   {
       CCTK_WARN (0, "Internal error in CCTK_FullName");
     }
     CCTK_VInfo (CCTK_THORNSTRING,
@@ -197,8 +189,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   /* Get and check group information */
   gi = CCTK_GroupIndexFromVarI (vi);
-  if (gi < 0 || gi > CCTK_NumGroups())
-  {
+  if (gi < 0 || gi > CCTK_NumGroups()) {
     CCTK_WARN (0, "Internal error in CCTK_GroupIndexFromVarI");
   }
   
@@ -258,11 +249,14 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   ttype = UNKNOWN;
   tcomponent = 0;
-  if (CCTK_EQUALS (tensortypealias, "scalar"))
-  {
+  if (CCTK_EQUALS (tensortypealias, "scalar")) {
     /* scalar */
     ttype = SCALAR;
     tcomponent = 0;
+  } else if (CCTK_EQUALS (tensortypealias, "4scalar")) {
+    /* 4-scalar */
+    ttype = SCALAR;
+    tcomponent = 0;
   } else if (CCTK_EQUALS (tensortypealias, "u")
              || CCTK_EQUALS (tensortypealias, "d"))
   {
@@ -270,12 +264,40 @@ BndReflectVI (cGH const * restrict const cctkGH,
     assert (numvars == 3);
     ttype = VECTOR;
     tcomponent = vi - firstvar;
+  } else if (CCTK_EQUALS (tensortypealias, "4u")
+             || CCTK_EQUALS (tensortypealias, "4d"))
+  {
+    /* 4-vector */
+    assert (numvars == 4);
+    if (vi == firstvar) {
+      ttype = SCALAR;
+      tcomponent = 0;
+    } else {
+      ttype = VECTOR;
+      tcomponent = vi - firstvar - 1;
+    }
   } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
-             || CCTK_EQUALS (tensortypealias, "dd_sym")) {
+             || CCTK_EQUALS (tensortypealias, "dd_sym"))
+  {
     /* symmetric tensor */
     assert (numvars == 6);
     ttype = TENSOR;
     tcomponent = vi - firstvar;
+  } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
+             || CCTK_EQUALS (tensortypealias, "4dd_sym"))
+  {
+    /* symmetric 4-tensor */
+    assert (numvars == 10);
+    if (vi == firstvar) {
+      ttype = SCALAR;
+      tcomponent = 0;
+    } else if (vi <= firstvar+3) {
+      ttype = VECTOR;
+      tcomponent = vi - firstvar - 1;
+    } else {
+      ttype = TENSOR;
+      tcomponent = vi - firstvar - 4;
+    }
   } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
     /* Weyl scalars, stored as 10 real values.  NOTE: This assumes
        that Psi_0 comes first, which is NOT the default with
@@ -292,8 +314,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
     free (groupname);
   }
   
-  switch (ttype)
-  {
+  switch (ttype) {
   case SCALAR:
     assert (tcomponent>=0 && tcomponent<1);
     break;
@@ -342,21 +363,16 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   
   /* Loop over all directions and faces */
-  for (dir=0; dir<3; ++dir)
-  {
-    for (face=0; face<2; ++face)
-    {
+  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 (do_reflection[2*dir+face]) {
         /* If we have the outer boundary of that face */
-        if (cctkGH->cctk_bbox[2*dir+face])
-        {
+        if (cctkGH->cctk_bbox[2*dir+face]) {
           
           /* Find parity */
           parity = tensorparity;
-          switch (ttype)
-          {
+          switch (ttype) {
           case SCALAR:
             parity *= +1;
             break;
@@ -364,8 +380,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
             parity *= dir == tcomponent ? -1 : +1;
             break;
           case TENSOR:
-            switch (tcomponent)
-            {
+            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;
@@ -375,44 +390,39 @@ BndReflectVI (cGH const * restrict const cctkGH,
             default: assert (0);
             }
             break;
-          case WEYLSCALARS_REAL:
-            {
-              static int const weylparities[10][3] =
-                {{+1,+1,+1},
-                 {-1,-1,-1},
-                 {+1,+1,+1},
-                 {-1,-1,-1},
-                 {+1,+1,+1},
-                 {-1,-1,-1},
-                 {+1,+1,+1},
-                 {-1,-1,-1},
+          case WEYLSCALARS_REAL: {
+            static int const weylparities[10][3] =
+              {{+1,+1,+1},
+               {-1,-1,-1},
+               {+1,+1,+1},
+               {-1,-1,-1},
+               {+1,+1,+1},
+               {-1,-1,-1},
+               {+1,+1,+1},
+               {-1,-1,-1},
                  {+1,+1,+1},
-                 {-1,-1,-1}};
-              parity *= weylparities[tcomponent][dir];
-            }
+               {-1,-1,-1}};
+            parity *= weylparities[tcomponent][dir];
             break;
+          }
           default:
             assert (0);
           }
           
           /* Find region extent */
-          for (d=0; d<3; ++d)
-          {
+          for (d=0; d<3; ++d) {
             lsh[d] = cctkGH->cctk_lsh[d];
             imin[d] = 0;
             imax[d] = cctkGH->cctk_lsh[d];
             ioff[d] = 0;
             idir[d] = 1;
           }
-          if (face == 0)
-          {
+          if (face == 0) {
             imax[dir] = cctkGH->cctk_nghostzones[dir];
             ioff[dir] = (+ 2*cctkGH->cctk_nghostzones[dir] - 1 
                        + (do_stagger[2*dir+face] ? 0 : 1));
             idir[dir] = -1;
-          }
-          else
-          {
+          } else {
             imin[dir] = cctkGH->cctk_lsh[dir] - cctkGH->cctk_nghostzones[dir];
             ioff[dir] = (- 2*cctkGH->cctk_nghostzones[dir]
                        - (do_stagger[2*dir+face] ? 0 : 1));
@@ -420,8 +430,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
           }
           
           /* Copy region */
-          switch (group.vartype)
-          {
+          switch (group.vartype) {
             
 #ifdef HAVE_CCTK_INT1
           case CCTK_VARIABLE_INT1:
@@ -725,86 +734,71 @@ ReflectionSymmetry_Apply (CCTK_ARGUMENTS)
   int istat;
   int ierr;
   
-  if (! cctkGH)
-  {
+  if (! cctkGH) {
     CCTK_WARN (0, "Argument cctkGH is NULL");
   }
   
   nvars = Boundary_SelectedGVs (cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
-  if (nvars < 0)
-  {
+  if (nvars < 0) {
     CCTK_WARN (0, "Internal error in Boundary_SelectedGVs");
   }
   
-  if (nvars == 0)
-  {
+  if (nvars == 0) {
     /* Nothing to do */
     return;
   }
   
   indices = malloc (nvars * sizeof *indices);
-  if (! indices)
-  {
+  if (! indices) {
     CCTK_WARN (0, "Out of memory");
   }
   faces = malloc (nvars * sizeof *faces);
-  if (! faces)
-  {
+  if (! faces) {
     CCTK_WARN (0, "Out of memory");
   }
   widths = malloc (nvars * sizeof *widths);
-  if (! widths)
-  {
+  if (! widths) {
     CCTK_WARN (0, "Out of memory");
   }
   tables = malloc (nvars * sizeof *tables);
-  if (! tables)
-  {
+  if (! tables) {
     CCTK_WARN (0, "Out of memory");
   }
   
   istat =  Boundary_SelectedGVs
     (cctkGH, nvars, indices, faces, widths, tables, 0);
-  if (istat != nvars)
-  {
+  if (istat != nvars) {
     CCTK_WARN (0, "Internal error in Boundary_SelectedGVs");
   }
   
-  for (i=0; i<nvars; ++i)
-  {
+  for (i=0; i<nvars; ++i) {
     vi = indices[i];
-    if (vi < 0 || vi >= CCTK_NumVars())
-    {
+    if (vi < 0 || vi >= CCTK_NumVars()) {
       CCTK_WARN (0, "Illegal variable index");
     }
     
-    if (widths[i] < 0)
-    {
+    if (widths[i] < 0) {
       CCTK_WARN (0, "Illegal boundary width");
     }
     
     dim = CCTK_GroupDimFromVarI (vi);
-    if (dim < 0)
-    {
+    if (dim < 0) {
       CCTK_WARN (0, "Illegal dimension");
     }
     
     stencil = malloc (dim * sizeof *stencil);
-    if (! stencil)
-    {
+    if (! stencil) {
       CCTK_WARN (0, "Out of memory");
     }
     ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi);
-    if (ierr)
-    {
+    if (ierr) {
       CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI");
     }
     
     CheckBoundaryParameters (cctkGH, vi, stencil);
     
     ierr = BndReflectVI (cctkGH, stencil, vi);
-    if (ierr)
-    {
+    if (ierr) {
       CCTK_WARN (0, "Internal error in BndReflectVI");
     }
     
diff --git a/src/interpolate.c b/src/interpolate.c
index e4c242d..eb60460 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -67,10 +67,8 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
   do_reflection[5] = reflection_upper_z;
   
   newfaces = faces;
-  for (d=0; d<6; ++d)
-  {
-    if (do_reflection[d])
-    {
+  for (d=0; d<6; ++d) {
+    if (do_reflection[d]) {
       assert (newfaces & (1 << d));
       newfaces &= ~ (1 << d);
     }
@@ -80,25 +78,20 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
   
   /* Fold coordinates */
   assert (interp_coords_type == CCTK_VARIABLE_REAL);
-  for (dir=0; dir<3; ++dir)
-  {
+  for (dir=0; dir<3; ++dir) {
     assert (! do_reflection[2*dir+1]);
     
-    if (do_reflection[2*dir])
-    {
+    if (do_reflection[2*dir]) {
       new_interp_coords[dir]
         = malloc (N_interp_points * sizeof (CCTK_REAL));
       assert (N_interp_points == 0 || new_interp_coords[dir]);
       
-      for (n=0; n<N_interp_points; ++n)
-      {
+      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
-    {
+    } else {
       new_interp_coords[dir] = interp_coords[dir];
     }
   }
@@ -117,15 +110,14 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
   
   
   /* Free coordinates */
-  for (dir=0; dir<3; ++dir)
-  {
-    if (do_reflection[2*dir])
-    {
+  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);
@@ -167,158 +159,185 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
   
   
   /* Unfold tensor types */
-  for (m=0; m<N_output_arrays; ++m)
-  {
-  if (output_array_indices[m]!=-1)
-  {
-    
-    int vi, gi;
-    int numvars, firstvar;
-    char * groupname;
-    
-    int table;
-    char tensortypealias[1000];
-    enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL };
-    enum tensortype ttype;
-    CCTK_INT tensorparity;
-    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);
-    table = CCTK_GroupTagsTableI(gi);
-    assert (table>=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) {
+  for (m=0; m<N_output_arrays; ++m) {
+    if (output_array_indices[m]!=-1) {
+      
+      int vi, gi;
+      int numvars, firstvar;
+      char * groupname;
+      
+      int table;
+      char tensortypealias[1000];
+      enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL };
+      enum tensortype ttype;
+      CCTK_INT tensorparity;
+      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);
+      table = CCTK_GroupTagsTableI(gi);
+      assert (table>=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 (2, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element",
+        CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Error in tensor type alias declaration for group \"%s\"",
                     groupname);
         free (groupname);
       }
-      ttype = SCALAR;
+      
+      ttype = UNKNOWN;
       tcomponent = 0;
-    }
-    else if (CCTK_EQUALS (tensortypealias, "u")
-             || CCTK_EQUALS (tensortypealias, "d")) {
-      /* vector */
-      assert (numvars == 3);
-      ttype = VECTOR;
-      tcomponent = vi - firstvar;
-    } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
-               || CCTK_EQUALS (tensortypealias, "dd_sym")) {
-      /* symmetric tensor */
-      assert (numvars == 6);
-      ttype = TENSOR;
-      tcomponent = vi - firstvar;
-    } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
-      /* Weyl scalars, stored as 10 real numbers */
-      assert (numvars == 10);
-      ttype = WEYLSCALARS_REAL;
-      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;
-    case WEYLSCALARS_REAL:
-      assert (tcomponent>=0 && tcomponent<10);
-      break;
-    default:
-      assert (0);
-    }
-    
-    ierr = Util_TableGetInt (table, & tensorparity, "tensorparity");
-    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
-      tensorparity = +1;
-    } else if (ierr<0) {
-      groupname = CCTK_GroupName(gi);
-      assert (groupname);
-      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Error in tensor parity declaration for group \"%s\"",
-                  groupname);
-      free (groupname);
-    }
-    
-    
-    
-    /* 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)
+      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, "4scalar")) {
+        /* 4-scalar */
+        assert (numvars == 1);
+        ttype = SCALAR;
+        tcomponent = 0;
+      } else if (CCTK_EQUALS (tensortypealias, "u")
+                 || CCTK_EQUALS (tensortypealias, "d"))
       {
-      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;
-    case WEYLSCALARS_REAL:
+        /* vector */
+        assert (numvars == 3);
+        ttype = VECTOR;
+        tcomponent = vi - firstvar;
+      } else if (CCTK_EQUALS (tensortypealias, "4u")
+                 || CCTK_EQUALS (tensortypealias, "4d"))
+      {
+        /* 4-vector */
+        assert (numvars == 4);
+        if (vi == firstvar) {
+          ttype = SCALAR;
+          tcomponent = 0;
+        } else {
+          ttype = VECTOR;
+          tcomponent = vi - firstvar - 1;
+        }
+      } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
+                 || CCTK_EQUALS (tensortypealias, "dd_sym"))
+      {
+        /* symmetric tensor */
+        assert (numvars == 6);
+        ttype = TENSOR;
+        tcomponent = vi - firstvar;
+      } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
+                 || CCTK_EQUALS (tensortypealias, "4dd_sym"))
       {
+        /* symmetric 4-tensor */
+        assert (numvars == 10);
+        if (vi == firstvar) {
+          ttype = SCALAR;
+          tcomponent = 0;
+        } else if (vi <= firstvar+3) {
+          ttype = VECTOR;
+          tcomponent = vi - firstvar - 1;
+        } else {
+          ttype = TENSOR;
+          tcomponent = vi - firstvar - 4;
+        }
+      } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
+        /* Weyl scalars, stored as 10 real numbers */
+        assert (numvars == 10);
+        ttype = WEYLSCALARS_REAL;
+        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;
+      case WEYLSCALARS_REAL:
+        assert (tcomponent>=0 && tcomponent<10);
+        break;
+      default:
+        assert (0);
+      }
+      
+      ierr = Util_TableGetInt (table, & tensorparity, "tensorparity");
+      if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+        tensorparity = +1;
+      } else if (ierr<0) {
+        groupname = CCTK_GroupName(gi);
+        assert (groupname);
+        CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Error in tensor parity declaration for group \"%s\"",
+                    groupname);
+        free (groupname);
+      }
+      
+      
+      
+      /* 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;
+      case WEYLSCALARS_REAL: {
         static int const weylparities[10][3] =
           {{+1,+1,+1},
            {-1,-1,-1},
@@ -333,62 +352,56 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
         for (dir=0; dir<3; ++dir) {
           parities[dir] = weylparities[tcomponent][dir];
         }
+        break;
       }
-      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;
+      default:
+        assert (0);
       }
-    }
-    
-    
-    
-    /* Are there negative parities? */
-    needs_checking = 0;
-    for (dir=0; dir<3; ++dir)
-    {
-      check_dir[dir] = do_reflection[2*dir] && parities[dir] < 0;
-      needs_checking |= check_dir[dir];
-    }
-    
-    
-    
-    /* Loop over all points and unfold */
-    if (needs_checking)
-    {
-      for (n=0; n<N_interp_points; ++n)
+      
+      
+      
+      /* Take derivatives into account */
       {
-        int parity = tensorparity;
-        /* 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;
+        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[dir] = do_reflection[2*dir] && parities[dir] < 0;
+        needs_checking |= check_dir[dir];
+      }
+      
+      
+      
+      /* Loop over all points and unfold */
+      if (needs_checking) {
+        for (n=0; n<N_interp_points; ++n) {
+          int parity = tensorparity;
+          /* 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;
         }
-        ((CCTK_REAL *)output_arrays[m])[n] *= parity;
       }
+      
     }
-    
-  }
-  } /* for m */
+    } /* for m */