← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty

« back to all changes in this revision

Viewing changes to intern/smoke/intern/FLUID_3D.cpp

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
intern/smoke/intern/FLUID_3D.cpp
34
34
#include "SPHERE.h"
35
35
#include <zlib.h>
36
36
 
 
37
#include "float.h"
 
38
 
37
39
#if PARALLEL==1
38
40
#include <omp.h>
39
41
#endif // PARALLEL 
42
44
// Construction/Destruction
43
45
//////////////////////////////////////////////////////////////////////
44
46
 
45
 
FLUID_3D::FLUID_3D(int *res, float *p0) :
 
47
FLUID_3D::FLUID_3D(int *res, float dx, float dtdef, int init_heat, int init_fire, int init_colors) :
46
48
        _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f)
47
49
{
48
50
        // set simulation consts
49
 
        _dt = DT_DEFAULT;       // just in case. set in step from a RNA factor
50
 
        
51
 
        // start point of array
52
 
        _p0[0] = p0[0];
53
 
        _p0[1] = p0[1];
54
 
        _p0[2] = p0[2];
 
51
        _dt = dtdef;    // just in case. set in step from a RNA factor
55
52
 
56
53
        _iterations = 100;
57
54
        _tempAmb = 0; 
70
67
        */
71
68
        
72
69
        // scale the constants according to the refinement of the grid
73
 
        _dx = 1.0f / (float)_maxRes;
 
70
        if (!dx)
 
71
                _dx = 1.0f / (float)_maxRes;
 
72
        else
 
73
                _dx = dx;
74
74
        _constantScaling = 64.0f / _maxRes;
75
75
        _constantScaling = (_constantScaling < 1.0f) ? 1.0f : _constantScaling;
76
76
        _vorticityEps = 2.0f / _constantScaling; // Just in case set a default value
81
81
        _xVelocity    = new float[_totalCells];
82
82
        _yVelocity    = new float[_totalCells];
83
83
        _zVelocity    = new float[_totalCells];
 
84
        _xVelocityOb  = new float[_totalCells];
 
85
        _yVelocityOb  = new float[_totalCells];
 
86
        _zVelocityOb  = new float[_totalCells];
84
87
        _xVelocityOld = new float[_totalCells];
85
88
        _yVelocityOld = new float[_totalCells];
86
89
        _zVelocityOld = new float[_totalCells];
89
92
        _zForce       = new float[_totalCells];
90
93
        _density      = new float[_totalCells];
91
94
        _densityOld   = new float[_totalCells];
92
 
        _heat         = new float[_totalCells];
93
 
        _heatOld      = new float[_totalCells];
94
95
        _obstacles    = new unsigned char[_totalCells]; // set 0 at end of step
95
96
 
96
97
        // For threaded version:
98
99
        _yVelocityTemp = new float[_totalCells];
99
100
        _zVelocityTemp = new float[_totalCells];
100
101
        _densityTemp   = new float[_totalCells];
101
 
        _heatTemp      = new float[_totalCells];
102
102
 
103
103
        // DG TODO: check if alloc went fine
104
104
 
106
106
        {
107
107
                _density[x]      = 0.0f;
108
108
                _densityOld[x]   = 0.0f;
109
 
                _heat[x]         = 0.0f;
110
 
                _heatOld[x]      = 0.0f;
111
109
                _xVelocity[x]    = 0.0f;
112
110
                _yVelocity[x]    = 0.0f;
113
111
                _zVelocity[x]    = 0.0f;
 
112
                _xVelocityOb[x]  = 0.0f;
 
113
                _yVelocityOb[x]  = 0.0f;
 
114
                _zVelocityOb[x]  = 0.0f;
114
115
                _xVelocityOld[x] = 0.0f;
115
116
                _yVelocityOld[x] = 0.0f;
116
117
                _zVelocityOld[x] = 0.0f;
120
121
                _obstacles[x]    = false;
121
122
        }
122
123
 
 
124
        /* heat */
 
125
        _heat = _heatOld = _heatTemp = NULL;
 
126
        if (init_heat) {
 
127
                initHeat();
 
128
        }
 
129
        // Fire simulation
 
130
        _flame = _fuel = _fuelTemp = _fuelOld = NULL;
 
131
        _react = _reactTemp = _reactOld = NULL;
 
132
        if (init_fire) {
 
133
                initFire();
 
134
        }
 
135
        // Smoke color
 
136
        _color_r = _color_rOld = _color_rTemp = NULL;
 
137
        _color_g = _color_gOld = _color_gTemp = NULL;
 
138
        _color_b = _color_bOld = _color_bTemp = NULL;
 
139
        if (init_colors) {
 
140
                initColors(0.0f, 0.0f, 0.0f);
 
141
        }
 
142
 
123
143
        // boundary conditions of the fluid domain
124
144
        // set default values -> vertically non-colliding
125
145
        _domainBcFront = true;
130
150
        _domainBcRight  = _domainBcLeft;
131
151
 
132
152
        _colloPrev = 1; // default value
133
 
 
134
 
 
 
153
}
 
154
 
 
155
void FLUID_3D::initHeat()
 
156
{
 
157
        if (!_heat) {
 
158
                _heat         = new float[_totalCells];
 
159
                _heatOld      = new float[_totalCells];
 
160
                _heatTemp      = new float[_totalCells];
 
161
 
 
162
                for (int x = 0; x < _totalCells; x++)
 
163
                {
 
164
                        _heat[x]         = 0.0f;
 
165
                        _heatOld[x]      = 0.0f;
 
166
                }
 
167
        }
 
168
}
 
169
 
 
170
void FLUID_3D::initFire()
 
171
{
 
172
        if (!_flame) {
 
173
                _flame          = new float[_totalCells];
 
174
                _fuel           = new float[_totalCells];
 
175
                _fuelTemp       = new float[_totalCells];
 
176
                _fuelOld        = new float[_totalCells];
 
177
                _react          = new float[_totalCells];
 
178
                _reactTemp      = new float[_totalCells];
 
179
                _reactOld       = new float[_totalCells];
 
180
 
 
181
                for (int x = 0; x < _totalCells; x++)
 
182
                {
 
183
                        _flame[x]               = 0.0f;
 
184
                        _fuel[x]                = 0.0f;
 
185
                        _fuelTemp[x]    = 0.0f;
 
186
                        _fuelOld[x]             = 0.0f;
 
187
                        _react[x]               = 0.0f;
 
188
                        _reactTemp[x]   = 0.0f;
 
189
                        _reactOld[x]    = 0.0f;
 
190
                }
 
191
        }
 
192
}
 
193
 
 
194
void FLUID_3D::initColors(float init_r, float init_g, float init_b)
 
