mantaflow  0.10
A framework for fluid simulation
interpol.h
Go to the documentation of this file.
1 
2 /******************************************************************************
3  *
4  * MantaFlow fluid solver framework
5  * Copyright 2011 Tobias Pfaff, Nils Thuerey
6  *
7  * This program is free software, distributed under the terms of the
8  * GNU General Public License (GPL)
9  * http://www.gnu.org/licenses
10  *
11  * Helper functions for interpolation
12  *
13  ******************************************************************************/
14 
15 #ifndef _INTERPOL_H
16 #define _INTERPOL_H
17 
18 #include "vectorbase.h"
19 
20 // Grid values are stored at i+0.5, j+0.5, k+0.5
21 // MAC grid values are stored at i,j+0.5,k+0.5 (for x) ...
22 
23 namespace Manta {
24 
25 inline Vec3 fdTangent(const Vec3& p0, const Vec3& p1, const Vec3& p2) {
26  return 0.5*(getNormalized(p2-p1) + getNormalized(p1-p0));
27 }
28 
29 inline Vec3 crTangent(const Vec3& p0, const Vec3& p1, const Vec3& p2) {
30  return 0.5*(p2-p0);
31 }
32 
33 inline Vec3 hermiteSpline(const Vec3& p0, const Vec3& p1, const Vec3& m0, const Vec3& m1, Real t) {
34  const Real t2=t*t, t3=t2*t;
35  return (2.0*t3 - 3.0*t2 + 1.0)*p0 + (t3 - 2.0*t2 + t)*m0 + (-2.0*t3 + 3.0*t2)*p1 + (t3 - t2)*m1;
36 }
37 
38 static inline void checkIndexInterpol(const Vec3i& size, IndexInt idx) {
39  if (idx<0 || idx > (IndexInt)size.x * size.y * size.z) {
40  std::ostringstream s;
41  s << "Grid interpol dim " << size << " : index " << idx << " out of bound ";
42  errMsg(s.str());
43  }
44 }
45 
46 
47 // ----------------------------------------------------------------------
48 // Grid interpolators
49 // ----------------------------------------------------------------------
50 
51 #define BUILD_INDEX \
52  Real px=pos.x-0.5f, py=pos.y-0.5f, pz=pos.z-0.5f; \
53  int xi = (int)px; \
54  int yi = (int)py; \
55  int zi = (int)pz; \
56  Real s1 = px-(Real)xi, s0 = 1.-s1; \
57  Real t1 = py-(Real)yi, t0 = 1.-t1; \
58  Real f1 = pz-(Real)zi, f0 = 1.-f1; \
59  /* clamp to border */ \
60  if (px < 0.) { xi = 0; s0 = 1.0; s1 = 0.0; } \
61  if (py < 0.) { yi = 0; t0 = 1.0; t1 = 0.0; } \
62  if (pz < 0.) { zi = 0; f0 = 1.0; f1 = 0.0; } \
63  if (xi >= size.x-1) { xi = size.x-2; s0 = 0.0; s1 = 1.0; } \
64  if (yi >= size.y-1) { yi = size.y-2; t0 = 0.0; t1 = 1.0; } \
65  if (size.z>1) { if (zi >= size.z-1) { zi = size.z-2; f0 = 0.0; f1 = 1.0; } } \
66  const int X = 1; \
67  const int Y = size.x;
68 
69 template <class T>
70 inline T interpol(const T* data, const Vec3i& size, const int Z, const Vec3& pos) {
71  BUILD_INDEX
72  IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi;
73  DEBUG_ONLY(checkIndexInterpol(size,idx)); DEBUG_ONLY(checkIndexInterpol(size,idx+X+Y+Z));
74 
75  return ((data[idx] *t0 + data[idx+Y] *t1) * s0
76  + (data[idx+X] *t0 + data[idx+X+Y] *t1) * s1) * f0
77  +((data[idx+Z] *t0 + data[idx+Y+Z] *t1) * s0
78  + (data[idx+X+Z]*t0 + data[idx+X+Y+Z]*t1) * s1) * f1;
79 }
80 
81 template <int c>
82 inline Real interpolComponent(const Vec3* data, const Vec3i& size, const int Z, const Vec3& pos) {
83  BUILD_INDEX
84  IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi;
85  DEBUG_ONLY(checkIndexInterpol(size,idx)); DEBUG_ONLY(checkIndexInterpol(size,idx+X+Y+Z));
86 
87  return ((data[idx][c] *t0 + data[idx+Y][c] *t1) * s0
88  + (data[idx+X][c] *t0 + data[idx+X+Y][c] *t1) * s1) * f0
89  +((data[idx+Z][c] *t0 + data[idx+Y+Z][c] *t1) * s0
90  + (data[idx+X+Z][c]*t0 + data[idx+X+Y+Z][c]*t1) * s1) * f1;
91 }
92 
93 template<class T>
94 inline void setInterpol(T* data, const Vec3i& size, const int Z, const Vec3& pos, const T& v, Real* sumBuffer)
95 {
96  BUILD_INDEX
97  IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi;
98  DEBUG_ONLY(checkIndexInterpol(size,idx)); DEBUG_ONLY(checkIndexInterpol(size,idx+X+Y+Z));
99 
100  T* ref = &data[idx];
101  Real* sum = &sumBuffer[idx];
102  Real s0f0=s0*f0, s1f0=s1*f0, s0f1=s0*f1, s1f1=s1*f1;
103  Real w0 = t0*s0f0, wx = t0*s1f0, wy = t1*s0f0, wxy = t1*s1f0;
104  Real wz = t0*s0f1, wxz = t0*s1f1, wyz = t1*s0f1, wxyz = t1*s1f1;
105 
106  sum[Z] += wz; sum[X+Z] += wxz; sum[Y+Z] += wyz; sum[X+Y+Z] += wxyz;
107  ref[Z] += wz*v; ref[X+Z] += wxz*v; ref[Y+Z] += wyz*v; ref[X+Y+Z] += wxyz*v;
108  sum[0] += w0; sum[X] += wx; sum[Y] += wy; sum[X+Y] += wxy;
109  ref[0] += w0*v; ref[X] += wx*v; ref[Y] += wy*v; ref[X+Y] += wxy*v;
110 }
111 
112 
113 #define BUILD_INDEX_SHIFT \
114  BUILD_INDEX \
115  /* shifted coords */ \
116  int s_xi = (int)pos.x, s_yi = (int)pos.y, s_zi = (int)pos.z; \
117  Real s_s1 = pos.x-(Real)s_xi, s_s0 = 1.-s_s1; \
118  Real s_t1 = pos.y-(Real)s_yi, s_t0 = 1.-s_t1; \
119  Real s_f1 = pos.z-(Real)s_zi, s_f0 = 1.-s_f1; \
120  /* clamp to border */ \
121  if (pos.x < 0) { s_xi = 0; s_s0 = 1.0; s_s1 = 0.0; } \
122  if (pos.y < 0) { s_yi = 0; s_t0 = 1.0; s_t1 = 0.0; } \
123  if (pos.z < 0) { s_zi = 0; s_f0 = 1.0; s_f1 = 0.0; } \
124  if (s_xi >= size.x-1) { s_xi = size.x-2; s_s0 = 0.0; s_s1 = 1.0; } \
125  if (s_yi >= size.y-1) { s_yi = size.y-2; s_t0 = 0.0; s_t1 = 1.0; } \
126  if (size.z>1) { if (s_zi >= size.z-1) { s_zi = size.z-2; s_f0 = 0.0; s_f1 = 1.0; } }
127 
128 inline Vec3 interpolMAC(const Vec3* data, const Vec3i& size, const int Z, const Vec3& pos)
129 {
130  BUILD_INDEX_SHIFT;
131  DEBUG_ONLY(checkIndexInterpol(size, (zi*(IndexInt)size.y + yi)*(IndexInt)size.x + xi));
132  DEBUG_ONLY(checkIndexInterpol(size, (s_zi*(IndexInt)size.y + s_yi)*(IndexInt)size.x + s_xi + X + Y + Z));
133 
134  // process individual components
135  Vec3 ret(0.);
136  { // X
137  const Vec3* ref = &data[((zi*size.y+yi)*size.x+s_xi)];
138  ret.x = f0 * ((ref[0].x *t0 + ref[Y].x *t1 )*s_s0 +
139  (ref[X].x *t0 + ref[X+Y].x *t1 )*s_s1) +
140  f1 * ((ref[Z].x *t0 + ref[Z+Y].x *t1 )*s_s0 +
141  (ref[X+Z].x*t0 + ref[X+Y+Z].x*t1 )*s_s1 );
142  }
143  { // Y
144  const Vec3* ref = &data[((zi*size.y+s_yi)*size.x+xi)];
145  ret.y = f0 * ((ref[0].y *s_t0 + ref[Y].y *s_t1 )*s0 +
146  (ref[X].y *s_t0 + ref[X+Y].y *s_t1 )*s1) +
147  f1 * ((ref[Z].y *s_t0 + ref[Z+Y].y *s_t1 )*s0 +
148  (ref[X+Z].y*s_t0 + ref[X+Y+Z].y*s_t1 )*s1 );
149  }
150  { // Z
151  const Vec3* ref = &data[((s_zi*size.y+yi)*size.x+xi)];
152  ret.z = s_f0 * ((ref[0].z *t0 + ref[Y].z *t1 )*s0 +
153  (ref[X].z *t0 + ref[X+Y].z *t1 )*s1) +
154  s_f1 * ((ref[Z].z *t0 + ref[Z+Y].z *t1 )*s0 +
155  (ref[X+Z].z*t0 + ref[X+Y+Z].z*t1 )*s1 );
156  }
157  return ret;
158 }
159 
160 inline void setInterpolMAC(Vec3* data, const Vec3i& size, const int Z, const Vec3& pos, const Vec3& val, Vec3* sumBuffer)
161 {
162  BUILD_INDEX_SHIFT;
163  DEBUG_ONLY(checkIndexInterpol(size, (zi*(IndexInt)size.y + yi)*(IndexInt)size.x + xi));
164  DEBUG_ONLY(checkIndexInterpol(size, (s_zi*(IndexInt)size.y + s_yi)*(IndexInt)size.x + s_xi + X + Y + Z));
165 
166  // process individual components
167  { // X
168  const IndexInt idx = (IndexInt)(zi*size.y+yi)*size.x+s_xi;
169  Vec3 *ref = &data[idx], *sum = &sumBuffer[idx];
170  Real s0f0=s_s0*f0, s1f0=s_s1*f0, s0f1=s_s0*f1, s1f1=s_s1*f1;
171  Real w0 = t0*s0f0, wx = t0*s1f0, wy = t1*s0f0, wxy = t1*s1f0;
172  Real wz = t0*s0f1, wxz = t0*s1f1, wyz = t1*s0f1, wxyz = t1*s1f1;
173 
174  sum[Z].x += wz; sum[X+Z].x += wxz; sum[Y+Z].x += wyz; sum[X+Y+Z].x += wxyz;
175  ref[Z].x += wz*val.x; ref[X+Z].x += wxz*val.x; ref[Y+Z].x += wyz*val.x; ref[X+Y+Z].x += wxyz*val.x;
176  sum[0].x += w0; sum[X].x += wx; sum[Y].x += wy; sum[X+Y].x += wxy;
177  ref[0].x += w0*val.x; ref[X].x += wx*val.x; ref[Y].x += wy*val.x; ref[X+Y].x += wxy*val.x;
178  }
179  { // Y
180  const IndexInt idx = (IndexInt)(zi*size.y+s_yi)*size.x+xi;
181  Vec3 *ref = &data[idx], *sum = &sumBuffer[idx];
182  Real s0f0=s0*f0, s1f0=s1*f0, s0f1=s0*f1, s1f1=s1*f1;
183  Real w0 = s_t0*s0f0, wx = s_t0*s1f0, wy = s_t1*s0f0, wxy = s_t1*s1f0;
184  Real wz = s_t0*s0f1, wxz = s_t0*s1f1, wyz = s_t1*s0f1, wxyz = s_t1*s1f1;
185 
186  sum[Z].y += wz; sum[X+Z].y += wxz; sum[Y+Z].y += wyz; sum[X+Y+Z].y += wxyz;
187  ref[Z].y += wz*val.y; ref[X+Z].y += wxz*val.y; ref[Y+Z].y += wyz*val.y; ref[X+Y+Z].y += wxyz*val.y;
188  sum[0].y += w0; sum[X].y += wx; sum[Y].y += wy; sum[X+Y].y += wxy;
189  ref[0].y += w0*val.y; ref[X].y += wx*val.y; ref[Y].y += wy*val.y; ref[X+Y].y += wxy*val.y;
190  }
191  { // Z
192  const IndexInt idx = (IndexInt)(s_zi*size.y+yi)*size.x+xi;
193  Vec3 *ref = &data[idx], *sum = &sumBuffer[idx];
194  Real s0f0=s0*s_f0, s1f0=s1*s_f0, s0f1=s0*s_f1, s1f1=s1*s_f1;
195  Real w0 = t0*s0f0, wx = t0*s1f0, wy = t1*s0f0, wxy = t1*s1f0;
196  Real wz = t0*s0f1, wxz = t0*s1f1, wyz = t1*s0f1, wxyz = t1*s1f1;
197 
198  sum[0].z += w0; sum[X].z += wx; sum[Y].z += wy; sum[X+Y].z += wxy;
199  sum[Z].z += wz; sum[X+Z].z += wxz; sum[Y+Z].z += wyz; sum[X+Y+Z].z += wxyz;
200  ref[0].z += w0*val.z; ref[X].z += wx*val.z; ref[Y].z += wy*val.z; ref[X+Y].z += wxy*val.z;
201  ref[Z].z += wz*val.z; ref[X+Z].z += wxz*val.z; ref[Y+Z].z += wyz*val.z; ref[X+Y+Z].z += wxyz*val.z;
202  }
203 }
204 
205 #undef BUILD_INDEX
206 #undef BUILD_INDEX_SHIFT
207 
208 } //namespace
209 
210 #endif
211 
212 
Manta
Definition: commonkernels.h:22
Manta::Manta::getNormalized
Vector3D< S > getNormalized(const Vector3D< S > &v)
Returns a normalized vector.
Definition: randomstream.h:392
Manta::Manta::Vec3
Vector3D< Real > Vec3
3D vector class of type Real (typically float)
Definition: randomstream.h:545
Generated on Tue Sep 19 2017 17:49:04 for mantaflow by   doxygen 1.8.13