Support non-symmetric tensors

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

2 files changed, 62 insertions, 8 deletions

diff --git a/src/apply.c b/src/apply.c
index 9375b56..c1736c6 100644
--- a/src/apply.c
+++ b/src/apply.c
@@ -179,7 +179,8 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   int table;
   char tensortypealias[1000];
-  enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL };
+  enum tensortype { UNKNOWN,
+                    SCALAR, VECTOR, SYMTENSOR, TENSOR, WEYLSCALARS_REAL };
   enum tensortype ttype;
   CCTK_INT tensorparity;
   int tcomponent;
@@ -317,6 +318,15 @@ BndReflectVI (cGH const * restrict const cctkGH,
   {
     /* symmetric tensor */
     assert (numvars == 6);
+    ttype = SYMTENSOR;
+    tcomponent = vi - firstvar;
+  } else if (CCTK_EQUALS (tensortypealias, "uu")
+             || CCTK_EQUALS (tensortypealias, "ud")
+             || CCTK_EQUALS (tensortypealias, "du")
+             || CCTK_EQUALS (tensortypealias, "dd"))
+  {
+    /* non-symmetric tensor */
+    assert (numvars == 9);
     ttype = TENSOR;
     tcomponent = vi - firstvar;
   } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
@@ -331,7 +341,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
       ttype = VECTOR;
       tcomponent = vi - firstvar - 1;
     } else {
-      ttype = TENSOR;
+      ttype = SYMTENSOR;
       tcomponent = vi - firstvar - 4;
     }
   } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
@@ -357,9 +367,12 @@ BndReflectVI (cGH const * restrict const cctkGH,
   case VECTOR:
     assert (tcomponent>=0 && tcomponent<3);
     break;
-  case TENSOR:
+  case SYMTENSOR:
     assert (tcomponent>=0 && tcomponent<6);
     break;
+  case TENSOR:
+    assert (tcomponent>=0 && tcomponent<9);
+    break;
   case WEYLSCALARS_REAL:
     assert (tcomponent>=0 && tcomponent<10);
     break;
@@ -415,7 +428,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
           case VECTOR:
             parity *= dir == tcomponent ? -1 : +1;
             break;
-          case TENSOR:
+          case SYMTENSOR:
             switch (tcomponent) {
             case 0: parity *= +1; break;
             case 1: parity *= (dir == 0 || dir == 1) ? -1 : +1; break;
@@ -426,6 +439,20 @@ BndReflectVI (cGH const * restrict const cctkGH,
             default: assert (0);
             }
             break;
+          case TENSOR:
+            switch (tcomponent) {
+            case 0: parity *= +1; break;
+            case 1: parity *= (dir == 0 || dir == 1) ? -1 : +1; break;
+            case 2: parity *= (dir == 0 || dir == 2) ? -1 : +1; break;
+            case 3: parity *= (dir == 1 || dir == 0) ? -1 : +1; break;
+            case 4: parity *= +1; break;
+            case 5: parity *= (dir == 1 || dir == 2) ? -1 : +1; break;
+            case 6: parity *= (dir == 2 || dir == 0) ? -1 : +1; break;
+            case 7: parity *= (dir == 2 || dir == 1) ? -1 : +1; break;
+            case 8: parity *= +1; break;
+            default: assert (0);
+            }
+            break;
           case WEYLSCALARS_REAL: {
             static int const weylparities[10][3] =
               {{+1,+1,+1},
diff --git a/src/interpolate.c b/src/interpolate.c
index 836fc8f..9ec8c17 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -168,7 +168,8 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       
       int table;
       char tensortypealias[1000];
-      enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL };
+      enum tensortype { UNKNOWN,
+                        SCALAR, VECTOR, SYMTENSOR, TENSOR, WEYLSCALARS_REAL };
       enum tensortype ttype;
       CCTK_INT tensorparity;
       int tcomponent;
@@ -255,6 +256,15 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       {
         /* symmetric tensor */
         assert (numvars == 6);
+        ttype = SYMTENSOR;
+        tcomponent = vi - firstvar;
+      } else if (CCTK_EQUALS (tensortypealias, "uu")
+                 || CCTK_EQUALS (tensortypealias, "ud")
+                 || CCTK_EQUALS (tensortypealias, "du")
+                 || CCTK_EQUALS (tensortypealias, "dd"))
+      {
+        /* non-symmetric tensor */
+        assert (numvars == 9);
         ttype = TENSOR;
         tcomponent = vi - firstvar;
       } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
@@ -269,7 +279,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
           ttype = VECTOR;
           tcomponent = vi - firstvar - 1;
         } else {
-          ttype = TENSOR;
+          ttype = SYMTENSOR;
           tcomponent = vi - firstvar - 4;
         }
       } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
@@ -293,9 +303,12 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       case VECTOR:
         assert (tcomponent>=0 && tcomponent<3);
         break;
-      case TENSOR:
+      case SYMTENSOR:
         assert (tcomponent>=0 && tcomponent<6);
         break;
+      case TENSOR:
+        assert (tcomponent>=0 && tcomponent<9);
+        break;
       case WEYLSCALARS_REAL:
         assert (tcomponent>=0 && tcomponent<10);
         break;
@@ -326,7 +339,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       case VECTOR:
         parities[tcomponent] = -1;
         break;
-      case TENSOR:
+      case SYMTENSOR:
         switch (tcomponent) {
         case 0: break;
         case 1: parities[0] = parities[1] = -1; break;
@@ -337,6 +350,20 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
         default: assert (0);
         }
         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: parities[1] = parities[0] = -1; break;
+        case 4: break;
+        case 5: parities[1] = parities[2] = -1; break;
+        case 6: parities[2] = parities[0] = -1; break;
+        case 7: parities[2] = parities[1] = -1; break;
+        case 8: break;
+        default: assert (0);
+        }
+        break;
       case WEYLSCALARS_REAL: {
         static int const weylparities[10][3] =
           {{+1,+1,+1},