Apply symmetry boundary conditions on all faces (not only on lower ones).

Now you can have negative octant, quadrant, etc.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/CartGrid3D/trunk@193 c78560ca-4b45-4335-b268-5f3340f3cb52

2 files changed, 155 insertions, 150 deletions

diff --git a/src/CartGrid3D.c b/src/CartGrid3D.c
index 3592adf..87b76ee 100644
--- a/src/CartGrid3D.c
+++ b/src/CartGrid3D.c
@@ -104,13 +104,16 @@ void CartGrid3D(CCTK_ARGUMENTS)
     /* Calculate physical indices, using symmetries and periodicity */
     for (i = 0; i < 3; i++)
     {
-      loweri[i] = domainsym[2*i] ? cctk_nghostzones[i] : 0;
+      loweri[i] = 0;
       upperi[i] = cctk_gsh[i] - 1;
 
-      if (do_periodic[i])
+      if (domainsym[2*i+0] || do_periodic[i])
       {
-        loweri[i] = cctk_nghostzones[i];
-        upperi[i] = cctk_gsh[i] - 1 - cctk_nghostzones[i];
+        loweri[i] += cctk_nghostzones[i];
+      }
+      if (domainsym[2*i+1] || do_periodic[i])
+      {
+        upperi[i] -= cctk_nghostzones[i];
       }
     }
 
@@ -162,19 +165,31 @@ void CartGrid3D(CCTK_ARGUMENTS)
       /* Grid spacing on coarsest grid */
       for (i = 0; i < 3; i++)
       {
-        if (domainsym[2*i])
+        if (domainsym[2*i+0])
         {
           if (cntstag[i])
           {
-            *coarse_delta[i] = max1[i] / (cctk_gsh[i] - cctk_nghostzones[i]-0.5);
+            *coarse_delta[i] = max1[i] / (cctk_gsh[i]-cctk_nghostzones[i]-0.5);
             origin[i] = - (cctk_nghostzones[i]-0.5) * *coarse_delta[i];
           }
           else
           {
-            *coarse_delta[i] = max1[i] / (cctk_gsh[i] - cctk_nghostzones[i]-1);
+            *coarse_delta[i] = max1[i] / (cctk_gsh[i]-cctk_nghostzones[i]-1);
             origin[i] = -  cctk_nghostzones[i]      * *coarse_delta[i];
           }
         }
+        else if (domainsym[2*i+1])
+        {
+          if (cntstag[i])
+          {
+            *coarse_delta[i] = fabs(min1[i]) / (cctk_gsh[i]-cctk_nghostzones[i]-0.5);
+          }
+          else
+          {
+            *coarse_delta[i] = fabs(min1[i]) / (cctk_gsh[i]-cctk_nghostzones[i]-1);
+          }
+          origin[i] = min1[i];
+        }
         else
         {
           if (cntstag[i])
@@ -215,15 +230,19 @@ void CartGrid3D(CCTK_ARGUMENTS)
         this_delta[i] = *coarse_delta[i];
 
         /* Set minimum values of coordinates */
-        if (domainsym[2*i])
+        if (domainsym[2*i+0])
+        {
+          origin[i] = -(cctk_nghostzones[i] - cntstag[i]*0.5);
+        }
+        else if (domainsym[2*i+1])
         {
-          origin[i] = (- cctk_nghostzones[i] + cntstag[i]*0.5) * this_delta[i];
+          origin[i] = -(cctk_gsh[i]-1 - cctk_nghostzones[i] + cntstag[i]*0.5);
         }
         else
         {
-          origin[i] = - 0.5 * (cctk_gsh[i]-1 - cntstag[i] * cctk_gsh[i]%2)
-            * this_delta[i];
+          origin[i] = -0.5 * (cctk_gsh[i]-1 - cntstag[i] * cctk_gsh[i]%2);
         }
+        origin[i] *= this_delta[i];
       }
     }
 
diff --git a/src/Symmetry.c b/src/Symmetry.c
index 8f3f832..0a2cb58 100644
--- a/src/Symmetry.c
+++ b/src/Symmetry.c
@@ -2,10 +2,10 @@
    @file      SymmetryWrappers.c
    @date      April 2000
    @author    Gerd Lanfermann
-   @desc 
+   @desc
               Routines to apply the 1/2/3D Symmetries for
               all symmetry domains (octant/bitant/quadrant).
-   @enddesc 
+   @enddesc
    @history
    @hdate     Sat 02 Nov 2002
    @hauthor   Thomas Radke
@@ -76,8 +76,8 @@ void DecodeSymParameters3D (int sym[6]);
                -1 if invalid group index was given
    @endreturndesc
 @@*/
