Nilorea Library
C utilities for networking, threading, graphics
n_fluids.c
Go to the documentation of this file.
1
7#include <math.h>
8#include <stdio.h>
9#include <stdlib.h>
10#include <stdint.h>
11#include <stdbool.h>
12#include <strings.h>
13#include "nilorea/n_fluids.h"
14#include "nilorea/n_common.h"
15#include "nilorea/n_thread_pool.h"
16
18#define N_FLUID_U_FIELD 0
20#define N_FLUID_V_FIELD 1
22#define N_FLUID_S_FIELD 2
23
24
32void n_memset( void *dst , void *val , size_t size , size_t count )
33{
34 void *ptr = dst;
35 while( count-- > 0 )
36 {
37 memcpy( ptr , val , size );
38 ptr += size;
39 }
40}
41
42
48int destroy_n_fluid( N_FLUID **fluid )
49{
50 __n_assert( (*fluid) , return FALSE );
51
52 FreeNoLog( (*fluid) -> u );
53 FreeNoLog( (*fluid) -> newU );
54 FreeNoLog( (*fluid) -> v );
55 FreeNoLog( (*fluid) -> newV );
56 FreeNoLog( (*fluid) -> p );
57 FreeNoLog( (*fluid) -> s );
58 FreeNoLog( (*fluid) -> m );
59 FreeNoLog( (*fluid) -> newM );
60 FreeNoLog( (*fluid) );
61
62 return TRUE ;
63}
64
65
66
78N_FLUID *new_n_fluid( double density , double gravity , size_t numIters , double dt , double overRelaxation , size_t sx , size_t sy )
79{
80 N_FLUID *fluid = NULL ;
81
82 Malloc( fluid , N_FLUID , 1 );
83 __n_assert( fluid , return NULL );
84
85 fluid -> density = density ;
86 fluid -> gravity = gravity ;
87 fluid -> numIters = numIters ;
88 fluid -> dt = dt ;
89 fluid -> h = 1.0 / 100.0 ;
90 fluid -> overRelaxation = overRelaxation ;
91 fluid -> numX = sx + 2 ;
92 fluid -> numY = sy + 2 ;
93 fluid -> numZ = 1 ; // TODO: expand Z
94 fluid -> numCells = fluid -> numX * fluid -> numY * fluid -> numZ;
95
96 Malloc( fluid -> u , double , fluid -> numCells );
97 Malloc( fluid -> newU , double , fluid -> numCells );
98 Malloc( fluid -> v , double , fluid -> numCells );
99 Malloc( fluid -> newV , double , fluid -> numCells );
100 Malloc( fluid -> p , double , fluid -> numCells );
101 Malloc( fluid -> s , double , fluid -> numCells );
102 Malloc( fluid -> m , double , fluid -> numCells );
103 Malloc( fluid -> newM , double , fluid -> numCells );
104
105 fluid -> showSmoke = 1 ;
106 fluid -> showPressure = 0 ;
107 fluid -> showPaint = 0 ;
108 fluid -> fluid_production_percentage = 0.1 ;
109 fluid -> cScale = 16.0 ;
110 // double precision. Taking 'value', if( fabs( value ) < float_tolerance ) value is considered as zero
111 fluid -> negative_float_tolerance = -0.00001 ;
112 fluid -> positive_float_tolerance = 0.00001 ;
113
114 double d_val = 1.0 ;
115 n_memset( fluid -> m , &d_val , sizeof( d_val ) , fluid -> numCells );
116
117 // precalculate and allocated N_PROC_PARAMS lists for threaded computing
118 fluid -> integrate_chunk_list = new_generic_list( 0 );
119 fluid -> solveIncompressibility_chunk_list = new_generic_list( 0 );
120 fluid -> advectVel_chunk_list = new_generic_list( 0 );
121 fluid -> advectSmoke_chunk_list = new_generic_list( 0 );
122
123 int nb_cores = get_nb_cpu_cores();
124 if( nb_cores == -1 )
125 {
126 nb_cores = 1 ;
127 }
128 size_t steps = fluid -> numX / nb_cores ;
129
130 N_FLUID_THREAD_PARAMS *params = NULL ;
131 // integrate
132 for( size_t i = 1 ; i < fluid -> numX ; i +=steps )
133 {
134 Malloc( params , N_FLUID_THREAD_PARAMS , 1 );
135 params -> ptr = fluid ;
136 params -> x_start = i ;
137 params -> x_end = i+steps ;
138 params -> y_start = 1 ;
139 params -> y_end = fluid -> numY - 1 ;
140 list_push( fluid -> integrate_chunk_list , params , &free );
141 }
142 // set the last batch at the end of the range
143 params = (N_FLUID_THREAD_PARAMS *)fluid -> integrate_chunk_list -> end -> ptr ;
144 params -> x_end = fluid -> numX ;
145
146 // solveIncompressibility
147 for( size_t i = 1 ; i < fluid -> numX - 1 ; i +=steps )
148 {
149 Malloc( params , N_FLUID_THREAD_PARAMS , 1 );
150 params -> ptr = fluid ;
151 params -> x_start = i ;
152 params -> x_end = i+steps ;
153 params -> y_start = 1 ;
154 params -> y_end = fluid -> numY - 1 ;
155 list_push( fluid -> solveIncompressibility_chunk_list , params , &free );
156 }
157 // set the last batch at the end of the range
158 params = (N_FLUID_THREAD_PARAMS *)fluid -> solveIncompressibility_chunk_list -> end -> ptr ;
159 params -> x_end = fluid -> numX - 1 ;
160
161 // advectVel
162 for( size_t i = 1 ; i < fluid -> numX ; i +=steps )
163 {
164 Malloc( params , N_FLUID_THREAD_PARAMS , 1 );
165 params -> ptr = fluid ;
166 params -> x_start = i ;
167 params -> x_end = i+steps ;
168 params -> y_start = 1 ;
169 params -> y_end = fluid -> numY ;
170 list_push( fluid -> advectVel_chunk_list , params , &free );
171 }
172 // set the last batch at the end of the range
173 params = (N_FLUID_THREAD_PARAMS *)fluid -> advectVel_chunk_list -> end -> ptr ;
174 params -> x_end = fluid -> numX ;
175
176
177 // adVectSmoke
178 for( size_t i = 1 ; i < fluid -> numX -1 ; i +=steps )
179 {
180 Malloc( params , N_FLUID_THREAD_PARAMS , 1 );
181 params -> ptr = fluid ;
182 params -> x_start = i ;
183 params -> x_end = i+steps ;
184 params -> y_start = 1 ;
185 params -> y_end = fluid -> numY - 1;
186 list_push( fluid -> advectSmoke_chunk_list , params , &free );
187 }
188 // set the last batch at the end of the range
189 params = (N_FLUID_THREAD_PARAMS *)fluid -> advectSmoke_chunk_list -> end -> ptr ;
190 params -> x_end = fluid -> numX -1 ;
191
192 return fluid ;
193} /* new_n_fluid */
194
195
196
202void *n_fluid_integrate_proc( void *ptr )
203{
204 N_FLUID_THREAD_PARAMS *params = (N_FLUID_THREAD_PARAMS *)ptr ;
205 N_FLUID *fluid = (N_FLUID *)params -> ptr ;
206
207 size_t n = fluid -> numY;
208 for( size_t i = params -> x_start ; i < params -> x_end ; i++ )
209 {
210 for( size_t j = params -> y_start ; j < params -> y_end ; j++ )
211 {
212 if( !_z( fluid , s[ i*n + j ] ) && !_z( fluid , s[ i*n + j-1 ] ) )
213 fluid -> v[ i*n + j ] += fluid -> gravity * fluid -> dt ;
214 }
215 }
216 return NULL ;
217}
218
219
225int n_fluid_integrate( N_FLUID *fluid )
226{
227 __n_assert( fluid , return FALSE );
228
229 size_t n = fluid -> numY;
230 for( size_t i = 1 ; i < fluid -> numX ; i++ )
231 {
232 for( size_t j = 1; j < fluid -> numY - 1 ; j++ )
233 {
234 if( !_z( fluid , s[ i*n + j ] ) && !_z( fluid , s[ i*n + j-1 ] ) )
235 fluid -> v[ i*n + j ] += fluid -> gravity * fluid -> dt ;
236 }
237 }
238 return TRUE ;
239}
240
241
242
248void *n_fluid_solveIncompressibility_proc( void *ptr )
249{
250 N_FLUID_THREAD_PARAMS *params = (N_FLUID_THREAD_PARAMS *)ptr ;
251 N_FLUID *fluid = (N_FLUID *)params -> ptr ;
252
253 double cp = ( fluid -> density * fluid -> h ) / fluid -> dt ;
254
255 size_t n = fluid -> numY;
256 for( size_t i = params -> x_start ; i < params -> x_end ; i++ )
257 {
258 for( size_t j = params -> y_start ; j < params -> y_end ; j++ )
259 {
260 if( _z( fluid , s[ i*n + j ] ) )
261 continue;
262
263 double sx0 = fluid ->s[ (i - 1)*n + j ];
264 double sx1 = fluid ->s[ (i + 1)*n + j ];
265 double sy0 = fluid ->s[ i*n + j - 1 ];
266 double sy1 = fluid ->s[ i*n + j + 1 ];
267 double s = sx0 + sx1 + sy0 + sy1;
268 if( _zd( fluid , s ) )
269 continue;
270
271 double div = fluid ->u[(i+1)*n + j] - fluid ->u[i*n + j] + fluid ->v[i*n + j+1] - fluid ->v[i*n + j];
272 double p = ( -div * fluid -> overRelaxation ) / s ;
273 fluid ->p[i*n + j] += cp * p;
274 fluid ->u[i*n + j] -= sx0 * p;
275 fluid ->u[(i+1)*n + j] += sx1 * p;
276 fluid ->v[i*n + j] -= sy0 * p;
277 fluid ->v[i*n + j+1] += sy1 * p;
278 }
279 }
280 return NULL ;
281}
282
283
289int n_fluid_solveIncompressibility( N_FLUID *fluid )
290{
291 __n_assert( fluid , return FALSE );
292 size_t n = fluid -> numY ;
293
294 double cp = ( fluid -> density * fluid -> h ) / fluid -> dt ;
295
296 for( size_t iter = 0 ; iter < fluid -> numIters ; iter++ )
297 {
298 for( size_t i = 1 ; i < fluid -> numX - 1 ; i++ )
299 {
300 for( size_t j = 1 ; j < fluid -> numY - 1 ; j++ )
301 {
302 if( _z( fluid , s[ i*n + j ] ) )
303 continue;
304
305 double sx0 = fluid ->s[ (i - 1)*n + j ];
306 double sx1 = fluid ->s[ (i + 1)*n + j ];
307 double sy0 = fluid ->s[ i*n + j - 1 ];
308 double sy1 = fluid ->s[ i*n + j + 1 ];
309 double s = sx0 + sx1 + sy0 + sy1;
310 if( _zd( fluid , s ) )
311 continue;
312
313 double div = fluid ->u[(i+1)*n + j] - fluid ->u[i*n + j] + fluid ->v[i*n + j+1] - fluid ->v[i*n + j];
314 double p = ( -div * fluid -> overRelaxation ) / s ;
315 fluid ->p[i*n + j] += cp * p;
316 fluid ->u[i*n + j] -= sx0 * p;
317 fluid ->u[(i+1)*n + j] += sx1 * p;
318 fluid ->v[i*n + j] -= sy0 * p;
319 fluid ->v[i*n + j+1] += sy1 * p;
320 }
321 }
322 }
323 return TRUE ;
324}
325
326
327
333int n_fluid_extrapolate( N_FLUID *fluid )
334{
335 __n_assert( fluid , return FALSE );
336 size_t n = fluid -> numY ;
337 for( size_t i = 0; i < fluid ->numX ; i++ )
338 {
339 fluid ->u[ i*n + 0 ] = fluid ->u[i*n + 1];
340 fluid ->u[ i*n + fluid ->numY-1 ] = fluid ->u[i*n + fluid ->numY-2];
341 }
342 for( size_t j = 0; j < fluid ->numY; j++)
343 {
344 fluid ->v[0*n + j] = fluid ->v[1*n + j];
345 fluid ->v[(fluid ->numX-1)*n + j] = fluid ->v[(fluid ->numX-2)*n + j] ;
346 }
347 return TRUE ;
348}
349
350
359double n_fluid_sampleField( N_FLUID *fluid , double x , double y , uint32_t field )
360{
361 __n_assert( fluid , return FALSE );
362 size_t n = fluid -> numY;
363 double h1 = 1.0 / fluid -> h ;
364 double h2 = 0.5 * fluid -> h ;
365
366 x = MAX( MIN( x , fluid ->numX * fluid -> h ) , fluid -> h );
367 y = MAX( MIN( y , fluid ->numY * fluid -> h ) , fluid -> h );
368
369 double dx = 0.0;
370 double dy = 0.0;
371
372 double *f = NULL ;
373 switch( field )
374 {
375 case N_FLUID_U_FIELD: f = fluid ->u ; dy = h2 ; break ;
376 case N_FLUID_V_FIELD: f = fluid ->v ; dx = h2 ; break ;
377 case N_FLUID_S_FIELD: f = fluid ->m ; dx = h2 ; dy = h2 ; break ;
378 }
379
380 double x0 = MIN( floor( (x-dx) * h1 ) , fluid ->numX - 1 );
381 double tx = ((x-dx) - x0 * fluid -> h ) * h1;
382 double x1 = MIN(x0 + 1, fluid ->numX-1);
383
384 double y0 = MIN( floor( (y-dy)*h1 ) , fluid ->numY - 1 );
385 double ty = ((y-dy) - y0 * fluid -> h ) * h1;
386 double y1 = MIN( y0 + 1 , fluid ->numY-1 );
387
388 double sx = 1.0 - tx;
389 double sy = 1.0 - ty;
390
391 double val = sx*sy * f[(size_t)(x0*n + y0)] +
392 tx*sy * f[(size_t)(x1*n + y0)] +
393 tx*ty * f[(size_t)(x1*n + y1)] +
394 sx*ty * f[(size_t)(x0*n + y1)];
395
396 return val;
397}
398
399
400
408double n_fluid_avgU( N_FLUID *fluid , size_t i , size_t j )
409{
410 __n_assert( fluid , return FALSE );
411 size_t n = fluid ->numY ;
412 double u = (fluid ->u[i*n + j-1] + fluid ->u[i*n + j] +
413 fluid ->u[(i+1)*n + j-1] + fluid ->u[(i+1)*n + j]) * 0.25;
414
415 return u;
416
417}
418
419
420
428double n_fluid_avgV( N_FLUID *fluid , size_t i , size_t j )
429{
430 __n_assert( fluid , return FALSE );
431 size_t n = fluid ->numY;
432 double v = (fluid ->v[(i-1)*n + j] + fluid ->v[i*n + j] +
433 fluid ->v[(i-1)*n + j+1] + fluid ->v[i*n + j+1]) * 0.25;
434 return v;
435}
436
437
438
444void *n_fluid_advectVel_proc( void *ptr )
445{
446 N_FLUID_THREAD_PARAMS *params = (N_FLUID_THREAD_PARAMS *)ptr ;
447 N_FLUID *fluid = (N_FLUID *)params -> ptr ;
448
449 size_t n = fluid ->numY;
450 double h2 = 0.5 * fluid -> h ;
451 for( size_t i = params -> x_start ; i < params -> x_end ; i++ )
452 {
453 for( size_t j = params -> y_start ; j < params -> y_end ; j++ )
454 {
455 size_t index = i*n + j ;
456 // u component
457 if( !_z( fluid , s[index] ) && !_z( fluid , s[ (i-1)*n + j ] ) && j < fluid -> numY - 1)
458 {
459 double x = i * fluid -> h;
460 double y = j * fluid -> h + h2;
461 double u = fluid -> u[index];
462 double v = n_fluid_avgV( fluid , i , j );
463 //double v = n_fluid_sampleField( fluid , x , y , N_FLUID_V_FIELD );
464 x = x - fluid -> dt * u;
465 y = y - fluid -> dt * v;
466 u = n_fluid_sampleField( fluid , x , y , N_FLUID_U_FIELD );
467 fluid ->newU[index] = u;
468 }
469 // v component
470 if( !_z( fluid , s[index] ) && !_z( fluid , s[index-1] ) && i < fluid ->numX - 1) {
471 double x = i * fluid -> h + h2;
472 double y = j * fluid -> h ;
473 double u = n_fluid_avgU( fluid , i, j );
474 // double u = n_fluid_sampleField( fluid , x , y , N_FLUID_U_FIELD );
475 double v = fluid ->v[index];
476 x = x - fluid -> dt*u;
477 y = y - fluid -> dt*v;
478 v = n_fluid_sampleField( fluid , x , y , N_FLUID_V_FIELD );
479 fluid ->newV[index] = v;
480 }
481 }
482 }
483 return NULL ;
484}
485
486
492int n_fluid_advectVel( N_FLUID *fluid )
493{
494 __n_assert( fluid , return FALSE );
495
496 memcpy( fluid -> newU , fluid -> u , fluid -> numCells * sizeof( double ) );
497 memcpy( fluid -> newV , fluid -> v , fluid -> numCells * sizeof( double ) );
498
499 size_t n = fluid ->numY;
500 double h2 = 0.5 * fluid -> h ;
501 for( size_t i = 1; i < fluid -> numX ; i++ )
502 {
503 for( size_t j = 1; j < fluid -> numY ; j++ )
504 {
505 size_t index = i*n + j ;
506 // u component
507 if( !_z( fluid , s[index] ) && !_z( fluid , s[ (i-1)*n + j ] ) && j < fluid -> numY - 1)
508 {
509 double x = i * fluid -> h;
510 double y = j * fluid -> h + h2;
511 double u = fluid -> u[index];
512 double v = n_fluid_avgV( fluid , i , j );
513 //double v = n_fluid_sampleField( fluid , x , y , N_FLUID_V_FIELD );
514 x = x - fluid -> dt * u;
515 y = y - fluid -> dt * v;
516 u = n_fluid_sampleField( fluid , x , y , N_FLUID_U_FIELD );
517 fluid ->newU[index] = u;
518 }
519 // v component
520 if( !_z( fluid , s[index] ) && !_z( fluid , s[index-1] ) && i < fluid ->numX - 1) {
521 double x = i * fluid -> h + h2;
522 double y = j * fluid -> h ;
523 double u = n_fluid_avgU( fluid , i, j );
524 // double u = n_fluid_sampleField( fluid , x , y , N_FLUID_U_FIELD );
525 double v = fluid ->v[index];
526 x = x - fluid -> dt*u;
527 y = y - fluid -> dt*v;
528 v = n_fluid_sampleField( fluid , x , y , N_FLUID_V_FIELD );
529 fluid ->newV[index] = v;
530 }
531 }
532 }
533 double *ptr = fluid -> u ;
534 fluid -> u = fluid -> newU ;
535 fluid -> newU = ptr ;
536
537 ptr = fluid -> v ;
538 fluid -> v = fluid -> newV ;
539 fluid -> newV = ptr ;
540
541 return TRUE ;
542}
543
544
550void *n_fluid_advectSmoke_proc( void *ptr )
551{
552 N_FLUID_THREAD_PARAMS *params = (N_FLUID_THREAD_PARAMS *)ptr ;
553 N_FLUID *fluid = (N_FLUID *)params -> ptr ;
554
555 size_t n = fluid ->numY;
556 double h2 = 0.5 * fluid -> h ;
557 for( size_t i = params -> x_start ; i < params -> x_end ; i++ )
558 {
559 for( size_t j = params -> y_start ; j < params -> y_end ; j++ )
560 {
561 size_t index = i*n + j ;
562 if( !_z( fluid , s[ index ] ) )
563 {
564 double u = (fluid ->u[ index ] + fluid ->u[ (i+1)*n + j ] ) * 0.5;
565 double v = (fluid ->v[ index ] + fluid ->v[ index + 1]) * 0.5;
566 double x = i * fluid -> h + h2 - fluid -> dt*u;
567 double y = j * fluid -> h + h2 - fluid -> dt*v;
568
569 fluid ->newM[index] = n_fluid_sampleField( fluid , x , y , N_FLUID_S_FIELD );
570 }
571 }
572 }
573 return NULL ;
574}
575
576
582int n_fluid_advectSmoke( N_FLUID *fluid )
583{
584 __n_assert( fluid , return FALSE );
585
586 size_t n = fluid ->numY;
587 double h2 = 0.5 * fluid -> h ;
588
589 memcpy( fluid -> newM , fluid -> m , fluid -> numCells * sizeof( double ) );
590
591 for( size_t i = 1; i < fluid -> numX - 1 ; i++ )
592 {
593 for( size_t j = 1; j < fluid -> numY - 1 ; j++ )
594 {
595 size_t index = i*n + j ;
596 if( !_z( fluid , s[ index ] ) )
597 {
598 double u = (fluid ->u[ index ] + fluid ->u[ (i+1)*n + j ] ) * 0.5;
599 double v = (fluid ->v[ index ] + fluid ->v[ index + 1]) * 0.5;
600 double x = i * fluid -> h + h2 - fluid -> dt*u;
601 double y = j * fluid -> h + h2 - fluid -> dt*v;
602
603 fluid ->newM[index] = n_fluid_sampleField( fluid , x , y , N_FLUID_S_FIELD );
604 }
605 }
606 }
607 double *ptr = fluid -> m ;
608 fluid -> m = fluid -> newM ;
609 fluid -> newM = ptr ;
610
611 return TRUE ;
612}
613
614
615
621int n_fluid_simulate( N_FLUID *fluid )
622{
623 __n_assert( fluid , return FALSE );
624 n_fluid_integrate( fluid );
625 memset( fluid -> p , 0 , fluid -> numCells * sizeof( double ) );
626 n_fluid_solveIncompressibility( fluid );
627 n_fluid_extrapolate( fluid );
628 n_fluid_advectVel( fluid );
629 n_fluid_advectSmoke( fluid );
630 return TRUE ;
631}
632
633
634
641int n_fluid_simulate_threaded( N_FLUID *fluid , THREAD_POOL *thread_pool )
642{
643 __n_assert( fluid , return FALSE );
644
645 //n_fluid_integrate( fluid );
646 list_foreach( node , fluid -> integrate_chunk_list )
647 {
648 add_threaded_process( thread_pool, &n_fluid_integrate_proc, (void *)node -> ptr , SYNCED_PROC);
649 }
650 start_threaded_pool( thread_pool );
651 wait_for_synced_threaded_pool( thread_pool );
652 refresh_thread_pool( thread_pool );
653
654 // set pressure to 0
655 memset( fluid -> p , 0 , fluid -> numCells * sizeof( double ) );
656
657 //n_fluid_solveIncompressibility( fluid );
658 for( size_t iter = 0 ; iter < fluid -> numIters ; iter++ )
659 {
660 list_foreach( node , fluid -> solveIncompressibility_chunk_list )
661 {
662 add_threaded_process( thread_pool, &n_fluid_solveIncompressibility_proc, (void *)node -> ptr , SYNCED_PROC);
663 }
664 start_threaded_pool( thread_pool );
665 wait_for_synced_threaded_pool( thread_pool );
666 refresh_thread_pool( thread_pool );
667 }
668
669 // extrapolate
670 n_fluid_extrapolate( fluid );
671
672 //n_fluid_advectVel( fluid );
673 memcpy( fluid -> newU , fluid -> u , fluid -> numCells * sizeof( double ) );
674 memcpy( fluid -> newV , fluid -> v , fluid -> numCells * sizeof( double ) );
675
676 list_foreach( node , fluid -> advectVel_chunk_list )
677 {
678 add_threaded_process( thread_pool, &n_fluid_advectVel_proc, (void *)node -> ptr , SYNCED_PROC);
679 }
680 start_threaded_pool( thread_pool );
681 wait_for_synced_threaded_pool( thread_pool );
682 refresh_thread_pool( thread_pool );
683
684 double *ptr = fluid -> u ;
685 fluid -> u = fluid -> newU ;
686 fluid -> newU = ptr ;
687
688 ptr = fluid -> v ;
689 fluid -> v = fluid -> newV ;
690 fluid -> newV = ptr ;
691
692
693 //n_fluid_advectSmoke( fluid );
694 memcpy( fluid -> newM , fluid -> m , fluid -> numCells * sizeof( double ) );
695 list_foreach( node , fluid -> advectSmoke_chunk_list )
696 {
697 add_threaded_process( thread_pool, &n_fluid_advectSmoke_proc, (void *)node -> ptr , SYNCED_PROC);
698 }
699 start_threaded_pool( thread_pool );
700 wait_for_synced_threaded_pool( thread_pool );
701 refresh_thread_pool( thread_pool );
702
703 ptr = fluid -> m ;
704 fluid -> m = fluid -> newM ;
705 fluid -> newM = ptr ;
706
707 return TRUE ;
708}
709
710
711
722int n_fluid_setObstacle( N_FLUID *fluid , double x , double y , double vx , double vy , double r )
723{
724 __n_assert( fluid , return FALSE );
725
726 size_t n = fluid -> numY;
727 for( size_t i = 1; i < fluid -> numX - 2 ; i++ )
728 {
729 for( size_t j = 1 ; j < fluid -> numY - 2 ; j++ )
730 {
731
732 double dx = (i + 0.5) - x;
733 double dy = (j + 0.5) - y;
734
735 if( i > 7 && ( dx * dx + dy * dy < r * r) )
736 {
737 fluid -> s[i*n + j] = 0.0;
738 fluid -> m[i*n + j] = 1.0;
739 fluid -> u[i*n + j] = vx;
740 fluid -> u[(i+1)*n + j] = vx;
741 fluid -> v[i*n + j] = vy;
742 fluid -> v[i*n + j+1] = vy;
743 }
744 }
745 }
746 return TRUE ;
747}
748
749
750
756int n_fluid_resetObstacles( N_FLUID *fluid )
757{
758 __n_assert( fluid , return FALSE );
759
760 size_t n = fluid -> numY;
761 for( size_t i = 1; i < fluid -> numX - 2 ; i++ )
762 {
763 for( size_t j = 1 ; j < fluid -> numY - 2 ; j++ )
764 {
765 fluid -> s[ i*n + j ] = 1.0 ;
766 }
767 }
768
769 return TRUE ;
770}
771
772
784int n_fluid_setObstacleFromBitmap( N_FLUID *fluid , ALLEGRO_BITMAP *bitmap , double x , double y , double vx , double vy , double r )
785{
786 __n_assert( fluid , return FALSE );
787
788 al_lock_bitmap( bitmap , al_get_bitmap_format( bitmap ) , ALLEGRO_LOCK_READONLY );
789
790 size_t n = fluid -> numY;
791 for( size_t i = 1; i < fluid -> numX - 2 ; i++ )
792 {
793 for( size_t j = 1 ; j < fluid -> numY - 2 ; j++ )
794 {
795
796 double dx = (i + 0.5) - x;
797 double dy = (j + 0.5) - y;
798
799 if( i > 7 && ( dx * dx + dy * dy < r * r) )
800 {
801 fluid -> s[i*n + j] = 0.0;
802 fluid -> m[i*n + j] = 1.0;
803 fluid -> u[i*n + j] = vx;
804 fluid -> v[i*n + j] = vy;
805
806 fluid -> u[(i+1)*n + j] = vx;
807 fluid -> v[i*n + j+1] = vy;
808 }
809 }
810 }
811
812 al_unlock_bitmap( bitmap );
813
814 return TRUE ;
815}
816
817
818
827ALLEGRO_COLOR n_fluid_getSciColor( N_FLUID *fluid , double val , double minVal , double maxVal )
828{
829 val = MIN( MAX( val , minVal ) , maxVal - 0.0001 );
830 double d = maxVal - minVal;
831 if( _zd( fluid , d ) )
832 {
833 val = 0.5 ;
834 }
835 else
836 {
837 val = (val - minVal) / d ;
838 }
839 double m = 0.25;
840 size_t num = floor( val / m );
841 double s = (val - num * m) / m;
842 double r = 0.0 , g = 0.0 , b = 0.0 ;
843 switch (num) {
844 case 0 : r = 0.0; g = s; b = 1.0; break;
845 case 1 : r = 0.0; g = 1.0; b = 1.0-s; break;
846 case 2 : r = s; g = 1.0; b = 0.0; break;
847 case 3 : r = 1.0; g = 1.0 - s; b = 0.0; break;
848 }
849 //return[255*r,255*g,255*b, 255]
850 return al_map_rgb_f( r , g , b );
851}
852
853
854
855
861int n_fluid_draw( N_FLUID *fluid )
862{
863 __n_assert( fluid , return FALSE );
864
865 size_t n = fluid -> numY ;
866
867 double minP = fluid -> p[0];
868 double maxP = fluid -> p[0];
869
870 if( fluid -> showPressure )
871 {
872 for( size_t i = 0; i < fluid -> numCells; i++)
873 {
874 minP = MIN( minP , fluid -> p[i] );
875 maxP = MAX( maxP , fluid -> p[i] );
876 }
877 }
878
879 ALLEGRO_COLOR color ;
880 double cScale = fluid -> cScale ;
881
882 for (size_t i = 0; i < fluid -> numX; i++)
883 {
884 for (size_t j = 0; j < fluid -> numY; j++)
885 {
886 int64_t x = i * cScale ;
887 int64_t y = j * cScale ;
888 int64_t cx = x + cScale ;
889 int64_t cy = y + cScale ;
890
891 double s = fluid -> m[i*n + j];
892
893 if( fluid -> showPaint )
894 {
895 color = n_fluid_getSciColor( fluid , s , 0.0 , 1.0 );
896 }
897 else if( fluid -> showPressure)
898 {
899 float color_vec_f[ 3 ] = { 0.0 , 0.0 , 0.0 };
900 double p = fluid -> p[i*n + j];
901 color = n_fluid_getSciColor( fluid , p , minP , maxP );
902 if( fluid -> showSmoke )
903 {
904 al_unmap_rgb_f( color , color_vec_f , color_vec_f + 1 , color_vec_f + 2 );
905 color = al_map_rgb_f( MAX( 0.0, color_vec_f[ 0 ] - s ) , MAX( 0.0 , color_vec_f[ 1 ] - s ) , MAX( 0.0, color_vec_f[ 2 ] - s) );
906 }
907 }
908 else if( fluid -> showSmoke )
909 {
910 color = al_map_rgb_f( 1.0 - s , 0.0 , 0.0 );
911 }
912 else
913 {
914 color = al_map_rgb_f( s , s , s );
915 }
916 al_draw_filled_rectangle( x , y , cx , cy , color );
917 }
918 }
919 return TRUE ;
920}
FreeNoLog
#define FreeNoLog(__ptr)
Free Handler without log.
Definition: n_common.h:268
Malloc
#define Malloc(__ptr, __struct, __size)
Malloc Handler to get errors and set to 0.
Definition: n_common.h:183
__n_assert
#define __n_assert(__ptr, __ret)
macro to assert things
Definition: n_common.h:276
list_push
int list_push(LIST *list, void *ptr, void(*destructor)(void *ptr))
Add a pointer to the end of the list.
Definition: n_list.c:244
list_foreach
#define list_foreach(__ITEM_, __LIST_)
ForEach macro helper.
Definition: n_list.h:70
new_generic_list
LIST * new_generic_list(int max_items)
Initialiaze a generic list container to max_items pointers.
Definition: n_list.c:20
N_FLUID::newU
double * newU
holder for newU arrays
Definition: n_fluids.h:95
N_FLUID::numY
size_t numY
number of cells in Y
Definition: n_fluids.h:58
N_FLUID::numX
size_t numX
number of cells in X
Definition: n_fluids.h:56
N_FLUID::u
double * u
holder for u arrays
Definition: n_fluids.h:93
N_FLUID::newV
double * newV
holder for newV arrays
Definition: n_fluids.h:100
N_FLUID::s
double * s
holder for s arrays
Definition: n_fluids.h:105
N_FLUID::v
double * v
holder for v arrays
Definition: n_fluids.h:98
N_FLUID::m
double * m
holder for m arrays
Definition: n_fluids.h:108
N_FLUID::newM
double * newM
holder for newM arrays
Definition: n_fluids.h:110
N_FLUID::p
double * p
holder for p arrays
Definition: n_fluids.h:103
n_fluid_avgV
double n_fluid_avgV(N_FLUID *fluid, size_t i, size_t j)
compute the average V value at a fluid position using it's surrounding
Definition: n_fluids.c:428
n_fluid_draw
int n_fluid_draw(N_FLUID *fluid)
draw a N_FLUID on screen / targert bitmap
Definition: n_fluids.c:861
n_fluid_simulate_threaded
int n_fluid_simulate_threaded(N_FLUID *fluid, THREAD_POOL *thread_pool)
a threaded version of N_FLUID global processing function
Definition: n_fluids.c:641
_z
#define _z(__fluid, __componant)
test if componant is near zero, according to fluid's precision
Definition: n_fluids.h:29
n_fluid_simulate
int n_fluid_simulate(N_FLUID *fluid)
non threaded version of N_FLUID global processing function
Definition: n_fluids.c:621
n_fluid_avgU
double n_fluid_avgU(N_FLUID *fluid, size_t i, size_t j)
compute the average U value at a fluid position using it's surrounding
Definition: n_fluids.c:408
n_fluid_sampleField
double n_fluid_sampleField(N_FLUID *fluid, double x, double y, uint32_t field)
compute a sample value at a field position
Definition: n_fluids.c:359
n_fluid_getSciColor
ALLEGRO_COLOR n_fluid_getSciColor(N_FLUID *fluid, double val, double minVal, double maxVal)
get fonky colors for the fluid
Definition: n_fluids.c:827
n_fluid_extrapolate
int n_fluid_extrapolate(N_FLUID *fluid)
non threaded extrapolation function
Definition: n_fluids.c:333
n_fluid_advectSmoke
int n_fluid_advectSmoke(N_FLUID *fluid)
non threaded version of add smoke function
Definition: n_fluids.c:582
n_fluid_solveIncompressibility
int n_fluid_solveIncompressibility(N_FLUID *fluid)
non threaded version of incompressibility solving function
Definition: n_fluids.c:289
n_fluid_resetObstacles
int n_fluid_resetObstacles(N_FLUID *fluid)
reset the obstacles set in a N_FLUID
Definition: n_fluids.c:756
_zd
#define _zd(__fluid, __value)
test if value is near zero, according to fluid's precision
Definition: n_fluids.h:34
n_fluid_integrate
int n_fluid_integrate(N_FLUID *fluid)
non threaded version of integration function
Definition: n_fluids.c:225
n_fluid_setObstacle
int n_fluid_setObstacle(N_FLUID *fluid, double x, double y, double vx, double vy, double r)
set an obstacle in the fluid grid
Definition: n_fluids.c:722
destroy_n_fluid
int destroy_n_fluid(N_FLUID **fluid)
ddestroy a fluid structure
Definition: n_fluids.c:48
n_fluid_advectVel
int n_fluid_advectVel(N_FLUID *fluid)
non threaded version of add velocities function
Definition: n_fluids.c:492
new_n_fluid
N_FLUID * new_n_fluid(double density, double gravity, size_t numIters, double dt, double overRelaxation, size_t sx, size_t sy)
!
Definition: n_fluids.c:78
N_FLUID
structure of a fluid
Definition: n_fluids.h:54
N_FLUID_THREAD_PARAMS
structure passed to a threaded fluid process
Definition: n_fluids.h:39
start_threaded_pool
int start_threaded_pool(THREAD_POOL *thread_pool)
Launch the process waiting for exectution in the thread pool.
Definition: n_thread_pool.c:286
get_nb_cpu_cores
int get_nb_cpu_cores()
get number of core of current system
Definition: n_thread_pool.c:27
SYNCED_PROC
#define SYNCED_PROC
processing mode for added func, synced start
Definition: n_thread_pool.h:27
add_threaded_process
int add_threaded_process(THREAD_POOL *thread_pool, void *(*func_ptr)(void *param), void *param, int mode)
add a function and params to a thread pool
Definition: n_thread_pool.c:205
refresh_thread_pool
int refresh_thread_pool(THREAD_POOL *thread_pool)
try to add some waiting DIRECT_PROCs on some free thread slots, else do nothing
Definition: n_thread_pool.c:469
wait_for_synced_threaded_pool
int wait_for_synced_threaded_pool(THREAD_POOL *thread_pool)
wait for all the launched process, blocking but light on the CPU as there is no polling
Definition: n_thread_pool.c:326
THREAD_POOL
Structure of a trhead pool.
Definition: n_thread_pool.h:81
n_common.h
Common headers and low-level hugly functions & define.
n_memset
void n_memset(void *dst, void *val, size_t size, size_t count)
memset bytes to a custom value
Definition: n_fluids.c:32
N_FLUID_S_FIELD
#define N_FLUID_S_FIELD
array number for s
Definition: n_fluids.c:22
N_FLUID_U_FIELD
#define N_FLUID_U_FIELD
array number for u
Definition: n_fluids.c:18
n_fluid_integrate_proc
void * n_fluid_integrate_proc(void *ptr)
ready to be threaded integration function
Definition: n_fluids.c:202
n_fluid_advectVel_proc
void * n_fluid_advectVel_proc(void *ptr)
ready to be threaded add velocities function
Definition: n_fluids.c:444
n_fluid_advectSmoke_proc
void * n_fluid_advectSmoke_proc(void *ptr)
ready to be threaded add smoke function
Definition: n_fluids.c:550
n_fluid_setObstacleFromBitmap
int n_fluid_setObstacleFromBitmap(N_FLUID *fluid, ALLEGRO_BITMAP *bitmap, double x, double y, double vx, double vy, double r)
set an obstacle in the fluid grid from a bitmap mask
Definition: n_fluids.c:784
N_FLUID_V_FIELD
#define N_FLUID_V_FIELD
array number for v
Definition: n_fluids.c:20
n_fluid_solveIncompressibility_proc
void * n_fluid_solveIncompressibility_proc(void *ptr)
ready to be threaded incompressibility solving function
Definition: n_fluids.c:248
n_fluids.h
fluid management port from "How to write an Eulerian fluid simulator with 200 lines of code",...
n_thread_pool.h
Thread pool declaration.