195
{
 
196
        if (!_color_r) {
 
197
                _color_r                = new float[_totalCells];
 
198
                _color_rOld             = new float[_totalCells];
 
199
                _color_rTemp    = new float[_totalCells];
 
200
                _color_g                = new float[_totalCells];
 
201
                _color_gOld             = new float[_totalCells];
 
202
                _color_gTemp    = new float[_totalCells];
 
203
                _color_b                = new float[_totalCells];
 
204
                _color_bOld             = new float[_totalCells];
 
205
                _color_bTemp    = new float[_totalCells];
 
206
 
 
207
                for (int x = 0; x < _totalCells; x++)
 
208
                {
 
209
                        _color_r[x]             = _density[x] * init_r;
 
210
                        _color_rOld[x]  = 0.0f;
 
211
                        _color_g[x]             = _density[x] * init_g;
 
212
                        _color_gOld[x]  = 0.0f;
 
213
                        _color_b[x]             = _density[x] * init_b;
 
214
                        _color_bOld[x]  = 0.0f;
 
215
                }
 
216
        }
 
217
}
 
218
 
 
219
void FLUID_3D::setBorderObstacles()
 
220
{
 
221
        
135
222
        // set side obstacles
136
 
        int index;
 
223
        unsigned int index;
137
224
        for (int y = 0; y < _yRes; y++)
138
225
        for (int x = 0; x < _xRes; x++)
139
226
        {
140
227
                // bottom slab
141
228
                index = x + y * _xRes;
142
 
                if(_domainBcBottom==1) _obstacles[index] = 1;
 
229
                if(_domainBcBottom) _obstacles[index] = 1;
143
230
 
144
231
                // top slab
145
232
                index += _totalCells - _slabSize;
146
 
                if(_domainBcTop==1) _obstacles[index] = 1;
 
233
                if(_domainBcTop) _obstacles[index] = 1;
147
234
        }
148
235
 
149
236
        for (int z = 0; z < _zRes; z++)
151
238
        {
152
239
                // front slab
153
240
                index = x + z * _slabSize;
154
 
                if(_domainBcFront==1) _obstacles[index] = 1;
 
241
                if(_domainBcFront) _obstacles[index] = 1;
155
242
 
156
243
                // back slab
157
244
                index += _slabSize - _xRes;
158
 
                if(_domainBcBack==1) _obstacles[index] = 1;
 
245
                if(_domainBcBack) _obstacles[index] = 1;
159
246
        }
160
247
 
161
248
        for (int z = 0; z < _zRes; z++)
163
250
        {
164
251
                // left slab
165
252
                index = y * _xRes + z * _slabSize;
166
 
                if(_domainBcLeft==1) _obstacles[index] = 1;
 
253
                if(_domainBcLeft) _obstacles[index] = 1;
167
254
 
168
255
                // right slab
169
256
                index += _xRes - 1;
170
 
                if(_domainBcRight==1) _obstacles[index] = 1;
 
257
                if(_domainBcRight) _obstacles[index] = 1;
171
258
        }
172
 
 
173
259
}
174
260
 
175
261
FLUID_3D::~FLUID_3D()
177
263
        if (_xVelocity) delete[] _xVelocity;
178
264
        if (_yVelocity) delete[] _yVelocity;
179
265
        if (_zVelocity) delete[] _zVelocity;
 
266
        if (_xVelocityOb) delete[] _xVelocityOb;
 
267
        if (_yVelocityOb) delete[] _yVelocityOb;
 
268
        if (_zVelocityOb) delete[] _zVelocityOb;
180
269
        if (_xVelocityOld) delete[] _xVelocityOld;
181
270
        if (_yVelocityOld) delete[] _yVelocityOld;
182
271
        if (_zVelocityOld) delete[] _zVelocityOld;
188
277
        if (_heat) delete[] _heat;
189
278
        if (_heatOld) delete[] _heatOld;
190
279
        if (_obstacles) delete[] _obstacles;
191
 
    // if (_wTurbulence) delete _wTurbulence;
192
280
 
193
281
        if (_xVelocityTemp) delete[] _xVelocityTemp;
194
282
        if (_yVelocityTemp) delete[] _yVelocityTemp;
196
284
        if (_densityTemp) delete[] _densityTemp;
197
285
        if (_heatTemp) delete[] _heatTemp;
198
286
 
 
287
        if (_flame) delete[] _flame;
 
288
        if (_fuel) delete[] _fuel;
 
289
        if (_fuelTemp) delete[] _fuelTemp;
 
290
        if (_fuelOld) delete[] _fuelOld;
 
291
        if (_react) delete[] _react;
 
292
        if (_reactTemp) delete[] _reactTemp;
 
293
        if (_reactOld) delete[] _reactOld;
 
294
 
 
295
        if (_color_r) delete[] _color_r;
 
296
        if (_color_rOld) delete[] _color_rOld;
 
297
        if (_color_rTemp) delete[] _color_rTemp;
 
298
        if (_color_g) delete[] _color_g;
 
299
        if (_color_gOld) delete[] _color_gOld;
 
300
        if (_color_gTemp) delete[] _color_gTemp;
 
301
        if (_color_b) delete[] _color_b;
 
302
        if (_color_bOld) delete[] _color_bOld;
 
303
        if (_color_bTemp) delete[] _color_bTemp;
 
304
 
199
305
    // printf("deleted fluid\n");
200
306
}
201
307
 
202
308
// init direct access functions from blender
203
 
void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision)
 
309
void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision, float *burning_rate,
 
310
                                                          float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp)
204
311
{
205
312
        _alpha = alpha;
206
313
        _beta = beta;
207
314
        _dtFactor = dt_factor;
208
315
        _vorticityRNA = vorticity;
209
316
        _borderColli = borderCollision;
 
317
        _burning_rate = burning_rate;
 
318
        _flame_smoke = flame_smoke;
 
319
        _flame_smoke_color = flame_smoke_color;
 
320
        _flame_vorticity = flame_vorticity;
 
321
        _ignition_temp = flame_ignition_temp;
 
322
        _max_temp = flame_max_temp;
210
323
}
211
324
 
212
325
//////////////////////////////////////////////////////////////////////
213
326
// step simulation once
214
327
//////////////////////////////////////////////////////////////////////
215
 
void FLUID_3D::step(float dt)
 
