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
|
/*@@
@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 <limits.h>
#include <stdlib.h>
#include <stdio.h>
#ifndef LOCALINTERP_STANDALONE_BUILD
#include "cctk.h"
#endif
#include "InterpLocalUniform.h"
/* the rcs ID and its dummy function to use it */
static const char *rcsid = "$Header$";
#ifndef LOCALINTERP_STANDALONE_BUILD
CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_molecule_posn_c)
#endif
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
/*@@
@routine LocalInterp_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
@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 out_of_range_tolerance_min, out_of_range_tolerance_max
@vdesc Specifies how out-of-range interpolation points should
be handled for the {minimum,maximum} ends of the grid
respectively.
If out_of_range_tolerance >= 0.0,
then an interpolation point is considered to be
"out of range" if and only if the interpolation
point is > out_of_range_tolerance * grid_delta
outside the grid in any coordinate.
If out_of_range_tolerance == -1.0,
then an interpolation point is considered to be
"out of range" if and only if a centered molecule
would require data from outside the grid.
Other values of out_of_range_tolerance are illegal.
@vtype fp out_of_range_tolerance_min, out_of_range_tolerance_max
@endvar
@var x
@vdesc The floating-point coordinate x at which we would like to center
the molecule.
@vtype fp x
@endvar
@var x_rel
@vdesc A pointer to a floating-point value where this function should
store the 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
@var molecule_m_min, molecule_m_max
@vdesc A pointer to an int where this function should store the
{minimum,maximum} molecule coordinate m of the molecule;
or NULL to skip storing this.
@vtype int *molecule_m_min, *molecule_m_max
@vio pointer to out
@endvar
@returntype int
@returndesc
This function returns the integer coordinate i_center at which
the molecule should be centered, or one of the error codes
(from <limits.h>)
INT_MIN if x is out-of-range on the min end of the grid
(i.e. x < the minimum allowable x)
INT_MAX if x is out-of-range on the max end of the grid
(i.e. x > the maximum allowable x)
@endreturndesc
@@*/
int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
int grid_i_min, int grid_i_max,
int molecule_size,
fp out_of_range_tolerance_min,
fp out_of_range_tolerance_max,
fp x,
fp *x_rel,
int *molecule_m_min, int *molecule_m_max)
{
/* min/max molecule coordinates */
const int m_max = (molecule_size >> 1); /* a positive number */
const int m_min = m_max - molecule_size + 1; /* a negative number */
/* integer coordinate i of point x, as a floating-point number */
const fp fp_i = (x - grid_origin) / grid_delta;
/* is point x out-of-range? */
if (out_of_range_tolerance_min >= 0.0)
then {
const fp fp_effective_grid_i_min
= ((fp) grid_i_min) - out_of_range_tolerance_min;
if (fp_i < fp_effective_grid_i_min)
then return INT_MIN; /*** ERROR RETURN ***/
}
if (out_of_range_tolerance_max >= 0.0)
then {
const fp fp_effective_grid_i_max
= ((fp) grid_i_max) + out_of_range_tolerance_max;
if (fp_i > fp_effective_grid_i_max)
then return INT_MAX; /*** ERROR RETURN ***/
}
/* integer coordinate i where we will center the molecule */
/* (see diagram in header comment for explanation) */
/* ... as a floating-point number */
{
const fp fp_i_center = is_even(molecule_size) ? floor(fp_i) : jt_round(fp_i);
/* ... as an integer */
int i_center = (int) fp_i_center;
/* min/max integer coordinates in molecule assuming it's fully centered */
const int centered_i_min = i_center + m_min;
const int centered_i_max = i_center + m_max;
/* check if off-centered molecules are forbidden? */
if ((out_of_range_tolerance_min == -1.0) && (centered_i_min < grid_i_min))
then return INT_MIN; /*** ERROR RETURN ***/
if ((out_of_range_tolerance_max == -1.0) && (centered_i_max > grid_i_max))
then return INT_MAX; /*** ERROR RETURN ***/
/* off-center as needed if we're close to the edge of the grid */
{
int shift = 0;
if (centered_i_min < grid_i_min)
then shift = grid_i_min - centered_i_min; /* a positive number */
if (centered_i_max > grid_i_max)
then shift = grid_i_max - centered_i_max; /* a negative number */
i_center += shift;
/* store the results */
if (x_rel != NULL)
then *x_rel = fp_i - ((fp) i_center);
if (molecule_m_min != NULL)
then *molecule_m_min = m_min;
if (molecule_m_max != NULL)
then *molecule_m_max = m_max;
return i_center;
}
}
}
|
