yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00710
[svn] r1539 - in trunk: core gui/py pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/StandAloneEngine scripts
Author: eudoxos
Date: 2008-10-09 14:22:33 +0200 (Thu, 09 Oct 2008)
New Revision: 1539
Modified:
trunk/core/Omega.cpp
trunk/core/Omega.hpp
trunk/gui/py/utils.py
trunk/gui/py/yadeControl.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/scripts/chain-distant-interactions.py
Log:
1. Add Omega().reload()
2. Create map<string,string> Omega::memSavedSimulations, where simulations are saved if they begin with special sequence ":memory:" and can be loaded again, within lifetime of Omega.
3. Add Omega().tmpLoad() and Omega().tmpSave() that use the above mentioned thing.
4. add utils.readParamsFromTable that reads any python variables from text file (useful for automatic parametric studies)
5. Add less expensive routine for sontact point relocation in SpheresContactGeometry
Modified: trunk/core/Omega.cpp
===================================================================
--- trunk/core/Omega.cpp 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/core/Omega.cpp 2008-10-09 12:22:33 UTC (rev 1539)
@@ -68,6 +68,7 @@
void Omega::init(){
simulationFileName="";
resetRootBody();
+ timeInit();
}
void Omega::timeInit(){
@@ -215,7 +216,6 @@
void Omega::loadSimulationFromStream(std::istream& stream){
LOG_DEBUG("Loading simulation from stream.");
resetRootBody();
- timeInit();
IOFormatManager::loadFromStream("XMLFormatManager",stream,"rootBody",rootBody);
}
void Omega::saveSimulationToStream(std::ostream& stream){
@@ -223,10 +223,10 @@
IOFormatManager::saveToStream("XMLFormatManager",stream,"rootBody",rootBody);
}
-void Omega::loadSimulation()
-{
+void Omega::loadSimulation(){
+
if(Omega::instance().getSimulationFileName().size()==0) throw yadeBadFile("Simulation filename to load has zero length");
- if(!filesystem::exists(simulationFileName)) throw yadeBadFile("Simulation file to load doesn't exist");
+ if(!filesystem::exists(simulationFileName) && !algorithm::starts_with(simulationFileName,":memory")) throw yadeBadFile("Simulation file to load doesn't exist");
// FIXME: should stop running simulation!!
LOG_INFO("Loading file " + simulationFileName);
@@ -244,6 +244,12 @@
resetRootBody();
IOFormatManager::loadFromFile("BINFormatManager",simulationFileName,"rootBody",rootBody);
}
+ else if(algorithm::starts_with(simulationFileName,":memory:")){
+ if(memSavedSimulations.count(simulationFileName)==0) throw yadeBadFile(("Cannot load nonexistent memory-saved simulation "+simulationFileName).c_str());
+ resetRootBody();
+ istringstream iss(memSavedSimulations[simulationFileName]);
+ IOFormatManager::loadFromStream("XMLFormatManager",iss,"rootBody",rootBody);
+ }
else throw (yadeBadFile("Extension of file not recognized."));
}
@@ -269,6 +275,12 @@
FormatChecker::format=FormatChecker::BIN;
IOFormatManager::saveToFile("BINFormatManager",name,"rootBody",rootBody);
}
+ else if(algorithm::starts_with(name,":memory:")){
+ if(memSavedSimulations.count(simulationFileName)>0) LOG_INFO("Overwriting in-memory saved simulation "<<name);
+ ostringstream oss;
+ IOFormatManager::saveToStream("XMLFormatManager",oss,"rootBody",rootBody);
+ memSavedSimulations[name]=oss.str();
+ }
else {
throw(yadeBadFile(("Filename extension not recognized in `"+name+"'").c_str()));
}
Modified: trunk/core/Omega.hpp
===================================================================
--- trunk/core/Omega.hpp 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/core/Omega.hpp 2008-10-09 12:22:33 UTC (rev 1539)
@@ -87,6 +87,8 @@
string simulationFileName;
void buildDynlibDatabase(const vector<string>& dynlibsList); // FIXME - maybe in ClassFactory ?
+ map<string,string> memSavedSimulations;
+
public :
shared_ptr<Preferences> preferences;
string yadeConfigPath; // FIXME - must be private and more clean
Modified: trunk/gui/py/utils.py
===================================================================
--- trunk/gui/py/utils.py 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/gui/py/utils.py 2008-10-09 12:22:33 UTC (rev 1539)
@@ -17,6 +17,8 @@
from yade._utils import *
+def typedEngine(name): return [e for e in Omega().engines if e.name==name][0]
+
def sphere(center,radius,density=1,young=30e9,poisson=.3,frictionAngle=0.5236,dynamic=True,wire=False,color=[1,1,1],physParamsClass='BodyMacroParameters'):
"""Create default sphere, with given parameters. Physical properties such as mass and inertia are calculated automatically."""
s=Body()
@@ -216,3 +218,50 @@
mainloop.run()
pipeline.set_state(gst.STATE_NULL); pipeline.get_state()
+def readParamsFromTable(tableFile=None,tableLine=None,noTableOk=False,**kw):
+ """
+ Read parameters from a file and assign them to __builtin__ variables.
+
+ tableFile is a text file (with one value per blank-separated columns)
+ tableLine is number of line where to get the values from
+
+ The format of the file is as follows (no comments, empty lines etc allowed…)
+
+ name1 name2 … # 0th line
+ val1 val2 … # 1st line
+ val2 val2 … # 2nd line
+ …
+
+ The name `description' is special and is assigned to Omega().tags['description']
+
+ assigns Omega().tags['params']="name1=val1,name2=val2,…"
+
+ assigns Omega().tags['defaultParams']="unassignedName1=defaultValue1,…"
+
+ return value is the number of assigned parameters.
+ """
+ o=Omega()
+ tagsParams=[]
+ import os, __builtin__
+ if not tableFile and not os.environ.has_key('PARAM_TABLE'):
+ if not noTableOk: raise EnvironmentError("PARAM_TABLE is not defined in the environment")
+ else:
+ if not tableFile: tableFile=os.environ['PARAM_TABLE']
+ if not tableLine: tableLine=int(os.environ['PARAM_LINE'])
+ ll=open(tableFile).readlines(); names=ll[0].split(); values=ll[tableLine].split()
+ assert(len(names)==len(values))
+ for i in range(len(names)):
+ if names[i]=='description': o.tags['description']=values[i]
+ else:
+ if names[i] not in kw.keys(): raise NameError("Parameter `%s' has no default value assigned"%names[i])
+ kw.pop(names[i])
+ eq="%s=%s"%(names[i],values[i])
+ exec('__builtin__.%s=%s'%(names[i],values[i])); tagsParams+=['%s=%s'%(names[i],values[i])]
+ defaults=[]
+ for k in kw.keys():
+ exec("__builtin__.%s=%s"%(k,kw[k]))
+ defaults+=["%s=%s"%(k,kw[k])]
+ o.tags['defaultParams']=",".join(defaults)
+ o.tags['params']=",".join(tagsParams)
+ return len(tagsParams)
+
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/gui/py/yadeControl.cpp 2008-10-09 12:22:33 UTC (rev 1539)
@@ -424,7 +424,11 @@
OMEGA.createSimulationLoop();
LOG_DEBUG("LOAD!");
}
+ void reload(){ load(OMEGA.getSimulationFileName());}
+ void saveTmp(string mark){ save(":memory:"+mark);}
+ void loadTmp(string mark){ load(":memory:"+mark);}
+
void reset(){Py_BEGIN_ALLOW_THREADS; OMEGA.reset(); Py_END_ALLOW_THREADS; }
void save(std::string fileName){
@@ -520,7 +524,10 @@
.add_property("dt",&pyOmega::dt_get,&pyOmega::dt_set)
.add_property("usesTimeStepper",&pyOmega::usesTimeStepper_get,&pyOmega::usesTimeStepper_set)
.def("load",&pyOmega::load)
+ .def("reload",&pyOmega::reload)
.def("save",&pyOmega::save)
+ .def("loadTmp",&pyOmega::loadTmp)
+ .def("saveTmp",&pyOmega::saveTmp)
.def("saveSpheres",&pyOmega::saveSpheres)
.def("run",&pyOmega::run,omega_run_overloads())
.def("pause",&pyOmega::pause)
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp 2008-10-09 12:22:33 UTC (rev 1539)
@@ -25,8 +25,15 @@
*/
void SpheresContactGeometry::relocateContactPoints(){
assert(hasShear);
- Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
- Vector3r midPt=.5*(p1+p2);
+ relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2());
+}
+
+/* Like SpheresContactGeometry::relocateContactPoints(), but use already computed tangent plane points.
+ *
+ *
+ */
+void SpheresContactGeometry::relocateContactPoints(const Vector3r& p1, const Vector3r& p2){
+ Vector3r midPt=(d1/(d1+d2))*(p1+p2); // proportionally to radii, so that angle would be the same
if((p1.SquaredLength()>4*d1 || p2.SquaredLength()>4*d2) && midPt.SquaredLength()>.5*min(d1,d2)){
//cerr<<"RELOCATION with displacementT="<<displacementT(); // should be the same before and after relocation
setTgPlanePts(p1-midPt,p2-midPt);
@@ -34,6 +41,7 @@
}
}
+
/*! Perform slip of the projected contact points so that their distance becomes equal (or remains smaller) than the given one.
* The slipped distance is returned.
*/
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2008-10-09 12:22:33 UTC (rev 1539)
@@ -22,7 +22,9 @@
*
* TODO: add torsion code, that will (if a flag is on) compute torsion angle on the contact axis.
*
+ * TODO: add bending code, that will compute bending angle of the contact axis
*
+ *
*/
class SpheresContactGeometry: public InteractionGeometry{
public:
@@ -58,7 +60,17 @@
Real displacementN(){assert(hasShear); return (pos2-pos1).Length()-d0;}
Real epsN(){return displacementN()*(1./d0);}
- Vector3r displacementT(){ assert(hasShear); return contPtInTgPlane2()-contPtInTgPlane1();}
+ Vector3r displacementT(){ assert(hasShear);
+ // enabling automatic relocation decreases overall simulation speed by about 3%
+ // perhaps: bool largeStrains ... ?
+ #if 0
+ Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
+ relocateContactPoints(p1,p2);
+ return p2-p1; // shear before relocation, but that should be OK
+ #else
+ return contPtInTgPlane2()-contPtInTgPlane1();
+ #endif
+ }
Vector3r epsT(){return displacementT()*(1./d0);}
Real slipToDisplacementTMax(Real displacementTMax);
@@ -66,6 +78,7 @@
Real slipToEpsTMax(Real epsTMax){ return (1/d0)*slipToDisplacementTMax(epsTMax*d0);}
void relocateContactPoints();
+ void relocateContactPoints(const Vector3r& tgPlanePt1, const Vector3r& tgPlanePt2);
SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO){createIndex();}
virtual ~SpheresContactGeometry();
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2008-10-09 12:22:33 UTC (rev 1539)
@@ -36,6 +36,7 @@
Real Fn=contPhys->kn*contGeom->displacementN(); // scalar normal force; displacementN()>=0 ≡ elongation of the contact
if(!isCohesive && contGeom->displacementN()>0){ cerr<<"deleting"<<endl; /* delete the interaction */ i->isReal=false; continue;}
contPhys->normalForce=Fn*contGeom->normal;
+ //contGeom->relocateContactPoints();
contGeom->slipToDisplacementTMax(max(0.,(-Fn*contPhys->tangensOfFrictionAngle)/contPhys->ks)); // limit shear displacement -- Coulomb criterion
contPhys->shearForce=contPhys->ks*contGeom->displacementT();
Vector3r force=contPhys->shearForce+contPhys->normalForce;
Modified: trunk/scripts/chain-distant-interactions.py
===================================================================
--- trunk/scripts/chain-distant-interactions.py 2008-10-07 07:02:48 UTC (rev 1538)
+++ trunk/scripts/chain-distant-interactions.py 2008-10-09 12:22:33 UTC (rev 1539)
@@ -15,7 +15,7 @@
MetaEngine('InteractionPhysicsMetaEngine',[EngineUnit('SimpleElasticRelationships')]),
StandAloneEngine('ElasticContactLaw2',{'isCohesive':True}),
DeusExMachina('GravityEngine',{'gravity':[0,0,-1e5]}),
- DeusExMachina('NewtonsDampedLaw',{'damping':0.2})
+ DeusExMachina('NewtonsDampedLaw',{'damping':0.1})
]
from yade import utils
for n in range(20):
@@ -25,4 +25,7 @@
o.dt=.04*utils.PWaveTimeStep()
-o.save('/tmp/a.xml.bz2')
+#o.save('/tmp/a.xml.bz2')
+#o.reload()
+#o.run(50000,True)
+#print o.iter/o.realtime
