handle vector groups of vectors correctly

this patch allows things like

cctk_real fnu[10]
{
  fnux,fnuy,fnuz
} "set of vectors"

to have symmetries applied to them. 


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

2 files changed, 79 insertions, 50 deletions

diff --git a/src/apply.c b/src/apply.c
index 0ef2674..0c09b8e 100644
--- a/src/apply.c
+++ b/src/apply.c
@@ -149,7 +149,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
   int gi;
   cGroup group;
   cGroupDynamicData data;
-  int firstvar, numvars;
+  int firstvar, numvars, vectorlength;
   char * restrict fullname;
   
   void * restrict varptr;
@@ -229,6 +229,9 @@ BndReflectVI (cGH const * restrict const cctkGH,
   assert (firstvar>=0 && firstvar<CCTK_NumVars());
   numvars = CCTK_NumVarsInGroupI (gi);
   assert (numvars>=0);
+  vectorlength = group.vectorlength;
+  assert (vectorlength>=0);
+  assert (vectorlength==1 || group.vectorgroup);
   
   ierr = CCTK_GroupDynamicData (cctkGH, gi, &data);
   assert (!ierr);
@@ -294,64 +297,73 @@ BndReflectVI (cGH const * restrict const cctkGH,
              || CCTK_EQUALS (tensortypealias, "d"))
   {
     /* vector */
-    assert (numvars == 3);
+    const int numcomps = 3;
+    /* special case to handle things like vel[3] */
+    assert (numvars % numcomps == 0 && 
+            (numvars == numcomps * vectorlength || numvars == vectorlength));
     ttype = VECTOR;
-    tcomponent = vi - firstvar;
+    tcomponent = (vi - firstvar) % numcomps;
   } else if (CCTK_EQUALS (tensortypealias, "4u")
              || CCTK_EQUALS (tensortypealias, "4d"))
   {
     /* 4-vector */
-    assert (numvars == 4);
-    if (vi == firstvar) {
+    const int numcomps = 4;
+    assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
+    if ((vi - firstvar) % numcomps == 0) {
       ttype = SCALAR;
       tcomponent = 0;
     } else {
       ttype = VECTOR;
-      tcomponent = vi - firstvar - 1;
+      tcomponent = (vi - firstvar) % numcomps - 1;
     }
   } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
              || CCTK_EQUALS (tensortypealias, "dd_sym"))
   {
     /* symmetric tensor */
-    assert (numvars == 6);
+    const int numcomps = 6;
+    assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
     ttype = SYMTENSOR;
-    tcomponent = vi - firstvar;
+    tcomponent = (vi - firstvar) % numcomps;
   } else if (CCTK_EQUALS (tensortypealias, "uu")
              || CCTK_EQUALS (tensortypealias, "ud")
              || CCTK_EQUALS (tensortypealias, "du")
              || CCTK_EQUALS (tensortypealias, "dd"))
   {
     /* non-symmetric tensor */
-    assert (numvars == 9);
+    const int numcomps = 9;
+    assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
     ttype = TENSOR;
-    tcomponent = vi - firstvar;
+    tcomponent = (vi - firstvar) % numcomps;
   } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
              || CCTK_EQUALS (tensortypealias, "4dd_sym"))
   {
     /* symmetric 4-tensor */
-    assert (numvars == 10);
-    if (vi == firstvar) {
+    const int numcomps = 10;
+    assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
+    if ((vi - firstvar) % numcomps == 0) {
       ttype = SCALAR;
       tcomponent = 0;
-    } else if (vi <= firstvar+3) {
+    } else if ((vi - firstvar) % numcomps <= 3) {
       ttype = VECTOR;
-      tcomponent = vi - firstvar - 1;
+      tcomponent = (vi - firstvar) % numcomps - 1;
     } else {
       ttype = SYMTENSOR;
-      tcomponent = vi - firstvar - 4;
+      tcomponent = (vi - firstvar) % numcomps - 4;
     }
   } else if (CCTK_EQUALS (tensortypealias, "ddd_sym")) {
     /* 3rd rank tensor, symmetric in last 2 indices */
-    assert (numvars == 18);
+    const int numcomps = 18;
+    assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
     ttype = SYMTENSOR3;
-    tcomponent = vi - firstvar;
+    tcomponent = (vi - firstvar) % numcomps;
   } 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
        PsiKadelia.  */
-    assert (numvars == 10);
+    const int numcomps = 10;
+    assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
     ttype = WEYLSCALARS_REAL;
-    tcomponent = vi - firstvar;
+    tcomponent = (vi - firstvar) % numcomps;
   } else if (CCTK_EQUALS (tensortypealias, "ManualCartesian")) {
       /* Reflection symmetries specified by hand */
       ttype = MANUALCARTESIAN;
diff --git a/src/interpolate.c b/src/interpolate.c
index 8362db4..571c048 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -76,15 +76,16 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
     assert (! do_reflection[2*dir+1]);
     
     if (do_reflection[2*dir]) {
-      new_interp_coords[dir]
+      CCTK_REAL * restrict coords_in_dir
         = malloc (N_interp_points * sizeof (CCTK_REAL));
-      assert (N_interp_points == 0 || new_interp_coords[dir]);
+      assert (N_interp_points == 0 || coords_in_dir);
       
       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;
+        coords_in_dir[n] = newpos;
       }
+      new_interp_coords[dir] = coords_in_dir;
     } else {
       new_interp_coords[dir] = interp_coords[dir];
     }
@@ -96,9 +97,9 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
   iret = SymmetryInterpolateFaces
     (cctkGH, 
      N_dims, local_interp_handle, param_table_handle, coord_system_handle,
-     N_interp_points, interp_coords_type, new_interp_coords,
+     N_interp_points, interp_coords_type, (CCTK_POINTER_TO_CONST)new_interp_coords,
      N_input_arrays, input_array_indices,
-     N_output_arrays, output_array_types, output_arrays,
+     N_output_arrays, output_array_types, (CCTK_POINTER_TO_CONST)output_arrays,
      newfaces);
   
   
