path: root/src/interpolate.c

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "cctk_Parameters.h"

#include "util_ErrorCodes.h"
#include "util_Table.h"

#include "reflection.h"



static const char rcsid[] = "$Header$";
CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_interpolate_c);



CCTK_INT
ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict 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 restrict 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 restrict const output_arrays[],
                                CCTK_INT const faces)
{
  cGH const * restrict const cctkGH = cctkGH_;
  DECLARE_CCTK_PARAMETERS;
  
  int do_reflection[6];
  
  CCTK_POINTER_TO_CONST restrict new_interp_coords[3];
  CCTK_INT newfaces;
  
  CCTK_INT * restrict operand_indices;
  CCTK_INT * restrict operation_codes;
  CCTK_INT * restrict output_array_indices;
  
  int dir;
  
  int iret;
  
  int m;
  int n;
  int d;
  
  int ierr;
  
  
  
  /* Get symmetry information */
  do_reflection[0] = reflection_x;
  do_reflection[1] = reflection_upper_x;
  do_reflection[2] = reflection_y;
  do_reflection[3] = reflection_upper_y;
  do_reflection[4] = reflection_z;
  do_reflection[5] = reflection_upper_z;
  
  newfaces = faces;
  for (d=0; d<6; ++d)
  {
    if (do_reflection[d])
    {
      assert (newfaces & (1 << d));
      newfaces &= ~ (1 << d);
    }
  }
  
  
  
  /* Fold coordinates */
  assert (interp_coords_type == CCTK_VARIABLE_REAL);
  for (dir=0; dir<3; ++dir)
  {
    assert (! do_reflection[2*dir+1]);
    
    if (do_reflection[2*dir])
    {
      new_interp_coords[dir]
        = malloc (N_interp_points * sizeof (CCTK_REAL));
      assert (N_interp_points == 0 || new_interp_coords[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;
      }
    }
    else
    {
      new_interp_coords[dir] = interp_coords[dir];
    }
  }
  
  
  
  /* 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,
     newfaces);
  
  
  
  /* Free coordinates */
  for (dir=0; dir<3; ++dir)
  {
    if (do_reflection[2*dir])
    {
      free (new_interp_coords[dir]);
    }
  }
  
  
  /* 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]==-1
             || (output_array_indices[m]>=0
                 && output_array_indices[m]<CCTK_NumVars()));
  }
  
  
  
  /* Unfold tensor types */
  for (m=0; m<N_output_arrays; ++m)
  {
  if (output_array_indices[m]!=-1)
  {
    
    int vi, gi;
    int numvars, firstvar;
    char * groupname;
    
    int table;
    char tensortypealias[1000];
    enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR };
    enum tensortype ttype;
    CCTK_INT tensorparity;
    int tcomponent;
    
    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());
    numvars = CCTK_NumVarsInGroupI(gi);
    assert (numvars>0);
    firstvar = CCTK_FirstVarIndexI(gi);
    assert (firstvar>=0);
    table = CCTK_GroupTagsTableI(gi);
    assert (table>=0);

    
    
    /* Get and check tensor type information */
    ierr = Util_TableGetString
      (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
      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) {
      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);
    }
    
    ttype = UNKNOWN;
    tcomponent = 0;
    if (CCTK_EQUALS (tensortypealias, "scalar"))
    {
      /* scalar */
      if (numvars != 1) {
        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);
      }
      ttype = SCALAR;
      tcomponent = 0;
    }
    else if (CCTK_EQUALS (tensortypealias, "u")
             || CCTK_EQUALS (tensortypealias, "d"))
    {
      /* vector */
      assert (numvars == 3);
      ttype = VECTOR;
      tcomponent = vi - firstvar;
    } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
               || CCTK_EQUALS (tensortypealias, "dd_sym")) {
      /* symmetric tensor */
      assert (numvars == 6);
      ttype = TENSOR;
      tcomponent = vi - firstvar;
    } 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);
    }
    
    switch (ttype)
    {
    case SCALAR:
      assert (tcomponent>=0 && tcomponent<1);
      break;
    case VECTOR:
      assert (tcomponent>=0 && tcomponent<3);
      break;
    case TENSOR:
      assert (tcomponent>=0 && tcomponent<6);
      break;
    default:
      assert (0);
    }
    
    ierr = Util_TableGetInt (table, & tensorparity, "tensorparity");
    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
      tensorparity = +1;
    } else if (ierr<0) {
      groupname = CCTK_GroupName(gi);
      assert (groupname);
      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Error in tensor parity declaration for group \"%s\"",
                  groupname);
      free (groupname);
    }
    
    
    
    /* Calculate parities */
    parities[0] = parities[1] = parities[2] = +1;
    switch (ttype)
    {
    case SCALAR:
      /* do nothing */
      break;
    case VECTOR:
      parities[tcomponent] = -1;
      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: break;
      case 4: parities[1] = parities[2] = -1; break;
      case 5: break;
      default: assert (0);
      }
      break;
    default:
      assert (0);
    }
    
    
    
    /* Take derivatives into account */
    {
      int code = operation_codes[m];
      while (code) {
        const int d = code % 10 - 1;
        code /= 10;
        assert (d>=0 && d<3);
        parities[d] *= -1;
      }
    }
    
    
    
    /* Are there negative parities? */
    needs_checking = 0;
    for (dir=0; dir<3; ++dir)
    {
      check_dir[dir] = do_reflection[2*dir] && parities[dir] < 0;
      needs_checking |= check_dir[dir];
    }
    
    
    
    /* Loop over all points and unfold */
    if (needs_checking)
    {
      for (n=0; n<N_interp_points; ++n)
      {
        int parity = tensorparity;
        /* Is the point outside the domain? */
        for (dir=0; dir<3; ++dir)
        {
          if (check_dir[dir])
          {
            CCTK_REAL const pos = ((CCTK_REAL const *)interp_coords[dir])[n];
            if (pos < 0)
            {
              /* Reflect the tensor component */
              parity *= -1;
            }
          }
        }
        ((CCTK_REAL *)output_arrays[m])[n] *= parity;
      }
    }
    
  }
  } /* for m */
  
  
  
  /* Free output variable indices */
  free (operand_indices);
  free (operation_codes);
  free (output_array_indices);
  
  
  
  /* Return */
  return iret;
}