forked from bartvdbraak/blender
This version now includes the fluid control sources, however the Blender
interface for it is still missing. Right now there is only a simple hard coded example, that moves a single control particle with strong attraction and velocity forces through the domain. I added more detailed information about the control code to the wiki http://wiki.blender.org/index.php/SoCFluidDevelDoc#The_Fluid-Control_Branch , together with some thoughts on how a Blender integration could be done.
This commit is contained in:
Nils Thuerey
parent
5606cb156e
commit
006d4d1176
@ -31,7 +31,7 @@ SET(INC ${PNG_INC} ${ZLIB_INC} ${SDL_INC})
|
||||
|
||||
FILE(GLOB SRC intern/*.cpp)
|
||||
|
||||
ADD_DEFINITIONS(-DNOGUI -DELBEEM_BLENDER=1)
|
||||
ADD_DEFINITIONS(-DNOGUI -DELBEEM_BLENDER=1 -DLBM_INCLUDE_CONTROL=1)
|
||||
IF(WINDOWS)
|
||||
ADD_DEFINITIONS(-DUSE_MSVC6FIXES)
|
||||
ENDIF(WINDOWS)
|
||||
|
||||
@ -5,7 +5,7 @@ Import('env')
|
||||
|
||||
sources = env.Glob('intern/*.cpp')
|
||||
|
||||
defs = ' NOGUI ELBEEM_BLENDER=1'
|
||||
defs = 'NOGUI ELBEEM_BLENDER=1 LBM_INCLUDE_CONTROL=1'
|
||||
|
||||
if env['WITH_BF_OPENMP'] == 1:
|
||||
defs += ' PARALLEL'
|
||||
|
||||
1320
intern/elbeem/intern/controlparticles.cpp
Normal file
1320
intern/elbeem/intern/controlparticles.cpp
Normal file
File diff suppressed because it is too large
Load Diff
297
intern/elbeem/intern/controlparticles.h
Normal file
297
intern/elbeem/intern/controlparticles.h
Normal file
@ -0,0 +1,297 @@
|
||||
// --------------------------------------------------------------------------
|
||||
//
|
||||
// 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
|
||||
|
||||
// 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();
|
||||
|
||||
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
|
||||
|
||||
@ -30,6 +30,7 @@ int guiRoiMaxLev=6, guiRoiMinLev=0;
|
||||
ntlWorld *gpWorld = NULL;
|
||||
|
||||
|
||||
|
||||
// API
|
||||
|
||||
// reset elbeemSimulationSettings struct with defaults
|
||||
|
||||
22
intern/elbeem/intern/elbeem_control.cpp
Normal file
22
intern/elbeem/intern/elbeem_control.cpp
Normal file
@ -0,0 +1,22 @@
|
||||
/******************************************************************************
|
||||
*
|
||||
* 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;
|
||||
}
|
||||
|
||||
62
intern/elbeem/intern/elbeem_control.h
Normal file
62
intern/elbeem/intern/elbeem_control.h
Normal file
@ -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
|
||||
@ -13,10 +13,6 @@
|
||||
#include <algorithm>
|
||||
#include <stdio.h>
|
||||
|
||||
#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
|
||||
#include <ieeefp.h>
|
||||
#endif
|
||||
|
||||
// just use default rounding for platforms where its not available
|
||||
#ifndef round
|
||||
#define round(x) (x)
|
||||
|
||||
191
intern/elbeem/intern/mvmcoords.cpp
Normal file
191
intern/elbeem/intern/mvmcoords.cpp
Normal file
@ -0,0 +1,191 @@
|
||||
/******************************************************************************
|
||||
*
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
77
intern/elbeem/intern/mvmcoords.h
Normal file
77
intern/elbeem/intern/mvmcoords.h
Normal file
@ -0,0 +1,77 @@
|
||||
/******************************************************************************
|
||||
*
|
||||
// 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
|
||||
|
||||
// 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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -105,6 +105,10 @@
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#if LBM_INCLUDE_CONTROL==1
|
||||
#include "solver_control.h"
|
||||
#endif
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
#include "solver_test.h"
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
@ -497,6 +501,14 @@ class LbmFsgrSolver :
|
||||
LbmFloat& debRAC(LbmFloat* s,int l);
|
||||
# endif // FSGR_STRICT_DEBUG==1
|
||||
|
||||
# if LBM_INCLUDE_CONTROL==1
|
||||
LbmControlData *mpControl;
|
||||
|
||||
void initCpdata();
|
||||
void handleCpdata();
|
||||
void cpDebugDisplay(int dispset);
|
||||
# endif // LBM_INCLUDE_CONTROL==1
|
||||
|
||||
bool mUseTestdata;
|
||||
# if LBM_INCLUDE_TESTSOLVERS==1
|
||||
// test functions
|
||||
@ -506,10 +518,6 @@ class LbmFsgrSolver :
|
||||
void handleTestdata();
|
||||
void set3dHeight(int ,int );
|
||||
|
||||
void initCpdata();
|
||||
void handleCpdata();
|
||||
void cpDebugDisplay(int dispset);
|
||||
|
||||
int mMpNum,mMpIndex;
|
||||
int mOrgSizeX;
|
||||
LbmFloat mOrgStartX;
|
||||
|
||||
961
intern/elbeem/intern/solver_control.cpp
Normal file
961
intern/elbeem/intern/solver_control.cpp
Normal file
@ -0,0 +1,961 @@
|
||||
/******************************************************************************
|
||||
*
|
||||
* 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"
|
||||
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
* 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(16), mCpOutfile(""),
|
||||
mCpForces(), mCpKernel(), mMdKernel(),
|
||||
mDiffVelCon(1.),
|
||||
mDebugCpscale(0.),
|
||||
mDebugVelScale(0.),
|
||||
mDebugCompavScale(0.),
|
||||
mDebugAttScale(0.),
|
||||
mDebugMaxdScale(0.),
|
||||
mDebugAvgVelScale(0.)
|
||||
{
|
||||
}
|
||||
|
||||
LbmControlData::~LbmControlData()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
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; }
|
||||
|
||||
// NT blender integration manual test setup
|
||||
if(1) {
|
||||
// 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 = 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 = mpControl->mCons[cpssi]->mCparts->getInfluenceAttraction();
|
||||
LbmFloat ivel2 = mpControl->mCons[cpssi]->mCparts->getInfluenceVelocity();
|
||||
LbmFloat imaxd2 = mpControl->mCons[cpssi]->mCparts->getInfluenceMaxdist();
|
||||
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;
|
||||
|
||||
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
|
||||
|
||||
|
||||
187
intern/elbeem/intern/solver_control.h
Normal file
187
intern/elbeem/intern/solver_control.h
Normal file
@ -0,0 +1,187 @@
|
||||
/******************************************************************************
|
||||
*
|
||||
* El'Beem - the visual lattice boltzmann freesurface simulator
|
||||
* All code distributed as part of El'Beem is covered by the version 2 of the
|
||||
* GNU General Public License. See the file COPYING for details.
|
||||
* Copyright 2003-2006 Nils Thuerey
|
||||
*
|
||||
* testing extensions
|
||||
*
|
||||
*****************************************************************************/
|
||||
|
||||
|
||||
#ifndef LBM_TESTCLASS_H
|
||||
#define LBM_TESTCLASS_H
|
||||
|
||||
//class IsoSurface;
|
||||
class ParticleObject;
|
||||
class ControlParticles;
|
||||
class ControlForces;
|
||||
|
||||
//#define NUMGRIDS 2
|
||||
//#define MAXNUMSWS 10
|
||||
|
||||
// farfield modes
|
||||
#define FARF_3DONLY -1
|
||||
#define FARF_BOTH 0
|
||||
#define FARF_SWEONLY 1
|
||||
// dont reuse 3d vars/init
|
||||
#define FARF_SEPSWE 2
|
||||
|
||||
// relaxation macros for solver_relax.h
|
||||
#if LBM_INCLUDE_CONTROL!=1
|
||||
|
||||
// defined in relax.h
|
||||
|
||||
#else // LBM_INCLUDE_TESTSOLVERS!=1
|
||||
|
||||
// WARNING has to match controlparts.h
|
||||
#define CPF_ENTRIES 12
|
||||
#define CPF_FORCE 0
|
||||
#define CPF_VELWEIGHT 3
|
||||
#define CPF_VELOCITY 4
|
||||
#define CPF_FORCEWEIGHT 7
|
||||
#define CPF_MINCPDIST 8
|
||||
#define CPF_MINCPDIR 9
|
||||
|
||||
#include "controlparticles.h"
|
||||
|
||||
// get force entry, set=0 is unused anyway
|
||||
#define LBMGET_FORCE(lev, i,j,k) mpControl->mCpForces[lev][ (LBMGI(lev,i,j,k,0)) ]
|
||||
|
||||
// debug mods off...
|
||||
// same as in src/solver_relax.h!
|
||||
#define __PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
|
||||
ux += (grav)[0]; \
|
||||
uy += (grav)[1]; \
|
||||
uz += (grav)[2];
|
||||
|
||||
//void testMaxdmod(int i, int j,int k, LbmFloat &ux,LbmFloat &uy,LbmFloat &uz,ControlForces &ff);
|
||||
#if LBMDIM==3
|
||||
#define MAXDGRAV \
|
||||
if(myforce->forceMaxd[0]*ux+myforce->forceMaxd[1]*uy<LBM_EPSILON) { \
|
||||
ux = v2w*myforce->forceVel[0]+ v2wi*ux; \
|
||||
uy = v2w*myforce->forceVel[1]+ v2wi*uy; } \
|
||||
/* movement inverse to g? */ \
|
||||
if((uz>LBM_EPSILON)&&(uz>myforce->forceVel[2])) { \
|
||||
uz = v2w*myforce->forceVel[2]+ v2wi*uz; }
|
||||
#else // LBMDIM==3
|
||||
#define MAXDGRAV \
|
||||
if(myforce->forceMaxd[0]*ux<LBM_EPSILON) { \
|
||||
ux = v2w*myforce->forceVel[0]+ v2wi*ux; } \
|
||||
/* movement inverse to g? */ \
|
||||
if((uy>LBM_EPSILON)&&(uy>myforce->forceVel[1])) { \
|
||||
uy = v2w*myforce->forceVel[1]+ v2wi*uy; }
|
||||
#endif // LBMDIM==3
|
||||
|
||||
// debug modifications of collide vars (testing)
|
||||
// requires: lev,i,j,k
|
||||
#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
|
||||
LbmFloat attforce = 1.; \
|
||||
if(this->mTForceStrength>0.) { \
|
||||
ControlForces* myforce = &LBMGET_FORCE(lev,i,j,k); \
|
||||
const LbmFloat vf = myforce->weightAtt;\
|
||||
const LbmFloat vw = myforce->weightVel;\
|
||||
if(vf!=0.) { attforce = MAX(0., 1.-vf); /* TODO FIXME? use ABS(vf) for repulsion force? */ \
|
||||
ux += myforce->forceAtt[0]; \
|
||||
uy += myforce->forceAtt[1]; \
|
||||
uz += myforce->forceAtt[2]; \
|
||||
\
|
||||
} else if(( myforce->maxDistance>0.) && ( myforce->maxDistance<CPF_MAXDINIT)) {\
|
||||
const LbmFloat v2w = mpControl->mCons[0]->mCparts->getInfluenceMaxdist() * \
|
||||
(myforce->maxDistance-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()) / (mpControl->mCons[0]->mCparts->getRadiusMaxd()-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()); \
|
||||
const LbmFloat v2wi = 1.-v2w; \
|
||||
if(v2w>0.){ MAXDGRAV; \
|
||||
/* errMsg("ERRMDTT","at "<<PRINT_IJK<<" maxd="<<myforce->maxDistance<<", newu"<<PRINT_VEC(ux,uy,uz)<<", org"<<PRINT_VEC(oux,ouy,ouz)<<", fv"<<myforce->forceVel<<" " ); */ \
|
||||
}\
|
||||
} \
|
||||
if(vw>0.) { \
|
||||
const LbmFloat vwi = 1.-vw;\
|
||||
const LbmFloat vwd = mpControl->mDiffVelCon;\
|
||||
ux += vw*(myforce->forceVel[0]-myforce->compAv[0] + vwd*(myforce->compAv[0]-ux) ); \
|
||||
uy += vw*(myforce->forceVel[1]-myforce->compAv[1] + vwd*(myforce->compAv[1]-uy) ); \
|
||||
uz += vw*(myforce->forceVel[2]-myforce->compAv[2] + vwd*(myforce->compAv[2]-uz) ); \
|
||||
/* TODO test!? modify smooth vel by influence of force for each lbm step, to account for force update only each N steps */ \
|
||||
myforce->compAv = (myforce->forceVel*vw+ myforce->compAv*vwi); \
|
||||
} \
|
||||
} \
|
||||
ux += (grav)[0]*attforce; \
|
||||
uy += (grav)[1]*attforce; \
|
||||
uz += (grav)[2]*attforce; \
|
||||
/* end PRECOLLIDE_MODS */
|
||||
|
||||
#define TEST_IF_CHECK \
|
||||
if((!iffilled)&&(LBMGET_FORCE(lev,i,j,k).weightAtt!=0.)) { \
|
||||
errMsg("TESTIFFILL"," at "<<PRINT_IJK<<" "<<mass<<" "<<rho); \
|
||||
iffilled = true; \
|
||||
if(mass<rho*1.0) mass = rho*1.0; myfrac = 1.0; \
|
||||
}
|
||||
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS!=1
|
||||
|
||||
|
||||
// a single set of control particles and params
|
||||
class LbmControlSet {
|
||||
public:
|
||||
LbmControlSet();
|
||||
~LbmControlSet();
|
||||
void initCparts();
|
||||
|
||||
// control particles
|
||||
ControlParticles *mCparts;
|
||||
// control particle overall motion (for easier manual generation)
|
||||
ControlParticles *mCpmotion;
|
||||
// cp data file
|
||||
string mContrPartFile;
|
||||
string mCpmotionFile;
|
||||
// cp debug displau
|
||||
LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
|
||||
|
||||
// params
|
||||
AnimChannel<float> mcForceAtt;
|
||||
AnimChannel<float> mcForceVel;
|
||||
AnimChannel<float> mcForceMaxd;
|
||||
|
||||
AnimChannel<float> mcRadiusAtt;
|
||||
AnimChannel<float> mcRadiusVel;
|
||||
AnimChannel<float> mcRadiusMind;
|
||||
AnimChannel<float> mcRadiusMaxd;
|
||||
|
||||
AnimChannel<ntlVec3f> mcCpScale;
|
||||
AnimChannel<ntlVec3f> mcCpOffset;
|
||||
};
|
||||
|
||||
|
||||
|
||||
// main control data storage
|
||||
class LbmControlData
|
||||
{
|
||||
public:
|
||||
LbmControlData();
|
||||
virtual ~LbmControlData();
|
||||
|
||||
// control data
|
||||
|
||||
// contorl params
|
||||
void parseControldataAttrList(AttributeList *attr);
|
||||
|
||||
// control strength, set for solver interface
|
||||
LbmFloat mSetForceStrength;
|
||||
// cp vars
|
||||
std::vector<LbmControlSet*> mCons;
|
||||
// update interval
|
||||
int mCpUpdateInterval;
|
||||
// output
|
||||
string mCpOutfile;
|
||||
// control particle precomputed influence
|
||||
std::vector< std::vector<ControlForces> > mCpForces;
|
||||
std::vector<ControlForces> mCpKernel;
|
||||
std::vector<ControlForces> mMdKernel;
|
||||
// activate differential velcon
|
||||
LbmFloat mDiffVelCon;
|
||||
|
||||
// cp debug displau
|
||||
LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
|
||||
};
|
||||
|
||||
#endif // LBM_TESTCLASS_H
|
||||
@ -328,7 +328,10 @@ LbmFsgrSolver::LbmFsgrSolver() :
|
||||
mInit2dYZ(false),
|
||||
mForceTadapRefine(-1), mCutoff(-1)
|
||||
{
|
||||
// not much to do here...
|
||||
#if LBM_INCLUDE_CONTROL==1
|
||||
mpControl = new LbmControlData();
|
||||
#endif
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
mpTest = new LbmTestdata();
|
||||
mMpNum = mMpIndex = 0;
|
||||
@ -488,6 +491,10 @@ void LbmFsgrSolver::parseAttrList()
|
||||
mSimulationTime += starttimeskip;
|
||||
if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
|
||||
|
||||
#if LBM_INCLUDE_CONTROL==1
|
||||
mpControl->parseControldataAttrList(mpSifAttrs);
|
||||
#endif
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
mUseTestdata = 0;
|
||||
mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
|
||||
@ -703,7 +710,7 @@ bool LbmFsgrSolver::initializeSolverMemory()
|
||||
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 +1271,11 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
|
||||
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
|
||||
mInitDone = 1;
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
#if LBM_INCLUDE_CONTROL==1
|
||||
initCpdata();
|
||||
#endif // LBM_INCLUDE_CONTROL==1
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
initTestdata();
|
||||
#endif // ELBEEM_PLUGIN!=1
|
||||
// not inited? dont use...
|
||||
|
||||
@ -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,7 +53,7 @@ void LbmFsgrSolver::stepMain() {
|
||||
|
||||
// init moving bc's, can change mMaxVlen
|
||||
initMovingObstacles(false);
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
#if LBM_INCLUDE_CONTROL==1
|
||||
handleCpdata();
|
||||
#endif
|
||||
|
||||
|
||||
@ -15,7 +15,8 @@
|
||||
#define CAUSE_PANIC { this->mPanic=1; } /*set flag*/
|
||||
#endif // FSGR_STRICT_DEBUG==1
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS!=1
|
||||
// #if LBM_INCLUDE_TESTSOLVERS!=1
|
||||
#if LBM_INCLUDE_CONTROL!=1
|
||||
|
||||
#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
|
||||
ux += (grav)[0]; \
|
||||
@ -24,42 +25,9 @@
|
||||
|
||||
#define TEST_IF_CHECK
|
||||
|
||||
#else // LBM_INCLUDE_TESTSOLVERS!=1
|
||||
// defined in test.h
|
||||
|
||||
#define NEWDIRVELMOTEST 0
|
||||
#if NEWDIRVELMOTEST==1
|
||||
// off for non testing
|
||||
#undef PRECOLLIDE_MODS
|
||||
#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
|
||||
ux += (grav)[0]; \
|
||||
uy += (grav)[1]; \
|
||||
uz += (grav)[2]; \
|
||||
{ \
|
||||
int lev = mMaxRefine, nomb=0; \
|
||||
LbmFloat bcnt = 0.,nux=0.,nuy=0.,nuz=0.; \
|
||||
for(int l=1; l<this->cDfNum; l++) { \
|
||||
if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { \
|
||||
if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBndMoving) { \
|
||||
nux += QCELL_NB(lev, i,j,k,SRCS(lev),l, dMass); \
|
||||
nuy += QCELL_NB(lev, i,j,k,SRCS(lev),l, dFfrac); \
|
||||
bcnt += 1.; \
|
||||
} else { \
|
||||
nomb++; \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
if((bcnt>0.)&&(nomb==0)) { \
|
||||
ux = nux/bcnt; \
|
||||
uy = nuy/bcnt; \
|
||||
uz = nuz/bcnt; \
|
||||
} \
|
||||
}
|
||||
#else // NEWDIRVELMOTEST==1
|
||||
// off for non testing
|
||||
#endif // NEWDIRVELMOTEST==1
|
||||
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS!=1
|
||||
#else // LBM_INCLUDE_CONTROL!=1
|
||||
// defined in solver_control.h
|
||||
#endif // LBM_INCLUDE_CONTROL!=1
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
|
||||
Loading…
Reference in New Issue
Block a user