diff options
Diffstat (limited to 'intern')
35 files changed, 3616 insertions, 738 deletions
diff --git a/intern/SoundSystem/SConscript b/intern/SoundSystem/SConscript index baf680f03f0..be19c4b7915 100644 --- a/intern/SoundSystem/SConscript +++ b/intern/SoundSystem/SConscript @@ -14,4 +14,7 @@ if env['WITH_BF_OPENAL']: else: defs = 'NO_SOUND' +if not env['WITH_BF_SDL']: + defs += ' DISABLE_SDL' + env.BlenderLib ('bf_soundsystem', sources, Split(incs), Split(defs), libtype=['core','player'], priority = [20,140] ) diff --git a/intern/SoundSystem/intern/SND_Utils.cpp b/intern/SoundSystem/intern/SND_Utils.cpp index 7d441964b24..05101a7e6e6 100644 --- a/intern/SoundSystem/intern/SND_Utils.cpp +++ b/intern/SoundSystem/intern/SND_Utils.cpp @@ -34,11 +34,6 @@ #include "SND_Utils.h" #include "SoundDefines.h" #include "SND_DependKludge.h" -/* -extern "C" { -#include "license_key.h" -} -*/ #include <stdio.h> #include <stdlib.h> #include <fcntl.h> diff --git a/intern/SoundSystem/sdl/SND_SDLCDDevice.cpp b/intern/SoundSystem/sdl/SND_SDLCDDevice.cpp index 8bb6c642b11..0ab0fa94c7b 100644 --- a/intern/SoundSystem/sdl/SND_SDLCDDevice.cpp +++ b/intern/SoundSystem/sdl/SND_SDLCDDevice.cpp @@ -55,6 +55,10 @@ SND_SDLCDDevice::SND_SDLCDDevice() : void SND_SDLCDDevice::init() { +#ifdef DISABLE_SDL + fprintf(stderr, "Blender compiled without SDL, no CDROM support\n"); + return; +#else if (SDL_InitSubSystem(SDL_INIT_CDROM)) { fprintf(stderr, "Error initializing CDROM\n"); @@ -75,19 +79,23 @@ void SND_SDLCDDevice::init() /* Did if open? Check if cdrom is NULL */ if(!m_cdrom) { - fprintf(stderr, "Couldn't open drive: %s", SDL_GetError()); + fprintf(stderr, "Couldn't open drive: %s\n", SDL_GetError()); return; } +#endif } SND_SDLCDDevice::~SND_SDLCDDevice() { +#ifndef DISABLE_SDL StopCD(); SDL_CDClose(m_cdrom); +#endif } void SND_SDLCDDevice::NextFrame() { +#ifndef DISABLE_SDL m_frame++; m_frame &= 127; @@ -111,20 +119,24 @@ void SND_SDLCDDevice::NextFrame() } } +#endif } void SND_SDLCDDevice::PlayCD(int track) { +#ifndef DISABLE_SDL if ( m_cdrom && CD_INDRIVE(SDL_CDStatus(m_cdrom)) ) { SDL_CDPlayTracks(m_cdrom, track-1, 0, track, 0); m_cdplaying = true; m_cdtrack = track; } +#endif } void SND_SDLCDDevice::PauseCD(bool pause) { +#ifndef DISABLE_SDL if (!m_cdrom) return; @@ -132,13 +144,16 @@ void SND_SDLCDDevice::PauseCD(bool pause) SDL_CDPause(m_cdrom); else SDL_CDResume(m_cdrom); +#endif } void SND_SDLCDDevice::StopCD() { +#ifndef DISABLE_SDL if (m_cdrom) SDL_CDStop(m_cdrom); m_cdplaying = false; +#endif } void SND_SDLCDDevice::SetCDPlaymode(int playmode) diff --git a/intern/elbeem/CMakeLists.txt b/intern/elbeem/CMakeLists.txt index bd88bc0c5f7..03fd4a3fefc 100644 --- a/intern/elbeem/CMakeLists.txt +++ b/intern/elbeem/CMakeLists.txt @@ -24,7 +24,7 @@ # # ***** 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) diff --git a/intern/elbeem/SConscript b/intern/elbeem/SConscript index bdcb0507987..68892ac261c 100644 --- a/intern/elbeem/SConscript +++ b/intern/elbeem/SConscript @@ -5,13 +5,14 @@ Import('env') sources = env.Glob('intern/*.cpp') -defs = ' NOGUI ELBEEM_BLENDER=1' +defs = 'NOGUI ELBEEM_BLENDER=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 ) diff --git a/intern/elbeem/extern/LBM_fluidsim.h b/intern/elbeem/extern/LBM_fluidsim.h index 6e7fc86213c..dd8c63d0d39 100644 --- a/intern/elbeem/extern/LBM_fluidsim.h +++ b/intern/elbeem/extern/LBM_fluidsim.h @@ -63,11 +63,6 @@ void loadFluidsimMesh(struct Object *srcob, int useRenderParams); // WARNING - implemented in intern/elbeem/blendercall.cpp int performElbeemSimulation(char *cfgfilename); -/* init axis aligned BB for mesh object */ -// implemented in source/blender/blenkernel/intern/DerivedMesh.c -void fluidsimGetAxisAlignedBB(struct Mesh *mesh, float obmat[][4], - /*RET*/ float start[3], /*RET*/ float size[3], /*RET*/ struct Mesh **bbmesh ); - #endif diff --git a/intern/elbeem/extern/elbeem.h b/intern/elbeem/extern/elbeem.h index 2a594dd07e6..e29890aba44 100644 --- a/intern/elbeem/extern/elbeem.h +++ b/intern/elbeem/extern/elbeem.h @@ -32,7 +32,7 @@ typedef struct elbeemSimulationSettings { 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]; @@ -97,10 +97,13 @@ typedef struct elbeemSimulationSettings { // 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 @@ -113,7 +116,7 @@ typedef struct elbeemSimulationSettings { // 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; @@ -155,6 +158,20 @@ typedef struct elbeemMesh { /* 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 @@ -170,6 +187,9 @@ void elbeemResetSettings(struct elbeemSimulationSettings*); // 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*); @@ -223,6 +243,7 @@ double elbeemEstimateMemreq(int res, // 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)) @@ -232,6 +253,7 @@ double elbeemEstimateMemreq(int res, #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 ) diff --git a/intern/elbeem/intern/controlparticles.cpp b/intern/elbeem/intern/controlparticles.cpp new file mode 100644 index 00000000000..4b496e4a3a2 --- /dev/null +++ b/intern/elbeem/intern/controlparticles.cpp @@ -0,0 +1,1461 @@ +// -------------------------------------------------------------------------- +// +// 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; + + // 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]); + + 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 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; + size_t len = 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; + size_t len = 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; +} + diff --git a/intern/elbeem/intern/controlparticles.h b/intern/elbeem/intern/controlparticles.h new file mode 100644 index 00000000000..712dfc40273 --- /dev/null +++ b/intern/elbeem/intern/controlparticles.h @@ -0,0 +1,301 @@ +// -------------------------------------------------------------------------- +// +// 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 + diff --git a/intern/elbeem/intern/elbeem.cpp b/intern/elbeem/intern/elbeem.cpp index b2779f51c3b..f7923322d5b 100644 --- a/intern/elbeem/intern/elbeem.cpp +++ b/intern/elbeem/intern/elbeem.cpp @@ -30,6 +30,7 @@ int guiRoiMaxLev=6, guiRoiMinLev=0; ntlWorld *gpWorld = NULL; + // API // reset elbeemSimulationSettings struct with defaults @@ -95,6 +96,13 @@ int elbeemInit() { return 0; } +// fluidsim end +extern "C" +int elbeemFree() { + + return 0; +} + // start fluidsim init extern "C" int elbeemAddDomain(elbeemSimulationSettings *settings) { @@ -158,13 +166,27 @@ void elbeemResetMesh(elbeemMesh *mesh) { /* 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) { @@ -176,9 +198,9 @@ int elbeemAddMesh(elbeemMesh *mesh) { 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 ); @@ -195,17 +217,34 @@ int elbeemAddMesh(elbeemMesh *mesh) { 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 ); @@ -227,6 +266,7 @@ int elbeemSimulate(void) { 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 { diff --git a/intern/elbeem/intern/elbeem.h b/intern/elbeem/intern/elbeem.h deleted file mode 100644 index 2a594dd07e6..00000000000 --- a/intern/elbeem/intern/elbeem.h +++ /dev/null @@ -1,240 +0,0 @@ -/****************************************************************************** - * - * 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 diff --git a/intern/elbeem/intern/elbeem_control.cpp b/intern/elbeem/intern/elbeem_control.cpp new file mode 100644 index 00000000000..800167baf73 --- /dev/null +++ b/intern/elbeem/intern/elbeem_control.cpp @@ -0,0 +1,25 @@ +/****************************************************************************** + * + * 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; +} + diff --git a/intern/elbeem/intern/elbeem_control.h b/intern/elbeem/intern/elbeem_control.h new file mode 100644 index 00000000000..70a58feda89 --- /dev/null +++ b/intern/elbeem/intern/elbeem_control.h @@ -0,0 +1,62 @@ +/****************************************************************************** + * + * 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 diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp index 1b0ba13c707..751a48cebba 100644 --- a/intern/elbeem/intern/isosurface.cpp +++ b/intern/elbeem/intern/isosurface.cpp @@ -13,8 +13,8 @@ #include <algorithm> #include <stdio.h> -#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__))) -#include <ieeefp.h> +#ifdef sun +#include "ieeefp.h" #endif // just use default rounding for platforms where its not available diff --git a/intern/elbeem/intern/mvmcoords.cpp b/intern/elbeem/intern/mvmcoords.cpp new file mode 100644 index 00000000000..ac954b956d0 --- /dev/null +++ b/intern/elbeem/intern/mvmcoords.cpp @@ -0,0 +1,193 @@ +/****************************************************************************** + * +// 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; + } +} + diff --git a/intern/elbeem/intern/mvmcoords.h b/intern/elbeem/intern/mvmcoords.h new file mode 100644 index 00000000000..889f5058a09 --- /dev/null +++ b/intern/elbeem/intern/mvmcoords.h @@ -0,0 +1,89 @@ +/****************************************************************************** + * +// 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 + +#ifdef sun +#include "ieeefp.h" +#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 + diff --git a/intern/elbeem/intern/ntl_geometryclass.h b/intern/elbeem/intern/ntl_geometryclass.h index 9282371d183..545f8c1d54f 100644 --- a/intern/elbeem/intern/ntl_geometryclass.h +++ b/intern/elbeem/intern/ntl_geometryclass.h @@ -37,6 +37,7 @@ class ntlGeometryClass //! Default destructor virtual ~ntlGeometryClass() { delete mpAttrs; + delete mpSwsAttrs; }; //! Return type id diff --git a/intern/elbeem/intern/ntl_geometryobject.cpp b/intern/elbeem/intern/ntl_geometryobject.cpp index f2ebd572682..c7a222af3e5 100644 --- a/intern/elbeem/intern/ntl_geometryobject.cpp +++ b/intern/elbeem/intern/ntl_geometryobject.cpp @@ -41,7 +41,9 @@ ntlGeometryObject::ntlGeometryObject() : 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.) { }; @@ -82,20 +84,21 @@ bool ntlGeometryObject::checkIsAnimated() { /*****************************************************************************/ /* 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) { @@ -330,7 +333,11 @@ void ntlGeometryObject::sceneAddTriangleNoVert(int *trips, 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); @@ -343,8 +350,15 @@ void ntlGeometryObject::initChannels( 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()<< @@ -565,7 +579,7 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) { } } - if( (this-getMeshAnimated()) + if( (this->getMeshAnimated()) || (mcTrans.accessValues().size()>1) // VALIDATE || (mcRot.accessValues().size()>1) || (mcScale.accessValues().size()>1) diff --git a/intern/elbeem/intern/ntl_geometryobject.h b/intern/elbeem/intern/ntl_geometryobject.h index a5afd6b2207..666798385f6 100644 --- a/intern/elbeem/intern/ntl_geometryobject.h +++ b/intern/elbeem/intern/ntl_geometryobject.h @@ -97,11 +97,37 @@ class ntlGeometryObject : public ntlGeometryClass /*! 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? */ @@ -202,6 +228,12 @@ class ntlGeometryObject : public ntlGeometryClass /*! 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: diff --git a/intern/elbeem/intern/ntl_ray.cpp b/intern/elbeem/intern/ntl_ray.cpp index d6593b6b335..242b82085be 100644 --- a/intern/elbeem/intern/ntl_ray.cpp +++ b/intern/elbeem/intern/ntl_ray.cpp @@ -688,7 +688,8 @@ ntlScene::~ntlScene() 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() ); diff --git a/intern/elbeem/intern/ntl_ray.h b/intern/elbeem/intern/ntl_ray.h index 096d5fd61d3..9b77fdcae28 100644 --- a/intern/elbeem/intern/ntl_ray.h +++ b/intern/elbeem/intern/ntl_ray.h @@ -20,6 +20,7 @@ class ntlRay; class ntlTree; class ntlScene; class ntlRenderGlobals; +class ntlGeometryObject; //! store data for an intersection of a ray and a triangle // NOT YET USED @@ -323,8 +324,8 @@ public: /* 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) { diff --git a/intern/elbeem/intern/ntl_world.cpp b/intern/elbeem/intern/ntl_world.cpp index 21d9d63aac0..c7aa3495c93 100644 --- a/intern/elbeem/intern/ntl_world.cpp +++ b/intern/elbeem/intern/ntl_world.cpp @@ -236,10 +236,15 @@ ntlWorld::~ntlWorld() { 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 @@ -895,6 +900,8 @@ ntlRenderGlobals::ntlRenderGlobals() : ntlRenderGlobals::~ntlRenderGlobals() { if(mpOpenGlAttr) delete mpOpenGlAttr; if(mpBlenderAttr) delete mpBlenderAttr; + + } diff --git a/intern/elbeem/intern/simulation_object.cpp b/intern/elbeem/intern/simulation_object.cpp index 9b47ae696af..81e056771fc 100644 --- a/intern/elbeem/intern/simulation_object.cpp +++ b/intern/elbeem/intern/simulation_object.cpp @@ -58,7 +58,8 @@ SimulationObject::~SimulationObject() 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); } diff --git a/intern/elbeem/intern/solver_adap.cpp b/intern/elbeem/intern/solver_adap.cpp index 5616d805232..ef516a578bd 100644 --- a/intern/elbeem/intern/solver_adap.cpp +++ b/intern/elbeem/intern/solver_adap.cpp @@ -11,9 +11,7 @@ #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 diff --git a/intern/elbeem/intern/solver_class.h b/intern/elbeem/intern/solver_class.h index d46f065adfd..f61ef72e8ab 100644 --- a/intern/elbeem/intern/solver_class.h +++ b/intern/elbeem/intern/solver_class.h @@ -100,11 +100,14 @@ // sirdude fix for solaris #if !defined(linux) && defined(sun) +#include "ieeefp.h" #ifndef expf #define expf(x) exp((double)(x)) #endif #endif +#include "solver_control.h" + #if LBM_INCLUDE_TESTSOLVERS==1 #include "solver_test.h" #endif // LBM_INCLUDE_TESTSOLVERS==1 @@ -497,6 +500,12 @@ class LbmFsgrSolver : LbmFloat& debRAC(LbmFloat* s,int l); # endif // FSGR_STRICT_DEBUG==1 + LbmControlData *mpControl; + + void initCpdata(); + void handleCpdata(); + void cpDebugDisplay(int dispset); + bool mUseTestdata; # if LBM_INCLUDE_TESTSOLVERS==1 // test functions @@ -506,10 +515,6 @@ class LbmFsgrSolver : void handleTestdata(); void set3dHeight(int ,int ); - void initCpdata(); - void handleCpdata(); - void cpDebugDisplay(int dispset); - int mMpNum,mMpIndex; int mOrgSizeX; LbmFloat mOrgStartX; diff --git a/intern/elbeem/intern/solver_control.cpp b/intern/elbeem/intern/solver_control.cpp new file mode 100644 index 00000000000..15b9468ed0d --- /dev/null +++ b/intern/elbeem/intern/solver_control.cpp @@ -0,0 +1,1016 @@ +/****************************************************************************** + * + * 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 + + diff --git a/intern/elbeem/intern/solver_control.h b/intern/elbeem/intern/solver_control.h new file mode 100644 index 00000000000..57112b365ce --- /dev/null +++ b/intern/elbeem/intern/solver_control.h @@ -0,0 +1,182 @@ +/****************************************************************************** + * + * 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 + +// 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; \ + } + + +// 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 diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp index a873f3c6987..270531d80be 100644 --- a/intern/elbeem/intern/solver_init.cpp +++ b/intern/elbeem/intern/solver_init.cpp @@ -328,7 +328,8 @@ LbmFsgrSolver::LbmFsgrSolver() : mInit2dYZ(false), mForceTadapRefine(-1), mCutoff(-1) { - // not much to do here... + mpControl = new LbmControlData(); + #if LBM_INCLUDE_TESTSOLVERS==1 mpTest = new LbmTestdata(); mMpNum = mMpIndex = 0; @@ -438,6 +439,8 @@ LbmFsgrSolver::~LbmFsgrSolver() delete mpIso; if(mpPreviewSurface) delete mpPreviewSurface; // cleanup done during scene deletion... + + if(mpControl) delete mpControl; // always output performance estimate debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5); @@ -488,6 +491,8 @@ void LbmFsgrSolver::parseAttrList() mSimulationTime += starttimeskip; if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1); + mpControl->parseControldataAttrList(mpSifAttrs); + #if LBM_INCLUDE_TESTSOLVERS==1 mUseTestdata = 0; mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false); @@ -689,21 +694,21 @@ bool LbmFsgrSolver::initializeSolverMemory() // 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; } @@ -1264,8 +1269,10 @@ bool LbmFsgrSolver::initializeSolverPostinit() { debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10); mInitDone = 1; -#if LBM_INCLUDE_TESTSOLVERS==1 + // init fluid control initCpdata(); + +#if LBM_INCLUDE_TESTSOLVERS==1 initTestdata(); #endif // ELBEEM_PLUGIN!=1 // not inited? dont use... diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp index 13ebf91b696..9df53a824d4 100644 --- a/intern/elbeem/intern/solver_main.cpp +++ b/intern/elbeem/intern/solver_main.cpp @@ -13,11 +13,6 @@ #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 */ /*****************************************************************************/ @@ -58,9 +53,9 @@ void LbmFsgrSolver::stepMain() { // init moving bc's, can change mMaxVlen initMovingObstacles(false); -#if LBM_INCLUDE_TESTSOLVERS==1 + + // handle fluid control handleCpdata(); -#endif // important - keep for tadap LbmFloat lastMass = mCurrentMass; diff --git a/intern/elbeem/intern/solver_relax.h b/intern/elbeem/intern/solver_relax.h index 287b73c77b2..c28e7aacd29 100644 --- a/intern/elbeem/intern/solver_relax.h +++ b/intern/elbeem/intern/solver_relax.h @@ -14,53 +14,6 @@ #else // FSGR_STRICT_DEBUG==1 #define CAUSE_PANIC { this->mPanic=1; } /*set flag*/ #endif // FSGR_STRICT_DEBUG==1 - -#if LBM_INCLUDE_TESTSOLVERS!=1 - -#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \ - ux += (grav)[0]; \ - uy += (grav)[1]; \ - uz += (grav)[2]; - -#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 - /****************************************************************************** * normal relaxation diff --git a/intern/elbeem/make/msvc_7_0/elbeem.vcproj b/intern/elbeem/make/msvc_7_0/elbeem.vcproj index 9de4482e2fd..958d74f8fa0 100644 --- a/intern/elbeem/make/msvc_7_0/elbeem.vcproj +++ b/intern/elbeem/make/msvc_7_0/elbeem.vcproj @@ -19,8 +19,8 @@ <Tool Name="VCCLCompilerTool" Optimization="0" - AdditionalIncludeDirectories="..\..\intern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" - PreprocessorDefinitions="WIN32;_DEBUG;_LIB;NOGUI;ELBEEM_BLENDER=1" + AdditionalIncludeDirectories="..\..\intern;..\..\extern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" + PreprocessorDefinitions="WIN32;_DEBUG;_LIB;NOGUI;ELBEEM_BLENDER=1;LBM_INCLUDE_CONTROL=1" MinimalRebuild="FALSE" BasicRuntimeChecks="3" RuntimeLibrary="1" @@ -72,8 +72,8 @@ ECHO Done CharacterSet="2"> <Tool Name="VCCLCompilerTool" - AdditionalIncludeDirectories="..\..\intern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" - PreprocessorDefinitions="WIN32;NDEBUG;_LIB;NOGUI;ELBEEM_BLENDER=1" + AdditionalIncludeDirectories="..\..\intern;..\..\extern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" + PreprocessorDefinitions="WIN32;NDEBUG;_LIB;NOGUI;ELBEEM_BLENDER=1;LBM_INCLUDE_CONTROL=1" StringPooling="TRUE" MinimalRebuild="FALSE" RuntimeLibrary="0" @@ -126,7 +126,7 @@ ECHO Done <Tool Name="VCCLCompilerTool" Optimization="0" - AdditionalIncludeDirectories="..\..\intern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" + AdditionalIncludeDirectories="..\..\intern;..\..\extern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" PreprocessorDefinitions="WIN32;_DEBUG;_LIB;NOGUI;ELBEEM_BLENDER=1" MinimalRebuild="FALSE" BasicRuntimeChecks="3" @@ -179,7 +179,7 @@ ECHO Done CharacterSet="2"> <Tool Name="VCCLCompilerTool" - AdditionalIncludeDirectories="..\..\intern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" + AdditionalIncludeDirectories="..\..\intern;..\..\extern;..\..\..\..\..\lib\windows\png\include;..\..\..\..\..\lib\windows\zlib\include;..\..\..\..\source\kernel\gen_messaging;..\..\..\..\..\lib\windows\sdl\include" PreprocessorDefinitions="WIN32;NDEBUG;_LIB;NOGUI;ELBEEM_BLENDER=1" StringPooling="TRUE" MinimalRebuild="FALSE" @@ -236,12 +236,21 @@ ECHO Done RelativePath="..\..\intern\attributes.cpp"> </File> <File + RelativePath="..\..\intern\controlparticles.cpp"> + </File> + <File RelativePath="..\..\intern\elbeem.cpp"> </File> <File + RelativePath="..\..\intern\elbeem_control.cpp"> + </File> + <File RelativePath="..\..\intern\isosurface.cpp"> </File> <File + RelativePath="..\..\intern\mvmcoords.cpp"> + </File> + <File RelativePath="..\..\intern\ntl_blenderdumper.cpp"> </File> <File @@ -275,6 +284,9 @@ ECHO Done RelativePath="..\..\intern\solver_adap.cpp"> </File> <File + RelativePath="..\..\intern\solver_control.cpp"> + </File> + <File RelativePath="..\..\intern\solver_init.cpp"> </File> <File @@ -298,7 +310,13 @@ ECHO Done RelativePath="..\..\intern\attributes.h"> </File> <File - RelativePath="..\..\intern\elbeem.h"> + RelativePath="..\..\intern\controlparticles.h"> + </File> + <File + RelativePath="..\..\extern\elbeem.h"> + </File> + <File + RelativePath="..\..\intern\elbeem_control.h"> </File> <File RelativePath="..\..\intern\isosurface.h"> @@ -307,6 +325,9 @@ ECHO Done RelativePath="..\..\intern\mcubes_tables.h"> </File> <File + RelativePath="..\..\intern\mvmcoords.h"> + </File> + <File RelativePath="..\..\intern\ntl_blenderdumper.h"> </File> <File @@ -352,6 +373,9 @@ ECHO Done RelativePath="..\..\intern\solver_class.h"> </File> <File + RelativePath="..\..\intern\solver_control.h"> + </File> + <File RelativePath="..\..\intern\solver_interface.h"> </File> <File diff --git a/intern/ghost/intern/GHOST_SystemX11.cpp b/intern/ghost/intern/GHOST_SystemX11.cpp index 047bc654559..8073756e453 100644 --- a/intern/ghost/intern/GHOST_SystemX11.cpp +++ b/intern/ghost/intern/GHOST_SystemX11.cpp @@ -96,25 +96,17 @@ GHOST_SystemX11( if (!m_display) return; #ifdef __sgi - m_delete_window_atom = XSGIFastInternAtom(m_display, + m_delete_window_atom + = XSGIFastInternAtom(m_display, "WM_DELETE_WINDOW", SGI_XA_WM_DELETE_WINDOW, False); #else - m_delete_window_atom = XInternAtom(m_display, "WM_DELETE_WINDOW", False); + m_delete_window_atom + = XInternAtom(m_display, "WM_DELETE_WINDOW", True); #endif m_wm_protocols= XInternAtom(m_display, "WM_PROTOCOLS", False); m_wm_take_focus= XInternAtom(m_display, "WM_TAKE_FOCUS", False); - m_wm_state= XInternAtom(m_display, "WM_STATE", False); - m_wm_change_state= XInternAtom(m_display, "WM_CHANGE_STATE", False); - m_net_state= XInternAtom(m_display, "_NET_WM_STATE", False); - m_net_max_horz= XInternAtom(m_display, - "_NET_WM_STATE_MAXIMIZED_HORZ", False); - m_net_max_vert= XInternAtom(m_display, - "_NET_WM_STATE_MAXIMIZED_VERT", False); - m_net_fullscreen= XInternAtom(m_display, - "_NET_WM_STATE_FULLSCREEN", False); - m_motif= XInternAtom(m_display, "_MOTIF_WM_HINTS", False); // compute the initial time timeval tv; @@ -522,12 +514,10 @@ GHOST_SystemX11::processEvent(XEvent *xe) GHOST_kEventNDOFButton, window, data); } - } - else if (((Atom)xcme.data.l[0]) == m_wm_take_focus) { + } else if (((Atom)xcme.data.l[0]) == m_wm_take_focus) { /* as ICCCM say, we need reply this event * with a SetInputFocus, the data[1] have - * the valid timestamp (send by the window - * manager). + * the valid timestamp (send by the wm). */ XSetInputFocus(m_display, xcme.window, RevertToParent, xcme.data.l[1]); } else { @@ -547,24 +537,6 @@ GHOST_SystemX11::processEvent(XEvent *xe) // XCrossingEvents pointer leave enter window. break; case MapNotify: - /* - * From ICCCM: - * [ Clients can select for StructureNotify on their - * top-level windows to track transition between - * Normal and Iconic states. Receipt of a MapNotify - * event will indicate a transition to the Normal - * state, and receipt of an UnmapNotify event will - * indicate a transition to the Iconic state. ] - */ - if (window->m_post_init == True) { - /* - * Now we are sure that the window is - * mapped, so only need change the state. - */ - window->setState(window->m_post_state); - window->m_post_init= False; - } - break; case UnmapNotify: break; case MappingNotify: diff --git a/intern/ghost/intern/GHOST_SystemX11.h b/intern/ghost/intern/GHOST_SystemX11.h index 683337b47e1..576577917ba 100644 --- a/intern/ghost/intern/GHOST_SystemX11.h +++ b/intern/ghost/intern/GHOST_SystemX11.h @@ -215,21 +215,9 @@ public: */ virtual void putClipboard(GHOST_TInt8 *buffer, int flag) const; - /** - * Atom used for ICCCM, WM-spec and Motif. - * We only need get this atom at the start, it's relative - * to the display not the window and are public for every - * window that need it. - */ - Atom m_wm_protocols; + /* Atom used for ICCCM. */ Atom m_wm_take_focus; - Atom m_wm_state; - Atom m_wm_change_state; - Atom m_net_state; - Atom m_net_max_horz; - Atom m_net_max_vert; - Atom m_net_fullscreen; - Atom m_motif; + Atom m_wm_protocols; Atom m_delete_window_atom; private : diff --git a/intern/ghost/intern/GHOST_WindowX11.cpp b/intern/ghost/intern/GHOST_WindowX11.cpp index 2f7aee95b82..73d61a30977 100644 --- a/intern/ghost/intern/GHOST_WindowX11.cpp +++ b/intern/ghost/intern/GHOST_WindowX11.cpp @@ -55,16 +55,6 @@ typedef struct { #define MWM_HINTS_DECORATIONS (1L << 1) /* - * A client can't change the window property, that is the - * work of the window manager. We send a ClientMessage - * event to the Root window with the property - * and the Action (WM-spec define this): - */ -#define _NET_WM_STATE_REMOVE 0 -#define _NET_WM_STATE_ADD 1 -#define _NET_WM_STATE_TOGGLE 2 - -/* import bpy I = bpy.data.images['blender.png'] # the 48x48 icon @@ -164,10 +154,10 @@ GHOST_WindowX11( // Set up the minimum atrributes that we require and see if // X can find us a visual matching those requirements. - Atom atoms[2]; - int natom; + int attributes[40], i = 0; - + Atom atoms[2]; + int natom; if(m_stereoVisual) attributes[i++] = GLX_STEREO; @@ -272,26 +262,46 @@ GHOST_WindowX11( // Are we in fullscreen mode - then include // some obscure blut code to remove decorations. - /* - * One of the problem with WM_spec is that can't set a property - * to a window that isn't mapped. That is why we can't "just - * call setState" here. - * - * To fix this, we first need know that the window is really - * mapped waiting for the MapNotify event. - * - * So, m_post_init indicate that we need wait for the MapNotify - * event and then set the window state to the m_post_state. - */ - if ((state != GHOST_kWindowStateNormal) && (state != GHOST_kWindowStateMinimized)) { - m_post_init = True; - m_post_state = state; - } - else { - m_post_init = False; - m_post_state = GHOST_kWindowStateNormal; - } + if (state == GHOST_kWindowStateFullScreen) { + MotifWmHints hints; + Atom atom; + + atom = XInternAtom(m_display, "_MOTIF_WM_HINTS", False); + + if (atom == None) { + GHOST_PRINT("Could not intern X atom for _MOTIF_WM_HINTS.\n"); + } else { + hints.flags = MWM_HINTS_DECORATIONS; + hints.decorations = 0; /* Absolutely no decorations. */ + // other hints.decorations make no sense + // you can't select individual decorations + + XChangeProperty(m_display, m_window, + atom, atom, 32, + PropModeReplace, (unsigned char *) &hints, 4); + } + } else if (state == GHOST_kWindowStateMaximized) { + // With this, xprop should report the following just after launch + // _NET_WM_STATE(ATOM) = _NET_WM_STATE_MAXIMIZED_VERT, _NET_WM_STATE_MAXIMIZED_HORZ + // After demaximization the right side is empty, though (maybe not the most correct then?) + Atom state, atomh, atomv; + + state = XInternAtom(m_display, "_NET_WM_STATE", False); + atomh = XInternAtom(m_display, "_NET_WM_STATE_MAXIMIZED_HORZ", False); + atomv = XInternAtom(m_display, "_NET_WM_STATE_MAXIMIZED_VERT", False); + if (state == None ) { + GHOST_PRINT("Atom _NET_WM_STATE requested but not avaliable nor created.\n"); + } else { + XChangeProperty(m_display, m_window, + state, XA_ATOM, 32, + PropModeAppend, (unsigned char *) &atomh, 1); + XChangeProperty(m_display, m_window, + state, XA_ATOM, 32, + PropModeAppend, (unsigned char *) &atomv, 1); + } + } + // Create some hints for the window manager on how // we want this window treated. @@ -654,298 +664,28 @@ clientToScreen( outY = ay; } -void GHOST_WindowX11::icccmSetState(int state) -{ - XEvent xev; - - if (state != IconicState) - return; - - xev.xclient.type = ClientMessage; - xev.xclient.serial = 0; - xev.xclient.send_event = True; - xev.xclient.display = m_display; - xev.xclient.window = m_window; - xev.xclient.format = 32; - xev.xclient.message_type = m_system->m_wm_change_state; - xev.xclient.data.l[0] = state; - XSendEvent (m_display, RootWindow(m_display, DefaultScreen(m_display)), - False, SubstructureNotifyMask | SubstructureRedirectMask, &xev); -} - -int GHOST_WindowX11::icccmGetState(void) const -{ - unsigned char *prop_ret; - unsigned long bytes_after, num_ret; - Atom type_ret; - int format_ret, st; - - prop_ret = NULL; - st = XGetWindowProperty(m_display, m_window, m_system->m_wm_state, 0, - 0x7fffffff, False, m_system->m_wm_state, &type_ret, - &format_ret, &num_ret, &bytes_after, &prop_ret); - - if ((st == Success) && (prop_ret) && (num_ret == 2)) - st = prop_ret[0]; - else - st = NormalState; - - if (prop_ret) - XFree(prop_ret); - return (st); -} - -void GHOST_WindowX11::netwmMaximized(bool set) -{ - XEvent xev; - - xev.xclient.type= ClientMessage; - xev.xclient.serial = 0; - xev.xclient.send_event = True; - xev.xclient.window = m_window; - xev.xclient.message_type = m_system->m_net_state; - xev.xclient.format = 32; - - if (set == True) - xev.xclient.data.l[0] = _NET_WM_STATE_ADD; - else - xev.xclient.data.l[0] = _NET_WM_STATE_REMOVE; - - xev.xclient.data.l[1] = m_system->m_net_max_horz; - xev.xclient.data.l[2] = m_system->m_net_max_vert; - xev.xclient.data.l[3] = 0; - xev.xclient.data.l[4] = 0; - XSendEvent (m_display, RootWindow(m_display, DefaultScreen(m_display)), - False, SubstructureRedirectMask | SubstructureNotifyMask, &xev); -} - -bool GHOST_WindowX11::netwmIsMaximized(void) const -{ - unsigned char *prop_ret; - unsigned long bytes_after, num_ret, i; - Atom type_ret; - bool st; - int format_ret, count; - - prop_ret = NULL; - st = False; - i = XGetWindowProperty(m_display, m_window, m_system->m_net_state, 0, - 0x7fffffff, False, XA_ATOM, &type_ret, &format_ret, - &num_ret, &bytes_after, &prop_ret); - if ((i == Success) && (prop_ret) && (format_ret == 32)) { - count = 0; - for (i = 0; i < num_ret; i++) { - if (((unsigned long *) prop_ret)[i] == m_system->m_net_max_horz) - count++; - if (((unsigned long *) prop_ret)[i] == m_system->m_net_max_vert) - count++; - if (count == 2) { - st = True; - break; - } - } - } - - if (prop_ret) - XFree(prop_ret); - return (st); -} - -void GHOST_WindowX11::netwmFullScreen(bool set) -{ - XEvent xev; - - xev.xclient.type = ClientMessage; - xev.xclient.serial = 0; - xev.xclient.send_event = True; - xev.xclient.window = m_window; - xev.xclient.message_type = m_system->m_net_state; - xev.xclient.format = 32; - - if (set == True) - xev.xclient.data.l[0] = _NET_WM_STATE_ADD; - else - xev.xclient.data.l[0] = _NET_WM_STATE_REMOVE; - - xev.xclient.data.l[1] = m_system->m_net_fullscreen; - xev.xclient.data.l[2] = 0; - xev.xclient.data.l[3] = 0; - xev.xclient.data.l[4] = 0; - XSendEvent (m_display, RootWindow(m_display, DefaultScreen(m_display)), - False, SubstructureRedirectMask | SubstructureNotifyMask, &xev); -} - -bool GHOST_WindowX11::netwmIsFullScreen(void) const -{ - unsigned char *prop_ret; - unsigned long bytes_after, num_ret, i; - Atom type_ret; - bool st; - int format_ret; - - prop_ret = NULL; - st = False; - i = XGetWindowProperty(m_display, m_window, m_system->m_net_state, 0, - 0x7fffffff, False, XA_ATOM, &type_ret, &format_ret, - &num_ret, &bytes_after, &prop_ret); - if ((i == Success) && (prop_ret) && (format_ret == 32)) { - for (i = 0; i < num_ret; i++) { - if (((unsigned long *)prop_ret)[i] == m_system->m_net_fullscreen) { - st = True; - break; - } - } - } - - if (prop_ret) - XFree(prop_ret); - return (st); -} - -void GHOST_WindowX11::motifFullScreen(bool set) -{ - MotifWmHints hints; - - hints.flags = MWM_HINTS_DECORATIONS; - if (set == True) - hints.decorations = 0; - else - hints.decorations = 1; - XChangeProperty(m_display, m_window, m_system->m_motif, - m_system->m_motif, 32, PropModeReplace, - (unsigned char *)&hints, 4); -} - -bool GHOST_WindowX11::motifIsFullScreen(void) const -{ - unsigned char *prop_ret; - unsigned long bytes_after, num_ret; - MotifWmHints *hints; - Atom type_ret; - bool state; - int format_ret, st; - - prop_ret = NULL; - state = False; - st = XGetWindowProperty(m_display, m_window, m_system->m_motif, 0, - 0x7fffffff, False, m_system->m_motif, - &type_ret, &format_ret, &num_ret, - &bytes_after, &prop_ret); - if ((st == Success) && (prop_ret)) { - hints = (MotifWmHints *)prop_ret; - if (hints->flags & MWM_HINTS_DECORATIONS) { - if (!hints->decorations) - state = True; - } - } - - if (prop_ret) - XFree(prop_ret); - return (state); -} - -GHOST_TWindowState GHOST_WindowX11::getState() const -{ - GHOST_TWindowState state_ret; - int state; - - state_ret = GHOST_kWindowStateNormal; - state = icccmGetState(); - /* - * In the Iconic and Withdrawn state, the window is - * unmaped, so only need return a Minimized state. - */ - if ((state == IconicState) || (state == WithdrawnState)) - state_ret = GHOST_kWindowStateMinimized; - else if (netwmIsMaximized() == True) - state_ret = GHOST_kWindowStateMaximized; - else if (netwmIsFullScreen() == True) - state_ret = GHOST_kWindowStateFullScreen; - else if (motifIsFullScreen() == True) - state_ret = GHOST_kWindowStateFullScreen; - return (state_ret); + GHOST_TWindowState +GHOST_WindowX11:: +getState( +) const { + //FIXME + return GHOST_kWindowStateNormal; } -GHOST_TSuccess GHOST_WindowX11::setState(GHOST_TWindowState state) -{ - GHOST_TWindowState cur_state; - bool is_max, is_full, is_motif_full; + GHOST_TSuccess +GHOST_WindowX11:: +setState( + GHOST_TWindowState state +){ + //TODO - cur_state = getState(); - if (state == (int)cur_state) + if (state == (int)getState()) { return GHOST_kSuccess; - - if (cur_state != GHOST_kWindowStateMinimized) { - /* - * The window don't have this property's - * if it's not mapped. - */ - is_max = netwmIsMaximized(); - is_full = netwmIsFullScreen(); - } - else { - is_max = False; - is_full = False; - } - - is_motif_full = motifIsFullScreen(); - - if (state == GHOST_kWindowStateNormal) { - if (is_max == True) - netwmMaximized(False); - if (is_full == True) - netwmFullScreen(False); - if (is_motif_full == True) - motifFullScreen(False); - icccmSetState(NormalState); - return (GHOST_kSuccess); - } - - if (state == GHOST_kWindowStateFullScreen) { - /* - * We can't change to full screen if the window - * isn't mapped. - */ - if (cur_state == GHOST_kWindowStateMinimized) - return (GHOST_kFailure); - - if (is_max == True) - netwmMaximized(False); - if (is_full == False) - netwmFullScreen(True); - if (is_motif_full == False) - motifFullScreen(True); - return (GHOST_kSuccess); - } - - if (state == GHOST_kWindowStateMaximized) { - /* - * We can't change to Maximized if the window - * isn't mapped. - */ - if (cur_state == GHOST_kWindowStateMinimized) - return (GHOST_kFailure); - - if (is_full == True) - netwmFullScreen(False); - if (is_motif_full == True) - motifFullScreen(False); - if (is_max == False) - netwmMaximized(True); - return (GHOST_kSuccess); - } - - if (state == GHOST_kWindowStateMinimized) { - /* - * The window manager need save the current state of - * the window (maximized, full screen, etc). - */ - icccmSetState(IconicState); - return (GHOST_kSuccess); + } else { + return GHOST_kFailure; } - return (GHOST_kFailure); } #include <iostream> diff --git a/intern/ghost/intern/GHOST_WindowX11.h b/intern/ghost/intern/GHOST_WindowX11.h index 1d73faaf3f6..abb5c131cb7 100644 --- a/intern/ghost/intern/GHOST_WindowX11.h +++ b/intern/ghost/intern/GHOST_WindowX11.h @@ -212,15 +212,6 @@ public: const GHOST_TabletData* GetTabletData() { return &m_xtablet.CommonData; } - - /* - * Need this in case that we want start the window - * in FullScreen or Maximized state. - * Check GHOST_WindowX11.cpp - */ - bool m_post_init; - GHOST_TWindowState m_post_state; - protected: /** * Tries to install a rendering context in this window. @@ -336,18 +327,6 @@ private : /* Tablet devices */ XTablet m_xtablet; - - void icccmSetState(int state); - int icccmGetState() const; - - void netwmMaximized(bool set); - bool netwmIsMaximized() const; - - void netwmFullScreen(bool set); - bool netwmIsFullScreen() const; - - void motifFullScreen(bool set); - bool motifIsFullScreen() const; }; |