aboutsummaryrefslogtreecommitdiff
path: root/src/GeneralizedPolynomial-Uniform/template.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/template.c')
-rw-r--r--src/GeneralizedPolynomial-Uniform/template.c200
1 files changed, 87 insertions, 113 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
index 0bb409b..a567119 100644
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -20,6 +20,13 @@
3-D molecules.
@endhdesc
+ @hdate 28.Jan.2003
+ @hauthor Jonathan Thornburg <jthorn@aei.mpg.de>
+ @hdesc Support parameter-table entries
+ @var boundary_off_centering_tolerance and
+ @var boundary_extrapolation_tolerance.
+ @endhdesc
+
@desc
The following header files must be #included before
#including this file:
@@ -301,6 +308,9 @@
only briefly; see the InterpLocalUniform() documentation
and/or the thorn guide for this thorn for further details.
+ If you change the arguments for this function, note that you
+ must also change the prototype in "template.h".
+
@var error_flags
@vdesc If we encounter an out-of-range point and this pointer
is non-NULL, we store a description of the out-of-range
@@ -346,13 +356,15 @@
@@*/
int FUNCTION_NAME(/***** coordinate system *****/
- const CCTK_REAL coord_system_origin[],
- const CCTK_REAL grid_spacing[],
+ const CCTK_REAL coord_origin[],
+ const CCTK_REAL coord_delta[],
/***** interpolation points *****/
int N_interp_points,
int interp_coords_type_code,
const void* const interp_coords[],
- const CCTK_REAL out_of_range_tolerance[],
+ const CCTK_INT N_boundary_points_to_omit[],
+ const CCTK_REAL boundary_off_centering_tolerance[],
+ const CCTK_REAL boundary_extrapolation_tolerance[],
/***** input arrays *****/
int N_input_arrays,
const CCTK_INT input_array_offsets[],
@@ -691,8 +703,8 @@ int out;
*/
{
#if N_DIMS >= 1
- const fp origin_x = coord_system_origin[X_AXIS];
- const fp delta_x = grid_spacing[X_AXIS];
+ const fp origin_x = coord_origin[X_AXIS];
+ const fp delta_x = coord_delta[X_AXIS];
#if defined(HAVE_OP_DX) \
|| defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \
|| defined(HAVE_OP_DXX)
@@ -700,8 +712,8 @@ int out;
#endif
#endif
#if N_DIMS >= 2
- const fp origin_y = coord_system_origin[Y_AXIS];
- const fp delta_y = grid_spacing[Y_AXIS];
+ const fp origin_y = coord_origin[Y_AXIS];
+ const fp delta_y = coord_delta[Y_AXIS];
#if defined(HAVE_OP_DY) \
|| defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \
|| defined(HAVE_OP_DYY)
@@ -709,8 +721,8 @@ int out;
#endif
#endif
#if N_DIMS >= 3
- const fp origin_z = coord_system_origin[Z_AXIS];
- const fp delta_z = grid_spacing[Z_AXIS];
+ const fp origin_z = coord_origin[Z_AXIS];
+ const fp delta_z = coord_delta[Z_AXIS];
#if defined(HAVE_OP_DZ) \
|| defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \
|| defined(HAVE_OP_DZZ)
@@ -764,6 +776,7 @@ int out;
/*
* interpolate at each point
*/
+int return_code = 0;
int pt;
for (pt = 0 ; pt < N_interp_points ; ++pt)
{
@@ -890,124 +903,84 @@ int pt;
*
* n.b. we need the final answers in variables with the magic
* names (x,y,z) (machine-generated code uses these names),
- * but we use temp variables as intermediates for (likely)
- * better performance: the temp variables have their addresses
- * taken and so may not be register-allocated, whereas the
- * final (x,y,z) are "clean" and thus more likely to be
- * register-allocated
+ * but we use temp variables as intermediates for these and
+ * for center_[ijk] for (likely) better performance:
+ * the temp variables have their addresses taken and so
+ * may not be register-allocated, whereas the final variables
+ * are "clean" and thus more likely to be register-allocated
*/
{
- #if (N_DIMS >= 1)
- fp x_temp;
- const int center_i
- = LocalInterp_molecule_posn(origin_x, delta_x,
- input_array_min_subscripts[X_AXIS],
- input_array_max_subscripts[X_AXIS],
- MOLECULE_SIZE,
- out_of_range_tolerance[X_AXIS_MIN],
- out_of_range_tolerance[X_AXIS_MAX],
- interp_coords_fp[X_AXIS],
- &x_temp,
- (int *) NULL, (int *) NULL);
- const fp x = x_temp;
- #endif
- #if (N_DIMS >= 2)
- fp y_temp;
- const int center_j
- = LocalInterp_molecule_posn(origin_y, delta_y,
- input_array_min_subscripts[Y_AXIS],
- input_array_max_subscripts[Y_AXIS],
- MOLECULE_SIZE,
- out_of_range_tolerance[Y_AXIS_MIN],
- out_of_range_tolerance[Y_AXIS_MAX],
- interp_coords_fp[Y_AXIS],
- &y_temp,
- (int *) NULL, (int *) NULL);
- const fp y = y_temp;
- #endif
- #if (N_DIMS >= 3)
- fp z_temp;
- const int center_k
- = LocalInterp_molecule_posn(origin_z, delta_z,
- input_array_min_subscripts[Z_AXIS],
- input_array_max_subscripts[Z_AXIS],
- MOLECULE_SIZE,
- out_of_range_tolerance[Z_AXIS_MIN],
- out_of_range_tolerance[Z_AXIS_MAX],
- interp_coords_fp[Z_AXIS],
- &z_temp,
- (int *) NULL, (int *) NULL);
- const fp z = z_temp;
- #endif
- #if (N_DIMS > 3)
- #error "N_DIMS may not be > 3!"
- #endif
-
- /* is the interpolation point out-of-range? */
- #if (N_DIMS >= 1)
- if (IS_MOLECULE_POSN_ERROR_CODE(center_i))
- then {
- if (center_i == MOLECULE_POSN_ERROR_GRID_TINY)
- then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
- if (error_flags != NULL)
- then {
- error_flags->error_pt = pt;
- error_flags->error_axis = X_AXIS;
- error_flags->error_end
- = (center_i == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1;
- }
- return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
- }
- #endif
- #if (N_DIMS >= 2)
- if (IS_MOLECULE_POSN_ERROR_CODE(center_j))
- then {
- if (center_j == MOLECULE_POSN_ERROR_GRID_TINY)
- then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
- if (error_flags != NULL)
- then {
- error_flags->error_pt = pt;
- error_flags->error_axis = Y_AXIS;
- error_flags->error_end
- = (center_j == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1;
- }
- return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
- }
- #endif
- #if (N_DIMS >= 3)
- if (IS_MOLECULE_POSN_ERROR_CODE(center_k))
- then {
- if (center_k == MOLECULE_POSN_ERROR_GRID_TINY)
- then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
- if (error_flags != NULL)
- then {
- error_flags->error_pt = pt;
- error_flags->error_axis = Z_AXIS;
- error_flags->error_end
- = (center_k == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1;
- }
- return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
- }
- #endif
+ int center_ijk_temp[MAX_N_DIMS];
+ fp xyz_temp[MAX_N_DIMS];
+ int axis;
+ for (axis = 0 ; axis < N_DIMS ; ++axis)
+ {
+ const int ibndry_min = 2*axis;
+ const int ibndry_max = 2*axis + 1;
+ const int status = LocalInterp_molecule_posn
+ (coord_origin[axis], coord_delta[axis],
+ input_array_min_subscripts[axis]
+ + N_boundary_points_to_omit[ibndry_min],
+ input_array_max_subscripts[axis]
+ - N_boundary_points_to_omit[ibndry_max],
+ MOLECULE_SIZE,
+ boundary_off_centering_tolerance[ibndry_min],
+ boundary_off_centering_tolerance[ibndry_max],
+ boundary_extrapolation_tolerance[ibndry_min],
+ boundary_extrapolation_tolerance[ibndry_max],
+ interp_coords_fp[axis],
+ &center_ijk_temp[axis], &xyz_temp[axis]);
+ if (status < 0)
+ then {
+ if (status == MOLECULE_POSN_ERROR_GRID_TINY)
+ then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
+ /*** ERROR RETURN ***/
+ if (error_flags != NULL)
+ then {
+ error_flags->error_pt = pt;
+ error_flags->error_axis = axis;
+ if (status == MOLECULE_POSN_ERROR_X_LT_MIN)
+ then {
+ error_flags->error_ibndry = ibndry_min;
+ error_flags->error_end = -1;
+ }
+ else {
+ error_flags->error_ibndry = ibndry_max;
+ error_flags->error_end = +1;
+ }
+ }
+ return CCTK_ERROR_INTERP_POINT_OUTSIDE; /*** ERROR RETURN ***/
+ }
+ }
+ {
+ #if (N_DIMS >= 1)
+ const int center_i = center_ijk_temp[X_AXIS];
+ const fp x = xyz_temp [X_AXIS];
+ #endif
+ #if (N_DIMS >= 2)
+ const int center_j = center_ijk_temp[Y_AXIS];
+ const fp y = xyz_temp [Y_AXIS];
+ #endif
+ #if (N_DIMS >= 3)
+ const int center_k = center_ijk_temp[Z_AXIS];
+ const fp z = xyz_temp [Z_AXIS];
+ #endif
/*
* compute 1-d position of molecule center in input arrays
*/
{
- #if (N_DIMS == 1)
+ #if (N_DIMS == 1)
const int molecule_center_posn = stride_i*center_i;
- #endif
- #if (N_DIMS == 2)
+ #elif (N_DIMS == 2)
const int molecule_center_posn = stride_i*center_i
+ stride_j*center_j;
- #endif
- #if (N_DIMS == 3)
+ #elif (N_DIMS == 3)
const int molecule_center_posn = stride_i*center_i
+ stride_j*center_j
+ stride_k*center_k;
- #endif
- #if (N_DIMS > 3)
+ #else
#error "N_DIMS may not be > 3!"
#endif
@@ -1566,6 +1539,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
}
}
}
+ }
/* end of for (pt = ...) loop */
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-20 16:10:53 +0000