aboutsummaryrefslogtreecommitdiff
path: root/src/interpolate.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/interpolate.c')
-rw-r--r--src/interpolate.c431
1 files changed, 222 insertions, 209 deletions
diff --git a/src/interpolate.c b/src/interpolate.c
index e4c242d..eb60460 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -67,10 +67,8 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
do_reflection[5] = reflection_upper_z;
newfaces = faces;
- for (d=0; d<6; ++d)
- {
- if (do_reflection[d])
- {
+ for (d=0; d<6; ++d) {
+ if (do_reflection[d]) {
assert (newfaces & (1 << d));
newfaces &= ~ (1 << d);
}
@@ -80,25 +78,20 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
/* Fold coordinates */
assert (interp_coords_type == CCTK_VARIABLE_REAL);
- for (dir=0; dir<3; ++dir)
- {
+ for (dir=0; dir<3; ++dir) {
assert (! do_reflection[2*dir+1]);
- if (do_reflection[2*dir])
- {
+ 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)
- {
+ 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
- {
+ } else {
new_interp_coords[dir] = interp_coords[dir];
}
}
@@ -117,15 +110,14 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
/* Free coordinates */
- for (dir=0; dir<3; ++dir)
- {
- if (do_reflection[2*dir])
- {
+ 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);
@@ -167,158 +159,185 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
/* 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, WEYLSCALARS_REAL };
- 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) {
+ 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, WEYLSCALARS_REAL };
+ 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 (2, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element",
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Error in tensor type alias declaration for group \"%s\"",
groupname);
free (groupname);
}
- ttype = SCALAR;
+
+ ttype = UNKNOWN;
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 if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
- /* Weyl scalars, stored as 10 real numbers */
- assert (numvars == 10);
- ttype = WEYLSCALARS_REAL;
- 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;
- case WEYLSCALARS_REAL:
- assert (tcomponent>=0 && tcomponent<10);
- 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)
+ 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, "4scalar")) {
+ /* 4-scalar */
+ assert (numvars == 1);
+ ttype = SCALAR;
+ tcomponent = 0;
+ } else if (CCTK_EQUALS (tensortypealias, "u")
+ || CCTK_EQUALS (tensortypealias, "d"))
{
- 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;
- case WEYLSCALARS_REAL:
+ /* vector */
+ assert (numvars == 3);
+ ttype = VECTOR;
+ tcomponent = vi - firstvar;
+ } else if (CCTK_EQUALS (tensortypealias, "4u")
+ || CCTK_EQUALS (tensortypealias, "4d"))
+ {
+ /* 4-vector */
+ assert (numvars == 4);
+ if (vi == firstvar) {
+ ttype = SCALAR;
+ tcomponent = 0;
+ } else {
+ ttype = VECTOR;
+ tcomponent = vi - firstvar - 1;
+ }
+ } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
+ || CCTK_EQUALS (tensortypealias, "dd_sym"))
+ {
+ /* symmetric tensor */
+ assert (numvars == 6);
+ ttype = TENSOR;
+ tcomponent = vi - firstvar;
+ } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
+ || CCTK_EQUALS (tensortypealias, "4dd_sym"))
{
+ /* symmetric 4-tensor */
+ assert (numvars == 10);
+ if (vi == firstvar) {
+ ttype = SCALAR;
+ tcomponent = 0;
+ } else if (vi <= firstvar+3) {
+ ttype = VECTOR;
+ tcomponent = vi - firstvar - 1;
+ } else {
+ ttype = TENSOR;
+ tcomponent = vi - firstvar - 4;
+ }
+ } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
+ /* Weyl scalars, stored as 10 real numbers */
+ assert (numvars == 10);
+ ttype = WEYLSCALARS_REAL;
+ 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;
+ case WEYLSCALARS_REAL:
+ assert (tcomponent>=0 && tcomponent<10);
+ 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;
+ case WEYLSCALARS_REAL: {
static int const weylparities[10][3] =
{{+1,+1,+1},
{-1,-1,-1},
@@ -333,62 +352,56 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
for (dir=0; dir<3; ++dir) {
parities[dir] = weylparities[tcomponent][dir];
}
+ break;
}
- 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;
+ default:
+ assert (0);
}
- }
-
-
-
- /* 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)
+
+
+
+ /* Take derivatives into account */
{
- 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;
+ 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;
}
- ((CCTK_REAL *)output_arrays[m])[n] *= parity;
}
+
}
-
- }
- } /* for m */
+ } /* for m */
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-18 11:33:18 +0000