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