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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
// 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 interpolator_order_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_edge_in.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_edge_in.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())
{
const linear_map<fp> &patch_par_map
= my_patch().ang_map(my_edge().par_is_rho());
// map defining subrange of patch's par map, from which we'll use data
// ... n.b. integer coordinate is ripar, not ipar;
// see setup_interpolation_for_gridfn() for discussion of this
ripar_map_ = new linear_map<fp>(min_ripar(), max_ripar(),
patch_par_map.min_fp(),
patch_par_map.delta_fp(),
patch_par_map.max_fp());
gridfn_buffer_
= new array3d<fp>(my_patch().min_gfn(), my_patch().max_gfn(),
min_iperp(), max_iperp(),
min_ipar(), max_ipar());
interpolator_
= new array2d<interpolator_t_ *>
(my_patch().min_gfn(), my_patch().max_gfn(),
min_iperp(), max_iperp());
// create the interpolator objects
// and set up the coordinates for the interpolators
for (int gfn = my_patch().min_gfn() ;
gfn <= my_patch().max_gfn() ;
++gfn)
{
for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
{
interpolator_t_* I = new interpolator_t_(interpolator_order_in);
(*interpolator_)(gfn,iperp) = I;
I->setup_data_x(*ripar_map_);
}
}
}
//*****************************************************************************
//
// This function destroys an patch_frontier object
//
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_;
}
//*****************************************************************************
//
// This function sets up the interpolation state for a specified gridfn,
// i.e. it copies the gridfn data to the contiguous-in-ipar buffer and
// sets up any interpolator state
//
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.
//
fp patch_frontier::interpolate(int gfn, int iperp, fp par) const
{
interpolator_t_* I = (*interpolator_)(gfn,iperp);
I->setup_interp_x(par, 'o');
return I->interpolate();
}
|
