diff options
Diffstat (limited to 'src/patch/patch_interp.cc')
| -rw-r--r-- | src/patch/patch_interp.cc | 96 |
1 files changed, 34 insertions, 62 deletions
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index 3335897..e21e09c 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -50,18 +50,33 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in, int ghosted_min_gfn_in, int ghosted_max_gfn_in, int interp_handle_in, int interp_par_table_handle_in) - : my_edge_(my_edge_in), - my_patch_(my_edge_in.my_patch()), + : my_patch_(my_edge_in.my_patch()), + my_edge_(my_edge_in), min_iperp_(min_iperp_in), max_iperp_(max_iperp_in), + min_ipar_(my_patch() + .ghost_zone_on_edge(my_edge_in.min_par_adjacent_edge()) + .is_symmetry() + ? my_edge_in.min_ipar_with_corners() + : my_edge_in.min_ipar_without_corners()), + max_ipar_(my_patch() + .ghost_zone_on_edge(my_edge_in.max_par_adjacent_edge()) + .is_symmetry() + ? my_edge_in.max_ipar_with_corners() + : my_edge_in.max_ipar_without_corners()) + min_parindex_used_(min_parindex_used_in), + max_parindex_used_(max_parindex_used_in), + interp_par_(interp_par_in), ghosted_min_gfn_(ghosted_min_gfn_in), ghosted_max_gfn_(ghosted_max_gfn_in), interp_handle_(interp_handle_in), - interp_par_table_handle_(interp_par_table_handle_in), + interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)), interp_coord_origin_(0.0), // computed later interp_coord_delta_(my_patch().delta_ang(my_edge().par_is_rho())), - interp_par_coords_(ghosted_min_gfn_, ghosted_max_gfn_, - min_iperp_, max_iperp_, - min_parindex_in, max_parindex_in) + gridfn_data_ptrs_(ghosted_min_gfn_in, ghosted_max_gfn_in, + min_iperp_in, max_iperp_in), + gridfn_offsets_(ghosted_min_gfn_in, ghosted_max_gfn_in, + min_iperp_in, max_iperp_in), + gridfn_par_stride_(0) // computed later { } @@ -72,67 +87,24 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in, // patch_frontier::~patch_frontier() { -// delete the interpolator objects - for (int gfn = my_patch().min_gfn() ; - gfn <= my_patch().max_gfn() ; - ++gfn) - { - for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) - { - delete (*interpolator_)(gfn,iperp); - } - } - -delete interpolator_; -delete gridfn_buffer_; -delete ripar_map_; +Util_TableDestroy(interp_par_table_handle_); } //***************************************************************************** // -// This function sets up a specified gridfn's interpolation state, i.e. -// it copies the gridfn data to the contiguous-in-ipar buffer and sets -// up any interpolator state +// This function interpolates the specified range of ghosted gridfns +// at all the coordinates specified when we were constructed. +// It stores the interpolated data in result_buffer(gfn, iperp,parindex) . +// It calls error_exit() if the point is out of range or any other +// interpolator error is detected. // -void patch_frontier::setup_interpolation_for_gridfn(int gfn) -{ - for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) - { - // - // copy this row of gridfn data from the (probably non-contiguous) - // gridfn array, into our (contiguous-in-ipar) scratch array - // - for (int ipar = min_ipar() ; ipar <= max_ipar() ; ++ipar) - { - int irho = my_edge(). irho_of_iperp_ipar(iperp, ipar); - int isigma = my_edge().isigma_of_iperp_ipar(iperp, ipar); - (*gridfn_buffer_)(gfn,iperp,ipar) - = my_patch().gridfn(gfn,irho,isigma); - } - - // - // Now we want to set up the interpolator for the just-copied row - // of data. The interpolator interface requires that we pass a - // pointer to the [0] data, where the [] subscripting refers to - // the interpolator's integer coordinate. Since ipar=0 may not - // even exist, we introduce ripar := ipar - min_ipar() and use - // this as the interpolator's integer coordinate, so that [0] - // refers to the start of this ipar row in the gridfn buffer. - // - const fp* row_ptr = & (*gridfn_buffer_)(gfn,iperp,min_ipar()); - interpolator_t_* I = (*interpolator_)(gfn,iperp); - I->setup_data_y(row_ptr); - } -} - -//****************************************************************************** - -// -// This function does the actual interpolation of data for a specified -// gridfn to a specified (iperp,par) point. It calls error_exit() -// if the point is out of range or any other interpolator error is found. -// -fp patch_frontier::interpolate(int gfn, int iperp, fp par) const +fp patch_frontier::interpolate(int ghosted_min_gfn, ghosted_max_gfn, + jtutil::array3d<fp>& result_buffer) const { const int N_dims = 1; + +const int status + = CCTK_InterpLocalUniform(N_dims, + interp_handle_, + interp_par_table_handle_, |