-int CartSymGI (const cGH *GH, int gindex) 
-{ 
+int CartSymGI (const cGH *GH, int gindex)
+{
   int numvars, first_vindex, retval;
 
 
@@ -184,7 +184,7 @@ void CCTK_FCALL CCTK_FNAME (CartSymGN)
                -1 if invalid variable index was given
    @endreturndesc
 @@*/
-int CartSymVI (const cGH *GH, int vindex) 
+int CartSymVI (const cGH *GH, int vindex)
 {
   int retval, gindex;
 
@@ -251,7 +251,7 @@ int CartSymVN (const cGH *GH, const char *vname)
     CCTK_VWarn (1,__LINE__,__FILE__,CCTK_THORNSTRING,
                 "Invalid variable name '%s' in CartSymVN", vname);
     retval = -1;
-  } 
+  }
 
   return (retval);
 }
@@ -288,55 +288,71 @@ void CCTK_FCALL CCTK_FNAME (CartSymVN)
    @vio        in
    @endvar
 @@*/
-#define SYMMETRY_BOUNDARY(cctk_type)                                          \
+#define APPLY_LOWER(dir, ii, jj, kk, jjend, kkend, iii, jjj, kkk)             \
 {                                                                             \
-  cctk_type *_var = (cctk_type *) GH->data[vindex][0];                        \
-                                                                              \
-                                                                              \
-  /* apply symmetry to the lower x face */                                    \
-  if (doSym[0] == GFSYM_REFLECTION)                                           \
+  for (kk = 0; kk < kkend; kk++)                                              \
   {                                                                           \
-    for (k = 0; k < lssh[2]; k++)                                             \
+    for (jj = 0; jj < jjend; jj++)                                            \
     {                                                                         \
-      for (j = 0; j < lssh[1]; j++)                                           \
+      for (ii = 0; ii < gdata.nghostzones[dir/2]; ii++)                       \
       {                                                                       \
-        for (i = 0; i < group_dynamic_data.nghostzones[0]; i++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][0] *      \
-            _var[INDEX_3D (lsh, offset[0]-i, j, k)] NUMBER_PART;              \
-        }                                                                     \
+        _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][dir] *      \
+             _var[INDEX_3D (lsh, iii, jjj, kkk)] NUMBER_PART;                 \
       }                                                                       \
     }                                                                         \
   }                                                                           \
-  else if (group_static_data.dim > 1)                                         \
+}
+
+#define APPLY_UPPER(dir, ii, jj, kk, jjend, kkend, iii, jjj, kkk)             \
+{                                                                             \
+  for (kk = 0; kk < kkend; kk++)                                              \
   {                                                                           \
-    if (doSym[0] == GFSYM_ROTATION_Z)                                         \
+    for (jj = 0; jj < jjend; jj++)                                            \
     {                                                                         \
-      for (k = 0; k < lssh[2]; k++)                                           \
+      for (ii = lsh[dir/2]-gdata.nghostzones[dir/2]; ii < lsh[dir/2]; ii++)   \
       {                                                                       \
-        for (j = 0; j < lssh[1]; j++)                                         \
-        {                                                                     \
-          for (i = 0; i < group_dynamic_data.nghostzones[0]; i++)             \
-          {                                                                   \
-            _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][0] *    \
-              _var[INDEX_3D (lsh, offset[0]-i, lssh[1]-j-1, k)] NUMBER_PART;  \
-          }                                                                   \
-        }                                                                     \
+        _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][dir] *      \
+             _var[INDEX_3D (lsh, iii, jjj, kkk)] NUMBER_PART;                 \
       }                                                                       \
     }                                                                         \
