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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
|
/*@@
@file molecule_posn.c
@date 23 Oct 2001
@author Jonathan Thornburg <jthorn@aei.mpg.de>
@desc Worker function to compute molecule positions for interpolator.
@enddesc
@version $Id$
@@*/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#ifndef AEILOCALINTERP_STANDALONE_TEST
#include "cctk.h"
#endif
#include "InterpLocalUniform.h"
/* the rcs ID and its dummy function to use it */
static const char *rcsid = "$Header$";
#ifndef AEILOCALINTERP_STANDALONE_TEST
CCTK_FILEVERSION(AEIThorns_AEILocalInterp_src_molecule_posn_c)
#endif
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
/*@@
@routine AEILocalInterp_molecule_posn
@date 23 Oct 2001
@author Jonathan Thornburg <jthorn@aei.mpg.de>
@desc
Given a uniformly-spaced grid, this function computes molecule
positions for an interpolation or similar operation:
(1) It computes the closest grid point to the input coordinate
(2) It does the slightly tricky odd/even computation to decide where
to center the molecule
(3) It checks for the molecule falling off the edge of the grid,
and off-centers it as appropriate.
For (1), we assume that grid points have floating-point
coordinates (we refer to these as "x") which are an arbitrary
linear function of the integer grid coordinates (we refer to
these as "i"), x = grid_origin + i*grid_delta.
For (2), suppose we have a data array indexed by [i], from which we
wish to select an N-point molecule [i_lo, i_hi] centered as close as
possible to the point x . The problem is just how to choose the
centering. The following diagram illustrates the issues:
N i_lo i_hi [i-3] [i-2] [i-1] [i] [i+1] [i+2] [i+3] [i+4]
- ---- ---- ----- ----- ----- ----- ----- ----- ----- -----
2 i i+1 *xxxxxx*
3 i-1 i+1 * xxx*xxx *
4 i-1 i+2 * *xxxxxx* *
5 i-2 i+2 * * xxx*xxx * *
6 i-2 i+3 * * *xxxxxx* * *
7 i-3 i+3 * * * xxx*xxx * * *
8 i-3 i+4 * * * *xxxxxx* * * *
The diagram shows the range of x values relative to the grid,
for which we'll choose each centering. Thus if N is even, the
centering is based on which grid zone contains x (take the floor
of the floating-point i coordinate), while if N is odd, the centering
is based on which grid point is closest to x (round the floating-
-point i coordinate to the nearest integer).
It's also convenient to introduce the integer molecule coordinate m;
this is the integer grid coordinate i relative to that of the molecule
center, i.e. the above diagram has columns labelled with [i+m].
@enddesc
@hdate 27.Jan.2003
@hauthor Jonathan Thornburg <jthorn@aei.mpg.de>
@hdesc Complete rewrite: now supports
@var boundary_off_centering_tolerance_min,
@var boundary_off_centering_tolerance_max,
@var boundary_extrapolation_tolerance_min,
@var boundary_extrapolation_tolerance_max,
also change to returning status code,
also drop returning @var min_m and @var max_m
since they were never used.
@endhdesc
@var grid_origin
@vdesc The floating-point coordinate x of the grid point i=0.
@vtype fp grid_origin
@endvar
@var grid_delta
@vdesc The grid spacing (in the floating-point coordinate x).
@vtype fp grid_delta
@endvar
@var grid_i_min
@vdesc The minimum integer coordinate i of the grid.
@vtype int grid_i_min
@endvar
@var grid_i_max
@vdesc The maximum integer coordinate i of the grid.
@vtype int grid_i_max
@endvar
@var molecule_size
@vdesc The size (number of points) of the molecule.
@vtype int molecule_size
@endvar
@var boundary_off_centering_tolerance_min,
boundary_off_centering_tolerance_max
@vdesc Specifies how many grid spacings the interpolation point
may be beyond the default-centering region before being
declared "out of range" on the {minimum,maximum} boundaries
of the grid respectively.
@vtype fp boundary_off_centering_tolerance_min,
boundary_off_centering_tolerance_max
@endvar
@var boundary_extrapolation_tolerance_min,
boundary_extrapolation_tolerance_max
@vdesc Specifies how many grid spacings the interpolation point
may be beyond the grid boundary before being declared
"out of range" on the {minimum,maximum} boundaries
of the grid respectively.
@vtype fp boundary_extrapolation_tolerance_min,
boundary_extrapolation_tolerance_max
@endvar
@var x
@vdesc The floating-point coordinate of the interpolation point.
@vtype fp x
@endvar
@var i_center
@vdesc A pointer to an value where this function should
store the integer coordinate of the molecule center,
or NULL to skip storing this
@vtype int *i_center
@vio pointer to out
@endvar
@var x_rel
@vdesc A pointer to where this function should store the
interpolation point's x coordinate relative to the
molecule center, measured in units of the grid spacing;
or NULL to skip storing this.
@vtype fp *x_rel
@vio pointer to out
@endvar
@returntype int
@returndesc
This function returns 0 if the interpolation point is "in range",
or one of the (negative) error codes defined in "InterpLocalUniform.h"
if the interpolation point is "out of range":
MOLECULE_POSN_ERROR_X_LT_MIN
if x is out-of-range on the min end of the grid
(i.e. x < the minimum allowable x)
MOLECULE_POSN_ERROR_X_GT_MAX
if x is out-of-range on the max end of the grid
(i.e. x > the maximum allowable x)
MOLECULE_POSN_ERROR_GRID_TINY
if the grid is smaller than the molecule
@endreturndesc
@@*/
int AEILocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
int grid_i_min, int grid_i_max,
int molecule_size,
fp boundary_off_centering_tolerance_min,
fp boundary_off_centering_tolerance_max,
fp boundary_extrapolation_tolerance_min,
fp boundary_extrapolation_tolerance_max,
fp x,
int* i_center, fp* x_rel)
{
/* molecule radia (inherently positive) in +/- directions */
const int mr_plus = (molecule_size >> 1);
const int mr_minus = molecule_size - mr_plus - 1;
/* range of x_rel for which this molecule is centered, cf. diagram in header comment */
const fp centered_min_x_rel = IS_EVEN(molecule_size) ? 0.0 : -0.5;
const fp centered_max_x_rel = IS_EVEN(molecule_size) ? 1.0 : +0.5;
/* range of i where we *could* center the molecule, as floating-point numbers */
const fp fp_centered_min_possible_i = grid_i_min + mr_minus + centered_min_x_rel;
const fp fp_centered_max_possible_i = grid_i_max - mr_plus + centered_max_x_rel;
/* integer coordinate i of interpolation point, as a floating-point number */
const fp fp_i = (x - grid_origin) / grid_delta;
/* is the molecule larger than the grid? */
if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max))
then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/
/* is interpolation point x beyond the extrapolation tolerance? */
if (fp_i < (fp)grid_i_min - boundary_extrapolation_tolerance_min)
then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/
if (fp_i > (fp)grid_i_max + boundary_extrapolation_tolerance_max)
then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/
/* is interpolation point x beyond the off-centering tolerance? */
if (fp_i < fp_centered_min_possible_i - boundary_off_centering_tolerance_min)
then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/
if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max)
then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/
/* now choose the actual molecule position/centering: */
/* first set up a centered molecule */
{
fp fp_i_center = IS_EVEN(molecule_size) /* ... as a floating-point number */
? floor(fp_i)
: JT_ROUND(fp_i);
int int_i_center = (int) fp_i_center; /* ... as an integer */
/* then clamp the molecule at the grid boundaries */
if (int_i_center < grid_i_min) int_i_center = grid_i_min;
if (int_i_center < grid_i_min + mr_minus)
then {
int_i_center = grid_i_min + mr_minus;
fp_i_center = (fp) int_i_center;
}
if (int_i_center > grid_i_max) int_i_center = grid_i_max;
if (int_i_center > grid_i_max - mr_plus)
then {
int_i_center = grid_i_max - mr_plus;
fp_i_center = (fp) int_i_center;
}
/* store the results */
if (i_center != NULL)
then *i_center = int_i_center;
if (x_rel != NULL)
then *x_rel = fp_i - fp_i_center;
return 0; /*** NORMAL RETURN ***/
}
}
|