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