A symmetry thorn with a 180 degree rotating symmetry. This will

replate PlanarPointSymmetry, which has a stupid name.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/RotatingSymmetry180/trunk@2 20f44201-0f4f-0410-9130-e5fc2714a787

3 files changed, 418 insertions, 0 deletions

diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..2b5b357
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn RotatingSymmetry180
+# $Header$
+
+# Source files in this directory
+SRCS = rotatingsymmetry180.c registersymmetry.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/registersymmetry.c b/src/registersymmetry.c
new file mode 100644
index 0000000..2916efd
--- /dev/null
+++ b/src/registersymmetry.c
@@ -0,0 +1,39 @@
+/* $Header$ */
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void Rot180_RegisterSymmetry (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int f;
+  CCTK_INT handle;
+  CCTK_INT faces[6];
+  CCTK_INT width[6];
+  CCTK_INT ierr;
+  
+  if (rotating_symmetry_180) {
+    
+    for (f=0; f<6; ++f) {
+      faces[f] = 0;
+      width[f] = 0;
+    }
+    
+    faces[0] = 1;
+    width[0] = cctk_nghostzones[0];
+    
+    handle = SymmetryRegister ("rotating_symmetry_180");
+    if (handle < 0) {
+      CCTK_WARN (0, "Could not register symmetry boundary condition");
+    }
+    
+    ierr = SymmetryRegisterGrid (cctkGH, handle, faces, width);
+    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");
+    }
+    
+  }
+}
diff --git a/src/rotatingsymmetry180.c b/src/rotatingsymmetry180.c
new file mode 100644
index 0000000..ebff1db
--- /dev/null
+++ b/src/rotatingsymmetry180.c
@@ -0,0 +1,371 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "Slab.h"
+
+
+
+int BndRot180VI (cGH const * restrict const cctkGH,
+                 int const * restrict const stencil,
+                 int const vi)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int gi;
+  cGroup group;
+  cGroupDynamicData data;
+  char * restrict fullname;
+  void * restrict varptr;
+  
+  int parities[3];
+  
+  int global_bbox[6];
+  int global_lbnd[3], global_ubnd[3];
+  int fake_bbox[6];
+  
+  CCTK_REAL x0[3], dx[3];
+  CCTK_REAL origin[3], dorigin[3];
+  int avoid_origin[3], iorigin[3];
+  int offset[3];                /* offset 0..1 due to avoid_origin */
+  
+  struct xferinfo * restrict xferinfo;
+  
+  int dirs[2];
+  int dir;                      /* direction of the symmetry face */
+  int otherdir;                 /* the other direction of the rotation */
+  int q;
+  int d;
+  int ierr;
+  
+  /* Check arguments */
+  assert (cctkGH);
+  assert (stencil);
+  assert (vi>=0 && vi<CCTK_NumVars());
+  
+  if (! rotating_symmetry_180) return;
+  
+  if (verbose) {
+    fullname = CCTK_FullName(vi);
+    assert (fullname);
+    CCTK_VInfo (CCTK_THORNSTRING,
+                "Applying 180 degree rotating symmetry boundary conditions to \"%s\"",
+                fullname);
+    free (fullname);
+  }
+  
+  for (d=0; d<group.dim; ++d) {
+    x0[d] = CCTK_ORIGIN_SPACE(d);
+    dx[d] = CCTK_DELTA_SPACE(d);
+  }
+  
+  /* Get and check group info */
+  gi = CCTK_GroupIndexFromVarI (vi);
+  assert (gi>=0 && gi<CCTK_NumGroups());
+  
+  ierr = CCTK_GroupData (gi, &group);
+  assert (!ierr);
+  assert (group.grouptype == CCTK_GF);
+  assert (group.disttype == CCTK_DISTRIB_DEFAULT);
+  assert (group.stagtype == 0);
+  
+  ierr = CCTK_GroupDynamicData (cctkGH, gi, &data);
+  assert (!ierr);
+  
+  varptr = CCTK_VarDataPtrI (cctkGH, 0, vi);
+  assert (varptr);
+  
+  /* find parity */
+  {
+    char tensortypealias[100];
+    int firstvar, numvars;
+    int table;
+    int index;
+    
+    numvars = CCTK_NumVarsInGroupI(gi);
+    assert (numvars>0);
+    firstvar = CCTK_FirstVarIndexI(gi);
+    assert (firstvar>=0);
+    index = vi - firstvar;
+    assert (index>=0 && index<numvars);
+    table = CCTK_GroupTagsTableI(gi);
+    assert (table>=0);
+    
+    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);
+    }
+    
+    if (CCTK_EQUALS (tensortypealias, "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);
+      }
+      parities[0] = parities[1] = parities[2] = +1;
+    } else if (CCTK_EQUALS (tensortypealias, "u")) {
+      assert (numvars == 3);
+      parities[0] = parities[1] = parities[2] = +1;
+      parities[index] = -1;
+    } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
+      assert (numvars == 6);
+      parities[0] = parities[1] = parities[2] = +1;
+      switch (index) {
+      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);
+      }
+    } 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);
+    }
+    
+    assert (abs(parities[0])==1 && abs(parities[1])==1
+            && abs(parities[2])==1);
+  }
+  
+  {
+    int min_handle, max_handle;
+    CCTK_REAL local[6], global[6];
+    min_handle = CCTK_ReductionArrayHandle ("minimum");
+    if (min_handle<0) CCTK_WARN (0, "Could not obtain reduction handle");
+    max_handle = CCTK_ReductionArrayHandle ("maximum");
+    if (max_handle<0) CCTK_WARN (0, "Could not obtain reduction handle");
+    
+    for (d=0; d<6; ++d) local[d] = cctkGH->cctk_bbox[d];
+    ierr = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, max_handle, local, global, 6, CCTK_VARIABLE_REAL);
+    for (d=0; d<6; ++d) global_bbox[d] = (int)global[d];
+    
+    for (d=0; d<3; ++d) local[d] = cctkGH->cctk_lbnd[d];
+    ierr = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, min_handle, local, global, 3, CCTK_VARIABLE_REAL);
+    for (d=0; d<3; ++d) global_lbnd[d] = (int)global[d];
+    
+    for (d=0; d<3; ++d) local[d] = cctkGH->cctk_ubnd[d];
+    ierr = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, max_handle, local, global, 3, CCTK_VARIABLE_REAL);
+    for (d=0; d<3; ++d) global_ubnd[d] = (int)global[d];
+    
+    for (d=0; d<3; ++d) {
+      fake_bbox[2*d  ] = data.lbnd[d] == global_lbnd[d];
+      fake_bbox[2*d+1] = data.ubnd[d] == global_ubnd[d];
+    }
+  }
+  
+  /* directions */
+  dirs[0] = 0;
+  dirs[1] = 1;
+  
+  /* Find grid point that corresponds to the origin */
+  for (q=0; q<2; ++q) {
+    d = dirs[q];
+    /* x0 + dx * origin == 0 */
+    origin[d] = - x0[d] / dx[d];
+    dorigin[d] = origin[d] - floor(origin[d]);
+    if (fabs(dorigin[d]) < 1.0e-6 || fabs(dorigin[d] - 1.0) < 1.0e-6) {
+      avoid_origin[d] = 0;
+    } else if (fabs(dorigin[d] - 0.5) < 1.0e-6) {
+      avoid_origin[d] = 1;
+    } else {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The coordinate origin in the %c-direction falls neither onto a grid point nor into the middle between two grid points.",
+                  "xyz"[d]);
+    }
+    iorigin[d] = floor(origin[d] + (avoid_origin[d] ? 0.5 : 0.0) + 0.5);
+    offset[d] = avoid_origin[d] ? 0 : 1;
+  }
+  
+  /* x direction */
+  dir = 0;
+  otherdir = 1;
+  
+  assert (stencil[dir] >= 0);
+  
+  if (global_bbox[2*dir]) {
+    
+    if (2*iorigin[otherdir] + offset[otherdir] != cctk_gsh[otherdir]) {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The coordinate origin in the %c-direction is not in the centre of the domain.  The boundary condition cannot be applied.",
+                  "xyz"[otherdir]);
+    }
+    
+    if (data.gsh[dir] < 2*stencil[dir] + offset[dir]) {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The group \"%s\" has in the %c-direction only %d grid points.  "
+                  "This is not large enough for a 180 degree rotating symmetry boundary that is %d grid points wide.  "
+                  "The group needs to have at least %d grid points in that direction.",
+                  CCTK_GroupNameFromVarI(vi), "xyz"[dir], data.gsh[dir],
+                  stencil[dir],
+                  2*stencil[dir] + offset);
+    }
+    
+    /* Allocate slab transfer description */
+    xferinfo = malloc(group.dim * sizeof *xferinfo);
+    assert (xferinfo);
+    
+    /* Fill in slab transfer description */
+    for (d=0; d<group.dim; ++d) {
+      xferinfo[d].src.gsh         = data.gsh        [d];
+      xferinfo[d].src.lbnd        = data.lbnd       [d];
+      xferinfo[d].src.lsh         = data.lsh        [d];
+      xferinfo[d].src.lbbox       = fake_bbox       [2*d  ];
+      xferinfo[d].src.ubbox       = fake_bbox       [2*d+1];
+      xferinfo[d].src.nghostzones = data.nghostzones[d];
+      xferinfo[d].src.off         = 0;
+      xferinfo[d].src.str         = 1;
+      xferinfo[d].src.len         = data.gsh        [d];
+      
+      xferinfo[d].dst.gsh         = data.gsh        [d];
+      xferinfo[d].dst.lbnd        = data.lbnd       [d];
+      xferinfo[d].dst.lsh         = data.lsh        [d];
+      xferinfo[d].dst.lbbox       = fake_bbox       [2*d  ];
+      xferinfo[d].dst.ubbox       = fake_bbox       [2*d+1];
+      xferinfo[d].dst.nghostzones = data.nghostzones[d];
+      xferinfo[d].dst.off         = 0;
+      xferinfo[d].dst.str         = 1;
+      xferinfo[d].dst.len         = data.gsh        [d];
+      
+      xferinfo[d].xpose = d;
+      xferinfo[d].flip  = 0;
+    }
+    
+    xferinfo[dir].src.off = stencil[dir] + offset[dir];
+    xferinfo[dir].src.len = stencil[dir];
+    xferinfo[dir].dst.off = 0;
+    xferinfo[dir].dst.len = stencil[dir];
+    
+    xferinfo[     dir].flip = 1;
+    xferinfo[otherdir].flip = 1;
+    
+    ierr = Slab_Transfer
+      (cctkGH, group.dim, xferinfo, -1,
+       group.vartype, varptr, group.vartype, varptr);
+    assert (!ierr);
+    
+    free (xferinfo);
+    
+  } /* if global_bbox */
+  
+  /* take parity into account */
+  if (cctkGH->cctk_bbox[2*dir]) {
+    int parity;
+    
+    parity = parities[0] * parities[1];
+    assert (abs(parity) == 1);
+    if (parity == -1) {
+      int i, j, k;
+      assert (dir == 0);
+      assert (group.dim == 3);
+      assert (group.vartype == CCTK_VARIABLE_REAL);
+      for (k=0; k<cctkGH->cctk_lsh[2]; ++k) {
+        for (j=0; j<cctkGH->cctk_lsh[1]; ++j) {
+          for (i=0; i<stencil[dir]; ++i) {
+            const int ind = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            ((CCTK_REAL *) varptr) [ind] *= -1;
+          }
+        }
+      }
+    }
+    
+  } /* if bbox */
+  
+  return 0;
+}
+
+
+
+void Rot180_ApplyBC (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  
+  int nvars;
+  CCTK_INT * restrict indices;
+  CCTK_INT * restrict faces;
+  CCTK_INT * restrict widths;
+  CCTK_INT * restrict tables;
+  int vi;
+  int dim;
+  int * restrict stencil;
+  int i;
+  int ierr;
+  
+  assert (cctkGH);
+  
+  nvars = Boundary_SelectedGVs (cctkGH, 0, 0, 0, 0, 0, 0);
+  assert (nvars>=0);
+  
+  indices = malloc (nvars * sizeof *indices);
+  assert (indices);
+  faces = malloc (nvars * sizeof *faces);
+  assert (faces);
+  widths = malloc (nvars * sizeof *widths);
+  assert (widths);
+  tables = malloc (nvars * sizeof *tables);
+  assert (tables);
+  
+  ierr =  Boundary_SelectedGVs
+    (cctkGH, nvars, indices, faces, widths, tables, 0);
+  assert (ierr == nvars);
+  
+  for (i=0; i<nvars; ++i) {
+    vi = indices[i];
+    assert (vi>=0 && vi<CCTK_NumVars());
+    
+    assert (widths[i] >= 0);
+    
+    dim = CCTK_GroupDimFromVarI (vi);
+    assert (dim>=0);
+    
+    stencil = malloc (dim * sizeof *stencil);
+    assert (stencil);
+    ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi);
+    assert (!ierr);
+    
+    ierr = BndRot180VI (cctkGH, stencil, vi);
+    assert (!ierr);
+    
+    free (stencil);
+  }
+  
+  free (indices);
+  free (faces);
+  free (widths);
+  free (tables);
+}