mantaflow  0.10
A framework for fluid simulation
noisefield.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  * Wavelet noise field
12  *
13  ******************************************************************************/
14 
15 #ifndef _NOISEFIELD_H_
16 #define _NOISEFIELD_H_
17 
18 #include "vectorbase.h"
19 #include "manta.h"
20 #include "grid.h"
21 
22 namespace Manta {
23 
24 #define NOISE_TILE_SIZE 128
25 
26 // wrapper for a parametrized field of wavelet noise
28 PYTHON class WaveletNoiseField : public PbClass {
29  public:
30  PYTHON() WaveletNoiseField( FluidSolver* parent, int fixedSeed=-1 , int loadFromFile=false );
31  ~WaveletNoiseField() {
32  if(mNoiseTile) { delete mNoiseTile; mNoiseTile=NULL; }
33  };
34 
36  inline Real evaluate(Vec3 pos, int tile=0);
38  inline Vec3 evaluateVec(Vec3 pos, int tile=0);
40  inline Vec3 evaluateCurl(Vec3 pos);
41 
43  Real* data() { return mNoiseTile; }
44 
46  static void computeCoefficients(Grid<Real>& input, Grid<Real>& tempIn1, Grid<Real>& tempIn2);
47 
48  // helper
49  std::string toString();
50 
51  // texcoord position and scale
52  PYTHON(name=posOffset) Vec3 mPosOffset;
53  PYTHON(name=posScale) Vec3 mPosScale;
54  // value offset & scale
55  PYTHON(name=valOffset) Real mValOffset;
56  PYTHON(name=valScale) Real mValScale;
57  // clamp? (default 0-1)
58  PYTHON(name=clamp) bool mClamp;
59  PYTHON(name=clampNeg) Real mClampNeg;
60  PYTHON(name=clampPos) Real mClampPos;
61  // animated over time
62  PYTHON(name=timeAnim) Real mTimeAnim;
63 
64  protected:
65  // noise evaluation functions
66  static inline Real WNoiseDx (const Vec3& p, Real *data);
67  static inline Vec3 WNoiseVec(const Vec3& p, Real *data);
68  static inline Real WNoise (const Vec3& p, Real *data);
69 
70  // helpers for tile generation , for periodic 128 grids only
71  static void downsample(Real *from, Real *to, int n, int stride);
72  static void upsample (Real *from, Real *to, int n, int stride);
73 
74  // for grids with arbitrary sizes, and neumann boundary conditions
75  static void downsampleNeumann(const Real *from, Real *to, int n, int stride);
76  static void upsampleNeumann (const Real *from, Real *to, int n, int stride);
77 
78  static inline int modSlow(int x, int n) { int m = x % n; return (m<0) ? m+n : m; }
79  // warning - noiseTileSize has to be 128^3!
80  #define modFast128(x) ((x) & 127)
81 
82  inline Real getTime() { return mParent->getTime() * mParent->getDx() * mTimeAnim; }
83 
84  // pre-compute tile data for wavelet noise
85  void generateTile( int loadFromFile );
86 
87  // animation over time
88  // grid size normalization (inverse size)
89  Real mGsInvX, mGsInvY, mGsInvZ;
90  // random offset into tile to simulate different random seeds
91  Vec3 mSeedOffset;
92 
93  static Real* mNoiseTile;
94  // global random seed storage
95  static int randomSeed;
96 };
97 
98 
99 
100 // **************************************************************************
101 // Implementation
102 
103 #define ADD_WEIGHTED(x,y,z)\
104  weight = 1.0f;\
105  xC = modFast128(midX + (x));\
106  weight *= w[0][(x) + 1];\
107  yC = modFast128(midY + (y));\
108  weight *= w[1][(y) + 1];\
109  zC = modFast128(midZ + (z));\
110  weight *= w[2][(z) + 1];\
111  result += weight * data[(zC * NOISE_TILE_SIZE + yC) * NOISE_TILE_SIZE + xC];
112 
114 // derivatives of 3D noise - unrolled for performance
116 inline Real WaveletNoiseField::WNoiseDx(const Vec3& p, Real *data) {
117  Real w[3][3], t, result = 0;
118 
119  // Evaluate quadratic B-spline basis functions
120  int midX = (int)ceil(p[0] - 0.5f);
121  t = midX - (p[0] - 0.5f);
122  w[0][0] = -t;
123  w[0][2] = (1.f - t);
124  w[0][1] = 2.0f * t - 1.0f;
125 
126  int midY = (int)ceil(p[1] - 0.5f);
127  t = midY - (p[1] - 0.5f);
128  w[1][0] = t * t * 0.5f;
129  w[1][2] = (1.f - t) * (1.f - t) *0.5f;
130  w[1][1] = 1.f - w[1][0] - w[1][2];
131 
132  int midZ = (int)ceil(p[2] - 0.5f);
133  t = midZ - (p[2] - 0.5f);
134  w[2][0] = t * t * 0.5f;
135  w[2][2] = (1.f - t) * (1.f - t) *0.5f;
136  w[2][1] = 1.f - w[2][0] - w[2][2];
137 
138  // Evaluate noise by weighting noise coefficients by basis function values
139  int xC, yC, zC;
140  Real weight = 1;
141 
142  ADD_WEIGHTED(-1,-1, -1); ADD_WEIGHTED( 0,-1, -1); ADD_WEIGHTED( 1,-1, -1);
143  ADD_WEIGHTED(-1, 0, -1); ADD_WEIGHTED( 0, 0, -1); ADD_WEIGHTED( 1, 0, -1);
144  ADD_WEIGHTED(-1, 1, -1); ADD_WEIGHTED( 0, 1, -1); ADD_WEIGHTED( 1, 1, -1);
145 
146  ADD_WEIGHTED(-1,-1, 0); ADD_WEIGHTED( 0,-1, 0); ADD_WEIGHTED( 1,-1, 0);
147  ADD_WEIGHTED(-1, 0, 0); ADD_WEIGHTED( 0, 0, 0); ADD_WEIGHTED( 1, 0, 0);
148  ADD_WEIGHTED(-1, 1, 0); ADD_WEIGHTED( 0, 1, 0); ADD_WEIGHTED( 1, 1, 0);
149 
150  ADD_WEIGHTED(-1,-1, 1); ADD_WEIGHTED( 0,-1, 1); ADD_WEIGHTED( 1,-1, 1);
151  ADD_WEIGHTED(-1, 0, 1); ADD_WEIGHTED( 0, 0, 1); ADD_WEIGHTED( 1, 0, 1);
152  ADD_WEIGHTED(-1, 1, 1); ADD_WEIGHTED( 0, 1, 1); ADD_WEIGHTED( 1, 1, 1);
153 
154  return result;
155 }
156 
157 inline Real WaveletNoiseField::WNoise(const Vec3& p, Real *data) {
158  Real w[3][3], t, result = 0;
159 
160  // Evaluate quadratic B-spline basis functions
161  int midX = (int)ceilf(p[0] - 0.5f);
162  t = midX - (p[0] - 0.5f);
163  w[0][0] = t * t * 0.5f;
164  w[0][2] = (1.f - t) * (1.f - t) *0.5f;
165  w[0][1] = 1.f - w[0][0] - w[0][2];
166 
167  int midY = (int)ceilf(p[1] - 0.5f);
168  t = midY - (p[1] - 0.5f);
169  w[1][0] = t * t * 0.5f;
170  w[1][2] = (1.f - t) * (1.f - t) *0.5f;
171  w[1][1] = 1.f - w[1][0] - w[1][2];
172 
173  int midZ = (int)ceilf(p[2] - 0.5f);
174  t = midZ - (p[2] - 0.5f);
175  w[2][0] = t * t * 0.5f;
176  w[2][2] = (1.f - t) * (1.f - t) *0.5f;
177  w[2][1] = 1.f - w[2][0] - w[2][2];
178 
179  // Evaluate noise by weighting noise coefficients by basis function values
180  int xC, yC, zC;
181  Real weight = 1;
182 
183  ADD_WEIGHTED(-1,-1, -1); ADD_WEIGHTED( 0,-1, -1); ADD_WEIGHTED( 1,-1, -1);
184  ADD_WEIGHTED(-1, 0, -1); ADD_WEIGHTED( 0, 0, -1); ADD_WEIGHTED( 1, 0, -1);
185  ADD_WEIGHTED(-1, 1, -1); ADD_WEIGHTED( 0, 1, -1); ADD_WEIGHTED( 1, 1, -1);
186 
187  ADD_WEIGHTED(-1,-1, 0); ADD_WEIGHTED( 0,-1, 0); ADD_WEIGHTED( 1,-1, 0);
188  ADD_WEIGHTED(-1, 0, 0); ADD_WEIGHTED( 0, 0, 0); ADD_WEIGHTED( 1, 0, 0);
189  ADD_WEIGHTED(-1, 1, 0); ADD_WEIGHTED( 0, 1, 0); ADD_WEIGHTED( 1, 1, 0);
190 
191  ADD_WEIGHTED(-1,-1, 1); ADD_WEIGHTED( 0,-1, 1); ADD_WEIGHTED( 1,-1, 1);
192  ADD_WEIGHTED(-1, 0, 1); ADD_WEIGHTED( 0, 0, 1); ADD_WEIGHTED( 1, 0, 1);
193  ADD_WEIGHTED(-1, 1, 1); ADD_WEIGHTED( 0, 1, 1); ADD_WEIGHTED( 1, 1, 1);
194 
195  return result;
196 }
197 
198 
199 
200 #define ADD_WEIGHTEDX(x,y,z)\
201  weight = dw[0][(x) + 1] * w[1][(y) + 1] * w[2][(z) + 1];\
202  result += weight * neighbors[x + 1][y + 1][z + 1];
203 
204 #define ADD_WEIGHTEDY(x,y,z)\
205  weight = w[0][(x) + 1] * dw[1][(y) + 1] * w[2][(z) + 1];\
206  result += weight * neighbors[x + 1][y + 1][z + 1];
207 
208 #define ADD_WEIGHTEDZ(x,y,z)\
209  weight = w[0][(x) + 1] * w[1][(y) + 1] * dw[2][(z) + 1];\
210  result += weight * neighbors[x + 1][y + 1][z + 1];
211 
213 // compute all derivatives in at once
215 inline Vec3 WaveletNoiseField::WNoiseVec(const Vec3& p, Real *data)
216 {
217  Vec3 final(0.);
218  Real w[3][3];
219  Real dw[3][3];
220  Real result = 0;
221  int xC, yC, zC;
222  Real weight;
223 
224  int midX = (int)ceil(p[0] - 0.5f);
225  int midY = (int)ceil(p[1] - 0.5f);
226  int midZ = (int)ceil(p[2] - 0.5f);
227 
228  Real t0 = midX - (p[0] - 0.5f);
229  Real t1 = midY - (p[1] - 0.5f);
230  Real t2 = midZ - (p[2] - 0.5f);
231 
232  // precache all the neighbors for fast access
233  Real neighbors[3][3][3];
234  for (int z = -1; z <=1; z++)
235  for (int y = -1; y <= 1; y++)
236  for (int x = -1; x <= 1; x++)
237  {
238  xC = modFast128(midX + (x));
239  yC = modFast128(midY + (y));
240  zC = modFast128(midZ + (z));
241  neighbors[x + 1][y + 1][z + 1] = data[zC * NOISE_TILE_SIZE * NOISE_TILE_SIZE + yC * NOISE_TILE_SIZE + xC];
242  }
243 
245  // evaluate splines
247  dw[0][0] = -t0;
248  dw[0][2] = (1.f - t0);
249  dw[0][1] = 2.0f * t0 - 1.0f;
250 
251  dw[1][0] = -t1;
252  dw[1][2] = (1.0f - t1);
253  dw[1][1] = 2.0f * t1 - 1.0f;
254 
255  dw[2][0] = -t2;
256  dw[2][2] = (1.0f - t2);
257  dw[2][1] = 2.0f * t2 - 1.0f;
258 
259  w[0][0] = t0 * t0 * 0.5f;
260  w[0][2] = (1.f - t0) * (1.f - t0) *0.5f;
261  w[0][1] = 1.f - w[0][0] - w[0][2];
262 
263  w[1][0] = t1 * t1 * 0.5f;
264  w[1][2] = (1.f - t1) * (1.f - t1) *0.5f;
265  w[1][1] = 1.f - w[1][0] - w[1][2];
266 
267  w[2][0] = t2 * t2 * 0.5f;
268  w[2][2] = (1.f - t2) * (1.f - t2) *0.5f;
269  w[2][1] = 1.f - w[2][0] - w[2][2];
270 
272  // x derivative
274  result = 0.0f;
275  ADD_WEIGHTEDX(-1,-1, -1); ADD_WEIGHTEDX( 0,-1, -1); ADD_WEIGHTEDX( 1,-1, -1);
276  ADD_WEIGHTEDX(-1, 0, -1); ADD_WEIGHTEDX( 0, 0, -1); ADD_WEIGHTEDX( 1, 0, -1);
277  ADD_WEIGHTEDX(-1, 1, -1); ADD_WEIGHTEDX( 0, 1, -1); ADD_WEIGHTEDX( 1, 1, -1);
278 
279  ADD_WEIGHTEDX(-1,-1, 0); ADD_WEIGHTEDX( 0,-1, 0); ADD_WEIGHTEDX( 1,-1, 0);
280  ADD_WEIGHTEDX(-1, 0, 0); ADD_WEIGHTEDX( 0, 0, 0); ADD_WEIGHTEDX( 1, 0, 0);
281  ADD_WEIGHTEDX(-1, 1, 0); ADD_WEIGHTEDX( 0, 1, 0); ADD_WEIGHTEDX( 1, 1, 0);
282 
283  ADD_WEIGHTEDX(-1,-1, 1); ADD_WEIGHTEDX( 0,-1, 1); ADD_WEIGHTEDX( 1,-1, 1);
284  ADD_WEIGHTEDX(-1, 0, 1); ADD_WEIGHTEDX( 0, 0, 1); ADD_WEIGHTEDX( 1, 0, 1);
285  ADD_WEIGHTEDX(-1, 1, 1); ADD_WEIGHTEDX( 0, 1, 1); ADD_WEIGHTEDX( 1, 1, 1);
286  final[0] = result;
287 
289  // y derivative
291  result = 0.0f;
292  ADD_WEIGHTEDY(-1,-1, -1); ADD_WEIGHTEDY( 0,-1, -1); ADD_WEIGHTEDY( 1,-1, -1);
293  ADD_WEIGHTEDY(-1, 0, -1); ADD_WEIGHTEDY( 0, 0, -1); ADD_WEIGHTEDY( 1, 0, -1);
294  ADD_WEIGHTEDY(-1, 1, -1); ADD_WEIGHTEDY( 0, 1, -1); ADD_WEIGHTEDY( 1, 1, -1);
295 
296  ADD_WEIGHTEDY(-1,-1, 0); ADD_WEIGHTEDY( 0,-1, 0); ADD_WEIGHTEDY( 1,-1, 0);
297  ADD_WEIGHTEDY(-1, 0, 0); ADD_WEIGHTEDY( 0, 0, 0); ADD_WEIGHTEDY( 1, 0, 0);
298  ADD_WEIGHTEDY(-1, 1, 0); ADD_WEIGHTEDY( 0, 1, 0); ADD_WEIGHTEDY( 1, 1, 0);
299 
300  ADD_WEIGHTEDY(-1,-1, 1); ADD_WEIGHTEDY( 0,-1, 1); ADD_WEIGHTEDY( 1,-1, 1);
301  ADD_WEIGHTEDY(-1, 0, 1); ADD_WEIGHTEDY( 0, 0, 1); ADD_WEIGHTEDY( 1, 0, 1);
302  ADD_WEIGHTEDY(-1, 1, 1); ADD_WEIGHTEDY( 0, 1, 1); ADD_WEIGHTEDY( 1, 1, 1);
303  final[1] = result;
304 
306  // z derivative
308  result = 0.0f;
309  ADD_WEIGHTEDZ(-1,-1, -1); ADD_WEIGHTEDZ( 0,-1, -1); ADD_WEIGHTEDZ( 1,-1, -1);
310  ADD_WEIGHTEDZ(-1, 0, -1); ADD_WEIGHTEDZ( 0, 0, -1); ADD_WEIGHTEDZ( 1, 0, -1);
311  ADD_WEIGHTEDZ(-1, 1, -1); ADD_WEIGHTEDZ( 0, 1, -1); ADD_WEIGHTEDZ( 1, 1, -1);
312 
313  ADD_WEIGHTEDZ(-1,-1, 0); ADD_WEIGHTEDZ( 0,-1, 0); ADD_WEIGHTEDZ( 1,-1, 0);
314  ADD_WEIGHTEDZ(-1, 0, 0); ADD_WEIGHTEDZ( 0, 0, 0); ADD_WEIGHTEDZ( 1, 0, 0);
315  ADD_WEIGHTEDZ(-1, 1, 0); ADD_WEIGHTEDZ( 0, 1, 0); ADD_WEIGHTEDZ( 1, 1, 0);
316 
317  ADD_WEIGHTEDZ(-1,-1, 1); ADD_WEIGHTEDZ( 0,-1, 1); ADD_WEIGHTEDZ( 1,-1, 1);
318  ADD_WEIGHTEDZ(-1, 0, 1); ADD_WEIGHTEDZ( 0, 0, 1); ADD_WEIGHTEDZ( 1, 0, 1);
319  ADD_WEIGHTEDZ(-1, 1, 1); ADD_WEIGHTEDZ( 0, 1, 1); ADD_WEIGHTEDZ( 1, 1, 1);
320  final[2] = result;
321 
322  //debMsg("FINAL","at "<<p<<" = "<<final); // DEBUG
323  return final;
324 }
325 #undef ADD_WEIGHTEDX
326 #undef ADD_WEIGHTEDY
327 #undef ADD_WEIGHTEDZ
328 
329 inline Real WaveletNoiseField::evaluate(Vec3 pos, int tile) {
330  pos[0] *= mGsInvX;
331  pos[1] *= mGsInvY;
332  pos[2] *= mGsInvZ;
333  pos += mSeedOffset;
334 
335  // time anim
336  pos += Vec3(getTime());
337 
338  pos[0] *= mPosScale[0];
339  pos[1] *= mPosScale[1];
340  pos[2] *= mPosScale[2];
341  pos += mPosOffset;
342 
343  const int n3 = square(NOISE_TILE_SIZE) * NOISE_TILE_SIZE;
344  Real v = WNoise(pos, &mNoiseTile[tile*n3]);
345 
346  v += mValOffset;
347  v *= mValScale;
348  if (mClamp) {
349  if (v< mClampNeg) v = mClampNeg;
350  if (v> mClampPos) v = mClampPos;
351  }
352  return v;
353 }
354 
355 inline Vec3 WaveletNoiseField::evaluateVec(Vec3 pos, int tile) {
356  pos[0] *= mGsInvX;
357  pos[1] *= mGsInvY;
358  pos[2] *= mGsInvZ;
359  pos += mSeedOffset;
360 
361  // time anim
362  pos += Vec3(getTime());
363 
364  pos[0] *= mPosScale[0];
365  pos[1] *= mPosScale[1];
366  pos[2] *= mPosScale[2];
367  pos += mPosOffset;
368 
369  const int n3 = square(NOISE_TILE_SIZE) * NOISE_TILE_SIZE;
370  Vec3 v = WNoiseVec(pos, &mNoiseTile[tile*n3]);
371 
372  v += Vec3(mValOffset);
373  v *= mValScale;
374 
375  if (mClamp) {
376  for(int i=0; i<3; i++) {
377  if (v[i]< mClampNeg) v[i] = mClampNeg;
378  if (v[i]> mClampPos) v[i] = mClampPos;
379  }
380  }
381  return v;
382 }
383 
384 inline Vec3 WaveletNoiseField::evaluateCurl(Vec3 pos) {
385  // gradients of w0-w2
386  Vec3 d0 = evaluateVec(pos,0),
387  d1 = evaluateVec(pos,1),
388  d2 = evaluateVec(pos,2);
389 
390  return Vec3(d0.y-d1.z, d2.z-d0.x, d1.x-d2.y);
391 }
392 
393 } // namespace
394 
395 #endif
396 
397 
Manta
Definition: commonkernels.h:22
Manta::WaveletNoiseField::evaluateVec
Vec3 evaluateVec(Vec3 pos, int tile=0)
evaluate noise as a vector
Definition: noisefield.h:355
Manta::WaveletNoiseField::computeCoefficients
static void computeCoefficients(Grid< Real > &input, Grid< Real > &tempIn1, Grid< Real > &tempIn2)
compute wavelet decomposition of an input grid (stores residual coefficients)
Definition: noisefield.cpp:231
Manta::WaveletNoiseField
Definition: noisefield.h:28
Manta::WaveletNoiseField::evaluateCurl
Vec3 evaluateCurl(Vec3 pos)
evaluate curl noise
Definition: noisefield.h:384
Manta::Vector3D
Basic inlined vector class.
Definition: vectorbase.h:71
Manta::WaveletNoiseField::evaluate
Real evaluate(Vec3 pos, int tile=0)
evaluate noise
Definition: noisefield.h:329
Manta::WaveletNoiseField::data
Real * data()
direct data access
Definition: noisefield.h:43
Manta::FluidSolver
Definition: fluidsolver.h:28
Generated on Tue Sep 19 2017 17:49:04 for mantaflow by   doxygen 1.8.13