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, float value)
49 {
50         const int slabSize = res[0]*res[1]; int maxRes = (int)MAX3V(res);
51         float dx = 1.0f / (float)maxRes;
52
53         float xTotal = dx * res[0];
54         float yTotal = dx * res[1];
55         float zTotal = dx * res[2];
56
57   float heighMin = 0.05;
58   float heighMax = 0.10;
59
60   for (int y = 0; y < res[1]; y++)
61     for (int z = (int)(heighMin*res[2]); z <= (int)(heighMax * res[1]); z++)
62       for (int x = 0; x < res[0]; x++)
63       {
64         float xLength = x * dx - xTotal * 0.4f;
65         float yLength = y * dx - zTotal * 0.5f;
66         float radius = sqrtf(xLength * xLength + yLength * yLength);
67
68         if (radius < 0.075f * xTotal)
69         {
70           int index = x + y * res[0] + z * slabSize;
71           field[index] = value;
72         }
73       }
74 }
75 */
76
77 //////////////////////////////////////////////////////////////////////
78 // set x direction to Neumann boundary conditions
79 //////////////////////////////////////////////////////////////////////
80 void FLUID_3D::setNeumannX(float* field, Vec3Int res)
81 {
82         const int slabSize = res[0] * res[1];
83         int index;
84         for (int z = 0; z < res[2]; z++)
85                 for (int y = 0; y < res[1]; y++)
86                 {
87                         // left slab
88                         index = y * res[0] + z * slabSize;
89                         field[index] = field[index + 2];
90
91                         // right slab
92                         index += res[0] - 1;
93                         field[index] = field[index - 2];
94                 }
95
96         // 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         static int hits = 0;
299         static int total = 0;
300         const int slabSize = res[0] * res[1];
301
302         // scale dt up to grid resolution
303 #if PARALLEL==1 && !_WIN32
304 #pragma omp parallel
305 #pragma omp for  schedule(static)
306 #endif
307         for (int z = 0; z < zres; z++)
308                 for (int y = 0; y < yres; y++)
309                         for (int x = 0; x < xres; x++)
310                         {
311                                 const int index = x + y * xres + z * xres*yres;
312                                 
313         // backtrace
314                                 float xTrace = x - dt * velx[index];
315                                 float yTrace = y - dt * vely[index];
316                                 float zTrace = z - dt * velz[index];
317
318                                 // clamp backtrace to grid boundaries
319                                 if (xTrace < 0.5) xTrace = 0.5;
320                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
321                                 if (yTrace < 0.5) yTrace = 0.5;
322                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
323                                 if (zTrace < 0.5) zTrace = 0.5;
324                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
325
326                                 // locate neighbors to interpolate
327                                 const int x0 = (int)xTrace;
328                                 const int x1 = x0 + 1;
329                                 const int y0 = (int)yTrace;
330                                 const int y1 = y0 + 1;
331                                 const int z0 = (int)zTrace;
332                                 const int z1 = z0 + 1;
333
334                                 // get interpolation weights
335                                 const float s1 = xTrace - x0;
336                                 const float s0 = 1.0f - s1;
337                                 const float t1 = yTrace - y0;
338                                 const float t0 = 1.0f - t1;
339                                 const float u1 = zTrace - z0;
340                                 const float u0 = 1.0f - u1;
341
342                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
343                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
344                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
345                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
346                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
347                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
348                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
349                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
350
351                                 // interpolate
352                                 // (indices could be computed once)
353                                 newField[index] = u0 * (s0 * (t0 * oldField[i000] +
354                                                         t1 * oldField[i010]) +
355                                                 s1 * (t0 * oldField[i100] +
356                                                         t1 * oldField[i110])) +
357                                         u1 * (s0 * (t0 * oldField[i001] +
358                                                                 t1 * oldField[i011]) +
359                                                         s1 * (t0 * oldField[i101] +
360                                                                 t1 * oldField[i111]));
361                         }
362 }
363
364 /////////////////////////////////////////////////////////////////////
365 // advect field with the maccormack method
366 //
367 // comments are the pseudocode from selle's paper
368 //////////////////////////////////////////////////////////////////////
369 void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
370                                 float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles)
371 {
372         float* phiHatN  = temp1;
373         float* phiHatN1 = temp2;
374         const int sx= res[0];
375         const int sy= res[1];
376         const int sz= res[2];
377
378         for (int x = 0; x < sx * sy * sz; x++)
379                 phiHatN[x] = phiHatN1[x] = oldField[x];
380
381         float*& phiN    = oldField;
382         float*& phiN1   = newField;
383
384         // phiHatN1 = A(phiN)
385         advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiHatN1, res);
386
387         // phiHatN = A^R(phiHatN1)
388         advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN1, phiHatN, res);
389
390         // phiN1 = phiHatN1 + (phiN - phiHatN) / 2
391         const int border = 0; 
392         for (int z = border; z < sz-border; z++)
393                 for (int y = border; y < sy-border; y++)
394                         for (int x = border; x < sx-border; x++) {
395                                 int index = x + y * sx + z * sx*sy;
396                                 phiN1[index] = phiHatN1[index] + (phiN[index] - phiHatN[index]) * 0.50f;
397                                 //phiN1[index] = phiHatN1[index]; // debug, correction off
398                         }
399         copyBorderX(phiN1, res);
400         copyBorderY(phiN1, res);
401         copyBorderZ(phiN1, res);
402
403         // clamp any newly created extrema
404         clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res);
405
406         // if the error estimate was bad, revert to first order
407         clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN1);
408
409
410
411 //////////////////////////////////////////////////////////////////////
412 // Clamp the extrema generated by the BFECC error correction
413 //////////////////////////////////////////////////////////////////////
414 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
415                 float* oldField, float* newField, Vec3Int res)
416 {
417         const int xres= res[0];
418         const int yres= res[1];
419         const int zres= res[2];
420         const int slabSize = res[0] * res[1];
421         for (int z = 1; z < zres-1; z++)
422                 for (int y = 1; y < yres-1; y++)
423                         for (int x = 1; x < xres-1; x++)
424                         {
425                                 const int index = x + y * xres+ z * xres*yres;
426                                 // backtrace
427                                 float xTrace = x - dt * velx[index];
428                                 float yTrace = y - dt * vely[index];
429                                 float zTrace = z - dt * velz[index];
430
431                                 // clamp backtrace to grid boundaries
432                                 if (xTrace < 0.5) xTrace = 0.5;
433                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
434                                 if (yTrace < 0.5) yTrace = 0.5;
435                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
436                                 if (zTrace < 0.5) zTrace = 0.5;
437                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
438
439                                 // locate neighbors to interpolate
440                                 const int x0 = (int)xTrace;
441                                 const int x1 = x0 + 1;
442                                 const int y0 = (int)yTrace;
443                                 const int y1 = y0 + 1;
444                                 const int z0 = (int)zTrace;
445                                 const int z1 = z0 + 1;
446
447                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
448                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
449                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
450                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
451                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
452                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
453                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
454                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
455
456                                 float minField = oldField[i000];
457                                 float maxField = oldField[i000];
458
459                                 minField = (oldField[i010] < minField) ? oldField[i010] : minField;
460                                 maxField = (oldField[i010] > maxField) ? oldField[i010] : maxField;
461
462                                 minField = (oldField[i100] < minField) ? oldField[i100] : minField;
463                                 maxField = (oldField[i100] > maxField) ? oldField[i100] : maxField;
464
465                                 minField = (oldField[i110] < minField) ? oldField[i110] : minField;
466                                 maxField = (oldField[i110] > maxField) ? oldField[i110] : maxField;
467
468                                 minField = (oldField[i001] < minField) ? oldField[i001] : minField;
469                                 maxField = (oldField[i001] > maxField) ? oldField[i001] : maxField;
470
471                                 minField = (oldField[i011] < minField) ? oldField[i011] : minField;
472                                 maxField = (oldField[i011] > maxField) ? oldField[i011] : maxField;
473
474                                 minField = (oldField[i101] < minField) ? oldField[i101] : minField;
475                                 maxField = (oldField[i101] > maxField) ? oldField[i101] : maxField;
476
477                                 minField = (oldField[i111] < minField) ? oldField[i111] : minField;
478                                 maxField = (oldField[i111] > maxField) ? oldField[i111] : maxField;
479
480                                 newField[index] = (newField[index] > maxField) ? maxField : newField[index];
481                                 newField[index] = (newField[index] < minField) ? minField : newField[index];
482                         }
483 }
484
485 //////////////////////////////////////////////////////////////////////
486 // Reverts any backtraces that go into boundaries back to first 
487 // order -- in this case the error correction term was totally
488 // incorrect
489 //////////////////////////////////////////////////////////////////////
490 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
491                                 float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection)
492 {
493         const int sx= res[0];
494         const int sy= res[1];
495         const int sz= res[2];
496         const int slabSize = res[0] * res[1];
497
498         for (int z = 1; z < sz-1; z++)
499                 for (int y = 1; y < sy-1; y++)
500                         for (int x = 1; x < sx-1; x++)
501                         {
502                                 const int index = x + y * sx+ z * slabSize;
503                                 // backtrace
504                                 float xBackward = x + dt * velx[index];
505                                 float yBackward = y + dt * vely[index];
506                                 float zBackward = z + dt * velz[index];
507                                 float xTrace    = x - dt * velx[index];
508                                 float yTrace    = y - dt * vely[index];
509                                 float zTrace    = z - dt * velz[index];
510
511                                 // see if it goes outside the boundaries
512                                 bool hasObstacle = 
513                                         (zTrace < 1.0f)    || (zTrace > sz - 2.0f) ||
514                                         (yTrace < 1.0f)    || (yTrace > sy - 2.0f) ||
515                                         (xTrace < 1.0f)    || (xTrace > sx - 2.0f) ||
516                                         (zBackward < 1.0f) || (zBackward > sz - 2.0f) ||
517                                         (yBackward < 1.0f) || (yBackward > sy - 2.0f) ||
518                                         (xBackward < 1.0f) || (xBackward > sx - 2.0f);
519                                 // reuse old advection instead of doing another one...
520                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
521
522                                 // clamp to prevent an out of bounds access when looking into
523                                 // the _obstacles array
524                                 zTrace = (zTrace < 0.5f) ? 0.5f : zTrace;
525                                 zTrace = (zTrace > sz - 1.5f) ? sz - 1.5f : zTrace;
526                                 yTrace = (yTrace < 0.5f) ? 0.5f : yTrace;
527                                 yTrace = (yTrace > sy - 1.5f) ? sy - 1.5f : yTrace;
528                                 xTrace = (xTrace < 0.5f) ? 0.5f : xTrace;
529                                 xTrace = (xTrace > sx - 1.5f) ? sx - 1.5f : xTrace;
530
531                                 // locate neighbors to interpolate,
532                                 // do backward first since we will use the forward indices if a
533                                 // reversion is actually necessary
534                                 zBackward = (zBackward < 0.5f) ? 0.5f : zBackward;
535                                 zBackward = (zBackward > sz - 1.5f) ? sz - 1.5f : zBackward;
536                                 yBackward = (yBackward < 0.5f) ? 0.5f : yBackward;
537                                 yBackward = (yBackward > sy - 1.5f) ? sy - 1.5f : yBackward;
538                                 xBackward = (xBackward < 0.5f) ? 0.5f : xBackward;
539                                 xBackward = (xBackward > sx - 1.5f) ? sx - 1.5f : xBackward;
540
541                                 int x0 = (int)xBackward;
542                                 int x1 = x0 + 1;
543                                 int y0 = (int)yBackward;
544                                 int y1 = y0 + 1;
545                                 int z0 = (int)zBackward;
546                                 int z1 = z0 + 1;
547                                 if(obstacles && !hasObstacle) {
548                                         hasObstacle = hasObstacle || 
549                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
550                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
551                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
552                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
553                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
554                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
555                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
556                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
557                                 }
558                                 // reuse old advection instead of doing another one...
559                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
560
561                                 x0 = (int)xTrace;
562                                 x1 = x0 + 1;
563                                 y0 = (int)yTrace;
564                                 y1 = y0 + 1;
565                                 z0 = (int)zTrace;
566                                 z1 = z0 + 1;
567                                 if(obstacles && !hasObstacle) {
568                                         hasObstacle = hasObstacle || 
569                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
570                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
571                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
572                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
573                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
574                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
575                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
576                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
577                                 } // obstacle array
578                                 // reuse old advection instead of doing another one...
579                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
580
581                                 // see if either the forward or backward ray went into
582                                 // a boundary
583                                 if (hasObstacle) {
584                                         // get interpolation weights
585                                         float s1 = xTrace - x0;
586                                         float s0 = 1.0f - s1;
587                                         float t1 = yTrace - y0;
588                                         float t0 = 1.0f - t1;
589                                         float u1 = zTrace - z0;
590                                         float u0 = 1.0f - u1;
591
592                                         const int i000 = x0 + y0 * sx + z0 * slabSize;
593                                         const int i010 = x0 + y1 * sx + z0 * slabSize;
594                                         const int i100 = x1 + y0 * sx + z0 * slabSize;
595                                         const int i110 = x1 + y1 * sx + z0 * slabSize;
596                                         const int i001 = x0 + y0 * sx + z1 * slabSize;
597                                         const int i011 = x0 + y1 * sx + z1 * slabSize;
598                                         const int i101 = x1 + y0 * sx + z1 * slabSize;
599                                         const int i111 = x1 + y1 * sx + z1 * slabSize;
600
601                                         // interpolate, (indices could be computed once)
602                                         newField[index] = u0 * (s0 * (
603                                                                 t0 * oldField[i000] +
604                                                                 t1 * oldField[i010]) +
605                                                         s1 * (t0 * oldField[i100] +
606                                                                 t1 * oldField[i110])) +
607                                                 u1 * (s0 * (t0 * oldField[i001] +
608                                                                         t1 * oldField[i011]) +
609                                                                 s1 * (t0 * oldField[i101] +
610                                                                         t1 * oldField[i111])); 
611                                 }
612                         } // xyz
613 }
614