ac485ad983a005d92748d9c5e94d7a28e6a87721
[blender-staging.git] / intern / smoke / intern / FLUID_3D_STATIC.cpp
1 /** \file smoke/intern/FLUID_3D_STATIC.cpp
2  *  \ingroup smoke
3  */
4 //////////////////////////////////////////////////////////////////////
5 // This file is part of Wavelet Turbulence.
6 //
7 // Wavelet Turbulence is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // Wavelet Turbulence is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
19 //
20 // Copyright 2008 Theodore Kim and Nils Thuerey
21 //
22 // FLUID_3D.cpp: implementation of the static functions of the FLUID_3D class.
23 //
24 //////////////////////////////////////////////////////////////////////
25 // Heavy parallel optimization done. Many of the old functions now
26 // take begin and end parameters and process only specified part of the data.
27 // Some functions were divided into multiple ones.
28 //              - MiikaH
29 //////////////////////////////////////////////////////////////////////
30
31 #include <zlib.h>
32 #include "FLUID_3D.h"
33 #include "IMAGE.h"
34 #include "WTURBULENCE.h"
35 #include "INTERPOLATE.h"
36
37 //////////////////////////////////////////////////////////////////////
38 // add a test cube of density to the center
39 //////////////////////////////////////////////////////////////////////
40 /*
41 void FLUID_3D::addSmokeColumn() {
42         addSmokeTestCase(_density, _res, 1.0);
43         // addSmokeTestCase(_zVelocity, _res, 1.0);
44         addSmokeTestCase(_heat, _res, 1.0);
45         if (_wTurbulence) {
46                 addSmokeTestCase(_wTurbulence->getDensityBig(), _wTurbulence->getResBig(), 1.0);
47         }
48 }
49 */
50
51 //////////////////////////////////////////////////////////////////////
52 // generic static version, so that it can be applied to the
53 // WTURBULENCE grid as well
54 //////////////////////////////////////////////////////////////////////
55
56 void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res)
57 {
58         const int slabSize = res[0]*res[1]; int maxRes = (int)MAX3V(res);
59         float dx = 1.0f / (float)maxRes;
60
61         float xTotal = dx * res[0];
62         float yTotal = dx * res[1];
63
64         float heighMin = 0.05;
65         float heighMax = 0.10;
66
67         for (int y = 0; y < res[2]; y++)
68                 for (int z = (int)(heighMin*res[2]); z <= (int)(heighMax * res[2]); z++)
69                         for (int x = 0; x < res[0]; x++) {
70                                 float xLength = x * dx - xTotal * 0.4f;
71                                 float yLength = y * dx - yTotal * 0.5f;
72                                 float radius = sqrtf(xLength * xLength + yLength * yLength);
73
74                                 if (radius < 0.075f * xTotal) {
75                                         int index = x + y * res[0] + z * slabSize;
76                                         field[index] = 1.0f;
77                                 }
78                         }
79 }
80
81
82 //////////////////////////////////////////////////////////////////////
83 // set x direction to Neumann boundary conditions
84 //////////////////////////////////////////////////////////////////////
85 void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd)
86 {
87         const int slabSize = res[0] * res[1];
88         int index;
89         for (int z = zBegin; z < zEnd; z++)
90                 for (int y = 0; y < res[1]; y++)
91                 {
92                         // left slab
93                         index = y * res[0] + z * slabSize;
94                         field[index] = field[index + 2];
95
96                         // right slab
97                         index = y * res[0] + z * slabSize + res[0] - 1;
98                         field[index] = field[index - 2];
99                 }
100  }
101
102 //////////////////////////////////////////////////////////////////////
103 // set y direction to Neumann boundary conditions
104 //////////////////////////////////////////////////////////////////////
105 void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd)
106 {
107         const int slabSize = res[0] * res[1];
108         int index;
109         for (int z = zBegin; z < zEnd; z++)
110                 for (int x = 0; x < res[0]; x++)
111                 {
112                         // front slab
113                         index = x + z * slabSize;
114                         field[index] = field[index + 2 * res[0]];
115
116                         // back slab
117                         index = x + z * slabSize + slabSize - res[0];
118                         field[index] = field[index - 2 * res[0]];
119                 }
120 }
121
122 //////////////////////////////////////////////////////////////////////
123 // set z direction to Neumann boundary conditions
124 //////////////////////////////////////////////////////////////////////
125 void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd)
126 {
127         const int slabSize = res[0] * res[1];
128         const int totalCells = res[0] * res[1] * res[2];
129         const int cellsslab = totalCells - slabSize;
130         int index;
131
132         if (zBegin == 0) {
133                 for (int y = 0; y < res[1]; y++)
134                         for (int x = 0; x < res[0]; x++)
135                         {
136                                 // front slab
137                                 index = x + y * res[0];
138                                 field[index] = field[index + 2 * slabSize];
139                         }
140         }
141
142         if (zEnd == res[2]) {
143                 for (int y = 0; y < res[1]; y++)
144                         for (int x = 0; x < res[0]; x++)
145                         {
146                                 // back slab
147                                 index = x + y * res[0] + cellsslab;
148                                 field[index] = field[index - 2 * slabSize];
149                         }
150         }
151                 
152 }
153
154 //////////////////////////////////////////////////////////////////////
155 // set x direction to zero
156 //////////////////////////////////////////////////////////////////////
157 void FLUID_3D::setZeroX(float* field, Vec3Int res, int zBegin, int zEnd)
158 {
159         const int slabSize = res[0] * res[1];
160         int index;
161         for (int z = zBegin; z < zEnd; z++)
162                 for (int y = 0; y < res[1]; y++)
163                 {
164                         // left slab
165                         index = y * res[0] + z * slabSize;
166                         field[index] = 0.0f;
167
168                         // right slab
169                         index += res[0] - 1;
170                         field[index] = 0.0f;
171                 }
172 }
173
174 //////////////////////////////////////////////////////////////////////
175 // set y direction to zero
176 //////////////////////////////////////////////////////////////////////
177 void FLUID_3D::setZeroY(float* field, Vec3Int res, int zBegin, int zEnd)
178 {
179         const int slabSize = res[0] * res[1];
180         int index;
181         for (int z = zBegin; z < zEnd; z++)
182                 for (int x = 0; x < res[0]; x++)
183                 {
184                         // bottom slab
185                         index = x + z * slabSize;
186                         field[index] = 0.0f;
187
188                         // top slab
189                         index += slabSize - res[0];
190                         field[index] = 0.0f;
191                 }
192 }
193
194 //////////////////////////////////////////////////////////////////////
195 // set z direction to zero
196 //////////////////////////////////////////////////////////////////////
197 void FLUID_3D::setZeroZ(float* field, Vec3Int res, int zBegin, int zEnd)
198 {
199         const int slabSize = res[0] * res[1];
200         const int totalCells = res[0] * res[1] * res[2];
201
202         int index = 0;
203         if ((zBegin == 0))
204         for (int y = 0; y < res[1]; y++)
205                 for (int x = 0; x < res[0]; x++, index++)
206                 {
207                         // front slab
208                         field[index] = 0.0f;
209     }
210
211         if (zEnd == res[2])
212         {
213                 index=0;
214                 int indexx=0;
215                 const int cellsslab = totalCells - slabSize;
216
217                 for (int y = 0; y < res[1]; y++)
218                         for (int x = 0; x < res[0]; x++, index++)
219                         {
220
221                                 // back slab
222                                 indexx = index + cellsslab;
223                                 field[indexx] = 0.0f;
224                         }
225         }
226  }
227 //////////////////////////////////////////////////////////////////////
228 // copy grid boundary
229 //////////////////////////////////////////////////////////////////////
230 void FLUID_3D::copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd)
231 {
232         const int slabSize = res[0] * res[1];
233         int index;
234         for (int z = zBegin; z < zEnd; z++)
235                 for (int y = 0; y < res[1]; y++)
236                 {
237                         // left slab
238                         index = y * res[0] + z * slabSize;
239                         field[index] = field[index + 1];
240
241                         // right slab
242                         index += res[0] - 1;
243                         field[index] = field[index - 1];
244                 }
245 }
246 void FLUID_3D::copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd)
247 {
248         const int slabSize = res[0] * res[1];
249         //const int totalCells = res[0] * res[1] * res[2];
250         int index;
251         for (int z = zBegin; z < zEnd; z++)
252                 for (int x = 0; x < res[0]; x++)
253                 {
254                         // bottom slab
255                         index = x + z * slabSize;
256                         field[index] = field[index + res[0]]; 
257                         // top slab
258                         index += slabSize - res[0];
259                         field[index] = field[index - res[0]];
260                 }
261 }
262 void FLUID_3D::copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd)
263 {
264         const int slabSize = res[0] * res[1];
265         const int totalCells = res[0] * res[1] * res[2];
266         int index=0;
267
268         if ((zBegin == 0))
269         for (int y = 0; y < res[1]; y++)
270                 for (int x = 0; x < res[0]; x++, index++)
271                 {
272                         field[index] = field[index + slabSize]; 
273                 }
274
275         if ((zEnd == res[2]))
276         {
277
278         index=0;
279         int indexx=0;
280         const int cellsslab = totalCells - slabSize;
281
282         for (int y = 0; y < res[1]; y++)
283                 for (int x = 0; x < res[0]; x++, index++)
284                 {
285                         // back slab
286                         indexx = index + cellsslab;
287                         field[indexx] = field[indexx - slabSize];
288                 }
289         }
290 }
291
292 /////////////////////////////////////////////////////////////////////
293 // advect field with the semi lagrangian method
294 //////////////////////////////////////////////////////////////////////
295 void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
296                 float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
297 {
298         const int xres = res[0];
299         const int yres = res[1];
300         const int zres = res[2];
301         const int slabSize = res[0] * res[1];
302
303
304         for (int z = zBegin; z < zEnd; z++)
305                 for (int y = 0; y < yres; y++)
306                         for (int x = 0; x < xres; x++)
307                         {
308                                 const int index = x + y * xres + z * xres*yres;
309                                 
310         // backtrace
311                                 float xTrace = x - dt * velx[index];
312                                 float yTrace = y - dt * vely[index];
313                                 float zTrace = z - dt * velz[index];
314
315                                 // clamp backtrace to grid boundaries
316                                 if (xTrace < 0.5) xTrace = 0.5;
317                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
318                                 if (yTrace < 0.5) yTrace = 0.5;
319                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
320                                 if (zTrace < 0.5) zTrace = 0.5;
321                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
322
323                                 // locate neighbors to interpolate
324                                 const int x0 = (int)xTrace;
325                                 const int x1 = x0 + 1;
326                                 const int y0 = (int)yTrace;
327                                 const int y1 = y0 + 1;
328                                 const int z0 = (int)zTrace;
329                                 const int z1 = z0 + 1;
330
331                                 // get interpolation weights
332                                 const float s1 = xTrace - x0;
333                                 const float s0 = 1.0f - s1;
334                                 const float t1 = yTrace - y0;
335                                 const float t0 = 1.0f - t1;
336                                 const float u1 = zTrace - z0;
337                                 const float u0 = 1.0f - u1;
338
339                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
340                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
341                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
342                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
343                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
344                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
345                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
346                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
347
348                                 // interpolate
349                                 // (indices could be computed once)
350                                 newField[index] = u0 * (s0 * (t0 * oldField[i000] +
351                                                         t1 * oldField[i010]) +
352                                                 s1 * (t0 * oldField[i100] +
353                                                         t1 * oldField[i110])) +
354                                         u1 * (s0 * (t0 * oldField[i001] +
355                                                                 t1 * oldField[i011]) +
356                                                         s1 * (t0 * oldField[i101] +
357                                                                 t1 * oldField[i111]));
358                         }
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::advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
368                                 float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd)
369 {
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];*/       // not needed as all the values are written first
376
377         float*& phiN    = oldField;
378         float*& phiN1   = tempResult;
379
380
381
382         // phiHatN1 = A(phiN)
383         advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiN1, res, zBegin, zEnd);         // uses wide data from old field and velocities (both are whole)
384 }
385
386
387
388 void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
389                                 float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd)
390 {
391         float* phiHatN  = tempResult;
392         float* t1  = temp1;
393         const int sx= res[0];
394         const int sy= res[1];
395
396         float*& phiN    = oldField;
397         float*& phiN1   = newField;
398
399
400
401         // phiHatN = A^R(phiHatN1)
402         advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN, t1, res, zBegin, zEnd);             // uses wide data from old field and velocities (both are whole)
403
404         // phiN1 = phiHatN1 + (phiN - phiHatN) / 2
405         const int border = 0; 
406         for (int z = zBegin+border; z < zEnd-border; z++)
407                 for (int y = border; y < sy-border; y++)
408                         for (int x = border; x < sx-border; x++) {
409                                 int index = x + y * sx + z * sx*sy;
410                                 phiN1[index] = phiHatN[index] + (phiN[index] - t1[index]) * 0.50f;
411                                 //phiN1[index] = phiHatN1[index]; // debug, correction off
412                         }
413         copyBorderX(phiN1, res, zBegin, zEnd);
414         copyBorderY(phiN1, res, zBegin, zEnd);
415         copyBorderZ(phiN1, res, zBegin, zEnd);
416
417         // clamp any newly created extrema
418         clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd);               // uses wide data from old field and velocities (both are whole)
419
420         // if the error estimate was bad, revert to first order
421         clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd);       // phiHatN is only used at cells within thread range, so its ok
422
423
424
425
426 //////////////////////////////////////////////////////////////////////
427 // Clamp the extrema generated by the BFECC error correction
428 //////////////////////////////////////////////////////////////////////
429 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
430                 float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
431 {
432         const int xres= res[0];
433         const int yres= res[1];
434         const int zres= res[2];
435         const int slabSize = res[0] * res[1];
436
437         int bb=0;
438         int bt=0;
439
440         if (zBegin == 0) {bb = 1;}
441         if (zEnd == res[2]) {bt = 1;}
442
443
444         for (int z = zBegin+bb; z < zEnd-bt; z++)
445                 for (int y = 1; y < yres-1; y++)
446                         for (int x = 1; x < xres-1; x++)
447                         {
448                                 const int index = x + y * xres+ z * xres*yres;
449                                 // backtrace
450                                 float xTrace = x - dt * velx[index];
451                                 float yTrace = y - dt * vely[index];
452                                 float zTrace = z - dt * velz[index];
453
454                                 // clamp backtrace to grid boundaries
455                                 if (xTrace < 0.5) xTrace = 0.5;
456                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
457                                 if (yTrace < 0.5) yTrace = 0.5;
458                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
459                                 if (zTrace < 0.5) zTrace = 0.5;
460                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
461
462                                 // locate neighbors to interpolate
463                                 const int x0 = (int)xTrace;
464                                 const int x1 = x0 + 1;
465                                 const int y0 = (int)yTrace;
466                                 const int y1 = y0 + 1;
467                                 const int z0 = (int)zTrace;
468                                 const int z1 = z0 + 1;
469
470                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
471                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
472                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
473                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
474                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
475                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
476                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
477                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
478
479                                 float minField = oldField[i000];
480                                 float maxField = oldField[i000];
481
482                                 minField = (oldField[i010] < minField) ? oldField[i010] : minField;
483                                 maxField = (oldField[i010] > maxField) ? oldField[i010] : maxField;
484
485                                 minField = (oldField[i100] < minField) ? oldField[i100] : minField;
486                                 maxField = (oldField[i100] > maxField) ? oldField[i100] : maxField;
487
488                                 minField = (oldField[i110] < minField) ? oldField[i110] : minField;
489                                 maxField = (oldField[i110] > maxField) ? oldField[i110] : maxField;
490
491                                 minField = (oldField[i001] < minField) ? oldField[i001] : minField;
492                                 maxField = (oldField[i001] > maxField) ? oldField[i001] : maxField;
493
494                                 minField = (oldField[i011] < minField) ? oldField[i011] : minField;
495                                 maxField = (oldField[i011] > maxField) ? oldField[i011] : maxField;
496
497                                 minField = (oldField[i101] < minField) ? oldField[i101] : minField;
498                                 maxField = (oldField[i101] > maxField) ? oldField[i101] : maxField;
499
500                                 minField = (oldField[i111] < minField) ? oldField[i111] : minField;
501                                 maxField = (oldField[i111] > maxField) ? oldField[i111] : maxField;
502
503                                 newField[index] = (newField[index] > maxField) ? maxField : newField[index];
504                                 newField[index] = (newField[index] < minField) ? minField : newField[index];
505                         }
506 }
507
508 //////////////////////////////////////////////////////////////////////
509 // Reverts any backtraces that go into boundaries back to first 
510 // order -- in this case the error correction term was totally
511 // incorrect
512 //////////////////////////////////////////////////////////////////////
513 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
514                                 float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd)
515 {
516         const int sx= res[0];
517         const int sy= res[1];
518         const int sz= res[2];
519         const int slabSize = res[0] * res[1];
520
521         int bb=0;
522         int bt=0;
523
524         if (zBegin == 0) {bb = 1;}
525         if (zEnd == res[2]) {bt = 1;}
526
527         for (int z = zBegin+bb; z < zEnd-bt; z++)
528                 for (int y = 1; y < sy-1; y++)
529                         for (int x = 1; x < sx-1; x++)
530                         {
531                                 const int index = x + y * sx+ z * slabSize;
532                                 // backtrace
533                                 float xBackward = x + dt * velx[index];
534                                 float yBackward = y + dt * vely[index];
535                                 float zBackward = z + dt * velz[index];
536                                 float xTrace    = x - dt * velx[index];
537                                 float yTrace    = y - dt * vely[index];
538                                 float zTrace    = z - dt * velz[index];
539
540                                 // see if it goes outside the boundaries
541                                 bool hasObstacle = 
542                                         (zTrace < 1.0f)    || (zTrace > sz - 2.0f) ||
543                                         (yTrace < 1.0f)    || (yTrace > sy - 2.0f) ||
544                                         (xTrace < 1.0f)    || (xTrace > sx - 2.0f) ||
545                                         (zBackward < 1.0f) || (zBackward > sz - 2.0f) ||
546                                         (yBackward < 1.0f) || (yBackward > sy - 2.0f) ||
547                                         (xBackward < 1.0f) || (xBackward > sx - 2.0f);
548                                 // reuse old advection instead of doing another one...
549                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
550
551                                 // clamp to prevent an out of bounds access when looking into
552                                 // the _obstacles array
553                                 zTrace = (zTrace < 0.5f) ? 0.5f : zTrace;
554                                 zTrace = (zTrace > sz - 1.5f) ? sz - 1.5f : zTrace;
555                                 yTrace = (yTrace < 0.5f) ? 0.5f : yTrace;
556                                 yTrace = (yTrace > sy - 1.5f) ? sy - 1.5f : yTrace;
557                                 xTrace = (xTrace < 0.5f) ? 0.5f : xTrace;
558                                 xTrace = (xTrace > sx - 1.5f) ? sx - 1.5f : xTrace;
559
560                                 // locate neighbors to interpolate,
561                                 // do backward first since we will use the forward indices if a
562                                 // reversion is actually necessary
563                                 zBackward = (zBackward < 0.5f) ? 0.5f : zBackward;
564                                 zBackward = (zBackward > sz - 1.5f) ? sz - 1.5f : zBackward;
565                                 yBackward = (yBackward < 0.5f) ? 0.5f : yBackward;
566                                 yBackward = (yBackward > sy - 1.5f) ? sy - 1.5f : yBackward;
567                                 xBackward = (xBackward < 0.5f) ? 0.5f : xBackward;
568                                 xBackward = (xBackward > sx - 1.5f) ? sx - 1.5f : xBackward;
569
570                                 int x0 = (int)xBackward;
571                                 int x1 = x0 + 1;
572                                 int y0 = (int)yBackward;
573                                 int y1 = y0 + 1;
574                                 int z0 = (int)zBackward;
575                                 int z1 = z0 + 1;
576                                 if(obstacles && !hasObstacle) {
577                                         hasObstacle = hasObstacle || 
578                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
579                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
580                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
581                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
582                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
583                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
584                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
585                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
586                                 }
587                                 // reuse old advection instead of doing another one...
588                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
589
590                                 x0 = (int)xTrace;
591                                 x1 = x0 + 1;
592                                 y0 = (int)yTrace;
593                                 y1 = y0 + 1;
594                                 z0 = (int)zTrace;
595                                 z1 = z0 + 1;
596                                 if(obstacles && !hasObstacle) {
597                                         hasObstacle = hasObstacle || 
598                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
599                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
600                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
601                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
602                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
603                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
604                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
605                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
606                                 } // obstacle array
607                                 // reuse old advection instead of doing another one...
608                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
609
610                                 // see if either the forward or backward ray went into
611                                 // a boundary
612                                 if (hasObstacle) {
613                                         // get interpolation weights
614                                         float s1 = xTrace - x0;
615                                         float s0 = 1.0f - s1;
616                                         float t1 = yTrace - y0;
617                                         float t0 = 1.0f - t1;
618                                         float u1 = zTrace - z0;
619                                         float u0 = 1.0f - u1;
620
621                                         const int i000 = x0 + y0 * sx + z0 * slabSize;
622                                         const int i010 = x0 + y1 * sx + z0 * slabSize;
623                                         const int i100 = x1 + y0 * sx + z0 * slabSize;
624                                         const int i110 = x1 + y1 * sx + z0 * slabSize;
625                                         const int i001 = x0 + y0 * sx + z1 * slabSize;
626                                         const int i011 = x0 + y1 * sx + z1 * slabSize;
627                                         const int i101 = x1 + y0 * sx + z1 * slabSize;
628                                         const int i111 = x1 + y1 * sx + z1 * slabSize;
629
630                                         // interpolate, (indices could be computed once)
631                                         newField[index] = u0 * (s0 * (
632                                                                 t0 * oldField[i000] +
633                                                                 t1 * oldField[i010]) +
634                                                         s1 * (t0 * oldField[i100] +
635                                                                 t1 * oldField[i110])) +
636                                                 u1 * (s0 * (t0 * oldField[i001] +
637                                                                         t1 * oldField[i011]) +
638                                                                 s1 * (t0 * oldField[i101] +
639                                                                         t1 * oldField[i111])); 
640                                 }
641                         } // xyz
642 }