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