328
void FLUID_3D::step(float dt, float gravity[3])
216
329
{
 
330
#if 0
217
331
        // If border rules have been changed
218
332
        if (_colloPrev != *_borderColli) {
 
333
                printf("Border collisions changed\n");
 
334
                
 
335
                // DG TODO: Need to check that no animated obstacle flags are overwritten
219
336
                setBorderCollisions();
220
337
        }
 
338
#endif
 
339
 
 
340
        // DG: TODO for the moment redo border for every timestep since it's been deleted every time by moving obstacles
 
341
        setBorderCollisions();
221
342
 
222
343
 
223
344
        // set delta time by dt_factor
225
346
        // set vorticity from RNA value
226
347
        _vorticityEps = (*_vorticityRNA)/_constantScaling;
227
348
 
228
 
 
229
349
#if PARALLEL==1
230
350
        int threadval = 1;
231
351
        threadval = omp_get_max_threads();
258
378
 
259
379
                wipeBoundariesSL(zBegin, zEnd);
260
380
                addVorticity(zBegin, zEnd);
261
 
                addBuoyancy(_heat, _density, zBegin, zEnd);
 
381
                addBuoyancy(_heat, _density, gravity, zBegin, zEnd);
262
382
                addForce(zBegin, zEnd);
263
383
 
264
384
#if PARALLEL==1
289
409
                if (i==0)
290
410
                {
291
411
#endif
292
 
                project();
 
412
                        project();
293
413
#if PARALLEL==1
294
414
                }
295
 
                else
 
415
                else if (i==1)
296
416
                {
297
417
#endif
298
 
                diffuseHeat();
 
418
                        if (_heat) {
 
419
                                diffuseHeat();
 
420
                        }
299
421
#if PARALLEL==1
300
422
                }
301
423
        }
315
437
        SWAP_POINTERS(_density, _densityOld);
316
438
        SWAP_POINTERS(_heat, _heatOld);
317
439
 
 
440
        SWAP_POINTERS(_fuel, _fuelOld);
 
441
        SWAP_POINTERS(_react, _reactOld);
 
442
 
 
443
        SWAP_POINTERS(_color_r, _color_rOld);
 
444
        SWAP_POINTERS(_color_g, _color_gOld);
 
445
        SWAP_POINTERS(_color_b, _color_bOld);
 
446
 
318
447
        advectMacCormackBegin(0, _zRes);
319
448
 
320
449
#if PARALLEL==1
375
504
        SWAP_POINTERS(_yVelocity, _yForce);
376
505
        SWAP_POINTERS(_zVelocity, _zForce);
377
506
 
378
 
 
379
 
 
380
 
 
381
507
        _totalTime += _dt;
382
 
        _totalSteps++;          
 
508
        _totalSteps++;
383
509
 
384
510
        for (int i = 0; i < _totalCells; i++)
385
511
        {
426
552
 
427
553
 
428
554
        // set side obstacles
429
 
        int index;
430
 
        for (int y = 0; y < _yRes; y++)
431
 
        for (int x = 0; x < _xRes; x++)
432
 
        {
433
 
                // front slab
434
 
                index = x + y * _xRes;
435
 
                if(_domainBcBottom==1) _obstacles[index] = 1;
436
 
                else _obstacles[index] = 0;
437
 
 
438
 
                // back slab
439
 
                index += _totalCells - _slabSize;
440
 
                if(_domainBcTop==1) _obstacles[index] = 1;
441
 
                else _obstacles[index] = 0;
442
 
        }
443
 
 
444
 
        for (int z = 0; z < _zRes; z++)
445
 
        for (int x = 0; x < _xRes; x++)
446
 
        {
447
 
                // bottom slab
448
 
                index = x + z * _slabSize;
449
 
                if(_domainBcFront==1) _obstacles[index] = 1;
450
 
                else _obstacles[index] = 0;
451
 
 
452
 
                // top slab
453
 
                index += _slabSize - _xRes;
454
 
                if(_domainBcBack==1) _obstacles[index] = 1;
455
 
                else _obstacles[index] = 0;
456
 
        }
457
 
 
458
 
        for (int z = 0; z < _zRes; z++)
459
 
        for (int y = 0; y < _yRes; y++)
460
 
        {
461
 
                // left slab
462
 
                index = y * _xRes + z * _slabSize;
463
 
                if(_domainBcLeft==1) _obstacles[index] = 1;
464
 
                else _obstacles[index] = 0;
465
 
 
466
 
                // right slab
467
 
                index += _xRes - 1;
468
 
                if(_domainBcRight==1) _obstacles[index] = 1;
469
 
                else _obstacles[index] = 0;
470
 
        }
 
555
        setBorderObstacles();
471
556
}
472
557
 
473
558
//////////////////////////////////////////////////////////////////////
661
746
        setZeroBorder(_yVelocity, _res, zBegin, zEnd);
662
747
        setZeroBorder(_zVelocity, _res, zBegin, zEnd);
663
748
        setZeroBorder(_density, _res, zBegin, zEnd);
 
749
        if (_fuel) {
 
750
                setZeroBorder(_fuel, _res, zBegin, zEnd);
 
751
                setZeroBorder(_react, _res, zBegin, zEnd);
 
752
        }
 
753
        if (_color_r) {
 
754
                setZeroBorder(_color_r, _res, zBegin, zEnd);
 
755
                setZeroBorder(_color_g, _res, zBegin, zEnd);
 
756
                setZeroBorder(_color_b, _res, zBegin, zEnd);
 
757
        }
664
758
}
665
759
 
666
760
void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
686
780
                        _yVelocity[index] = 0.0f;
687
781
                        _zVelocity[index] = 0.0f;
688
782
                        _density[index] = 0.0f;
 
783
                        if (_fuel) {
 
784
                                _fuel[index] = 0.0f;
 
785
                                _react[index] = 0.0f;
 
786
                        }
 
787
                        if (_color_r) {
 
788
                                _color_r[index] = 0.0f;
 
789
                                _color_g[index] = 0.0f;
 
790
                                _color_b[index] = 0.0f;
 
791
                        }
689
792
 
690
793
                        // right slab
691
794
                        index += _xRes - 1;
693
796
                        _yVelocity[index] = 0.0f;
694
797
                        _zVelocity[index] = 0.0f;
695
798
                        _density[index] = 0.0f;
 
799
                        if (_fuel) {
 
800
                                _fuel[index] = 0.0f;
 
801
                                _react[index] = 0.0f;
 
802
                        }
 
803
                        if (_color_r) {
 
804
                                _color_r[index] = 0.0f;
 
805
                                _color_g[index] = 0.0f;
 
806
                                _color_b[index] = 0.0f;
 
807
                        }
696
808
                }
697
809
 
698
810
        /////////////////////////////////////
708
820
                        _yVelocity[index] = 0.0f;
709
821
                        _zVelocity[index] = 0.0f;
710
822
                        _density[index] = 0.0f;
 
823
                        if (_fuel) {
 
824
                                _fuel[index] = 0.0f;
 
825
                                _react[index] = 0.0f;
 
826
                        }
 
827
                        if (_color_r) {
 
828
                                _color_r[index] = 0.0f;
 
829
                                _color_g[index] = 0.0f;
 
830
                                _color_b[index] = 0.0f;
 
831
                        }
711
832
 
712
833
                        // top slab
713
834
                        index += slabSize - _xRes;
715
836
                        _yVelocity[index] = 0.0f;
716
837
                        _zVelocity[index] = 0.0f;
717
838
                        _density[index] = 0.0f;
 
