diff options
-rw-r--r-- | src/patch/patch_interp.cc | 51 | ||||
-rw-r--r-- | src/patch/patch_interp.hh | 25 |
2 files changed, 19 insertions, 57 deletions
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index ee12e10..421405e 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -61,7 +61,9 @@ patch_interp::patch_interp(const patch_edge& my_edge_in, interp_par_(interp_par_in), interp_handle_(interp_handle_in), interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)), - gridfn_coord_origin_(my_edge().par_map().origin()), + // note interpolator coordinate origin = par of min_ipar() point + // cf. comments in "patch_interp.hh" on gridfn array subscripting + gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar())), gridfn_coord_delta_ (my_edge().par_map().delta_fp()), gridfn_type_codes_ (min_gfn_, max_gfn_), gridfn_data_ptrs_ (min_gfn_, max_gfn_), @@ -81,8 +83,8 @@ if (status < 0) status); /*NOTREACHED*/ // -// set up the offset, stride, and min/max subscripts in the -// parameter table so the interpolator can access the gridfn data arrays +// set up the gridfn array stride in the parameter table +// so the interpolator can access the gridfn data arrays // const int N_dims = 1; @@ -100,49 +102,6 @@ if (status < 0) , interp_par_table_handle_, status); /*NOTREACHED*/ -// offsets (same for all gfn) -jtutil::array1d offsets(min_gfn_, max_gfn_); - for (int gfn = min_gfn_ ; gfn <= max_gfn_ ; ++gfn) - { - offsets(gfn) = - min_ipar_ * stride; - } -const int N_gridfns = jtutil::how_many_in_range(min_gfn_, max_gfn_); -// n.b. Util_TableSetIntArray() *copies* our array, -// so it doesn't matter that the array is a local variable -// in this function and will disappear when this function returns -status = Util_TableSetIntArray(interp_par_table_handle_, - N_gridfns, &offsets(min_gfn_), - "input_array_offsets"); -if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::patch_interp():\n" -" can't set gridfn offsets in interpolator parmameter table!\n" -" table handle=%d error status=%d\n" -, - interp_par_table_handle_, status); /*NOTREACHED*/ - -// min/max subscripts -status = Util_TableSetIntArray(interp_par_table_handle_, - N_dims, &min_ipar_, - "input_array_min_subscripts"); -if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::patch_interp():\n" -" can't set min gridfn subscript in interpolator parmameter table!\n" -" table handle=%d error status=%d\n" -, - interp_par_table_handle_, status); /*NOTREACHED*/ -status = Util_TableSetIntArray(interp_par_table_handle_, - N_dims, &max_ipar_, - "input_array_max_subscripts"); -if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::patch_interp():\n" -" can't set max gridfn subscript in interpolator parmameter table!\n" -" table handle=%d error status=%d\n" -, - interp_par_table_handle_, status); /*NOTREACHED*/ - // // set up the gridfn type codes for the interpolator // diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index a213b43..c8cf8a0 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -89,17 +89,20 @@ // the gridfn data is a bit tricky: The interpolator accesses the // gridfn data using the generic (1-D) subscripting expression // data[offset + i*stride] -// where i is the data array index. It's convenient to take i=ipar -// (this way the mapping from i to par is just given by the usual -// origin and delta values for the par coordinate). The stride is -// fixed by the gridfn storage mapping (we get this from the patch). -// We then calculate the offset such that offset + min_ipar*stride == 0 , -// and then pass the data pointer for each interpolator call as a pointer -// to the min_ipar grid point at that iperp. (We also need to set the -// min and max i in the interpolator parameter table.) With this strategy -// we can share the interpolator parameter table across all the iperp -// values, and we don't need to modify the parameter table after the -// initial setup. +// where i is the data array index. However, we don't want to use +// offset, because it has to be supplied in the parameter table as +// an array subscripted by gfn, and so would require changing the +// parameter table for each call on interpolate() (with potentially +// different numbers of gridfns being interpolated). Instead, at each +// iperp we use i = ipar-min_ipar, so the default offset=0 makes the +// subscripting expression zero for ipar = min_ipar. This also makes +// the interpolator's min_i = 0 and max_i be dims-1 (both the defaults), +// so those also don't have to be set in the parameter table. We set +// the interpolator's data coordinate origin to the par coordinate for +// min_ipar, so it correctly maps i --> par. With this strategy we can +// share the interpolator parameter table across all the iperp values, +// and we don't need to modify the parameter table at all after the +// initial setup in our constructor. // class patch_interp |