-    else if (group_static_data.dim > 2 && doSym[0] == GFSYM_ROTATION_Y)       \
+  }                                                                           \
+}
+
+
+#define SYMMETRY_BOUNDARY(cctk_type)                                          \
+{                                                                             \
+  cctk_type *_var = GH->data[vindex][0];                                      \
+                                                                              \
+                                                                              \
+  /* apply symmetry to the x faces */                                         \
+  if (doSym[0] == GFSYM_REFLECTION)                                           \
+  {                                                                           \
+    APPLY_LOWER (0, i, j, k, lssh[1], lssh[2], offset[0]-i, j, k);            \
+  }                                                                           \
+  if (doSym[1] == GFSYM_REFLECTION)                                           \
+  {                                                                           \
+    APPLY_UPPER (1, i, j, k, lssh[1], lssh[2], offset[1]-i, j, k);            \
+  }                                                                           \
+  if (group_static_data.dim > 1)                                              \
+  {                                                                           \
+    if (doSym[0] == GFSYM_ROTATION_Z)                                         \
     {                                                                         \
-      for (k = 0; k < lssh[2]; k++)                                           \
-      {                                                                       \
-        for (j = 0; j < lssh[1]; j++)                                         \
-        {                                                                     \
-          for (i = 0; i < group_dynamic_data.nghostzones[0]; i++)             \
-          {                                                                   \
-            _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][0] *    \
-              _var[INDEX_3D (lsh, offset[0]-i, j, lssh[2]-k-1)] NUMBER_PART;  \
-          }                                                                   \
-        }                                                                     \
-      }                                                                       \
+      APPLY_LOWER (0, i, j, k, lssh[1], lssh[2], offset[0]-i, lssh[1]-j-1, k);\
+    }                                                                         \
+    if (doSym[1] == GFSYM_ROTATION_Z)                                         \
+    {                                                                         \
+      APPLY_UPPER (1, i, j, k, lssh[1], lssh[2], offset[1]-i, lssh[1]-j-1, k);\
+    }                                                                         \
+  }                                                                           \
+  if (group_static_data.dim > 2)                                              \
+  {                                                                           \
+    if (doSym[0] == GFSYM_ROTATION_Y)                                         \
+    {                                                                         \
+      APPLY_LOWER (0, i, j, k, lssh[1], lssh[2], offset[0]-i, j, lssh[2]-k-1);\
+    }                                                                         \
+    if (doSym[1] == GFSYM_ROTATION_Y)                                         \
+    {                                                                         \
+      APPLY_UPPER (1, i, j, k, lssh[1], lssh[2], offset[1]-i, j, lssh[2]-k-1);\
     }                                                                         \
   }                                                                           \
                                                                               \
@@ -345,47 +361,32 @@ void CCTK_FCALL CCTK_FNAME (CartSymVN)
     break;                                                                    \
   }                                                                           \
                                                                               \
-  /* apply symmetry to the lower y face */                                    \
+  /* apply symmetry to the y faces */                                         \
   if (doSym[2] == GFSYM_REFLECTION)                                           \
   {                                                                           \
-    for (i = 0; i < lssh[0]; i++)                                             \
-    {                                                                         \
-      for (k = 0; k < lssh[2]; k++)                                           \
-      {                                                                       \
-        for (j = 0; j < group_dynamic_data.nghostzones[1]; j++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][2] *      \
-            _var[INDEX_3D (lsh, i, offset[1]-j, k)] NUMBER_PART;              \
-        }                                                                     \
-      }                                                                       \
-    }                                                                         \
+    APPLY_LOWER (2, j, k, i, lssh[2], lssh[0], i, offset[2]-j, k);            \
   }                                                                           \
-  else if (doSym[2] == GFSYM_ROTATION_Z)                                      \
+  if (doSym[3] == GFSYM_REFLECTION)                                           \
   {                                                                           \
-    for (i = 0; i < lssh[0]; i++)                                             \
-    {                                                                         \
-      for (k = 0; k < lssh[2]; k++)                                           \
-      {                                                                       \
-        for (j = 0; j < group_dynamic_data.nghostzones[1]; j++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][2] *      \
-            _var[INDEX_3D (lsh, lssh[0]-i-1, offset[1]-j, k)] NUMBER_PART;    \
-        }                                                                     \
-      }                                                                       \
-    }                                                                         \
+    APPLY_UPPER (3, j, k, i, lssh[2], lssh[0], i, offset[3]-j, k);            \
+  }                                                                           \
+  if (doSym[2] == GFSYM_ROTATION_Z)                                           \
+  {                                                                           \
+    APPLY_LOWER (2, j, k, i, lssh[2], lssh[0], lssh[0]-i-1, offset[2]-j, k);  \
+  }                                                                           \
+  if (doSym[3] == GFSYM_ROTATION_Z)                                           \
+  {                                                                           \
+    APPLY_UPPER (3, j, k, i, lssh[2], lssh[0], lssh[0]-i-1, offset[3]-j, k);  \
   }                                                                           \