@@ -106,7 +107,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
   /* Free coordinates */
   for (dir=0; dir<3; ++dir) {
     if (do_reflection[2*dir]) {
-      free (new_interp_coords[dir]);
+      free ((void*)new_interp_coords[dir]);
     }
   }
   
@@ -157,7 +158,8 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
     if (output_array_indices[m]!=-1) {
       
       int vi, gi;
-      int numvars, firstvar;
+      int numvars, firstvar, vectorlength;
+      cGroup group;
       char * groupname;
       
       int table;
@@ -172,17 +174,23 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
       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());
+
+      ierr = CCTK_GroupData (gi, &group);
+      assert (!ierr);
+      
       numvars = CCTK_NumVarsInGroupI(gi);
       assert (numvars>0);
       firstvar = CCTK_FirstVarIndexI(gi);
       assert (firstvar>=0);
+      vectorlength = group.vectorlength;
+      assert (vectorlength>=0);
+      assert (vectorlength==1 || group.vectorgroup);
       table = CCTK_GroupTagsTableI(gi);
       assert (table>=0);
   
@@ -212,7 +220,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
       tcomponent = 0;
       if (CCTK_EQUALS (tensortypealias, "scalar")) {
         /* scalar */
-        if (numvars != 1) {
+        if (numvars != vectorlength) {
           groupname = CCTK_GroupName(gi);
           assert (groupname);
           CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -224,72 +232,81 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH,
         tcomponent = 0;
       } else if (CCTK_EQUALS (tensortypealias, "4scalar")) {
         /* 4-scalar */
-        assert (numvars == 1);
+        assert (numvars == vectorlength);
         ttype = SCALAR;
         tcomponent = 0;
       } else if (CCTK_EQUALS (tensortypealias, "u")
                  || CCTK_EQUALS (tensortypealias, "d"))
       {
         /* vector */
-        assert (numvars == 3);
+        const int numcomps = 3;
+        /* special case to handle things like vel[3] */
+        assert (numvars % numcomps == 0 && 
+                (numvars == numcomps * vectorlength || numvars == vectorlength));
         ttype = VECTOR;
-        tcomponent = vi - firstvar;
+        tcomponent = (vi - firstvar) % numcomps;
       } else if (CCTK_EQUALS (tensortypealias, "4u")
                  || CCTK_EQUALS (tensortypealias, "4d"))
       {
         /* 4-vector */
-        assert (numvars == 4);
-        if (vi == firstvar) {
+        const int numcomps = 4;
+        assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
+        if ((vi - firstvar) % numcomps == 0) {
           ttype = SCALAR;
           tcomponent = 0;
         } else {
           ttype = VECTOR;
-          tcomponent = vi - firstvar - 1;
+          tcomponent = (vi - firstvar) % numcomps - 1;
         }
       } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
                  || CCTK_EQUALS (tensortypealias, "dd_sym"))
       {
         /* symmetric tensor */
-        assert (numvars == 6);
+        const int numcomps = 6;
+        assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
         ttype = SYMTENSOR;
-        tcomponent = vi - firstvar;
+        tcomponent = (vi - firstvar) % numcomps;
       } else if (CCTK_EQUALS (tensortypealias, "uu")
                  || CCTK_EQUALS (tensortypealias, "ud")
                  || CCTK_EQUALS (tensortypealias, "du")
                  || CCTK_EQUALS (tensortypealias, "dd"))
       {
         /* non-symmetric tensor */
-        assert (numvars == 9);
+        const int numcomps = 9;
+        assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
         ttype = TENSOR;
-        tcomponent = vi - firstvar;
+        tcomponent = (vi - firstvar) % numcomps;
       } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
                  || CCTK_EQUALS (tensortypealias, "4dd_sym"))
       {
         /* symmetric 4-tensor */
-        assert (numvars == 10);
-        if (vi == firstvar) {
+        const int numcomps = 10;
+        assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
+        if ((vi - firstvar) % numcomps == 0) {
           ttype = SCALAR;
           tcomponent = 0;
-        } else if (vi <= firstvar+3) {
+        } else if ((vi - firstvar) % numcomps <= 3) {
           ttype = VECTOR;
-          tcomponent = vi - firstvar - 1;
+          tcomponent = (vi - firstvar) % numcomps - 1;
         } else {
           ttype = SYMTENSOR;
-          tcomponent = vi - firstvar - 4;
+          tcomponent = (vi - firstvar) % numcomps - 4;
         }
       } else if (CCTK_EQUALS (tensortypealias, "ddd_sym")) {
         /* 3rd rank tensor, symmetric in last 2 indices */
-        assert (numvars == 18);
+        const int numcomps = 18;
+        assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
         ttype = SYMTENSOR3;
-        tcomponent = vi - firstvar;
+        tcomponent = (vi - firstvar) % numcomps;
       } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
         /* Weyl scalars, stored as 10 real numbers */
-        assert (numvars == 10);
+        const int numcomps = 10;
+        assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength);
         ttype = WEYLSCALARS_REAL;
-        tcomponent = vi - firstvar;
+        tcomponent = (vi - firstvar) % numcomps;
       } else if (CCTK_EQUALS (tensortypealias, "ManualCartesian")) {
-      /* Reflection symmetries specified by hand */
-      ttype = MANUALCARTESIAN;
+        /* Reflection symmetries specified by hand */
+        ttype = MANUALCARTESIAN;
       tcomponent = vi - firstvar;
       } else {
         groupname = CCTK_GroupName(gi);