839
                        if (_fuel) {
 
840
                                _fuel[index] = 0.0f;
 
841
                                _react[index] = 0.0f;
 
842
                        }
 
843
                        if (_color_r) {
 
844
                                _color_r[index] = 0.0f;
 
845
                                _color_g[index] = 0.0f;
 
846
                                _color_b[index] = 0.0f;
 
847
                        }
718
848
 
719
849
                }
720
850
 
735
865
                        _yVelocity[index] = 0.0f;
736
866
                        _zVelocity[index] = 0.0f;
737
867
                        _density[index] = 0.0f;
 
868
                        if (_fuel) {
 
869
                                _fuel[index] = 0.0f;
 
870
                                _react[index] = 0.0f;
 
871
                        }
 
872
                        if (_color_r) {
 
873
                                _color_r[index] = 0.0f;
 
874
                                _color_g[index] = 0.0f;
 
875
                                _color_b[index] = 0.0f;
 
876
                        }
738
877
    }
739
878
 
740
879
        if (zEnd == _zRes)
753
892
                                _yVelocity[indexx] = 0.0f;
754
893
                                _zVelocity[indexx] = 0.0f;
755
894
                                _density[indexx] = 0.0f;
 
895
                                if (_fuel) {
 
896
                                        _fuel[index] = 0.0f;
 
897
                                        _react[index] = 0.0f;
 
898
                                }
 
899
                                if (_color_r) {
 
900
                                        _color_r[index] = 0.0f;
 
901
                                        _color_g[index] = 0.0f;
 
902
                                        _color_b[index] = 0.0f;
 
903
                                }
756
904
                        }
757
905
        }
758
906
 
786
934
        memset(_pressure, 0, sizeof(float)*_totalCells);
787
935
        memset(_divergence, 0, sizeof(float)*_totalCells);
788
936
        
 
937
        // set velocity and pressure inside of obstacles to zero
789
938
        setObstacleBoundaries(_pressure, 0, _zRes);
790
939
        
791
940
        // copy out the boundaries
792
 
        if(_domainBcLeft == 0)  setNeumannX(_xVelocity, _res, 0, _zRes);
 
941
        if(!_domainBcLeft)  setNeumannX(_xVelocity, _res, 0, _zRes);
793
942
        else setZeroX(_xVelocity, _res, 0, _zRes); 
794
943
 
795
 
        if(_domainBcFront == 0)   setNeumannY(_yVelocity, _res, 0, _zRes);
 
944
        if(!_domainBcFront)   setNeumannY(_yVelocity, _res, 0, _zRes);
796
945
        else setZeroY(_yVelocity, _res, 0, _zRes); 
797
946
 
798
 
        if(_domainBcTop == 0) setNeumannZ(_zVelocity, _res, 0, _zRes);
 
947
        if(!_domainBcTop) setNeumannZ(_zVelocity, _res, 0, _zRes);
799
948
        else setZeroZ(_zVelocity, _res, 0, _zRes);
800
949
 
801
950
        // calculate divergence
804
953
                for (y = 1; y < _yRes - 1; y++, index += 2)
805
954
                        for (x = 1; x < _xRes - 1; x++, index++)
806
955
                        {
 
956
                                
 
957
                                if(_obstacles[index])
 
958
                                {
 
959
                                        _divergence[index] = 0.0f;
 
960
                                        continue;
 
961
                                }
 
962
                                
 
963
 
807
964
                                float xright = _xVelocity[index + 1];
808
965
                                float xleft  = _xVelocity[index - 1];
809
966
                                float yup    = _yVelocity[index + _xRes];
811
968
                                float ztop   = _zVelocity[index + _slabSize];
812
969
                                float zbottom = _zVelocity[index - _slabSize];
813
970
 
814
 
                                if(_obstacles[index+1]) xright = - _xVelocity[index];
 
971
                                if(_obstacles[index+1]) xright = - _xVelocity[index]; // DG: +=
815
972
                                if(_obstacles[index-1]) xleft  = - _xVelocity[index];
816
973
                                if(_obstacles[index+_xRes]) yup    = - _yVelocity[index];
817
974
                                if(_obstacles[index-_xRes]) ydown  = - _yVelocity[index];
818
975
                                if(_obstacles[index+_slabSize]) ztop    = - _zVelocity[index];
819
976
                                if(_obstacles[index-_slabSize]) zbottom = - _zVelocity[index];
820
977
 
 
978
                                if(_obstacles[index+1] & 8)                     xright  += _xVelocityOb[index + 1];
 
979
                                if(_obstacles[index-1] & 8)                     xleft   += _xVelocityOb[index - 1];
 
980
                                if(_obstacles[index+_xRes] & 8)         yup             += _yVelocityOb[index + _xRes];
 
981
                                if(_obstacles[index-_xRes] & 8)         ydown   += _yVelocityOb[index - _xRes];
 
982
                                if(_obstacles[index+_slabSize] & 8) ztop    += _zVelocityOb[index + _slabSize];
 
983
                                if(_obstacles[index-_slabSize] & 8) zbottom += _zVelocityOb[index - _slabSize];
 
984
 
821
985
                                _divergence[index] = -_dx * 0.5f * (
822
986
                                                xright - xleft +
823
987
                                                yup - ydown +
824
988
                                                ztop - zbottom );
825
989
 
826
 
                                // DG: commenting this helps CG to get a better start, 10-20% speed improvement
827
 
                                // _pressure[index] = 0.0f;
 
990
                                // Pressure is zero anyway since now a local array is used
 
991
                                _pressure[index] = 0.0f;
828
992
                        }
 
993
 
 
994
 
 
995
        /*
 
996
        {
 
997
                float maxvalue = 0;
 
998
                for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
 
999
                {
 
1000
                        if(_divergence[i] > maxvalue)
 
1001
                                maxvalue = _divergence[i];
 
1002
 
 
1003
                }
 
1004
                printf("Max divergence: %f\n", maxvalue);
 
1005
        }
 
1006
        */
 
1007
 
829
1008
        copyBorderAll(_pressure, 0, _zRes);
830
1009
 
 
1010
        /*
 
1011
        {
 
1012
                float maxvalue = 0;
 
1013
                for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
 
1014
                {
 
1015
                        if(_pressure[i] > maxvalue)
 
1016
                                maxvalue = _pressure[i];
 
1017
                }
 
1018
                printf("Max pressure BEFORE: %f\n", maxvalue);
 
1019
        }
 
1020
        */
 
1021
 
831
1022
        // solve Poisson equation
832
1023
        solvePressurePre(_pressure, _divergence, _obstacles);