-  else if (group_static_data.dim > 2 && doSym[2] == GFSYM_ROTATION_X)         \
+  if (group_static_data.dim > 2)                                              \
   {                                                                           \
-    for (i = 0; i < lssh[0]; i++)                                             \
+    if (doSym[2] == GFSYM_ROTATION_X)                                         \
     {                                                                         \
-      for (k = 0; k < lssh[2]; k++)                                           \
-      {                                                                       \
-        for (j = 0; j < group_dynamic_data.nghostzones[1]; j++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][2] *      \
-            _var[INDEX_3D (lsh, i, offset[1]-j, lssh[2]-k-1)] NUMBER_PART;    \
-        }                                                                     \
-      }                                                                       \
+      APPLY_LOWER (2, j, k, i, lssh[2], lssh[0], i, offset[2]-j, lssh[2]-k-1);\
+    }                                                                         \
+    if (doSym[3] == GFSYM_ROTATION_X)                                         \
+    {                                                                         \
+      APPLY_UPPER (3, j, k, i, lssh[2], lssh[0], i, offset[3]-j, lssh[2]-k-1);\
     }                                                                         \
   }                                                                           \
                                                                               \
@@ -394,48 +395,30 @@ void CCTK_FCALL CCTK_FNAME (CartSymVN)
     break;                                                                    \
   }                                                                           \
                                                                               \
-  /* apply symmetry to the lower z face */                                    \
+  /* apply symmetry to the z faces */                                         \
   if (doSym[4] == GFSYM_REFLECTION)                                           \
   {                                                                           \
-    for (i = 0; i < lssh[0]; i++)                                             \
-    {                                                                         \
-      for (j = 0; j < lssh[1]; j++)                                           \
-      {                                                                       \
-        for (k = 0; k < group_dynamic_data.nghostzones[2]; k++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][4] *      \
-            _var[INDEX_3D (lsh, i, j, offset[2]-k)] NUMBER_PART;              \
-        }                                                                     \
-      }                                                                       \
-    }                                                                         \
+    APPLY_LOWER (4, k, j, i, lssh[1], lssh[0], i, j, offset[4]-k);            \
   }                                                                           \
-  else if (doSym[4] == GFSYM_ROTATION_X)                                      \
+  if (doSym[5] == GFSYM_REFLECTION)                                           \
   {                                                                           \
-    for (i = 0; i < lssh[0]; i++)                                             \
-    {                                                                         \
-      for (j = 0; j < lssh[1]; j++)                                           \
-      {                                                                       \
-        for (k = 0; k < group_dynamic_data.nghostzones[2]; k++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][4] *      \
-            _var[INDEX_3D (lsh, i, lssh[1]-j-1, offset[2]-k)] NUMBER_PART;    \
-        }                                                                     \
-      }                                                                       \
-    }                                                                         \
+    APPLY_UPPER (5, k, j, i, lssh[1], lssh[0], i, j, offset[5]-k);            \
   }                                                                           \
-  else if (doSym[4] == GFSYM_ROTATION_Y)                                      \
+  if (doSym[4] == GFSYM_ROTATION_X)                                           \
   {                                                                           \
-    for (i = 0; i < lssh[0]; i++)                                             \
-    {                                                                         \
-      for (j = 0; j < lssh[1]; j++)                                           \
-      {                                                                       \
-        for (k = 0; k < group_dynamic_data.nghostzones[2]; k++)               \
-        {                                                                     \
-          _var[INDEX_3D (lsh, i, j, k)] NUMBER_PART = GFSym[vindex][4] *      \
-            _var[INDEX_3D (lsh, lssh[0]-i-1, j, offset[2]-k)] NUMBER_PART;    \
-        }                                                                     \
-      }                                                                       \
-    }                                                                         \
+    APPLY_LOWER (4, k, j, i, lssh[1], lssh[0], i, lssh[1]-j-1, offset[4]-k);  \
+  }                                                                           \
+  if (doSym[5] == GFSYM_ROTATION_X)                                           \
+  {                                                                           \
+    APPLY_UPPER (5, k, j, i, lssh[1], lssh[0], i, lssh[1]-j-1, offset[5]-k);  \
+  }                                                                           \
+  if (doSym[4] == GFSYM_ROTATION_Y)                                           \
+  {                                                                           \
+    APPLY_LOWER (4, k, j, i, lssh[1], lssh[0], lssh[0]-i-1, j, offset[4]-k);  \
+  }                                                                           \
+  if (doSym[5] == GFSYM_ROTATION_Y)                                           \
+  {                                                                           \
+    APPLY_UPPER (5, k, j, i, lssh[1], lssh[0], lssh[0]-i-1, j, offset[5]-k);  \
   }                                                                           \
 }
 
