diff options
Diffstat (limited to 'src/interpolate.c')
| -rw-r--r-- | src/interpolate.c | 79 |
1 files changed, 48 insertions, 31 deletions
diff --git a/src/interpolate.c b/src/interpolate.c index 8362db4..571c048 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -76,15 +76,16 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, assert (! do_reflection[2*dir+1]); if (do_reflection[2*dir]) { - new_interp_coords[dir] + CCTK_REAL * restrict coords_in_dir = malloc (N_interp_points * sizeof (CCTK_REAL)); - assert (N_interp_points == 0 || new_interp_coords[dir]); + assert (N_interp_points == 0 || coords_in_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; + coords_in_dir[n] = newpos; } + new_interp_coords[dir] = coords_in_dir; } else { new_interp_coords[dir] = interp_coords[dir]; } @@ -96,9 +97,9 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, iret = SymmetryInterpolateFaces (cctkGH, N_dims, local_interp_handle, param_table_handle, coord_system_handle, - N_interp_points, interp_coords_type, new_interp_coords, + N_interp_points, interp_coords_type, (CCTK_POINTER_TO_CONST)new_interp_coords, N_input_arrays, input_array_indices, - N_output_arrays, output_array_types, output_arrays, + N_output_arrays, output_array_types, (CCTK_POINTER_TO_CONST)output_arrays, newfaces); @@ -106,7 +107,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, /* Free coordinates */ for (dir=0; dir<3; ++dir) { if (do_reflection[2*dir]) { - free (new_interp_coords[dir]); + free ((void*)new_interp_coords[dir]); } } @@ -157,7 +158,8 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, if (output_array_indices[m]!=-1) { int vi, gi; - int numvars, firstvar; + int numvars, firstvar, vectorlength; + cGroup group; char * groupname; int table; @@ -172,17 +174,23 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, 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()); + + ierr = CCTK_GroupData (gi, &group); + assert (!ierr); + numvars = CCTK_NumVarsInGroupI(gi); assert (numvars>0); firstvar = CCTK_FirstVarIndexI(gi); assert (firstvar>=0); + vectorlength = group.vectorlength; + assert (vectorlength>=0); + assert (vectorlength==1 || group.vectorgroup); table = CCTK_GroupTagsTableI(gi); assert (table>=0); @@ -212,7 +220,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, tcomponent = 0; if (CCTK_EQUALS (tensortypealias, "scalar")) { /* scalar */ - if (numvars != 1) { + if (numvars != vectorlength) { groupname = CCTK_GroupName(gi); assert (groupname); CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING, @@ -224,72 +232,81 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, tcomponent = 0; } else if (CCTK_EQUALS (tensortypealias, "4scalar")) { /* 4-scalar */ - assert (numvars == 1); + assert (numvars == vectorlength); ttype = SCALAR; tcomponent = 0; } else if (CCTK_EQUALS (tensortypealias, "u") || CCTK_EQUALS (tensortypealias, "d")) { /* vector */ - assert (numvars == 3); + const int numcomps = 3; + /* special case to handle things like vel[3] */ + assert (numvars % numcomps == 0 && + (numvars == numcomps * vectorlength || numvars == vectorlength)); ttype = VECTOR; - tcomponent = vi - firstvar; + tcomponent = (vi - firstvar) % numcomps; } else if (CCTK_EQUALS (tensortypealias, "4u") || CCTK_EQUALS (tensortypealias, "4d")) { /* 4-vector */ - assert (numvars == 4); - if (vi == firstvar) { + const int numcomps = 4; + assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength); + if ((vi - firstvar) % numcomps == 0) { ttype = SCALAR; tcomponent = 0; } else { ttype = VECTOR; - tcomponent = vi - firstvar - 1; + tcomponent = (vi - firstvar) % numcomps - 1; } } else if (CCTK_EQUALS (tensortypealias, "uu_sym") || CCTK_EQUALS (tensortypealias, "dd_sym")) { /* symmetric tensor */ - assert (numvars == 6); + const int numcomps = 6; + assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength); ttype = SYMTENSOR; - tcomponent = vi - firstvar; + tcomponent = (vi - firstvar) % numcomps; } else if (CCTK_EQUALS (tensortypealias, "uu") || CCTK_EQUALS (tensortypealias, "ud") || CCTK_EQUALS (tensortypealias, "du") || CCTK_EQUALS (tensortypealias, "dd")) { /* non-symmetric tensor */ - assert (numvars == 9); + const int numcomps = 9; + assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength); ttype = TENSOR; - tcomponent = vi - firstvar; + tcomponent = (vi - firstvar) % numcomps; } else if (CCTK_EQUALS (tensortypealias, "4uu_sym") || CCTK_EQUALS (tensortypealias, "4dd_sym")) { /* symmetric 4-tensor */ - assert (numvars == 10); - if (vi == firstvar) { + const int numcomps = 10; + assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength); + if ((vi - firstvar) % numcomps == 0) { ttype = SCALAR; tcomponent = 0; - } else if (vi <= firstvar+3) { + } else if ((vi - firstvar) % numcomps <= 3) { ttype = VECTOR; - tcomponent = vi - firstvar - 1; + tcomponent = (vi - firstvar) % numcomps - 1; } else { ttype = SYMTENSOR; - tcomponent = vi - firstvar - 4; + tcomponent = (vi - firstvar) % numcomps - 4; } } else if (CCTK_EQUALS (tensortypealias, "ddd_sym")) { /* 3rd rank tensor, symmetric in last 2 indices */ - assert (numvars == 18); + const int numcomps = 18; + assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength); ttype = SYMTENSOR3; - tcomponent = vi - firstvar; + tcomponent = (vi - firstvar) % numcomps; } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) { /* Weyl scalars, stored as 10 real numbers */ - assert (numvars == 10); + const int numcomps = 10; + assert (numvars % numcomps == 0 && numvars == numcomps * vectorlength); ttype = WEYLSCALARS_REAL; - tcomponent = vi - firstvar; + tcomponent = (vi - firstvar) % numcomps; } else if (CCTK_EQUALS (tensortypealias, "ManualCartesian")) { - /* Reflection symmetries specified by hand */ - ttype = MANUALCARTESIAN; + /* Reflection symmetries specified by hand */ + ttype = MANUALCARTESIAN; tcomponent = vi - firstvar; } else { groupname = CCTK_GroupName(gi); |