diff options
Diffstat (limited to 'src/interpolate.c')
| -rw-r--r-- | src/interpolate.c | 359 |
1 files changed, 359 insertions, 0 deletions
diff --git a/src/interpolate.c b/src/interpolate.c new file mode 100644 index 0000000..8591fd7 --- /dev/null +++ b/src/interpolate.c @@ -0,0 +1,359 @@ +#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 *new_interp_coords[dir]); + assert (new_interp_coords[d]); + + 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]>=0 + && output_array_indices[m]<CCTK_NumVars()); + } + + + + /* Unfold tensor types */ + for (m=0; m<N_output_arrays; ++m) + { + + int vi, gi; + int numvars, firstvar; + char * groupname; + + int table; + char tensortypealias[1000]; + enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR }; + enum tensortype ttype; + 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); + + + + /* 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")) + { + /* vector */ + assert (numvars == 3); + ttype = VECTOR; + tcomponent = vi - firstvar; + } else if (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); + } + + + + /* 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[d] = do_reflection[2*dir] && parities[dir] < 0; + needs_checking |= check_dir[d]; + } + + + + /* Loop over all points and unfold */ + if (needs_checking) + { + for (n=0; n<N_interp_points; ++n) + { + int parity = +1; + /* 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; +} |