|
Blender
V2.59
|
00001 00004 00005 // This file is part of Wavelet Turbulence. 00006 // 00007 // Wavelet Turbulence is free software: you can redistribute it and/or modify 00008 // it under the terms of the GNU General Public License as published by 00009 // the Free Software Foundation, either version 3 of the License, or 00010 // (at your option) any later version. 00011 // 00012 // Wavelet Turbulence is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU General Public License for more details. 00016 // 00017 // You should have received a copy of the GNU General Public License 00018 // along with Wavelet Turbulence. If not, see <http://www.gnu.org/licenses/>. 00019 // 00020 // Copyright 2008 Theodore Kim and Nils Thuerey 00021 // 00022 // FLUID_3D.cpp: implementation of the static functions of the FLUID_3D class. 00023 // 00025 // Heavy parallel optimization done. Many of the old functions now 00026 // take begin and end parameters and process only specified part of the data. 00027 // Some functions were divided into multiple ones. 00028 // - MiikaH 00030 00031 #include <zlib.h> 00032 #include "FLUID_3D.h" 00033 #include "IMAGE.h" 00034 #include "WTURBULENCE.h" 00035 #include "INTERPOLATE.h" 00036 00038 // add a test cube of density to the center 00040 /* 00041 void FLUID_3D::addSmokeColumn() { 00042 addSmokeTestCase(_density, _res, 1.0); 00043 // addSmokeTestCase(_zVelocity, _res, 1.0); 00044 addSmokeTestCase(_heat, _res, 1.0); 00045 if (_wTurbulence) { 00046 addSmokeTestCase(_wTurbulence->getDensityBig(), _wTurbulence->getResBig(), 1.0); 00047 } 00048 } 00049 */ 00050 00052 // generic static version, so that it can be applied to the 00053 // WTURBULENCE grid as well 00055 00056 void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res) 00057 { 00058 const int slabSize = res[0]*res[1]; int maxRes = (int)MAX3V(res); 00059 float dx = 1.0f / (float)maxRes; 00060 00061 float xTotal = dx * res[0]; 00062 float yTotal = dx * res[1]; 00063 00064 float heighMin = 0.05; 00065 float heighMax = 0.10; 00066 00067 for (int y = 0; y < res[2]; y++) 00068 for (int z = (int)(heighMin*res[2]); z <= (int)(heighMax * res[2]); z++) 00069 for (int x = 0; x < res[0]; x++) { 00070 float xLength = x * dx - xTotal * 0.4f; 00071 float yLength = y * dx - yTotal * 0.5f; 00072 float radius = sqrtf(xLength * xLength + yLength * yLength); 00073 00074 if (radius < 0.075f * xTotal) { 00075 int index = x + y * res[0] + z * slabSize; 00076 field[index] = 1.0f; 00077 } 00078 } 00079 } 00080 00081 00083 // set x direction to Neumann boundary conditions 00085 void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd) 00086 { 00087 const int slabSize = res[0] * res[1]; 00088 int index; 00089 for (int z = zBegin; z < zEnd; z++) 00090 for (int y = 0; y < res[1]; y++) 00091 { 00092 // left slab 00093 index = y * res[0] + z * slabSize; 00094 field[index] = field[index + 2]; 00095 00096 // right slab 00097 index += res[0] - 1; 00098 field[index] = field[index - 2]; 00099 } 00100 00101 // fix, force top slab to only allow outwards flux 00102 for (int y = 0; y < res[1]; y++) 00103 for (int z = zBegin; z < zEnd; z++) 00104 { 00105 // top slab 00106 index = y * res[0] + z * slabSize; 00107 index += res[0] - 1; 00108 if(field[index]<0.) field[index] = 0.; 00109 index -= 1; 00110 if(field[index]<0.) field[index] = 0.; 00111 } 00112 } 00113 00115 // set y direction to Neumann boundary conditions 00117 void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd) 00118 { 00119 const int slabSize = res[0] * res[1]; 00120 int index; 00121 for (int z = zBegin; z < zEnd; z++) 00122 for (int x = 0; x < res[0]; x++) 00123 { 00124 // bottom slab 00125 index = x + z * slabSize; 00126 field[index] = field[index + 2 * res[0]]; 00127 00128 // top slab 00129 index += slabSize - res[0]; 00130 field[index] = field[index - 2 * res[0]]; 00131 } 00132 00133 // fix, force top slab to only allow outwards flux 00134 for (int z = zBegin; z < zEnd; z++) 00135 for (int x = 0; x < res[0]; x++) 00136 { 00137 // top slab 00138 index = x + z * slabSize; 00139 index += slabSize - res[0]; 00140 if(field[index]<0.) field[index] = 0.; 00141 index -= res[0]; 00142 if(field[index]<0.) field[index] = 0.; 00143 } 00144 00145 } 00146 00148 // set z direction to Neumann boundary conditions 00150 void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd) 00151 { 00152 const int slabSize = res[0] * res[1]; 00153 const int totalCells = res[0] * res[1] * res[2]; 00154 const int cellsslab = totalCells - slabSize; 00155 int index; 00156 00157 index = 0; 00158 if (zBegin == 0) 00159 for (int y = 0; y < res[1]; y++) 00160 for (int x = 0; x < res[0]; x++, index++) 00161 { 00162 // front slab 00163 field[index] = field[index + 2 * slabSize]; 00164 } 00165 00166 if (zEnd == res[2]) 00167 { 00168 index = 0; 00169 int indexx = 0; 00170 00171 for (int y = 0; y < res[1]; y++) 00172 for (int x = 0; x < res[0]; x++, index++) 00173 { 00174 00175 // back slab 00176 indexx = index + cellsslab; 00177 field[indexx] = field[indexx - 2 * slabSize]; 00178 } 00179 00180 00181 // fix, force top slab to only allow outwards flux 00182 for (int y = 0; y < res[1]; y++) 00183 for (int x = 0; x < res[0]; x++) 00184 { 00185 // top slab 00186 index = x + y * res[0]; 00187 index += cellsslab; 00188 if(field[index]<0.) field[index] = 0.; 00189 index -= slabSize; 00190 if(field[index]<0.) field[index] = 0.; 00191 } 00192 00193 } // zEnd == res[2] 00194 00195 } 00196 00198 // set x direction to zero 00200 void FLUID_3D::setZeroX(float* field, Vec3Int res, int zBegin, int zEnd) 00201 { 00202 const int slabSize = res[0] * res[1]; 00203 int index; 00204 for (int z = zBegin; z < zEnd; z++) 00205 for (int y = 0; y < res[1]; y++) 00206 { 00207 // left slab 00208 index = y * res[0] + z * slabSize; 00209 field[index] = 0.0f; 00210 00211 // right slab 00212 index += res[0] - 1; 00213 field[index] = 0.0f; 00214 } 00215 } 00216 00218 // set y direction to zero 00220 void FLUID_3D::setZeroY(float* field, Vec3Int res, int zBegin, int zEnd) 00221 { 00222 const int slabSize = res[0] * res[1]; 00223 int index; 00224 for (int z = zBegin; z < zEnd; z++) 00225 for (int x = 0; x < res[0]; x++) 00226 { 00227 // bottom slab 00228 index = x + z * slabSize; 00229 field[index] = 0.0f; 00230 00231 // top slab 00232 index += slabSize - res[0]; 00233 field[index] = 0.0f; 00234 } 00235 } 00236 00238 // set z direction to zero 00240 void FLUID_3D::setZeroZ(float* field, Vec3Int res, int zBegin, int zEnd) 00241 { 00242 const int slabSize = res[0] * res[1]; 00243 const int totalCells = res[0] * res[1] * res[2]; 00244 00245 int index = 0; 00246 if ((zBegin == 0)) 00247 for (int y = 0; y < res[1]; y++) 00248 for (int x = 0; x < res[0]; x++, index++) 00249 { 00250 // front slab 00251 field[index] = 0.0f; 00252 } 00253 00254 if (zEnd == res[2]) 00255 { 00256 index=0; 00257 int indexx=0; 00258 const int cellsslab = totalCells - slabSize; 00259 00260 for (int y = 0; y < res[1]; y++) 00261 for (int x = 0; x < res[0]; x++, index++) 00262 { 00263 00264 // back slab 00265 indexx = index + cellsslab; 00266 field[indexx] = 0.0f; 00267 } 00268 } 00269 } 00271 // copy grid boundary 00273 void FLUID_3D::copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd) 00274 { 00275 const int slabSize = res[0] * res[1]; 00276 int index; 00277 for (int z = zBegin; z < zEnd; z++) 00278 for (int y = 0; y < res[1]; y++) 00279 { 00280 // left slab 00281 index = y * res[0] + z * slabSize; 00282 field[index] = field[index + 1]; 00283 00284 // right slab 00285 index += res[0] - 1; 00286 field[index] = field[index - 1]; 00287 } 00288 } 00289 void FLUID_3D::copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd) 00290 { 00291 const int slabSize = res[0] * res[1]; 00292 //const int totalCells = res[0] * res[1] * res[2]; 00293 int index; 00294 for (int z = zBegin; z < zEnd; z++) 00295 for (int x = 0; x < res[0]; x++) 00296 { 00297 // bottom slab 00298 index = x + z * slabSize; 00299 field[index] = field[index + res[0]]; 00300 // top slab 00301 index += slabSize - res[0]; 00302 field[index] = field[index - res[0]]; 00303 } 00304 } 00305 void FLUID_3D::copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd) 00306 { 00307 const int slabSize = res[0] * res[1]; 00308 const int totalCells = res[0] * res[1] * res[2]; 00309 int index=0; 00310 00311 if ((zBegin == 0)) 00312 for (int y = 0; y < res[1]; y++) 00313 for (int x = 0; x < res[0]; x++, index++) 00314 { 00315 field[index] = field[index + slabSize]; 00316 } 00317 00318 if ((zEnd == res[2])) 00319 { 00320 00321 index=0; 00322 int indexx=0; 00323 const int cellsslab = totalCells - slabSize; 00324 00325 for (int y = 0; y < res[1]; y++) 00326 for (int x = 0; x < res[0]; x++, index++) 00327 { 00328 // back slab 00329 indexx = index + cellsslab; 00330 field[indexx] = field[indexx - slabSize]; 00331 } 00332 } 00333 } 00334 00336 // advect field with the semi lagrangian method 00338 void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely, const float* velz, 00339 float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd) 00340 { 00341 const int xres = res[0]; 00342 const int yres = res[1]; 00343 const int zres = res[2]; 00344 const int slabSize = res[0] * res[1]; 00345 00346 00347 for (int z = zBegin; z < zEnd; z++) 00348 for (int y = 0; y < yres; y++) 00349 for (int x = 0; x < xres; x++) 00350 { 00351 const int index = x + y * xres + z * xres*yres; 00352 00353 // backtrace 00354 float xTrace = x - dt * velx[index]; 00355 float yTrace = y - dt * vely[index]; 00356 float zTrace = z - dt * velz[index]; 00357 00358 // clamp backtrace to grid boundaries 00359 if (xTrace < 0.5) xTrace = 0.5; 00360 if (xTrace > xres - 1.5) xTrace = xres - 1.5; 00361 if (yTrace < 0.5) yTrace = 0.5; 00362 if (yTrace > yres - 1.5) yTrace = yres - 1.5; 00363 if (zTrace < 0.5) zTrace = 0.5; 00364 if (zTrace > zres - 1.5) zTrace = zres - 1.5; 00365 00366 // locate neighbors to interpolate 00367 const int x0 = (int)xTrace; 00368 const int x1 = x0 + 1; 00369 const int y0 = (int)yTrace; 00370 const int y1 = y0 + 1; 00371 const int z0 = (int)zTrace; 00372 const int z1 = z0 + 1; 00373 00374 // get interpolation weights 00375 const float s1 = xTrace - x0; 00376 const float s0 = 1.0f - s1; 00377 const float t1 = yTrace - y0; 00378 const float t0 = 1.0f - t1; 00379 const float u1 = zTrace - z0; 00380 const float u0 = 1.0f - u1; 00381 00382 const int i000 = x0 + y0 * xres + z0 * slabSize; 00383 const int i010 = x0 + y1 * xres + z0 * slabSize; 00384 const int i100 = x1 + y0 * xres + z0 * slabSize; 00385 const int i110 = x1 + y1 * xres + z0 * slabSize; 00386 const int i001 = x0 + y0 * xres + z1 * slabSize; 00387 const int i011 = x0 + y1 * xres + z1 * slabSize; 00388 const int i101 = x1 + y0 * xres + z1 * slabSize; 00389 const int i111 = x1 + y1 * xres + z1 * slabSize; 00390 00391 // interpolate 00392 // (indices could be computed once) 00393 newField[index] = u0 * (s0 * (t0 * oldField[i000] + 00394 t1 * oldField[i010]) + 00395 s1 * (t0 * oldField[i100] + 00396 t1 * oldField[i110])) + 00397 u1 * (s0 * (t0 * oldField[i001] + 00398 t1 * oldField[i011]) + 00399 s1 * (t0 * oldField[i101] + 00400 t1 * oldField[i111])); 00401 } 00402 } 00403 00404 00406 // advect field with the maccormack method 00407 // 00408 // comments are the pseudocode from selle's paper 00410 void FLUID_3D::advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 00411 float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd) 00412 { 00413 /*const int sx= res[0]; 00414 const int sy= res[1]; 00415 const int sz= res[2]; 00416 00417 for (int x = 0; x < sx * sy * sz; x++) 00418 phiHatN[x] = phiHatN1[x] = oldField[x];*/ // not needed as all the values are written first 00419 00420 float*& phiN = oldField; 00421 float*& phiN1 = tempResult; 00422 00423 00424 00425 // phiHatN1 = A(phiN) 00426 advectFieldSemiLagrange( dt, xVelocity, yVelocity, zVelocity, phiN, phiN1, res, zBegin, zEnd); // uses wide data from old field and velocities (both are whole) 00427 } 00428 00429 00430 00431 void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 00432 float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd) 00433 { 00434 float* phiHatN = tempResult; 00435 float* t1 = temp1; 00436 const int sx= res[0]; 00437 const int sy= res[1]; 00438 00439 float*& phiN = oldField; 00440 float*& phiN1 = newField; 00441 00442 00443 00444 // phiHatN = A^R(phiHatN1) 00445 advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN, t1, res, zBegin, zEnd); // uses wide data from old field and velocities (both are whole) 00446 00447 // phiN1 = phiHatN1 + (phiN - phiHatN) / 2 00448 const int border = 0; 00449 for (int z = zBegin+border; z < zEnd-border; z++) 00450 for (int y = border; y < sy-border; y++) 00451 for (int x = border; x < sx-border; x++) { 00452 int index = x + y * sx + z * sx*sy; 00453 phiN1[index] = phiHatN[index] + (phiN[index] - t1[index]) * 0.50f; 00454 //phiN1[index] = phiHatN1[index]; // debug, correction off 00455 } 00456 copyBorderX(phiN1, res, zBegin, zEnd); 00457 copyBorderY(phiN1, res, zBegin, zEnd); 00458 copyBorderZ(phiN1, res, zBegin, zEnd); 00459 00460 // clamp any newly created extrema 00461 clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd); // uses wide data from old field and velocities (both are whole) 00462 00463 // if the error estimate was bad, revert to first order 00464 clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd); // phiHatN is only used at cells within thread range, so its ok 00465 00466 } 00467 00468 00470 // Clamp the extrema generated by the BFECC error correction 00472 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely, const float* velz, 00473 float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd) 00474 { 00475 const int xres= res[0]; 00476 const int yres= res[1]; 00477 const int zres= res[2]; 00478 const int slabSize = res[0] * res[1]; 00479 00480 int bb=0; 00481 int bt=0; 00482 00483 if (zBegin == 0) {bb = 1;} 00484 if (zEnd == res[2]) {bt = 1;} 00485 00486 00487 for (int z = zBegin+bb; z < zEnd-bt; z++) 00488 for (int y = 1; y < yres-1; y++) 00489 for (int x = 1; x < xres-1; x++) 00490 { 00491 const int index = x + y * xres+ z * xres*yres; 00492 // backtrace 00493 float xTrace = x - dt * velx[index]; 00494 float yTrace = y - dt * vely[index]; 00495 float zTrace = z - dt * velz[index]; 00496 00497 // clamp backtrace to grid boundaries 00498 if (xTrace < 0.5) xTrace = 0.5; 00499 if (xTrace > xres - 1.5) xTrace = xres - 1.5; 00500 if (yTrace < 0.5) yTrace = 0.5; 00501 if (yTrace > yres - 1.5) yTrace = yres - 1.5; 00502 if (zTrace < 0.5) zTrace = 0.5; 00503 if (zTrace > zres - 1.5) zTrace = zres - 1.5; 00504 00505 // locate neighbors to interpolate 00506 const int x0 = (int)xTrace; 00507 const int x1 = x0 + 1; 00508 const int y0 = (int)yTrace; 00509 const int y1 = y0 + 1; 00510 const int z0 = (int)zTrace; 00511 const int z1 = z0 + 1; 00512 00513 const int i000 = x0 + y0 * xres + z0 * slabSize; 00514 const int i010 = x0 + y1 * xres + z0 * slabSize; 00515 const int i100 = x1 + y0 * xres + z0 * slabSize; 00516 const int i110 = x1 + y1 * xres + z0 * slabSize; 00517 const int i001 = x0 + y0 * xres + z1 * slabSize; 00518 const int i011 = x0 + y1 * xres + z1 * slabSize; 00519 const int i101 = x1 + y0 * xres + z1 * slabSize; 00520 const int i111 = x1 + y1 * xres + z1 * slabSize; 00521 00522 float minField = oldField[i000]; 00523 float maxField = oldField[i000]; 00524 00525 minField = (oldField[i010] < minField) ? oldField[i010] : minField; 00526 maxField = (oldField[i010] > maxField) ? oldField[i010] : maxField; 00527 00528 minField = (oldField[i100] < minField) ? oldField[i100] : minField; 00529 maxField = (oldField[i100] > maxField) ? oldField[i100] : maxField; 00530 00531 minField = (oldField[i110] < minField) ? oldField[i110] : minField; 00532 maxField = (oldField[i110] > maxField) ? oldField[i110] : maxField; 00533 00534 minField = (oldField[i001] < minField) ? oldField[i001] : minField; 00535 maxField = (oldField[i001] > maxField) ? oldField[i001] : maxField; 00536 00537 minField = (oldField[i011] < minField) ? oldField[i011] : minField; 00538 maxField = (oldField[i011] > maxField) ? oldField[i011] : maxField; 00539 00540 minField = (oldField[i101] < minField) ? oldField[i101] : minField; 00541 maxField = (oldField[i101] > maxField) ? oldField[i101] : maxField; 00542 00543 minField = (oldField[i111] < minField) ? oldField[i111] : minField; 00544 maxField = (oldField[i111] > maxField) ? oldField[i111] : maxField; 00545 00546 newField[index] = (newField[index] > maxField) ? maxField : newField[index]; 00547 newField[index] = (newField[index] < minField) ? minField : newField[index]; 00548 } 00549 } 00550 00552 // Reverts any backtraces that go into boundaries back to first 00553 // order -- in this case the error correction term was totally 00554 // incorrect 00556 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely, const float* velz, 00557 float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd) 00558 { 00559 const int sx= res[0]; 00560 const int sy= res[1]; 00561 const int sz= res[2]; 00562 const int slabSize = res[0] * res[1]; 00563 00564 int bb=0; 00565 int bt=0; 00566 00567 if (zBegin == 0) {bb = 1;} 00568 if (zEnd == res[2]) {bt = 1;} 00569 00570 for (int z = zBegin+bb; z < zEnd-bt; z++) 00571 for (int y = 1; y < sy-1; y++) 00572 for (int x = 1; x < sx-1; x++) 00573 { 00574 const int index = x + y * sx+ z * slabSize; 00575 // backtrace 00576 float xBackward = x + dt * velx[index]; 00577 float yBackward = y + dt * vely[index]; 00578 float zBackward = z + dt * velz[index]; 00579 float xTrace = x - dt * velx[index]; 00580 float yTrace = y - dt * vely[index]; 00581 float zTrace = z - dt * velz[index]; 00582 00583 // see if it goes outside the boundaries 00584 bool hasObstacle = 00585 (zTrace < 1.0f) || (zTrace > sz - 2.0f) || 00586 (yTrace < 1.0f) || (yTrace > sy - 2.0f) || 00587 (xTrace < 1.0f) || (xTrace > sx - 2.0f) || 00588 (zBackward < 1.0f) || (zBackward > sz - 2.0f) || 00589 (yBackward < 1.0f) || (yBackward > sy - 2.0f) || 00590 (xBackward < 1.0f) || (xBackward > sx - 2.0f); 00591 // reuse old advection instead of doing another one... 00592 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; } 00593 00594 // clamp to prevent an out of bounds access when looking into 00595 // the _obstacles array 00596 zTrace = (zTrace < 0.5f) ? 0.5f : zTrace; 00597 zTrace = (zTrace > sz - 1.5f) ? sz - 1.5f : zTrace; 00598 yTrace = (yTrace < 0.5f) ? 0.5f : yTrace; 00599 yTrace = (yTrace > sy - 1.5f) ? sy - 1.5f : yTrace; 00600 xTrace = (xTrace < 0.5f) ? 0.5f : xTrace; 00601 xTrace = (xTrace > sx - 1.5f) ? sx - 1.5f : xTrace; 00602 00603 // locate neighbors to interpolate, 00604 // do backward first since we will use the forward indices if a 00605 // reversion is actually necessary 00606 zBackward = (zBackward < 0.5f) ? 0.5f : zBackward; 00607 zBackward = (zBackward > sz - 1.5f) ? sz - 1.5f : zBackward; 00608 yBackward = (yBackward < 0.5f) ? 0.5f : yBackward; 00609 yBackward = (yBackward > sy - 1.5f) ? sy - 1.5f : yBackward; 00610 xBackward = (xBackward < 0.5f) ? 0.5f : xBackward; 00611 xBackward = (xBackward > sx - 1.5f) ? sx - 1.5f : xBackward; 00612 00613 int x0 = (int)xBackward; 00614 int x1 = x0 + 1; 00615 int y0 = (int)yBackward; 00616 int y1 = y0 + 1; 00617 int z0 = (int)zBackward; 00618 int z1 = z0 + 1; 00619 if(obstacles && !hasObstacle) { 00620 hasObstacle = hasObstacle || 00621 obstacles[x0 + y0 * sx + z0*slabSize] || 00622 obstacles[x0 + y1 * sx + z0*slabSize] || 00623 obstacles[x1 + y0 * sx + z0*slabSize] || 00624 obstacles[x1 + y1 * sx + z0*slabSize] || 00625 obstacles[x0 + y0 * sx + z1*slabSize] || 00626 obstacles[x0 + y1 * sx + z1*slabSize] || 00627 obstacles[x1 + y0 * sx + z1*slabSize] || 00628 obstacles[x1 + y1 * sx + z1*slabSize] ; 00629 } 00630 // reuse old advection instead of doing another one... 00631 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; } 00632 00633 x0 = (int)xTrace; 00634 x1 = x0 + 1; 00635 y0 = (int)yTrace; 00636 y1 = y0 + 1; 00637 z0 = (int)zTrace; 00638 z1 = z0 + 1; 00639 if(obstacles && !hasObstacle) { 00640 hasObstacle = hasObstacle || 00641 obstacles[x0 + y0 * sx + z0*slabSize] || 00642 obstacles[x0 + y1 * sx + z0*slabSize] || 00643 obstacles[x1 + y0 * sx + z0*slabSize] || 00644 obstacles[x1 + y1 * sx + z0*slabSize] || 00645 obstacles[x0 + y0 * sx + z1*slabSize] || 00646 obstacles[x0 + y1 * sx + z1*slabSize] || 00647 obstacles[x1 + y0 * sx + z1*slabSize] || 00648 obstacles[x1 + y1 * sx + z1*slabSize] ; 00649 } // obstacle array 00650 // reuse old advection instead of doing another one... 00651 if(hasObstacle) { newField[index] = oldAdvection[index]; continue; } 00652 00653 // see if either the forward or backward ray went into 00654 // a boundary 00655 if (hasObstacle) { 00656 // get interpolation weights 00657 float s1 = xTrace - x0; 00658 float s0 = 1.0f - s1; 00659 float t1 = yTrace - y0; 00660 float t0 = 1.0f - t1; 00661 float u1 = zTrace - z0; 00662 float u0 = 1.0f - u1; 00663 00664 const int i000 = x0 + y0 * sx + z0 * slabSize; 00665 const int i010 = x0 + y1 * sx + z0 * slabSize; 00666 const int i100 = x1 + y0 * sx + z0 * slabSize; 00667 const int i110 = x1 + y1 * sx + z0 * slabSize; 00668 const int i001 = x0 + y0 * sx + z1 * slabSize; 00669 const int i011 = x0 + y1 * sx + z1 * slabSize; 00670 const int i101 = x1 + y0 * sx + z1 * slabSize; 00671 const int i111 = x1 + y1 * sx + z1 * slabSize; 00672 00673 // interpolate, (indices could be computed once) 00674 newField[index] = u0 * (s0 * ( 00675 t0 * oldField[i000] + 00676 t1 * oldField[i010]) + 00677 s1 * (t0 * oldField[i100] + 00678 t1 * oldField[i110])) + 00679 u1 * (s0 * (t0 * oldField[i001] + 00680 t1 * oldField[i011]) + 00681 s1 * (t0 * oldField[i101] + 00682 t1 * oldField[i111])); 00683 } 00684 } // xyz 00685 }