forked from bartvdbraak/blender
818bfcca63
compilation. Currently, it is not yet active by default, but Genscher wanted to do some tests. It can be used to distribute the computation load onto multiple shared- memory CPUs by splitting the domain along the y-axis (assuming a gravity force along z). However, there is no load balancing: so if there's fluid only in one of the y-axis halves you will not get a speedup for 2 CPUs. - Added a fix for the memory allocation bugs #7120 and #6775. In solver_init.cpp there are now several variables max___MemChunk (line 692+), that set upper limits for various systems. The same problem existed for mac & linux, but the limit is higher, so it probably went by undetected. The windows limit is currently 1GB, if the strange 700MB limit problems mentioned in the bug regports the bugs persist, this could be further reduced. For 64bit compilations this problem shouldn't exist anyway. What's still missing is a display of how much the resolution was reduced to fit into memory... - And some minor solver code cleanup.
742 lines
27 KiB
C++
742 lines
27 KiB
C++
/******************************************************************************
|
|
*
|
|
* 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
|
|
*
|
|
* Combined 2D/3D Lattice Boltzmann Interface Class
|
|
* contains stuff to be statically compiled
|
|
*
|
|
*****************************************************************************/
|
|
|
|
/* LBM Files */
|
|
#include "ntl_matrices.h"
|
|
#include "solver_interface.h"
|
|
#include "ntl_ray.h"
|
|
#include "ntl_world.h"
|
|
#include "elbeem.h"
|
|
|
|
|
|
|
|
|
|
/******************************************************************************
|
|
* Interface Constructor
|
|
*****************************************************************************/
|
|
LbmSolverInterface::LbmSolverInterface() :
|
|
mPanic( false ),
|
|
mSizex(10), mSizey(10), mSizez(10),
|
|
mAllfluid(false), mStepCnt( 0 ),
|
|
mFixMass( 0.0 ),
|
|
mOmega( 1.0 ),
|
|
mGravity(0.0),
|
|
mSurfaceTension( 0.0 ),
|
|
mBoundaryEast( (CellFlagType)(CFBnd) ),mBoundaryWest( (CellFlagType)(CFBnd) ),mBoundaryNorth( (CellFlagType)(CFBnd) ),
|
|
mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop( (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ),
|
|
mInitDone( false ),
|
|
mInitDensityGradient( false ),
|
|
mpSifAttrs( NULL ), mpSifSwsAttrs(NULL), mpParam( NULL ), mpParticles(NULL),
|
|
mNumParticlesLost(0),
|
|
mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0),
|
|
mDebugVelScale( 0.01 ), mNodeInfoString("+"),
|
|
mvGeoStart(-1.0), mvGeoEnd(1.0), mpSimTrafo(NULL),
|
|
mAccurateGeoinit(0),
|
|
mName("lbm_default") ,
|
|
mpIso( NULL ), mIsoValue(0.499),
|
|
mSilent(false) ,
|
|
mLbmInitId(1) ,
|
|
mpGiTree( NULL ),
|
|
mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ),
|
|
mRefinementDesired(0),
|
|
mOutputSurfacePreview(0), mPreviewFactor(0.25),
|
|
mSmoothSurface(1.0), mSmoothNormals(1.0),
|
|
mIsoSubdivs(1), mPartGenProb(0.),
|
|
mDumpVelocities(false),
|
|
mMarkedCells(), mMarkedCellIndex(0),
|
|
mDomainBound("noslip"), mDomainPartSlipValue(0.1),
|
|
mFarFieldSize(0.),
|
|
mPartDropMassSub(0.1), // good default
|
|
mPartUsePhysModel(false),
|
|
mTForceStrength(0.0),
|
|
mCppfStage(0),
|
|
mDumpRawText(false),
|
|
mDumpRawBinary(false),
|
|
mDumpRawBinaryZip(true)
|
|
{
|
|
#if ELBEEM_PLUGIN==1
|
|
if(gDebugLevel<=1) setSilent(true);
|
|
#endif
|
|
mpSimTrafo = new ntlMat4Gfx(0.0);
|
|
mpSimTrafo->initId();
|
|
|
|
if(getenv("ELBEEM_RAWDEBUGDUMP")) {
|
|
debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2);
|
|
mDumpRawText = true;
|
|
}
|
|
|
|
if(getenv("ELBEEM_BINDEBUGDUMP")) {
|
|
debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2);
|
|
mDumpRawBinary = true;
|
|
}
|
|
}
|
|
|
|
LbmSolverInterface::~LbmSolverInterface()
|
|
{
|
|
if(mpSimTrafo) delete mpSimTrafo;
|
|
}
|
|
|
|
|
|
/******************************************************************************
|
|
* initialize correct grid sizes given a geometric bounding box
|
|
* and desired grid resolutions, all params except maxrefine
|
|
* will be modified
|
|
*****************************************************************************/
|
|
void initGridSizes(int &sizex, int &sizey, int &sizez,
|
|
ntlVec3Gfx &geoStart, ntlVec3Gfx &geoEnd,
|
|
int mMaxRefine, bool parallel)
|
|
{
|
|
// fix size inits to force cubic cells and mult4 level dimensions
|
|
const int debugGridsizeInit = 1;
|
|
if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Called - size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" "<<geoStart<<","<<geoEnd ,10);
|
|
|
|
int maxGridSize = sizex; // get max size
|
|
if(sizey>maxGridSize) maxGridSize = sizey;
|
|
if(sizez>maxGridSize) maxGridSize = sizez;
|
|
LbmFloat maxGeoSize = (geoEnd[0]-geoStart[0]); // get max size
|
|
if((geoEnd[1]-geoStart[1])>maxGeoSize) maxGeoSize = (geoEnd[1]-geoStart[1]);
|
|
if((geoEnd[2]-geoStart[2])>maxGeoSize) maxGeoSize = (geoEnd[2]-geoStart[2]);
|
|
// FIXME better divide max geo size by corresponding resolution rather than max? no prob for rx==ry==rz though
|
|
LbmFloat cellSize = (maxGeoSize / (LbmFloat)maxGridSize);
|
|
if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Start:"<<geoStart<<" End:"<<geoEnd<<" maxS:"<<maxGeoSize<<" maxG:"<<maxGridSize<<" cs:"<<cellSize, 10);
|
|
// force grid sizes according to geom. size, rounded
|
|
sizex = (int) ((geoEnd[0]-geoStart[0]) / cellSize +0.5);
|
|
sizey = (int) ((geoEnd[1]-geoStart[1]) / cellSize +0.5);
|
|
sizez = (int) ((geoEnd[2]-geoStart[2]) / cellSize +0.5);
|
|
// match refinement sizes, round downwards to multiple of 4
|
|
int sizeMask = 0;
|
|
int maskBits = mMaxRefine;
|
|
if(parallel==1) maskBits+=2;
|
|
for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); }
|
|
|
|
// at least size 4 on coarsest level
|
|
int minSize = 4<<mMaxRefine; //(maskBits+2);
|
|
if(sizex<minSize) sizex = minSize;
|
|
if(sizey<minSize) sizey = minSize;
|
|
if(sizez<minSize) sizez = minSize;
|
|
|
|
sizeMask = ~sizeMask;
|
|
if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" m"<<convertFlags2String(sizeMask)<<", maskBits:"<<maskBits<<",minSize:"<<minSize ,10);
|
|
sizex &= sizeMask;
|
|
sizey &= sizeMask;
|
|
sizez &= sizeMask;
|
|
// debug
|
|
int nextinc = (~sizeMask)+1;
|
|
debMsgStd("initGridSizes",DM_MSG,"MPTeffMLtest, curr size:"<<PRINT_VEC(sizex,sizey,sizez)<<", next size:"<<PRINT_VEC(sizex+nextinc,sizey+nextinc,sizez+nextinc) ,10);
|
|
|
|
// force geom size to match rounded/modified grid sizes
|
|
geoEnd[0] = geoStart[0] + cellSize*(LbmFloat)sizex;
|
|
geoEnd[1] = geoStart[1] + cellSize*(LbmFloat)sizey;
|
|
geoEnd[2] = geoStart[2] + cellSize*(LbmFloat)sizez;
|
|
}
|
|
|
|
void calculateMemreqEstimate( int resx,int resy,int resz,
|
|
int refine, float farfield,
|
|
double *reqret, double *reqretFine, string *reqstr) {
|
|
// debug estimation?
|
|
const bool debugMemEst = true;
|
|
// COMPRESSGRIDS define is not available here, make sure it matches
|
|
const bool useGridComp = true;
|
|
// make sure we can handle bid numbers here... all double
|
|
double memCnt = 0.0;
|
|
double ddTotalNum = (double)dTotalNum;
|
|
if(reqretFine) *reqretFine = -1.;
|
|
|
|
double currResx = (double)resx;
|
|
double currResy = (double)resy;
|
|
double currResz = (double)resz;
|
|
double rcellSize = ((currResx*currResy*currResz) *ddTotalNum);
|
|
memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0);
|
|
if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"res:"<<PRINT_VEC(currResx,currResy,currResz)<<" rcellSize:"<<rcellSize<<" mc:"<<memCnt, 10);
|
|
if(!useGridComp) {
|
|
memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
|
|
if(reqretFine) *reqretFine = (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
|
|
if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," no-comp, mc:"<<memCnt, 10);
|
|
} else {
|
|
double compressOffset = (double)(currResx*currResy*ddTotalNum*2.0);
|
|
memCnt += (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0));
|
|
if(reqretFine) *reqretFine = (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0));
|
|
if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," w-comp, mc:"<<memCnt, 10);
|
|
}
|
|
for(int i=refine-1; i>=0; i--) {
|
|
currResx /= 2.0;
|
|
currResy /= 2.0;
|
|
currResz /= 2.0;
|
|
rcellSize = ((currResz*currResy*currResx) *ddTotalNum);
|
|
memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0);
|
|
memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
|
|
if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"refine "<<i<<", mc:"<<memCnt, 10);
|
|
}
|
|
|
|
// isosurface memory, use orig res values
|
|
if(farfield>0.) {
|
|
memCnt += (double)( (3*sizeof(int)+sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) );
|
|
} else {
|
|
// ignore 3 int slices...
|
|
memCnt += (double)( ( sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) );
|
|
}
|
|
if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"iso, mc:"<<memCnt, 10);
|
|
|
|
// cpdata init check missing...
|
|
|
|
double memd = memCnt;
|
|
const char *sizeStr = "";
|
|
const double sfac = 1024.0;
|
|
if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
|
|
if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
|
|
if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
|
|
if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
|
|
|
|
// return values
|
|
std::ostringstream ret;
|
|
if(memCnt< 1024.0*1024.0) {
|
|
// show full MBs
|
|
ret << (ceil(memd));
|
|
} else {
|
|
// two digits for anything larger than MB
|
|
ret << (ceil(memd*100.0)/100.0);
|
|
}
|
|
ret << " "<< sizeStr;
|
|
*reqret = memCnt;
|
|
*reqstr = ret.str();
|
|
//debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Required Grid memory: "<< memd <<" "<< sizeStr<<" ",4);
|
|
}
|
|
|
|
void LbmSolverInterface::initDomainTrafo(float *mat) {
|
|
mpSimTrafo->initArrayCheck(mat);
|
|
}
|
|
|
|
/*******************************************************************************/
|
|
/*! parse a boundary flag string */
|
|
CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed) {
|
|
string val = mpSifAttrs->readString(name, "", source, target, needed);
|
|
if(!strcmp(val.c_str(),"")) {
|
|
// no value given...
|
|
return CFEmpty;
|
|
}
|
|
if(!strcmp(val.c_str(),"bnd_no")) {
|
|
return (CellFlagType)( CFBnd );
|
|
}
|
|
if(!strcmp(val.c_str(),"bnd_free")) {
|
|
return (CellFlagType)( CFBnd );
|
|
}
|
|
if(!strcmp(val.c_str(),"fluid")) {
|
|
/* might be used for some in/out flow cases */
|
|
return (CellFlagType)( CFFluid );
|
|
}
|
|
errMsg("LbmSolverInterface::readBoundaryFlagInt","Invalid value '"<<val<<"' " );
|
|
# if FSGR_STRICT_DEBUG==1
|
|
errFatal("readBoundaryFlagInt","Strict abort..."<<val, SIMWORLD_INITERROR);
|
|
# endif
|
|
return defaultValue;
|
|
}
|
|
|
|
/*******************************************************************************/
|
|
/*! parse standard attributes */
|
|
void LbmSolverInterface::parseStdAttrList() {
|
|
if(!mpSifAttrs) {
|
|
errFatal("LbmSolverInterface::parseAttrList","mpSifAttrs pointer not initialized!",SIMWORLD_INITERROR);
|
|
return; }
|
|
|
|
// st currently unused
|
|
mBoundaryEast = readBoundaryFlagInt("boundary_east", mBoundaryEast, "LbmSolverInterface", "mBoundaryEast", false);
|
|
mBoundaryWest = readBoundaryFlagInt("boundary_west", mBoundaryWest, "LbmSolverInterface", "mBoundaryWest", false);
|
|
mBoundaryNorth = readBoundaryFlagInt("boundary_north", mBoundaryNorth,"LbmSolverInterface", "mBoundaryNorth", false);
|
|
mBoundarySouth = readBoundaryFlagInt("boundary_south", mBoundarySouth,"LbmSolverInterface", "mBoundarySouth", false);
|
|
mBoundaryTop = readBoundaryFlagInt("boundary_top", mBoundaryTop,"LbmSolverInterface", "mBoundaryTop", false);
|
|
mBoundaryBottom= readBoundaryFlagInt("boundary_bottom", mBoundaryBottom,"LbmSolverInterface", "mBoundaryBottom", false);
|
|
|
|
mpSifAttrs->readMat4Gfx("domain_trafo" , (*mpSimTrafo), "ntlBlenderDumper","mpSimTrafo", false, mpSimTrafo );
|
|
|
|
LbmVec sizeVec(mSizex,mSizey,mSizez);
|
|
sizeVec = vec2L( mpSifAttrs->readVec3d("size", vec2P(sizeVec), "LbmSolverInterface", "sizeVec", false) );
|
|
mSizex = (int)sizeVec[0];
|
|
mSizey = (int)sizeVec[1];
|
|
mSizez = (int)sizeVec[2];
|
|
// param needs size in any case
|
|
// test solvers might not have mpPara, though
|
|
if(mpParam) mpParam->setSize(mSizex, mSizey, mSizez );
|
|
|
|
mInitDensityGradient = mpSifAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmSolverInterface", "mInitDensityGradient", false);
|
|
mIsoValue = mpSifAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false );
|
|
|
|
mDebugVelScale = mpSifAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmSolverInterface", "mDebugVelScale", false);
|
|
mNodeInfoString = mpSifAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false );
|
|
|
|
mDumpVelocities = mpSifAttrs->readBool("dump_velocities", mDumpVelocities, "SimulationLbm","mDumpVelocities", false );
|
|
if(getenv("ELBEEM_DUMPVELOCITIES")) {
|
|
int get = atoi(getenv("ELBEEM_DUMPVELOCITIES"));
|
|
if((get==0)||(get==1)) {
|
|
mDumpVelocities = get;
|
|
debMsgStd("LbmSolverInterface::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_DUMPVELOCITIES to set mDumpVelocities to "<<get<<","<<mDumpVelocities,8);
|
|
}
|
|
}
|
|
|
|
// new test vars
|
|
// mTForceStrength set from test solver
|
|
mFarFieldSize = mpSifAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false);
|
|
// old compat
|
|
float sizeScale = mpSifAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false);
|
|
if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); }
|
|
|
|
mPartDropMassSub = mpSifAttrs->readFloat("part_dropmass", mPartDropMassSub,"LbmSolverInterface", "mPartDropMassSub", false);
|
|
|
|
mCppfStage = mpSifAttrs->readInt("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false);
|
|
mPartGenProb = mpSifAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false);
|
|
mPartUsePhysModel = mpSifAttrs->readBool("partusephysmodel", mPartUsePhysModel,"LbmFsgrSolver", "mPartUsePhysModel", false);
|
|
mIsoSubdivs = mpSifAttrs->readInt("isosubdivs", mIsoSubdivs,"LbmFsgrSolver", "mIsoSubdivs", false);
|
|
}
|
|
|
|
|
|
/*******************************************************************************/
|
|
/*! geometry initialization */
|
|
/*******************************************************************************/
|
|
|
|
/*****************************************************************************/
|
|
/*! init tree for certain geometry init */
|
|
void LbmSolverInterface::initGeoTree() {
|
|
if(mpGlob == NULL) { errFatal("LbmSolverInterface::initGeoTree","Requires globals!",SIMWORLD_INITERROR); return; }
|
|
ntlScene *scene = mpGlob->getSimScene();
|
|
mpGiObjects = scene->getObjects();
|
|
mGiObjInside.resize( mpGiObjects->size() );
|
|
mGiObjDistance.resize( mpGiObjects->size() );
|
|
mGiObjSecondDist.resize( mpGiObjects->size() );
|
|
for(size_t i=0; i<mpGiObjects->size(); i++) {
|
|
if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true;
|
|
//debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<<mLbmInitId<<" obj:"<< (*mpGiObjects)[i]->getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 );
|
|
}
|
|
debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<<mAccurateGeoinit, 9)
|
|
debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Objects: "<<mpGiObjects->size() ,10);
|
|
|
|
if(mpGiTree != NULL) delete mpGiTree;
|
|
char treeFlag = (1<<(this->mLbmInitId+4));
|
|
mpGiTree = new ntlTree(
|
|
# if FSGR_STRICT_DEBUG!=1
|
|
15, 8, // TREEwarning - fixed values for depth & maxtriangles here...
|
|
# else // FSGR_STRICT_DEBUG==1
|
|
10, 20, // reduced/slower debugging values
|
|
# endif // FSGR_STRICT_DEBUG==1
|
|
scene, treeFlag );
|
|
}
|
|
|
|
/*****************************************************************************/
|
|
/*! destroy tree etc. when geometry init done */
|
|
void LbmSolverInterface::freeGeoTree() {
|
|
if(mpGiTree != NULL) {
|
|
delete mpGiTree;
|
|
mpGiTree = NULL;
|
|
}
|
|
}
|
|
|
|
|
|
int globGeoInitDebug = 0;
|
|
int globGICPIProblems = 0;
|
|
/*****************************************************************************/
|
|
/*! check for a certain flag type at position org */
|
|
bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance, int shootDir) {
|
|
// shift ve ctors to avoid rounding errors
|
|
org += ntlVec3Gfx(0.0001);
|
|
OId = -1;
|
|
|
|
// select shooting direction
|
|
ntlVec3Gfx dir = ntlVec3Gfx(1., 0., 0.);
|
|
if(shootDir==1) dir = ntlVec3Gfx(0., 1., 0.);
|
|
else if(shootDir==2) dir = ntlVec3Gfx(0., 0., 1.);
|
|
else if(shootDir==-2) dir = ntlVec3Gfx(0., 0., -1.);
|
|
else if(shootDir==-1) dir = ntlVec3Gfx(0., -1., 0.);
|
|
|
|
ntlRay ray(org, dir, 0, 1.0, mpGlob);
|
|
bool done = false;
|
|
bool inside = false;
|
|
|
|
vector<int> giObjFirstHistSide;
|
|
giObjFirstHistSide.resize( mpGiObjects->size() );
|
|
for(size_t i=0; i<mGiObjInside.size(); i++) {
|
|
mGiObjInside[i] = 0;
|
|
mGiObjDistance[i] = -1.0;
|
|
mGiObjSecondDist[i] = -1.0;
|
|
giObjFirstHistSide[i] = 0;
|
|
}
|
|
// if not inside, return distance to first hit
|
|
gfxReal firstHit=-1.0;
|
|
int firstOId = -1;
|
|
if(globGeoInitDebug) errMsg("IIIstart"," isect "<<org<<" f"<<flags<<" acc"<<mAccurateGeoinit);
|
|
|
|
if(mAccurateGeoinit) {
|
|
while(!done) {
|
|
// find first inside intersection
|
|
ntlTriangle *triIns = NULL;
|
|
distance = -1.0;
|
|
ntlVec3Gfx normal(0.0);
|
|
if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true);
|
|
else mpGiTree->intersect(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[OId]++;
|
|
if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = 1;
|
|
if(globGeoInitDebug) errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
|
|
} else {
|
|
// inside hit
|
|
mGiObjInside[OId]++;
|
|
if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
|
|
if(globGeoInitDebug) errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
|
|
if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = -1;
|
|
}
|
|
norg += normal * getVecEpsilon();
|
|
ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
|
|
// 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;
|
|
for(size_t i=0; i<mGiObjInside.size(); i++) {
|
|
if(mGiObjInside[i]>0) {
|
|
bool mess = false;
|
|
if((mGiObjInside[i]%2)==1) {
|
|
if(giObjFirstHistSide[i] != -1) mess=true;
|
|
} else {
|
|
if(giObjFirstHistSide[i] != 1) mess=true;
|
|
}
|
|
if(mess) {
|
|
//errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
|
|
globGICPIProblems++;
|
|
mGiObjInside[i]++; // believe first hit side...
|
|
}
|
|
}
|
|
}
|
|
for(size_t i=0; i<mGiObjInside.size(); i++) {
|
|
if(globGeoInitDebug) errMsg("CHIII","i"<<i<<" ins="<<mGiObjInside[i]<<" t"<<mGiObjDistance[i]<<" d"<<distance);
|
|
if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
|
|
if( (distance<0.0) || // first intersection -> good
|
|
((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
|
|
) {
|
|
distance = mGiObjDistance[i];
|
|
OId = i;
|
|
inside = true;
|
|
}
|
|
}
|
|
}
|
|
if(!inside) {
|
|
distance = firstHit;
|
|
OId = firstOId;
|
|
}
|
|
if(globGeoInitDebug) errMsg("CHIII","i"<<inside<<" fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
|
|
|
|
return inside;
|
|
} else {
|
|
|
|
// find first inside intersection
|
|
ntlTriangle *triIns = NULL;
|
|
distance = -1.0;
|
|
ntlVec3Gfx normal(0.0);
|
|
if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true);
|
|
else mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
|
|
if(triIns) {
|
|
// check outside intersect
|
|
LbmFloat orientation = dot(normal, dir);
|
|
if(orientation<=0.0) return false;
|
|
|
|
OId = triIns->getObjectId();
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
}
|
|
bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, const gfxReal halfCellsize, bool &thinHit, bool recurse) {
|
|
// shift ve ctors to avoid rounding errors
|
|
org += ntlVec3Gfx(0.0001); //?
|
|
OId = -1;
|
|
ntlRay ray(org, dir, 0, 1.0, mpGlob);
|
|
//int insCnt = 0;
|
|
bool done = false;
|
|
bool inside = false;
|
|
for(size_t i=0; i<mGiObjInside.size(); i++) {
|
|
mGiObjInside[i] = 0;
|
|
mGiObjDistance[i] = -1.0;
|
|
mGiObjSecondDist[i] = -1.0;
|
|
}
|
|
// if not inside, return distance to first hit
|
|
gfxReal firstHit=-1.0;
|
|
int firstOId = -1;
|
|
thinHit = false;
|
|
|
|
if(mAccurateGeoinit) {
|
|
while(!done) {
|
|
// find first inside intersection
|
|
ntlTriangle *triIns = NULL;
|
|
distance = -1.0;
|
|
ntlVec3Gfx normal(0.0);
|
|
mpGiTree->intersect(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;
|
|
//mGiObjDistance[OId] = -1.0;
|
|
//errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
|
|
} else {
|
|
// inside hit
|
|
//if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
|
|
//errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
|
|
if(mGiObjInside[OId]==1) {
|
|
// second inside hit
|
|
if(mGiObjSecondDist[OId]<0.0) mGiObjSecondDist[OId] = distance;
|
|
}
|
|
}
|
|
mGiObjInside[OId]++;
|
|
// always store first hit for thin obj init
|
|
if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
|
|
|
|
norg += normal * getVecEpsilon();
|
|
ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
|
|
// 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;
|
|
//if(insCnt%2) inside=true;
|
|
}
|
|
}
|
|
|
|
distance = -1.0;
|
|
// standard inside check
|
|
for(size_t i=0; i<mGiObjInside.size(); i++) {
|
|
if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
|
|
if( (distance<0.0) || // first intersection -> good
|
|
((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
|
|
) {
|
|
distance = mGiObjDistance[i];
|
|
OId = i;
|
|
inside = true;
|
|
}
|
|
}
|
|
}
|
|
// now check for thin hits
|
|
if(!inside) {
|
|
distance = -1.0;
|
|
for(size_t i=0; i<mGiObjInside.size(); i++) {
|
|
if((mGiObjInside[i]>=2)&&(mGiObjDistance[i]>0.0)&&(mGiObjSecondDist[i]>0.0)&&
|
|
(mGiObjDistance[i]<1.0*halfCellsize)&&(mGiObjSecondDist[i]<2.0*halfCellsize) ) {
|
|
if( (distance<0.0) || // first intersection -> good
|
|
((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
|
|
) {
|
|
distance = mGiObjDistance[i];
|
|
OId = i;
|
|
inside = true;
|
|
thinHit = true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if(!inside) {
|
|
// check for hit in this cell, opposite to current dir (only recurse once)
|
|
if(recurse) {
|
|
gfxReal r_distance;
|
|
int r_OId;
|
|
bool ret = geoInitCheckPointInside(org, dir*-1.0, flags, r_OId, r_distance, halfCellsize, thinHit, false);
|
|
if((ret)&&(thinHit)) {
|
|
OId = r_OId;
|
|
distance = 0.0;
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
// really no hit...
|
|
if(!inside) {
|
|
distance = firstHit;
|
|
OId = firstOId;
|
|
/*if((mGiObjDistance[OId]>0.0)&&(mGiObjSecondDist[OId]>0.0)) {
|
|
const gfxReal thisdist = mGiObjSecondDist[OId]-mGiObjDistance[OId];
|
|
// dont walk over this cell...
|
|
if(thisdist<halfCellsize) distance-=2.0*halfCellsize;
|
|
} // ? */
|
|
}
|
|
//errMsg("CHIII","i"<<inside<<" fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
|
|
|
|
return inside;
|
|
} else {
|
|
|
|
// find first inside intersection
|
|
ntlTriangle *triIns = NULL;
|
|
distance = -1.0;
|
|
ntlVec3Gfx normal(0.0);
|
|
mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
|
|
if(triIns) {
|
|
// check outside intersect
|
|
LbmFloat orientation = dot(normal, dir);
|
|
if(orientation<=0.0) return false;
|
|
|
|
OId = triIns->getObjectId();
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
/*****************************************************************************/
|
|
ntlVec3Gfx LbmSolverInterface::getGeoMaxMovementVelocity(LbmFloat simtime, LbmFloat stepsize) {
|
|
ntlVec3Gfx max(0.0);
|
|
if(mpGlob == NULL) return max;
|
|
// mpGiObjects has to be inited here...
|
|
|
|
for(int i=0; i< (int)mpGiObjects->size(); i++) {
|
|
//errMsg("LbmSolverInterface::getGeoMaxMovementVelocity","i="<<i<<" "<< (*mpGiObjects)[i]->getName() ); // DEBUG
|
|
if( (*mpGiObjects)[i]->getGeoInitType() & (FGI_FLUID|FGI_MBNDINFLOW) ){
|
|
//ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
|
|
ntlVec3Gfx orgvel = (*mpGiObjects)[i]->calculateMaxVel( simtime, simtime+stepsize );
|
|
if( normNoSqrt(orgvel) > normNoSqrt(max) ) { max = orgvel; }
|
|
//errMsg("MVT","i"<<i<<" a"<< (*mpGiObjects)[i]->getInitialVelocity(simtime)<<" o"<<orgvel ); // DEBUG
|
|
// TODO check if inflow and simtime
|
|
ntlVec3Gfx inivel = (*mpGiObjects)[i]->getInitialVelocity(simtime);
|
|
if( normNoSqrt(inivel) > normNoSqrt(max) ) { max = inivel; }
|
|
}
|
|
}
|
|
errMsg("LbmSolverInterface::getGeoMaxMovementVelocity", "max="<< max ); // DEBUG
|
|
return max;
|
|
}
|
|
|
|
|
|
|
|
|
|
/*******************************************************************************/
|
|
/*! cell iteration functions */
|
|
/*******************************************************************************/
|
|
|
|
|
|
|
|
|
|
/*****************************************************************************/
|
|
//! add cell to mMarkedCells list
|
|
void
|
|
LbmSolverInterface::addCellToMarkedList( CellIdentifierInterface *cid ) {
|
|
for(size_t i=0; i<mMarkedCells.size(); i++) {
|
|
// check if cids alreay in
|
|
if( mMarkedCells[i]->equal(cid) ) return;
|
|
//mMarkedCells[i]->setEnd(false);
|
|
}
|
|
mMarkedCells.push_back( cid );
|
|
//cid->setEnd(true);
|
|
}
|
|
|
|
/*****************************************************************************/
|
|
//! marked cell iteration methods
|
|
CellIdentifierInterface*
|
|
LbmSolverInterface::markedGetFirstCell( ) {
|
|
if(mMarkedCells.size() > 0){ return mMarkedCells[0]; }
|
|
return NULL;
|
|
}
|
|
|
|
CellIdentifierInterface*
|
|
LbmSolverInterface::markedAdvanceCell() {
|
|
mMarkedCellIndex++;
|
|
if(mMarkedCellIndex>=(int)mMarkedCells.size()) return NULL;
|
|
return mMarkedCells[mMarkedCellIndex];
|
|
}
|
|
|
|
void LbmSolverInterface::markedClearList() {
|
|
// FIXME free cids?
|
|
mMarkedCells.clear();
|
|
}
|
|
|
|
/*******************************************************************************/
|
|
/*! string helper functions */
|
|
/*******************************************************************************/
|
|
|
|
|
|
|
|
// 32k
|
|
string convertSingleFlag2String(CellFlagType cflag) {
|
|
CellFlagType flag = cflag;
|
|
if(flag == CFUnused ) return string("cCFUnused");
|
|
if(flag == CFEmpty ) return string("cCFEmpty");
|
|
if(flag == CFBnd ) return string("cCFBnd");
|
|
if(flag == CFBndNoslip ) return string("cCFBndNoSlip");
|
|
if(flag == CFBndFreeslip ) return string("cCFBndFreeSlip");
|
|
if(flag == CFBndPartslip ) return string("cCFBndPartSlip");
|
|
if(flag == CFNoInterpolSrc ) return string("cCFNoInterpolSrc");
|
|
if(flag == CFFluid ) return string("cCFFluid");
|
|
if(flag == CFInter ) return string("cCFInter");
|
|
if(flag == CFNoNbFluid ) return string("cCFNoNbFluid");
|
|
if(flag == CFNoNbEmpty ) return string("cCFNoNbEmpty");
|
|
if(flag == CFNoDelete ) return string("cCFNoDelete");
|
|
if(flag == CFNoBndFluid ) return string("cCFNoBndFluid");
|
|
if(flag == CFGrNorm ) return string("cCFGrNorm");
|
|
if(flag == CFGrFromFine ) return string("cCFGrFromFine");
|
|
if(flag == CFGrFromCoarse ) return string("cCFGrFromCoarse");
|
|
if(flag == CFGrCoarseInited ) return string("cCFGrCoarseInited");
|
|
if(flag == CFMbndInflow ) return string("cCFMbndInflow");
|
|
if(flag == CFMbndOutflow ) return string("cCFMbndOutflow");
|
|
if(flag == CFInvalid ) return string("cfINVALID");
|
|
|
|
std::ostringstream mult;
|
|
int val = 0;
|
|
if(flag != 0) {
|
|
while(! (flag&1) ) {
|
|
flag = flag>>1;
|
|
val++;
|
|
}
|
|
} else {
|
|
val = -1;
|
|
}
|
|
if(val>=24) {
|
|
mult << "cfOID_" << (flag>>24) <<"_TYPE";
|
|
} else {
|
|
mult << "cfUNKNOWN_" << val <<"_TYPE";
|
|
}
|
|
return mult.str();
|
|
}
|
|
|
|
//! helper function to convert flag to string (for debuggin)
|
|
string convertCellFlagType2String( CellFlagType cflag ) {
|
|
int flag = cflag;
|
|
|
|
const int jmax = sizeof(CellFlagType)*8;
|
|
bool somefound = false;
|
|
std::ostringstream mult;
|
|
mult << "[";
|
|
for(int j=0; j<jmax ; j++) {
|
|
if(flag& (1<<j)) {
|
|
if(somefound) mult << "|";
|
|
mult << j<<"<"<< convertSingleFlag2String( (CellFlagType)(1<<j) ); // this call should always be _non_-recursive
|
|
somefound = true;
|
|
}
|
|
};
|
|
mult << "]";
|
|
|
|
// return concatenated string
|
|
if(somefound) return mult.str();
|
|
|
|
// empty?
|
|
return string("[emptyCFT]");
|
|
}
|
|
|