@@ -493,24 +476,31 @@ static int ApplySymmetry (const cGH *GH, int gindex, int first_vindex,
 {
   int i, j, k, dim, vindex, retval;
   int **GFSym;
-  int domainsym[MAX_FACE], doSym[MAX_FACE];
-  int dstag[MAX_DIM], lsh[MAX_DIM], lssh[MAX_DIM],
-      cntstag[MAX_DIM], offset[MAX_DIM];
+  int domainsym[2*MAX_DIM], doSym[2*MAX_DIM], offset[2*MAX_DIM];
+  int dstag[MAX_DIM], lsh[MAX_DIM], lssh[MAX_DIM], cntstag[MAX_DIM];
   cGroup group_static_data;
-  cGroupDynamicData group_dynamic_data;
+  cGroupDynamicData gdata;
   DECLARE_CCTK_PARAMETERS
 
 
-  /* get out if we are sure no symmetries should be applied */
-  /* FIXME: There has to be a better way early bailout! */
-  if (CCTK_Equals (domain, "full"))
+  DecodeSymParameters3D (domainsym);
+
+  /* check if any symmetries are to be applied */
+  for (i = 0; i < 2*MAX_DIM; i++)
+  {
+    if (domainsym[i])
+    {
+      break;
+    }
+  }
+  if (i == 2*MAX_DIM)
   {
     return (0);
   }
 
   /* get the group's static and dynamic data structure */
   CCTK_GroupData (gindex, &group_static_data);
-  CCTK_GroupDynamicData (GH, gindex, &group_dynamic_data);
+  CCTK_GroupDynamicData (GH, gindex, &gdata);
   if (group_static_data.dim <= 0 || group_static_data.dim > MAX_DIM)
   {
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -526,52 +516,48 @@ static int ApplySymmetry (const cGH *GH, int gindex, int first_vindex,
   /* get the directional staggering of the group */
   CCTK_StaggerDirArray (dstag, group_static_data.dim,
                         group_static_data.stagtype);
-  
+
   /* initialize array for variables with less dimensions than MAX_DIM
      so that we can use the INDEX_3D macro later on */
   for (i = 0; i < MAX_DIM; i++)
   {
     if (i < group_static_data.dim)
     {
-      lsh[i] = group_dynamic_data.lsh[i];
+      lsh[i] = gdata.lsh[i];
       lssh[i] = GH->cctk_lssh[CCTK_LSSH_IDX (dstag[i], i)];
     }
     else
     {
       lsh[i] = lssh[i] = 1;
     }
-    offset[i] = 2*group_dynamic_data.nghostzones[i] - cntstag[i];
+    offset[2*i+0] = 2*gdata.nghostzones[i] - cntstag[i];
+    offset[2*i+1] = 2*(lsh[i]-1) - offset[2*i+0];
   }
 
-  DecodeSymParameters3D (domainsym);
-
   GFSym = ((SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry"))->GFSym;
 
-  /* Apply Symmetries to lower sides [0,2,4] if:
-     + if the Symmetry is activated (== NOT NOSYM)
-     + if the Symmetry is set (== NOT UNSET)
-     + if the length in the direction is more than 1 grid point
-     + if the processor has a lower physical boundary.
+  /* Apply Symmetries if:
+     + the Symmetry is activated (== NOT NOSYM)
+     + the Symmetry is set (== NOT UNSET)
+     + the length in the direction is more than 1 grid point
+     + the processor has a lower/upper physical boundary.
      Whether a grid allows a symmetry along a direction (e.g. octant=all)
      is part if the Symmetry Setup process.
-
-     No Symmetries for "upper" sides : [1,3,5]
   */
   retval = 0;
-  for (vindex = first_vindex; vindex < first_vindex + numvars; vindex++) 
+  for (vindex = first_vindex; vindex < first_vindex + numvars; vindex++)
   {
-    for (dim = 0; dim < group_static_data.dim; dim++)
+    for (dim = 0; dim < 2*group_static_data.dim; dim++)
     {
-      if (GFSym[vindex][2*dim] == GFSYM_UNSET)
+      if (GFSym[vindex][dim] == GFSYM_UNSET)
       {
         CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "Symmetries unspecified for '%s'", CCTK_VarName (vindex));
       }
 
-      doSym[2*dim] = (GFSym[vindex][2*dim] != GFSYM_NOSYM &&
-                      GFSym[vindex][2*dim] != GFSYM_UNSET &&
-                      lssh[dim] > 1 && GH->cctk_bbox[2*dim]) ?
-                      domainsym[2*dim+0] : 0;
+      doSym[dim] = (GFSym[vindex][dim] != GFSYM_NOSYM &&
+                    GFSym[vindex][dim] != GFSYM_UNSET &&
+                    lssh[dim/2] > 1 && GH->cctk_bbox[dim]) ? domainsym[dim] : 0;
     }
 
 #define NUMBER_PART