Provide symmetry interpolation.

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

4 files changed, 286 insertions, 2 deletions

diff --git a/src/interpolate.c b/src/interpolate.c
new file mode 100644
index 0000000..f48aa1d
--- /dev/null
+++ b/src/interpolate.c
@@ -0,0 +1,274 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "rotatingsymmetry180.h"
+
+
+
+/* Symmetry interpolation */
+CCTK_INT
+Rot180_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
+                            CCTK_INT const N_dims,
+                            CCTK_INT const local_interp_handle,
+                            CCTK_INT const param_table_handle,
+                            CCTK_INT const coord_system_handle,
+                            CCTK_INT const N_interp_points,
+                            CCTK_INT const interp_coords_type,
+                            CCTK_POINTER_TO_CONST const interp_coords[],
+                            CCTK_INT const N_input_arrays,
+                            CCTK_INT const input_array_indices[],
+                            CCTK_INT const N_output_arrays,
+                            CCTK_INT const output_array_types[],
+                            CCTK_POINTER const output_arrays[],
+                            CCTK_INT const faces)
+{
+  cGH const * restrict const cctkGH = cctkGH_;
+  CCTK_POINTER new_interp_coords[3];
+  CCTK_INT new_faces;
+  CCTK_INT * restrict operand_indices;
+  CCTK_INT * restrict operation_codes;
+  CCTK_INT * restrict output_array_indices;
+  int parities[3];              /* variable parities */
+  int m;                        /* output array */
+  int n;                        /* point */
+  int f;                        /* face */
+  int d;                        /* dimension */
+  int iret;                     /* interpolator return value */
+  int ierr;
+  
+  /* Check arguments */
+  assert (N_dims == 3);
+  assert (N_interp_points >= 0);
+  assert (interp_coords_type >= 0);
+  for (d=0; d<3; ++d) {
+    assert (N_interp_points == 0 || interp_coords[d]);
+  }
+  assert (N_output_arrays >= 0);
+  
+  /* Coordinates must be CCTK_REAL */
+  assert (interp_coords_type == CCTK_VARIABLE_REAL);
+  for (m=0; m<N_output_arrays; ++m) {
+    assert (output_array_types[m] == CCTK_VARIABLE_REAL);
+  }
+  
+  
+  
+  /* Claim faces */
+  assert (faces & (1 << 0));
+  new_faces = faces;
+  new_faces &= ~ (1 << 0);
+  
+  
+  
+  /* Copy coordinates */
+  for (d=0; d<3; ++d) {
+    new_interp_coords[d] = malloc (N_interp_points * sizeof(CCTK_REAL));
+    assert (new_interp_coords[d]);
+    memcpy (new_interp_coords[d], interp_coords[d],
+            N_interp_points * sizeof(CCTK_REAL));
+  }
+  
+  /* Fold coordinates */
+  for (n=0; n<N_interp_points; ++n) {
+    /* Is the point outside the domain? */
+    if (((CCTK_REAL *)new_interp_coords[0])[n] < 0) {
+      /* Rotate the point by 180 degrees */
+      ((CCTK_REAL *)new_interp_coords[0])[n] *= -1;
+      ((CCTK_REAL *)new_interp_coords[1])[n] *= -1;
+    }
+  }
+  
+  
+  
+  /* Recursive call */
+  iret = SymmetryInterpolateFaces
+    (cctkGH_, N_dims,
+     local_interp_handle, param_table_handle, coord_system_handle,
+     N_interp_points, interp_coords_type, new_interp_coords,
+     N_input_arrays, input_array_indices,
+     N_output_arrays, output_array_types, output_arrays,
+     new_faces);
+  
+  
+  
+  /* Free coordinates */
+  for (d=0; d<3; ++d) {
+    free (new_interp_coords[d]);
+  }
+  
+  
+  
+  /* Find output variable indices */
+  operand_indices = malloc (N_output_arrays * sizeof *operand_indices);
+  assert (operand_indices);
+  ierr = Util_TableGetIntArray
+    (param_table_handle, N_output_arrays, operand_indices, "operand_indices");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    assert (N_output_arrays == N_input_arrays);
+    for (m=0; m<N_output_arrays; ++m) {
+      operand_indices[m] = m;   /* set output index to input index */
+    }
+  } else {
+    assert (ierr == N_output_arrays);
+  }
+  
+  operation_codes = malloc (N_output_arrays * sizeof *operation_codes);
+  assert (operation_codes);
+  ierr = Util_TableGetIntArray
+    (param_table_handle, N_output_arrays, operation_codes, "operation_codes");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    assert (N_output_arrays == N_input_arrays);
+    for (m=0; m<N_output_arrays; ++m) {
+      operation_codes[m] = 0;     /* do not take derivatives */
+    }
+  } else {
+    assert (ierr == N_output_arrays);
+  }
+  
+  output_array_indices
+    = malloc (N_output_arrays * sizeof *output_array_indices);
+  assert (output_array_indices);
+  for (m=0; m<N_output_arrays; ++m) {
+    assert (operand_indices[m]>=0 && operand_indices[m]<N_input_arrays);
+    output_array_indices[m] = input_array_indices[operand_indices[m]];
+     assert (output_array_indices[m]>=0
+             && output_array_indices[m]<CCTK_NumVars());
+  }
+  
+  
+  
+  /* Loop over all output arrays */
+  for (m=0; m<N_output_arrays; ++m) {
+    
+    
+    
+    /* Find variable parity */
+    {
+      int vi, gi;
+      char tensortypealias[100];
+      int firstvar, numvars;
+      int table;
+      int index;
+      
+      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);
+      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);
+      }
+      
+      for (d=0; d<3; ++d) {
+        assert (abs(parities[d])==1);
+      }
+    } /* Find parity */
+    
+    
+    
+    /* Modify parity */
+    {
+      int code = operation_codes[m];
+      while (code) {
+        const int dir = code % 10;
+        code /= 10;
+        assert (dir>=1 && dir<=3);
+        parities[dir-1] *= -1;
+      }
+    }
+    
+    
+    
+    /* Loop over all points and unfold */
+    if (parities[0] * parities[1] == -1) {
+      for (n=0; n<N_interp_points; ++n) {
+        /* Is the point outside the domain? */
+        if (((const CCTK_REAL *)interp_coords[0])[n] < 0) {
+          /* Rotate the tensor components back by 180 degrees */
+          ((CCTK_REAL *)output_arrays[m])[n] *= -1;
+        }
+      }
+    }
+    
+    
+    
+  } /* for m */
+  
+  
+  
+  /* Free output variable descriptions */
+  free (operand_indices);
+  free (operation_codes);
+  free (output_array_indices);
+  
+  
+
+  return iret;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index 2b5b357..03f55f4 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = rotatingsymmetry180.c registersymmetry.c
+SRCS = interpolate.c rotatingsymmetry180.c registersymmetry.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/registersymmetry.c b/src/registersymmetry.c
index 50a286c..b326181 100644
--- a/src/registersymmetry.c
+++ b/src/registersymmetry.c
@@ -4,6 +4,10 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
+#include "rotatingsymmetry180.h"
+
+
+
 void Rot180_RegisterSymmetry (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
@@ -32,4 +36,10 @@ void Rot180_RegisterSymmetry (CCTK_ARGUMENTS)
   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");
   }
+  
+  ierr = SymmetryRegisterGridInterpolator
+    (cctkGH, handle, Rot180_SymmetryInterpolate);
+  if (ierr < 0) {
+    CCTK_WARN (0, "Could not register the symmetry interpolator");
+  }
 }
diff --git a/src/rotatingsymmetry180.h b/src/rotatingsymmetry180.h
index d795e8c..e560f20 100644
--- a/src/rotatingsymmetry180.h
+++ b/src/rotatingsymmetry180.h
@@ -27,6 +27,6 @@ Rot180_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
                             CCTK_INT const N_output_arrays,
                             CCTK_INT const output_array_types[],
                             CCTK_POINTER const output_arrays[],
-                            CCTK_INT const faces[]);
+                            CCTK_INT const faces);
 
 #endif /* ! defined ROTATINGSYMMETRY180_H */