31 namespace interpolation {
34 : x_m(2, 0), y_m(2, 0), z_m(2, 0), xSize_m(0), ySize_m(0), zSize_m(0),
35 maps_m(0), constantSpacing_m(false) {
41 int ySize,
const double *y,
42 int zSize,
const double *z)
43 : x_m(x, x+xSize), y_m(y, y+ySize), z_m(z, z+zSize),
44 xSize_m(xSize), ySize_m(ySize), zSize_m(zSize),
45 maps_m(), constantSpacing_m(false) {
49 "ThreeDGrid::ThreeDGrid(...)",
50 "3D Grid must be at least 2x2x2 grid"
56 std::vector<double> y,
57 std::vector<double> z)
58 : x_m(x), y_m(y), z_m(z),
59 xSize_m(x.size()), ySize_m(y.size()), zSize_m(z.size()),
60 maps_m(), constantSpacing_m(false) {
64 "ThreeDGrid::ThreeDGrid(...)",
65 "3D Grid must be at least 2x2x2 grid"
71 double minX,
double minY,
double minZ,
72 int numberOfXCoords,
int numberOfYCoords,
74 : x_m(numberOfXCoords), y_m(numberOfYCoords), z_m(numberOfZCoords),
75 xSize_m(numberOfXCoords), ySize_m(numberOfYCoords),
76 zSize_m(numberOfZCoords), maps_m(), constantSpacing_m(true) {
77 for (
int i = 0; i < numberOfXCoords; i++)
79 for (
int j = 0; j < numberOfYCoords; j++)
81 for (
int k = 0; k < numberOfZCoords; k++)
96 std::vector<int>
end(3, 1);
102 double * position)
const {
109 std::vector<double> new_x(
x_m.size()-1);
110 std::vector<double> new_y(
y_m.size()-1);
111 std::vector<double> new_z(
z_m.size()-1);
112 for (
size_t i = 0; i <
x_m.size()-1; ++i) {
113 new_x[i] = (
x_m[i]+
x_m[i+1])/2.;
115 for (
size_t i = 0; i <
y_m.size()-1; ++i) {
116 new_y[i] = (
y_m[i]+
y_m[i+1])/2.;
118 for (
size_t i = 0; i <
z_m.size()-1; ++i) {
119 new_z[i] = (
z_m[i]+
z_m[i+1])/2.;
127 return subEquals(lhs, -difference);
129 lhs.
state_m[0] += difference/(ySize_m*zSize_m);
130 difference = difference%(ySize_m*zSize_m);
131 lhs.
state_m[1] += difference/(zSize_m);
132 lhs.
state_m[2] += difference%(zSize_m);
134 if (lhs.
state_m[1] > ySize_m) {
138 if (lhs.
state_m[2] > zSize_m) {
149 return addEquals(lhs, -difference);
151 lhs.
state_m[0] -= difference/(ySize_m*zSize_m);
152 difference = difference%(ySize_m*zSize_m);
153 lhs.
state_m[1] -= difference/(zSize_m);
154 lhs.
state_m[2] -= difference%(zSize_m);
175 if (x >= vec.back()) {
176 index = vec.size()-1;
180 size_t xUpper = vec.size()-1;
181 while (xUpper - xLower > 1) {
182 index = (xUpper+xLower)/2;
183 if (x >= vec[index]) {
220 }
else if (lhs.
state_m[2] == 1) {
267 for (
unsigned int i = 0; i <
x_m.size()-1; i++)
272 for (
unsigned int i = 0; i <
y_m.size()-1; i++)
277 for (
unsigned int i = 0; i <
z_m.size()-1; i++)
285 std::vector<int> index(3);
287 position[1], index[1],
288 position[2], index[2]);
289 if (index[0] <
xSize_m-1 && index[0] >= 0)
290 index[0] += (2*(position[0] -
x_m[index[0]]) >
291 x_m[index[0]+1]-
x_m[index[0]] ? 2 : 1);
294 if (index[1] <
ySize_m-1 && index[1] >= 0)
295 index[1] += (2*(position[1] -
y_m[index[1]]) >
296 y_m[index[1]+1]-
y_m[index[1]] ? 2 : 1);
299 if (index[2] <
zSize_m-1 && index[2] >= 0)
300 index[2] += (2*(position[2] -
z_m[index[2]]) >
301 z_m[index[2]+1]-
z_m[index[2]] ? 2 : 1);
Mesh::Iterator getNearest(const double *position) const
Mesh::Iterator end() const
constexpr double e
The value of .
std::vector< int > state_m
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
std::vector< double > z_m
void remove(VectorMap *map)
virtual Mesh::Iterator & subOne(Mesh::Iterator &lhs) const
std::vector< double > y_m
virtual Mesh::Iterator & subEquals(Mesh::Iterator &lhs, int difference) const
virtual void getPosition(const Mesh::Iterator &it, double *position) const
void setConstantSpacing()
Mesh::Iterator begin() const
virtual bool isGreater(const Mesh::Iterator &lhs, const Mesh::Iterator &rhs) const
static void vectorLowerBound(std::vector< double > vec, double x, int &index)
virtual Mesh::Iterator & addEquals(Mesh::Iterator &lhs, int difference) const
const T * find(const T table[], const std::string &name)
Look up name.
Base class for meshing routines.
std::vector< VectorMap * > maps_m
virtual Mesh::Iterator & addOne(Mesh::Iterator &lhs) const
std::string::iterator iterator
void lowerBound(const double &x, int &xIndex, const double &y, int &yIndex, const double &z, int &zIndex) const
std::vector< double > x_m