#
# ***** END GPL LICENSE BLOCK *****
-SET(INC ${PNG_INC} ${ZLIB_INC} ${SDL_INC})
+SET(INC ${PNG_INC} ${ZLIB_INC} ${SDL_INC} extern)
FILE(GLOB SRC intern/*.cpp)
-ADD_DEFINITIONS(-DNOGUI -DELBEEM_BLENDER=1)
+ADD_DEFINITIONS(-DNOGUI -DELBEEM_BLENDER=1 -DLBM_INCLUDE_CONTROL=1)
IF(WINDOWS)
ADD_DEFINITIONS(-DUSE_MSVC6FIXES)
ENDIF(WINDOWS)
sources = env.Glob('intern/*.cpp')
-defs = ' NOGUI ELBEEM_BLENDER=1'
+defs = 'NOGUI ELBEEM_BLENDER=1 LBM_INCLUDE_CONTROL=1'
if env['WITH_BF_OPENMP'] == 1:
defs += ' PARALLEL'
if env['OURPLATFORM']=='win32-vc':
defs += ' USE_MSVC6FIXES'
-incs = env['BF_PNG_INC'] + ' ' + env['BF_ZLIB_INC'] + ' ' +env['BF_SDL_INC']
+incs = env['BF_PNG_INC'] + ' ' + env['BF_ZLIB_INC'] + ' ' +env['BF_SDL_INC']
+incs += ' extern '
env.BlenderLib ('bf_elbeem', sources, Split(incs), Split(defs), libtype='blender', priority=0 )
short version;
/* id number of simulation domain, needed if more than a
* single domain should be simulated */
- short domainId;
+ short domainId; // unused within blender
/* geometrical extent */
float geoStart[3], geoSize[3];
// defines for elbeemMesh->type below
+/* please keep in sync with DNA_object_fluidsim.h */
#define OB_FLUIDSIM_FLUID 4
#define OB_FLUIDSIM_OBSTACLE 8
#define OB_FLUIDSIM_INFLOW 16
#define OB_FLUIDSIM_OUTFLOW 32
+#define OB_FLUIDSIM_PARTICLE 64
+#define OB_FLUIDSIM_CONTROL 128
// defines for elbeemMesh->obstacleType below
#define FLUIDSIM_OBSTACLE_NOSLIP 1
// a single mesh object
typedef struct elbeemMesh {
- /* obstacle,fluid or inflow... */
+ /* obstacle,fluid or inflow or control ... */
short type;
/* id of simulation domain it belongs to */
short parentDomainId;
/* name of the mesh, mostly for debugging */
const char *name;
+
+ /* fluid control settings */
+ float cpsTimeStart;
+ float cpsTimeEnd;
+ float cpsQuality;
+
+ int channelSizeAttractforceStrength;
+ float *channelAttractforceStrength;
+ int channelSizeAttractforceRadius;
+ float *channelAttractforceRadius;
+ int channelSizeVelocityforceStrength;
+ float *channelVelocityforceStrength;
+ int channelSizeVelocityforceRadius;
+ float *channelVelocityforceRadius;
} elbeemMesh;
// API functions
// start fluidsim init (returns !=0 upon failure)
int elbeemInit(void);
+// frees fluidsim
+int elbeemFree(void);
+
// start fluidsim init (returns !=0 upon failure)
int elbeemAddDomain(struct elbeemSimulationSettings*);
// structs, for these use OB_xxx defines above
/*! fluid geometry init types */
+// type "int" used, so max is 8
#define FGI_FLAGSTART 16
#define FGI_FLUID (1<<(FGI_FLAGSTART+ 0))
#define FGI_NO_FLUID (1<<(FGI_FLAGSTART+ 1))
#define FGI_NO_BND (1<<(FGI_FLAGSTART+ 5))
#define FGI_MBNDINFLOW (1<<(FGI_FLAGSTART+ 6))
#define FGI_MBNDOUTFLOW (1<<(FGI_FLAGSTART+ 7))
+#define FGI_CONTROL (1<<(FGI_FLAGSTART+ 8))
// all boundary types at once
#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
--- /dev/null
+// --------------------------------------------------------------------------
+//
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the
+// GNU General Public License. See the file COPYING for details.
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+// implementation of control particle handling
+//
+// --------------------------------------------------------------------------
+
+// indicator for LBM inclusion
+#include "ntl_geometrymodel.h"
+#include "ntl_world.h"
+#include "solver_class.h"
+#include "controlparticles.h"
+#include "mvmcoords.h"
+#include <zlib.h>
+
+#ifndef sqrtf
+#define sqrtf sqrt
+#endif
+
+// brute force circle test init in initTimeArray
+// replaced by mDebugInit
+//#define CP_FORCECIRCLEINIT 0
+
+
+void ControlParticles::initBlenderTest() {
+ mPartSets.clear();
+
+ ControlParticleSet cps;
+ mPartSets.push_back(cps);
+ int setCnt = mPartSets.size()-1;
+ ControlParticle p;
+
+ // set for time zero
+ mPartSets[setCnt].time = 0.;
+
+ // add single particle
+ p.reset();
+ p.pos = LbmVec(0.5, 0.5, -0.5);
+ mPartSets[setCnt].particles.push_back(p);
+
+ // add second set for animation
+ mPartSets.push_back(cps);
+ setCnt = mPartSets.size()-1;
+ mPartSets[setCnt].time = 0.15;
+
+ // insert new position
+ p.reset();
+ p.pos = LbmVec(-0.5, -0.5, 0.5);
+ mPartSets[setCnt].particles.push_back(p);
+
+ // applyTrafos();
+ initTime(0. , 1.);
+}
+
+// blender control object gets converted to mvm flui control object
+int ControlParticles::initFromObject(ntlGeometryObjModel *model) {
+ vector<ntlTriangle> triangles;
+ vector<ntlVec3Gfx> vertices;
+ vector<ntlVec3Gfx> normals;
+
+ /*
+ model->loadBobjModel(string(infile));
+
+ model->setLoaded(true);
+
+ model->setGeoInitId(gid);
+
+
+ printf("a animated? %d\n", model->getIsAnimated());
+ printf("b animated? %d\n", model->getMeshAnimated());
+ */
+
+ model->setGeoInitType(FGI_FLUID);
+
+ model->getTriangles(mCPSTimeStart, &triangles, &vertices, &normals, 1 );
+ // model->applyTransformation(mCPSTimeStart, &vertices, &normals, 0, vertices.size(), true);
+
+ // valid mesh?
+ if(triangles.size() <= 0) {
+ return 0;
+ }
+
+ ntlRenderGlobals *glob = new ntlRenderGlobals;
+ ntlScene *genscene = new ntlScene( glob, false );
+ genscene->addGeoClass(model);
+ genscene->addGeoObject(model);
+ genscene->buildScene(0., false);
+ char treeFlag = (1<<(4+model->getGeoInitId()));
+
+ ntlTree *tree = new ntlTree(
+ 15, 8, // TREEwarning - fixed values for depth & maxtriangles here...
+ genscene, treeFlag );
+
+ // TODO? use params
+ ntlVec3Gfx start,end;
+ model->getExtends(start,end);
+ /*
+ printf("start - x: %f, y: %f, z: %f\n", start[0], start[1], start[2]);
+ printf("end - x: %f, y: %f, z: %f\n", end[0], end[1], end[2]);
+ printf("mCPSWidth: %f\n");
+*/
+ LbmFloat width = mCPSWidth;
+ if(width<=LBM_EPSILON) { errMsg("ControlParticles::initFromMVMCMesh","Invalid mCPSWidth! "<<mCPSWidth); width=mCPSWidth=0.1; }
+ ntlVec3Gfx org = start+ntlVec3Gfx(width*0.5);
+ gfxReal distance = -1.;
+ vector<ntlVec3Gfx> inspos;
+ int approxmax = (int)( ((end[0]-start[0])/width)*((end[1]-start[1])/width)*((end[2]-start[2])/width) );
+
+ // printf("distance: %f, width: %f\n", distance, width);
+
+ while(org[2]<end[2]) {
+ while(org[1]<end[1]) {
+ while(org[0]<end[0]) {
+ if(checkPointInside(tree, org, distance)) {
+ inspos.push_back(org);
+ }
+ // TODO optimize, use distance
+ org[0] += width;
+ }
+ org[1] += width;
+ org[0] = start[0];
+ }
+ org[2] += width;
+ org[1] = start[1];
+ }
+
+ // printf("inspos.size(): %d\n", inspos.size());
+
+ MeanValueMeshCoords mvm;
+ mvm.calculateMVMCs(vertices,triangles, inspos, mCPSWeightFac);
+ vector<ntlVec3Gfx> ninspos;
+ mvm.transfer(vertices, ninspos);
+
+ // init first set, check dist
+ ControlParticleSet firstcps; //T
+ mPartSets.push_back(firstcps);
+ mPartSets[mPartSets.size()-1].time = mCPSTimeStart;
+ vector<bool> useCP;
+
+ for(int i=0; i<(int)inspos.size(); i++) {
+ ControlParticle p; p.reset();
+ p.pos = vec2L(inspos[i]);
+
+ double cpdist = norm(inspos[i]-ninspos[i]);
+ bool usecpv = true;
+
+ mPartSets[mPartSets.size()-1].particles.push_back(p);
+ useCP.push_back(usecpv);
+ }
+
+ // init further sets, temporal mesh sampling
+ double tsampling = mCPSTimestep;
+ // printf("tsampling: %f, ninspos.size(): %d, mCPSTimeEnd: %f\n", tsampling, ninspos.size(), mCPSTimeEnd);
+
+ int totcnt = (int)( (mCPSTimeEnd-mCPSTimeStart)/tsampling ), tcnt=0;
+ for(double t=mCPSTimeStart+tsampling; ((t<mCPSTimeEnd) && (ninspos.size()>0.)); t+=tsampling) {
+ ControlParticleSet nextcps; //T
+ mPartSets.push_back(nextcps);
+ mPartSets[mPartSets.size()-1].time = (gfxReal)t;
+
+ vertices.clear(); triangles.clear(); normals.clear();
+ model->getTriangles(t, &triangles, &vertices, &normals, 1 );
+ mvm.transfer(vertices, ninspos);
+
+ tcnt++;
+ for(size_t i=0; i < ninspos.size(); i++) {
+
+ if(useCP[i]) {
+ ControlParticle p; p.reset();
+ p.pos = vec2L(ninspos[i]);
+ mPartSets[mPartSets.size()-1].particles.push_back(p);
+ }
+ }
+ }
+
+ model->setGeoInitType(FGI_CONTROL);
+
+ delete tree;
+ delete genscene;
+ delete glob;
+
+ // do reverse here
+ if(model->getGeoPartSlipValue())
+ {
+ mirrorTime();
+ }
+
+ return 1;
+}
+
+
+// init all zero / defaults for a single particle
+void ControlParticle::reset() {
+ pos = LbmVec(0.,0.,0.);
+ vel = LbmVec(0.,0.,0.);
+ influence = 1.;
+ size = 1.;
+#ifndef LBMDIM
+#ifdef MAIN_2D
+ rotaxis = LbmVec(0.,1.,0.); // SPH xz
+#else // MAIN_2D
+ // 3d - roate in xy plane, vortex
+ rotaxis = LbmVec(0.,0.,1.);
+ // 3d - rotate for wave
+ //rotaxis = LbmVec(0.,1.,0.);
+#endif // MAIN_2D
+#else // LBMDIM
+ rotaxis = LbmVec(0.,1.,0.); // LBM xy , is swapped afterwards
+#endif // LBMDIM
+
+ density = 0.;
+ densityWeight = 0.;
+ avgVelAcc = avgVel = LbmVec(0.);
+ avgVelWeight = 0.;
+}
+
+
+// default preset/empty init
+ControlParticles::ControlParticles() :
+ _influenceTangential(0.f),
+ _influenceAttraction(0.f),
+ _influenceVelocity(0.f),
+ _influenceMaxdist(0.f),
+ _radiusAtt(1.0f),
+ _radiusVel(1.0f),
+ _radiusMinMaxd(2.0f),
+ _radiusMaxd(3.0f),
+ _currTime(-1.0), _currTimestep(1.),
+ _initTimeScale(1.),
+ _initPartOffset(0.), _initPartScale(1.),
+ _initLastPartOffset(0.), _initLastPartScale(1.),
+ _initMirror(""),
+ _fluidSpacing(1.), _kernelWeight(-1.),
+ _charLength(1.), _charLengthInv(1.),
+ mvCPSStart(-10000.), mvCPSEnd(10000.),
+ mCPSWidth(0.1), mCPSTimestep(0.02), // was 0.05
+ mCPSTimeStart(0.), mCPSTimeEnd(0.5), mCPSWeightFac(1.),
+ mDebugInit(0)
+{
+ _radiusAtt = 0.15f;
+ _radiusVel = 0.15f;
+ _radiusMinMaxd = 0.16f;
+ _radiusMaxd = 0.3;
+
+ _influenceAttraction = 0.f;
+ _influenceTangential = 0.f;
+ _influenceVelocity = 0.f;
+ // 3d tests */
+}
+
+
+
+ControlParticles::~ControlParticles() {
+ // nothing to do...
+}
+
+LbmFloat ControlParticles::getControlTimStart() {
+ if(mPartSets.size()>0) { return mPartSets[0].time; }
+ return -1000.;
+}
+LbmFloat ControlParticles::getControlTimEnd() {
+ if(mPartSets.size()>0) { return mPartSets[mPartSets.size()-1].time; }
+ return -1000.;
+}
+
+// calculate for delta t
+void ControlParticles::setInfluenceVelocity(LbmFloat set, LbmFloat dt) {
+ const LbmFloat dtInter = 0.01;
+ LbmFloat facFv = 1.-set; //cparts->getInfluenceVelocity();
+ // mLevel[mMaxRefine].timestep
+ LbmFloat facNv = (LbmFloat)( 1.-pow( (double)facFv, (double)(dt/dtInter)) );
+ //errMsg("vwcalc","ts:"<<dt<< " its:"<<(dt/dtInter) <<" fv"<<facFv<<" nv"<<facNv<<" test:"<< pow( (double)(1.-facNv),(double)(dtInter/dt)) );
+ _influenceVelocity = facNv;
+}
+
+int ControlParticles::initExampleSet()
+{
+ // unused
+ return 0;
+}
+
+int ControlParticles::getTotalSize()
+{
+ int s=0;
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ s+= mPartSets[i].particles.size();
+ }
+ return s;
+}
+
+// --------------------------------------------------------------------------
+// load positions & timing from text file
+// WARNING - make sure file has unix format, no win/dos linefeeds...
+#define LINE_LEN 100
+int ControlParticles::initFromTextFile(string filename)
+{
+ const bool debugRead = false;
+ char line[LINE_LEN];
+ line[LINE_LEN-1] = '\0';
+ mPartSets.clear();
+ if(filename.size()<2) return 0;
+
+ // HACK , use "cparts" suffix as old
+ // e.g. "cpart2" as new
+ if(filename[ filename.size()-1 ]=='s') {
+ return initFromTextFileOld(filename);
+ }
+
+ FILE *infile = fopen(filename.c_str(), "r");
+ if(!infile) {
+ errMsg("ControlParticles::initFromTextFile","unable to open '"<<filename<<"' " );
+ // try to open as gz sequence
+ if(initFromBinaryFile(filename)) { return 1; }
+ // try mesh MVCM generation
+ if(initFromMVCMesh(filename)) { return 1; }
+ // failed...
+ return 0;
+ }
+
+ int haveNo = false;
+ int haveScale = false;
+ int haveTime = false;
+ int noParts = -1;
+ int partCnt = 0;
+ int setCnt = 0;
+ //ControlParticle p; p.reset();
+ // scale times by constant factor while reading
+ LbmFloat timeScale= 1.0;
+ int lineCnt = 0;
+ bool abortParse = false;
+#define LASTCP mPartSets[setCnt].particles[ mPartSets[setCnt].particles.size()-1 ]
+
+ while( (!feof(infile)) && (!abortParse)) {
+ lineCnt++;
+ fgets(line, LINE_LEN, infile);
+
+ //if(debugRead) printf("\nDEBUG%d r '%s'\n",lineCnt, line);
+ if(!line) continue;
+ int len = (int)strlen(line);
+
+ // skip empty lines and comments (#,//)
+ if(len<1) continue;
+ if( (line[0]=='#') || (line[0]=='\n') ) continue;
+ if((len>1) && (line[0]=='/' && line[1]=='/')) continue;
+
+ // debug remove newline
+ if((len>=1)&&(line[len-1]=='\n')) line[len-1]='\0';
+
+ switch(line[0]) {
+
+ case 'N': { // total number of particles, more for debugging...
+ noParts = atoi(line+2);
+ if(noParts<=0) {
+ errMsg("ControlParticles::initFromTextFile","file '"<<filename<<"' - invalid no of particles "<<noParts);
+ mPartSets.clear(); fclose(infile); return 0;
+ }
+ if(debugRead) printf("CPDEBUG%d no parts '%d'\n",lineCnt, noParts );
+ haveNo = true;
+ } break;
+
+ case 'T': { // global time scale
+ timeScale *= (LbmFloat)atof(line+2);
+ if(debugRead) printf("ControlParticles::initFromTextFile - line %d , set timescale '%f', org %f\n",lineCnt, timeScale , _initTimeScale);
+ if(timeScale==0.) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error: timescale = 0.! reseting to 1 ...\n",lineCnt); timeScale=1.; }
+ haveScale = true;
+ } break;
+
+ case 'I': { // influence settings, overrides others as of now...
+ float val = (LbmFloat)atof(line+3);
+ const char *setvar = "[invalid]";
+ switch(line[1]) {
+ //case 'f': { _influenceFalloff = val; setvar = "falloff"; } break;
+ case 't': { _influenceTangential = val; setvar = "tangential"; } break;
+ case 'a': { _influenceAttraction = val; setvar = "attraction"; } break;
+ case 'v': { _influenceVelocity = val; setvar = "velocity"; } break;
+ case 'm': { _influenceMaxdist = val; setvar = "maxdist"; } break;
+ default:
+ fprintf(stdout,"ControlParticles::initFromTextFile (%s) - line %d , invalid influence setting %c, %f\n",filename.c_str() ,lineCnt, line[1], val);
+ }
+ if(debugRead) printf("CPDEBUG%d set influence '%s'=%f \n",lineCnt, setvar, val);
+ } break;
+
+ case 'R': { // radius settings, overrides others as of now...
+ float val = (LbmFloat)atof(line+3);
+ const char *setvar = "[invalid]";
+ switch(line[1]) {
+ case 'a': { _radiusAtt = val; setvar = "r_attraction"; } break;
+ case 'v': { _radiusVel = val; setvar = "r_velocity"; } break;
+ case 'm': { _radiusMaxd = val; setvar = "r_maxdist"; } break;
+ default:
+ fprintf(stdout,"ControlParticles::initFromTextFile (%s) - line %d , invalid influence setting %c, %f\n",filename.c_str() ,lineCnt, line[1], val);
+ }
+ if(debugRead) printf("CPDEBUG%d set influence '%s'=%f \n",lineCnt, setvar, val);
+ } break;
+
+ case 'S': { // new particle set at time T
+ ControlParticleSet cps;
+ mPartSets.push_back(cps);
+ setCnt = (int)mPartSets.size()-1;
+
+ LbmFloat val = (LbmFloat)atof(line+2);
+ mPartSets[setCnt].time = val * timeScale;
+ if(debugRead) printf("CPDEBUG%d new set, time '%f', %d\n",lineCnt, mPartSets[setCnt].time, setCnt );
+ haveTime = true;
+ partCnt = -1;
+ } break;
+
+ case 'P': // new particle with pos
+ case 'n': { // new particle without pos
+ if((!haveTime)||(setCnt<0)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error: set missing!\n",lineCnt); abortParse=true; break; }
+ partCnt++;
+ if(partCnt>=noParts) {
+ if(debugRead) printf("CPDEBUG%d partset done \n",lineCnt);
+ haveTime = false;
+ } else {
+ ControlParticle p; p.reset();
+ mPartSets[setCnt].particles.push_back(p);
+ }
+ }
+ // only new part, or new with pos?
+ if(line[0] == 'n') break;
+
+ // particle properties
+
+ case 'p': { // new particle set at time T
+ if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|p: particle missing!\n",lineCnt); abortParse=true; break; }
+ float px=0.,py=0.,pz=0.;
+ if( sscanf(line+2,"%f %f %f",&px,&py,&pz) != 3) {
+ fprintf(stdout,"CPDEBUG%d, unable to parse position!\n",lineCnt); abortParse=true; break;
+ }
+ if(!(finite(px)&&finite(py)&&finite(pz))) { px=py=pz=0.; }
+ LASTCP.pos[0] = px;
+ LASTCP.pos[1] = py;
+ LASTCP.pos[2] = pz;
+ if(debugRead) printf("CPDEBUG%d part%d,%d: position %f,%f,%f \n",lineCnt,setCnt,partCnt, px,py,pz);
+ } break;
+
+ case 's': { // particle size
+ if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|s: particle missing!\n",lineCnt); abortParse=true; break; }
+ float ps=1.;
+ if( sscanf(line+2,"%f",&ps) != 1) {
+ fprintf(stdout,"CPDEBUG%d, unable to parse size!\n",lineCnt); abortParse=true; break;
+ }
+ if(!(finite(ps))) { ps=0.; }
+ LASTCP.size = ps;
+ if(debugRead) printf("CPDEBUG%d part%d,%d: size %f \n",lineCnt,setCnt,partCnt, ps);
+ } break;
+
+ case 'i': { // particle influence
+ if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|i: particle missing!\n",lineCnt); abortParse=true; break; }
+ float pinf=1.;
+ if( sscanf(line+2,"%f",&pinf) != 1) {
+ fprintf(stdout,"CPDEBUG%d, unable to parse size!\n",lineCnt); abortParse=true; break;
+ }
+ if(!(finite(pinf))) { pinf=0.; }
+ LASTCP.influence = pinf;
+ if(debugRead) printf("CPDEBUG%d part%d,%d: influence %f \n",lineCnt,setCnt,partCnt, pinf);
+ } break;
+
+ case 'a': { // rotation axis
+ if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|a: particle missing!\n",lineCnt); abortParse=true; break; }
+ float px=0.,py=0.,pz=0.;
+ if( sscanf(line+2,"%f %f %f",&px,&py,&pz) != 3) {
+ fprintf(stdout,"CPDEBUG%d, unable to parse rotaxis!\n",lineCnt); abortParse=true; break;
+ }
+ if(!(finite(px)&&finite(py)&&finite(pz))) { px=py=pz=0.; }
+ LASTCP.rotaxis[0] = px;
+ LASTCP.rotaxis[1] = py;
+ LASTCP.rotaxis[2] = pz;
+ if(debugRead) printf("CPDEBUG%d part%d,%d: rotaxis %f,%f,%f \n",lineCnt,setCnt,partCnt, px,py,pz);
+ } break;
+
+
+ default:
+ if(debugRead) printf("CPDEBUG%d ignored: '%s'\n",lineCnt, line );
+ break;
+ }
+ }
+ if(debugRead && abortParse) printf("CPDEBUG aborted parsing after set... %d\n",(int)mPartSets.size() );
+
+ // sanity check
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ if( (int)mPartSets[i].particles.size()!=noParts) {
+ fprintf(stdout,"ControlParticles::initFromTextFile (%s) - invalid no of particles in set %d, is:%d, shouldbe:%d \n",filename.c_str() ,i,(int)mPartSets[i].particles.size(), noParts);
+ mPartSets.clear();
+ fclose(infile);
+ return 0;
+ }
+ }
+
+ // print stats
+ printf("ControlParticles::initFromTextFile (%s): Read %d sets, each %d particles\n",filename.c_str() ,
+ (int)mPartSets.size(), noParts );
+ if(mPartSets.size()>0) {
+ printf("ControlParticles::initFromTextFile (%s): Time: %f,%f\n",filename.c_str() ,mPartSets[0].time, mPartSets[mPartSets.size()-1].time );
+ }
+
+ // done...
+ fclose(infile);
+ applyTrafos();
+ return 1;
+}
+
+
+int ControlParticles::initFromTextFileOld(string filename)
+{
+ const bool debugRead = false;
+ char line[LINE_LEN];
+ line[LINE_LEN-1] = '\0';
+ mPartSets.clear();
+ if(filename.size()<1) return 0;
+
+ FILE *infile = fopen(filename.c_str(), "r");
+ if(!infile) {
+ fprintf(stdout,"ControlParticles::initFromTextFileOld - unable to open '%s'\n",filename.c_str() );
+ return 0;
+ }
+
+ int haveNo = false;
+ int haveScale = false;
+ int haveTime = false;
+ int noParts = -1;
+ int coordCnt = 0;
+ int partCnt = 0;
+ int setCnt = 0;
+ ControlParticle p; p.reset();
+ // scale times by constant factor while reading
+ LbmFloat timeScale= 1.0;
+ int lineCnt = 0;
+
+ while(!feof(infile)) {
+ lineCnt++;
+ fgets(line, LINE_LEN, infile);
+
+ if(debugRead) printf("\nDEBUG%d r '%s'\n",lineCnt, line);
+
+ if(!line) continue;
+ int len = (int)strlen(line);
+
+ // skip empty lines and comments (#,//)
+ if(len<1) continue;
+ if( (line[0]=='#') || (line[0]=='\n') ) continue;
+ if((len>1) && (line[0]=='/' && line[1]=='/')) continue;
+
+ // debug remove newline
+ if((len>=1)&&(line[len-1]=='\n')) line[len-1]='\0';
+
+ // first read no. of particles
+ if(!haveNo) {
+ noParts = atoi(line);
+ if(noParts<=0) {
+ fprintf(stdout,"ControlParticles::initFromTextFileOld - invalid no of particles %d\n",noParts);
+ mPartSets.clear();
+ fclose(infile);
+ return 0;
+ }
+ if(debugRead) printf("DEBUG%d noparts '%d'\n",lineCnt, noParts );
+ haveNo = true;
+ }
+
+ // then read time scale
+ else if(!haveScale) {
+ timeScale *= (LbmFloat)atof(line);
+ if(debugRead) printf("DEBUG%d tsc '%f', org %f\n",lineCnt, timeScale , _initTimeScale);
+ haveScale = true;
+ }
+
+ // then get set time
+ else if(!haveTime) {
+ ControlParticleSet cps;
+ mPartSets.push_back(cps);
+ setCnt = (int)mPartSets.size()-1;
+
+ LbmFloat val = (LbmFloat)atof(line);
+ mPartSets[setCnt].time = val * timeScale;
+ if(debugRead) printf("DEBUG%d time '%f', %d\n",lineCnt, mPartSets[setCnt].time, setCnt );
+ haveTime = true;
+ }
+
+ // default read all parts
+ else {
+ LbmFloat val = (LbmFloat)atof(line);
+ if(debugRead) printf("DEBUG: l%d s%d,particle%d '%f' %d,%d/%d\n",lineCnt,(int)mPartSets.size(),(int)mPartSets[setCnt].particles.size(), val ,coordCnt,partCnt,noParts);
+ p.pos[coordCnt] = val;
+ coordCnt++;
+ if(coordCnt>=3) {
+ mPartSets[setCnt].particles.push_back(p);
+ p.reset();
+ coordCnt=0;
+ partCnt++;
+ }
+ if(partCnt>=noParts) {
+ partCnt = 0;
+ haveTime = false;
+ }
+ //if(debugRead) printf("DEBUG%d par2 %d,%d/%d\n",lineCnt, coordCnt,partCnt,noParts);
+ }
+ //read pos, vel ...
+ }
+
+ // sanity check
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ if( (int)mPartSets[i].particles.size()!=noParts) {
+ fprintf(stdout,"ControlParticles::initFromTextFileOld - invalid no of particles in set %d, is:%d, shouldbe:%d \n",i,(int)mPartSets[i].particles.size(), noParts);
+ mPartSets.clear();
+ fclose(infile);
+ return 0;
+ }
+ }
+ // print stats
+ printf("ControlParticles::initFromTextFileOld: Read %d sets, each %d particles\n",
+ (int)mPartSets.size(), noParts );
+ if(mPartSets.size()>0) {
+ printf("ControlParticles::initFromTextFileOld: Time: %f,%f\n",mPartSets[0].time, mPartSets[mPartSets.size()-1].time );
+ }
+
+ // done...
+ fclose(infile);
+ applyTrafos();
+ return 1;
+}
+
+// load positions & timing from gzipped binary file
+int ControlParticles::initFromBinaryFile(string filename) {
+ mPartSets.clear();
+ if(filename.size()<1) return 0;
+ int fileNotFound=0;
+ int fileFound=0;
+ char ofile[256];
+
+ for(int set=0; ((set<10000)&&(fileNotFound<10)); set++) {
+ snprintf(ofile,256,"%s%04d.gz",filename.c_str(),set);
+ //errMsg("ControlParticle::initFromBinaryFile","set"<<set<<" notf"<<fileNotFound<<" ff"<<fileFound);
+
+ gzFile gzf;
+ gzf = gzopen(ofile, "rb");
+ if (!gzf) {
+ //errMsg("ControlParticles::initFromBinaryFile","Unable to open file for reading '"<<ofile<<"' ");
+ fileNotFound++;
+ continue;
+ }
+ fileNotFound=0;
+ fileFound++;
+
+ ControlParticleSet cps;
+ mPartSets.push_back(cps);
+ int setCnt = (int)mPartSets.size()-1;
+ //LbmFloat val = (LbmFloat)atof(line+2);
+ mPartSets[setCnt].time = (gfxReal)set;
+
+ int totpart = 0;
+ gzread(gzf, &totpart, sizeof(totpart));
+
+ for(int a=0; a<totpart; a++) {
+ int ptype=0;
+ float psize=0.0;
+ ntlVec3Gfx ppos,pvel;
+ gzread(gzf, &ptype, sizeof( ptype ));
+ gzread(gzf, &psize, sizeof( float ));
+
+ for(int j=0; j<3; j++) { gzread(gzf, &ppos[j], sizeof( float )); }
+ for(int j=0; j<3; j++) { gzread(gzf, &pvel[j], sizeof( float )); }
+
+ ControlParticle p;
+ p.reset();
+ p.pos = vec2L(ppos);
+ mPartSets[setCnt].particles.push_back(p);
+ }
+
+ gzclose(gzf);
+ //errMsg("ControlParticle::initFromBinaryFile","Read set "<<ofile<<", #"<<mPartSets[setCnt].particles.size() ); // DEBUG
+ } // sets
+
+ if(fileFound==0) return 0;
+ applyTrafos();
+ return 1;
+}
+
+int globCPIProblems =0;
+bool ControlParticles::checkPointInside(ntlTree *tree, ntlVec3Gfx org, gfxReal &distance) {
+ // warning - stripped down version of geoInitCheckPointInside
+ const int globGeoInitDebug = 0;
+ const int flags = FGI_FLUID;
+ org += ntlVec3Gfx(0.0001);
+ ntlVec3Gfx dir = ntlVec3Gfx(1.0, 0.0, 0.0);
+ int OId = -1;
+ ntlRay ray(org, dir, 0, 1.0, NULL);
+ bool done = false;
+ bool inside = false;
+ int mGiObjInside = 0;
+ LbmFloat mGiObjDistance = -1.0;
+ LbmFloat giObjFirstHistSide = 0;
+
+ // if not inside, return distance to first hit
+ gfxReal firstHit=-1.0;
+ int firstOId = -1;
+ if(globGeoInitDebug) errMsg("IIIstart"," isect "<<org);
+
+ while(!done) {
+ // find first inside intersection
+ ntlTriangle *triIns = NULL;
+ distance = -1.0;
+ ntlVec3Gfx normal(0.0);
+ tree->intersectX(ray,distance,normal, triIns, flags, true);
+ if(triIns) {
+ ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
+ LbmFloat orientation = dot(normal, dir);
+ OId = triIns->getObjectId();
+ if(orientation<=0.0) {
+ // outside hit
+ normal *= -1.0;
+ mGiObjInside++;
+ if(giObjFirstHistSide==0) giObjFirstHistSide = 1;
+ if(globGeoInitDebug) errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
+ } else {
+ // inside hit
+ mGiObjInside++;
+ if(mGiObjDistance<0.0) mGiObjDistance = distance;
+ if(globGeoInitDebug) errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
+ if(giObjFirstHistSide==0) giObjFirstHistSide = -1;
+ }
+ norg += normal * getVecEpsilon();
+ ray = ntlRay(norg, dir, 0, 1.0, NULL);
+ // remember first hit distance, in case we're not
+ // inside anything
+ if(firstHit<0.0) {
+ firstHit = distance;
+ firstOId = OId;
+ }
+ } else {
+ // no more intersections... return false
+ done = true;
+ }
+ }
+
+ distance = -1.0;
+ if(mGiObjInside>0) {
+ bool mess = false;
+ if((mGiObjInside%2)==1) {
+ if(giObjFirstHistSide != -1) mess=true;
+ } else {
+ if(giObjFirstHistSide != 1) mess=true;
+ }
+ if(mess) {
+ // ?
+ //errMsg("IIIproblem","At "<<org<<" obj inside:"<<mGiObjInside<<" firstside:"<<giObjFirstHistSide );
+ globCPIProblems++;
+ mGiObjInside++; // believe first hit side...
+ }
+ }
+
+ if(globGeoInitDebug) errMsg("CHIII"," ins="<<mGiObjInside<<" t"<<mGiObjDistance<<" d"<<distance);
+ if(((mGiObjInside%2)==1)&&(mGiObjDistance>0.0)) {
+ if( (distance<0.0) || // first intersection -> good
+ ((distance>0.0)&&(distance>mGiObjDistance)) // more than one intersection -> use closest one
+ ) {
+ distance = mGiObjDistance;
+ OId = 0;
+ inside = true;
+ }
+ }
+
+ if(!inside) {
+ distance = firstHit;
+ OId = firstOId;
+ }
+ if(globGeoInitDebug) errMsg("CHIII","ins"<<inside<<" fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
+
+ return inside;
+}
+int ControlParticles::initFromMVCMesh(string filename) {
+ myTime_t mvmstart = getTime();
+ ntlGeometryObjModel *model = new ntlGeometryObjModel();
+ int gid=1;
+ char infile[256];
+ vector<ntlTriangle> triangles;
+ vector<ntlVec3Gfx> vertices;
+ vector<ntlVec3Gfx> normals;
+ snprintf(infile,256,"%s.bobj.gz", filename.c_str() );
+ model->loadBobjModel(string(infile));
+ model->setLoaded(true);
+ model->setGeoInitId(gid);
+ model->setGeoInitType(FGI_FLUID);
+ debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG,"infile:"<<string(infile) ,4);
+
+ //getTriangles(double t, vector<ntlTriangle> *triangles, vector<ntlVec3Gfx> *vertices, vector<ntlVec3Gfx> *normals, int objectId );
+ model->getTriangles(mCPSTimeStart, &triangles, &vertices, &normals, 1 );
+ debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG," tris:"<<triangles.size()<<" verts:"<<vertices.size()<<" norms:"<<normals.size() , 2);
+
+ // valid mesh?
+ if(triangles.size() <= 0) {
+ return 0;
+ }
+
+ ntlRenderGlobals *glob = new ntlRenderGlobals;
+ ntlScene *genscene = new ntlScene( glob, false );
+ genscene->addGeoClass(model);
+ genscene->addGeoObject(model);
+ genscene->buildScene(0., false);
+ char treeFlag = (1<<(4+gid));
+
+ ntlTree *tree = new ntlTree(
+ 15, 8, // TREEwarning - fixed values for depth & maxtriangles here...
+ genscene, treeFlag );
+
+ // TODO? use params
+ ntlVec3Gfx start,end;
+ model->getExtends(start,end);
+
+ LbmFloat width = mCPSWidth;
+ if(width<=LBM_EPSILON) { errMsg("ControlParticles::initFromMVMCMesh","Invalid mCPSWidth! "<<mCPSWidth); width=mCPSWidth=0.1; }
+ ntlVec3Gfx org = start+ntlVec3Gfx(width*0.5);
+ gfxReal distance = -1.;
+ vector<ntlVec3Gfx> inspos;
+ int approxmax = (int)( ((end[0]-start[0])/width)*((end[1]-start[1])/width)*((end[2]-start[2])/width) );
+
+ debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG,"start"<<start<<" end"<<end<<" w="<<width<<" maxp:"<<approxmax, 5);
+ while(org[2]<end[2]) {
+ while(org[1]<end[1]) {
+ while(org[0]<end[0]) {
+ if(checkPointInside(tree, org, distance)) {
+ inspos.push_back(org);
+ //inspos.push_back(org+ntlVec3Gfx(width));
+ //inspos.push_back(start+end*0.5);
+ }
+ // TODO optimize, use distance
+ org[0] += width;
+ }
+ org[1] += width;
+ org[0] = start[0];
+ }
+ org[2] += width;
+ org[1] = start[1];
+ }
+ debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG,"points: "<<inspos.size()<<" initproblems: "<<globCPIProblems,5 );
+
+ MeanValueMeshCoords mvm;
+ mvm.calculateMVMCs(vertices,triangles, inspos, mCPSWeightFac);
+ vector<ntlVec3Gfx> ninspos;
+ mvm.transfer(vertices, ninspos);
+
+ // init first set, check dist
+ ControlParticleSet firstcps; //T
+ mPartSets.push_back(firstcps);
+ mPartSets[mPartSets.size()-1].time = (gfxReal)0.;
+ vector<bool> useCP;
+ bool debugPos=false;
+
+ for(int i=0; i<(int)inspos.size(); i++) {
+ ControlParticle p; p.reset();
+ p.pos = vec2L(inspos[i]);
+ //errMsg("COMP "," "<<inspos[i]<<" vs "<<ninspos[i] );
+ double cpdist = norm(inspos[i]-ninspos[i]);
+ bool usecpv = true;
+ if(debugPos) errMsg("COMP "," "<<cpdist<<usecpv);
+
+ mPartSets[mPartSets.size()-1].particles.push_back(p);
+ useCP.push_back(usecpv);
+ }
+
+ // init further sets, temporal mesh sampling
+ double tsampling = mCPSTimestep;
+ int totcnt = (int)( (mCPSTimeEnd-mCPSTimeStart)/tsampling ), tcnt=0;
+ for(double t=mCPSTimeStart+tsampling; ((t<mCPSTimeEnd) && (ninspos.size()>0.)); t+=tsampling) {
+ ControlParticleSet nextcps; //T
+ mPartSets.push_back(nextcps);
+ mPartSets[mPartSets.size()-1].time = (gfxReal)t;
+
+ vertices.clear(); triangles.clear(); normals.clear();
+ model->getTriangles(t, &triangles, &vertices, &normals, 1 );
+ mvm.transfer(vertices, ninspos);
+ if(tcnt%(totcnt/10)==1) debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"Transferring animation, frame: "<<tcnt<<"/"<<totcnt,5 );
+ tcnt++;
+ for(int i=0; i<(int)ninspos.size(); i++) {
+ if(debugPos) errMsg("COMP "," "<<norm(inspos[i]-ninspos[i]) );
+ if(useCP[i]) {
+ ControlParticle p; p.reset();
+ p.pos = vec2L(ninspos[i]);
+ mPartSets[mPartSets.size()-1].particles.push_back(p);
+ }
+ }
+ }
+
+ applyTrafos();
+
+ myTime_t mvmend = getTime();
+ debMsgStd("ControlParticle::initFromMVMCMesh",DM_MSG,"t:"<<getTimeString(mvmend-mvmstart)<<" ",7 );
+ delete tree;
+ delete genscene;
+ delete glob;
+//exit(1); // DEBUG
+ return 1;
+}
+
+#define TRISWAP(v,a,b) { LbmFloat tmp = (v)[b]; (v)[b]=(v)[a]; (v)[a]=tmp; }
+#define TRISWAPALL(v,a,b) { \
+ TRISWAP( (v).pos ,a,b ); \
+ TRISWAP( (v).vel ,a,b ); \
+ TRISWAP( (v).rotaxis ,a,b ); }
+
+// helper function for LBM 2D -> swap Y and Z components everywhere
+void ControlParticles::swapCoords(int a, int b) {
+ //return;
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ for(int j=0; j<(int)mPartSets[i].particles.size(); j++) {
+ TRISWAPALL( mPartSets[i].particles[j],a,b );
+ }
+ }
+}
+
+// helper function for LBM 2D -> mirror time
+void ControlParticles::mirrorTime() {
+ LbmFloat maxtime = mPartSets[mPartSets.size()-1].time;
+ const bool debugTimeswap = false;
+
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ mPartSets[i].time = maxtime - mPartSets[i].time;
+ }
+
+ for(int i=0; i<(int)mPartSets.size()/2; i++) {
+ ControlParticleSet cps = mPartSets[i];
+ if(debugTimeswap) errMsg("TIMESWAP", " s"<<i<<","<<mPartSets[i].time<<" and s"<<(mPartSets.size()-1-i)<<","<< mPartSets[mPartSets.size()-1-i].time <<" mt:"<<maxtime );
+ mPartSets[i] = mPartSets[mPartSets.size()-1-i];
+ mPartSets[mPartSets.size()-1-i] = cps;
+ }
+
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ if(debugTimeswap) errMsg("TIMESWAP", "done: s"<<i<<","<<mPartSets[i].time<<" "<<mPartSets[i].particles.size() );
+ }
+}
+
+// apply init transformations
+void ControlParticles::applyTrafos() {
+ // apply trafos
+ for(int i=0; i<(int)mPartSets.size(); i++) {
+ mPartSets[i].time *= _initTimeScale;
+ /*for(int j=0; j<(int)mPartSets[i].particles.size(); j++) {
+ for(int k=0; k<3; k++) {
+ mPartSets[i].particles[j].pos[k] *= _initPartScale[k];
+ mPartSets[i].particles[j].pos[k] += _initPartOffset[k];
+ }
+ } now done in initarray */
+ }
+
+ // mirror coords...
+ for(int l=0; l<(int)_initMirror.length(); l++) {
+ switch(_initMirror[l]) {
+ case 'X':
+ case 'x':
+ //printf("ControlParticles::applyTrafos - mirror x\n");
+ swapCoords(1,2);
+ break;
+ case 'Y':
+ case 'y':
+ //printf("ControlParticles::applyTrafos - mirror y\n");
+ swapCoords(0,2);
+ break;
+ case 'Z':
+ case 'z':
+ //printf("ControlParticles::applyTrafos - mirror z\n");
+ swapCoords(0,1);
+ break;
+ case 'T':
+ case 't':
+ //printf("ControlParticles::applyTrafos - mirror time\n");
+ mirrorTime();
+ break;
+ case ' ':
+ case '-':
+ case '\n':
+ break;
+ default:
+ //printf("ControlParticles::applyTrafos - mirror unknown %c !?\n", _initMirror[l] );
+ break;
+ }
+ }
+
+ // reset 2d positions
+#if (CP_PROJECT2D==1) && ( defined(MAIN_2D) || LBMDIM==2 )
+ for(size_t j=0; j<mPartSets.size(); j++)
+ for(size_t i=0; i<mPartSets[j].particles.size(); i++) {
+ // DEBUG
+ mPartSets[j].particles[i].pos[1] = 0.f;
+ }
+#endif
+
+#if defined(LBMDIM)
+ //? if( (getenv("ELBEEM_CPINFILE")) || (getenv("ELBEEM_CPOUTFILE")) ){
+ // gui control test, don swap...
+ //? } else {
+ //? swapCoords(1,2); // LBM 2D -> swap Y and Z components everywhere
+ //? }
+#endif
+
+ initTime(0.f, 0.f);
+}
+
+#undef TRISWAP
+
+// --------------------------------------------------------------------------
+// init for a given time
+void ControlParticles::initTime(LbmFloat t, LbmFloat dt)
+{
+ //fprintf(stdout, "CPINITTIME init %f\n",t);
+ _currTime = t;
+ if(mPartSets.size()<1) return;
+
+ // init zero velocities
+ initTimeArray(t, _particles);
+
+ // calculate velocities from prev. timestep?
+ if(dt>0.) {
+ _currTimestep = dt;
+ std::vector<ControlParticle> prevparts;
+ initTimeArray(t-dt, prevparts);
+ LbmFloat invdt = 1.0/dt;
+ for(size_t j=0; j<_particles.size(); j++) {
+ ControlParticle &p = _particles[j];
+ ControlParticle &prevp = prevparts[j];
+ for(int k=0; k<3; k++) {
+ p.pos[k] *= _initPartScale[k];
+ p.pos[k] += _initPartOffset[k];
+ prevp.pos[k] *= _initLastPartScale[k];
+ prevp.pos[k] += _initLastPartOffset[k];
+ }
+ p.vel = (p.pos - prevp.pos)*invdt;
+ }
+
+ if(0) {
+ LbmVec avgvel(0.);
+ for(size_t j=0; j<_particles.size(); j++) {
+ avgvel += _particles[j].vel;
+ }
+ avgvel /= (LbmFloat)_particles.size();
+ //fprintf(stdout," AVGVEL %f,%f,%f \n",avgvel[0],avgvel[1],avgvel[2]); // DEBUG
+ }
+ }
+}
+
+// helper, init given array
+void ControlParticles::initTimeArray(LbmFloat t, std::vector<ControlParticle> &parts) {
+ if(mPartSets.size()<1) return;
+
+ if(parts.size()!=mPartSets[0].particles.size()) {
+ //fprintf(stdout,"PRES \n");
+ parts.resize(mPartSets[0].particles.size());
+ // TODO reset all?
+ for(size_t j=0; j<parts.size(); j++) {
+ parts[j].reset();
+ }
+ }
+ if(parts.size()<1) return;
+
+ // debug inits
+ if(mDebugInit==1) {
+ // hard coded circle init
+ for(size_t j=0; j<mPartSets[0].particles.size(); j++) {
+ ControlParticle p = mPartSets[0].particles[j];
+ // remember old
+ p.density = parts[j].density;
+ p.densityWeight = parts[j].densityWeight;
+ p.avgVel = parts[j].avgVel;
+ p.avgVelAcc = parts[j].avgVelAcc;
+ p.avgVelWeight = parts[j].avgVelWeight;
+ LbmVec ppos(0.); { // DEBUG
+ const float tscale=10.;
+ const float tprevo = 0.33;
+ const LbmVec toff(50,50,0);
+ const LbmVec oscale(30,30,0);
+ ppos[0] = cos(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[0] + toff[0];
+ ppos[1] = -sin(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[1] + toff[1];
+ ppos[2] = toff[2]; } // DEBUG
+ p.pos = ppos;
+ parts[j] = p;
+ //errMsg("ControlParticle::initTimeArray","j:"<<j<<" p:"<<parts[j].pos );
+ }
+ return;
+ }
+ else if(mDebugInit==2) {
+ // hard coded spiral init
+ const float tscale=-10.;
+ const float tprevo = 0.33;
+ LbmVec toff(50,0,-50);
+ const LbmVec oscale(20,20,0);
+ toff[2] += 30. * t +30.;
+ for(size_t j=0; j<mPartSets[0].particles.size(); j++) {
+ ControlParticle p = mPartSets[0].particles[j];
+ // remember old
+ p.density = parts[j].density;
+ p.densityWeight = parts[j].densityWeight;
+ p.avgVel = parts[j].avgVel;
+ p.avgVelAcc = parts[j].avgVelAcc;
+ p.avgVelWeight = parts[j].avgVelWeight;
+ LbmVec ppos(0.);
+ ppos[1] = toff[2];
+ LbmFloat zscal = (ppos[1]+100.)/200.;
+ ppos[0] = cos(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[0]*zscal + toff[0];
+ ppos[2] = -sin(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[1]*zscal + toff[1];
+ p.pos = ppos;
+ parts[j] = p;
+
+ toff[2] += 0.25;
+ }
+ return;
+ }
+
+ // use first set
+ if((t<=mPartSets[0].time)||(mPartSets.size()==1)) {
+ //fprintf(stdout,"PINI %f \n", t);
+ //parts = mPartSets[0].particles;
+ const int i=0;
+ for(size_t j=0; j<mPartSets[i].particles.size(); j++) {
+ ControlParticle p = mPartSets[i].particles[j];
+ // remember old
+ p.density = parts[j].density;
+ p.densityWeight = parts[j].densityWeight;
+ p.avgVel = parts[j].avgVel;
+ p.avgVelAcc = parts[j].avgVelAcc;
+ p.avgVelWeight = parts[j].avgVelWeight;
+ parts[j] = p;
+ }
+ return;
+ }
+
+ for(int i=0; i<(int)mPartSets.size()-1; i++) {
+ if((mPartSets[i].time<=t) && (mPartSets[i+1].time>t)) {
+ LbmFloat d = mPartSets[i+1].time-mPartSets[i].time;
+ LbmFloat f = (t-mPartSets[i].time)/d;
+ LbmFloat omf = 1.0f - f;
+
+ for(size_t j=0; j<mPartSets[i].particles.size(); j++) {
+ ControlParticle *src1=&mPartSets[i ].particles[j];
+ ControlParticle *src2=&mPartSets[i+1].particles[j];
+ ControlParticle &p = parts[j];
+ // do linear interpolation
+ p.pos = src1->pos * omf + src2->pos *f;
+ p.vel = LbmVec(0.); // reset, calculated later on src1->vel * omf + src2->vel *f;
+ p.rotaxis = src1->rotaxis * omf + src2->rotaxis *f;
+ p.influence = src1->influence * omf + src2->influence *f;
+ p.size = src1->size * omf + src2->size *f;
+ // dont modify: density, densityWeight
+ }
+ }
+ }
+
+ // after last?
+ if(t>=mPartSets[ mPartSets.size() -1 ].time) {
+ //parts = mPartSets[ mPartSets.size() -1 ].particles;
+ const int i= (int)mPartSets.size() -1;
+ for(size_t j=0; j<mPartSets[i].particles.size(); j++) {
+ ControlParticle p = mPartSets[i].particles[j];
+ // restore
+ p.density = parts[j].density;
+ p.densityWeight = parts[j].densityWeight;
+ p.avgVel = parts[j].avgVel;
+ p.avgVelAcc = parts[j].avgVelAcc;
+ p.avgVelWeight = parts[j].avgVelWeight;
+ parts[j] = p;
+ }
+ }
+}
+
+
+
+
+// --------------------------------------------------------------------------
+
+#define DEBUG_MODVEL 0
+
+// recalculate
+void ControlParticles::calculateKernelWeight() {
+ const bool debugKernel = true;
+
+ // calculate kernel area with respect to particlesize/cellsize
+ LbmFloat kernelw = -1.;
+ LbmFloat kernelnorm = -1.;
+ LbmFloat krad = (_radiusAtt*0.75); // FIXME use real cone approximation...?
+ //krad = (_influenceFalloff*1.);
+#if (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+ kernelw = CP_PI*krad*krad;
+ kernelnorm = 1.0 / (_fluidSpacing * _fluidSpacing);
+#else // 2D
+ kernelw = CP_PI*krad*krad*krad* (4./3.);
+ kernelnorm = 1.0 / (_fluidSpacing * _fluidSpacing * _fluidSpacing);
+#endif // MAIN_2D
+
+ if(debugKernel) debMsgStd("ControlParticles::calculateKernelWeight",DM_MSG,"kw"<<kernelw<<", norm"<<
+ kernelnorm<<", w*n="<<(kernelw*kernelnorm)<<", rad"<<krad<<", sp"<<_fluidSpacing<<" ", 7);
+ LbmFloat kernelws = kernelw*kernelnorm;
+ _kernelWeight = kernelws;
+ if(debugKernel) debMsgStd("ControlParticles::calculateKernelWeight",DM_MSG,"influence f="<<_radiusAtt<<" t="<<
+ _influenceTangential<<" a="<<_influenceAttraction<<" v="<<_influenceVelocity<<" kweight="<<_kernelWeight, 7);
+ if(_kernelWeight<=0.) {
+ errMsg("ControlParticles::calculateKernelWeight", "invalid kernel! "<<_kernelWeight<<", resetting");
+ _kernelWeight = 1.;
+ }
+}
+
+void
+ControlParticles::prepareControl(LbmFloat simtime, LbmFloat dt, ControlParticles *motion) {
+ debMsgStd("ControlParticle::prepareControl",DM_MSG," simtime="<<simtime<<" dt="<<dt<<" ", 5);
+
+ //fprintf(stdout,"PREPARE \n");
+ LbmFloat avgdw = 0.;
+ for(size_t i=0; i<_particles.size(); i++) {
+ ControlParticle *cp = &_particles[i];
+
+ if(this->getInfluenceAttraction()<0.) {
+ cp->density=
+ cp->densityWeight = 1.0;
+ continue;
+ }
+
+ // normalize by kernel
+ //cp->densityWeight = (1.0 - (cp->density / _kernelWeight)); // store last
+#if (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+ cp->densityWeight = (1.0 - (cp->density / (_kernelWeight*cp->size*cp->size) )); // store last
+#else // 2D
+ cp->densityWeight = (1.0 - (cp->density / (_kernelWeight*cp->size*cp->size*cp->size) )); // store last
+#endif // MAIN_2D
+
+ if(i<10) debMsgStd("ControlParticle::prepareControl",DM_MSG,"kernelDebug i="<<i<<" densWei="<<cp->densityWeight<<" 1/kw"<<(1.0/_kernelWeight)<<" cpdensity="<<cp->density, 9 );
+ if(cp->densityWeight<0.) cp->densityWeight=0.;
+ if(cp->densityWeight>1.) cp->densityWeight=1.;
+
+ avgdw += cp->densityWeight;
+ // reset for next step
+ cp->density = 0.;
+
+ if(cp->avgVelWeight>0.) {
+ cp->avgVel = cp->avgVelAcc/cp->avgVelWeight;
+ cp->avgVelWeight = 0.;
+ cp->avgVelAcc = LbmVec(0.,0.,0.);
+ }
+ }
+ //if(debugKernel) for(size_t i=0; i<_particles.size(); i++) { ControlParticle *cp = &_particles[i]; fprintf(stdout,"A %f,%f \n",cp->density,cp->densityWeight); }
+ avgdw /= (LbmFloat)(_particles.size());
+ //if(motion) { printf("ControlParticle::kernel: avgdw:%f, kw%f, sp%f \n", avgdw, _kernelWeight, _fluidSpacing); }
+
+ //if((simtime>=0.) && (simtime != _currTime))
+ initTime(simtime, dt);
+
+ if((motion) && (motion->getSize()>0)){
+ ControlParticle *motionp = motion->getParticle(0);
+ //printf("ControlParticle::prepareControl motion: pos[%f,%f,%f] vel[%f,%f,%f] \n", motionp->pos[0], motionp->pos[1], motionp->pos[2], motionp->vel[0], motionp->vel[1], motionp->vel[2] );
+ for(size_t i=0; i<_particles.size(); i++) {
+ ControlParticle *cp = &_particles[i];
+ cp->pos = cp->pos + motionp->pos;
+ cp->vel = cp->vel + motionp->vel;
+ cp->size = cp->size * motionp->size;
+ cp->influence = cp->size * motionp->influence;
+ }
+ }
+
+ // reset to radiusAtt by default
+ if(_radiusVel==0.) _radiusVel = _radiusAtt;
+ if(_radiusMinMaxd==0.) _radiusMinMaxd = _radiusAtt;
+ if(_radiusMaxd==0.) _radiusMaxd = 2.*_radiusAtt;
+ // has to be radiusVel<radiusAtt<radiusMinMaxd<radiusMaxd
+ if(_radiusVel>_radiusAtt) _radiusVel = _radiusAtt;
+ if(_radiusAtt>_radiusMinMaxd) _radiusAtt = _radiusMinMaxd;
+ if(_radiusMinMaxd>_radiusMaxd) _radiusMinMaxd = _radiusMaxd;
+
+ //printf("ControlParticle::radii vel:%f att:%f min:%f max:%f \n", _radiusVel,_radiusAtt,_radiusMinMaxd,_radiusMaxd);
+ // prepareControl done
+}
+
+void ControlParticles::finishControl(std::vector<ControlForces> &forces, LbmFloat iatt, LbmFloat ivel, LbmFloat imaxd) {
+
+ //const LbmFloat iatt = this->getInfluenceAttraction() * this->getCurrTimestep();
+ //const LbmFloat ivel = this->getInfluenceVelocity();
+ //const LbmFloat imaxd = this->getInfluenceMaxdist() * this->getCurrTimestep();
+ // prepare for usage
+ iatt *= this->getCurrTimestep();
+ ivel *= 1.; // not necessary!
+ imaxd *= this->getCurrTimestep();
+
+ // skip when size=0
+ for(int i=0; i<(int)forces.size(); i++) {
+ if(DEBUG_MODVEL) fprintf(stdout, "CPFORGF %d , wf:%f,f:%f,%f,%f , v:%f,%f,%f \n",i, forces[i].weightAtt, forces[i].forceAtt[0],forces[i].forceAtt[1],forces[i].forceAtt[2], forces[i].forceVel[0], forces[i].forceVel[1], forces[i].forceVel[2] );
+ LbmFloat cfweight = forces[i].weightAtt; // always normalize
+ if((cfweight!=0.)&&(iatt!=0.)) {
+ // multiple kernels, normalize - note this does not normalize in d>r/2 region
+ if(ABS(cfweight)>1.) { cfweight = 1.0/cfweight; }
+ // multiply iatt afterwards to allow stronger force
+ cfweight *= iatt;
+ forces[i].forceAtt *= cfweight;
+ } else {
+ forces[i].weightAtt = 0.;
+ forces[i].forceAtt = LbmVec(0.);
+ }
+
+ if( (cfweight==0.) && (imaxd>0.) && (forces[i].maxDistance>0.) ) {
+ forces[i].forceMaxd *= imaxd;
+ } else {
+ forces[i].maxDistance= 0.;
+ forces[i].forceMaxd = LbmVec(0.);
+ }
+
+ LbmFloat cvweight = forces[i].weightVel; // always normalize
+ if(cvweight>0.) {
+ forces[i].forceVel /= cvweight;
+ forces[i].compAv /= cvweight;
+ // now modify cvweight, and write back
+ // important, cut at 1 - otherwise strong vel. influences...
+ if(cvweight>1.) { cvweight = 1.; }
+ // thus cvweight is in the range of 0..influenceVelocity, currently not normalized by numCParts
+ cvweight *= ivel;
+ if(cvweight<0.) cvweight=0.; if(cvweight>1.) cvweight=1.;
+ // LBM, FIXME todo use relaxation factor
+ //pvel = (cvel*0.5 * cvweight) + (pvel * (1.0-cvweight));
+ forces[i].weightVel = cvweight;
+
+ //errMsg("COMPAV","i"<<i<<" compav"<<forces[i].compAv<<" forcevel"<<forces[i].forceVel<<" ");
+ } else {
+ forces[i].weightVel = 0.;
+ if(forces[i].maxDistance==0.) forces[i].forceVel = LbmVec(0.);
+ forces[i].compAvWeight = 0.;
+ forces[i].compAv = LbmVec(0.);
+ }
+ if(DEBUG_MODVEL) fprintf(stdout, "CPFINIF %d , wf:%f,f:%f,%f,%f , v:%f,%f,%f \n",i, forces[i].weightAtt, forces[i].forceAtt[0],forces[i].forceAtt[1],forces[i].forceAtt[2], forces[i].forceVel[0],forces[i].forceVel[1],forces[i].forceVel[2] );
+ }
+
+ // unused...
+ if(DEBUG_MODVEL) fprintf(stdout,"MFC iatt:%f,%f ivel:%f,%f ifmd:%f,%f \n", iatt,_radiusAtt, ivel,_radiusVel, imaxd, _radiusMaxd);
+ //for(size_t i=0; i<_particles.size(); i++) { ControlParticle *cp = &_particles[i]; fprintf(stdout," %f,%f,%f ",cp->density,cp->densityWeight, (1.0 - (12.0*cp->densityWeight))); }
+ //fprintf(stdout,"\n\nCP DONE \n\n\n");
+}
+
+
+// --------------------------------------------------------------------------
+// calculate forces at given position, and modify velocity
+// according to timestep
+void ControlParticles::calculateCpInfluenceOpt(ControlParticle *cp, LbmVec fluidpos, LbmVec fluidvel, ControlForces *force, LbmFloat fillFactor) {
+ // dont reset, only add...
+ // test distance, simple squared distance reject
+ const LbmFloat cpfo = _radiusAtt*cp->size;
+
+ LbmVec posDelta;
+ if(DEBUG_MODVEL) fprintf(stdout, "CP at %f,%f,%f bef fw:%f, f:%f,%f,%f , vw:%f, v:%f,%f,%f \n",fluidpos[0],fluidpos[1],fluidpos[2], force->weightAtt, force->forceAtt[0], force->forceAtt[1], force->forceAtt[2], force->weightVel, force->forceVel[0], force->forceVel[1], force->forceVel[2]);
+ posDelta = cp->pos - fluidpos;
+#if LBMDIM==2 && (CP_PROJECT2D==1)
+ posDelta[2] = 0.; // project to xy plane, z-velocity should already be gone...
+#endif
+
+ const LbmFloat distsqr = posDelta[0]*posDelta[0]+posDelta[1]*posDelta[1]+posDelta[2]*posDelta[2];
+ if(DEBUG_MODVEL) fprintf(stdout, " Pd at %f,%f,%f d%f \n",posDelta[0],posDelta[1],posDelta[2], distsqr);
+ // cut at influence=0.5 , scaling not really makes sense
+ if(cpfo*cpfo < distsqr) {
+ /*if(cp->influence>0.5) {
+ if(force->weightAtt == 0.) {
+ if(force->maxDistance*force->maxDistance > distsqr) {
+ const LbmFloat dis = sqrtf((float)distsqr);
+ const LbmFloat sc = dis-cpfo;
+ force->maxDistance = dis;
+ force->forceMaxd = (posDelta)*(sc/dis);
+ }
+ } } */
+ return;
+ }
+ force->weightAtt += 1e-6; // for distance
+ force->maxDistance = 0.; // necessary for SPH?
+
+ const LbmFloat pdistance = MAGNITUDE(posDelta);
+ LbmFloat pdistinv = 0.;
+ if(ABS(pdistance)>0.) pdistinv = 1./pdistance;
+ posDelta *= pdistinv;
+
+ LbmFloat falloffAtt = 0.; //CPKernel::kernel(cpfo * 1.0, pdistance);
+ const LbmFloat qac = pdistance / cpfo ;
+ if (qac < 1.0){ // return 0.;
+ if(qac < 0.5) falloffAtt = 1.0f;
+ else falloffAtt = (1.0f - qac) * 2.0f;
+ }
+
+ // vorticity force:
+ // - //LbmVec forceVort;
+ // - //CROSS(forceVort, posDelta, cp->rotaxis);
+ // - //NORMALIZE(forceVort);
+ // - if(falloffAtt>1.0) falloffAtt=1.0;
+
+#if (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+ // fillFactor *= 2.0 *0.75 * pdistance; // 2d>3d sampling
+#endif // (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+
+ LbmFloat signum = getInfluenceAttraction() > 0.0 ? 1.0 : -1.0;
+ cp->density += falloffAtt * fillFactor;
+ force->forceAtt += posDelta *cp->densityWeight *cp->influence *signum;
+ force->weightAtt += falloffAtt*cp->densityWeight *cp->influence;
+
+ LbmFloat falloffVel = 0.; //CPKernel::kernel(cpfo * 1.0, pdistance);
+ const LbmFloat cpfv = _radiusVel*cp->size;
+ if(cpfv*cpfv < distsqr) { return; }
+ const LbmFloat qvc = pdistance / cpfo ;
+ //if (qvc < 1.0){
+ //if(qvc < 0.5) falloffVel = 1.0f;
+ //else falloffVel = (1.0f - qvc) * 2.0f;
+ //}
+ falloffVel = 1.-qvc;
+
+ LbmFloat pvWeight; // = (1.0-cp->densityWeight) * _currTimestep * falloffVel;
+ pvWeight = falloffVel *cp->influence; // std, without density influence
+ //pvWeight *= (1.0-cp->densityWeight); // use inverse density weight
+ //pvWeight *= cp->densityWeight; // test, use density weight
+ LbmVec modvel(0.);
+ modvel += cp->vel * pvWeight;
+ //pvWeight = 1.; modvel = partVel; // DEBUG!?
+
+ if(pvWeight>0.) {
+ force->forceVel += modvel;
+ force->weightVel += pvWeight;
+
+ cp->avgVelWeight += falloffVel;
+ cp->avgVel += fluidvel;
+ }
+ if(DEBUG_MODVEL) fprintf(stdout, "CP at %f,%f,%f aft fw:%f, f:%f,%f,%f , vw:%f, v:%f,%f,%f \n",fluidpos[0],fluidpos[1],fluidpos[2], force->weightAtt, force->forceAtt[0], force->forceAtt[1], force->forceAtt[2], force->weightVel, force->forceVel[0], force->forceVel[1], force->forceVel[2]);
+ return;
+}
+
+void ControlParticles::calculateMaxdForce(ControlParticle *cp, LbmVec fluidpos, ControlForces *force) {
+ if(force->weightAtt != 0.) return; // maxd force off
+ if(cp->influence <= 0.5) return; // ignore
+
+ LbmVec posDelta;
+ //if(DEBUG_MODVEL) fprintf(stdout, "CP at %f,%f,%f bef fw:%f, f:%f,%f,%f , vw:%f, v:%f,%f,%f \n",fluidpos[0],fluidpos[1],fluidpos[2], force->weightAtt, force->forceAtt[0], force->forceAtt[1], force->forceAtt[2], force->weightVel, force->forceVel[0], force->forceVel[1], force->forceVel[2]);
+ posDelta = cp->pos - fluidpos;
+#if LBMDIM==2 && (CP_PROJECT2D==1)
+ posDelta[2] = 0.; // project to xy plane, z-velocity should already be gone...
+#endif
+
+ // dont reset, only add...
+ // test distance, simple squared distance reject
+ const LbmFloat distsqr = posDelta[0]*posDelta[0]+posDelta[1]*posDelta[1]+posDelta[2]*posDelta[2];
+
+ // closer cp found
+ if(force->maxDistance*force->maxDistance < distsqr) return;
+
+ const LbmFloat dmin = _radiusMinMaxd*cp->size;
+ if(distsqr<dmin*dmin) return; // inside min
+ const LbmFloat dmax = _radiusMaxd*cp->size;
+ if(distsqr>dmax*dmax) return; // outside
+
+
+ if(DEBUG_MODVEL) fprintf(stdout, " Pd at %f,%f,%f d%f \n",posDelta[0],posDelta[1],posDelta[2], distsqr);
+ // cut at influence=0.5 , scaling not really makes sense
+ const LbmFloat dis = sqrtf((float)distsqr);
+ //const LbmFloat sc = dis - dmin;
+ const LbmFloat sc = (dis-dmin)/(dmax-dmin); // scale from 0-1
+ force->maxDistance = dis;
+ force->forceMaxd = (posDelta/dis) * sc;
+ //debug errMsg("calculateMaxdForce","pos"<<fluidpos<<" dis"<<dis<<" sc"<<sc<<" dmin"<<dmin<<" maxd"<< force->maxDistance <<" fmd"<<force->forceMaxd );
+ return;
+}
+
--- /dev/null
+// --------------------------------------------------------------------------
+//
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the
+// GNU General Public License. See the file COPYING for details.
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+// control particle classes
+//
+// --------------------------------------------------------------------------
+
+#ifndef CONTROLPARTICLES_H
+#define CONTROLPARTICLES_H
+
+#include "ntl_geometrymodel.h"
+
+// indicator for LBM inclusion
+//#ifndef LBMDIM
+
+//#include <NxFoundation.h>
+//#include <vector>
+//class MultisphGUI;
+//#define NORMALIZE(a) a.normalize()
+//#define MAGNITUDE(a) a.magnitude()
+//#define CROSS(a,b,c) a.cross(b,c)
+//#define ABS(a) (a>0. ? (a) : -(a))
+//#include "cpdefines.h"
+
+//#else // LBMDIM
+
+// use compatibility defines
+//#define NORMALIZE(a) normalize(a)
+//#define MAGNITUDE(a) norm(a)
+//#define CROSS(a,b,c) a=cross(b,c)
+
+//#endif // LBMDIM
+
+#define MAGNITUDE(a) norm(a)
+
+// math.h compatibility
+#define CP_PI ((LbmFloat)3.14159265358979323846)
+
+// project 2d test cases onto plane?
+// if not, 3d distance is used for 2d sim as well
+#define CP_PROJECT2D 1
+
+
+// default init for mincpdist, ControlForces::maxDistance
+#define CPF_MAXDINIT 10000.
+
+// storage of influence for a fluid cell/particle in lbm/sph
+class ControlForces
+{
+public:
+ ControlForces() { };
+ ~ControlForces() {};
+
+ // attraction force
+ LbmFloat weightAtt;
+ LbmVec forceAtt;
+ // velocity influence
+ LbmFloat weightVel;
+ LbmVec forceVel;
+ // maximal distance influence,
+ // first is max. distance to first control particle
+ // second attraction strength
+ LbmFloat maxDistance;
+ LbmVec forceMaxd;
+
+ LbmFloat compAvWeight;
+ LbmVec compAv;
+
+ void resetForces() {
+ weightAtt = weightVel = 0.;
+ maxDistance = CPF_MAXDINIT;
+ forceAtt = forceVel = forceMaxd = LbmVec(0.,0.,0.);
+ compAvWeight=0.; compAv=LbmVec(0.);
+ };
+};
+
+
+// single control particle
+class ControlParticle
+{
+public:
+ ControlParticle() { reset(); };
+ ~ControlParticle() {};
+
+ // control parameters
+
+ // position
+ LbmVec pos;
+ // size (influences influence radius)
+ LbmFloat size;
+ // overall strength of influence
+ LbmFloat influence;
+ // rotation axis
+ LbmVec rotaxis;
+
+ // computed values
+
+ // velocity
+ LbmVec vel;
+ // computed density
+ LbmFloat density;
+ LbmFloat densityWeight;
+
+ LbmVec avgVel;
+ LbmVec avgVelAcc;
+ LbmFloat avgVelWeight;
+
+ // init all zero / defaults
+ void reset();
+};
+
+
+// container for a particle configuration at time t
+class ControlParticleSet
+{
+public:
+
+ // time of particle set
+ LbmFloat time;
+ // particle positions
+ std::vector<ControlParticle> particles;
+
+};
+
+
+// container & management of control particles
+class ControlParticles
+{
+public:
+ ControlParticles();
+ ~ControlParticles();
+
+ // reset datastructures for next influence step
+ // if motion object is given, particle 1 of second system is used for overall
+ // position and speed offset
+ void prepareControl(LbmFloat simtime, LbmFloat dt, ControlParticles *motion);
+ // post control operations
+ void finishControl(std::vector<ControlForces> &forces, LbmFloat iatt, LbmFloat ivel, LbmFloat imaxd);
+ // recalculate
+ void calculateKernelWeight();
+
+ // calculate forces at given position, and modify velocity
+ // according to timestep (from initControl)
+ void calculateCpInfluenceOpt (ControlParticle *cp, LbmVec fluidpos, LbmVec fluidvel, ControlForces *force, LbmFloat fillFactor);
+ void calculateMaxdForce (ControlParticle *cp, LbmVec fluidpos, ControlForces *force);
+
+ // no. of particles
+ inline int getSize() { return (int)_particles.size(); }
+ int getTotalSize();
+ // get particle [i]
+ inline ControlParticle* getParticle(int i){ return &_particles[i]; }
+
+ // set influence parameters
+ void setInfluenceTangential(LbmFloat set) { _influenceTangential=set; }
+ void setInfluenceAttraction(LbmFloat set) { _influenceAttraction=set; }
+ void setInfluenceMaxdist(LbmFloat set) { _influenceMaxdist=set; }
+ // calculate for delta t
+ void setInfluenceVelocity(LbmFloat set, LbmFloat dt);
+ // get influence parameters
+ inline LbmFloat getInfluenceAttraction() { return _influenceAttraction; }
+ inline LbmFloat getInfluenceTangential() { return _influenceTangential; }
+ inline LbmFloat getInfluenceVelocity() { return _influenceVelocity; }
+ inline LbmFloat getInfluenceMaxdist() { return _influenceMaxdist; }
+ inline LbmFloat getCurrTimestep() { return _currTimestep; }
+
+ void setRadiusAtt(LbmFloat set) { _radiusAtt=set; }
+ inline LbmFloat getRadiusAtt() { return _radiusAtt; }
+ void setRadiusVel(LbmFloat set) { _radiusVel=set; }
+ inline LbmFloat getRadiusVel() { return _radiusVel; }
+ void setRadiusMaxd(LbmFloat set) { _radiusMaxd=set; }
+ inline LbmFloat getRadiusMaxd() { return _radiusMaxd; }
+ void setRadiusMinMaxd(LbmFloat set) { _radiusMinMaxd=set; }
+ inline LbmFloat getRadiusMinMaxd() { return _radiusMinMaxd; }
+
+ LbmFloat getControlTimStart();
+ LbmFloat getControlTimEnd();
+
+ // set/get characteristic length (and inverse)
+ void setCharLength(LbmFloat set) { _charLength=set; _charLengthInv=1./_charLength; }
+ inline LbmFloat getCharLength() { return _charLength;}
+ inline LbmFloat getCharLengthInv() { return _charLengthInv;}
+
+ // set init parameters
+ void setInitTimeScale(LbmFloat set) { _initTimeScale = set; };
+ void setInitMirror(string set) { _initMirror = set; };
+ string getInitMirror() { return _initMirror; };
+
+ void setLastOffset(LbmVec set) { _initLastPartOffset = set; };
+ void setLastScale(LbmVec set) { _initLastPartScale = set; };
+ void setOffset(LbmVec set) { _initPartOffset = set; };
+ void setScale(LbmVec set) { _initPartScale = set; };
+
+ // set/get cps params
+ void setCPSWith(LbmFloat set) { mCPSWidth = set; };
+ void setCPSTimestep(LbmFloat set) { mCPSTimestep = set; };
+ void setCPSTimeStart(LbmFloat set) { mCPSTimeStart = set; };
+ void setCPSTimeEnd(LbmFloat set) { mCPSTimeEnd = set; };
+ void setCPSMvmWeightFac(LbmFloat set) { mCPSWeightFac = set; };
+
+ LbmFloat getCPSWith() { return mCPSWidth; };
+ LbmFloat getCPSTimestep() { return mCPSTimestep; };
+ LbmFloat getCPSTimeStart() { return mCPSTimeStart; };
+ LbmFloat getCPSTimeEnd() { return mCPSTimeEnd; };
+ LbmFloat getCPSMvmWeightFac() { return mCPSWeightFac; };
+
+ void setDebugInit(int set) { mDebugInit = set; };
+
+ // set init parameters
+ void setFluidSpacing(LbmFloat set) { _fluidSpacing = set; };
+
+ // load positions & timing from text file
+ int initFromTextFile(string filename);
+ int initFromTextFileOld(string filename);
+ // load positions & timing from gzipped binary file
+ int initFromBinaryFile(string filename);
+ int initFromMVCMesh(string filename);
+ // init an example test case
+ int initExampleSet();
+
+ // init for a given time
+ void initTime(LbmFloat t, LbmFloat dt);
+
+ // blender test init
+ void initBlenderTest();
+
+ int initFromObject(ntlGeometryObjModel *model);
+
+protected:
+ // sets influence params
+ friend class MultisphGUI;
+
+ // tangential and attraction influence
+ LbmFloat _influenceTangential, _influenceAttraction;
+ // direct velocity influence
+ LbmFloat _influenceVelocity;
+ // maximal distance influence
+ LbmFloat _influenceMaxdist;
+
+ // influence radii
+ LbmFloat _radiusAtt, _radiusVel, _radiusMinMaxd, _radiusMaxd;
+
+ // currently valid time & timestep
+ LbmFloat _currTime, _currTimestep;
+ // all particles
+ std::vector<ControlParticle> _particles;
+
+ // particle sets
+ std::vector<ControlParticleSet> mPartSets;
+
+ // additional parameters for initing particles
+ LbmFloat _initTimeScale;
+ LbmVec _initPartOffset;
+ LbmVec _initPartScale;
+ LbmVec _initLastPartOffset;
+ LbmVec _initLastPartScale;
+ // mirror particles for loading?
+ string _initMirror;
+
+ // row spacing paramter, e.g. use for approximation of kernel area/volume
+ LbmFloat _fluidSpacing;
+ // save current kernel weight
+ LbmFloat _kernelWeight;
+ // charateristic length in world coordinates for normalizatioon of forces
+ LbmFloat _charLength, _charLengthInv;
+
+
+ /*! do ani mesh CPS */
+ void calculateCPS(string filename);
+ //! ani mesh cps params
+ ntlVec3Gfx mvCPSStart, mvCPSEnd;
+ gfxReal mCPSWidth, mCPSTimestep;
+ gfxReal mCPSTimeStart, mCPSTimeEnd;
+ gfxReal mCPSWeightFac;
+
+ int mDebugInit;
+
+
+protected:
+ // apply init transformations
+ void applyTrafos();
+
+ // helper function for init -> swap components everywhere
+ void swapCoords(int a,int b);
+ // helper function for init -> mirror time
+ void mirrorTime();
+
+ // helper, init given array
+ void initTimeArray(LbmFloat t, std::vector<ControlParticle> &parts);
+
+ bool checkPointInside(ntlTree *tree, ntlVec3Gfx org, gfxReal &distance);
+};
+
+
+
+#endif
+
ntlWorld *gpWorld = NULL;
+
// API
// reset elbeemSimulationSettings struct with defaults
return 0;
}
+// fluidsim end
+extern "C"
+int elbeemFree() {
+
+ return 0;
+}
+
// start fluidsim init
extern "C"
int elbeemAddDomain(elbeemSimulationSettings *settings) {
/* name of the mesh, mostly for debugging */
mesh->name = "[unnamed]";
+
+ /* fluid control settings */
+ mesh->cpsTimeStart = 0;
+ mesh->cpsTimeEnd = 0;
+ mesh->cpsQuality = 0;
+
+ mesh->channelSizeAttractforceStrength = 0;
+ mesh->channelAttractforceStrength = NULL;
+ mesh->channelSizeAttractforceRadius = 0;
+ mesh->channelAttractforceRadius = NULL;
+ mesh->channelSizeVelocityforceStrength = 0;
+ mesh->channelVelocityforceStrength = NULL;
+ mesh->channelSizeVelocityforceRadius = 0;
+ mesh->channelVelocityforceRadius = NULL;
}
int globalMeshCounter = 1;
// add mesh as fluidsim object
extern "C"
int elbeemAddMesh(elbeemMesh *mesh) {
- int initType = -1;
+ int initType;
if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddMesh","World and domain not initialized, call elbeemInit and elbeemAddDomain before...", SIMWORLD_INITERROR); }
switch(mesh->type) {
case OB_FLUIDSIM_FLUID: initType = FGI_FLUID; break;
case OB_FLUIDSIM_INFLOW: initType = FGI_MBNDINFLOW; break;
case OB_FLUIDSIM_OUTFLOW: initType = FGI_MBNDOUTFLOW; break;
+ case OB_FLUIDSIM_CONTROL: initType = FGI_CONTROL; break;
+ default: return 1; // invalid type
}
- // invalid type?
- if(initType<0) return 1;
ntlGeometryObjModel *obj = new ntlGeometryObjModel( );
gpWorld->getRenderGlobals()->getSimScene()->addGeoClass( obj );
obj->setGeoInitId( mesh->parentDomainId+1 );
obj->setGeoInitIntersect(true);
obj->setGeoInitType(initType);
- obj->setGeoPartSlipValue(mesh->obstaclePartslip);
+
+ // abuse partslip value for control fluid: reverse control keys or not
+ if(initType == FGI_CONTROL)
+ obj->setGeoPartSlipValue(mesh->obstacleType);
+ else
+ obj->setGeoPartSlipValue(mesh->obstaclePartslip);
+
obj->setGeoImpactFactor(mesh->obstacleImpactFactor);
+
+ /* fluid control features */
+ obj->setCpsTimeStart(mesh->cpsTimeStart);
+ obj->setCpsTimeEnd(mesh->cpsTimeEnd);
+ obj->setCpsQuality(mesh->cpsQuality);
+
if((mesh->volumeInitType<VOLUMEINIT_VOLUME)||(mesh->volumeInitType>VOLUMEINIT_BOTH)) mesh->volumeInitType = VOLUMEINIT_VOLUME;
obj->setVolumeInit(mesh->volumeInitType);
// use channel instead, obj->setInitialVelocity( ntlVec3Gfx(mesh->iniVelocity[0], mesh->iniVelocity[1], mesh->iniVelocity[2]) );
+
obj->initChannels(
mesh->channelSizeTranslation, mesh->channelTranslation,
mesh->channelSizeRotation, mesh->channelRotation,
mesh->channelSizeScale, mesh->channelScale,
mesh->channelSizeActive, mesh->channelActive,
- mesh->channelSizeInitialVel, mesh->channelInitialVel
+ mesh->channelSizeInitialVel, mesh->channelInitialVel,
+ mesh->channelSizeAttractforceStrength, mesh->channelAttractforceStrength,
+ mesh->channelSizeAttractforceRadius, mesh->channelAttractforceRadius,
+ mesh->channelSizeVelocityforceStrength, mesh->channelVelocityforceStrength,
+ mesh->channelSizeVelocityforceRadius, mesh->channelVelocityforceRadius
);
obj->setLocalCoordInivel( mesh->localInivelCoords );
if(getElbeemState() != SIMWORLD_STOP) {
// ok, we're done...
delete gpWorld;
+
gpWorld = NULL;
debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 );
} else {
+++ /dev/null
-/******************************************************************************
- *
- * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * All code distributed as part of El'Beem is covered by the version 2 of the
- * GNU General Public License. See the file COPYING for details.
- * Copyright 2003-2006 Nils Thuerey
- *
- * API header
- */
-#ifndef ELBEEM_API_H
-#define ELBEEM_API_H
-
-
-// simulation run callback function type (elbeemSimulationSettings->runsimCallback)
-// best use with FLUIDSIM_CBxxx defines below.
-// >parameters
-// return values: 0=continue, 1=stop, 2=abort
-// data pointer: user data pointer from elbeemSimulationSettings->runsimUserData
-// status integer: 1=running simulation, 2=new frame saved
-// frame integer: if status is 1, contains current frame number
-typedef int (*elbeemRunSimulationCallback)(void *data, int status, int frame);
-#define FLUIDSIM_CBRET_CONTINUE 0
-#define FLUIDSIM_CBRET_STOP 1
-#define FLUIDSIM_CBRET_ABORT 2
-#define FLUIDSIM_CBSTATUS_STEP 1
-#define FLUIDSIM_CBSTATUS_NEWFRAME 2
-
-
-// global settings for the simulation
-typedef struct elbeemSimulationSettings {
- /* version number */
- short version;
- /* id number of simulation domain, needed if more than a
- * single domain should be simulated */
- short domainId;
-
- /* geometrical extent */
- float geoStart[3], geoSize[3];
-
- /* resolutions */
- short resolutionxyz;
- short previewresxyz;
- /* size of the domain in real units (meters along largest resolution x,y,z extent) */
- float realsize;
-
- /* fluid properties */
- double viscosity;
- /* gravity strength */
- float gravity[3];
- /* anim start end time */
- float animStart, aniFrameTime;
- /* no. of frames to simulate & output */
- short noOfFrames;
- /* g star param (LBM compressibility) */
- float gstar;
- /* activate refinement? */
- short maxRefine;
- /* probability for surface particle generation (0.0=off) */
- float generateParticles;
- /* amount of tracer particles to generate (0=off) */
- int numTracerParticles;
-
- /* store output path, and file prefix for baked fluid surface */
- char outputPath[160+80];
-
- /* channel for frame time, visc & gravity animations */
- int channelSizeFrameTime;
- float *channelFrameTime;
- int channelSizeViscosity;
- float *channelViscosity;
- int channelSizeGravity;
- float *channelGravity; // vector
-
- /* boundary types and settings for domain walls */
- short domainobsType;
- float domainobsPartslip;
- /* generate speed vectors for vertices (e.g. for image based motion blur)*/
- short generateVertexVectors;
- /* strength of surface smoothing */
- float surfaceSmoothing;
- /* no. of surface subdivisions */
- int surfaceSubdivs;
-
- /* global transformation to apply to fluidsim mesh */
- float surfaceTrafo[4*4];
-
- /* development variables, testing for upcoming releases...*/
- float farFieldSize;
-
- /* callback function to notify calling program of performed simulation steps
- * or newly available frame data, if NULL it is ignored */
- elbeemRunSimulationCallback runsimCallback;
- /* pointer passed to runsimCallback for user data storage */
- void* runsimUserData;
-
-} elbeemSimulationSettings;
-
-
-// defines for elbeemMesh->type below
-#define OB_FLUIDSIM_FLUID 4
-#define OB_FLUIDSIM_OBSTACLE 8
-#define OB_FLUIDSIM_INFLOW 16
-#define OB_FLUIDSIM_OUTFLOW 32
-
-// defines for elbeemMesh->obstacleType below
-#define FLUIDSIM_OBSTACLE_NOSLIP 1
-#define FLUIDSIM_OBSTACLE_PARTSLIP 2
-#define FLUIDSIM_OBSTACLE_FREESLIP 3
-
-#define OB_VOLUMEINIT_VOLUME 1
-#define OB_VOLUMEINIT_SHELL 2
-#define OB_VOLUMEINIT_BOTH (OB_VOLUMEINIT_SHELL|OB_VOLUMEINIT_VOLUME)
-
-// a single mesh object
-typedef struct elbeemMesh {
- /* obstacle,fluid or inflow... */
- short type;
- /* id of simulation domain it belongs to */
- short parentDomainId;
-
- /* vertices */
- int numVertices;
- float *vertices; // = float[n][3];
- /* animated vertices */
- int channelSizeVertices;
- float *channelVertices; // = float[channelSizeVertices* (n*3+1) ];
-
- /* triangles */
- int numTriangles;
- int *triangles; // = int[][3];
-
- /* animation channels */
- int channelSizeTranslation;
- float *channelTranslation;
- int channelSizeRotation;
- float *channelRotation;
- int channelSizeScale;
- float *channelScale;
-
- /* active channel */
- int channelSizeActive;
- float *channelActive;
- /* initial velocity channel (e.g. for inflow) */
- int channelSizeInitialVel;
- float *channelInitialVel; // vector
- /* use initial velocity in object coordinates? (e.g. for rotation) */
- short localInivelCoords;
- /* boundary types and settings */
- short obstacleType;
- float obstaclePartslip;
- /* amount of force transfer from fluid to obj, 0=off, 1=normal */
- float obstacleImpactFactor;
- /* init volume, shell or both? use OB_VOLUMEINIT_xxx defines above */
- short volumeInitType;
-
- /* name of the mesh, mostly for debugging */
- const char *name;
-} elbeemMesh;
-
-// API functions
-
-#ifdef __cplusplus
-extern "C" {
-#endif // __cplusplus
-
-
-// reset elbeemSimulationSettings struct with defaults
-void elbeemResetSettings(struct elbeemSimulationSettings*);
-
-// start fluidsim init (returns !=0 upon failure)
-int elbeemInit(void);
-
-// start fluidsim init (returns !=0 upon failure)
-int elbeemAddDomain(struct elbeemSimulationSettings*);
-
-// get failure message during simulation or init
-// if an error occured (the string is copied into buffer,
-// max. length = 256 chars )
-void elbeemGetErrorString(char *buffer);
-
-// reset elbeemMesh struct with zeroes
-void elbeemResetMesh(struct elbeemMesh*);
-
-// add mesh as fluidsim object
-int elbeemAddMesh(struct elbeemMesh*);
-
-// do the actual simulation
-int elbeemSimulate(void);
-
-// continue a previously stopped simulation
-int elbeemContinueSimulation(void);
-
-
-// helper functions
-
-// simplify animation channels
-// returns if the channel and its size changed
-int elbeemSimplifyChannelFloat(float *channel, int *size);
-int elbeemSimplifyChannelVec3(float *channel, int *size);
-
-// helper functions implemented in utilities.cpp
-
-/* set elbeem debug output level (0=off to 10=full on) */
-void elbeemSetDebugLevel(int level);
-/* elbeem debug output function, prints if debug level >0 */
-void elbeemDebugOut(char *msg);
-
-/* estimate how much memory a given setup will require */
-double elbeemEstimateMemreq(int res,
- float sx, float sy, float sz,
- int refine, char *retstr);
-
-
-
-#ifdef __cplusplus
-}
-#endif // __cplusplus
-
-
-
-/******************************************************************************/
-// internal defines, do not use for initializing elbeemMesh
-// structs, for these use OB_xxx defines above
-
-/*! fluid geometry init types */
-#define FGI_FLAGSTART 16
-#define FGI_FLUID (1<<(FGI_FLAGSTART+ 0))
-#define FGI_NO_FLUID (1<<(FGI_FLAGSTART+ 1))
-#define FGI_BNDNO (1<<(FGI_FLAGSTART+ 2))
-#define FGI_BNDFREE (1<<(FGI_FLAGSTART+ 3))
-#define FGI_BNDPART (1<<(FGI_FLAGSTART+ 4))
-#define FGI_NO_BND (1<<(FGI_FLAGSTART+ 5))
-#define FGI_MBNDINFLOW (1<<(FGI_FLAGSTART+ 6))
-#define FGI_MBNDOUTFLOW (1<<(FGI_FLAGSTART+ 7))
-
-// all boundary types at once
-#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
-
-
-#endif // ELBEEM_API_H
--- /dev/null
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * All code distributed as part of El'Beem is covered by the version 2 of the
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2006 Nils Thuerey
+ *
+ * Control API header
+ */
+
+#include "elbeem.h"
+#include "elbeem_control.h"
+
+// add mesh as fluidsim object
+int elbeemControlAddSet(struct elbeemControl*) {
+
+ return 0;
+}
+
+int elbeemControlComputeMesh(struct elbeemMesh*) {
+
+
+ return 0;
+}
+
--- /dev/null
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * All code distributed as part of El'Beem is covered by the version 2 of the
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2006 Nils Thuerey
+ *
+ * Control API header
+ */
+#ifndef ELBEEMCONTROL_API_H
+#define ELBEEMCONTROL_API_H
+
+// a single control particle set
+typedef struct elbeemControl {
+ /* influence forces */
+ float influenceAttraction;
+ float *channelInfluenceAttraction;
+ float channelSizeInfluenceAttraction;
+
+ float influenceVelocity;
+ float *channelInfluenceVelocity;
+ float channelSizeInfluenceVelocity;
+
+ float influenceMaxdist;
+ float *channelInfluenceMaxdist;
+ float channelSizeInfluenceMaxdist;
+
+ /* influence force radii */
+ float radiusAttraction;
+ float *channelRadiusAttraction;
+ float channelSizeRadiusAttraction;
+
+ float radiusVelocity;
+ float *channelRadiusVelocity;
+ float channelSizeRadiusVelocity;
+
+ float radiusMindist;
+ float *channelRadiusMindist;
+ float channelSizeRadiusMindist;
+ float radiusMaxdist;
+ float *channelRadiusMaxdist;
+ float channelSizeRadiusMaxdist;
+
+ /* control particle positions/scale */
+ float offset[3];
+ float *channelOffset;
+ float channelSizeOffset;
+
+ float scale[3];
+ float *channelScale;
+ float channelSizeScale;
+
+} elbeemControl;
+
+
+// add mesh as fluidsim object
+int elbeemControlAddSet(struct elbeemControl*);
+
+// sample & track mesh control particles, TODO add return type...
+int elbeemControlComputeMesh(struct elbeemMesh*);
+
+#endif // ELBEEMCONTROL_API_H
#include <algorithm>
#include <stdio.h>
-#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
-#include <ieeefp.h>
-#endif
-
// just use default rounding for platforms where its not available
#ifndef round
#define round(x) (x)
--- /dev/null
+/******************************************************************************
+ *
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the
+// GNU General Public License. See the file COPYING for details.
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+ *
+ * Mean Value Mesh Coords class
+ *
+ *****************************************************************************/
+
+#include "mvmcoords.h"
+#include <algorithm>
+using std::vector;
+
+void MeanValueMeshCoords::clear()
+{
+ mVertices.resize(0);
+ mNumVerts = 0;
+}
+
+void MeanValueMeshCoords::calculateMVMCs(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle> &tris,
+ vector<ntlVec3Gfx> &points, gfxReal numweights)
+{
+ clear();
+ mvmTransferPoint tds;
+ int mem = 0;
+ int i = 0;
+
+ mNumVerts = (int)reference_vertices.size();
+
+ for (vector<ntlVec3Gfx>::iterator iter = points.begin(); iter != points.end(); ++iter, ++i) {
+ /*
+ if(i%(points.size()/10)==1) debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"Computing weights, points: "<<i<<"/"<<points.size(),5 );
+ */
+ tds.lastpos = *iter;
+ tds.weights.resize(0); // clear
+ computeWeights(reference_vertices, tris, tds, numweights);
+ mem += (int)tds.weights.size();
+ mVertices.push_back(tds);
+ }
+ int mbmem = mem * sizeof(mvmFloat) / (1024*1024);
+ debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"vertices:"<<mNumVerts<<" points:"<<points.size()<<" weights:"<<mem<<", wmem:"<<mbmem<<"MB ",7 );
+}
+
+// from: mean value coordinates for closed triangular meshes
+// attention: fails if a point is exactly (or very close) to a vertex
+void MeanValueMeshCoords::computeWeights(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle>& tris,
+ mvmTransferPoint& tds, gfxReal numweights)
+{
+ const bool mvmFullDebug=false;
+ //const ntlVec3Gfx cEPS = 1.0e-6;
+ const mvmFloat cEPS = 1.0e-14;
+
+ //mvmFloat d[3], s[3], phi[3],c[3];
+ ntlVec3d u[3],c,d,s,phi;
+ int indices[3];
+
+ for (int i = 0; i < (int)reference_vertices.size(); ++i) {
+ tds.weights.push_back(mvmIndexWeight(i, 0.0));
+ }
+
+ // for each triangle
+ //for (vector<ntlTriangle>::iterator iter = tris.begin(); iter != tris.end();) {
+ for(int t=0; t<(int)tris.size(); t++) {
+
+ for (int i = 0; i < 3; ++i) { //, ++iter) {
+ indices[i] = tris[t].getPoints()[i];
+ u[i] = vec2D(reference_vertices[ indices[i] ]-tds.lastpos);
+ d[i] = normalize(u[i]); //.normalize();
+ //assert(d[i] != 0.);
+ if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","t"<<t<<" i"<<indices[i] //<<" lp"<<tds.lastpos
+ <<" v"<<reference_vertices[indices[i]]<<" u"<<u[i]<<" ");
+ // on vertex!
+ //? if(d[i]<=0.) continue;
+ }
+ //for (int i = 0; i < 3; ++i) { errMsg("III"," "<<i <<" i"<<indices[i]<<reference_vertices[ indices[i] ] ); }
+
+ // arcsin is not needed, see paper
+ phi[0] = 2.*asin( (mvmFloat)(0.5* norm(u[1]-u[2]) ) );
+ phi[1] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[2]) ) );
+ phi[2] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[1]) ) );
+ mvmFloat h = (phi[0] + phi[1] + phi[2])*0.5;
+ if (M_PI-h < cEPS) {
+ if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","point on triangle");
+ tds.weights.resize(0);
+ tds.weights.push_back( mvmIndexWeight(indices[0], sin(phi[0])*d[1]*d[2]));
+ tds.weights.push_back( mvmIndexWeight(indices[1], sin(phi[1])*d[0]*d[2]));
+ tds.weights.push_back( mvmIndexWeight(indices[2], sin(phi[2])*d[1]*d[0]));
+ break;
+ }
+ mvmFloat sinh = 2.*sin(h);
+ c[0] = (sinh*sin(h-phi[0]))/(sin(phi[1])*sin(phi[2]))-1.;
+ c[1] = (sinh*sin(h-phi[1]))/(sin(phi[0])*sin(phi[2]))-1.;
+ c[2] = (sinh*sin(h-phi[2]))/(sin(phi[0])*sin(phi[1]))-1.;
+ if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","c="<<c<<" phi="<<phi<<" d="<<d);
+ //if (c[0] > 1. || c[0] < 0. || c[1] > 1. || c[1] < 0. || c[2] > 1. || c[2] < 0.) continue;
+
+ s[0] = sqrtf((float)(1.-c[0]*c[0]));
+ s[1] = sqrtf((float)(1.-c[1]*c[1]));
+ s[2] = sqrtf((float)(1.-c[2]*c[2]));
+
+ if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","s");
+ if (s[0] <= cEPS || s[1] <= cEPS || s[2] <= cEPS) {
+ //MSG("position lies outside the triangle on the same plane -> ignore it");
+ continue;
+ }
+ const mvmFloat u0x = u[0][0];
+ const mvmFloat u0y = u[0][1];
+ const mvmFloat u0z = u[0][2];
+ const mvmFloat u1x = u[1][0];
+ const mvmFloat u1y = u[1][1];
+ const mvmFloat u1z = u[1][2];
+ const mvmFloat u2x = u[2][0];
+ const mvmFloat u2y = u[2][1];
+ const mvmFloat u2z = u[2][2];
+ mvmFloat det = u0x*u1y*u2z - u0x*u1z*u2y + u0y*u1z*u2x - u0y*u1x*u2z + u0z*u1x*u2y - u0z*u1y*u2x;
+ //assert(det != 0.);
+ if (det < 0.) {
+ s[0] = -s[0];
+ s[1] = -s[1];
+ s[2] = -s[2];
+ }
+
+ tds.weights[indices[0]].weight += (phi[0]-c[1]*phi[2]-c[2]*phi[1])/(d[0]*sin(phi[1])*s[2]);
+ tds.weights[indices[1]].weight += (phi[1]-c[2]*phi[0]-c[0]*phi[2])/(d[1]*sin(phi[2])*s[0]);
+ tds.weights[indices[2]].weight += (phi[2]-c[0]*phi[1]-c[1]*phi[0])/(d[2]*sin(phi[0])*s[1]);
+ if(mvmFullDebug) { errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[0]<<" o"<<tds.weights[indices[0]].weight);
+ errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[1]<<" o"<<tds.weights[indices[1]].weight);
+ errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[2]<<" o"<<tds.weights[indices[2]].weight);
+ errMsg("MeanValueMeshCoords::computeWeights","\n\n\n"); }
+ }
+
+ //sort weights
+ if((numweights>0.)&& (numweights<1.) ) {
+ //if( ((int)tds.weights.size() > maxNumWeights) && (maxNumWeights > 0) ) {
+ int maxNumWeights = (int)(tds.weights.size()*numweights);
+ if(maxNumWeights<=0) maxNumWeights = 1;
+ std::sort(tds.weights.begin(), tds.weights.end(), std::greater<mvmIndexWeight>());
+ // only use maxNumWeights-th largest weights
+ tds.weights.resize(maxNumWeights);
+ }
+
+ // normalize weights
+ mvmFloat totalWeight = 0.;
+ for (vector<mvmIndexWeight>::const_iterator witer = tds.weights.begin();
+ witer != tds.weights.end(); ++witer) {
+ totalWeight += witer->weight;
+ }
+ mvmFloat invTotalWeight;
+ if (totalWeight == 0.) {
+ if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","totalWeight == 0");
+ invTotalWeight = 0.0;
+ } else {
+ invTotalWeight = 1.0/totalWeight;
+ }
+
+ for (vector<mvmIndexWeight>::iterator viter = tds.weights.begin();
+ viter != tds.weights.end(); ++viter) {
+ viter->weight *= invTotalWeight;
+ //assert(finite(viter->weight) != 0);
+ if(!finite(viter->weight)) viter->weight=0.;
+ }
+}
+
+void MeanValueMeshCoords::transfer(vector<ntlVec3Gfx> &vertices, vector<ntlVec3Gfx>& displacements)
+{
+ displacements.resize(0);
+
+ //debMsgStd("MeanValueMeshCoords::transfer",DM_MSG,"vertices:"<<mNumVerts<<" curr_verts:"<<vertices.size()<<" ",7 );
+ if((int)vertices.size() != mNumVerts) {
+ errMsg("MeanValueMeshCoords::transfer","Different no of verts: "<<vertices.size()<<" vs "<<mNumVerts);
+ return;
+ }
+
+ for (vector<mvmTransferPoint>::iterator titer = mVertices.begin(); titer != mVertices.end(); ++titer) {
+ mvmTransferPoint &tds = *titer;
+ ntlVec3Gfx newpos(0.0);
+
+ for (vector<mvmIndexWeight>::iterator witer = tds.weights.begin();
+ witer != tds.weights.end(); ++witer) {
+ newpos += vertices[witer->index] * witer->weight;
+ //errMsg("transfer","np"<<newpos<<" v"<<vertices[witer->index]<<" w"<< witer->weight);
+ }
+
+ displacements.push_back(newpos);
+ //displacements.push_back(newpos - tds.lastpos);
+ //tds.lastpos = newpos;
+ }
+}
+
--- /dev/null
+/******************************************************************************
+ *
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the
+// GNU General Public License. See the file COPYING for details.
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+ *
+ * Mean Value Mesh Coords class
+ *
+ *****************************************************************************/
+
+#ifndef MVMCOORDS_H
+#define MVMCOORDS_H
+
+#include "utilities.h"
+#include "ntl_ray.h"
+#include <vector>
+#define mvmFloat double
+
+#ifdef WIN32
+#ifndef FREE_WINDOWS
+#include "float.h"
+#define isnan(n) _isnan(n)
+#define finite _finite
+#endif
+#endif
+
+// weight and triangle index
+class mvmIndexWeight {
+ public:
+
+ mvmIndexWeight() : weight(0.0) {}
+
+ mvmIndexWeight(int const& i, mvmFloat const& w) :
+ weight(w), index(i) {}
+
+ // for sorting
+ bool operator> (mvmIndexWeight const& w) const { return this->weight > w.weight; }
+ bool operator< (mvmIndexWeight const& w) const { return this->weight < w.weight; }
+
+ mvmFloat weight;
+ int index;
+};
+
+// transfer point with weights
+class mvmTransferPoint {
+ public:
+ //! position of transfer point
+ ntlVec3Gfx lastpos;
+ //! triangle weights
+ std::vector<mvmIndexWeight> weights;
+};
+
+
+//! compute mvmcs
+class MeanValueMeshCoords {
+
+ public:
+
+ MeanValueMeshCoords() {}
+ ~MeanValueMeshCoords() {
+ clear();
+ }
+
+ void clear();
+
+ void calculateMVMCs(std::vector<ntlVec3Gfx> &reference_vertices,
+ std::vector<ntlTriangle> &tris, std::vector<ntlVec3Gfx> &points, gfxReal numweights);
+
+ void transfer(std::vector<ntlVec3Gfx> &vertices, std::vector<ntlVec3Gfx>& displacements);
+
+ protected:
+
+ void computeWeights(std::vector<ntlVec3Gfx> &reference_vertices,
+ std::vector<ntlTriangle> &tris, mvmTransferPoint& tds, gfxReal numweights);
+
+ std::vector<mvmTransferPoint> mVertices;
+ int mNumVerts;
+
+};
+
+#endif
+
//! Default destructor
virtual ~ntlGeometryClass() {
delete mpAttrs;
+ delete mpSwsAttrs;
};
//! Return type id
mCachedMovPoints(), mCachedMovNormals(),
mTriangleDivs1(), mTriangleDivs2(), mTriangleDivs3(),
mMovPntsInited(-100.0), mMaxMovPnt(-1),
- mcGeoActive(1.)
+ mcGeoActive(1.),
+ mCpsTimeStart(0.), mCpsTimeEnd(1.0), mCpsQuality(10.),
+ mcAttrFStr(0.),mcAttrFRad(0.), mcVelFStr(0.), mcVelFRad(0.)
{
};
/*****************************************************************************/
/* Init attributes etc. of this object */
/*****************************************************************************/
-#define GEOINIT_STRINGS 9
+#define GEOINIT_STRINGS 10
static const char *initStringStrs[GEOINIT_STRINGS] = {
"fluid",
"bnd_no","bnd_noslip",
"bnd_free","bnd_freeslip",
"bnd_part","bnd_partslip",
- "inflow", "outflow"
+ "inflow", "outflow", "control",
};
static int initStringTypes[GEOINIT_STRINGS] = {
FGI_FLUID,
FGI_BNDNO, FGI_BNDNO,
FGI_BNDFREE, FGI_BNDFREE,
FGI_BNDPART, FGI_BNDPART,
- FGI_MBNDINFLOW, FGI_MBNDOUTFLOW
+ FGI_MBNDINFLOW, FGI_MBNDOUTFLOW,
+ FGI_CONTROL
};
void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
{
void ntlGeometryObject::initChannels(
int nTrans, float *trans, int nRot, float *rot, int nScale, float *scale,
- int nAct, float *act, int nIvel, float *ivel
+ int nAct, float *act, int nIvel, float *ivel,
+ int nAttrFStr, float *attrFStr,
+ int nAttrFRad, float *attrFRad,
+ int nVelFStr, float *velFStr,
+ int nVelFRad, float *velFRad
) {
const bool debugInitc=true;
if(debugInitc) { debMsgStd("ntlGeometryObject::initChannels",DM_MSG,"nt:"<<nTrans<<" nr:"<<nRot<<" ns:"<<nScale, 10);
if((scale)&&(nScale>0)) { ADD_CHANNEL_VEC(mcScale, nScale, scale); }
if((act)&&(nAct>0)) { ADD_CHANNEL_FLOAT(mcGeoActive, nAct, act); }
if((ivel)&&(nIvel>0)) { ADD_CHANNEL_VEC(mcInitialVelocity, nIvel, ivel); }
+
+ /* fluid control channels */
+ if((attrFStr)&&(nAttrFStr>0)) { ADD_CHANNEL_FLOAT(mcAttrFStr, nAttrFStr, attrFStr); }
+ if((attrFRad)&&(nAttrFRad>0)) { ADD_CHANNEL_FLOAT(mcAttrFRad, nAttrFRad, attrFRad); }
+ if((velFStr)&&(nVelFStr>0)) { ADD_CHANNEL_FLOAT(mcVelFStr, nAct, velFStr); }
+ if((velFRad)&&(nVelFRad>0)) { ADD_CHANNEL_FLOAT(mcVelFRad, nVelFRad, velFRad); }
checkIsAnimated();
+
if(debugInitc) {
debMsgStd("ntlGeometryObject::initChannels",DM_MSG,getName()<<
" nt:"<<mcTrans.accessValues().size()<<" nr:"<<mcRot.accessValues().size()<<
}
}
- if( (this-getMeshAnimated())
+ if( (this->getMeshAnimated())
|| (mcTrans.accessValues().size()>1) // VALIDATE
|| (mcRot.accessValues().size()>1)
|| (mcScale.accessValues().size()>1)
/*! Set/get the local inivel coords flag */
inline bool getLocalCoordInivel() const { return mLocalCoordInivel; }
inline void setLocalCoordInivel(bool set) { mLocalCoordInivel=set; }
-
+
+ /****************************************/
+ /* fluid control features */
+ /****************************************/
+ /*! Set/get the particle control set attract force strength */
+ inline float getCpsTimeStart() const { return mCpsTimeStart; }
+ inline void setCpsTimeStart(float set) { mCpsTimeStart=set; }
+
+ /*! Set/get the particle control set attract force strength */
+ inline float getCpsTimeEnd() const { return mCpsTimeEnd; }
+ inline void setCpsTimeEnd(float set) { mCpsTimeEnd=set; }
+
+ /*! Set/get the particle control set quality */
+ inline float getCpsQuality() const { return mCpsQuality; }
+ inline void setCpsQuality(float set) { mCpsQuality=set; }
+
+ inline AnimChannel<float> getCpsAttrFStr() const { return mcAttrFStr; }
+ inline AnimChannel<float> getCpsAttrFRad() const { return mcAttrFRad; }
+ inline AnimChannel<float> getCpsVelFStr() const { return mcVelFStr; }
+ inline AnimChannel<float> getCpsVelFRad() const { return mcVelFRad; }
+
+ /****************************************/
+
/*! Init channels from float arrays (for elbeem API) */
void initChannels(
int nTrans, float *trans, int nRot, float *rot, int nScale, float *scale,
- int nAct, float *act, int nIvel, float *ivel
+ int nAct, float *act, int nIvel, float *ivel,
+ int nAttrFStr, float *attrFStr,
+ int nAttrFRad, float *attrFRad,
+ int nVelFStr, float *velFStr,
+ int nVelFRad, float *velFRad
);
/*! is the object animated? */
/*! animated channels for in/outflow on/off */
AnimChannel<float> mcGeoActive;
+
+ /* fluid control settings */
+ float mCpsTimeStart;
+ float mCpsTimeEnd;
+ float mCpsQuality;
+ AnimChannel<float> mcAttrFStr, mcAttrFRad, mcVelFStr, mcVelFRad;
public:
if(mpTree != NULL) delete mpTree;
// cleanup lists, only if this is the rendering cleanup scene
- if(mSceneDel) {
+ if(mSceneDel)
+ {
for (vector<ntlGeometryClass*>::iterator iter = mGeos.begin();
iter != mGeos.end(); iter++) {
//errMsg("ntlScene::~ntlScene","Deleting obj "<<(*iter)->getName() );
class ntlTree;
class ntlScene;
class ntlRenderGlobals;
+class ntlGeometryObject;
//! store data for an intersection of a ray and a triangle
// NOT YET USED
/* CONSTRUCTORS */
/*! Default constructor */
ntlScene( ntlRenderGlobals *glob, bool del=true );
- /*! Default destructor */
- ~ntlScene();
+ /*! Default destructor */
+ ~ntlScene();
/*! Add an object to the scene */
inline void addGeoClass(ntlGeometryClass *geo) {
{
delete mpGlob->getRenderScene();
delete mpGlob->getSimScene();
- delete mpGlob;
- delete mpLightList;
- delete mpPropList;
- delete mpSims;
+
+ delete mpGlob;
+
+
+ // these get assigned to mpGlob but not freed there
+ delete mpLightList;
+ delete mpPropList; // materials
+ delete mpSims;
+
#ifndef NOGUI
if(mpOpenGLRenderer) delete mpOpenGLRenderer;
#endif // NOGUI
ntlRenderGlobals::~ntlRenderGlobals() {
if(mpOpenGlAttr) delete mpOpenGlAttr;
if(mpBlenderAttr) delete mpBlenderAttr;
+
+
}
if(mpGiTree) delete mpGiTree;
if(mpElbeemSettings) delete mpElbeemSettings;
if(mpLbm) delete mpLbm;
- if(mpParam) delete mpParam;
+ if(mpParam) delete mpParam;
+ if(mpParts) delete mpParts;
debMsgStd("SimulationObject",DM_MSG,"El'Beem Done!\n",10);
}
#include "solver_relax.h"
#include "particletracer.h"
-#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
-#include <ieeefp.h>
-#endif
+
/*****************************************************************************/
//! coarse step functions
#endif
#endif
+#if LBM_INCLUDE_CONTROL==1
+#include "solver_control.h"
+#endif
+
#if LBM_INCLUDE_TESTSOLVERS==1
#include "solver_test.h"
#endif // LBM_INCLUDE_TESTSOLVERS==1
LbmFloat& debRAC(LbmFloat* s,int l);
# endif // FSGR_STRICT_DEBUG==1
+# if LBM_INCLUDE_CONTROL==1
+ LbmControlData *mpControl;
+
+ void initCpdata();
+ void handleCpdata();
+ void cpDebugDisplay(int dispset);
+# endif // LBM_INCLUDE_CONTROL==1
+
bool mUseTestdata;
# if LBM_INCLUDE_TESTSOLVERS==1
// test functions
void handleTestdata();
void set3dHeight(int ,int );
- void initCpdata();
- void handleCpdata();
- void cpDebugDisplay(int dispset);
-
int mMpNum,mMpIndex;
int mOrgSizeX;
LbmFloat mOrgStartX;
--- /dev/null
+/******************************************************************************
+ *
+ * El'Beem - the visual lattice boltzmann freesurface simulator
+ * All code distributed as part of El'Beem is covered by the version 2 of the
+ * GNU General Public License. See the file COPYING for details.
+ *
+ * Copyright 2003-2008 Nils Thuerey
+ *
+ * control extensions
+ *
+ *****************************************************************************/
+#include "solver_class.h"
+#include "solver_relax.h"
+#include "particletracer.h"
+
+#include "solver_control.h"
+
+#include "controlparticles.h"
+
+#include "elbeem.h"
+
+#include "ntl_geometrymodel.h"
+
+/******************************************************************************
+ * LbmControlData control set
+ *****************************************************************************/
+
+LbmControlSet::LbmControlSet() :
+ mCparts(NULL), mCpmotion(NULL), mContrPartFile(""), mCpmotionFile(""),
+ mcForceAtt(0.), mcForceVel(0.), mcForceMaxd(0.),
+ mcRadiusAtt(0.), mcRadiusVel(0.), mcRadiusMind(0.), mcRadiusMaxd(0.),
+ mcCpScale(1.), mcCpOffset(0.)
+{
+}
+LbmControlSet::~LbmControlSet() {
+ if(mCparts) delete mCparts;
+ if(mCpmotion) delete mCpmotion;
+}
+void LbmControlSet::initCparts() {
+ mCparts = new ControlParticles();
+ mCpmotion = new ControlParticles();
+}
+
+
+
+/******************************************************************************
+ * LbmControlData control
+ *****************************************************************************/
+
+LbmControlData::LbmControlData() :
+ mSetForceStrength(0.),
+ mCons(),
+ mCpUpdateInterval(8), // DG: was 16 --> causes problems (big sphere after some time), unstable
+ mCpOutfile(""),
+ mCpForces(), mCpKernel(), mMdKernel(),
+ mDiffVelCon(1.),
+ mDebugCpscale(0.),
+ mDebugVelScale(0.),
+ mDebugCompavScale(0.),
+ mDebugAttScale(0.),
+ mDebugMaxdScale(0.),
+ mDebugAvgVelScale(0.)
+{
+}
+
+LbmControlData::~LbmControlData()
+{
+ while (!mCons.empty()) {
+ delete mCons.back(); mCons.pop_back();
+ }
+}
+
+
+void LbmControlData::parseControldataAttrList(AttributeList *attr) {
+ // controlpart vars
+ mSetForceStrength = attr->readFloat("tforcestrength", mSetForceStrength,"LbmControlData", "mSetForceStrength", false);
+ //errMsg("tforcestrength set to "," "<<mSetForceStrength);
+ mCpUpdateInterval = attr->readInt("controlparticle_updateinterval", mCpUpdateInterval,"LbmControlData","mCpUpdateInterval", false);
+ // tracer output file
+ mCpOutfile = attr->readString("controlparticle_outfile",mCpOutfile,"LbmControlData","mCpOutfile", false);
+ if(getenv("ELBEEM_CPOUTFILE")) {
+ string outfile(getenv("ELBEEM_CPOUTFILE"));
+ mCpOutfile = outfile;
+ debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPOUTFILE to set mCpOutfile to "<<outfile<<","<<mCpOutfile,7);
+ }
+
+ for(int cpii=0; cpii<10; cpii++) {
+ string suffix("");
+ //if(cpii>0)
+ { suffix = string("0"); suffix[0]+=cpii; }
+ LbmControlSet *cset;
+ cset = new LbmControlSet();
+ cset->initCparts();
+
+ cset->mContrPartFile = attr->readString("controlparticle"+suffix+"_file",cset->mContrPartFile,"LbmControlData","cset->mContrPartFile", false);
+ if((cpii==0) && (getenv("ELBEEM_CPINFILE")) ) {
+ string infile(getenv("ELBEEM_CPINFILE"));
+ cset->mContrPartFile = infile;
+ debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPINFILE to set mContrPartFile to "<<infile<<","<<cset->mContrPartFile,7);
+ }
+
+ LbmFloat cpvort=0.;
+ cset->mcRadiusAtt = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusatt", 0., "LbmControlData","mcRadiusAtt" );
+ cset->mcRadiusVel = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
+ cset->mcRadiusVel = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
+ cset->mCparts->setRadiusAtt(cset->mcRadiusAtt.get(0.));
+ cset->mCparts->setRadiusVel(cset->mcRadiusVel.get(0.));
+
+ // WARNING currently only for first set
+ //if(cpii==0) {
+ cset->mcForceAtt = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_attraction", 0. , "LbmControlData","cset->mcForceAtt", false);
+ cset->mcForceVel = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_velocity", 0. , "LbmControlData","mcForceVel", false);
+ cset->mcForceMaxd = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_maxdist", 0. , "LbmControlData","mcForceMaxd", false);
+ cset->mCparts->setInfluenceAttraction(cset->mcForceAtt.get(0.) );
+ // warning - stores temprorarily, value converted to dt dep. factor
+ cset->mCparts->setInfluenceVelocity(cset->mcForceVel.get(0.) , 0.01 ); // dummy dt
+ cset->mCparts->setInfluenceMaxdist(cset->mcForceMaxd.get(0.) );
+ cpvort = attr->readFloat("controlparticle"+suffix+"_vorticity", cpvort, "LbmControlData","cpvort", false);
+ cset->mCparts->setInfluenceTangential(cpvort);
+
+ cset->mcRadiusMind = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmin", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMind", false);
+ cset->mcRadiusMaxd = attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmax", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMaxd", false);
+ cset->mCparts->setRadiusMinMaxd(cset->mcRadiusMind.get(0.));
+ cset->mCparts->setRadiusMaxd(cset->mcRadiusMaxd.get(0.));
+ //}
+
+ // now local...
+ //LbmVec cpOffset(0.), cpScale(1.);
+ LbmFloat cpTimescale = 1.;
+ string cpMirroring("");
+
+ //cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
+ //cset->mcCpScale = attr->readChannelVec3f("controlparticle"+suffix+"_scale", ntlVec3f(1.), "LbmControlData","mcCpScale", false);
+ cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
+ cset->mcCpScale = attr->readChannelVec3f("controlparticle"+suffix+"_scale", ntlVec3f(1.), "LbmControlData","mcCpScale", false);
+ cpTimescale = attr->readFloat("controlparticle"+suffix+"_timescale", cpTimescale, "LbmControlData","cpTimescale", false);
+ cpMirroring = attr->readString("controlparticle"+suffix+"_mirror", cpMirroring, "LbmControlData","cpMirroring", false);
+
+ LbmFloat cpsWidth = cset->mCparts->getCPSWith();
+ cpsWidth = attr->readFloat("controlparticle"+suffix+"_cpswidth", cpsWidth, "LbmControlData","cpsWidth", false);
+ LbmFloat cpsDt = cset->mCparts->getCPSTimestep();
+ cpsDt = attr->readFloat("controlparticle"+suffix+"_cpstimestep", cpsDt, "LbmControlData","cpsDt", false);
+ LbmFloat cpsTstart = cset->mCparts->getCPSTimeStart();
+ cpsTstart = attr->readFloat("controlparticle"+suffix+"_cpststart", cpsTstart, "LbmControlData","cpsTstart", false);
+ LbmFloat cpsTend = cset->mCparts->getCPSTimeEnd();
+ cpsTend = attr->readFloat("controlparticle"+suffix+"_cpstend", cpsTend, "LbmControlData","cpsTend", false);
+ LbmFloat cpsMvmfac = cset->mCparts->getCPSMvmWeightFac();
+ cpsMvmfac = attr->readFloat("controlparticle"+suffix+"_cpsmvmfac", cpsMvmfac, "LbmControlData","cpsMvmfac", false);
+ cset->mCparts->setCPSWith(cpsWidth);
+ cset->mCparts->setCPSTimestep(cpsDt);
+ cset->mCparts->setCPSTimeStart(cpsTstart);
+ cset->mCparts->setCPSTimeEnd(cpsTend);
+ cset->mCparts->setCPSMvmWeightFac(cpsMvmfac);
+
+ cset->mCparts->setOffset( vec2L(cset->mcCpOffset.get(0.)) );
+ cset->mCparts->setScale( vec2L(cset->mcCpScale.get(0.)) );
+ cset->mCparts->setInitTimeScale( cpTimescale );
+ cset->mCparts->setInitMirror( cpMirroring );
+
+ int mDebugInit = 0;
+ mDebugInit = attr->readInt("controlparticle"+suffix+"_debuginit", mDebugInit,"LbmControlData","mDebugInit", false);
+ cset->mCparts->setDebugInit(mDebugInit);
+
+ // motion particle settings
+ LbmVec mcpOffset(0.), mcpScale(1.);
+ LbmFloat mcpTimescale = 1.;
+ string mcpMirroring("");
+
+ cset->mCpmotionFile = attr->readString("cpmotion"+suffix+"_file",cset->mCpmotionFile,"LbmControlData","mCpmotionFile", false);
+ mcpTimescale = attr->readFloat("cpmotion"+suffix+"_timescale", mcpTimescale, "LbmControlData","mcpTimescale", false);
+ mcpMirroring = attr->readString("cpmotion"+suffix+"_mirror", mcpMirroring, "LbmControlData","mcpMirroring", false);
+ mcpOffset = vec2L( attr->readVec3d("cpmotion"+suffix+"_offset", vec2P(mcpOffset),"LbmControlData","cpOffset", false) );
+ mcpScale = vec2L( attr->readVec3d("cpmotion"+suffix+"_scale", vec2P(mcpScale), "LbmControlData","cpScale", false) );
+
+ cset->mCpmotion->setOffset( vec2L(mcpOffset) );
+ cset->mCpmotion->setScale( vec2L(mcpScale) );
+ cset->mCpmotion->setInitTimeScale( mcpTimescale );
+ cset->mCpmotion->setInitMirror( mcpMirroring );
+
+ if(cset->mContrPartFile.length()>1) {
+ errMsg("LbmControlData","Using control particle set "<<cpii<<" file:"<<cset->mContrPartFile<<" cpmfile:"<<cset->mCpmotionFile<<" mirr:'"<<cset->mCpmotion->getInitMirror()<<"' " );
+ mCons.push_back( cset );
+ } else {
+ delete cset;
+ }
+ }
+
+ // debug, testing - make sure theres at least an empty set
+ if(mCons.size()<1) {
+ mCons.push_back( new LbmControlSet() );
+ mCons[0]->initCparts();
+ }
+
+ // take from first set
+ for(int cpii=1; cpii<(int)mCons.size(); cpii++) {
+ mCons[cpii]->mCparts->setRadiusMinMaxd( mCons[0]->mCparts->getRadiusMinMaxd() );
+ mCons[cpii]->mCparts->setRadiusMaxd( mCons[0]->mCparts->getRadiusMaxd() );
+ mCons[cpii]->mCparts->setInfluenceAttraction( mCons[0]->mCparts->getInfluenceAttraction() );
+ mCons[cpii]->mCparts->setInfluenceTangential( mCons[0]->mCparts->getInfluenceTangential() );
+ mCons[cpii]->mCparts->setInfluenceVelocity( mCons[0]->mCparts->getInfluenceVelocity() , 0.01 ); // dummy dt
+ mCons[cpii]->mCparts->setInfluenceMaxdist( mCons[0]->mCparts->getInfluenceMaxdist() );
+ }
+
+ // invert for usage in relax macro
+ mDiffVelCon = 1.-attr->readFloat("cpdiffvelcon", mDiffVelCon, "LbmControlData","mDiffVelCon", false);
+
+ mDebugCpscale = attr->readFloat("cpdebug_cpscale", mDebugCpscale, "LbmControlData","mDebugCpscale", false);
+ mDebugMaxdScale = attr->readFloat("cpdebug_maxdscale", mDebugMaxdScale, "LbmControlData","mDebugMaxdScale", false);
+ mDebugAttScale = attr->readFloat("cpdebug_attscale", mDebugAttScale, "LbmControlData","mDebugAttScale", false);
+ mDebugVelScale = attr->readFloat("cpdebug_velscale", mDebugVelScale, "LbmControlData","mDebugVelScale", false);
+ mDebugCompavScale = attr->readFloat("cpdebug_compavscale", mDebugCompavScale, "LbmControlData","mDebugCompavScale", false);
+ mDebugAvgVelScale = attr->readFloat("cpdebug_avgvelsc", mDebugAvgVelScale, "LbmControlData","mDebugAvgVelScale", false);
+}
+
+
+void
+LbmFsgrSolver::initCpdata()
+{
+ // enable for cps via env. vars
+ //if( (getenv("ELBEEM_CPINFILE")) || (getenv("ELBEEM_CPOUTFILE")) ){ mUseTestdata=1; }
+
+
+ // manually switch on! if this is zero, nothing is done...
+ mpControl->mSetForceStrength = this->mTForceStrength = 1.;
+ mpControl->mCons.clear();
+
+ // init all control fluid objects
+ int numobjs = (int)(mpGiObjects->size());
+ for(int o=0; o<numobjs; o++) {
+ ntlGeometryObjModel *obj = (ntlGeometryObjModel *)(*mpGiObjects)[o];
+ if(obj->getGeoInitType() & FGI_CONTROL) {
+ // add new control set per object
+ LbmControlSet *cset;
+
+ cset = new LbmControlSet();
+ cset->initCparts();
+
+ // dont load any file
+ cset->mContrPartFile = string("");
+
+ cset->mcForceAtt = obj->getCpsAttrFStr();
+ cset->mcRadiusAtt = obj->getCpsAttrFRad();
+ cset->mcForceVel = obj->getCpsVelFStr();
+ cset->mcRadiusVel = obj->getCpsVelFRad();
+
+ cset->mCparts->setCPSTimeStart(obj->getCpsTimeStart());
+ cset->mCparts->setCPSTimeEnd(obj->getCpsTimeEnd());
+
+ if(obj->getCpsQuality() > LBM_EPSILON)
+ cset->mCparts->setCPSWith(1.0 / obj->getCpsQuality());
+
+ // this value can be left at 0.5:
+ cset->mCparts->setCPSMvmWeightFac(0.5);
+
+ mpControl->mCons.push_back( cset );
+ mpControl->mCons[mpControl->mCons.size()-1]->mCparts->initFromObject(obj);
+ }
+ }
+
+ // NT blender integration manual test setup
+ if(0) {
+ // manually switch on! if this is zero, nothing is done...
+ mpControl->mSetForceStrength = this->mTForceStrength = 1.;
+ mpControl->mCons.clear();
+
+ // add new set
+ LbmControlSet *cset;
+
+ cset = new LbmControlSet();
+ cset->initCparts();
+ // dont load any file
+ cset->mContrPartFile = string("");
+
+ // set radii for attraction & velocity forces
+ // set strength of the forces
+ // don't set directly! but use channels:
+ // mcForceAtt, mcForceVel, mcForceMaxd, mcRadiusAtt, mcRadiusVel, mcRadiusMind, mcRadiusMaxd etc.
+
+ // wrong: cset->mCparts->setInfluenceAttraction(1.15); cset->mCparts->setRadiusAtt(1.5);
+ // right, e.g., to init some constant values:
+ cset->mcForceAtt = AnimChannel<float>(0.2);
+ cset->mcRadiusAtt = AnimChannel<float>(0.75);
+ cset->mcForceVel = AnimChannel<float>(0.2);
+ cset->mcRadiusVel = AnimChannel<float>(0.75);
+
+ // this value can be left at 0.5:
+ cset->mCparts->setCPSMvmWeightFac(0.5);
+
+ mpControl->mCons.push_back( cset );
+
+ // instead of reading from file (cset->mContrPartFile), manually init some particles
+ mpControl->mCons[0]->mCparts->initBlenderTest();
+
+ // other values that might be interesting to change:
+ //cset->mCparts->setCPSTimestep(0.02);
+ //cset->mCparts->setCPSTimeStart(0.);
+ //cset->mCparts->setCPSTimeEnd(1.);
+
+ //mpControl->mDiffVelCon = 1.; // more rigid velocity control, 0 (default) allows more turbulence
+ }
+
+ // control particle -------------------------------------------------------------------------------------
+
+ // init cppf stage, use set 0!
+ if(mCppfStage>0) {
+ if(mpControl->mCpOutfile.length()<1) mpControl->mCpOutfile = string("cpout"); // use getOutFilename !?
+ char strbuf[100];
+ const char *cpFormat = "_d%dcppf%d";
+
+ // initial coarse stage, no input
+ if(mCppfStage==1) {
+ mpControl->mCons[0]->mContrPartFile = "";
+ } else {
+ // read from prev stage
+ snprintf(strbuf,100, cpFormat ,LBMDIM,mCppfStage-1);
+ mpControl->mCons[0]->mContrPartFile = mpControl->mCpOutfile;
+ mpControl->mCons[0]->mContrPartFile += strbuf;
+ mpControl->mCons[0]->mContrPartFile += ".cpart2";
+ }
+
+ snprintf(strbuf,100, cpFormat ,LBMDIM,mCppfStage);
+ mpControl->mCpOutfile += strbuf;
+ } // */
+
+ for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+ ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+
+ // now set with real dt
+ cparts->setInfluenceVelocity( mpControl->mCons[cpssi]->mcForceVel.get(0.), mLevel[mMaxRefine].timestep);
+ cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
+ cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
+ errMsg("LbmControlData","CppfStage "<<mCppfStage<<" in:"<<mpControl->mCons[cpssi]->mContrPartFile<<
+ " out:"<<mpControl->mCpOutfile<<" cl:"<< cparts->getCharLength() );
+
+ // control particle test init
+ if(mpControl->mCons[cpssi]->mCpmotionFile.length()>=1) cpmotion->initFromTextFile(mpControl->mCons[cpssi]->mCpmotionFile);
+ // not really necessary...
+ //? cparts->setFluidSpacing( mLevel[mMaxRefine].nodeSize ); // use grid coords!?
+ //? cparts->calculateKernelWeight();
+ //? debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - motion inited: "<<cparts->getSize() ,10);
+
+ // ensure both are on for env. var settings
+ // when no particles, but outfile enabled, initialize
+ const int lev = mMaxRefine;
+ if((mpParticles) && (mpControl->mCpOutfile.length()>=1) && (cpssi==0)) {
+ // check if auto num
+ if( (mpParticles->getNumInitialParticles()<=1) &&
+ (mpParticles->getNumParticles()<=1) ) { // initParticles done afterwards anyway
+ int tracers = 0;
+ const int workSet = mLevel[lev].setCurr;
+ FSGR_FORIJK_BOUNDS(lev) {
+ if(RFLAG(lev,i,j,k, workSet)&(CFFluid)) tracers++;
+ }
+ if(LBMDIM==3) tracers /= 8;
+ else tracers /= 4;
+ mpParticles->setNumInitialParticles(tracers);
+ mpParticles->setDumpTextFile(mpControl->mCpOutfile);
+ debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - set tracers #"<<tracers<<", actual #"<<mpParticles->getNumParticles() ,10);
+ }
+ if(mpParticles->getDumpTextInterval()<=0.) {
+ mpParticles->setDumpTextInterval(mLevel[lev].timestep * mLevel[lev].lSizex);
+ debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - dump delta t not set, using dti="<< mpParticles->getDumpTextInterval()<<", sim dt="<<mLevel[lev].timestep, 5 );
+ }
+ mpParticles->setDumpParts(true); // DEBUG? also dump as particle system
+ }
+
+ if(mpControl->mCons[cpssi]->mContrPartFile.length()>=1) cparts->initFromTextFile(mpControl->mCons[cpssi]->mContrPartFile);
+ cparts->setFluidSpacing( mLevel[lev].nodeSize ); // use grid coords!?
+ cparts->calculateKernelWeight();
+ debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles mCons"<<cpssi<<" - inited, parts:"<<cparts->getTotalSize()<<","<<cparts->getSize()<<" dt:"<<mpParam->getTimestep()<<" control time:"<<cparts->getControlTimStart()<<" to "<<cparts->getControlTimEnd() ,10);
+ } // cpssi
+
+ if(getenv("ELBEEM_CPINFILE")) {
+ this->mTForceStrength = 1.0;
+ }
+ this->mTForceStrength = mpControl->mSetForceStrength;
+ if(mpControl->mCpOutfile.length()>=1) mpParticles->setDumpTextFile(mpControl->mCpOutfile);
+
+ // control particle init end -------------------------------------------------------------------------------------
+
+ // make sure equiv to solver init
+ if(this->mTForceStrength>0.) { \
+ mpControl->mCpForces.resize( mMaxRefine+1 );
+ for(int lev = 0; lev<=mMaxRefine; lev++) {
+ LONGINT rcellSize = (mLevel[lev].lSizex*mLevel[lev].lSizey*mLevel[lev].lSizez);
+ debMsgStd("LbmFsgrSolver::initControl",DM_MSG,"mCpForces init, lev="<<lev<<" rcs:"<<(int)(rcellSize+4)<<","<<(rcellSize*sizeof(ControlForces)/(1024*1024)), 9 );
+ mpControl->mCpForces[lev].resize( (int)(rcellSize+4) );
+ //for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces.push_back( ControlForces() );
+ for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces[lev][i].resetForces();
+ }
+ } // on?
+
+ debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles #mCons "<<mpControl->mCons.size()<<" done", 6);
+}
+
+
+#define CPODEBUG 0
+//define CPINTER ((int)(mpControl->mCpUpdateInterval))
+
+#define KERN(x,y,z) mpControl->mCpKernel[ (((z)*cpkarWidth + (y))*cpkarWidth + (x)) ]
+#define MDKERN(x,y,z) mpControl->mMdKernel[ (((z)*mdkarWidth + (y))*mdkarWidth + (x)) ]
+
+#define BOUNDCHECK(x,low,high) ( ((x)<low) ? low : (((x)>high) ? high : (x) ) )
+#define BOUNDSKIP(x,low,high) ( ((x)<low) || ((x)>high) )
+
+void
+LbmFsgrSolver::handleCpdata()
+{
+ myTime_t cpstart = getTime();
+ int cpChecks=0;
+ int cpInfs=0;
+ //debMsgStd("ControlData::handleCpdata",DM_MSG,"called... "<<this->mTForceStrength,1);
+
+ // add cp influence
+ if((true) && (this->mTForceStrength>0.)) {
+ // ok continue...
+ } // on off
+ else {
+ return;
+ }
+
+ if((mpControl->mCpUpdateInterval<1) || (this->mStepCnt%mpControl->mCpUpdateInterval==0)) {
+ // do full reinit later on...
+ }
+ else if(this->mStepCnt>mpControl->mCpUpdateInterval) {
+ // only reinit new cells
+ // TODO !? remove loop dependance!?
+#define NOFORCEENTRY(lev, i,j,k) (LBMGET_FORCE(lev, i,j,k).maxDistance==CPF_MAXDINIT)
+ // interpolate missing
+ for(int lev=0; lev<=mMaxRefine; lev++) {
+ FSGR_FORIJK_BOUNDS(lev) {
+ if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFFluid|CFInter) )
+ //if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFInter) )
+ //if(0)
+ { // only check new inter? RFLAG?check
+ if(NOFORCEENTRY(lev, i,j,k)) {
+ //errMsg("CP","FE_MISSING at "<<PRINT_IJK<<" f"<<LBMGET_FORCE(lev, i,j,k).weightAtt<<" md"<<LBMGET_FORCE(lev, i,j,k).maxDistance );
+
+ LbmFloat nbs=0.;
+ ControlForces vals;
+ vals.resetForces(); vals.maxDistance = 0.;
+ for(int l=1; l<this->cDirNum; l++) {
+ int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+ //errMsg("CP","FE_MISSING check "<<PRINT_VEC(ni,nj,nk)<<" f"<<LBMGET_FORCE(lev, ni,nj,nk).weightAtt<<" md"<<LBMGET_FORCE(lev, ni,nj,nk).maxDistance );
+ if(!NOFORCEENTRY(lev, ni,nj,nk)) {
+ //? vals.weightAtt += LBMGET_FORCE(lev, ni,nj,nk).weightAtt;
+ //? vals.forceAtt += LBMGET_FORCE(lev, ni,nj,nk).forceAtt;
+ vals.maxDistance += LBMGET_FORCE(lev, ni,nj,nk).maxDistance;
+ vals.forceMaxd += LBMGET_FORCE(lev, ni,nj,nk).forceMaxd;
+ vals.weightVel += LBMGET_FORCE(lev, ni,nj,nk).weightVel;
+ vals.forceVel += LBMGET_FORCE(lev, ni,nj,nk).forceVel;
+ // ignore att/compAv/avgVel here for now
+ nbs += 1.;
+ }
+ }
+ if(nbs>0.) {
+ nbs = 1./nbs;
+ //? LBMGET_FORCE(lev, i,j,k).weightAtt = vals.weightAtt*nbs;
+ //? LBMGET_FORCE(lev, i,j,k).forceAtt = vals.forceAtt*nbs;
+ LBMGET_FORCE(lev, i,j,k).maxDistance = vals.maxDistance*nbs;
+ LBMGET_FORCE(lev, i,j,k).forceMaxd = vals.forceMaxd*nbs;
+ LBMGET_FORCE(lev, i,j,k).weightVel = vals.weightVel*nbs;
+ LBMGET_FORCE(lev, i,j,k).forceVel = vals.forceVel*nbs;
+ }
+ /*ControlForces *ff = &LBMGET_FORCE(lev, i,j,k); // DEBUG
+ errMsg("CP","FE_MISSING rec at "<<PRINT_IJK // DEBUG
+ <<" w:"<<ff->weightAtt<<" wa:" <<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )
+ <<" v:"<<ff->weightVel<<" wv:" <<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )
+ <<" v:"<<ff->maxDistance<<" wv:" <<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] ) ); // DEBUG */
+ // else errMsg("CP","FE_MISSING rec at "<<PRINT_IJK<<" failed!"); // DEBUG
+
+ }
+ }
+ }} // ijk, lev
+
+ // mStepCnt > mpControl->mCpUpdateInterval
+ return;
+ } else {
+ // nothing to do ...
+ return;
+ }
+
+ // reset
+ for(int lev=0; lev<=mMaxRefine; lev++) {
+ FSGR_FORIJK_BOUNDS(lev) { LBMGET_FORCE(lev,i,j,k).resetForces(); }
+ }
+ // do setup for coarsest level
+ const int coarseLev = 0;
+ const int fineLev = mMaxRefine;
+
+ // init for current time
+ for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+ LbmControlSet *cset = mpControl->mCons[cpssi];
+
+ cparts->setRadiusAtt(cset->mcRadiusAtt.get(mSimulationTime));
+ cparts->setRadiusVel(cset->mcRadiusVel.get(mSimulationTime));
+ cparts->setInfluenceAttraction(cset->mcForceAtt.get(mSimulationTime) );
+ cparts->setInfluenceMaxdist(cset->mcForceMaxd.get(mSimulationTime) );
+ cparts->setRadiusMinMaxd(cset->mcRadiusMind.get(mSimulationTime));
+ cparts->setRadiusMaxd(cset->mcRadiusMaxd.get(mSimulationTime));
+ cparts->calculateKernelWeight(); // always necessary!?
+ cparts->setOffset( vec2L(cset->mcCpOffset.get(mSimulationTime)) );
+ cparts->setScale( vec2L(cset->mcCpScale.get(mSimulationTime)) );
+
+ cparts->setInfluenceVelocity( cset->mcForceVel.get(mSimulationTime), mLevel[fineLev].timestep );
+ cparts->setLastOffset( vec2L(cset->mcCpOffset.get(mSimulationTime-mLevel[fineLev].timestep)) );
+ cparts->setLastScale( vec2L(cset->mcCpScale.get(mSimulationTime-mLevel[fineLev].timestep)) );
+
+ }
+
+ // check actual values
+ LbmFloat iatt = ABS(mpControl->mCons[0]->mCparts->getInfluenceAttraction());
+ LbmFloat ivel = mpControl->mCons[0]->mCparts->getInfluenceVelocity();
+ LbmFloat imaxd = mpControl->mCons[0]->mCparts->getInfluenceMaxdist();
+ //errMsg("FINCIT","iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
+ for(int cpssi=1; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ LbmFloat iatt2 = ABS(mpControl->mCons[cpssi]->mCparts->getInfluenceAttraction());
+ LbmFloat ivel2 = mpControl->mCons[cpssi]->mCparts->getInfluenceVelocity();
+ LbmFloat imaxd2 = mpControl->mCons[cpssi]->mCparts->getInfluenceMaxdist();
+
+ // we allow negative attraction force here!
+ if(iatt2 > iatt) iatt = iatt2;
+
+ if(ivel2 >ivel) ivel = ivel2;
+ if(imaxd2>imaxd) imaxd= imaxd2;
+ //errMsg("FINCIT"," "<<cpssi<<" iatt2="<<iatt2<<" ivel2="<<ivel2<<" imaxd2="<<imaxd<<" NEW "<<" iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
+ }
+
+ if(iatt==0. && ivel==0. && imaxd==0.) {
+ debMsgStd("ControlData::initControl",DM_MSG,"Skipped, all zero...",4);
+ return;
+ }
+ //iatt = mpControl->mCons[1]->mCparts->getInfluenceAttraction(); //ivel = mpControl->mCons[1]->mCparts->getInfluenceVelocity(); //imaxd = mpControl->mCons[1]->mCparts->getInfluenceMaxdist(); // TTTTTT
+
+ // do control setup
+ for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+ ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+
+ // TEST!?
+ bool radmod = false;
+ const LbmFloat minRadSize = mLevel[coarseLev].nodeSize * 1.5;
+ if((cparts->getRadiusAtt()>0.) && (cparts->getRadiusAtt()<minRadSize) && (!radmod) ) {
+ LbmFloat radfac = minRadSize / cparts->getRadiusAtt(); radmod=true;
+ debMsgStd("ControlData::initControl",DM_MSG,"Modified radii att, fac="<<radfac, 7);
+ cparts->setRadiusAtt(cparts->getRadiusAtt()*radfac);
+ cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
+ cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
+ cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
+ } else if((cparts->getRadiusVel()>0.) && (cparts->getRadiusVel()<minRadSize) && (!radmod) ) {
+ LbmFloat radfac = minRadSize / cparts->getRadiusVel();
+ debMsgStd("ControlData::initControl",DM_MSG,"Modified radii vel, fac="<<radfac, 7);
+ cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
+ cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
+ cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
+ } else if((cparts->getRadiusMaxd()>0.) && (cparts->getRadiusMaxd()<minRadSize) && (!radmod) ) {
+ LbmFloat radfac = minRadSize / cparts->getRadiusMaxd();
+ debMsgStd("ControlData::initControl",DM_MSG,"Modified radii maxd, fac="<<radfac, 7);
+ cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
+ cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
+ }
+ if(radmod) {
+ debMsgStd("ControlData::initControl",DM_MSG,"Modified radii: att="<<
+ cparts->getRadiusAtt()<<", vel=" << cparts->getRadiusVel()<<", maxd=" <<
+ cparts->getRadiusMaxd()<<", mind=" << cparts->getRadiusMinMaxd() ,5);
+ }
+
+ cpmotion->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), NULL );
+ cparts->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), cpmotion );
+ }
+
+ // do control...
+ for(int lev=0; lev<=mMaxRefine; lev++) {
+ LbmFloat levVolume = 1.;
+ LbmFloat levForceScale = 1.;
+ for(int ll=lev; ll<mMaxRefine; ll++) {
+ if(LBMDIM==3) levVolume *= 8.;
+ else levVolume *= 4.;
+ levForceScale *= 2.;
+ }
+ errMsg("LbmFsgrSolver::handleCpdata","levVolume="<<levVolume<<" levForceScale="<<levForceScale );
+ //todo: scale velocity, att by level timestep!?
+
+ for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+ // ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+
+ // if control set is not active skip it
+ if((cparts->getControlTimStart() > mSimulationTime) || (cparts->getControlTimEnd() < mLastSimTime))
+ {
+ continue;
+ }
+
+ const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
+ LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
+ LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
+ LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
+#if LBMDIM==2
+ gsz = gsx;
+#endif
+ LbmFloat goffx = mvGeoStart[0];
+ LbmFloat goffy = mvGeoStart[1];
+ LbmFloat goffz = mvGeoStart[2];
+
+ //const LbmFloat cpwIncFac = 2.0;
+ // max to two thirds of domain size
+ const int cpw = MIN( mLevel[lev].lSizex/3, MAX( (int)( cparts->getRadiusAtt() /gsx) +1 , 2) ); // normal kernel, att,vel
+ const int cpkarWidth = 2*cpw+1;
+ mpControl->mCpKernel.resize(cpkarWidth* cpkarWidth* cpkarWidth);
+ ControlParticle cpt; cpt.reset();
+ cpt.pos = LbmVec( (gsx*(LbmFloat)cpw)+goffx, (gsy*(LbmFloat)cpw)+goffy, (gsz*(LbmFloat)cpw)+goffz ); // optimize?
+ cpt.density = 0.5; cpt.densityWeight = 0.5;
+#if LBMDIM==3
+ for(int k= 0; k<cpkarWidth; ++k) {
+#else // LBMDIM==3
+ { int k = cpw;
+#endif
+ for(int j= 0; j<cpkarWidth; ++j)
+ for(int i= 0; i<cpkarWidth; ++i) {
+ KERN(i,j,k).resetForces();
+ //LbmFloat dx = i-cpw; LbmFloat dy = j-cpw; LbmFloat dz = k-cpw;
+ //LbmVec dv = ( LbmVec(dx,dy,dz) );
+ //LbmFloat dl = norm( dv ); //LbmVec dir = dv / dl;
+ LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz ); // optimize?
+ cparts->calculateCpInfluenceOpt( &cpt, pos, LbmVec(0,0,0), &KERN(i,j,k) ,1. );
+ /*if((CPODEBUG)&&(k==cpw)) errMsg("kern"," at "<<PRINT_IJK<<" pos"<<pos<<" cpp"<<cpt.pos
+ <<" wf:"<<KERN(i,j,k).weightAtt<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceAtt[0],KERN(i,j,k).forceAtt[1],KERN(i,j,k).forceAtt[2] )
+ <<" wf:"<<KERN(i,j,k).weightVel<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceVel[0],KERN(i,j,k).forceVel[1],KERN(i,j,k).forceVel[2] )
+ <<" wf:"<<KERN(i,j,k).maxDistance<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceMaxd[0],KERN(i,j,k).forceMaxd[1],KERN(i,j,k).forceMaxd[2] ) ); // */
+ KERN(i,j,k).weightAtt *= 2.;
+ KERN(i,j,k).forceAtt *= 2.;
+ //KERN(i,j,k).forceAtt[1] *= 2.; KERN(i,j,k).forceAtt[2] *= 2.;
+ KERN(i,j,k).weightVel *= 2.;
+ KERN(i,j,k).forceVel *= 2.;
+ //KERN(i,j,k).forceVel[1] *= 2.; KERN(i,j,k).forceVel[2] *= 2.;
+ }
+ }
+
+ if(CPODEBUG) errMsg("cpw"," = "<<cpw<<" f"<< cparts->getRadiusAtt()<<" gsx"<<gsx<<" kpw"<<cpkarWidth); // DEBUG
+ // first cp loop - add att and vel forces
+ for(int cppi=0; cppi<cparts->getSize(); cppi++) {
+ ControlParticle *cp = cparts->getParticle(cppi);
+ if(cp->influence<=0.) continue;
+ const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
+ const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
+ int cpk = (int)( (cp->pos[2]-goffz)/gsz );
+ /*if( ((LBMDIM==3)&&(BOUNDSKIP(cpk - cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
+ ((LBMDIM==3)&&(BOUNDSKIP(cpk + cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
+ BOUNDSKIP(cpj - cpwsm, 0, mLevel[lev].lSizey ) ||
+ BOUNDSKIP(cpj + cpwsm, 0, mLevel[lev].lSizey ) ||
+ BOUNDSKIP(cpi - cpwsm, 0, mLevel[lev].lSizex ) ||
+ BOUNDSKIP(cpi + cpwsm, 0, mLevel[lev].lSizex ) ) {
+ continue;
+ } // */
+ int is,ie,js,je,ks,ke;
+ ks = BOUNDCHECK(cpk - cpw, getForZMinBnd(), getForZMaxBnd(lev) );
+ ke = BOUNDCHECK(cpk + cpw, getForZMinBnd(), getForZMaxBnd(lev) );
+ js = BOUNDCHECK(cpj - cpw, 0, mLevel[lev].lSizey );
+ je = BOUNDCHECK(cpj + cpw, 0, mLevel[lev].lSizey );
+ is = BOUNDCHECK(cpi - cpw, 0, mLevel[lev].lSizex );
+ ie = BOUNDCHECK(cpi + cpw, 0, mLevel[lev].lSizex );
+ if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
+ if(CPODEBUG) errMsg("cppft","i"<<cppi<<" cpw"<<cpw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<" i:"<<is<<","<<ie<<" j:"<<js<<","<<je<<" k:"<<ks<<","<<ke<<" "); // DEBUG
+ cpInfs++;
+
+ for(int k= ks; k<ke; ++k) {
+ for(int j= js; j<je; ++j) {
+
+ CellFlagType *pflag = &RFLAG(lev,is,j,k, mLevel[lev].setCurr);
+ ControlForces *kk = &KERN( is-cpi+cpw, j-cpj+cpw, k-cpk+cpw);
+ ControlForces *ff = &LBMGET_FORCE(lev,is,j,k);
+ pflag--; kk--; ff--;
+
+ for(int i= is; i<ie; ++i) {
+ // first cp loop (att,vel)
+ pflag++; kk++; ff++;
+
+ //add weight for bnd cells
+ const LbmFloat pwforce = kk->weightAtt;
+ // control particle mod,
+ // dont add multiple CFFluid fsgr boundaries
+ if(lev==mMaxRefine) {
+ //if( ( ((*pflag)&(CFFluid )) && (lev==mMaxRefine) ) ||
+ //( ((*pflag)&(CFGrNorm)) && (lev <mMaxRefine) ) ) {
+ if((*pflag)&(CFFluid|CFUnused)) {
+ // check not fromcoarse?
+ cp->density += levVolume* kk->weightAtt; // old CFFluid
+ } else if( (*pflag) & (CFEmpty) ) {
+ cp->density -= levVolume* 0.5;
+ } else { //if( ((*pflag) & (CFBnd)) ) {
+ cp->density -= levVolume* 0.2; // penalty
+ }
+ } else {
+ //if((*pflag)&(CFGrNorm)) {
+ //cp->density += levVolume* kk->weightAtt; // old CFFluid
+ //}
+ }
+ //else if(!((*pflag) & (CFUnused)) ) { cp->density -= levVolume* 0.2; } // penalty
+
+ if( (*pflag) & (CFFluid|CFInter) ) // RFLAG_check
+ {
+
+ cpChecks++;
+ //const LbmFloat pwforce = kk->weightAtt;
+ LbmFloat pwvel = kk->weightVel;
+ if((pwforce==0.)&&(pwvel==0.)) { continue; }
+ ff->weightAtt += 1e-6; // for distance
+
+ if(pwforce>0.) {
+ ff->weightAtt += pwforce *cp->densityWeight *cp->influence;
+ ff->forceAtt += kk->forceAtt *levForceScale *cp->densityWeight *cp->influence;
+
+ // old fill handling here
+ const int workSet =mLevel[lev].setCurr;
+ LbmFloat ux=0., uy=0., uz=0.;
+ FORDF1{
+ const LbmFloat dfn = QCELL(lev, i,j,k, workSet, l);
+ ux += (this->dfDvecX[l]*dfn);
+ uy += (this->dfDvecY[l]*dfn);
+ uz += (this->dfDvecZ[l]*dfn);
+ }
+ // control particle mod
+ cp->avgVelWeight += levVolume*pwforce;
+ cp->avgVelAcc += LbmVec(ux,uy,uz) * levVolume*pwforce;
+ }
+
+ if(pwvel>0.) {
+ // TODO make switch? vel.influence depends on density weight...
+ // (reduced lowering with 0.75 factor)
+ pwvel *= cp->influence *(1.-0.75*cp->densityWeight);
+ // control particle mod
+ // todo use Omega instead!?
+ ff->forceVel += cp->vel*levVolume*pwvel * velLatticeScale; // levVolume?
+ ff->weightVel += levVolume*pwvel; // levVolume?
+ ff->compAv += cp->avgVel*levVolume*pwvel; // levVolume?
+ ff->compAvWeight += levVolume*pwvel; // levVolume?
+ }
+
+ if(CPODEBUG) errMsg("cppft","i"<<cppi<<" at "<<PRINT_IJK<<" kern:"<<
+ PRINT_VEC(i-cpi+cpw, j-cpj+cpw, k-cpk+cpw )
+ //<<" w:"<<ff->weightAtt<<" wa:"
+ //<<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )
+ //<<" v:"<<ff->weightVel<<" wv:"
+ //<<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )
+ //<<" v:"<<ff->maxDistance<<" wv:"
+ //<<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )
+ );
+ } // celltype
+
+ } // ijk
+ } // ijk
+ } // ijk
+ } // cpi, end first cp loop (att,vel)
+ debMsgStd("LbmFsgrSolver::handleCpdata",DM_MSG,"Force cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
+ } //cpssi
+ } // lev
+
+ // second loop
+ for(int lev=0; lev<=mMaxRefine; lev++) {
+ LbmFloat levVolume = 1.;
+ LbmFloat levForceScale = 1.;
+ for(int ll=lev; ll<mMaxRefine; ll++) {
+ if(LBMDIM==3) levVolume *= 8.;
+ else levVolume *= 4.;
+ levForceScale *= 2.;
+ }
+ // prepare maxd forces
+ for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+
+ // WARNING copied from above!
+ const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
+ LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
+ LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
+ LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
+#if LBMDIM==2
+ gsz = gsx;
+#endif
+ LbmFloat goffx = mvGeoStart[0];
+ LbmFloat goffy = mvGeoStart[1];
+ LbmFloat goffz = mvGeoStart[2];
+
+ //const LbmFloat cpwIncFac = 2.0;
+ const int mdw = MIN( mLevel[lev].lSizex/2, MAX( (int)( cparts->getRadiusMaxd() /gsx) +1 , 2) ); // wide kernel, md
+ const int mdkarWidth = 2*mdw+1;
+ mpControl->mMdKernel.resize(mdkarWidth* mdkarWidth* mdkarWidth);
+ ControlParticle cpt; cpt.reset();
+ cpt.density = 0.5; cpt.densityWeight = 0.5;
+ cpt.pos = LbmVec( (gsx*(LbmFloat)mdw)+goffx, (gsy*(LbmFloat)mdw)+goffy, (gsz*(LbmFloat)mdw)+goffz ); // optimize?
+#if LBMDIM==3
+ for(int k= 0; k<mdkarWidth; ++k) {
+#else // LBMDIM==3
+ { int k = mdw;
+#endif
+ for(int j= 0; j<mdkarWidth; ++j)
+ for(int i= 0; i<mdkarWidth; ++i) {
+ MDKERN(i,j,k).resetForces();
+ LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz ); // optimize?
+ cparts->calculateMaxdForce( &cpt, pos, &MDKERN(i,j,k) );
+ }
+ }
+
+ // second cpi loop, maxd forces
+ if(cparts->getInfluenceMaxdist()>0.) {
+ for(int cppi=0; cppi<cparts->getSize(); cppi++) {
+ ControlParticle *cp = cparts->getParticle(cppi);
+ if(cp->influence<=0.) continue;
+ const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
+ const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
+ int cpk = (int)( (cp->pos[2]-goffz)/gsz );
+
+ int is,ie,js,je,ks,ke;
+ ks = BOUNDCHECK(cpk - mdw, getForZMinBnd(), getForZMaxBnd(lev) );
+ ke = BOUNDCHECK(cpk + mdw, getForZMinBnd(), getForZMaxBnd(lev) );
+ js = BOUNDCHECK(cpj - mdw, 0, mLevel[lev].lSizey );
+ je = BOUNDCHECK(cpj + mdw, 0, mLevel[lev].lSizey );
+ is = BOUNDCHECK(cpi - mdw, 0, mLevel[lev].lSizex );
+ ie = BOUNDCHECK(cpi + mdw, 0, mLevel[lev].lSizex );
+ if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
+ if(CPODEBUG) errMsg("cppft","i"<<cppi<<" mdw"<<mdw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<" i:"<<is<<","<<ie<<" j:"<<js<<","<<je<<" k:"<<ks<<","<<ke<<" "); // DEBUG
+ cpInfs++;
+
+ for(int k= ks; k<ke; ++k)
+ for(int j= js; j<je; ++j) {
+ CellFlagType *pflag = &RFLAG(lev,is-1,j,k, mLevel[lev].setCurr);
+ for(int i= is; i<ie; ++i) {
+ // second cpi loop, maxd forces
+ pflag++;
+ if( (*pflag) & (CFFluid|CFInter) ) // RFLAG_check
+ {
+ cpChecks++;
+ ControlForces *ff = &LBMGET_FORCE(lev,i,j,k);
+ if(ff->weightAtt == 0.) {
+ ControlForces *kk = &MDKERN( i-cpi+mdw, j-cpj+mdw, k-cpk+mdw);
+ const LbmFloat pmdf = kk->maxDistance;
+ if((ff->maxDistance > pmdf) || (ff->maxDistance<0.))
+ ff->maxDistance = pmdf;
+ ff->forceMaxd = kk->forceMaxd;
+ // todo use Omega instead!?
+ ff->forceVel = cp->vel* velLatticeScale;
+ }
+ } // celltype
+ } } // ijk
+ } // cpi, md loop
+ } // maxd inf>0 */
+
+
+ debMsgStd("ControlData::initControl",DM_MSG,"Maxd cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
+ } //cpssi
+
+ // normalize, only done once for the whole array
+ mpControl->mCons[0]->mCparts->finishControl( mpControl->mCpForces[lev], iatt,ivel,imaxd );
+
+ } // lev loop
+
+ myTime_t cpend = getTime();
+ debMsgStd("ControlData::handleCpdata",DM_MSG,"Time for cpgrid generation:"<< getTimeString(cpend-cpstart)<<", checks:"<<cpChecks<<" infs:"<<cpInfs<<" " ,8);
+
+ // warning, may return before
+}
+
+#if LBM_USE_GUI==1
+
+#define USE_GLUTILITIES
+#include "../gui/gui_utilities.h"
+
+void LbmFsgrSolver::cpDebugDisplay(int dispset)
+{
+ for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+ ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+ //ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+ // display cp parts
+ const bool cpCubes = false;
+ const bool cpDots = true;
+ const bool cpCpdist = true;
+ const bool cpHideIna = true;
+ glShadeModel(GL_FLAT);
+ glDisable( GL_LIGHTING ); // dont light lines
+
+ // dot influence
+ if((mpControl->mDebugCpscale>0.) && cpDots) {
+ glPointSize(mpControl->mDebugCpscale * 8.);
+ glBegin(GL_POINTS);
+ for(int i=0; i<cparts->getSize(); i++) {
+ if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
+ ntlVec3Gfx org( vec2G(cparts->getParticle(i)->pos ) );
+ //LbmFloat halfsize = 0.5;
+ LbmFloat scale = cparts->getParticle(i)->densityWeight;
+ //glColor4f( scale,scale,scale,scale );
+ glColor4f( 0.,scale,0.,scale );
+ glVertex3f( org[0],org[1],org[2] );
+ //errMsg("lbmDebugDisplay","CP "<<i<<" at "<<org); // DEBUG
+ }
+ glEnd();
+ }
+
+ // cp positions
+ if((mpControl->mDebugCpscale>0.) && cpDots) {
+ glPointSize(mpControl->mDebugCpscale * 3.);
+ glBegin(GL_POINTS);
+ glColor3f( 0,1,0 );
+ }
+ for(int i=0; i<cparts->getSize(); i++) {
+ if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
+ ntlVec3Gfx org( vec2G(cparts->getParticle(i)->pos ) );
+ LbmFloat halfsize = 0.5;
+ LbmFloat scale = cparts->getRadiusAtt() * cparts->getParticle(i)->densityWeight;
+ if(cpCubes){ glLineWidth( 1 );
+ glColor3f( 1,1,1 );
+ ntlVec3Gfx s = org-(halfsize * (scale));
+ ntlVec3Gfx e = org+(halfsize * (scale));
+ drawCubeWire( s,e ); }
+ if((mpControl->mDebugCpscale>0.) && cpDots) {
+ glVertex3f( org[0],org[1],org[2] );
+ }
+ }
+ if(cpDots) glEnd();
+
+ if(mpControl->mDebugAvgVelScale>0.) {
+ const float scale = mpControl->mDebugAvgVelScale;
+
+ glColor3f( 1.0,1.0,1 );
+ glBegin(GL_LINES);
+ for(int i=0; i<cparts->getSize(); i++) {
+ if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
+ ntlVec3Gfx org( vec2G(cparts->getParticle(i)->pos ) );
+
+ //errMsg("CPAVGVEL","i"<<i<<" pos"<<org<<" av"<<cparts->getParticle(i)->avgVel);// DEBUG
+ float dx = cparts->getParticle(i)->avgVel[0];
+ float dy = cparts->getParticle(i)->avgVel[1];
+ float dz = cparts->getParticle(i)->avgVel[2];
+ dx *= scale; dy *= scale; dz *= scale;
+ glVertex3f( org[0],org[1],org[2] );
+ glVertex3f( org[0]+dx,org[1]+dy,org[2]+dz );
+ }
+ glEnd();
+ } // */
+
+ if( (LBMDIM==2) && (cpCpdist) ) {
+
+ // debug, for use of e.g. LBMGET_FORCE LbmControlData *mpControl = this;
+# define TESTGET_FORCE(lev,i,j,k) mpControl->mCpForces[lev][ ((k*mLevel[lev].lSizey)+j)*mLevel[lev].lSizex+i ]
+
+ glBegin(GL_LINES);
+ //const int lev=0;
+ for(int lev=0; lev<=mMaxRefine; lev++) {
+ FSGR_FORIJK_BOUNDS(lev) {
+ LbmVec pos = LbmVec(
+ ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex) * ((LbmFloat)i+0.5) + mvGeoStart[0],
+ ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey) * ((LbmFloat)j+0.5) + mvGeoStart[1],
+ ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez) * ((LbmFloat)k+0.5) + mvGeoStart[2] );
+ if(LBMDIM==2) pos[2] = ((mvGeoEnd[2]-mvGeoStart[2])*0.5 + mvGeoStart[2]);
+
+ if((mpControl->mDebugMaxdScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt<=0.) )
+ if(TESTGET_FORCE(lev,i,j,k).maxDistance>=0.)
+ if(TESTGET_FORCE(lev,i,j,k).maxDistance<CPF_MAXDINIT ) {
+ const float scale = mpControl->mDebugMaxdScale*10001.;
+ float dx = TESTGET_FORCE(lev,i,j,k).forceMaxd[0];
+ float dy = TESTGET_FORCE(lev,i,j,k).forceMaxd[1];
+ float dz = TESTGET_FORCE(lev,i,j,k).forceMaxd[2];
+ dx *= scale; dy *= scale; dz *= scale;
+ glColor3f( 0,1,0 );
+ glVertex3f( pos[0],pos[1],pos[2] );
+ glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+ } // */
+ if((mpControl->mDebugAttScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt>0.)) {
+ const float scale = mpControl->mDebugAttScale*100011.;
+ float dx = TESTGET_FORCE(lev,i,j,k).forceAtt[0];
+ float dy = TESTGET_FORCE(lev,i,j,k).forceAtt[1];
+ float dz = TESTGET_FORCE(lev,i,j,k).forceAtt[2];
+ dx *= scale; dy *= scale; dz *= scale;
+ glColor3f( 1,0,0 );
+ glVertex3f( pos[0],pos[1],pos[2] );
+ glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+ } // */
+ // why check maxDistance?
+ if((mpControl->mDebugVelScale>0.) && (TESTGET_FORCE(lev,i,j,k).maxDistance+TESTGET_FORCE(lev,i,j,k).weightVel>0.)) {
+ float scale = mpControl->mDebugVelScale*1.;
+ float wvscale = TESTGET_FORCE(lev,i,j,k).weightVel;
+ float dx = TESTGET_FORCE(lev,i,j,k).forceVel[0];
+ float dy = TESTGET_FORCE(lev,i,j,k).forceVel[1];
+ float dz = TESTGET_FORCE(lev,i,j,k).forceVel[2];
+ scale *= wvscale;
+ dx *= scale; dy *= scale; dz *= scale;
+ glColor3f( 0.2,0.2,1 );
+ glVertex3f( pos[0],pos[1],pos[2] );
+ glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+ } // */
+ if((mpControl->mDebugCompavScale>0.) && (TESTGET_FORCE(lev,i,j,k).compAvWeight>0.)) {
+ const float scale = mpControl->mDebugCompavScale*1.;
+ float dx = TESTGET_FORCE(lev,i,j,k).compAv[0];
+ float dy = TESTGET_FORCE(lev,i,j,k).compAv[1];
+ float dz = TESTGET_FORCE(lev,i,j,k).compAv[2];
+ dx *= scale; dy *= scale; dz *= scale;
+ glColor3f( 0.2,0.2,1 );
+ glVertex3f( pos[0],pos[1],pos[2] );
+ glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+ } // */
+ } // att,maxd
+ }
+ glEnd();
+ }
+ } // cpssi
+
+ //fprintf(stderr,"BLA\n");
+ glEnable( GL_LIGHTING ); // dont light lines
+ glShadeModel(GL_SMOOTH);
+}
+
+#else // LBM_USE_GUI==1
+void LbmFsgrSolver::cpDebugDisplay(int dispset) { }
+#endif // LBM_USE_GUI==1
+
+
--- /dev/null
+/******************************************************************************
+ *
+ * El'Beem - the visual lattice boltzmann freesurface simulator
+ * All code distributed as part of El'Beem is covered by the version 2 of the
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2006 Nils Thuerey
+ *
+ * testing extensions
+ *
+ *****************************************************************************/
+
+
+#ifndef LBM_TESTCLASS_H
+#define LBM_TESTCLASS_H
+
+//class IsoSurface;
+class ParticleObject;
+class ControlParticles;
+class ControlForces;
+
+//#define NUMGRIDS 2
+//#define MAXNUMSWS 10
+
+// farfield modes
+#define FARF_3DONLY -1
+#define FARF_BOTH 0
+#define FARF_SWEONLY 1
+// dont reuse 3d vars/init
+#define FARF_SEPSWE 2
+
+// relaxation macros for solver_relax.h
+#if LBM_INCLUDE_CONTROL!=1
+
+// defined in relax.h
+
+#else // LBM_INCLUDE_TESTSOLVERS!=1
+
+// WARNING has to match controlparts.h
+#define CPF_ENTRIES 12
+#define CPF_FORCE 0
+#define CPF_VELWEIGHT 3
+#define CPF_VELOCITY 4
+#define CPF_FORCEWEIGHT 7
+#define CPF_MINCPDIST 8
+#define CPF_MINCPDIR 9
+
+#include "controlparticles.h"
+
+#include "ntl_geometrymodel.h"
+
+// get force entry, set=0 is unused anyway
+#define LBMGET_FORCE(lev, i,j,k) mpControl->mCpForces[lev][ (LBMGI(lev,i,j,k,0)) ]
+
+// debug mods off...
+// same as in src/solver_relax.h!
+#define __PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
+ ux += (grav)[0]; \
+ uy += (grav)[1]; \
+ uz += (grav)[2];
+
+//void testMaxdmod(int i, int j,int k, LbmFloat &ux,LbmFloat &uy,LbmFloat &uz,ControlForces &ff);
+#if LBMDIM==3
+#define MAXDGRAV \
+ if(myforce->forceMaxd[0]*ux+myforce->forceMaxd[1]*uy<LBM_EPSILON) { \
+ ux = v2w*myforce->forceVel[0]+ v2wi*ux; \
+ uy = v2w*myforce->forceVel[1]+ v2wi*uy; } \
+ /* movement inverse to g? */ \
+ if((uz>LBM_EPSILON)&&(uz>myforce->forceVel[2])) { \
+ uz = v2w*myforce->forceVel[2]+ v2wi*uz; }
+#else // LBMDIM==3
+#define MAXDGRAV \
+ if(myforce->forceMaxd[0]*ux<LBM_EPSILON) { \
+ ux = v2w*myforce->forceVel[0]+ v2wi*ux; } \
+ /* movement inverse to g? */ \
+ if((uy>LBM_EPSILON)&&(uy>myforce->forceVel[1])) { \
+ uy = v2w*myforce->forceVel[1]+ v2wi*uy; }
+#endif // LBMDIM==3
+
+// debug modifications of collide vars (testing)
+// requires: lev,i,j,k
+#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
+ LbmFloat attforce = 1.; \
+ if(this->mTForceStrength>0.) { \
+ ControlForces* myforce = &LBMGET_FORCE(lev,i,j,k); \
+ const LbmFloat vf = myforce->weightAtt;\
+ const LbmFloat vw = myforce->weightVel;\
+ if(vf!=0.) { attforce = MAX(0., 1.-vf); /* TODO FIXME? use ABS(vf) for repulsion force? */ \
+ ux += myforce->forceAtt[0]; \
+ uy += myforce->forceAtt[1]; \
+ uz += myforce->forceAtt[2]; \
+ \
+ } else if(( myforce->maxDistance>0.) && ( myforce->maxDistance<CPF_MAXDINIT)) {\
+ const LbmFloat v2w = mpControl->mCons[0]->mCparts->getInfluenceMaxdist() * \
+ (myforce->maxDistance-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()) / (mpControl->mCons[0]->mCparts->getRadiusMaxd()-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()); \
+ const LbmFloat v2wi = 1.-v2w; \
+ if(v2w>0.){ MAXDGRAV; \
+ /* errMsg("ERRMDTT","at "<<PRINT_IJK<<" maxd="<<myforce->maxDistance<<", newu"<<PRINT_VEC(ux,uy,uz)<<", org"<<PRINT_VEC(oux,ouy,ouz)<<", fv"<<myforce->forceVel<<" " ); */ \
+ }\
+ } \
+ if(vw>0.) { \
+ const LbmFloat vwi = 1.-vw;\
+ const LbmFloat vwd = mpControl->mDiffVelCon;\
+ ux += vw*(myforce->forceVel[0]-myforce->compAv[0] + vwd*(myforce->compAv[0]-ux) ); \
+ uy += vw*(myforce->forceVel[1]-myforce->compAv[1] + vwd*(myforce->compAv[1]-uy) ); \
+ uz += vw*(myforce->forceVel[2]-myforce->compAv[2] + vwd*(myforce->compAv[2]-uz) ); \
+ /* TODO test!? modify smooth vel by influence of force for each lbm step, to account for force update only each N steps */ \
+ myforce->compAv = (myforce->forceVel*vw+ myforce->compAv*vwi); \
+ } \
+ } \
+ ux += (grav)[0]*attforce; \
+ uy += (grav)[1]*attforce; \
+ uz += (grav)[2]*attforce; \
+ /* end PRECOLLIDE_MODS */
+
+#define TEST_IF_CHECK \
+ if((!iffilled)&&(LBMGET_FORCE(lev,i,j,k).weightAtt!=0.)) { \
+ errMsg("TESTIFFILL"," at "<<PRINT_IJK<<" "<<mass<<" "<<rho); \
+ iffilled = true; \
+ if(mass<rho*1.0) mass = rho*1.0; myfrac = 1.0; \
+ }
+
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
+
+
+// a single set of control particles and params
+class LbmControlSet {
+ public:
+ LbmControlSet();
+ ~LbmControlSet();
+ void initCparts();
+
+ // control particles
+ ControlParticles *mCparts;
+ // control particle overall motion (for easier manual generation)
+ ControlParticles *mCpmotion;
+ // cp data file
+ string mContrPartFile;
+ string mCpmotionFile;
+ // cp debug displau
+ LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
+
+ // params
+ AnimChannel<float> mcForceAtt;
+ AnimChannel<float> mcForceVel;
+ AnimChannel<float> mcForceMaxd;
+
+ AnimChannel<float> mcRadiusAtt;
+ AnimChannel<float> mcRadiusVel;
+ AnimChannel<float> mcRadiusMind;
+ AnimChannel<float> mcRadiusMaxd;
+
+ AnimChannel<ntlVec3f> mcCpScale;
+ AnimChannel<ntlVec3f> mcCpOffset;
+};
+
+
+
+// main control data storage
+class LbmControlData
+{
+ public:
+ LbmControlData();
+ virtual ~LbmControlData();
+
+ // control data
+
+ // contorl params
+ void parseControldataAttrList(AttributeList *attr);
+
+ // control strength, set for solver interface
+ LbmFloat mSetForceStrength;
+ // cp vars
+ std::vector<LbmControlSet*> mCons;
+ // update interval
+ int mCpUpdateInterval;
+ // output
+ string mCpOutfile;
+ // control particle precomputed influence
+ std::vector< std::vector<ControlForces> > mCpForces;
+ std::vector<ControlForces> mCpKernel;
+ std::vector<ControlForces> mMdKernel;
+ // activate differential velcon
+ LbmFloat mDiffVelCon;
+
+ // cp debug displau
+ LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
+};
+
+#endif // LBM_TESTCLASS_H
mInit2dYZ(false),
mForceTadapRefine(-1), mCutoff(-1)
{
- // not much to do here...
+#if LBM_INCLUDE_CONTROL==1
+ mpControl = new LbmControlData();
+#endif
+
#if LBM_INCLUDE_TESTSOLVERS==1
mpTest = new LbmTestdata();
mMpNum = mMpIndex = 0;
delete mpIso;
if(mpPreviewSurface) delete mpPreviewSurface;
// cleanup done during scene deletion...
+
+#if LBM_INCLUDE_CONTROL==1
+ if(mpControl) delete mpControl;
+#endif
// always output performance estimate
debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
mSimulationTime += starttimeskip;
if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
+#if LBM_INCLUDE_CONTROL==1
+ mpControl->parseControldataAttrList(mpSifAttrs);
+#endif
+
#if LBM_INCLUDE_TESTSOLVERS==1
mUseTestdata = 0;
mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
// restrict max. chunk of 1 mem block to 1GB for windos
bool memBlockAllocProblem = false;
- double maxWinMemChunk = 1100.*1024.*1024.;
- double maxMacMemChunk = 1200.*1024.*1024.;
double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
//std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
#ifdef WIN32
+ double maxWinMemChunk = 1100.*1024.*1024.;
if(sizeof(void *)==4 && memEstFine>maxWinMemChunk) {
memBlockAllocProblem = true;
}
#endif // WIN32
#ifdef __APPLE__
+ double maxMacMemChunk = 1200.*1024.*1024.;
if(memEstFine> maxMacMemChunk) {
memBlockAllocProblem = true;
}
#endif // Mac
- if(sizeof(void *)==4 && memEstFine>maxDefaultMemChunk) {
+ if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
// max memory chunk for 32bit systems 2gig
memBlockAllocProblem = true;
}
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
mInitDone = 1;
-#if LBM_INCLUDE_TESTSOLVERS==1
+#if LBM_INCLUDE_CONTROL==1
initCpdata();
+#endif // LBM_INCLUDE_CONTROL==1
+
+#if LBM_INCLUDE_TESTSOLVERS==1
initTestdata();
#endif // ELBEEM_PLUGIN!=1
// not inited? dont use...
#include "loop_tools.h"
#include <stdlib.h>
-#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
-#include <ieeefp.h>
-#endif
-
-
/*****************************************************************************/
/*! perform a single LBM step */
/*****************************************************************************/
// init moving bc's, can change mMaxVlen
initMovingObstacles(false);
-#if LBM_INCLUDE_TESTSOLVERS==1
+#if LBM_INCLUDE_CONTROL==1
handleCpdata();
#endif
#define CAUSE_PANIC { this->mPanic=1; } /*set flag*/
#endif // FSGR_STRICT_DEBUG==1
-#if LBM_INCLUDE_TESTSOLVERS!=1
+// #if LBM_INCLUDE_TESTSOLVERS!=1
+#if LBM_INCLUDE_CONTROL!=1
#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
ux += (grav)[0]; \
#define TEST_IF_CHECK
-#else // LBM_INCLUDE_TESTSOLVERS!=1
-// defined in test.h
-
-#define NEWDIRVELMOTEST 0
-#if NEWDIRVELMOTEST==1
-// off for non testing
-#undef PRECOLLIDE_MODS
-#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
- ux += (grav)[0]; \
- uy += (grav)[1]; \
- uz += (grav)[2]; \
- { \
- int lev = mMaxRefine, nomb=0; \
- LbmFloat bcnt = 0.,nux=0.,nuy=0.,nuz=0.; \
- for(int l=1; l<this->cDfNum; l++) { \
- if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { \
- if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBndMoving) { \
- nux += QCELL_NB(lev, i,j,k,SRCS(lev),l, dMass); \
- nuy += QCELL_NB(lev, i,j,k,SRCS(lev),l, dFfrac); \
- bcnt += 1.; \
- } else { \
- nomb++; \
- } \
- } \
- } \
- if((bcnt>0.)&&(nomb==0)) { \
- ux = nux/bcnt; \
- uy = nuy/bcnt; \
- uz = nuz/bcnt; \
- } \
- }
-#else // NEWDIRVELMOTEST==1
-// off for non testing
-#endif // NEWDIRVELMOTEST==1
-
-#endif // LBM_INCLUDE_TESTSOLVERS!=1
+#else // LBM_INCLUDE_CONTROL!=1
+// defined in solver_control.h
+#endif // LBM_INCLUDE_CONTROL!=1
/******************************************************************************
DerivedMesh *mesh_create_derived_render(struct Object *ob,
CustomDataMask dataMask);
-/* same as above but wont use render settings */
+
+DerivedMesh *mesh_create_derived_index_render(struct Object *ob, CustomDataMask dataMask, int index);
+
+ /* same as above but wont use render settings */
DerivedMesh *mesh_create_derived_view(struct Object *ob,
CustomDataMask dataMask);
DerivedMesh *mesh_create_derived_no_deform(struct Object *ob,
--- /dev/null
+/**
+ * BKE_fluidsim.h
+ *
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "DNA_modifier_types.h"
+#include "DNA_object_fluidsim.h" // N_T
+#include "DNA_object_types.h"
+
+#include "BKE_DerivedMesh.h"
+
+/* old interface */
+FluidsimSettings *fluidsimSettingsNew(Object *srcob);
+
+void initElbeemMesh(Object *ob, int *numVertices, float **vertices, int *numTriangles, int **triangles, int useGlobalCoords, int modifierIndex);
+
+
+/* new fluid-modifier interface */
+void fluidsim_init(FluidsimModifierData *fluidmd);
+void fluidsim_free(FluidsimModifierData *fluidmd);
+
+DerivedMesh *fluidsim_read_cache(Object *ob, DerivedMesh *orgdm, FluidsimModifierData *fluidmd, int framenr, int useRenderParams);
+DerivedMesh *fluidsimModifier_do(FluidsimModifierData *fluidmd, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc);
+
+// get bounding box of mesh
+void fluid_get_bb(MVert *mvert, int totvert, float obmat[][4],
+ /*RET*/ float start[3], /*RET*/ float size[3] );
+
+
+
#include <config.h>
#endif
-#include <zlib.h>
-
#include "PIL_time.h"
#include "MEM_guardedalloc.h"
#include "BKE_deform.h"
#include "BKE_displist.h"
#include "BKE_effect.h"
+#include "BKE_fluidsim.h"
#include "BKE_global.h"
#include "BKE_key.h"
#include "BKE_material.h"
#include "GPU_extensions.h"
#include "GPU_material.h"
-// headers for fluidsim bobj meshes
-#include <stdlib.h>
-#include "LBM_fluidsim.h"
-#include "elbeem.h"
-
///////////////////////////////////
///////////////////////////////////
static DerivedMesh *getMeshDerivedMesh(Mesh *me, Object *ob, float (*vertCos)[3])
{
DerivedMesh *dm = CDDM_from_mesh(me, ob);
- int i, dofluidsim;
-
- dofluidsim = ((ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) &&
- (ob->fluidsimSettings->type & OB_FLUIDSIM_DOMAIN)&&
- (ob->fluidsimSettings->meshSurface) &&
- (me->totvert == ((Mesh *)(ob->fluidsimSettings->meshSurface))->totvert));
-
- if (vertCos && !dofluidsim)
+
+ if(!dm)
+ return NULL;
+
+ if (vertCos)
CDDM_apply_vert_coords(dm, vertCos);
CDDM_calc_normals(dm);
- /* apply fluidsim normals */
- if (dofluidsim) {
- // use normals from readBobjgz
- // TODO? check for modifiers!?
- MVert *fsvert = ob->fluidsimSettings->meshSurfNormals;
- short (*normals)[3] = MEM_mallocN(sizeof(short)*3*me->totvert, "fluidsim nor");
-
- for (i=0; i<me->totvert; i++) {
- VECCOPY(normals[i], fsvert[i].no);
- //mv->no[0]= 30000; mv->no[1]= mv->no[2]= 0; // DEBUG fixed test normals
- }
-
- CDDM_apply_vert_normals(dm, normals);
-
- MEM_freeN(normals);
- }
-
return dm;
}
static void mesh_calc_modifiers(Object *ob, float (*inputVertexCos)[3],
DerivedMesh **deform_r, DerivedMesh **final_r,
int useRenderParams, int useDeform,
- int needMapping, CustomDataMask dataMask)
+ int needMapping, CustomDataMask dataMask, int index)
{
Mesh *me = ob->data;
ModifierData *firstmd, *md;
float (*deformedVerts)[3] = NULL;
DerivedMesh *dm, *orcodm, *finaldm;
int numVerts = me->totvert;
- int fluidsimMeshUsed = 0;
int required_mode;
md = firstmd = modifiers_getVirtualModifierList(ob);
if(deform_r) *deform_r = NULL;
*final_r = NULL;
- /* replace original mesh by fluidsim surface mesh for fluidsim
- * domain objects
- */
- if((G.obedit!=ob) && !needMapping) {
- if((ob->fluidsimFlag & OB_FLUIDSIM_ENABLE)) {
- if(ob->fluidsimSettings->type & OB_FLUIDSIM_DOMAIN) {
- loadFluidsimMesh(ob,useRenderParams);
- fluidsimMeshUsed = 1;
- /* might have changed... */
- me = ob->data;
- numVerts = me->totvert;
- }
- }
- }
-
if(useRenderParams) required_mode = eModifierMode_Render;
else required_mode = eModifierMode_Realtime;
deformedVerts = mesh_getVertexCos(me, &numVerts);
/* Apply all leading deforming modifiers */
- for(; md; md = md->next, curr = curr->next) {
+ for(;md; md = md->next, curr = curr->next) {
ModifierTypeInfo *mti = modifierType_getInfo(md->type);
if((md->mode & required_mode) != required_mode) continue;
} else {
break;
}
+
+ /* grab modifiers until index i */
+ if((index >= 0) && (modifiers_indexInObject(ob, md) >= index))
+ break;
}
/* Result of all leading deforming modifiers is cached for
#endif
}
} else {
- if(!fluidsimMeshUsed) {
- /* default behaviour for meshes */
- if(inputVertexCos)
- deformedVerts = inputVertexCos;
- else
- deformedVerts = mesh_getRefKeyCos(me, &numVerts);
- } else {
- /* the fluid sim mesh might have more vertices than the original
- * one, so inputVertexCos shouldnt be used
- */
- deformedVerts = mesh_getVertexCos(me, &numVerts);
- }
+ /* default behaviour for meshes */
+ if(inputVertexCos)
+ deformedVerts = inputVertexCos;
+ else
+ deformedVerts = mesh_getRefKeyCos(me, &numVerts);
}
if(me->vnode) dm = derivedmesh_from_versemesh(me->vnode, deformedVerts);
#endif
- for(; md; md = md->next, curr = curr->next) {
+ for(;md; md = md->next, curr = curr->next) {
ModifierTypeInfo *mti = modifierType_getInfo(md->type);
if((md->mode & required_mode) != required_mode) continue;
}
}
}
+
+ /* grab modifiers until index i */
+ if((index >= 0) && (modifiers_indexInObject(ob, md) >= index))
+ break;
}
for(md=firstmd; md; md=md->next)
MEM_freeN(deformedVerts);
BLI_linklist_free(datamasks, NULL);
-
- /* restore mesh in any case */
- if(fluidsimMeshUsed) ob->data = ob->fluidsimSettings->orgMesh;
}
static float (*editmesh_getVertexCos(EditMesh *em, int *numVerts_r))[3]
mesh_calc_modifiers(ob, NULL, &ob->derivedDeform,
&ob->derivedFinal, 0, 1,
- needMapping, dataMask);
+ needMapping, dataMask, -1);
CustomData_free_layer_active(&me->fdata, CD_MCOL, me->totface);
} else {
mesh_calc_modifiers(ob, NULL, &ob->derivedDeform,
&ob->derivedFinal, G.rendering, 1,
- needMapping, dataMask);
+ needMapping, dataMask, -1);
}
INIT_MINMAX(min, max);
int orig_lvl= 0;
vert_copy= multires_render_pin(ob, me, &orig_lvl);
- mesh_calc_modifiers(ob, NULL, NULL, &final, 1, 1, 0, dataMask);
+ mesh_calc_modifiers(ob, NULL, NULL, &final, 1, 1, 0, dataMask, -1);
+ multires_render_final(ob, me, &final, vert_copy, orig_lvl, dataMask);
+
+ return final;
+}
+
+DerivedMesh *mesh_create_derived_index_render(Object *ob, CustomDataMask dataMask, int index)
+{
+ DerivedMesh *final;
+ Mesh *me= get_mesh(ob);
+ float *vert_copy= NULL;
+ int orig_lvl= 0;
+
+ vert_copy= multires_render_pin(ob, me, &orig_lvl);
+ mesh_calc_modifiers(ob, NULL, NULL, &final, 1, 1, 0, dataMask, index);
multires_render_final(ob, me, &final, vert_copy, orig_lvl, dataMask);
return final;
{
DerivedMesh *final;
- mesh_calc_modifiers(ob, NULL, NULL, &final, 0, 1, 0, dataMask);
+ mesh_calc_modifiers(ob, NULL, NULL, &final, 0, 1, 0, dataMask, -1);
return final;
}
{
DerivedMesh *final;
- mesh_calc_modifiers(ob, vertCos, NULL, &final, 0, 0, 0, dataMask);
+ mesh_calc_modifiers(ob, vertCos, NULL, &final, 0, 0, 0, dataMask, -1);
return final;
}
int orig_lvl= 0;
vert_copy= multires_render_pin(ob, me, &orig_lvl);
- mesh_calc_modifiers(ob, vertCos, NULL, &final, 1, 0, 0, dataMask);
+ mesh_calc_modifiers(ob, vertCos, NULL, &final, 1, 0, 0, dataMask, -1);
multires_render_final(ob, me, &final, vert_copy, orig_lvl, dataMask);
return final;
}
}
-/* ************************* fluidsim bobj file handling **************************** */
-
-#ifndef DISABLE_ELBEEM
-
-#ifdef WIN32
-#ifndef snprintf
-#define snprintf _snprintf
-#endif
-#endif
-
-/* write .bobj.gz file for a mesh object */
-void writeBobjgz(char *filename, struct Object *ob, int useGlobalCoords, int append, float time)
-{
- char debugStrBuffer[256];
- int wri=0,i,j,totvert,totface;
- float wrf;
- gzFile gzf;
- DerivedMesh *dm;
- float vec[3];
- float rotmat[3][3];
- MVert *mvert;
- MFace *mface;
- //if(append)return; // DEBUG
-
- if(!ob->data || (ob->type!=OB_MESH)) {
- snprintf(debugStrBuffer,256,"Writing GZ_BOBJ Invalid object %s ...\n", ob->id.name);
- elbeemDebugOut(debugStrBuffer);
- return;
- }
- if((ob->size[0]<0.0) || (ob->size[0]<0.0) || (ob->size[0]<0.0) ) {
- snprintf(debugStrBuffer,256,"\nfluidSim::writeBobjgz:: Warning object %s has negative scaling - check triangle ordering...?\n\n", ob->id.name);
- elbeemDebugOut(debugStrBuffer);
- }
-
- snprintf(debugStrBuffer,256,"Writing GZ_BOBJ '%s' ... ",filename); elbeemDebugOut(debugStrBuffer);
- if(append) gzf = gzopen(filename, "a+b9");
- else gzf = gzopen(filename, "wb9");
- if (!gzf) {
- snprintf(debugStrBuffer,256,"writeBobjgz::error - Unable to open file for writing '%s'\n", filename);
- elbeemDebugOut(debugStrBuffer);
- return;
- }
-
- dm = mesh_create_derived_render(ob, CD_MASK_BAREMESH);
- //dm = mesh_create_derived_no_deform(ob,NULL);
-
- mvert = dm->getVertArray(dm);
- mface = dm->getFaceArray(dm);
- totvert = dm->getNumVerts(dm);
- totface = dm->getNumFaces(dm);
-
- // write time value for appended anim mesh
- if(append) {
- gzwrite(gzf, &time, sizeof(time));
- }
-
- // continue with verts/norms
- if(sizeof(wri)!=4) { snprintf(debugStrBuffer,256,"Writing GZ_BOBJ, Invalid int size %d...\n", wri); elbeemDebugOut(debugStrBuffer); return; } // paranoia check
- wri = dm->getNumVerts(dm);
- mvert = dm->getVertArray(dm);
- gzwrite(gzf, &wri, sizeof(wri));
- for(i=0; i<wri;i++) {
- VECCOPY(vec, mvert[i].co);
- if(useGlobalCoords) { Mat4MulVecfl(ob->obmat, vec); }
- for(j=0; j<3; j++) {
- wrf = vec[j];
- gzwrite(gzf, &wrf, sizeof( wrf ));
- }
- }
-
- // should be the same as Vertices.size
- wri = totvert;
- gzwrite(gzf, &wri, sizeof(wri));
- EulToMat3(ob->rot, rotmat);
- for(i=0; i<wri;i++) {
- VECCOPY(vec, mvert[i].no);
- Normalize(vec);
- if(useGlobalCoords) { Mat3MulVecfl(rotmat, vec); }
- for(j=0; j<3; j++) {
- wrf = vec[j];
- gzwrite(gzf, &wrf, sizeof( wrf ));
- }
- }
-
- // append only writes verts&norms
- if(!append) {
- //float side1[3],side2[3],norm1[3],norm2[3];
- //float inpf;
-
- // compute no. of triangles
- wri = 0;
- for(i=0; i<totface; i++) {
- wri++;
- if(mface[i].v4) { wri++; }
- }
- gzwrite(gzf, &wri, sizeof(wri));
- for(i=0; i<totface; i++) {
-
- int face[4];
- face[0] = mface[i].v1;
- face[1] = mface[i].v2;
- face[2] = mface[i].v3;
- face[3] = mface[i].v4;
- //snprintf(debugStrBuffer,256,"F %s %d = %d,%d,%d,%d \n",ob->id.name, i, face[0],face[1],face[2],face[3] ); elbeemDebugOut(debugStrBuffer);
- //VecSubf(side1, mvert[face[1]].co,mvert[face[0]].co);
- //VecSubf(side2, mvert[face[2]].co,mvert[face[0]].co);
- //Crossf(norm1,side1,side2);
- gzwrite(gzf, &(face[0]), sizeof( face[0] ));
- gzwrite(gzf, &(face[1]), sizeof( face[1] ));
- gzwrite(gzf, &(face[2]), sizeof( face[2] ));
- if(face[3]) {
- //VecSubf(side1, mvert[face[2]].co,mvert[face[0]].co);
- //VecSubf(side2, mvert[face[3]].co,mvert[face[0]].co);
- //Crossf(norm2,side1,side2);
- //inpf = Inpf(norm1,norm2);
- //if(inpf>0.) {
- gzwrite(gzf, &(face[0]), sizeof( face[0] ));
- gzwrite(gzf, &(face[2]), sizeof( face[2] ));
- gzwrite(gzf, &(face[3]), sizeof( face[3] ));
- //} else {
- //gzwrite(gzf, &(face[0]), sizeof( face[0] ));
- //gzwrite(gzf, &(face[3]), sizeof( face[3] ));
- //gzwrite(gzf, &(face[2]), sizeof( face[2] ));
- //}
- } // quad
- }
- }
-
- snprintf(debugStrBuffer,256,"Done. #Vertices: %d, #Triangles: %d\n", totvert, totface );
- elbeemDebugOut(debugStrBuffer);
-
- gzclose( gzf );
- dm->release(dm);
-}
-
-void initElbeemMesh(struct Object *ob,
- int *numVertices, float **vertices,
- int *numTriangles, int **triangles,
- int useGlobalCoords)
-{
- DerivedMesh *dm = NULL;
- MVert *mvert;
- MFace *mface;
- int countTris=0, i, totvert, totface;
- float *verts;
- int *tris;
-
- dm = mesh_create_derived_render(ob, CD_MASK_BAREMESH);
- //dm = mesh_create_derived_no_deform(ob,NULL);
-
- mvert = dm->getVertArray(dm);
- mface = dm->getFaceArray(dm);
- totvert = dm->getNumVerts(dm);
- totface = dm->getNumFaces(dm);
-
- *numVertices = totvert;
- verts = MEM_callocN( totvert*3*sizeof(float), "elbeemmesh_vertices");
- for(i=0; i<totvert; i++) {
- VECCOPY( &verts[i*3], mvert[i].co);
- if(useGlobalCoords) { Mat4MulVecfl(ob->obmat, &verts[i*3]); }
- }
- *vertices = verts;
-
- for(i=0; i<totface; i++) {
- countTris++;
- if(mface[i].v4) { countTris++; }
- }
- *numTriangles = countTris;
- tris = MEM_callocN( countTris*3*sizeof(int), "elbeemmesh_triangles");
- countTris = 0;
- for(i=0; i<totface; i++) {
- int face[4];
- face[0] = mface[i].v1;
- face[1] = mface[i].v2;
- face[2] = mface[i].v3;
- face[3] = mface[i].v4;
-
- tris[countTris*3+0] = face[0];
- tris[countTris*3+1] = face[1];
- tris[countTris*3+2] = face[2];
- countTris++;
- if(face[3]) {
- tris[countTris*3+0] = face[0];
- tris[countTris*3+1] = face[2];
- tris[countTris*3+2] = face[3];
- countTris++;
- }
- }
- *triangles = tris;
-
- dm->release(dm);
-}
-
-/* read .bobj.gz file into a fluidsimDerivedMesh struct */
-Mesh* readBobjgz(char *filename, Mesh *orgmesh, float* bbstart, float *bbsize) //, fluidsimDerivedMesh *fsdm)
-{
- int wri,i,j;
- char debugStrBuffer[256];
- float wrf;
- Mesh *newmesh;
- const int debugBobjRead = 1;
- // init data from old mesh (materials,flags)
- MFace *origMFace = &((MFace*) orgmesh->mface)[0];
- int mat_nr = -1;
- int flag = -1;
- MFace *fsface = NULL;
- int gotBytes;
- gzFile gzf;
-
- if(!orgmesh) return NULL;
- if(!origMFace) return NULL;
- mat_nr = origMFace->mat_nr;
- flag = origMFace->flag;
-
- // similar to copy_mesh
- newmesh = MEM_dupallocN(orgmesh);
- newmesh->mat= orgmesh->mat;
-
- newmesh->mvert= NULL;
- newmesh->medge= NULL;
- newmesh->mface= NULL;
- newmesh->mtface= NULL;
-
- newmesh->dvert = NULL;
-
- newmesh->mcol= NULL;
- newmesh->msticky= NULL;
- newmesh->texcomesh= NULL;
- memset(&newmesh->vdata, 0, sizeof(newmesh->vdata));
- memset(&newmesh->edata, 0, sizeof(newmesh->edata));
- memset(&newmesh->fdata, 0, sizeof(newmesh->fdata));
-
- newmesh->key= NULL;
- newmesh->totface = 0;
- newmesh->totvert = 0;
- newmesh->totedge = 0;
- newmesh->medge = NULL;
-
-
- snprintf(debugStrBuffer,256,"Reading '%s' GZ_BOBJ... ",filename); elbeemDebugOut(debugStrBuffer);
- gzf = gzopen(filename, "rb");
- // gzf = fopen(filename, "rb");
- // debug: fread(b,c,1,a) = gzread(a,b,c)
- if (!gzf) {
- //snprintf(debugStrBuffer,256,"readBobjgz::error - Unable to open file for reading '%s'\n", filename); // DEBUG
- MEM_freeN(newmesh);
- return NULL;
- }
-
- //if(sizeof(wri)!=4) { snprintf(debugStrBuffer,256,"Reading GZ_BOBJ, Invalid int size %d...\n", wri); return NULL; } // paranoia check
- gotBytes = gzread(gzf, &wri, sizeof(wri));
- newmesh->totvert = wri;
- newmesh->mvert = CustomData_add_layer(&newmesh->vdata, CD_MVERT, CD_CALLOC, NULL, newmesh->totvert);
- if(debugBobjRead){ snprintf(debugStrBuffer,256,"#vertices %d ", newmesh->totvert); elbeemDebugOut(debugStrBuffer); } //DEBUG
- for(i=0; i<newmesh->totvert;i++) {
- //if(debugBobjRead) snprintf(debugStrBuffer,256,"V %d = ",i);
- for(j=0; j<3; j++) {
- gotBytes = gzread(gzf, &wrf, sizeof( wrf ));
- newmesh->mvert[i].co[j] = wrf;
- //if(debugBobjRead) snprintf(debugStrBuffer,256,"%25.20f ", wrf);
- }
- //if(debugBobjRead) snprintf(debugStrBuffer,256,"\n");
- }
-
- // should be the same as Vertices.size
- gotBytes = gzread(gzf, &wri, sizeof(wri));
- if(wri != newmesh->totvert) {
- // complain #vertices has to be equal to #normals, reset&abort
- CustomData_free_layer_active(&newmesh->vdata, CD_MVERT, newmesh->totvert);
- MEM_freeN(newmesh);
- snprintf(debugStrBuffer,256,"Reading GZ_BOBJ, #normals=%d, #vertices=%d, aborting...\n", wri,newmesh->totvert );
- return NULL;
- }
- for(i=0; i<newmesh->totvert;i++) {
- for(j=0; j<3; j++) {
- gotBytes = gzread(gzf, &wrf, sizeof( wrf ));
- newmesh->mvert[i].no[j] = (short)(wrf*32767.0f);
- //newmesh->mvert[i].no[j] = 0.5; // DEBUG tst
- }
- //fprintf(stderr," DEBDPCN nm%d, %d = %d,%d,%d \n",
- //(int)(newmesh->mvert), i, newmesh->mvert[i].no[0], newmesh->mvert[i].no[1], newmesh->mvert[i].no[2]);
- }
- //fprintf(stderr," DPCN 0 = %d,%d,%d \n", newmesh->mvert[0].no[0], newmesh->mvert[0].no[1], newmesh->mvert[0].no[2]);
-
-
- /* compute no. of triangles */
- gotBytes = gzread(gzf, &wri, sizeof(wri));
- newmesh->totface = wri;
- newmesh->mface = CustomData_add_layer(&newmesh->fdata, CD_MFACE, CD_CALLOC, NULL, newmesh->totface);
- if(debugBobjRead){ snprintf(debugStrBuffer,256,"#faces %d ", newmesh->totface); elbeemDebugOut(debugStrBuffer); } //DEBUG
- fsface = newmesh->mface;
- for(i=0; i<newmesh->totface; i++) {
- int face[4];
-
- gotBytes = gzread(gzf, &(face[0]), sizeof( face[0] ));
- gotBytes = gzread(gzf, &(face[1]), sizeof( face[1] ));
- gotBytes = gzread(gzf, &(face[2]), sizeof( face[2] ));
- face[3] = 0;
-
- fsface[i].v1 = face[0];
- fsface[i].v2 = face[1];
- fsface[i].v3 = face[2];
- fsface[i].v4 = face[3];
- }
-
- // correct triangles with v3==0 for blender, cycle verts
- for(i=0; i<newmesh->totface; i++) {
- if(!fsface[i].v3) {
- int temp = fsface[i].v1;
- fsface[i].v1 = fsface[i].v2;
- fsface[i].v2 = fsface[i].v3;
- fsface[i].v3 = temp;
- }
- }
-
- gzclose( gzf );
- for(i=0;i<newmesh->totface;i++) {
- fsface[i].mat_nr = mat_nr;
- fsface[i].flag = flag;
- fsface[i].edcode = ME_V1V2 | ME_V2V3 | ME_V3V1;
- //snprintf(debugStrBuffer,256,"%d : %d,%d,%d\n", i,fsface[i].mat_nr, fsface[i].flag, fsface[i].edcode );
- }
-
- snprintf(debugStrBuffer,256," (%d,%d) done\n", newmesh->totvert,newmesh->totface); elbeemDebugOut(debugStrBuffer); //DEBUG
- return newmesh;
-}
-
-/* read zipped fluidsim velocities into the co's of the fluidsimsettings normals struct */
-void readVelgz(char *filename, Object *srcob)
-{
- char debugStrBuffer[256];
- int wri, i, j;
- float wrf;
- gzFile gzf;
- MVert *vverts = srcob->fluidsimSettings->meshSurfNormals;
- int len = strlen(filename);
- Mesh *mesh = srcob->data;
- // mesh and vverts have to be valid from loading...
-
- // clean up in any case
- for(i=0; i<mesh->totvert;i++) {
- for(j=0; j<3; j++) {
- vverts[i].co[j] = 0.;
- }
- }
- if(srcob->fluidsimSettings->domainNovecgen>0) return;
-
- if(len<7) {
- //printf("readVelgz Eror: invalid filename '%s'\n",filename); // DEBUG
- return;
- }
-
- // .bobj.gz , correct filename
- // 87654321
- filename[len-6] = 'v';
- filename[len-5] = 'e';
- filename[len-4] = 'l';
-
- snprintf(debugStrBuffer,256,"Reading '%s' GZ_VEL... ",filename); elbeemDebugOut(debugStrBuffer);
- gzf = gzopen(filename, "rb");
- if (!gzf) {
- //printf("readVelgz Eror: unable to open file '%s'\n",filename); // DEBUG
- return;
- }
-
- gzread(gzf, &wri, sizeof( wri ));
- if(wri != mesh->totvert) {
- //printf("readVelgz Eror: invalid no. of velocities %d vs. %d aborting.\n" ,wri ,mesh->totvert ); // DEBUG
- return;
- }
-
- for(i=0; i<mesh->totvert;i++) {
- for(j=0; j<3; j++) {
- gzread(gzf, &wrf, sizeof( wrf ));
- vverts[i].co[j] = wrf;
- }
- //if(i<20) fprintf(stderr, "GZ_VELload %d = %f,%f,%f \n",i,vverts[i].co[0],vverts[i].co[1],vverts[i].co[2]); // DEBUG
- }
-
- gzclose(gzf);
-}
-
-
-/* ***************************** fluidsim derived mesh ***************************** */
-
-/* check which file to load, and replace old mesh of the object with it */
-/* this replacement is undone at the end of mesh_calc_modifiers */
-void loadFluidsimMesh(Object *srcob, int useRenderParams)
-{
- Mesh *mesh = NULL;
- float *bbStart = NULL, *bbSize = NULL;
- float lastBB[3];
- int displaymode = 0;
- int curFrame = G.scene->r.cfra - 1 /*G.scene->r.sfra*/; /* start with 0 at start frame */
- char targetDir[FILE_MAXFILE+FILE_MAXDIR], targetFile[FILE_MAXFILE+FILE_MAXDIR];
- char debugStrBuffer[256];
- //snprintf(debugStrBuffer,256,"loadFluidsimMesh call (obid '%s', rp %d)\n", srcob->id.name, useRenderParams); // debug
-
- if((!srcob)||(!srcob->fluidsimSettings)) {
- snprintf(debugStrBuffer,256,"DEBUG - Invalid loadFluidsimMesh call, rp %d, dm %d)\n", useRenderParams, displaymode); // debug
- elbeemDebugOut(debugStrBuffer); // debug
- return;
- }
- // make sure the original mesh data pointer is stored
- if(!srcob->fluidsimSettings->orgMesh) {
- srcob->fluidsimSettings->orgMesh = srcob->data;
- }
-
- // free old mesh, if there is one (todo, check if it's still valid?)
- if(srcob->fluidsimSettings->meshSurface) {
- Mesh *freeFsMesh = srcob->fluidsimSettings->meshSurface;
-
- // similar to free_mesh(...) , but no things like unlink...
- CustomData_free(&freeFsMesh->vdata, freeFsMesh->totvert);
- CustomData_free(&freeFsMesh->edata, freeFsMesh->totedge);
- CustomData_free(&freeFsMesh->fdata, freeFsMesh->totface);
- MEM_freeN(freeFsMesh);
-
- if(srcob->data == srcob->fluidsimSettings->meshSurface)
- srcob->data = srcob->fluidsimSettings->orgMesh;
- srcob->fluidsimSettings->meshSurface = NULL;
-
- if(srcob->fluidsimSettings->meshSurfNormals) MEM_freeN(srcob->fluidsimSettings->meshSurfNormals);
- srcob->fluidsimSettings->meshSurfNormals = NULL;
- }
-
- // init bounding box
- bbStart = srcob->fluidsimSettings->bbStart;
- bbSize = srcob->fluidsimSettings->bbSize;
- lastBB[0] = bbSize[0]; // TEST
- lastBB[1] = bbSize[1];
- lastBB[2] = bbSize[2];
- fluidsimGetAxisAlignedBB(srcob->fluidsimSettings->orgMesh, srcob->obmat, bbStart, bbSize, &srcob->fluidsimSettings->meshBB);
- // check free fsmesh... TODO
-
- if(!useRenderParams) {
- displaymode = srcob->fluidsimSettings->guiDisplayMode;
- } else {
- displaymode = srcob->fluidsimSettings->renderDisplayMode;
- }
-
- snprintf(debugStrBuffer,256,"loadFluidsimMesh call (obid '%s', rp %d, dm %d), curFra=%d, sFra=%d #=%d \n",
- srcob->id.name, useRenderParams, displaymode, G.scene->r.cfra, G.scene->r.sfra, curFrame ); // debug
- elbeemDebugOut(debugStrBuffer); // debug
-
- strncpy(targetDir, srcob->fluidsimSettings->surfdataPath, FILE_MAXDIR);
- // use preview or final mesh?
- if(displaymode==1) {
- // just display original object
- srcob->data = srcob->fluidsimSettings->orgMesh;
- return;
- } else if(displaymode==2) {
- strcat(targetDir,"fluidsurface_preview_####");
- } else { // 3
- strcat(targetDir,"fluidsurface_final_####");
- }
- BLI_convertstringcode(targetDir, G.sce);
- BLI_convertstringframe(targetDir, curFrame); // fixed #frame-no
-
- strcpy(targetFile,targetDir);
- strcat(targetFile, ".bobj.gz");
-
- snprintf(debugStrBuffer,256,"loadFluidsimMesh call (obid '%s', rp %d, dm %d) '%s' \n", srcob->id.name, useRenderParams, displaymode, targetFile); // debug
- elbeemDebugOut(debugStrBuffer); // debug
-
- if(displaymode!=2) { // dont add bounding box for final
- mesh = readBobjgz(targetFile, srcob->fluidsimSettings->orgMesh ,NULL,NULL);
- } else {
- mesh = readBobjgz(targetFile, srcob->fluidsimSettings->orgMesh, bbSize,bbSize );
- }
- if(!mesh) {
- // switch, abort background rendering when fluidsim mesh is missing
- const char *strEnvName2 = "BLENDER_ELBEEMBOBJABORT"; // from blendercall.cpp
- if(G.background==1) {
- if(getenv(strEnvName2)) {
- int elevel = atoi(getenv(strEnvName2));
- if(elevel>0) {
- printf("Env. var %s set, fluid sim mesh '%s' not found, aborting render...\n",strEnvName2, targetFile);
- exit(1);
- }
- }
- }
-
- // display org. object upon failure
- srcob->data = srcob->fluidsimSettings->orgMesh;
- return;
- }
-
- if((mesh)&&(mesh->totvert>0)) {
- make_edges(mesh, 0); // 0 = make all edges draw
- }
- srcob->fluidsimSettings->meshSurface = mesh;
- srcob->data = mesh;
- srcob->fluidsimSettings->meshSurfNormals = MEM_dupallocN(mesh->mvert);
-
- // load vertex velocities, if they exist...
- // TODO? use generate flag as loading flag as well?
- // warning, needs original .bobj.gz mesh loading filename
- if(displaymode==3) {
- readVelgz(targetFile, srcob);
- } else {
- // no data for preview, only clear...
- int i,j;
- for(i=0; i<mesh->totvert;i++) { for(j=0; j<3; j++) { srcob->fluidsimSettings->meshSurfNormals[i].co[j] = 0.; }}
- }
-
- //fprintf(stderr,"LOADFLM DEBXHCH fs=%d 3:%d,%d,%d \n", (int)mesh, ((Mesh *)(srcob->fluidsimSettings->meshSurface))->mvert[3].no[0], ((Mesh *)(srcob->fluidsimSettings->meshSurface))->mvert[3].no[1], ((Mesh *)(srcob->fluidsimSettings->meshSurface))->mvert[3].no[2]);
- return;
-}
-
-/* helper function */
-/* init axis aligned BB for mesh object */
-void fluidsimGetAxisAlignedBB(struct Mesh *mesh, float obmat[][4],
- /*RET*/ float start[3], /*RET*/ float size[3], /*RET*/ struct Mesh **bbmesh )
-{
- float bbsx=0.0, bbsy=0.0, bbsz=0.0;
- float bbex=1.0, bbey=1.0, bbez=1.0;
- int i;
- float vec[3];
-
- VECCOPY(vec, mesh->mvert[0].co);
- Mat4MulVecfl(obmat, vec);
- bbsx = vec[0]; bbsy = vec[1]; bbsz = vec[2];
- bbex = vec[0]; bbey = vec[1]; bbez = vec[2];
-
- for(i=1; i<mesh->totvert;i++) {
- VECCOPY(vec, mesh->mvert[i].co);
- Mat4MulVecfl(obmat, vec);
-
- if(vec[0] < bbsx){ bbsx= vec[0]; }
- if(vec[1] < bbsy){ bbsy= vec[1]; }
- if(vec[2] < bbsz){ bbsz= vec[2]; }
- if(vec[0] > bbex){ bbex= vec[0]; }
- if(vec[1] > bbey){ bbey= vec[1]; }
- if(vec[2] > bbez){ bbez= vec[2]; }
- }
-
- // return values...
- if(start) {
- start[0] = bbsx;
- start[1] = bbsy;
- start[2] = bbsz;
- }
- if(size) {
- size[0] = bbex-bbsx;
- size[1] = bbey-bbsy;
- size[2] = bbez-bbsz;
- }
-
- // init bounding box mesh?
- if(bbmesh) {
- int i,j;
- Mesh *newmesh = NULL;
- if(!(*bbmesh)) { newmesh = MEM_callocN(sizeof(Mesh), "fluidsimGetAxisAlignedBB_meshbb"); }
- else { newmesh = *bbmesh; }
-
- newmesh->totvert = 8;
- if(!newmesh->mvert)
- newmesh->mvert = CustomData_add_layer(&newmesh->vdata, CD_MVERT, CD_CALLOC, NULL, newmesh->totvert);
- for(i=0; i<8; i++) {
- for(j=0; j<3; j++) newmesh->mvert[i].co[j] = start[j];
- }
-
- newmesh->totface = 6;
- if(!newmesh->mface)
- newmesh->mface = CustomData_add_layer(&newmesh->fdata, CD_MFACE, CD_CALLOC, NULL, newmesh->totface);
-
- *bbmesh = newmesh;
- }
-}
-
-#else // DISABLE_ELBEEM
-
-/* dummy for mesh_calc_modifiers */
-void loadFluidsimMesh(Object *srcob, int useRenderParams) {
-}
-
-#endif // DISABLE_ELBEEM
-
dm->deformedOnly = 1;
- if(ob && ob->fluidsimSettings && ob->fluidsimSettings->meshSurface)
- alloctype= CD_DUPLICATE;
- else
- alloctype= CD_REFERENCE;
+ alloctype= CD_REFERENCE;
CustomData_merge(&mesh->vdata, &dm->vertData, CD_MASK_MESH, alloctype,
mesh->totvert);
ob->shapeflag &= ~OB_SHAPE_TEMPLOCK;
}
}
- if((ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && (ob->fluidsimSettings)) {
- // fluidsimSettings might not be initialized during load...
- if(ob->fluidsimSettings->type & (OB_FLUIDSIM_DOMAIN|OB_FLUIDSIM_PARTICLE)) {
- ob->recalc |= OB_RECALC_DATA; // NT FSPARTICLE
- }
- }
if(ob->particlesystem.first)
ob->recalc |= OB_RECALC_DATA;
break;
--- /dev/null
+/**
+ * fluidsim.c
+ *
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_object_force.h" // for pointcache
+#include "DNA_particle_types.h"
+#include "DNA_scene_types.h" // N_T
+
+#include "BLI_arithb.h"
+#include "BLI_blenlib.h"
+
+#include "BKE_cdderivedmesh.h"
+#include "BKE_customdata.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_fluidsim.h"
+#include "BKE_global.h"
+#include "BKE_modifier.h"
+#include "BKE_mesh.h"
+#include "BKE_pointcache.h"
+#include "BKE_utildefines.h"
+
+// headers for fluidsim bobj meshes
+#include <stdlib.h>
+#include "LBM_fluidsim.h"
+#include "elbeem.h"
+#include <zlib.h>
+#include <string.h>
+#include <stdio.h>
+
+/* ************************* fluidsim bobj file handling **************************** */
+
+#ifndef DISABLE_ELBEEM
+
+// -----------------------------------------
+// forward decleration
+// -----------------------------------------
+
+// -----------------------------------------
+
+void fluidsim_init(FluidsimModifierData *fluidmd)
+{
+ if(fluidmd)
+ {
+ FluidsimSettings *fss = MEM_callocN(sizeof(FluidsimSettings), "fluidsimsettings");
+
+ fluidmd->fss = fss;
+
+ if(!fss)
+ return;
+
+ fss->type = 0;
+ fss->show_advancedoptions = 0;
+
+ fss->resolutionxyz = 50;
+ fss->previewresxyz = 25;
+ fss->realsize = 0.03;
+ fss->guiDisplayMode = 2; // preview
+ fss->renderDisplayMode = 3; // render
+
+ fss->viscosityMode = 2; // default to water
+ fss->viscosityValue = 1.0;
+ fss->viscosityExponent = 6;
+
+ // dg TODO: change this to []
+ fss->gravx = 0.0;
+ fss->gravy = 0.0;
+ fss->gravz = -9.81;
+ fss->animStart = 0.0;
+ fss->animEnd = 0.30;
+ fss->gstar = 0.005; // used as normgstar
+ fss->maxRefine = -1;
+ // maxRefine is set according to resolutionxyz during bake
+
+ // fluid/inflow settings
+ // fss->iniVel --> automatically set to 0
+
+ /* elubie: changed this to default to the same dir as the render output
+ to prevent saving to C:\ on Windows */
+ BLI_strncpy(fss->surfdataPath, btempdir, FILE_MAX);
+
+ // first init of bounding box
+ // no bounding box needed
+
+ // todo - reuse default init from elbeem!
+ fss->typeFlags = 0;
+ fss->domainNovecgen = 0;
+ fss->volumeInitType = 1; // volume
+ fss->partSlipValue = 0.0;
+
+ fss->generateTracers = 0;
+ fss->generateParticles = 0.0;
+ fss->surfaceSmoothing = 1.0;
+ fss->surfaceSubdivs = 1.0;
+ fss->particleInfSize = 0.0;
+ fss->particleInfAlpha = 0.0;
+
+ // init fluid control settings
+ fss->attractforceStrength = 0.2;
+ fss->attractforceRadius = 0.75;
+ fss->velocityforceStrength = 0.2;
+ fss->velocityforceRadius = 0.75;
+ fss->cpsTimeStart = fss->animStart;
+ fss->cpsTimeEnd = fss->animEnd;
+ fss->cpsQuality = 10.0; // 1.0 / 10.0 => means 0.1 width
+
+ /*
+ BAD TODO: this is done in buttons_object.c in the moment
+ Mesh *mesh = ob->data;
+ // calculate bounding box
+ fluid_get_bb(mesh->mvert, mesh->totvert, ob->obmat, fss->bbStart, fss->bbSize);
+ */
+
+ fss->lastgoodframe = -1;
+
+ fss->flag = 0;
+
+ }
+
+ return;
+}
+
+void fluidsim_free(FluidsimModifierData *fluidmd)
+{
+ if(fluidmd)
+ {
+ MEM_freeN(fluidmd->fss);
+ }
+
+ return;
+}
+
+DerivedMesh *fluidsimModifier_do(FluidsimModifierData *fluidmd, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
+{
+ DerivedMesh *result = NULL;
+ int framenr;
+ FluidsimSettings *fss = NULL;
+
+ framenr= (int)G.scene->r.cfra;
+
+ // only handle fluidsim domains
+ if(fluidmd && fluidmd->fss && (fluidmd->fss->type != OB_FLUIDSIM_DOMAIN))
+ return dm;
+
+ // sanity check
+ if(!fluidmd || (fluidmd && !fluidmd->fss))
+ return dm;
+
+ fss = fluidmd->fss;
+
+ // timescale not supported yet
+ // clmd->sim_parms->timescale= timescale;
+
+ // support reversing of baked fluid frames here
+ if((fss->flag & OB_FLUIDSIM_REVERSE) && (fss->lastgoodframe >= 0))
+ {
+ framenr = fss->lastgoodframe - framenr + 1;
+ CLAMP(framenr, 1, fss->lastgoodframe);
+ }
+
+ /* try to read from cache */
+ if(((fss->lastgoodframe >= framenr) || (fss->lastgoodframe < 0)) && (result = fluidsim_read_cache(ob, dm, fluidmd, framenr, useRenderParams)))
+ {
+ // fss->lastgoodframe = framenr; // set also in src/fluidsim.c
+ return result;
+ }
+ else
+ {
+ // display last known good frame
+ if(fss->lastgoodframe >= 0)
+ {
+ if((result = fluidsim_read_cache(ob, dm, fluidmd, fss->lastgoodframe, useRenderParams)))
+ {
+ return result;
+ }
+
+ // it was supposed to be a valid frame but it isn't!
+ fss->lastgoodframe = framenr - 1;
+
+
+ // this could be likely the case when you load an old fluidsim
+ if((result = fluidsim_read_cache(ob, dm, fluidmd, fss->lastgoodframe, useRenderParams)))
+ {
+ return result;
+ }
+ }
+
+ result = CDDM_copy(dm);
+
+ if(result)
+ {
+ return result;
+ }
+ }
+
+ return dm;
+}
+
+/* read .bobj.gz file into a fluidsimDerivedMesh struct */
+static DerivedMesh *fluidsim_read_obj(char *filename)
+{
+ int wri,i,j;
+ float wrf;
+ int gotBytes;
+ gzFile gzf;
+ int numverts = 0, numfaces = 0;
+ DerivedMesh *dm = NULL;
+ MFace *mface;
+ MVert *mvert;
+ short *normals;
+
+ // ------------------------------------------------
+ // get numverts + numfaces first
+ // ------------------------------------------------
+ gzf = gzopen(filename, "rb");
+ if (!gzf)
+ {
+ return NULL;
+ }
+
+ // read numverts
+ gotBytes = gzread(gzf, &wri, sizeof(wri));
+ numverts = wri;
+
+ // skip verts
+ for(i=0; i<numverts*3; i++)
+ {
+ gotBytes = gzread(gzf, &wrf, sizeof( wrf ));
+ }
+
+ // read number of normals
+ gotBytes = gzread(gzf, &wri, sizeof(wri));
+
+ // skip normals
+ for(i=0; i<numverts*3; i++)
+ {
+ gotBytes = gzread(gzf, &wrf, sizeof( wrf ));
+ }
+
+ /* get no. of triangles */
+ gotBytes = gzread(gzf, &wri, sizeof(wri));
+ numfaces = wri;
+
+ gzclose( gzf );
+ // ------------------------------------------------
+
+ if(!numfaces || !numverts)
+ return NULL;
+
+ gzf = gzopen(filename, "rb");
+ if (!gzf)
+ {
+ return NULL;
+ }
+
+ dm = CDDM_new(numverts, 0, numfaces);
+
+ if(!dm)
+ {
+ gzclose( gzf );
+ return NULL;
+ }
+
+ // read numverts
+ gotBytes = gzread(gzf, &wri, sizeof(wri));
+
+ // read vertex position from file
+ mvert = CDDM_get_verts(dm);
+ for(i=0; i<numverts; i++)
+ {
+ MVert *mv = &mvert[i];
+
+ for(j=0; j<3; j++)
+ {
+ gotBytes = gzread(gzf, &wrf, sizeof( wrf ));
+ mv->co[j] = wrf;
+ }
+ }
+
+ // should be the same as numverts
+ gotBytes = gzread(gzf, &wri, sizeof(wri));
+ if(wri != numverts)
+ {
+ if(dm)
+ dm->release(dm);
+ gzclose( gzf );
+ return NULL;
+ }
+
+ normals = MEM_callocN(sizeof(short) * numverts * 3, "fluid_tmp_normals" );
+ if(!normals)
+ {
+ if(dm)
+ dm->release(dm);
+ gzclose( gzf );
+ return NULL;
+ }
+
+ // read normals from file (but don't save them yet)
+ for(i=0; i<numverts*3; i++)
+ {
+ gotBytes = gzread(gzf, &wrf, sizeof( wrf ));
+ normals[i] = (short)(wrf*32767.0f);
+ }
+
+ /* read no. of triangles */
+ gotBytes = gzread(gzf, &wri, sizeof(wri));
+
+ if(wri!=numfaces)
+ printf("Fluidsim: error in reading data from file.\n");
+
+ // read triangles from file
+ mface = CDDM_get_faces(dm);
+ for(i=0; i<numfaces; i++)
+ {
+ int face[4];
+ MFace *mf = &mface[i];
+
+ gotBytes = gzread(gzf, &(face[0]), sizeof( face[0] ));
+ gotBytes = gzread(gzf, &(face[1]), sizeof( face[1] ));
+ gotBytes = gzread(gzf, &(face[2]), sizeof( face[2] ));
+ face[3] = 0;
+
+ // check if 3rd vertex has index 0 (not allowed in blender)
+ if(face[2])
+ {
+ mf->v1 = face[0];
+ mf->v2 = face[1];
+ mf->v3 = face[2];
+ }
+ else
+ {
+ mf->v1 = face[1];
+ mf->v2 = face[2];
+ mf->v3 = face[0];
+ }
+ mf->v4 = face[3];
+
+ test_index_face(mf, NULL, 0, 3);
+ }
+
+ gzclose( gzf );
+
+ CDDM_calc_edges(dm);
+
+ CDDM_apply_vert_normals(dm, (short (*)[3])normals);
+ MEM_freeN(normals);
+
+ // CDDM_calc_normals(result);
+
+ return dm;
+}
+
+DerivedMesh *fluidsim_read_cache(Object *ob, DerivedMesh *orgdm, FluidsimModifierData *fluidmd, int framenr, int useRenderParams)
+{
+ int displaymode = 0;
+ int curFrame = framenr - 1 /*G.scene->r.sfra*/; /* start with 0 at start frame */
+ char targetDir[FILE_MAXFILE+FILE_MAXDIR], targetFile[FILE_MAXFILE+FILE_MAXDIR];
+ FluidsimSettings *fss = fluidmd->fss;
+ DerivedMesh *dm = NULL;
+ MFace *mface;
+ int numfaces;
+ int mat_nr, flag, i;
+
+ if(!useRenderParams) {
+ displaymode = fss->guiDisplayMode;
+ } else {
+ displaymode = fss->renderDisplayMode;
+ }
+
+ strncpy(targetDir, fss->surfdataPath, FILE_MAXDIR);
+
+ // use preview or final mesh?
+ if(displaymode==1)
+ {
+ // just display original object
+ return NULL;
+ }
+ else if(displaymode==2)
+ {
+ strcat(targetDir,"fluidsurface_preview_####");
+ }
+ else
+ { // 3
+ strcat(targetDir,"fluidsurface_final_####");
+ }
+
+ BLI_convertstringcode(targetDir, G.sce);
+ BLI_convertstringframe(targetDir, curFrame); // fixed #frame-no
+
+ strcpy(targetFile,targetDir);
+ strcat(targetFile, ".bobj.gz");
+
+ dm = fluidsim_read_obj(targetFile);
+
+ if(!dm)
+ {
+ // switch, abort background rendering when fluidsim mesh is missing
+ const char *strEnvName2 = "BLENDER_ELBEEMBOBJABORT"; // from blendercall.cpp
+
+ if(G.background==1) {
+ if(getenv(strEnvName2)) {
+ int elevel = atoi(getenv(strEnvName2));
+ if(elevel>0) {
+ printf("Env. var %s set, fluid sim mesh '%s' not found, aborting render...\n",strEnvName2, targetFile);
+ exit(1);
+ }
+ }
+ }
+
+ // display org. object upon failure which is in dm
+ return NULL;
+ }
+
+ // assign material + flags to new dm
+ mface = orgdm->getFaceArray(orgdm);
+ mat_nr = mface[0].mat_nr;
+ flag = mface[0].flag;
+
+ mface = dm->getFaceArray(dm);
+ numfaces = dm->getNumFaces(dm);
+ for(i=0; i<numfaces; i++)
+ {
+ mface[i].mat_nr = mat_nr;
+ mface[i].flag = flag;
+ }
+
+ // load vertex velocities, if they exist...
+ // TODO? use generate flag as loading flag as well?
+ // warning, needs original .bobj.gz mesh loading filename
+ /*
+ if(displaymode==3)
+ {
+ readVelgz(targetFile, srcob);
+ }
+ else
+ {
+ // no data for preview, only clear...
+ int i,j;
+ for(i=0; i<mesh->totvert;i++) { for(j=0; j<3; j++) { srcob->fluidsimSettings->meshSurfNormals[i].co[j] = 0.; }}
+ }*/
+
+ return dm;
+}
+
+void fluid_get_bb(MVert *mvert, int totvert, float obmat[][4],
+ /*RET*/ float start[3], /*RET*/ float size[3] )
+{
+ float bbsx=0.0, bbsy=0.0, bbsz=0.0;
+ float bbex=1.0, bbey=1.0, bbez=1.0;
+ int i;
+ float vec[3];
+
+ VECCOPY(vec, mvert[0].co);
+ Mat4MulVecfl(obmat, vec);
+ bbsx = vec[0]; bbsy = vec[1]; bbsz = vec[2];
+ bbex = vec[0]; bbey = vec[1]; bbez = vec[2];
+
+ for(i = 1; i < totvert; i++) {
+ VECCOPY(vec, mvert[i].co);
+ Mat4MulVecfl(obmat, vec);
+
+ if(vec[0] < bbsx){ bbsx= vec[0]; }
+ if(vec[1] < bbsy){ bbsy= vec[1]; }
+ if(vec[2] < bbsz){ bbsz= vec[2]; }
+ if(vec[0] > bbex){ bbex= vec[0]; }
+ if(vec[1] > bbey){ bbey= vec[1]; }
+ if(vec[2] > bbez){ bbez= vec[2]; }
+ }
+
+ // return values...
+ if(start) {
+ start[0] = bbsx;
+ start[1] = bbsy;
+ start[2] = bbsz;
+ }
+ if(size) {
+ size[0] = bbex-bbsx;
+ size[1] = bbey-bbsy;
+ size[2] = bbez-bbsz;
+ }
+}
+
+//-------------------------------------------------------------------------------
+// old interface
+//-------------------------------------------------------------------------------
+
+
+
+//-------------------------------------------------------------------------------
+// file handling
+//-------------------------------------------------------------------------------
+
+void initElbeemMesh(struct Object *ob,
+ int *numVertices, float **vertices,
+ int *numTriangles, int **triangles,
+ int useGlobalCoords, int modifierIndex)
+{
+ DerivedMesh *dm = NULL;
+ MVert *mvert;
+ MFace *mface;
+ int countTris=0, i, totvert, totface;
+ float *verts;
+ int *tris;
+
+ dm = mesh_create_derived_index_render(ob, CD_MASK_BAREMESH, modifierIndex);
+ //dm = mesh_create_derived_no_deform(ob,NULL);
+
+ mvert = dm->getVertArray(dm);
+ mface = dm->getFaceArray(dm);
+ totvert = dm->getNumVerts(dm);
+ totface = dm->getNumFaces(dm);
+
+ *numVertices = totvert;
+ verts = MEM_callocN( totvert*3*sizeof(float), "elbeemmesh_vertices");
+ for(i=0; i<totvert; i++) {
+ VECCOPY( &verts[i*3], mvert[i].co);
+ if(useGlobalCoords) { Mat4MulVecfl(ob->obmat, &verts[i*3]); }
+ }
+ *vertices = verts;
+
+ for(i=0; i<totface; i++) {
+ countTris++;
+ if(mface[i].v4) { countTris++; }
+ }
+ *numTriangles = countTris;
+ tris = MEM_callocN( countTris*3*sizeof(int), "elbeemmesh_triangles");
+ countTris = 0;
+ for(i=0; i<totface; i++) {
+ int face[4];
+ face[0] = mface[i].v1;
+ face[1] = mface[i].v2;
+ face[2] = mface[i].v3;
+ face[3] = mface[i].v4;
+
+ tris[countTris*3+0] = face[0];
+ tris[countTris*3+1] = face[1];
+ tris[countTris*3+2] = face[2];
+ countTris++;
+ if(face[3]) {
+ tris[countTris*3+0] = face[0];
+ tris[countTris*3+1] = face[2];
+ tris[countTris*3+2] = face[3];
+ countTris++;
+ }
+ }
+ *triangles = tris;
+
+ dm->release(dm);
+}
+
+/* read zipped fluidsim velocities into the co's of the fluidsimsettings normals struct */
+void readVelgz(char *filename, Object *srcob)
+{
+ int wri, i, j;
+ float wrf;
+ gzFile gzf;
+ MVert *vverts = srcob->fluidsimSettings->meshSurfNormals;
+ int len = strlen(filename);
+ Mesh *mesh = srcob->data;
+ // mesh and vverts have to be valid from loading...
+
+ // clean up in any case
+ for(i=0; i<mesh->totvert;i++)
+ {
+ for(j=0; j<3; j++)
+ {
+ vverts[i].co[j] = 0.;
+ }
+ }
+ if(srcob->fluidsimSettings->domainNovecgen>0) return;
+
+ if(len<7)
+ {
+ return;
+ }
+
+ // .bobj.gz , correct filename
+ // 87654321
+ filename[len-6] = 'v';
+ filename[len-5] = 'e';
+ filename[len-4] = 'l';
+
+ gzf = gzopen(filename, "rb");
+ if (!gzf)
+ return;
+
+ gzread(gzf, &wri, sizeof( wri ));
+ if(wri != mesh->totvert)
+ {
+ return;
+ }
+
+ for(i=0; i<mesh->totvert;i++)
+ {
+ for(j=0; j<3; j++)
+ {
+ gzread(gzf, &wrf, sizeof( wrf ));
+ vverts[i].co[j] = wrf;
+ }
+ }
+
+ gzclose(gzf);
+}
+
+
+#endif // DISABLE_ELBEEM
+
FLUIDSIM_VISC, FLUIDSIM_TIME,
FLUIDSIM_GRAV_X , FLUIDSIM_GRAV_Y , FLUIDSIM_GRAV_Z ,
FLUIDSIM_VEL_X , FLUIDSIM_VEL_Y , FLUIDSIM_VEL_Z ,
- FLUIDSIM_ACTIVE
+ FLUIDSIM_ACTIVE,
+ FLUIDSIM_ATTR_FORCE_STR, FLUIDSIM_ATTR_FORCE_RADIUS,
+ FLUIDSIM_VEL_FORCE_STR, FLUIDSIM_VEL_FORCE_RADIUS,
};
int part_ar[PART_TOTIPO]= {
#include "MEM_guardedalloc.h"
#include "DNA_armature_types.h"
+#include "DNA_camera_types.h"
#include "DNA_cloth_types.h"
+#include "DNA_curve_types.h"
#include "DNA_effect_types.h"
#include "DNA_material_types.h"
#include "DNA_mesh_types.h"
#include "DNA_particle_types.h"
#include "DNA_scene_types.h"
#include "DNA_texture_types.h"
-#include "DNA_curve_types.h"
-#include "DNA_camera_types.h"
#include "BLI_editVert.h"
#include "BKE_main.h"
#include "BKE_anim.h"
#include "BKE_bad_level_calls.h"
+#include "BKE_bmesh.h"
+#include "BKE_booleanops.h"
#include "BKE_cloth.h"
#include "BKE_collision.h"
+#include "BKE_cdderivedmesh.h"
#include "BKE_curve.h"
#include "BKE_customdata.h"
-#include "BKE_global.h"
-#include "BKE_cdderivedmesh.h"
#include "BKE_DerivedMesh.h"
-#include "BKE_booleanops.h"
#include "BKE_displist.h"
-#include "BKE_modifier.h"
+#include "BKE_fluidsim.h"
+#include "BKE_global.h"
#include "BKE_lattice.h"
#include "BKE_library.h"
-#include "BKE_subsurf.h"
-#include "BKE_object.h"
-#include "BKE_mesh.h"
-#include "BKE_softbody.h"
-#include "BKE_cloth.h"
#include "BKE_material.h"
+#include "BKE_mesh.h"
+#include "BKE_modifier.h"
+#include "BKE_object.h"
#include "BKE_particle.h"
#include "BKE_pointcache.h"
+#include "BKE_softbody.h"
+#include "BKE_subsurf.h"
#include "BKE_texture.h"
#include "BKE_utildefines.h"
+
#include "depsgraph_private.h"
-#include "BKE_bmesh.h"
#include "BKE_deform.h"
#include "BKE_shrinkwrap.h"
}
return derivedData;
}
+
+/* Fluidsim */
+static void fluidsimModifier_initData(ModifierData *md)
+{
+ FluidsimModifierData *fluidmd= (FluidsimModifierData*) md;
+
+ fluidsim_init(fluidmd);
+}
+static void fluidsimModifier_freeData(ModifierData *md)
+{
+ FluidsimModifierData *fluidmd= (FluidsimModifierData*) md;
+
+ fluidsim_free(fluidmd);
+}
+
+static void fluidsimModifier_copyData(ModifierData *md, ModifierData *target)
+{
+ FluidsimModifierData *fluidmd= (FluidsimModifierData*) md;
+ FluidsimModifierData *tfluidmd= (FluidsimModifierData*) target;
+
+ if(tfluidmd->fss)
+ MEM_freeN(tfluidmd->fss);
+
+ tfluidmd->fss = MEM_dupallocN(fluidmd->fss);
+}
+
+static DerivedMesh * fluidsimModifier_applyModifier(
+ ModifierData *md, Object *ob, DerivedMesh *derivedData,
+ int useRenderParams, int isFinalCalc)
+{
+ FluidsimModifierData *fluidmd= (FluidsimModifierData*) md;
+ DerivedMesh *result = NULL;
+
+ /* check for alloc failing */
+ if(!fluidmd->fss)
+ {
+ fluidsimModifier_initData(md);
+
+ if(!fluidmd->fss)
+ return derivedData;
+ }
+
+ result = fluidsimModifier_do(fluidmd, ob, derivedData, useRenderParams, isFinalCalc);
+
+ if(result)
+ {
+ return result;
+ }
+
+ return derivedData;
+}
+
+static void fluidsimModifier_updateDepgraph(
+ ModifierData *md, DagForest *forest,
+ Object *ob, DagNode *obNode)
+{
+ FluidsimModifierData *fluidmd= (FluidsimModifierData*) md;
+ Base *base;
+
+ if(fluidmd && fluidmd->fss)
+ {
+ if(fluidmd->fss->type == OB_FLUIDSIM_DOMAIN)
+ {
+ for(base = G.scene->base.first; base; base= base->next)
+ {
+ Object *ob1= base->object;
+ if(ob1 != ob)
+ {
+ FluidsimModifierData *fluidmdtmp = (FluidsimModifierData *)modifiers_findByType(ob1, eModifierType_Fluidsim);
+
+ // only put dependancies from NON-DOMAIN fluids in here
+ if(fluidmdtmp && fluidmdtmp->fss && (fluidmdtmp->fss->type!=OB_FLUIDSIM_DOMAIN))
+ {
+ DagNode *curNode = dag_get_node(forest, ob1);
+ dag_add_relation(forest, curNode, obNode, DAG_RL_DATA_DATA|DAG_RL_OB_DATA, "Fluidsim Object");
+ }
+ }
+ }
+ }
+ }
+}
+
+static int fluidsimModifier_dependsOnTime(ModifierData *md)
+{
+ return 1;
+}
+
/* MeshDeform */
static void meshdeformModifier_initData(ModifierData *md)
mti->dependsOnTime = explodeModifier_dependsOnTime;
mti->requiredDataMask = explodeModifier_requiredDataMask;
mti->applyModifier = explodeModifier_applyModifier;
+
+ mti = INIT_TYPE(Fluidsim);
+ mti->type = eModifierTypeType_Nonconstructive
+ | eModifierTypeFlag_RequiresOriginalData;
+ mti->flags = eModifierTypeFlag_AcceptsMesh;
+ mti->initData = fluidsimModifier_initData;
+ mti->freeData = fluidsimModifier_freeData;
+ mti->copyData = fluidsimModifier_copyData;
+ mti->dependsOnTime = fluidsimModifier_dependsOnTime;
+ mti->applyModifier = fluidsimModifier_applyModifier;
+ mti->updateDepgraph = fluidsimModifier_updateDepgraph;
mti = INIT_TYPE(Shrinkwrap);
mti->type = eModifierTypeType_OnlyDeform;
MEM_freeN(ob->pd);
}
if(ob->soft) sbFree(ob->soft);
- if(ob->fluidsimSettings) fluidsimSettingsFree(ob->fluidsimSettings);
if(ob->gpulamp.first) GPU_lamp_free(ob);
}
}
obn->soft= copy_softbody(ob->soft);
- /* NT copy fluid sim setting memory */
- if(obn->fluidsimSettings) {
- obn->fluidsimSettings = fluidsimSettingsCopy(ob->fluidsimSettings);
- /* copying might fail... */
- if(obn->fluidsimSettings) {
- obn->fluidsimSettings->orgMesh = (Mesh *)obn->data;
- }
- }
-
copy_object_particlesystems(obn, ob);
obn->derivedDeform = NULL;
}
static void particles_fluid_step(Object *ob, ParticleSystem *psys, int cfra)
-{
+{
if(psys->particles){
MEM_freeN(psys->particles);
psys->particles = 0;
/* fluid sim particle import handling, actual loading of particles from file */
#ifndef DISABLE_ELBEEM
- if( (1) && (ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && // broken, disabled for now!
- (ob->fluidsimSettings)) {
- ParticleSettings *part = psys->part;
- ParticleData *pa=0;
- char *suffix = "fluidsurface_particles_####";
- char *suffix2 = ".gz";
- char filename[256];
- char debugStrBuffer[256];
- int curFrame = G.scene->r.cfra -1; // warning - sync with derived mesh fsmesh loading
- int p, j, numFileParts, totpart;
- int readMask, activeParts = 0, fileParts = 0;
- gzFile gzf;
-
- if(ob==G.obedit) // off...
- return;
-
- // ok, start loading
- strcpy(filename, ob->fluidsimSettings->surfdataPath);
- strcat(filename, suffix);
- BLI_convertstringcode(filename, G.sce);
- BLI_convertstringframe(filename, curFrame); // fixed #frame-no
- strcat(filename, suffix2);
-
- gzf = gzopen(filename, "rb");
- if (!gzf) {
- snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s' \n", filename);
- //elbeemDebugOut(debugStrBuffer);
- return;
- }
-
- gzread(gzf, &totpart, sizeof(totpart));
- numFileParts = totpart;
- totpart = (G.rendering)?totpart:(part->disp*totpart)/100;
+ {
+ FluidsimModifierData *fluidmd = (FluidsimModifierData *)modifiers_findByType(ob, eModifierType_Fluidsim);
- part->totpart= totpart;
- part->sta=part->end = 1.0f;
- part->lifetime = G.scene->r.efra + 1;
-
- /* initialize particles */
- realloc_particles(ob, psys, part->totpart);
- initialize_all_particles(ob, psys, 0);
-
- // set up reading mask
- readMask = ob->fluidsimSettings->typeFlags;
-
- for(p=0, pa=psys->particles; p<totpart; p++, pa++) {
- int ptype=0;
-
- gzread(gzf, &ptype, sizeof( ptype ));
- if(ptype&readMask) {
- activeParts++;
-
- gzread(gzf, &(pa->size), sizeof( float ));
-
- pa->size /= 10.0f;
-
- for(j=0; j<3; j++) {
- float wrf;
- gzread(gzf, &wrf, sizeof( wrf ));
- pa->state.co[j] = wrf;
- //fprintf(stderr,"Rj%d ",j);
- }
- for(j=0; j<3; j++) {
- float wrf;
- gzread(gzf, &wrf, sizeof( wrf ));
- pa->state.vel[j] = wrf;
- }
-
- pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0f;
- pa->state.rot[0] = 1.0;
- pa->state.rot[1] = pa->state.rot[2] = pa->state.rot[3] = 0.0;
-
- pa->alive = PARS_ALIVE;
- //if(a<25) fprintf(stderr,"FSPARTICLE debug set %s , a%d = %f,%f,%f , life=%f \n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
- } else {
- // skip...
- for(j=0; j<2*3+1; j++) {
- float wrf; gzread(gzf, &wrf, sizeof( wrf ));
+ if( fluidmd && fluidmd->fss) {
+ FluidsimSettings *fss= fluidmd->fss;
+ ParticleSettings *part = psys->part;
+ ParticleData *pa=0;
+ char *suffix = "fluidsurface_particles_####";
+ char *suffix2 = ".gz";
+ char filename[256];
+ char debugStrBuffer[256];
+ int curFrame = G.scene->r.cfra -1; // warning - sync with derived mesh fsmesh loading
+ int p, j, numFileParts, totpart;
+ int readMask, activeParts = 0, fileParts = 0;
+ gzFile gzf;
+
+ if(ob==G.obedit) // off...
+ return;
+
+ // ok, start loading
+ strcpy(filename, fss->surfdataPath);
+ strcat(filename, suffix);
+ BLI_convertstringcode(filename, G.sce);
+ BLI_convertstringframe(filename, curFrame); // fixed #frame-no
+ strcat(filename, suffix2);
+
+ gzf = gzopen(filename, "rb");
+ if (!gzf) {
+ snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s' \n", filename);
+ //elbeemDebugOut(debugStrBuffer);
+ return;
+ }
+
+ gzread(gzf, &totpart, sizeof(totpart));
+ numFileParts = totpart;
+ totpart = (G.rendering)?totpart:(part->disp*totpart)/100;
+
+ part->totpart= totpart;
+ part->sta=part->end = 1.0f;
+ part->lifetime = G.scene->r.efra + 1;
+
+ /* initialize particles */
+ realloc_particles(ob, psys, part->totpart);
+ initialize_all_particles(ob, psys, 0);
+
+ // set up reading mask
+ readMask = fss->typeFlags;
+
+ for(p=0, pa=psys->particles; p<totpart; p++, pa++) {
+ int ptype=0;
+
+ gzread(gzf, &ptype, sizeof( ptype ));
+ if(ptype&readMask) {
+ activeParts++;
+
+ gzread(gzf, &(pa->size), sizeof( float ));
+
+ pa->size /= 10.0f;
+
+ for(j=0; j<3; j++) {
+ float wrf;
+ gzread(gzf, &wrf, sizeof( wrf ));
+ pa->state.co[j] = wrf;
+ //fprintf(stderr,"Rj%d ",j);
+ }
+ for(j=0; j<3; j++) {
+ float wrf;
+ gzread(gzf, &wrf, sizeof( wrf ));
+ pa->state.vel[j] = wrf;
+ }
+
+ pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0f;
+ pa->state.rot[0] = 1.0;
+ pa->state.rot[1] = pa->state.rot[2] = pa->state.rot[3] = 0.0;
+
+ pa->alive = PARS_ALIVE;
+ //if(a<25) fprintf(stderr,"FSPARTICLE debug set %s , a%d = %f,%f,%f , life=%f \n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
+ } else {
+ // skip...
+ for(j=0; j<2*3+1; j++) {
+ float wrf; gzread(gzf, &wrf, sizeof( wrf ));
+ }
}
+ fileParts++;
}
- fileParts++;
- }
- gzclose( gzf );
-
- totpart = psys->totpart = activeParts;
- snprintf(debugStrBuffer,256,"readFsPartData::done - particles:%d, active:%d, file:%d, mask:%d \n", psys->totpart,activeParts,fileParts,readMask);
- elbeemDebugOut(debugStrBuffer);
- } // fluid sim particles done
+ gzclose( gzf );
+
+ totpart = psys->totpart = activeParts;
+ snprintf(debugStrBuffer,256,"readFsPartData::done - particles:%d, active:%d, file:%d, mask:%d \n", psys->totpart,activeParts,fileParts,readMask);
+ elbeemDebugOut(debugStrBuffer);
+ } // fluid sim particles done
+ }
#endif // DISABLE_ELBEEM
}
#ifdef WIN32
#include "winsock2.h"
#include "BLI_winstuff.h"
+#ifndef INT_MAX
+#include "limits.h"
+#endif
#endif
#include <stdio.h> // for printf fopen fwrite fclose sprintf FILE
}
act= act->next;
}
-
- if(ob->fluidsimSettings) {
- ob->fluidsimSettings->ipo = newlibadr_us(fd, ob->id.lib, ob->fluidsimSettings->ipo);
+
+ {
+ FluidsimModifierData *fluidmd = (FluidsimModifierData *)modifiers_findByType(ob, eModifierType_Fluidsim);
+
+ if(fluidmd && fluidmd->fss)
+ fluidmd->fss->ipo = newlibadr_us(fd, ob->id.lib, fluidmd->fss->ipo);
}
/* texture field */
}
}
+ else if (md->type==eModifierType_Fluidsim) {
+ FluidsimModifierData *fluidmd = (FluidsimModifierData*) md;
+
+ fluidmd->fss= newdataadr(fd, fluidmd->fss);
+ }
else if (md->type==eModifierType_Collision) {
CollisionModifierData *collmd = (CollisionModifierData*) md;
direct_link_pointcache(fd, sb->pointcache);
}
ob->fluidsimSettings= newdataadr(fd, ob->fluidsimSettings); /* NT */
- if(ob->fluidsimSettings) {
- // reinit mesh pointers
- ob->fluidsimSettings->orgMesh = NULL; //ob->data;
- ob->fluidsimSettings->meshSurface = NULL;
- ob->fluidsimSettings->meshBB = NULL;
- ob->fluidsimSettings->meshSurfNormals = NULL;
- }
link_list(fd, &ob->particlesystem);
direct_link_particlesystems(fd,&ob->particlesystem);
}
}
- if(ob->fluidsimSettings && ob->fluidsimSettings->type == OB_FLUIDSIM_PARTICLE)
- part->type = PART_FLUID;
+
+ {
+ FluidsimModifierData *fluidmd = (FluidsimModifierData *)modifiers_findByType(ob, eModifierType_Fluidsim);
+ if(fluidmd && fluidmd->fss && fluidmd->fss->type == OB_FLUIDSIM_PARTICLE)
+ part->type = PART_FLUID;
+ }
free_effects(&ob->effect);
la->sun_intensity = 1.0;
}
}
+
+ // convert fluids to modifier
+ if(main->versionfile <= 246 && main->subversionfile < 1)
+ {
+ Object *ob;
+
+ for(ob = main->object.first; ob; ob= ob->id.next) {
+ if(ob->fluidsimSettings)
+ {
+ FluidsimModifierData *fluidmd = (FluidsimModifierData *)modifier_new(eModifierType_Fluidsim);
+ BLI_addhead(&ob->modifiers, (ModifierData *)fluidmd);
+
+ MEM_freeN(fluidmd->fss);
+ fluidmd->fss = MEM_dupallocN(ob->fluidsimSettings);
+ fluidmd->fss->ipo = newlibadr_us(fd, ob->id.lib, ob->fluidsimSettings->ipo);
+ MEM_freeN(ob->fluidsimSettings);
+
+ fluidmd->fss->lastgoodframe = INT_MAX;
+ fluidmd->fss->flag = 0;
+ }
+ }
+ }
+
if(main->versionfile < 246 || (main->versionfile == 246 && main->subversionfile < 1)) {
Mesh *me;
writestruct(wd, DATA, "ClothCollSettings", 1, clmd->coll_parms);
writestruct(wd, DATA, "PointCache", 1, clmd->point_cache);
}
+ else if(md->type==eModifierType_Fluidsim) {
+ FluidsimModifierData *fluidmd = (FluidsimModifierData*) md;
+
+ writestruct(wd, DATA, "FluidsimSettings", 1, fluidmd->fss);
+ }
else if (md->type==eModifierType_Collision) {
/*
writestruct(wd, DATA, "PartDeflect", 1, ob->pd);
writestruct(wd, DATA, "SoftBody", 1, ob->soft);
if(ob->soft) writestruct(wd, DATA, "PointCache", 1, ob->soft->pointcache);
- writestruct(wd, DATA, "FluidsimSettings", 1, ob->fluidsimSettings); // NT
write_particlesystems(wd, &ob->particlesystem);
write_modifiers(wd, &ob->modifiers);
/* ****** FluidSim (ID_FLUIDSIM) ****** */
-#define FLUIDSIM_TOTIPO 9
-#define FLUIDSIM_TOTNAM 9
+#define FLUIDSIM_TOTIPO 13
+#define FLUIDSIM_TOTNAM 13
#define FLUIDSIM_VISC 1
#define FLUIDSIM_TIME 2
#define FLUIDSIM_ACTIVE 9
-/* ******* Particle (ID_PA) ******** */
+#define FLUIDSIM_ATTR_FORCE_STR 10
+#define FLUIDSIM_ATTR_FORCE_RADIUS 11
+#define FLUIDSIM_VEL_FORCE_STR 12
+#define FLUIDSIM_VEL_FORCE_RADIUS 13
+
+/* ******************** */
+/* particle ipos */
+/* ******* Particle (ID_PA) ******** */
#define PART_TOTIPO 25
#define PART_TOTNAM 25
eModifierType_Collision,
eModifierType_Bevel,
eModifierType_Shrinkwrap,
+ eModifierType_Fluidsim,
NUM_MODIFIER_TYPES
} ModifierType;
float protect;
} ExplodeModifierData;
+typedef struct FluidsimModifierData {
+ ModifierData modifier;
+
+ struct FluidsimSettings *fss; /* definition is is DNA_object_fluidsim.h */
+ struct PointCache *point_cache; /* definition is in DNA_object_force.h */
+} FluidsimModifierData;
+
typedef struct ShrinkwrapModifierData {
ModifierData modifier;
float surfaceSmoothing;
/* number of surface subdivisions*/
int surfaceSubdivs;
- int unusedDNADummy;
+ int flag; /* GUI flags */
/* particle display - size scaling, and alpha influence */
float particleInfSize, particleInfAlpha;
/* save fluidsurface normals in mvert.no, and surface vertex velocities (if available) in mvert.co */
struct MVert *meshSurfNormals;
+
+ /* Fluid control settings */
+ float cpsTimeStart;
+ float cpsTimeEnd;
+ float cpsQuality;
+
+ float attractforceStrength;
+ float attractforceRadius;
+ float velocityforceStrength;
+ float velocityforceRadius;
+
+ int lastgoodframe;
} FluidsimSettings;
#define OB_FLUIDSIM_INFLOW 16
#define OB_FLUIDSIM_OUTFLOW 32
#define OB_FLUIDSIM_PARTICLE 64
+#define OB_FLUIDSIM_CONTROL 128
#define OB_TYPEFLAG_START 0
#define OB_FSGEO_THIN (1<<(OB_TYPEFLAG_START+1))
#define OB_FSPART_NEWPART (1<<3)
#define OB_FSPART_FLOAT (1<<4)
+// new fluid bit flags for fss->flags - dg
+#define OB_FLUIDSIM_REVERSE (1 << 0)
+
#ifdef __cplusplus
}
#endif
MEM_freeN(vtangents);
}
-// NT same as calc_vertexnormals, but dont modify the existing vertex normals
-// only recalculate other render data. If this is at some point used for other things than fluidsim,
-// this could be made on option for the normal calc_vertexnormals
-static void calc_fluidsimnormals(Render *re, ObjectRen *obr, int do_nmap_tangent)
-{
- int a;
-
- /* dont clear vertex normals here */
- // OFF for(a=0; a<obr->totvert; a++) { VertRen *ver= RE_findOrAddVert(obr, a); ver->n[0]=ver->n[1]=ver->n[2]= 0.0; }
- /* calculate cos of angles and point-masses, use as weight factor to add face normal to vertex */
- for(a=0; a<obr->totvlak; a++) {
- VlakRen *vlr= RE_findOrAddVlak(obr, a);
- if(vlr->flag & ME_SMOOTH) {
- VertRen *v1= vlr->v1;
- VertRen *v2= vlr->v2;
- VertRen *v3= vlr->v3;
- VertRen *v4= vlr->v4;
- float n1[3], n2[3], n3[3], n4[3];
- float fac1, fac2, fac3, fac4=0.0f;
-
- if(re->flag & R_GLOB_NOPUNOFLIP)
- vlr->flag |= R_NOPUNOFLIP;
-
- VecSubf(n1, v2->co, v1->co);
- Normalize(n1);
- VecSubf(n2, v3->co, v2->co);
- Normalize(n2);
- if(v4==NULL) {
- VecSubf(n3, v1->co, v3->co);
- Normalize(n3);
- fac1= saacos(-n1[0]*n3[0]-n1[1]*n3[1]-n1[2]*n3[2]);
- fac2= saacos(-n1[0]*n2[0]-n1[1]*n2[1]-n1[2]*n2[2]);
- fac3= saacos(-n2[0]*n3[0]-n2[1]*n3[1]-n2[2]*n3[2]);
- }
- else {
- VecSubf(n3, v4->co, v3->co);
- Normalize(n3);
- VecSubf(n4, v1->co, v4->co);
- Normalize(n4);
-
- fac1= saacos(-n4[0]*n1[0]-n4[1]*n1[1]-n4[2]*n1[2]);
- &
projects / blender-staging.git / commitdiff