833
1024
 
 
1025
        {
 
1026
                float maxvalue = 0;
 
1027
                for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
 
1028
                {
 
1029
                        if(_pressure[i] > maxvalue)
 
1030
                                maxvalue = _pressure[i];
 
1031
 
 
1032
                        /* HACK: Animated collision object sometimes result in a non converging solvePressurePre() */ 
 
1033
                        if(_pressure[i] > _dx * _dt)
 
1034
                                _pressure[i] = _dx * _dt;
 
1035
                        else if(_pressure[i] < -_dx * _dt)
 
1036
                                _pressure[i] = -_dx * _dt;
 
1037
 
 
1038
                        // if(_obstacle[i] && _pressure[i] != 0.0)
 
1039
                        //      printf("BAD PRESSURE i\n");
 
1040
 
 
1041
                        // if(_pressure[i]>1)
 
1042
                        //      printf("index: %d\n", i);
 
1043
                }
 
1044
                // printf("Max pressure: %f, dx: %f\n", maxvalue, _dx);
 
1045
        }
 
1046
 
834
1047
        setObstaclePressure(_pressure, 0, _zRes);
835
1048
 
836
1049
        // project out solution
 
1050
        // New idea for code from NVIDIA graphic gems 3 - DG
837
1051
        float invDx = 1.0f / _dx;
838
1052
        index = _slabSize + _xRes + 1;
839
1053
        for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
840
1054
                for (y = 1; y < _yRes - 1; y++, index += 2)
841
1055
                        for (x = 1; x < _xRes - 1; x++, index++)
842
1056
                        {
 
1057
                                float vMask[3] = {1.0f, 1.0f, 1.0f}, vObst[3] = {0, 0, 0};
 
1058
                                // float vR = 0.0f, vL = 0.0f, vT = 0.0f, vB = 0.0f, vD = 0.0f, vU = 0.0f;  // UNUSED
 
1059
 
 
1060
                                float pC = _pressure[index]; // center
 
1061
                                float pR = _pressure[index + 1]; // right
 
1062
                                float pL = _pressure[index - 1]; // left
 
1063
                                float pU = _pressure[index + _xRes]; // Up
 
1064
                                float pD = _pressure[index - _xRes]; // Down
 
1065
                                float pT = _pressure[index + _slabSize]; // top
 
1066
                                float pB = _pressure[index - _slabSize]; // bottom
 
1067
 
843
1068
                                if(!_obstacles[index])
844
1069
                                {
845
 
                                        _xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
846
 
                                        _yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
847
 
                                        _zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
 
1070
                                        // DG TODO: What if obstacle is left + right and one of them is moving?
 
1071
                                        if(_obstacles[index+1])                 { pR = pC; vObst[0] = _xVelocityOb[index + 1];                  vMask[0] = 0; }
 
1072
                                        if(_obstacles[index-1])                 { pL = pC; vObst[0]     = _xVelocityOb[index - 1];                      vMask[0] = 0; }
 
1073
                                        if(_obstacles[index+_xRes])             { pU = pC; vObst[1]     = _yVelocityOb[index + _xRes];          vMask[1] = 0; }
 
1074
                                        if(_obstacles[index-_xRes])             { pD = pC; vObst[1]     = _yVelocityOb[index - _xRes];          vMask[1] = 0; }
 
1075
                                        if(_obstacles[index+_slabSize]) { pT = pC; vObst[2] = _zVelocityOb[index + _slabSize];  vMask[2] = 0; }
 
1076
                                        if(_obstacles[index-_slabSize]) { pB = pC; vObst[2]     = _zVelocityOb[index - _slabSize];      vMask[2] = 0; }
 
1077
 
 
1078
                                        _xVelocity[index] -= 0.5f * (pR - pL) * invDx;
 
1079
                                        _yVelocity[index] -= 0.5f * (pU - pD) * invDx;
 
1080
                                        _zVelocity[index] -= 0.5f * (pT - pB) * invDx;
 
1081
 
 
1082
                                        _xVelocity[index] = (vMask[0] * _xVelocity[index]) + vObst[0];
 
1083
                                        _yVelocity[index] = (vMask[1] * _yVelocity[index]) + vObst[1];
 
1084
                                        _zVelocity[index] = (vMask[2] * _zVelocity[index]) + vObst[2];
 
1085
                                }
 
1086
                                else
 
1087
                                {
 
1088
                                        _xVelocity[index] = _xVelocityOb[index];
 
1089
                                        _yVelocity[index] = _yVelocityOb[index];
 
1090
                                        _zVelocity[index] = _zVelocityOb[index];
848
1091
                                }
849
1092
                        }
850
1093
 
 
1094
        // DG: was enabled in original code but now we do this later
 
1095
        // setObstacleVelocity(0, _zRes);
 
1096
 
851
1097
        if (_pressure) delete[] _pressure;
852
1098
        if (_divergence) delete[] _divergence;
853
1099
}
854
1100
 
855
 
 
856
 
 
 
1101
//////////////////////////////////////////////////////////////////////
 
1102
// calculate the obstacle velocity at boundary
 
1103
//////////////////////////////////////////////////////////////////////
 
1104
void FLUID_3D::setObstacleVelocity(int zBegin, int zEnd)
 
1105
{
 
1106
        
 
1107
        // completely TODO <-- who wrote this and what is here TODO? DG
 
1108
 
 
1109
        const size_t index_ = _slabSize + _xRes + 1;
 
1110
 
 
1111
        //int vIndex=_slabSize + _xRes + 1;
 
1112
 
 
1113
        int bb=0;
 
1114
        int bt=0;
 
1115
 
 
1116
        if (zBegin == 0) {bb = 1;}
 
1117
        if (zEnd == _zRes) {bt = 1;}
 
1118
 
 
1119
        // tag remaining obstacle blocks
 
1120
        for (int z = zBegin + bb; z < zEnd - bt; z++)
 
1121
        {
 
1122
                size_t index = index_ +(z-1)*_slabSize;
 
1123
 
 
1124
                for (int y = 1; y < _yRes - 1; y++, index += 2)
 
1125
                {
 
1126
                        for (int x = 1; x < _xRes - 1; x++, index++)
 
1127
                {
 
1128
                        if (!_obstacles[index])
 
1129
                        {
 
1130
                                // if(_obstacles[index+1]) xright = - _xVelocityOb[index]; 
 
1131
                                if((_obstacles[index - 1] & 8) && abs(_xVelocityOb[index - 1]) > FLT_EPSILON )
 
1132
                                {
 
1133
                                        // printf("velocity x!\n");
 
1134
                                        _xVelocity[index]  = _xVelocityOb[index - 1];
 
1135
                                        _xVelocity[index - 1]  = _xVelocityOb[index - 1];
 
1136
                                }
 
1137
                                // if(_obstacles[index+_xRes]) yup    = - _yVelocityOb[index];
 
1138
                                if((_obstacles[index - _xRes] & 8) && abs(_yVelocityOb[index - _xRes]) > FLT_EPSILON)
 
1139
                                {
 
1140
                                        // printf("velocity y!\n");
 
1141
                                        _yVelocity[index]  = _yVelocityOb[index - _xRes];
 
1142
                                        _yVelocity[index - _xRes]  = _yVelocityOb[index - _xRes];
 
1143
                                }
 
1144
                                // if(_obstacles[index+_slabSize]) ztop    = - _zVelocityOb[index];
 
1145
                                if((_obstacles[index - _slabSize] & 8) && abs(_zVelocityOb[index - _slabSize]) > FLT_EPSILON)
 
1146
                                {
 
1147
                                        // printf("velocity z!\n");
 
1148
                                        _zVelocity[index] = _zVelocityOb[index - _slabSize];
 
1149
                                        _zVelocity[index - _slabSize] = _zVelocityOb[index - _slabSize];
 
1150
                                }
 
1151
                        }
 
1152
                        else
 
1153
                        {
 
1154
                                _density[index] = 0;
 
1155
                        }
 
1156
                        //vIndex++;
 
1157
                }       // x loop
 
1158
                //vIndex += 2;
 
1159
                }       // y loop
 
1160
                //vIndex += 2 * _xRes;
 
1161
        }       // z loop
 
1162
}
857
1163
 
