下面是用osg3.6.5可视化的烟雾模拟,smoke simulation.这里fluid solver来自" Jos Stam, Real-Time Fluid Dynamics for Games"，
按D键，添加烟雾，按G，T，H分别添加x,y,z方向的力，添加的烟雾过一阵自动会消散。
screenshot为
solver3d.h
#ifndef SOLVER3D_H #define SOLVER3D_H #ifdef __cplusplus extern "C"{ #endif void dens_step(int,int,int,float *,float *,float *,float *,float *,float,float); void vel_step(int,int,int,float *,float *,float *,float *,float *,float *,float,float); #ifdef __cplusplus } #endif #endif
solver3d.c
/*A 3D Real Time Fluid Solver based on Jos Stam's fluid solver (stable N-S solver). Reference: Jos Stam, "Real-Time Fluid Dynamics for Games". Proceedings of the Game Developer Conference, March 2003. *http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf] *The code is only for research. *By SoRoMan. */ #define IX(i,j,k) ((i)+(M+2)*(j) + (M+2)*(N+2)*(k)) /* 3d arrary access index*/ #define SWAP(x0,x) {float * tmp=x0;x0=x;x=tmp;} #define FOR_EACH_CELL for ( i=1 ; i<=M ; i++ ) { for ( j=1 ; j<=N ; j++ ) { for ( k=1 ; k<=O ; k++ ) { #define END_FOR }}} #define MAX(a,b) (((a) > (b)) ? (a) : (b)) #define LINEARSOLVERTIMES 10 void add_source ( int M, int N, int O, float * x, float * s, float dt ) { int i, size=(M+2)*(N+2)*(O+2); for ( i=0 ; i<size ; i++ ) x[i] += dt*s[i]; } void set_bnd ( int M, int N, int O, int b, float * x ) { int i,j; for ( i=1 ; i<=M ; i++ ) { for ( j=1 ; j<=N ; j++ ) { x[IX(i,j,0 )] = b==3 ? -x[IX(i,j,1)] : x[IX(i,j,1)]; x[IX(i,j,O+1)] = b==3 ? -x[IX(i,j,O)] : x[IX(i,j,O)]; } } for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=O ; j++ ) { x[IX(0 ,i, j)] = b==1 ? -x[IX(1,i,j)] : x[IX(1,i,j)]; x[IX(M+1,i, j)] = b==1 ? -x[IX(M,i,j)] : x[IX(M,i,j)]; } } for ( i=1 ; i<=M ; i++ ) { for ( j=1 ; j<=O ; j++ ) { x[IX(i,0,j )] = b==2 ? -x[IX(i,1,j)] : x[IX(i,1,j)]; x[IX(i,N+1,j)] = b==2 ? -x[IX(i,N,j)] : x[IX(i,N,j)]; } } x[IX(0 ,0, 0 )] = 1.0/3.0*(x[IX(1,0,0 )]+x[IX(0 ,1,0)]+x[IX(0 ,0,1)]); x[IX(0 ,N+1, 0)] = 1.0/3.0*(x[IX(1,N+1, 0)]+x[IX(0 ,N, 0)] + x[IX(0 ,N+1, 1)]); x[IX(M+1,0, 0 )] = 1.0/3.0*(x[IX(M,0,0 )]+x[IX(M+1,1,0)] + x[IX(M+1,0,1)]) ; x[IX(M+1,N+1,0)] = 1.0/3.0*(x[IX(M,N+1,0)]+x[IX(M+1,N,0)]+x[IX(M+1,N+1,1)]); x[IX(0 ,0, O+1 )] = 1.0/3.0*(x[IX(1,0,O+1 )]+x[IX(0 ,1,O+1)]+x[IX(0 ,0,O)]); x[IX(0 ,N+1, O+1)] = 1.0/3.0*(x[IX(1,N+1, O+1)]+x[IX(0 ,N, O+1)] + x[IX(0 ,N+1, O)]); x[IX(M+1,0, O+1 )] = 1.0/3.0*(x[IX(M,0,O+1 )]+x[IX(M+1,1,O+1)] + x[IX(M+1,0,O)]) ; x[IX(M+1,N+1,O+1)] = 1.0/3.0*(x[IX(M,N+1,O+1)]+x[IX(M+1,N,O+1)]+x[IX(M+1,N+1,O)]); } void lin_solve ( int M, int N, int O, int b, float * x, float * x0, float a, float c ) { int i, j, k,l; for ( l=0 ; l<LINEARSOLVERTIMES ; l++ ) { FOR_EACH_CELL x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c; END_FOR set_bnd ( M, N, O, b, x ); } } void diffuse ( int M, int N, int O, int b, float * x, float * x0, float diff, float dt ) { /*float a=dt*diff*M*N*O;*/ int max = MAX(MAX(M, N), MAX(N, O)); float a=dt*diff*max*max*max; lin_solve ( M, N, O, b, x, x0, a, 1+6*a ); } void advect ( int M, int N, int O, int b, float * d, float * d0, float * u, float * v, float * w, float dt ) { int i, j, k, i0, j0, k0, i1, j1, k1; float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz; dtx=dty=dtz=dt*MAX(MAX(M, N), MAX(N, O)); /*dtx = dt*M; dty = dt*N; dtz = dt*O;*/ FOR_EACH_CELL x = i-dtx*u[IX(i,j,k)]; y = j-dty*v[IX(i,j,k)]; z = k-dtz*w[IX(i,j,k)]; if (x<0.5f) x=0.5f; if (x>M+0.5f) x=M+0.5f; i0=(int)x; i1=i0+1; if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1; if (z<0.5f) z=0.5f; if (z>O+0.5f) z=O+0.5f; k0=(int)z; k1=k0+1; s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; u1 = z-k0; u0 = 1-u1; d[IX(i,j,k)] = s0*(t0*u0*d0[IX(i0,j0,k0)]+t1*u0*d0[IX(i0,j1,k0)]+t0*u1*d0[IX(i0,j0,k1)]+t1*u1*d0[IX(i0,j1,k1)])+ s1*(t0*u0*d0[IX(i1,j0,k0)]+t1*u0*d0[IX(i1,j1,k0)]+t0*u1*d0[IX(i1,j0,k1)]+t1*u1*d0[IX(i1,j1,k1)]); END_FOR set_bnd (M, N, O, b, d ); } void project ( int M, int N, int O, float * u, float * v, float * w, float * p, float * div ) { int i, j, k; FOR_EACH_CELL div[IX(i,j,k)] = -1.0/3.0*((u[IX(i+1,j,k)]-u[IX(i-1,j,k)])/M+(v[IX(i,j+1,k)]-v[IX(i,j-1,k)])/M+(w[IX(i,j,k+1)]-w[IX(i,j,k-1)])/M); p[IX(i,j,k)] = 0; END_FOR set_bnd ( M, N, O, 0, div ); set_bnd (M, N, O, 0, p ); lin_solve ( M, N, O, 0, p, div, 1, 6 ); FOR_EACH_CELL u[IX(i,j,k)] -= 0.5f*M*(p[IX(i+1,j,k)]-p[IX(i-1,j,k)]); v[IX(i,j,k)] -= 0.5f*M*(p[IX(i,j+1,k)]-p[IX(i,j-1,k)]); w[IX(i,j,k)] -= 0.5f*M*(p[IX(i,j,k+1)]-p[IX(i,j,k-1)]); END_FOR set_bnd ( M, N, O, 1, u ); set_bnd ( M, N, O, 2, v );set_bnd ( M, N, O, 3, w); } void dens_step ( int M, int N, int O, float * x, float * x0, float * u, float * v, float * w, float diff, float dt ) { add_source ( M, N, O, x, x0, dt ); SWAP ( x0, x ); diffuse ( M, N, O, 0, x, x0, diff, dt ); SWAP ( x0, x ); advect ( M, N, O, 0, x, x0, u, v, w, dt ); } void vel_step ( int M, int N, int O, float * u, float * v, float * w, float * u0, float * v0, float * w0, float visc, float dt ) { add_source ( M, N, O, u, u0, dt ); add_source ( M, N, O, v, v0, dt );add_source ( M, N, O, w, w0, dt ); SWAP ( u0, u ); diffuse ( M, N, O, 1, u, u0, visc, dt ); SWAP ( v0, v ); diffuse ( M, N, O, 2, v, v0, visc, dt ); SWAP ( w0, w ); diffuse ( M, N, O, 3, w, w0, visc, dt ); project ( M, N, O, u, v, w, u0, v0 ); SWAP ( u0, u ); SWAP ( v0, v );SWAP ( w0, w ); advect ( M, N, O, 1, u, u0, u0, v0, w0, dt ); advect ( M, N, O, 2, v, v0, u0, v0, w0, dt );advect ( M, N, O, 3, w, w0, u0, v0, w0, dt ); project ( M, N, O, u, v, w, u0, v0 ); }
main.cpp
// OpenSceneGraph library. #include <osgDB/ReadFile> #include <osgViewer/Viewer> #include <osgGA/StateSetManipulator> #include <osgGA/TrackballManipulator> #include <osgViewer/ViewerEventHandlers> #include <osgUtil/SmoothingVisitor> #include <osg/BlendFunc> #include <osg/PositionAttitudeTransform> #include <stdio.h> #include "solver3d.h" #include <time.h> //#pragma comment(lib, "osg.lib") //#pragma comment(lib, "osgDB.lib") //#pragma comment(lib, "osgGA.lib") //#pragma comment(lib, "osgViewer.lib") //#pragma comment(lib, "osgUtil.lib") typedef unsigned char boolean; #define FALSE 0 #define TRUE 1 /* global variables */ /*fluid field information*/ static int M = 30; /*grid size x*/ static int N = 30; /*grid size y*/ static int O = 30; /*grid size z*/ static float dt = 0.4f; /*time step*/ static float diff = 0.0f; /*diffuse factor*/ static float visc = 0.0f; /*viscous factor*/ static float force = 10.0f; /*force added each time*/ static float source = 200.0f; /*source added each time*/ static int vert_sum = M*N*O; //#define VID(i,j,k) ((i)+(M)*(j) + (M)*(N)*(k)) //#define VNEXT_I(id) (id+1) //#define VNEXT_J(id) (id+M) //#define VNEXT_K(id) (id+M*N) #define VID(i,j,k) ((i)+(M+2)*(j) + (M+2)*(N+2)*(k)) #define VNEXT_I(id) (id+1) #define VNEXT_J(id) (id+M+2) #define VNEXT_K(id) (id+(M+2)*(N+2)) #define IX(i,j,k) ((i)+(M+2)*(j) + (M+2)*(N+2)*(k)) #define MAX(a,b) (((a) > (b)) ? (a) : (b)) static boolean addforce[3] = {FALSE, FALSE, FALSE}; static boolean addsource = FALSE; static float * u, * v, *w, * u_prev, * v_prev, * w_prev; static float * dens, * dens_prev; osg::ref_ptr<osg::Geometry> s_volmesh; static void free_data ( void ) { if ( u ) free ( u ); if ( v ) free ( v ); if ( w ) free ( w ); if ( u_prev ) free ( u_prev ); if ( v_prev ) free ( v_prev ); if ( w_prev ) free ( w_prev ); if ( dens ) free ( dens ); if ( dens_prev ) free ( dens_prev ); } static void clear_data ( void ) { int i, size=(M+2)*(N+2)*(O+2); for ( i=0 ; i<size ; i++ ) { u[i] = v[i] = w[i] = u_prev[i] = v_prev[i] =w_prev[i] = dens[i] = dens_prev[i] = 0.0f; } addforce[0] = addforce[1] = addforce[2] = FALSE; } static int allocate_data ( void ) { int size = (M+2)*(N+2)*(O+2); u = (float *) malloc ( size*sizeof(float) ); v = (float *) malloc ( size*sizeof(float) ); w = (float *) malloc ( size*sizeof(float) ); u_prev = (float *) malloc ( size*sizeof(float) ); v_prev = (float *) malloc ( size*sizeof(float) ); w_prev = (float *) malloc ( size*sizeof(float) ); dens = (float *) malloc ( size*sizeof(float) ); dens_prev = (float *) malloc ( size*sizeof(float) ); if ( !u || !v || !w || !u_prev || !v_prev || !w_prev || !dens || !dens_prev ) { fprintf ( stderr, "cannot allocate data\n" ); return ( 0 ); } return ( 1 ); } osg::Geode* buildVolumeMesh(int m, int n, int o, float cellh) { osg::Geode* geod = new osg::Geode; s_volmesh = new osg::Geometry; geod->addDrawable(s_volmesh); s_volmesh->setUseDisplayList(false); s_volmesh->setUseVertexBufferObjects(true); osg::ref_ptr<osg::Vec3Array> pointsVec = new osg::Vec3Array(); osg::ref_ptr<osg::Vec4Array> colorVec = new osg::Vec4Array(); pointsVec->reserve(o*n*m); colorVec->assign(o*n*m, osg::Vec4(0.0,0.0,0.0,0.05)); for (int k = 0; k < o; k++) { for (int j = 0; j < n; j++) { for (int i = 0; i < m; i++) { pointsVec->push_back(osg::Vec3(i*cellh,j*cellh,k*cellh)); } } } s_volmesh->setVertexArray(pointsVec); s_volmesh->setColorArray(colorVec, osg::Array::BIND_PER_VERTEX); std::vector<unsigned int> indices; for (int k = 0; k < o; k++) { for (int j = 0; j < (n-1); j++) { for (int i = 0; i < (m-1); i++) { //indices.push_back(); // 3 2 // 0 1 unsigned int id0 = VID(i,j,k); unsigned int id1 = VNEXT_I(id0); unsigned int id3 = VNEXT_J(id0); unsigned int id2 = VNEXT_I(id3); indices.push_back(id0); indices.push_back(id1); indices.push_back(id2); indices.push_back(id0); indices.push_back(id2); indices.push_back(id3); } } } for (int i = 0; i < m; i++) { for (int k = 0; k < (o-1); k++) { for (int j = 0; j < (n-1); j++) { unsigned int id0 = VID(i,j,k); unsigned int id1 = VNEXT_J(id0); unsigned int id3 = VNEXT_K(id0); unsigned int id2 = VNEXT_J(id3); indices.push_back(id0); indices.push_back(id1); indices.push_back(id2); indices.push_back(id0); indices.push_back(id2); indices.push_back(id3); } } } for (int j = 0; j < n; j++) { for (int i = 0; i < (m-1); i++) { for (int k = 0; k < (o-1); k++) { unsigned int id0 = VID(i,j,k); unsigned int id1 = VNEXT_K(id0); unsigned int id3 = VNEXT_I(id0); unsigned int id2 = VNEXT_K(id3); indices.push_back(id0); indices.push_back(id1); indices.push_back(id2); indices.push_back(id0); indices.push_back(id2); indices.push_back(id3); } } } s_volmesh->addPrimitiveSet(new osg::DrawElementsUInt(osg::PrimitiveSet::TRIANGLES, indices.size(), (unsigned int*)&indices[0])); geod->getOrCreateStateSet()->setMode(GL_LIGHTING, osg::StateAttribute::OFF | osg::StateAttribute::OVERRIDE); geod->getOrCreateStateSet()->setMode(GL_CULL_FACE, osg::StateAttribute::OFF | osg::StateAttribute::OVERRIDE); geod->getOrCreateStateSet()->setMode(GL_DEPTH_TEST, osg::StateAttribute::OFF | osg::StateAttribute::OVERRIDE); geod->getOrCreateStateSet()->setMode(GL_BLEND, osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE); geod->getOrCreateStateSet()->setAttribute(new osg::BlendFunc(GL_SRC_ALPHA, GL_DST_ALPHA),osg::StateAttribute::ON|osg::StateAttribute::OVERRIDE); return geod; } static void get_force_source ( float * d, float * u, float * v, float * w ) { int i, j, k, size=(M+2)*(N+2)*(O+2);; for ( i=0 ; i<size ; i++ ) { u[i] = v[i] = w[i]= d[i] = 0.0f; } if(addforce[0]) /*x axis*/ { i=2, j=N/2; k=O/2; if ( i<1 || i>M || j<1 || j>N || k <1 || k>O) return; u[IX(i,j,k)] = force*10; addforce[0] = FALSE; } if(addforce[1]) { i=M/2, j=2; k=O/2; if ( i<1 || i>M || j<1 || j>N || k <1 || k>O) return; v[IX(i,j,k)] = force*10; addforce[1] = FALSE; } if(addforce[2]) /* z-axis*/ { i=M/2, j=N/2; k=2; if ( i<1 || i>M || j<1 || j>N || k <1 || k>O) return; w[IX(i,j,k)] = force*10; addforce[2] = FALSE; } if(addsource) { i=M/2, j=N/2; k=O/2; d[IX(i,j,k)] = source; addsource = FALSE; } return; } int Game_Init(void) { /*init grid size, temp step, diffuse factor, viscous factor, force and source added each time*/ /*M = 16; N = 16; O= 16; dt = 0.4f; diff = 0.0f; visc = 0.0f; force = 10.0f; source = 200.0f;*/ /*init rot*/ if ( !allocate_data () ) return (0); clear_data (); //init_dens_offs(); return (1); } void display() { osg::Vec4Array* cols = dynamic_cast<osg::Vec4Array*>(s_volmesh->getColorArray()); for (int i = 0; i < (int)cols->size(); i++) { float d = dens[i]; cols->at(i)._v[0] = d; cols->at(i)._v[1] = d; cols->at(i)._v[2] = d; if (d >0.f) { d *= 0.999f; } else { d = 0.f; } dens[i] = d; } cols->dirty(); } int Game_Main(void) { /* compute fluids*/ get_force_source( dens_prev, u_prev, v_prev, w_prev ); vel_step ( M,N,O, u, v, w, u_prev, v_prev,w_prev, visc, dt ); dens_step ( M,N,O, dens, dens_prev, u, v, w, diff, dt ); /* render fluids*/ display(); return (1); } int Game_Shutdown(void) { /* this function is where you shutdown your game and release all resources that you allocated*/ free_data (); /* return success*/ return(1); } void Game_Reset() { clear_data(); } class KeyMouseHandler : public osgGA::GUIEventHandler { public: KeyMouseHandler(){}; ~KeyMouseHandler(){}; bool handle(const osgGA::GUIEventAdapter& ea,osgGA::GUIActionAdapter& aa) { if (osgGA::GUIEventAdapter::KEYDOWN == ea.getEventType()) { if (osgGA::GUIEventAdapter::KEY_D == ea.getKey()) { addsource = TRUE; } else if (osgGA::GUIEventAdapter::KEY_G == ea.getKey()) { addforce[1] = TRUE; } else if (osgGA::GUIEventAdapter::KEY_H == ea.getKey()) { addforce[0] = TRUE; } else if (osgGA::GUIEventAdapter::KEY_T == ea.getKey()) { addforce[2] = TRUE; } else if (osgGA::GUIEventAdapter::KEY_C == ea.getKey()) { Game_Reset(); } } else if (osgGA::GUIEventAdapter::FRAME == ea.getEventType()) { Game_Main(); } return false; } }; float GetMidHeight(int i, int stride, std::vector<osg::Vec3>& points) { return 0.5f * (points[i - stride].y() + points[i + stride].y()); } int CalculateArrayLenght(int iterateTime) { unsigned long base = 0x1; unsigned int c = base << iterateTime; return c+1; } float getRand(float lowb, float upb) { float t = rand() / double(RAND_MAX); return lowb*(1.f-t)+upb*t; } std::vector<osg::Vec3> Generater1DLineIterative(osg::Vec3 start, osg::Vec3 end, int iterateTime, float heightScale, float h) { srand((unsigned)time(0)); int length = CalculateArrayLenght(iterateTime); std::vector<osg::Vec3> points; points.assign(length, osg::Vec3(0,0,0)); points[0] = start; points[length - 1] = end; float gap = fabs(end.x() - start.x()) / length; //Debug.Log("gap: " + gap); for(int i=0; i< length; i++) { points[i] = osg::Vec3(start.x() + i*gap, 0.f, 0.f); } float ratio, scale; ratio = (float)pow(2.0f, -h); scale = heightScale;// * ratio; int stride = length / 2; while(stride != 0) { for (int i = stride; i < length; i += stride) { float rd = getRand(-h, h); points[i].y() = scale * rd + GetMidHeight(i, stride, points); i += stride; } scale *= ratio; //Debug.Log("Stride: " + stride); stride >>= 1; } // DumpAllPoint(points); return points; } void Generate1DLineIter(std::vector<osg::Vec3>& points, osg::Vec3 start, osg::Vec3 end, int iterateTime, float roughless) { if (iterateTime > 0) { float rand = getRand(-2.999f, 2.999f) * roughless; osg::Vec3 mid = osg::Vec3(0.5f * start.x() + 0.5f * end.x(), 0.5f * start.y() + 0.5f * end.y() + rand, 0); --iterateTime; Generate1DLineIter(points, start, mid, iterateTime, roughless*0.5/* * iterateTime * iterateTime*/); points.push_back(mid); Generate1DLineIter(points, mid, end, iterateTime, roughless*0.5/* * iterateTime * iterateTime*/); } } std::vector<osg::Vec3> Generate1DLine(osg::Vec3 start, osg::Vec3 end, int iterateTime, float roughless) { srand((unsigned)time(0)); std::vector<osg::Vec3> pts; pts.push_back(start); Generate1DLineIter(pts, start, end, iterateTime, roughless); pts.push_back(end); return pts; } osg::Geode* testRandStrip() { //std::vector<osg::Vec3> pts = Generater1DLineIterative(osg::Vec3(0,0,0),osg::Vec3(10,0,0),3,2.0,5); std::vector<osg::Vec3> pts = Generate1DLine(osg::Vec3(0,0,0),osg::Vec3(10,0,0),12,1.5); osg::Geode* geod = new osg::Geode; osg::Geometry* geo = new osg::Geometry; geod->addDrawable(geo); geo->setUseDisplayList(false); geo->setUseVertexBufferObjects(true); geo->setVertexArray(new osg::Vec3Array(pts.begin(),pts.end())); osg::Vec4Array* colorVec = new osg::Vec4Array; colorVec->push_back(osg::Vec4(1,1,1,1)); geo->setColorArray(colorVec, osg::Array::BIND_OVERALL); std::vector<unsigned int> indices; for (int i = 1; i < pts.size(); i++) { indices.push_back(i-1); indices.push_back(i); } geo->addPrimitiveSet(new osg::DrawElementsUInt(osg::PrimitiveSet::LINES, indices.size(), (unsigned int*)&indices[0])); return geod; } int main(void) { osgViewer::Viewer myViewer; //myViewer.setSceneData(BuildScene()); osg::Group* root = new osg::Group; float h = 1.0f/MAX(MAX((M+2), (N+2)), MAX((N+2), (O+2))); root->addChild(buildVolumeMesh(M+2,N+2,O+2,h)); //root->addChild(testRandStrip()); myViewer.setSceneData(root); myViewer.addEventHandler(new KeyMouseHandler); myViewer.addEventHandler(new osgGA::StateSetManipulator(myViewer.getCamera()->getOrCreateStateSet())); myViewer.addEventHandler(new osgViewer::StatsHandler); myViewer.addEventHandler(new osgViewer::WindowSizeHandler); myViewer.setCameraManipulator(new osgGA::TrackballManipulator); Game_Init(); myViewer.realize(); while (!myViewer.done()) { Game_Main(); myViewer.frame(); } //return myViewer.run(); return 0; }