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