858
1164
//////////////////////////////////////////////////////////////////////
859
1165
// diffuse heat
869
1175
        for (int x = 0; x < _totalCells; x++)
870
1176
                if (_obstacles[x])
871
1177
                        _heat[x] = 0.0f;
872
 
 
873
1178
}
874
1179
 
875
1180
//////////////////////////////////////////////////////////////////////
892
1197
void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd)
893
1198
{
894
1199
 
895
 
        // compleately TODO
 
1200
        // completely TODO <-- who wrote this and what is here TODO? DG
896
1201
 
897
1202
        const size_t index_ = _slabSize + _xRes + 1;
898
1203
 
914
1219
                        for (int x = 1; x < _xRes - 1; x++, index++)
915
1220
                {
916
1221
                        // could do cascade of ifs, but they are a pain
917
 
                        if (_obstacles[index])
 
1222
                        if (_obstacles[index] /* && !(_obstacles[index] & 8) DG TODO TEST THIS CONDITION */)
918
1223
                        {
919
1224
                                const int top   = _obstacles[index + _slabSize];
920
1225
                                const int bottom= _obstacles[index - _slabSize];
928
1233
                                // const bool fully = (up && down);
929
1234
                                //const bool fullx = (left && right);
930
1235
 
 
1236
                                /*
931
1237
                                _xVelocity[index] =
932
1238
                                _yVelocity[index] =
933
1239
                                _zVelocity[index] = 0.0f;
 
1240
                                */
934
1241
                                _pressure[index] = 0.0f;
935
1242
 
936
1243
                                // average pressure neighbors
1035
1342
//////////////////////////////////////////////////////////////////////
1036
1343
// add buoyancy forces
1037
1344
//////////////////////////////////////////////////////////////////////
1038
 
void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
 
1345
void FLUID_3D::addBuoyancy(float *heat, float *density, float gravity[3], int zBegin, int zEnd)
1039
1346
{
1040
1347
        int index = zBegin*_slabSize;
1041
1348
 
1043
1350
                for (int y = 0; y < _yRes; y++)
1044
1351
                        for (int x = 0; x < _xRes; x++, index++)
1045
1352
                        {
1046
 
                                _zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender
 
1353
                                float buoyancy = *_alpha * density[index] + (*_beta * (((heat) ? heat[index] : 0.0f) - _tempAmb));
 
1354
                                _xForce[index] -= gravity[0] * buoyancy;
 
1355
                                _yForce[index] -= gravity[1] * buoyancy;
 
1356
                                _zForce[index] -= gravity[2] * buoyancy;
1047
1357
                        }
1048
1358
}
1049
1359
 
 
1360
 
1050
1361
//////////////////////////////////////////////////////////////////////
1051
1362
// add vorticity to the force field
1052
1363
//////////////////////////////////////////////////////////////////////
 
1364
#define VORT_VEL(i, j) \
 
1365
        ((_obstacles[obpos[(i)]] & 8) ? ((abs(objvelocity[(j)][obpos[(i)]]) > FLT_EPSILON) ? objvelocity[(j)][obpos[(i)]] : velocity[(j)][index]) : velocity[(j)][obpos[(i)]])
 
1366
 
1053
1367
void FLUID_3D::addVorticity(int zBegin, int zEnd)
1054
1368
{
 
1369
        // set flame vorticity from RNA value
 
1370
        float flame_vorticity = (*_flame_vorticity)/_constantScaling;
1055
1371
        //int x,y,z,index;
1056
 
        if(_vorticityEps<=0.) return;
 
1372
        if(_vorticityEps+flame_vorticity<=0.) return;
1057
1373
 
1058
1374
        int _blockSize=zEnd-zBegin;
1059
1375
        int _blockTotalCells = _slabSize * (_blockSize+2);
1085
1401
        float gridSize = 0.5f / _dx;
1086
1402
        //index = _slabSize + _xRes + 1;
1087
1403
 
 
1404
        float *velocity[3];
 
1405
        float *objvelocity[3];
 
1406
 
 
1407
        velocity[0] = _xVelocity;
 
1408
        velocity[1] = _yVelocity;
 
1409
        velocity[2] = _zVelocity;
 
1410
 
 
1411
        objvelocity[0] = _xVelocityOb;
 
1412
        objvelocity[1] = _yVelocityOb;
 
1413
        objvelocity[2] = _zVelocityOb;
1088
1414
 
1089
1415
        size_t vIndex=_xRes + 1;
1090
 
 
1091
1416
        for (int z = zBegin + bb1; z < (zEnd - bt1); z++)
1092
1417
        {
1093
1418
                size_t index = index_ +(z-1)*_slabSize;
1097
1422
                {
1098
1423
                        for (int x = 1; x < _xRes - 1; x++, index++)
1099
1424
                        {
1100
 
 
1101
 
                                int up    = _obstacles[index + _xRes] ? index : index + _xRes;
1102
 
                                int down  = _obstacles[index - _xRes] ? index : index - _xRes;
1103
 
                                float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
1104
 
                                int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
1105
 
                                int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
1106
 
                                float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
1107
 
                                int right = _obstacles[index + 1] ? index : index + 1;
1108
 
                                int left  = _obstacles[index - 1] ? index : index - 1;
1109
 
                                float dx  = (right == index || left == index) ? 1.0f / _dx : gridSize;
1110
 
 
1111
 
                                _xVorticity[vIndex] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
1112
 
                                _yVorticity[vIndex] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
1113
 
                                _zVorticity[vIndex] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
1114
 
 
1115
 
                                _vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
 
1425
                                if (!_obstacles[index])
 
1426
                                {
 
1427
                                        int obpos[6];
 
1428
 
 
1429
                                        obpos[0] = (_obstacles[index + _xRes] == 1) ? index : index + _xRes; // up
 
1430
                                        obpos[1] = (_obstacles[index - _xRes] == 1) ? index : index - _xRes; // down
 
1431
                                        float dy = (obpos[0] == index || obpos[1] == index) ? 1.0f / _dx : gridSize;
 
1432
 
 
1433
                                        obpos[2]  = (_obstacles[index + _slabSize] == 1) ? index : index + _slabSize; // out
 
1434
                                        obpos[3]  = (_obstacles[index - _slabSize] == 1) ? index : index - _slabSize; // in
 
1435
                                        float dz  = (obpos[2] == index || obpos[3] == index) ? 1.0f / _dx : gridSize;
 
1436
 
 
1437
                                        obpos[4] = (_obstacles[index + 1] == 1) ? index : index + 1; // right
 
1438
                                        obpos[5] = (_obstacles[index - 1] == 1) ? index : index - 1; // left
 
1439
                                        float dx = (obpos[4] == index || obpos[5] == index) ? 1.0f / _dx : gridSize;
 
1440
 
 
1441
                                        float xV[2], yV[2], zV[2];
 
1442
 
 
1443
                                        zV[1] = VORT_VEL(0, 2);
 
1444
                                        zV[0] = VORT_VEL(1, 2);
 
1445
                                        yV[1] = VORT_VEL(2, 1);
 
1446
                                        yV[0] = VORT_VEL(3, 1);
 
1447
                                        _xVorticity[vIndex] = (zV[1] - zV[0]) * dy + (-yV[1] + yV[0]) * dz;
 
1448
 
 
1449
                                        xV[1] = VORT_VEL(2, 0);
 
1450
                                        xV[0] = VORT_VEL(3, 0);
 
1451
                                        zV[1] = VORT_VEL(4, 2);
 
1452
                                        zV[0] = VORT_VEL(5, 2);
 
1453
                                        _yVorticity[vIndex] = (xV[1] - xV[0]) * dz + (-zV[1] + zV[0]) * dx;
 
1454
 
 
1455
                                        yV[1] = VORT_VEL(4, 1);
 
1456
                                        yV[0] = VORT_VEL(5, 1);
 
1457
                                        xV[1] = VORT_VEL(0, 0);
 
1458
                                        xV[0] = VORT_VEL(1, 0);
 
1459
                                        _zVorticity[vIndex] = (yV[1] - yV[0]) * dx + (-xV[1] + xV[0])* dy;
 
1460
 
 
1461
                                        _vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
1116
1462
                                                _yVorticity[vIndex] * _yVorticity[vIndex] +
1117
1463
                                                _zVorticity[vIndex] * _zVorticity[vIndex]);
1118
1464
 
 
1465
                                }
1119
1466
                                vIndex++;
1120
1467
                        }
1121
1468
                        vIndex+=2;
1144
1491
                                {
1145
1492
                                        float N[3];
1146
1493
 
1147
 
                                        int up    = _obstacles[index + _xRes] ? vIndex : vIndex + _xRes;
1148
 
                                        int down  = _obstacles[index - _xRes] ? vIndex : vIndex - _xRes;
 
1494
                                        int up    = (_obstacles[index + _xRes] == 1) ? vIndex : vIndex + _xRes;
 
1495
                                        int down  = (_obstacles[index - _xRes] == 1) ? vIndex : vIndex - _xRes;
1149
1496
                                        float dy  = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
1150
 
                                        int out   = _obstacles[index + _slabSize] ? vIndex : vIndex + _slabSize;
1151
 
                                        int in    = _obstacles[index - _slabSize] ? vIndex : vIndex - _slabSize;
 
1497
 
 
1498
                                        int out   = (_obstacles[index + _slabSize] == 1) ? vIndex : vIndex + _slabSize;
 
1499
                                        int in    = (_obstacles[index - _slabSize] == 1) ? vIndex : vIndex - _slabSize;
1152
1500
                                        float dz  = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
1153
 
                                        int right = _obstacles[index + 1] ? vIndex : vIndex + 1;
1154
 
                                        int left  = _obstacles[index - 1] ? vIndex : vIndex - 1;
 
1501
 
 
1502
                                        int right = (_obstacles[index + 1] == 1) ? vIndex : vIndex + 1;
 
1503
                                        int left  = (_obstacles[index - 1] == 1) ? vIndex : vIndex - 1;
1155
1504
                                        float dx  = (right == vIndex || left == vIndex) ? 1.0f / _dx : gridSize;
 
1505
 
1156
1506
                                        N[0] = (_vorticity[right] - _vorticity[left]) * dx;
1157
1507
                                        N[1] = (_vorticity[up] - _vorticity[down]) * dy;
1158
1508
                                        N[2] = (_vorticity[out] - _vorticity[in]) * dz;
1159
1509
 
1160
1510
                                        float magnitude = sqrtf(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1161
 
                                        if (magnitude > 0.0f)
 
1511
                                        if (magnitude > FLT_EPSILON)
1162
1512
                                        {
 
1513
                                                float flame_vort = (_fuel) ? _fuel[index]*flame_vorticity : 0.0f;
1163
1514
                                                magnitude = 1.0f / magnitude;
1164
1515
                                                N[0] *= magnitude;
1165
1516
                                                N[1] *= magnitude;
1166
1517
                                                N[2] *= magnitude;
1167
1518
 
1168
 
                                                _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * eps;
1169
 
                                                _yForce[index] -= (N[0] * _zVorticity[vIndex] - N[2] * _xVorticity[vIndex]) * _dx * eps;
1170
 
                                                _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
 
1519
                                                _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * (eps + flame_vort);
 
1520
                                                _yForce[index] += (N[2] * _xVorticity[vIndex] - N[0] * _zVorticity[vIndex]) * _dx * (eps + flame_vort);
 
1521
                                                _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * (eps + flame_vort);
1171
1522
                                        }
1172
1523
                                        }       // if
1173
1524
                                        vIndex++;
1188
1539
{
1189
1540
        Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1190
1541
 
1191
 
        if(_domainBcLeft == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
1192
 
        else setZeroX(_xVelocityOld, res, zBegin, zEnd);
1193
 
 
1194
 
        if(_domainBcFront == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
1195
 
        else setZeroY(_yVelocityOld, res, zBegin, zEnd); 
1196
 
 
1197
 
        if(_domainBcTop == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
1198
 
        else setZeroZ(_zVelocityOld, res, zBegin, zEnd);
 
1542
        setZeroX(_xVelocityOld, res, zBegin, zEnd);
 
1543
        setZeroY(_yVelocityOld, res, zBegin, zEnd);
 
1544
        setZeroZ(_zVelocityOld, res, zBegin, zEnd);
1199
1545
}
1200
1546
 
1201
1547
//////////////////////////////////////////////////////////////////////
1215
1561
        // advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res)
1216
1562
 
1217
1563
        advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, res, zBegin, zEnd);
1218
 
        advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
 
1564
        if (_heat) {
 
1565
                advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
 
1566
        }
 
1567
        if (_fuel) {
 
1568
                advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _fuelOld, _fuelTemp, res, zBegin, zEnd);
 
1569
                advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _reactOld, _reactTemp, res, zBegin, zEnd);
 
1570
        }
 
1571
        if (_color_r) {
 
1572
                advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_rOld, _color_rTemp, res, zBegin, zEnd);
 
1573
                advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_gOld, _color_gTemp, res, zBegin, zEnd);
 
1574
                advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_bOld, _color_bTemp, res, zBegin, zEnd);
 
1575
        }
1219
1576
        advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, res, zBegin, zEnd);
1220
1577
        advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, res, zBegin, zEnd);
1221
1578
        advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, res, zBegin, zEnd);
