aboutsummaryrefslogtreecommitdiff
path: root/src/patch/patch_interp.cc
blob: e21e09c544d8b041b016dededae48bcad3208661 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
// patch_frontier.cc -- interpolate gridfn data for another patch's ghost zone
// $Id$

//
// patch_frontier::patch_frontier
// patch_frontier::~patch_frontier
// patch_frontier::setup_interpolation_for_gridfn
// patch_frontier::interpolate
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
#include "jt/interpolate.hh"

#include "fp.hh"
#include "coords.hh"
#include "grid.hh"
#include "fd_grid.hh"
#include "patch.hh"
#include "patch_edge.hh"
#include "ghost_zone.hh"
#include "patch_frontier.hh"

//*****************************************************************************
//*****************************************************************************
// patch_frontier - interpolation to get data for another patch's ghost zone
//*****************************************************************************
//*****************************************************************************

//
// This function constructs an  patch_frontier  object.
//
// Note this requires that the adjacent-side  ghost_zone  objects already
// exist, since they're used in determining the ipar range over which the
// interpolation will be done.
//
patch_frontier::patch_frontier(const patch_edge& my_edge_in,
			       int min_iperp_in, int max_iperp_in,
			       int min_parindex_in, int max_parindex_in,
			       const jtutil::array1d<fp>& min_parindex_used_in,
			       const jtutil::array1d<fp>& max_parindex_used_in,
			       const jtutil::array2d<fp>& interp_par_in,
			       int ghosted_min_gfn_in, int ghosted_max_gfn_in,
			       int interp_handle_in,
			       int interp_par_table_handle_in)
	: 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_(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())),
	  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
{
}

//*****************************************************************************

//
// This function destroys an  patch_frontier  object
//
patch_frontier::~patch_frontier()
{
Util_TableDestroy(interp_par_table_handle_);
}

//*****************************************************************************

//
// 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.
//
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_,