- another minor solver update to fix
authorNils Thuerey <nils@thuerey.de>
Mon, 12 Jun 2006 12:55:51 +0000 (12:55 +0000)
committerNils Thuerey <nils@thuerey.de>
Mon, 12 Jun 2006 12:55:51 +0000 (12:55 +0000)
  obstacle fluid surface generation bug
- also contains some code clean ups
  and safer initializations

13 files changed:
intern/elbeem/intern/attributes.cpp
intern/elbeem/intern/elbeem.cpp
intern/elbeem/intern/ntl_geometrymodel.cpp
intern/elbeem/intern/ntl_geometrymodel.h
intern/elbeem/intern/ntl_ray.h
intern/elbeem/intern/particletracer.cpp
intern/elbeem/intern/simulation_object.cpp
intern/elbeem/intern/solver_class.h
intern/elbeem/intern/solver_init.cpp
intern/elbeem/intern/solver_interface.cpp
intern/elbeem/intern/solver_main.cpp
intern/elbeem/intern/solver_relax.h
intern/elbeem/intern/solver_util.cpp

index f970f6e..b7f1945 100644 (file)
@@ -472,35 +472,45 @@ bool AttributeList::ignoreParameter(string name, string source) {
                
 // read channels
 AnimChannel<int> AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) {
-       if(!exists(name)) { return AnimChannel<int>(defaultValue); } 
+       if(!exists(name)) { 
+               if(needed) { errFatal("AttributeList::readChannelInt","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+               return AnimChannel<int>(defaultValue); } 
        AnimChannel<int> ret = find(name)->getChannelInt(); 
        find(name)->setUsed(true);
        channelSimplifyi(ret);
        return ret;
 }
 AnimChannel<double> AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) {
-       if(!exists(name)) { return AnimChannel<double>(defaultValue); } 
+       if(!exists(name)) {
+               if(needed) { errFatal("AttributeList::readChannelFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+               return AnimChannel<double>(defaultValue); } 
        AnimChannel<double> ret = find(name)->getChannelFloat(); 
        find(name)->setUsed(true);
        channelSimplifyd(ret);
        return ret;
 }
 AnimChannel<ntlVec3d> AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) {
-       if(!exists(name)) { return AnimChannel<ntlVec3d>(defaultValue); } 
+       if(!exists(name)) { 
+               if(needed) { errFatal("AttributeList::readChannelVec3d","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+               return AnimChannel<ntlVec3d>(defaultValue); } 
        AnimChannel<ntlVec3d> ret = find(name)->getChannelVec3d(); 
        find(name)->setUsed(true);
        channelSimplifyVd(ret);
        return ret;
 }
 AnimChannel<ntlSetVec3f> AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) {
-       if(!exists(name)) { return AnimChannel<ntlSetVec3f>(defaultValue); } 
+       if(!exists(name)) { 
+               if(needed) { errFatal("AttributeList::readChannelSetVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+               return AnimChannel<ntlSetVec3f>(defaultValue); } 
        AnimChannel<ntlSetVec3f> ret = find(name)->getChannelSetVec3f(); 
        find(name)->setUsed(true);
        //channelSimplifyVf(ret);
        return ret;
 }
 AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) {
-       if(!exists(name)) { return AnimChannel<float>(defaultValue); } 
+       if(!exists(name)) { 
+               if(needed) { errFatal("AttributeList::readChannelSinglePrecFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+               return AnimChannel<float>(defaultValue); } 
        AnimChannel<double> convert = find(name)->getChannelFloat(); 
        find(name)->setUsed(true);
        channelSimplifyd(convert);
@@ -514,7 +524,9 @@ AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float
        return ret;
 }
 AnimChannel<ntlVec3f> AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) {
-       if(!exists(name)) { return AnimChannel<ntlVec3f>(defaultValue); } 
+       if(!exists(name)) { 
+               if(needed) { errFatal("AttributeList::readChannelVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+               return AnimChannel<ntlVec3f>(defaultValue); } 
 
        AnimChannel<ntlVec3d> convert = find(name)->getChannelVec3d(); 
        // convert to float
@@ -737,6 +749,19 @@ void AnimChannel<Scalar>::debugPrintChannel() {
 }
 
 
+// hack to force instantiation
+void __forceAnimChannelInstantiation() {
+       AnimChannel< float > tmp1;
+       AnimChannel< double > tmp2;
+       AnimChannel< string > tmp3;
+       AnimChannel< ntlVector3Dim<float> > tmp4;
+       tmp1.debugPrintChannel();
+       tmp2.debugPrintChannel();
+       tmp3.debugPrintChannel();
+       tmp4.debugPrintChannel();
+}
+
+
 ntlSetVec3f::ntlSetVec3f(double v ) {
        mVerts.clear();
        mVerts.push_back( ntlVec3f(v) );
index 9074013..d4242da 100644 (file)
@@ -115,32 +115,45 @@ void elbeemGetErrorString(char *buffer) {
 extern "C" 
 void elbeemResetMesh(elbeemMesh *mesh) {
        if(!mesh) return;
+       // init typedef struct elbeemMesh
   mesh->type = 0;
-  mesh->parentDomainId = 0;
+
+       mesh->parentDomainId = 0;
+
+       /* vertices */
   mesh->numVertices = 0;
        mesh->vertices = NULL;
+
+       mesh->channelSizeVertices = 0;
+       mesh->channelVertices = NULL;
+
+       /* triangles */
        mesh->numTriangles = 0;
   mesh->triangles = NULL;
 
+       /* animation channels */
        mesh->channelSizeTranslation = 0;
        mesh->channelTranslation = NULL;
        mesh->channelSizeRotation = 0;
        mesh->channelRotation = NULL;
        mesh->channelSizeScale = 0;
        mesh->channelScale = NULL;
+
+       /* active channel */
        mesh->channelSizeActive = 0;
        mesh->channelActive = NULL;
+
        mesh->channelSizeInitialVel = 0;
        mesh->channelInitialVel = NULL;
+
        mesh->localInivelCoords = 0;
 
        mesh->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP;
-       mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
        mesh->obstaclePartslip= 0.;
 
-       mesh->channelSizeVertices = 0;
-       mesh->channelVertices = NULL;
+       mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
 
+       /* name of the mesh, mostly for debugging */
        mesh->name = "[unnamed]";
 }
 
index 16d18c4..ba6e638 100644 (file)
@@ -28,9 +28,7 @@ ntlGeometryObjModel::ntlGeometryObjModel( void ) :
        mLoaded( false ),
        mTriangles(), mVertices(), mNormals(),
        mcAniVerts(), mcAniNorms(), 
-       mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.),
-       mvCPSStart(-10000.), mvCPSEnd(10000.),
-       mCPSFilename(""), mCPSWidth(0.1), mCPSTimestep(1.)
+       mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.)
 {
 }
 
@@ -52,6 +50,39 @@ bool ntlGeometryObjModel::getMeshAnimated() {
        return ret;
 }
 
+/*! calculate max extends of (ani) mesh */
+void ntlGeometryObjModel::getExtends(ntlVec3Gfx &sstart, ntlVec3Gfx &send) {
+       bool ini=false;
+       ntlVec3Gfx start(0.),end(0.);
+       for(int s=0; s<=(int)mcAniVerts.accessValues().size(); s++) {
+               vector<ntlVec3f> *sverts;
+               if(mcAniVerts.accessValues().size()>0) {
+                       if(s==(int)mcAniVerts.accessValues().size()) continue;
+                       sverts  = &(mcAniVerts.accessValues()[s].mVerts);
+               } else sverts = &mVertices;
+
+               for(int i=0; i<(int)sverts->size(); i++) {
+
+                       if(!ini) {
+                               start=(*sverts)[i];
+                               end=(*sverts)[i];
+                               //errMsg("getExtends","ini "<<s<<","<<i<<" "<<start<<","<<end);
+                               ini=true;
+                       } else {
+                               for(int j=0; j<3; j++) {
+                                       if(start[j] > (*sverts)[i][j]) { start[j]= (*sverts)[i][j]; }
+                                       if(end[j]   < (*sverts)[i][j]) { end[j]  = (*sverts)[i][j]; }
+                               }
+                               //errMsg("getExtends","check "<<s<<","<<i<<" "<<start<<","<<end<<" "<< (*sverts)[i]);
+                       }
+
+               }
+       }
+       sstart=start;
+       send=end;
+}
+
+
 /*****************************************************************************/
 /* Init attributes etc. of this object */
 /*****************************************************************************/
@@ -90,11 +121,6 @@ void ntlGeometryObjModel::initialize(ntlRenderGlobals *glob)
        mAniTimeScale = mpAttrs->readFloat("ani_timescale", mAniTimeScale,"ntlGeometryObjModel", "mAniTimeScale", false);
        mAniTimeOffset = mpAttrs->readFloat("ani_timeoffset", mAniTimeOffset,"ntlGeometryObjModel", "mAniTimeOffset", false);
 
-       mCPSWidth = mpAttrs->readFloat("cps_width", mCPSWidth,"ntlGeometryObjModel", "mCPSWidth", false);
-       mCPSTimestep = mpAttrs->readFloat("cps_timestep", mCPSTimestep,"ntlGeometryObjModel", "mCPSTimestep", false);
-       mvCPSStart = vec2G( mpAttrs->readVec3d("cps_start", vec2D(mvCPSStart),"ntlGeometryObjModel", "mvCPSStart", false));
-       mvCPSEnd = vec2G( mpAttrs->readVec3d("cps_end", vec2D(mvCPSEnd),"ntlGeometryObjModel", "mvCPSEnd", false));
-
        // continue with standard obj
        if(loadBobjModel(mFilename)==0) mLoaded=1;
        if(!mLoaded) {
@@ -322,7 +348,7 @@ int ntlGeometryObjModel::loadBobjModel(string filename)
                bytesRead += gzread(gzf, &frameTime, sizeof(frameTime) );
                //if(bytesRead!=3*sizeof(float)) {
                if(bytesRead!=sizeof(float)) {
-                       debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' no ani sets. ", 10 );
+                       debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' end of gzfile. ", 10 );
                        if(anitimes.size()>0) {
                                // finally init channels and stop reading file
                                mcAniVerts = AnimChannel<ntlSetVec3f>(aniverts,anitimes);
@@ -428,7 +454,7 @@ ntlGeometryObjModel::getTriangles(double t, vector<ntlTriangle> *triangles,
                (*normals)[startvert+i] = mNormals[i];
        }
 
-       triangles->reserve(triangles->size() + mTriangles.size() );
+       triangles->reserve(triangles->size() + mTriangles.size()/3 );
        for(int i=0; i<(int)mTriangles.size(); i+=3) {
                int trip[3];
                trip[0] = startvert+mTriangles[i+0];
index f911279..95cadc4 100644 (file)
@@ -44,8 +44,8 @@ class ntlGeometryObjModel : public ntlGeometryObject
                /*! init triangle divisions */
                virtual void calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri);
 
-               /*! do ani mesh CPS */
-               void calculateCPS(string filename);
+               /*! calculate max extends of (ani) mesh */
+               void getExtends(ntlVec3Gfx &start, ntlVec3Gfx &end);
 
        private:
 
@@ -71,12 +71,6 @@ class ntlGeometryObjModel : public ntlGeometryObject
                /*! timing mapping & offset for config files */
                double mAniTimeScale, mAniTimeOffset;
 
-               /*! ani mesh cps params */
-               ntlVec3Gfx mvCPSStart, mvCPSEnd;
-               string mCPSFilename;
-               gfxReal mCPSWidth, mCPSTimestep;
-
-
        public:
 
                /* Access methods */
@@ -87,6 +81,9 @@ class ntlGeometryObjModel : public ntlGeometryObject
                inline ntlVec3Gfx getEnd( void ){ return mvEnd; }
                inline void setEnd( const ntlVec3Gfx &set ){ mvEnd = set; }
 
+               inline bool getLoaded( void ){ return mLoaded; }
+               inline void setLoaded( bool set ){ mLoaded = set; }
+
                /*! set data file name */
                inline void setFilename(string set) { mFilename = set; }
 };
index 2321cee..8c0188a 100644 (file)
@@ -331,6 +331,8 @@ public:
                mGeos.push_back( geo ); 
                geo->setObjectId(mGeos.size());
        }
+       /*! Add a geo object to the scene, warning - only needed for hand init */
+       inline void addGeoObject(ntlGeometryObject *geo) { mObjects.push_back( geo ); }
 
        /*! Acces a certain object */
        inline ntlGeometryObject *getObject(int id) { 
index c49c1b1..1d15cec 100644 (file)
@@ -154,7 +154,7 @@ void ParticleTracer::cleanup() {
  *! dump particles if desired 
  *****************************************************************************/
 void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
-       debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts, 10); // DEBUG
+       debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts<<" t"<<simtime, 10); // DEBUG
 
        if(
                        (dumptype==DUMP_FULLGEOMETRY)&&
@@ -305,7 +305,7 @@ void ParticleTracer::checkTrails(double time) {
 /******************************************************************************
  * Get triangles for rendering
  *****************************************************************************/
-void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles, 
+void ParticleTracer::getTriangles(double time, vector<ntlTriangle> *triangles, 
                                                                                                         vector<ntlVec3Gfx> *vertices, 
                                                                                                         vector<ntlVec3Gfx> *normals, int objectId )
 {
@@ -326,7 +326,7 @@ void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles,
        int segments = mPartSegments;
        ntlVec3Gfx scale = ntlVec3Gfx( (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2]));
        ntlVec3Gfx trans = mStart;
-       t = 0.; // doesnt matter
+       time = 0.; // doesnt matter
 
        for(size_t t=0; t<mPrevs.size()+1; t++) {
                vector<ParticleObject> *dparts;
index 5d8a73e..e340301 100644 (file)
@@ -96,6 +96,7 @@ void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) {
        *mpElbeemSettings = *settings;
 
        mGeoInitId = settings->domainId+1;
+       debMsgStd("SimulationObject",DM_MSG,"mGeoInitId="<<mGeoInitId<<", domainId="<<settings->domainId, 8);
 }
 
 /******************************************************************************
@@ -115,7 +116,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        }
 
 
-       this->mGeoInitId = mpAttrs->readInt("geoinitid", this->mGeoInitId,"LbmSolverInterface", "mGeoInitId", false);
+       mGeoInitId = mpAttrs->readInt("geoinitid", mGeoInitId,"LbmSolverInterface", "mGeoInitId", false);
        //mDimension, mSolverType are deprecated
        string mSolverType(""); 
        mSolverType = mpAttrs->readString("solver", mSolverType, "SimulationObject","mSolverType", false ); 
@@ -297,6 +298,7 @@ void SimulationObject::step( void )
                // dont advance for stopped time
                mpLbm->step();
                mTime += mpParam->getTimestep();
+//if(mTime>0.001) { errMsg("DEBUG!!!!!!!!","quit mlsu..."); exit(1); } // PROFILE DEBUG TEST!
        }
        if(mpLbm->getPanic()) mPanic = true;
 
index 10b172d..46462ec 100644 (file)
 #endif
 #endif
 
-#if LBM_INCLUDE_TESTSOLVERS==1
-#include "solver_test.h"
-#endif // LBM_INCLUDE_TESTSOLVERS==1
-
 /*****************************************************************************/
 /*! cell access classes */
 class UniformFsgrCellIdentifier : 
@@ -467,20 +463,6 @@ class LbmFsgrSolver :
 #              endif // FSGR_STRICT_DEBUG==1
 
                bool mUseTestdata;
-#if LBM_INCLUDE_TESTSOLVERS==1
-               // test functions
-               LbmTestdata *mpTest;
-               void initTestdata();
-               void destroyTestdata();
-               void handleTestdata();
-               void set3dHeight(int ,int );
-
-               void initCpdata();
-               void handleCpdata();
-       public:
-               // needed for testdata
-               void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz);
-#endif // LBM_INCLUDE_TESTSOLVERS==1
 
        public: // former LbmModelLBGK  functions
                // relaxation funtions - implemented together with relax macros
index 908355f..226af4e 100644 (file)
@@ -411,7 +411,7 @@ LbmFsgrSolver::~LbmFsgrSolver()
 {
   if(!this->mInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; }
 #if COMPRESSGRIDS==1
-       delete mLevel[mMaxRefine].mprsCells[1];
+       delete [] mLevel[mMaxRefine].mprsCells[1];
        mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
 #endif // COMPRESSGRIDS==1
 
@@ -469,6 +469,12 @@ void LbmFsgrSolver::parseAttrList()
        // FIXME check needed?
        mFVArea   = this->mpAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
 
+       // debugging - skip some time...
+       double starttimeskip = 0.;
+       starttimeskip = this->mpAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
+       mSimulationTime += starttimeskip;
+       if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
+
 #if LBM_INCLUDE_TESTSOLVERS==1
        mUseTestdata = 0;
        mUseTestdata = this->mpAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
@@ -483,7 +489,7 @@ void LbmFsgrSolver::parseAttrList()
        mUseTestdata = 0;
        if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
 #endif // LBM_INCLUDE_TESTSOLVERS!=1
-       if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
+  if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
 
        if(mMaxRefine==0) mInitialCsmago=0.02;
        mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
@@ -1184,6 +1190,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
                                        }
                                        break; 
+                                       // off */
                                        /*
                                case FGI_BNDPART: rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndPartslip;
@@ -1191,7 +1198,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                case FGI_BNDFREE: rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndFreeslip;
                                        break;
-                               // */
+                                       // off */
                                case FGI_BNDNO:   rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndNoslip;
                                        break;
@@ -1876,7 +1883,6 @@ void LbmFsgrSolver::initFreeSurfaces() {
                mLevel[lev].setCurr ^= 1;
        }
        // copy back...?
-
 }
 
 /*****************************************************************************/
index d682445..c4e8d22 100644 (file)
@@ -44,6 +44,7 @@ LbmSolverInterface::LbmSolverInterface() :
        mName("lbm_default") ,
        mpIso( NULL ), mIsoValue(0.499),
        mSilent(false) , 
+       mLbmInitId(1) ,
        mpGiTree( NULL ),
        mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ),
        mRefinementDesired(0),
@@ -284,7 +285,7 @@ void LbmSolverInterface::initGeoTree() {
        if(mpGiTree != NULL) delete mpGiTree;
        char treeFlag = (1<<(this->mLbmInitId+4));
        mpGiTree = new ntlTree( 
-                       15, 8,  // warning - fixed values for depth & maxtriangles here...
+                       15, 8,  // TREEwarning - fixed values for depth & maxtriangles here...
                        scene, treeFlag );
 }
 
@@ -299,6 +300,7 @@ void LbmSolverInterface::freeGeoTree() {
 
 
 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) {
@@ -371,7 +373,8 @@ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int
                                        if(giObjFirstHistSide[i] !=  1) mess=true;
                                }
                                if(mess) {
-                                       errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
+                                       //errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
+                                       globGICPIProblems++;
                                        mGiObjInside[i]++; // believe first hit side...
                                }
                        }
index 2888c7e..8fa3b8d 100644 (file)
@@ -294,6 +294,8 @@ LbmFsgrSolver::mainLoop(int lev)
        int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
        //{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIMcDimension<<" mlsz:"<< mLevel[mMaxRefine].lSizez<<" zmax1:"<<getForZMax1(mMaxRefine) ); } // DEBUG
 #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
+#define LIST_EMPTY(x) mListEmpty.push_back( x );
+#define LIST_FULL(x)  mListFull.push_back( x );
 #endif // PARALLEL==1
 
        // local to loop
@@ -360,22 +362,15 @@ LbmFsgrSolver::mainLoop(int lev)
        } // COMPRT
 
 #if PARALLEL==0
-  const int id = 0, Nthrds = 1;
-#endif // PARALLEL==1
-  const int Nj = mLevel[mMaxRefine].lSizey;
-       int jstart = 0+( id * (Nj / Nthrds) );
-       int jend   = 0+( (id+1) * (Nj / Nthrds) );
-  if( ((Nj/Nthrds) *Nthrds) != Nj) {
-    errMsg("LbmFsgrSolver","Invalid domain size Nj="<<Nj<<" Nthrds="<<Nthrds);
-  }
-       // cutoff obstacle boundary
-       if(jstart<1) jstart = 1;
-       if(jend>mLevel[mMaxRefine].lSizey-1) jend = mLevel[mMaxRefine].lSizey-1;
-
-#if PARALLEL==1
+  //const int id = 0, Nthrds = 1;
+       const int jstart = 0;
+       const int jend   = mLevel[mMaxRefine].lSizey;
+//#endif // PARALLEL==1
+#else // PARALLEL==1
        PARA_INITIALIZE();
-       //errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
+       errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
 #endif // PARALLEL==1
+
   for(int k=kstart;k!=kend;k+=kdir) {
 
        //errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir<<" x*y="<<mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum ); // debug
@@ -459,13 +454,14 @@ LbmFsgrSolver::mainLoop(int lev)
                                changeFlag(lev, i,j,k, TSET(lev), CFInter);
 
                                // same as ifemptied for if below
-                               LbmPoint emptyp; emptyp.flag = 0;
-                               emptyp.x = i; emptyp.y = j; emptyp.z = k;
+                               LbmPoint oemptyp; oemptyp.flag = 0;
+                               oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
 #if PARALLEL==1
-                               calcListEmpty[id].push_back( emptyp );
+                               //calcListEmpty[id].push_back( oemptyp );
 #else // PARALLEL==1
-                               mListEmpty.push_back( emptyp );
+                               //mListEmpty.push_back( oemptyp );
 #endif // PARALLEL==1
+                               LIST_EMPTY(oemptyp);
                                calcCellsEmptied++;
                                continue;
                        }
@@ -601,6 +597,7 @@ LbmFsgrSolver::mainLoop(int lev)
                // for fluid cells - just the f_i difference from streaming to empty cells  ----
                numRecons = 0;
                bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
+       //onlyBndnb = false; // DEBUG test off
 
                FORDF1 { // dfl loop
                        recons[l] = 0;
@@ -1008,6 +1005,7 @@ LbmFsgrSolver::mainLoop(int lev)
 
                // looks much nicer... LISTTRICK
 #if FSGR_LISTTRICK==1
+               if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
                if(newFlag&CFNoBndFluid) { // test NEW TEST
                        if(!iffilled) {
                                // remove cells independent from amount of change...
@@ -1035,10 +1033,11 @@ LbmFsgrSolver::mainLoop(int lev)
                        if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
                        filledp.x = i; filledp.y = j; filledp.z = k;
 #if PARALLEL==1
-                       calcListFull[id].push_back( filledp );
+                       //calcListFull[id].push_back( filledp );
 #else // PARALLEL==1
-                       mListFull.push_back( filledp );
+                       //mListFull.push_back( filledp );
 #endif // PARALLEL==1
+                       LIST_FULL(filledp);
                        //this->mNumFilledCells++; // DEBUG
                        calcCellsFilled++;
                }
@@ -1047,10 +1046,11 @@ LbmFsgrSolver::mainLoop(int lev)
                        if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
                        emptyp.x = i; emptyp.y = j; emptyp.z = k;
 #if PARALLEL==1
-                       calcListEmpty[id].push_back( emptyp );
+                       //calcListEmpty[id].push_back( emptyp );
 #else // PARALLEL==1
-                       mListEmpty.push_back( emptyp );
+                       //mListEmpty.push_back( emptyp );
 #endif // PARALLEL==1
+                       LIST_EMPTY(emptyp);
                        //this->mNumEmptiedCells++; // DEBUG
                        calcCellsEmptied++;
                } 
index eba6742..1a1e2c0 100644 (file)
@@ -17,7 +17,6 @@
 
 
 
-#if LBM_INCLUDE_TESTSOLVERS!=1
 
 // off for non testing
 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
        uy += (grav)[1]; \
        uz += (grav)[2]; \
 
-#else // LBM_INCLUDE_TESTSOLVERS!=1
+#define TEST_IF_CHECK 
 
-// defined in test.h
-
-#endif // LBM_INCLUDE_TESTSOLVERS!=1
 
        
 /******************************************************************************
@@ -1125,7 +1121,7 @@ inline void LbmFsgrSolver::collideArrays(
        for(l=0; l<this->cDfNum; l++) { 
                df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 
        }  
-       if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df);
+       //if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df);
 
        mux = ux;
        muy = uy;
index 8346299..05cfcaa 100644 (file)
@@ -53,33 +53,60 @@ void LbmFsgrSolver::prepareVisualization( void ) {
    for(int j=1;j<mLevel[lev].lSizey-1;j++) 
     for(int i=1;i<mLevel[lev].lSizex-1;i++) {
                        const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
-                       //continue; // OFF DEBUG
-                       if(cflag&(CFBnd|CFEmpty)) {
+                       //if(cflag&(CFBnd|CFEmpty)) {
+                       if(cflag&(CFBnd)) {
                                continue;
 
-                       } else if( (cflag&CFInter) ) {
-                               //} else if( (cflag&CFInter) && (!(cflag&CFNoBndFluid)) && (cflag&CFNoNbFluid) ) {
-                               //} else if( (cflag&CFInter) && (!(cflag&CFNoBndFluid)) ) {
+                       } else if( (cflag&CFEmpty) ) {
+                               //continue; // OFF DEBUG
+                               int noslipbnd = 0;
+                               int intercnt = 0;
+                               CellFlagType nbored;
+                               FORDF1 { 
+                                       const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
+                                       if((nbflag&CFBnd)&&(nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; }
+                                       if(nbflag&CFInter){ intercnt++; }
+                                       nbored |= nbflag;
+                               }
+                               //? val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
+                               if((noslipbnd)&&(intercnt>6)) {
+                                       //if(val<minval) val = minval; 
+                                       //*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
+                                       *this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
+                               } else if((noslipbnd)&&(intercnt>0)) {
+                                       *this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.95;
+                               } else {
+                               }
+                               continue;
+
+                       //} else if( (cflag&CFInter) ) {
+
+                       } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) {
+                               // optimized fluid
+                               val = 1.;
+
+                       } else if( (cflag&(CFInter|CFFluid)) ) {
                                int noslipbnd = 0;
                                FORDF1 { 
                                        const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
                                        if((nbflag&CFBnd)&&(nbflag&CFBnd)&&(CFBndNoslip)){ noslipbnd=1; l=100; } 
                                }
+                               // no empty nb interface cells are treated as full
+                               if(cflag&(CFNoNbEmpty|CFFluid)) {
+                                       val=1.0;
+                               }
                                val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
+                               
                                if(noslipbnd) {
                                        //errMsg("NEWVAL", PRINT_IJK<<" val"<<val <<" f"<< convertCellFlagType2String(cflag)<<" "<<noslipbnd); //" nbfl"<<convertCellFlagType2String(nbored) );
                                        if(val<minval) val = minval; 
-                                       *this->mpIso->lbmGetData( i   , j    ,ZKOFF  ) += minval-( val * mIsoWeight[13] ); 
+                                       *this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
                                }
                                // */
-                               
-                               // no empty nb interface cells are treated as full
-                               if(cflag&CFNoNbEmpty) {
-                                       val=1.0;
-                               }
 
-                       } else { // fluid?
-                               val = 1.0; 
+                       } else { // unused?
+                               continue;
+                               // old fluid val = 1.0; 
                        } // */
 
                        *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
@@ -885,10 +912,11 @@ void LbmFsgrSolver::advanceParticles() {
 }
 
 void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
+       int workSet = mLevel[mMaxRefine].setCurr;
+       std::ostringstream name;
+
        // debug - raw dump of ffrac values
        if(getenv("ELBEEM_RAWDEBUGDUMP")) {
-               int workSet = mLevel[mMaxRefine].setCurr;
-               std::ostringstream name;
                //name <<"fill_" << this->mStepCnt <<".dump";
                name << outfilename<< frameNrStr <<".dump";
                FILE *file = fopen(name.str().c_str(),"w");
@@ -898,9 +926,12 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt
                                for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
                                        for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
                                                float val = 0.;
-                                               if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+                                               if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
+                                                       val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+                                                       if(val<0.) val=0.;
+                                                       if(val>1.) val=1.;
+                                               }
                                                if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
-                                               //fwrite( &val, sizeof(val), 1, file); // binary
                                                fprintf(file, "%f ",val); // text
                                                //errMsg("W", PRINT_IJK<<" val:"<<val);
                                        }
@@ -912,6 +943,27 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt
 
                } // file
        } // */
+       if(getenv("ELBEEM_BINDEBUGDUMP")) {
+               name << outfilename<< frameNrStr <<".bdump";
+               FILE *file = fopen(name.str().c_str(),"w");
+               if(file) {
+                       for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
+                               for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
+                                       for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
+                                               float val = 0.;
+                                               if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
+                                                       val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+                                                       if(val<0.) val=0.;
+                                                       if(val>1.) val=1.;
+                                               }
+                                               if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
+                                               fwrite( &val, sizeof(val), 1, file); // binary
+                                       }
+                               }
+                       }
+                       fclose(file);
+               } // file
+       }
 
        dumptype = 0; frameNr = 0; // get rid of warning
 }