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