BF2.5: First commit of smoke code.
[blender.git] / intern / smoke / intern / FLUID_3D_STATIC.cpp
1 //////////////////////////////////////////////////////////////////////
2 // This file is part of Wavelet Turbulence.
3 //
4 // Wavelet Turbulence is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // Wavelet Turbulence is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
16 //
17 // Copyright 2008 Theodore Kim and Nils Thuerey
18 //
19 // FLUID_3D.cpp: implementation of the static functions of the FLUID_3D class.
20 //
21 //////////////////////////////////////////////////////////////////////
22
23 #include <zlib.h>
24 #include "FLUID_3D.h"
25 #include "IMAGE.h"
26 #include "WTURBULENCE.h"
27 #include "INTERPOLATE.h"
28
29 //////////////////////////////////////////////////////////////////////
30 // add a test cube of density to the center
31 //////////////////////////////////////////////////////////////////////
32 void FLUID_3D::addSmokeColumn() {
33         addSmokeTestCase(_density, _res, 1.0);
34         // addSmokeTestCase(_zVelocity, _res, 1.0);
35         addSmokeTestCase(_heat, _res, 1.0);
36         if (_wTurbulence) {
37                 addSmokeTestCase(_wTurbulence->getDensityBig(), _wTurbulence->getResBig(), 1.0);
38         }
39 }
40
41 //////////////////////////////////////////////////////////////////////
42 // generic static version, so that it can be applied to the
43 // WTURBULENCE grid as well
44 //////////////////////////////////////////////////////////////////////
45 void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res, float value)
46 {
47         const int slabSize = res[0]*res[1]; int maxRes = (int)MAX3V(res);
48         float dx = 1.0f / (float)maxRes;
49
50         float xTotal = dx * res[0];
51         float yTotal = dx * res[1];
52         float zTotal = dx * res[2];
53
54   float heighMin = 0.05;
55   float heighMax = 0.10;
56
57   for (int y = 0; y < res[1]; y++)
58     for (int z = (int)(heighMin*res[2]); z <= (int)(heighMax * res[1]); z++)
59       for (int x = 0; x < res[0]; x++)
60       {
61         float xLength = x * dx - xTotal * 0.4f;
62         float yLength = y * dx - zTotal * 0.5f;
63         float radius = sqrtf(xLength * xLength + yLength * yLength);
64
65         if (radius < 0.075f * xTotal)
66         {
67           int index = x + y * res[0] + z * slabSize;
68           field[index] = value;
69         }
70       }
71 }
72
73 //////////////////////////////////////////////////////////////////////
74 // set x direction to Neumann boundary conditions
75 //////////////////////////////////////////////////////////////////////
76 void FLUID_3D::setNeumannX(float* field, Vec3Int res)
77 {
78         const int slabSize = res[0] * res[1];
79         int index;
80         for (int z = 0; z < res[2]; z++)
81                 for (int y = 0; y < res[1]; y++)
82                 {
83                         // left slab
84                         index = y * res[0] + z * slabSize;
85                         field[index] = field[index + 2];
86
87                         // right slab
88                         index += res[0] - 1;
89                         field[index] = field[index - 2];
90                 }
91  }
92
93 //////////////////////////////////////////////////////////////////////
94 // set y direction to Neumann boundary conditions
95 //////////////////////////////////////////////////////////////////////
96 void FLUID_3D::setNeumannY(float* field, Vec3Int res)
97 {
98         const int slabSize = res[0] * res[1];
99         int index;
100         for (int z = 0; z < res[2]; z++)
101                 for (int x = 0; x < res[0]; x++)
102                 {
103                         // bottom slab
104                         index = x + z * slabSize;
105                         field[index] = field[index + 2 * res[0]];
106
107                         // top slab
108                         index += slabSize - res[0];
109                         field[index] = field[index - 2 * res[0]];
110                 }
111
112         // fix, force top slab to only allow outwards flux
113         for (int z = 0; z < res[2]; z++)
114                 for (int x = 0; x < res[0]; x++)
115                 {
116                         // top slab
117                         int index = x + z * slabSize;
118                         index += slabSize - res[0];
119                         if(field[index]<0.) field[index] = 0.;
120                         index -= res[0];
121                         if(field[index]<0.) field[index] = 0.;
122                 }
123 }
124
125 //////////////////////////////////////////////////////////////////////
126 // set z direction to Neumann boundary conditions
127 //////////////////////////////////////////////////////////////////////
128 void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
129 {
130         const int slabSize = res[0] * res[1];
131         const int totalCells = res[0] * res[1] * res[2];
132         int index;
133         for (int y = 0; y < res[1]; y++)
134                 for (int x = 0; x < res[0]; x++)
135                 {
136                         // front slab
137                         index = x + y * res[0];
138                         field[index] = field[index + 2 * slabSize];
139
140                         // back slab
141                         index += totalCells - slabSize;
142                         field[index] = field[index - 2 * slabSize];
143                 }
144 }
145
146 //////////////////////////////////////////////////////////////////////
147 // set x direction to zero
148 //////////////////////////////////////////////////////////////////////
149 void FLUID_3D::setZeroX(float* field, Vec3Int res)
150 {
151         const int slabSize = res[0] * res[1];
152         int index;
153         for (int z = 0; z < res[2]; z++)
154                 for (int y = 0; y < res[1]; y++)
155                 {
156                         // left slab
157                         index = y * res[0] + z * slabSize;
158                         field[index] = 0.0f;
159
160                         // right slab
161                         index += res[0] - 1;
162                         field[index] = 0.0f;
163                 }
164 }
165
166 //////////////////////////////////////////////////////////////////////
167 // set y direction to zero
168 //////////////////////////////////////////////////////////////////////
169 void FLUID_3D::setZeroY(float* field, Vec3Int res)
170 {
171         const int slabSize = res[0] * res[1];
172         int index;
173         for (int z = 0; z < res[2]; z++)
174                 for (int x = 0; x < res[0]; x++)
175                 {
176                         // bottom slab
177                         index = x + z * slabSize;
178                         field[index] = 0.0f;
179
180                         // top slab
181                         index += slabSize - res[0];
182                         field[index] = 0.0f;
183                 }
184 }
185
186 //////////////////////////////////////////////////////////////////////
187 // set z direction to zero
188 //////////////////////////////////////////////////////////////////////
189 void FLUID_3D::setZeroZ(float* field, Vec3Int res)
190 {
191         const int slabSize = res[0] * res[1];
192         const int totalCells = res[0] * res[1] * res[2];
193         int index;
194         for (int y = 0; y < res[1]; y++)
195                 for (int x = 0; x < res[0]; x++)
196                 {
197                         // front slab
198                         index = x + y * res[0];
199                         field[index] = 0.0f;
200
201                         // back slab
202                         index += totalCells - slabSize;
203                         field[index] = 0.0f;
204                 }
205  }
206
207 //////////////////////////////////////////////////////////////////////
208 // copy grid boundary
209 //////////////////////////////////////////////////////////////////////
210 void FLUID_3D::copyBorderX(float* field, Vec3Int res)
211 {
212         const int slabSize = res[0] * res[1];
213         int index;
214         for (int z = 0; z < res[2]; z++)
215                 for (int y = 0; y < res[1]; y++)
216                 {
217                         // left slab
218                         index = y * res[0] + z * slabSize;
219                         field[index] = field[index + 1];
220
221                         // right slab
222                         index += res[0] - 1;
223                         field[index] = field[index - 1];
224                 }
225 }
226 void FLUID_3D::copyBorderY(float* field, Vec3Int res)
227 {
228         const int slabSize = res[0] * res[1];
229         const int totalCells = res[0] * res[1] * res[2];
230         int index;
231         for (int z = 0; z < res[2]; z++)
232                 for (int x = 0; x < res[0]; x++)
233                 {
234                         // bottom slab
235                         index = x + z * slabSize;
236                         field[index] = field[index + res[0]];
237                         // top slab
238                         index += slabSize - res[0];
239                         field[index] = field[index - res[0]];
240                 }
241 }
242 void FLUID_3D::copyBorderZ(float* field, Vec3Int res)
243 {
244         const int slabSize = res[0] * res[1];
245         const int totalCells = res[0] * res[1] * res[2];
246         int index;
247         for (int y = 0; y < res[1]; y++)
248                 for (int x = 0; x < res[0]; x++)
249                 {
250                         // front slab
251                         index = x + y * res[0];
252                         field[index] = field[index + slabSize];
253                         // back slab
254                         index += totalCells - slabSize;
255                         field[index] = field[index - slabSize];
256                 }
257 }
258
259 /////////////////////////////////////////////////////////////////////
260 // advect field with the semi lagrangian method
261 //////////////////////////////////////////////////////////////////////
262 void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
263                 float* oldField, float* newField, Vec3Int res)
264 {
265         const int xres = res[0];
266         const int yres = res[1];
267         const int zres = res[2];
268         static int hits = 0;
269         static int total = 0;
270         const int slabSize = res[0] * res[1];
271
272         // scale dt up to grid resolution
273 #if PARALLEL==1
274 #pragma omp parallel for schedule(static)
275 #endif
276         for (int z = 0; z < zres; z++)
277                 for (int y = 0; y < yres; y++)
278                         for (int x = 0; x < xres; x++)
279                         {
280                                 const int index = x + y * xres + z * xres*yres;
281
282         // backtrace
283                                 float xTrace = x - dt * velx[index];
284                                 float yTrace = y - dt * vely[index];
285                                 float zTrace = z - dt * velz[index];
286
287                                 // clamp backtrace to grid boundaries
288                                 if (xTrace < 0.5) xTrace = 0.5;
289                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
290                                 if (yTrace < 0.5) yTrace = 0.5;
291                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
292                                 if (zTrace < 0.5) zTrace = 0.5;
293                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
294
295                                 // locate neighbors to interpolate
296                                 const int x0 = (int)xTrace;
297                                 const int x1 = x0 + 1;
298                                 const int y0 = (int)yTrace;
299                                 const int y1 = y0 + 1;
300                                 const int z0 = (int)zTrace;
301                                 const int z1 = z0 + 1;
302
303                                 // get interpolation weights
304                                 const float s1 = xTrace - x0;
305                                 const float s0 = 1.0f - s1;
306                                 const float t1 = yTrace - y0;
307                                 const float t0 = 1.0f - t1;
308                                 const float u1 = zTrace - z0;
309                                 const float u0 = 1.0f - u1;
310
311                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
312                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
313                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
314                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
315                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
316                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
317                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
318                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
319
320                                 // interpolate
321                                 // (indices could be computed once)
322                                 newField[index] = u0 * (s0 * (t0 * oldField[i000] +
323                                                         t1 * oldField[i010]) +
324                                                 s1 * (t0 * oldField[i100] +
325                                                         t1 * oldField[i110])) +
326                                         u1 * (s0 * (t0 * oldField[i001] +
327                                                                 t1 * oldField[i011]) +
328                                                         s1 * (t0 * oldField[i101] +
329                                                                 t1 * oldField[i111]));
330                         }
331 }
332
333 /////////////////////////////////////////////////////////////////////
334 // advect field with the maccormack method
335 //
336 // comments are the pseudocode from selle's paper
337 //////////////////////////////////////////////////////////////////////
338 void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity,
339                 float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const float* obstacles)
340 {
341         float* phiHatN  = temp1;
342         float* phiHatN1 = temp2;
343         const int sx= res[0];
344         const int sy= res[1];
345         const int sz= res[2];
346
347         for (int x = 0; x < sx * sy * sz; x++)
348                 phiHatN[x] = phiHatN1[x] = oldField[x];
349
350         float*& phiN    = oldField;
351         float*& phiN1   = newField;
352
353         // phiHatN1 = A(phiN)
354         advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiHatN1, res);
355
356         // phiHatN = A^R(phiHatN1)
357         advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN1, phiHatN, res);
358
359         // phiN1 = phiHatN1 + (phiN - phiHatN) / 2
360         const int border = 0;
361         for (int z = border; z < sz-border; z++)
362                 for (int y = border; y < sy-border; y++)
363                         for (int x = border; x < sx-border; x++) {
364                                 int index = x + y * sx + z * sx*sy;
365                                 phiN1[index] = phiHatN1[index] + (phiN[index] - phiHatN[index]) * 0.50f;
366                                 //phiN1[index] = phiHatN1[index]; // debug, correction off
367                         }
368         copyBorderX(phiN1, res);
369         copyBorderY(phiN1, res);
370         copyBorderZ(phiN1, res);
371
372         // clamp any newly created extrema
373         clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res);
374
375         // if the error estimate was bad, revert to first order
376         clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN1);
377 }
378
379
380 //////////////////////////////////////////////////////////////////////
381 // Clamp the extrema generated by the BFECC error correction
382 //////////////////////////////////////////////////////////////////////
383 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
384                 float* oldField, float* newField, Vec3Int res)
385 {
386         const int xres= res[0];
387         const int yres= res[1];
388         const int zres= res[2];
389         const int slabSize = res[0] * res[1];
390         for (int z = 1; z < zres-1; z++)
391                 for (int y = 1; y < yres-1; y++)
392                         for (int x = 1; x < xres-1; x++)
393                         {
394                                 const int index = x + y * xres+ z * xres*yres;
395                                 // backtrace
396                                 float xTrace = x - dt * velx[index];
397                                 float yTrace = y - dt * vely[index];
398                                 float zTrace = z - dt * velz[index];
399
400                                 // clamp backtrace to grid boundaries
401                                 if (xTrace < 0.5) xTrace = 0.5;
402                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
403                                 if (yTrace < 0.5) yTrace = 0.5;
404                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
405                                 if (zTrace < 0.5) zTrace = 0.5;
406                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
407
408                                 // locate neighbors to interpolate
409                                 const int x0 = (int)xTrace;
410                                 const int x1 = x0 + 1;
411                                 const int y0 = (int)yTrace;
412                                 const int y1 = y0 + 1;
413                                 const int z0 = (int)zTrace;
414                                 const int z1 = z0 + 1;
415
416                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
417                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
418                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
419                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
420                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
421                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
422                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
423                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
424
425                                 float minField = oldField[i000];
426                                 float maxField = oldField[i000];
427
428                                 minField = (oldField[i010] < minField) ? oldField[i010] : minField;
429                                 maxField = (oldField[i010] > maxField) ? oldField[i010] : maxField;
430
431                                 minField = (oldField[i100] < minField) ? oldField[i100] : minField;
432                                 maxField = (oldField[i100] > maxField) ? oldField[i100] : maxField;
433
434                                 minField = (oldField[i110] < minField) ? oldField[i110] : minField;
435                                 maxField = (oldField[i110] > maxField) ? oldField[i110] : maxField;
436
437                                 minField = (oldField[i001] < minField) ? oldField[i001] : minField;
438                                 maxField = (oldField[i001] > maxField) ? oldField[i001] : maxField;
439
440                                 minField = (oldField[i011] < minField) ? oldField[i011] : minField;
441                                 maxField = (oldField[i011] > maxField) ? oldField[i011] : maxField;
442
443                                 minField = (oldField[i101] < minField) ? oldField[i101] : minField;
444                                 maxField = (oldField[i101] > maxField) ? oldField[i101] : maxField;
445
446                                 minField = (oldField[i111] < minField) ? oldField[i111] : minField;
447                                 maxField = (oldField[i111] > maxField) ? oldField[i111] : maxField;
448
449                                 newField[index] = (newField[index] > maxField) ? maxField : newField[index];
450                                 newField[index] = (newField[index] < minField) ? minField : newField[index];
451                         }
452 }
453
454 //////////////////////////////////////////////////////////////////////
455 // Reverts any backtraces that go into boundaries back to first
456 // order -- in this case the error correction term was totally
457 // incorrect
458 //////////////////////////////////////////////////////////////////////
459 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely, const float* velz,
460                 float* oldField, float* newField, Vec3Int res, const float* obstacles, const float *oldAdvection)
461 {
462         const int sx= res[0];
463         const int sy= res[1];
464         const int sz= res[2];
465         const int slabSize = res[0] * res[1];
466         for (int z = 1; z < sz-1; z++)
467                 for (int y = 1; y < sy-1; y++)
468                         for (int x = 1; x < sx-1; x++)
469                         {
470                                 const int index = x + y * sx+ z * slabSize;
471                                 // backtrace
472                                 float xBackward = x + dt * velx[index];
473                                 float yBackward = y + dt * vely[index];
474                                 float zBackward = z + dt * velz[index];
475                                 float xTrace    = x - dt * velx[index];
476                                 float yTrace    = y - dt * vely[index];
477                                 float zTrace    = z - dt * velz[index];
478
479                                 // see if it goes outside the boundaries
480                                 bool hasObstacle =
481                                         (zTrace < 1.0f)    || (zTrace > sz - 2.0f) ||
482                                         (yTrace < 1.0f)    || (yTrace > sy - 2.0f) ||
483                                         (xTrace < 1.0f)    || (xTrace > sx - 2.0f) ||
484                                         (zBackward < 1.0f) || (zBackward > sz - 2.0f) ||
485                                         (yBackward < 1.0f) || (yBackward > sy - 2.0f) ||
486                                         (xBackward < 1.0f) || (xBackward > sx - 2.0f);
487                                 // reuse old advection instead of doing another one...
488                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
489
490                                 // clamp to prevent an out of bounds access when looking into
491                                 // the _obstacles array
492                                 zTrace = (zTrace < 0.5f) ? 0.5f : zTrace;
493                                 zTrace = (zTrace > sz - 1.5f) ? sz - 1.5f : zTrace;
494                                 yTrace = (yTrace < 0.5f) ? 0.5f : yTrace;
495                                 yTrace = (yTrace > sy - 1.5f) ? sy - 1.5f : yTrace;
496                                 xTrace = (xTrace < 0.5f) ? 0.5f : xTrace;
497                                 xTrace = (xTrace > sx - 1.5f) ? sx - 1.5f : xTrace;
498
499                                 // locate neighbors to interpolate,
500                                 // do backward first since we will use the forward indices if a
501                                 // reversion is actually necessary
502                                 zBackward = (zBackward < 0.5f) ? 0.5f : zBackward;
503                                 zBackward = (zBackward > sz - 1.5f) ? sz - 1.5f : zBackward;
504                                 yBackward = (yBackward < 0.5f) ? 0.5f : yBackward;
505                                 yBackward = (yBackward > sy - 1.5f) ? sy - 1.5f : yBackward;
506                                 xBackward = (xBackward < 0.5f) ? 0.5f : xBackward;
507                                 xBackward = (xBackward > sx - 1.5f) ? sx - 1.5f : xBackward;
508
509                                 int x0 = (int)xBackward;
510                                 int x1 = x0 + 1;
511                                 int y0 = (int)yBackward;
512                                 int y1 = y0 + 1;
513                                 int z0 = (int)zBackward;
514                                 int z1 = z0 + 1;
515                                 if(obstacles && !hasObstacle) {
516                                         hasObstacle = hasObstacle ||
517                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
518                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
519                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
520                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
521                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
522                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
523                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
524                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
525                                 }
526                                 // reuse old advection instead of doing another one...
527                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
528
529                                 x0 = (int)xTrace;
530                                 x1 = x0 + 1;
531                                 y0 = (int)yTrace;
532                                 y1 = y0 + 1;
533                                 z0 = (int)zTrace;
534                                 z1 = z0 + 1;
535                                 if(obstacles && !hasObstacle) {
536                                         hasObstacle = hasObstacle ||
537                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
538                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
539                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
540                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
541                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
542                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
543                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
544                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
545                                 } // obstacle array
546                                 // reuse old advection instead of doing another one...
547                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
548
549                                 // see if either the forward or backward ray went into
550                                 // a boundary
551                                 if (hasObstacle) {
552                                         // get interpolation weights
553                                         float s1 = xTrace - x0;
554                                         float s0 = 1.0f - s1;
555                                         float t1 = yTrace - y0;
556                                         float t0 = 1.0f - t1;
557                                         float u1 = zTrace - z0;
558                                         float u0 = 1.0f - u1;
559
560                                         const int i000 = x0 + y0 * sx + z0 * slabSize;
561                                         const int i010 = x0 + y1 * sx + z0 * slabSize;
562                                         const int i100 = x1 + y0 * sx + z0 * slabSize;
563                                         const int i110 = x1 + y1 * sx + z0 * slabSize;
564                                         const int i001 = x0 + y0 * sx + z1 * slabSize;
565                                         const int i011 = x0 + y1 * sx + z1 * slabSize;
566                                         const int i101 = x1 + y0 * sx + z1 * slabSize;
567                                         const int i111 = x1 + y1 * sx + z1 * slabSize;
568
569                                         // interpolate, (indices could be computed once)
570                                         newField[index] = u0 * (s0 * (
571                                                                 t0 * oldField[i000] +
572                                                                 t1 * oldField[i010]) +
573                                                         s1 * (t0 * oldField[i100] +
574                                                                 t1 * oldField[i110])) +
575                                                 u1 * (s0 * (t0 * oldField[i001] +
576                                                                         t1 * oldField[i011]) +
577                                                                 s1 * (t0 * oldField[i101] +
578                                                                         t1 * oldField[i111]));
579                                 }
580                         } // xyz
581 }
582
583 //////////////////////////////////////////////////////////////////////
584 // image output
585 //////////////////////////////////////////////////////////////////////
586 /*
587 void FLUID_3D::writeImageSliceXY(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) {
588   writeProjectedIntern(field, res, 0,1, prefix, picCnt, scale);
589 }
590 void FLUID_3D::writeImageSliceYZ(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) {
591   writeProjectedIntern(field, res, 1,2, prefix, picCnt, scale);
592 }
593 void FLUID_3D::writeImageSliceXZ(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) {
594   writeProjectedIntern(field, res, 0,2, prefix, picCnt, scale);
595 }
596 */
597
598 //////////////////////////////////////////////////////////////////////
599 // Helper function for projecting densities along a dimension
600 //////////////////////////////////////////////////////////////////////
601 static int getOtherDir(int dir1, int dir2) {
602         switch(dir1) {
603                 case 0:
604                         switch(dir2) {
605                                 case 1: return 2;
606                                 case 2: return 1; }
607                         break;
608                 case 1:
609                         switch(dir2) {
610                                 case 0: return 2;
611                                 case 2: return 0; }
612                         break;
613                 case 2:
614                         switch(dir2) {
615                                 case 0: return 1;
616                                 case 1: return 0; }
617                         break;
618         }
619 }
620
621 //////////////////////////////////////////////////////////////////////
622 // average densities along third spatial direction
623 //////////////////////////////////////////////////////////////////////
624 void FLUID_3D::writeProjectedIntern(const float *field, Vec3Int res,
625                 int dir1, int dir2, string prefix, int picCnt, float scale) {
626         const int nitems = res[dir1]*res[dir2];
627         const int otherDir = getOtherDir(dir1,dir2);
628         float *buf = new float[nitems];
629         Vec3Int min = Vec3Int(0);
630         Vec3Int max = res;
631
632         min[otherDir] = 0;
633         max[otherDir] = res[otherDir];
634         float div = 1./(float)MIN3V(res); // normalize for shorter sides, old: res[otherDir];
635         div *= 4.; //slightly increase contrast
636         for(int i=0; i<nitems; i++) buf[i] = 0.;
637
638         Vec3Int cnt = 0;
639         for (cnt[2] = min[2]; cnt[2] < max[2]; cnt[2]++) {
640                 for (cnt[1] = min[1]; cnt[1] < max[1]; cnt[1]++)
641                         for (cnt[0] = min[0]; cnt[0] < max[0]; cnt[0]++)
642                         {
643                                 const int index = cnt[0] + cnt[1] * res[0] + cnt[2] * res[0]*res[1];
644                                 const int bufindex = cnt[dir1] + cnt[dir2] * res[dir1];
645                                 buf[bufindex] += field[index] * scale *div;
646                         }
647         }
648         IMAGE::dumpNumberedPNG(picCnt, prefix, buf, res[dir1], res[dir2]);
649         delete[] buf;
650 }