1231
1588
        const float dt0 = _dt / _dx;
1232
1589
        Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1233
1590
 
1234
 
        // use force array as temp arrays
 
1591
        // use force array as temp array
1235
1592
        float* t1 = _xForce;
1236
1593
 
1237
1594
        // advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles)
1238
1595
 
 
1596
        /* finish advection */
1239
1597
        advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd);
1240
 
        advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
 
1598
        if (_heat) {
 
1599
                advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
 
1600
        }
 
1601
        if (_fuel) {
 
1602
                advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _fuelOld, _fuel, _fuelTemp, t1, res, _obstacles, zBegin, zEnd);
 
1603
                advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _reactOld, _react, _reactTemp, t1, res, _obstacles, zBegin, zEnd);
 
1604
        }
 
1605
        if (_color_r) {
 
1606
                advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_rOld, _color_r, _color_rTemp, t1, res, _obstacles, zBegin, zEnd);
 
1607
                advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_gOld, _color_g, _color_gTemp, t1, res, _obstacles, zBegin, zEnd);
 
1608
                advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_bOld, _color_b, _color_bTemp, t1, res, _obstacles, zBegin, zEnd);
 
1609
        }
1241
1610
        advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocityTemp, _xVelocity, t1, res, _obstacles, zBegin, zEnd);
