Fix various warnings with clang build, and adjust cmake clang warnings flags
[blender.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         {
205                 for (int y = 0; y < res[1]; y++)
206                         for (int x = 0; x < res[0]; x++, index++)
207                         {
208                                 // front slab
209                                 field[index] = 0.0f;
210                 }
211         }
212
213         if (zEnd == res[2])
214         {
215                 index=0;
216                 int indexx=0;
217                 const int cellsslab = totalCells - slabSize;
218
219                 for (int y = 0; y < res[1]; y++)
220                         for (int x = 0; x < res[0]; x++, index++)
221                         {
222
223                                 // back slab
224                                 indexx = index + cellsslab;
225                                 field[indexx] = 0.0f;
226                         }
227         }
228  }
229 //////////////////////////////////////////////////////////////////////
230 // copy grid boundary
231 //////////////////////////////////////////////////////////////////////
232 void FLUID_3D::copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd)
233 {
234         const int slabSize = res[0] * res[1];
235         int index;
236         for (int z = zBegin; z < zEnd; z++)
237                 for (int y = 0; y < res[1]; y++)
238                 {
239                         // left slab
240                         index = y * res[0] + z * slabSize;
241                         field[index] = field[index + 1];
242
243                         // right slab
244                         index += res[0] - 1;
245                         field[index] = field[index - 1];
246                 }
247 }
248 void FLUID_3D::copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd)
249 {
250         const int slabSize = res[0] * res[1];
251         //const int totalCells = res[0] * res[1] * res[2];
252         int index;
253         for (int z = zBegin; z < zEnd; z++)
254                 for (int x = 0; x < res[0]; x++)
255                 {
256                         // bottom slab
257                         index = x + z * slabSize;
258                         field[index] = field[index + res[0]]; 
259                         // top slab
260                         index += slabSize - res[0];
261                         field[index] = field[index - res[0]];
262                 }
263 }
264 void FLUID_3D::copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd)
265 {
266         const int slabSize = res[0] * res[1];
267         const int totalCells = res[0] * res[1] * res[2];
268         int index=0;
269
270         if (zBegin == 0)
271         {
272                 for (int y = 0; y < res[1]; y++)
273                         for (int x = 0; x < res[0]; x++, index++)
274                         {
275                                 field[index] = field[index + slabSize]; 
276                         }
277         }
278
279         if (zEnd == res[2])
280         {
281
282                 index=0;
283                 int indexx=0;
284                 const int cellsslab = totalCells - slabSize;
285
286                 for (int y = 0; y < res[1]; y++)
287                         for (int x = 0; x < res[0]; x++, index++)
288                         {
289                                 // back slab
290                                 indexx = index + cellsslab;
291                                 field[indexx] = field[indexx - slabSize];
292                         }
293         }
294 }
295
296 /////////////////////////////////////////////////////////////////////
297 // advect field with the semi lagrangian method
298 //////////////////////////////////////////////////////////////////////
299 void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
300                 float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
301 {
302         const int xres = res[0];
303         const int yres = res[1];
304         const int zres = res[2];
305         const int slabSize = res[0] * res[1];
306
307
308         for (int z = zBegin; z < zEnd; z++)
309                 for (int y = 0; y < yres; y++)
310                         for (int x = 0; x < xres; x++)
311                         {
312                                 const int index = x + y * xres + z * xres*yres;
313                                 
314         // backtrace
315                                 float xTrace = x - dt * velx[index];
316                                 float yTrace = y - dt * vely[index];
317                                 float zTrace = z - dt * velz[index];
318
319                                 // clamp backtrace to grid boundaries
320                                 if (xTrace < 0.5) xTrace = 0.5;
321                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
322                                 if (yTrace < 0.5) yTrace = 0.5;
323                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
324                                 if (zTrace < 0.5) zTrace = 0.5;
325                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
326
327                                 // locate neighbors to interpolate
328                                 const int x0 = (int)xTrace;
329                                 const int x1 = x0 + 1;
330                                 const int y0 = (int)yTrace;
331                                 const int y1 = y0 + 1;
332                                 const int z0 = (int)zTrace;
333                                 const int z1 = z0 + 1;
334
335                                 // get interpolation weights
336                                 const float s1 = xTrace - x0;
337                                 const float s0 = 1.0f - s1;
338                                 const float t1 = yTrace - y0;
339                                 const float t0 = 1.0f - t1;
340                                 const float u1 = zTrace - z0;
341                                 const float u0 = 1.0f - u1;
342
343                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
344                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
345                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
346                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
347                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
348                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
349                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
350                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
351
352                                 // interpolate
353                                 // (indices could be computed once)
354                                 newField[index] = u0 * (s0 * (t0 * oldField[i000] +
355                                                         t1 * oldField[i010]) +
356                                                 s1 * (t0 * oldField[i100] +
357                                                         t1 * oldField[i110])) +
358                                         u1 * (s0 * (t0 * oldField[i001] +
359                                                                 t1 * oldField[i011]) +
360                                                         s1 * (t0 * oldField[i101] +
361                                                                 t1 * oldField[i111]));
362                         }
363 }
364
365
366 /////////////////////////////////////////////////////////////////////
367 // advect field with the maccormack method
368 //
369 // comments are the pseudocode from selle's paper
370 //////////////////////////////////////////////////////////////////////
371 void FLUID_3D::advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
372                                 float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd)
373 {
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];*/       // not needed as all the values are written first
380
381         float*& phiN    = oldField;
382         float*& phiN1   = tempResult;
383
384
385
386         // phiHatN1 = A(phiN)
387         advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiN1, res, zBegin, zEnd);         // uses wide data from old field and velocities (both are whole)
388 }
389
390
391
392 void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
393                                 float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd)
394 {
395         float* phiHatN  = tempResult;
396         float* t1  = temp1;
397         const int sx= res[0];
398         const int sy= res[1];
399
400         float*& phiN    = oldField;
401         float*& phiN1   = newField;
402
403
404
405         // phiHatN = A^R(phiHatN1)
406         advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN, t1, res, zBegin, zEnd);             // uses wide data from old field and velocities (both are whole)
407
408         // phiN1 = phiHatN1 + (phiN - phiHatN) / 2
409         const int border = 0; 
410         for (int z = zBegin+border; z < zEnd-border; z++)
411                 for (int y = border; y < sy-border; y++)
412                         for (int x = border; x < sx-border; x++) {
413                                 int index = x + y * sx + z * sx*sy;
414                                 phiN1[index] = phiHatN[index] + (phiN[index] - t1[index]) * 0.50f;
415                                 //phiN1[index] = phiHatN1[index]; // debug, correction off
416                         }
417         copyBorderX(phiN1, res, zBegin, zEnd);
418         copyBorderY(phiN1, res, zBegin, zEnd);
419         copyBorderZ(phiN1, res, zBegin, zEnd);
420
421         // clamp any newly created extrema
422         clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd);               // uses wide data from old field and velocities (both are whole)
423
424         // if the error estimate was bad, revert to first order
425         clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd);       // phiHatN is only used at cells within thread range, so its ok
426
427
428
429
430 //////////////////////////////////////////////////////////////////////
431 // Clamp the extrema generated by the BFECC error correction
432 //////////////////////////////////////////////////////////////////////
433 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
434                 float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
435 {
436         const int xres= res[0];
437         const int yres= res[1];
438         const int zres= res[2];
439         const int slabSize = res[0] * res[1];
440
441         int bb=0;
442         int bt=0;
443
444         if (zBegin == 0) {bb = 1;}
445         if (zEnd == res[2]) {bt = 1;}
446
447
448         for (int z = zBegin+bb; z < zEnd-bt; z++)
449                 for (int y = 1; y < yres-1; y++)
450                         for (int x = 1; x < xres-1; x++)
451                         {
452                                 const int index = x + y * xres+ z * xres*yres;
453                                 // backtrace
454                                 float xTrace = x - dt * velx[index];
455                                 float yTrace = y - dt * vely[index];
456                                 float zTrace = z - dt * velz[index];
457
458                                 // clamp backtrace to grid boundaries
459                                 if (xTrace < 0.5) xTrace = 0.5;
460                                 if (xTrace > xres - 1.5) xTrace = xres - 1.5;
461                                 if (yTrace < 0.5) yTrace = 0.5;
462                                 if (yTrace > yres - 1.5) yTrace = yres - 1.5;
463                                 if (zTrace < 0.5) zTrace = 0.5;
464                                 if (zTrace > zres - 1.5) zTrace = zres - 1.5;
465
466                                 // locate neighbors to interpolate
467                                 const int x0 = (int)xTrace;
468                                 const int x1 = x0 + 1;
469                                 const int y0 = (int)yTrace;
470                                 const int y1 = y0 + 1;
471                                 const int z0 = (int)zTrace;
472                                 const int z1 = z0 + 1;
473
474                                 const int i000 = x0 + y0 * xres + z0 * slabSize;
475                                 const int i010 = x0 + y1 * xres + z0 * slabSize;
476                                 const int i100 = x1 + y0 * xres + z0 * slabSize;
477                                 const int i110 = x1 + y1 * xres + z0 * slabSize;
478                                 const int i001 = x0 + y0 * xres + z1 * slabSize;
479                                 const int i011 = x0 + y1 * xres + z1 * slabSize;
480                                 const int i101 = x1 + y0 * xres + z1 * slabSize;
481                                 const int i111 = x1 + y1 * xres + z1 * slabSize;
482
483                                 float minField = oldField[i000];
484                                 float maxField = oldField[i000];
485
486                                 minField = (oldField[i010] < minField) ? oldField[i010] : minField;
487                                 maxField = (oldField[i010] > maxField) ? oldField[i010] : maxField;
488
489                                 minField = (oldField[i100] < minField) ? oldField[i100] : minField;
490                                 maxField = (oldField[i100] > maxField) ? oldField[i100] : maxField;
491
492                                 minField = (oldField[i110] < minField) ? oldField[i110] : minField;
493                                 maxField = (oldField[i110] > maxField) ? oldField[i110] : maxField;
494
495                                 minField = (oldField[i001] < minField) ? oldField[i001] : minField;
496                                 maxField = (oldField[i001] > maxField) ? oldField[i001] : maxField;
497
498                                 minField = (oldField[i011] < minField) ? oldField[i011] : minField;
499                                 maxField = (oldField[i011] > maxField) ? oldField[i011] : maxField;
500
501                                 minField = (oldField[i101] < minField) ? oldField[i101] : minField;
502                                 maxField = (oldField[i101] > maxField) ? oldField[i101] : maxField;
503
504                                 minField = (oldField[i111] < minField) ? oldField[i111] : minField;
505                                 maxField = (oldField[i111] > maxField) ? oldField[i111] : maxField;
506
507                                 newField[index] = (newField[index] > maxField) ? maxField : newField[index];
508                                 newField[index] = (newField[index] < minField) ? minField : newField[index];
509                         }
510 }
511
512 //////////////////////////////////////////////////////////////////////
513 // Reverts any backtraces that go into boundaries back to first 
514 // order -- in this case the error correction term was totally
515 // incorrect
516 //////////////////////////////////////////////////////////////////////
517 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
518                                 float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd)
519 {
520         const int sx= res[0];
521         const int sy= res[1];
522         const int sz= res[2];
523         const int slabSize = res[0] * res[1];
524
525         int bb=0;
526         int bt=0;
527
528         if (zBegin == 0) {bb = 1;}
529         if (zEnd == res[2]) {bt = 1;}
530
531         for (int z = zBegin+bb; z < zEnd-bt; z++)
532                 for (int y = 1; y < sy-1; y++)
533                         for (int x = 1; x < sx-1; x++)
534                         {
535                                 const int index = x + y * sx+ z * slabSize;
536                                 // backtrace
537                                 float xBackward = x + dt * velx[index];
538                                 float yBackward = y + dt * vely[index];
539                                 float zBackward = z + dt * velz[index];
540                                 float xTrace    = x - dt * velx[index];
541                                 float yTrace    = y - dt * vely[index];
542                                 float zTrace    = z - dt * velz[index];
543
544                                 // see if it goes outside the boundaries
545                                 bool hasObstacle = 
546                                         (zTrace < 1.0f)    || (zTrace > sz - 2.0f) ||
547                                         (yTrace < 1.0f)    || (yTrace > sy - 2.0f) ||
548                                         (xTrace < 1.0f)    || (xTrace > sx - 2.0f) ||
549                                         (zBackward < 1.0f) || (zBackward > sz - 2.0f) ||
550                                         (yBackward < 1.0f) || (yBackward > sy - 2.0f) ||
551                                         (xBackward < 1.0f) || (xBackward > sx - 2.0f);
552                                 // reuse old advection instead of doing another one...
553                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
554
555                                 // clamp to prevent an out of bounds access when looking into
556                                 // the _obstacles array
557                                 zTrace = (zTrace < 0.5f) ? 0.5f : zTrace;
558                                 zTrace = (zTrace > sz - 1.5f) ? sz - 1.5f : zTrace;
559                                 yTrace = (yTrace < 0.5f) ? 0.5f : yTrace;
560                                 yTrace = (yTrace > sy - 1.5f) ? sy - 1.5f : yTrace;
561                                 xTrace = (xTrace < 0.5f) ? 0.5f : xTrace;
562                                 xTrace = (xTrace > sx - 1.5f) ? sx - 1.5f : xTrace;
563
564                                 // locate neighbors to interpolate,
565                                 // do backward first since we will use the forward indices if a
566                                 // reversion is actually necessary
567                                 zBackward = (zBackward < 0.5f) ? 0.5f : zBackward;
568                                 zBackward = (zBackward > sz - 1.5f) ? sz - 1.5f : zBackward;
569                                 yBackward = (yBackward < 0.5f) ? 0.5f : yBackward;
570                                 yBackward = (yBackward > sy - 1.5f) ? sy - 1.5f : yBackward;
571                                 xBackward = (xBackward < 0.5f) ? 0.5f : xBackward;
572                                 xBackward = (xBackward > sx - 1.5f) ? sx - 1.5f : xBackward;
573
574                                 int x0 = (int)xBackward;
575                                 int x1 = x0 + 1;
576                                 int y0 = (int)yBackward;
577                                 int y1 = y0 + 1;
578                                 int z0 = (int)zBackward;
579                                 int z1 = z0 + 1;
580                                 if(obstacles && !hasObstacle) {
581                                         hasObstacle = hasObstacle || 
582                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
583                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
584                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
585                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
586                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
587                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
588                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
589                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
590                                 }
591                                 // reuse old advection instead of doing another one...
592                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
593
594                                 x0 = (int)xTrace;
595                                 x1 = x0 + 1;
596                                 y0 = (int)yTrace;
597                                 y1 = y0 + 1;
598                                 z0 = (int)zTrace;
599                                 z1 = z0 + 1;
600                                 if(obstacles && !hasObstacle) {
601                                         hasObstacle = hasObstacle || 
602                                                 obstacles[x0 + y0 * sx + z0*slabSize] ||
603                                                 obstacles[x0 + y1 * sx + z0*slabSize] ||
604                                                 obstacles[x1 + y0 * sx + z0*slabSize] ||
605                                                 obstacles[x1 + y1 * sx + z0*slabSize] ||
606                                                 obstacles[x0 + y0 * sx + z1*slabSize] ||
607                                                 obstacles[x0 + y1 * sx + z1*slabSize] ||
608                                                 obstacles[x1 + y0 * sx + z1*slabSize] ||
609                                                 obstacles[x1 + y1 * sx + z1*slabSize] ;
610                                 } // obstacle array
611                                 // reuse old advection instead of doing another one...
612                                 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }
613
614                                 // see if either the forward or backward ray went into
615                                 // a boundary
616                                 if (hasObstacle) {
617                                         // get interpolation weights
618                                         float s1 = xTrace - x0;
619                                         float s0 = 1.0f - s1;
620                                         float t1 = yTrace - y0;
621                                         float t0 = 1.0f - t1;
622                                         float u1 = zTrace - z0;
623                                         float u0 = 1.0f - u1;
624
625                                         const int i000 = x0 + y0 * sx + z0 * slabSize;
626                                         const int i010 = x0 + y1 * sx + z0 * slabSize;
627                                         const int i100 = x1 + y0 * sx + z0 * slabSize;
628                                         const int i110 = x1 + y1 * sx + z0 * slabSize;
629                                         const int i001 = x0 + y0 * sx + z1 * slabSize;
630                                         const int i011 = x0 + y1 * sx + z1 * slabSize;
631                                         const int i101 = x1 + y0 * sx + z1 * slabSize;
632                                         const int i111 = x1 + y1 * sx + z1 * slabSize;
633
634                                         // interpolate, (indices could be computed once)
635                                         newField[index] = u0 * (s0 * (
636                                                                 t0 * oldField[i000] +
637                                                                 t1 * oldField[i010]) +
638                                                         s1 * (t0 * oldField[i100] +
639                                                                 t1 * oldField[i110])) +
640                                                 u1 * (s0 * (t0 * oldField[i001] +
641                                                                         t1 * oldField[i011]) +
642                                                                 s1 * (t0 * oldField[i101] +
643                                                                         t1 * oldField[i111])); 
644                                 }
645                         } // xyz
646 }