1242
1611
        advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd);
1243
1612
        advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd);
1244
1613
 
1245
 
        if(_domainBcLeft == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
 
1614
        /* set boundary conditions for velocity */
 
1615
        if(!_domainBcLeft) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
1246
1616
        else setZeroX(_xVelocityTemp, res, zBegin, zEnd);                               
1247
1617
 
1248
 
        if(_domainBcFront == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
 
1618
        if(!_domainBcFront) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
1249
1619
        else setZeroY(_yVelocityTemp, res, zBegin, zEnd); 
1250
1620
 
1251
 
        if(_domainBcTop == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
 
1621
        if(!_domainBcTop) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
1252
1622
        else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
1253
1623
 
 
1624
        /* clear data boundaries */
1254
1625
        setZeroBorder(_density, res, zBegin, zEnd);
1255
 
        setZeroBorder(_heat, res, zBegin, zEnd);
1256
 
 
1257
 
 
1258
 
        /*int begin=zBegin * _slabSize;
1259
 
        int end=begin + (zEnd - zBegin) * _slabSize;
1260
 
  for (int x = begin; x < end; x++)
1261
 
    _xForce[x] = _yForce[x] = 0.0f;*/
 
1626
        if (_fuel) {
 
1627
                setZeroBorder(_fuel, res, zBegin, zEnd);
 
1628
                setZeroBorder(_react, res, zBegin, zEnd);
 
1629
        }
 
1630
        if (_color_r) {
 
1631
                setZeroBorder(_color_r, res, zBegin, zEnd);
 
1632
                setZeroBorder(_color_g, res, zBegin, zEnd);
 
1633
                setZeroBorder(_color_b, res, zBegin, zEnd);
 
1634
        }
 
1635
}
 
1636
 
 
1637
 
 
1638
void FLUID_3D::processBurn(float *fuel, float *smoke, float *react, float *flame, float *heat,
 
1639
                                                   float *r, float *g, float *b, int total_cells, float dt)
 
1640
{
 
1641
        float burning_rate = *_burning_rate;
 
1642
        float flame_smoke = *_flame_smoke;
 
1643
        float ignition_point = *_ignition_temp;
 
1644
        float temp_max = *_max_temp;
 
1645
 
 
1646
        for (int index = 0; index < total_cells; index++)
 
1647
        {
 
1648
                float orig_fuel = fuel[index];
 
1649
                float orig_smoke = smoke[index];
 
1650
                float smoke_emit = 0.0f;
 
1651
                float react_coord = 0.0f;
 
1652
 
 
1653
                /* process fuel */
 
1654
                fuel[index] -= burning_rate * dt;
 
1655
                if (fuel[index] < 0.0f) fuel[index] = 0.0f;
 
1656
                /* process reaction coordinate */
 
1657
                if (orig_fuel > FLT_EPSILON) {
 
1658
                        react[index] *= fuel[index]/orig_fuel;
 
1659
                        react_coord = react[index];
 
1660
                }
 
1661
                else {
 
1662
                        react[index] = 0.0f;
 
1663
                }
 
1664
 
 
1665
                /* emit smoke based on fuel burn rate and "flame_smoke" factor */
 
1666
                smoke_emit = (orig_fuel < 1.0f) ? (1.0f - orig_fuel)*0.5f : 0.0f;
 
1667
                smoke_emit = (smoke_emit + 0.5f) * (orig_fuel-fuel[index]) * 0.1f * flame_smoke;
 
1668
                smoke[index] += smoke_emit;
 
1669
                CLAMP(smoke[index], 0.0f, 1.0f);
 
1670
 
 
1671
                /* model flame temperature curve from the reaction coordinate (fuel) */
 
1672
                if (react_coord>0.0f) {
 
1673
                        /* do a smooth falloff for rest of the values */
 
1674
                        flame[index] = pow(react_coord, 0.5f);
 
1675
                }
 
1676
                else
 
1677
                        flame[index] = 0.0f;
 
1678
 
 
1679
                /* set fluid temperature from the flame temperature profile */
 
1680
                if (heat && flame[index])
 
1681
                        heat[index] = (1.0f-flame[index])*ignition_point + flame[index]*temp_max;
 
1682
 
 
1683
                /* mix new color */
 
1684
                if (r && smoke_emit > FLT_EPSILON) {
 
1685
                        float smoke_factor = smoke[index]/(orig_smoke+smoke_emit);
 
1686
                        r[index] = (r[index] + _flame_smoke_color[0] * smoke_emit) * smoke_factor;
 
1687
                        g[index] = (g[index] + _flame_smoke_color[1] * smoke_emit) * smoke_factor;
 
1688
                        b[index] = (b[index] + _flame_smoke_color[2] * smoke_emit) * smoke_factor;
 
1689
                }
 
1690
        }
1262
1691
}