yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00824
[svn] r1596 - in trunk: extra/clump gui/py gui/qt3 pkg/common/Engine/DeusExMachina
Author: eudoxos
Date: 2008-12-09 12:43:29 +0100 (Tue, 09 Dec 2008)
New Revision: 1596
Modified:
trunk/extra/clump/Shop.cpp
trunk/extra/clump/Shop.hpp
trunk/gui/py/_eudoxos.cpp
trunk/gui/py/_utils.cpp
trunk/gui/py/eudoxos.py
trunk/gui/py/utils.py
trunk/gui/qt3/GLViewer.cpp
trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp
Log:
1. Fix warnings in GLViewer
2. Add (35.4%)-style counter if stopAtIteration is set to the OSD in gl view
3. Add utils.downCast to downcast object to its derived class and copying all attributes.
4. Add utils.SpherePWaveTimeStep to compute critical dt from a few parameters only
5. Add routine to compute overall kinetic energy to Shop (indirectly to utils as well)
6. Add routine to sum torque (momentum) with respect to an axis from forces and torques on given set of bodies.
Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/extra/clump/Shop.cpp 2008-12-09 11:43:29 UTC (rev 1596)
@@ -135,10 +135,23 @@
return (useMaxForce?maxF:meanF)/maxContactF;
}
+Real Shop::kineticEnergy(MetaBody* _rb){
+ MetaBody* rb=_rb ? _rb : Omega::instance().getRootBody().get();
+ Real ret=0.;
+ FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ if(!b->isDynamic) continue;
+ shared_ptr<RigidBodyParameters> rbp=YADE_PTR_CAST<RigidBodyParameters>(b->physicalParameters); assert(pp);
+ Matrix3r inertiaMatrix(rbp->inertia[0],rbp->inertia[1],rbp->inertia[2]);
+ // ½(mv²+ωIω)
+ ret+=.5*(rbp->mass*rbp->velocity.SquaredLength()+rbp->angularVelocity.Dot(diagMult(rbp->inertia,rbp->angularVelocity)));
+ }
+ return ret;
+}
+
template <typename valType> valType Shop::getDefault(const string& key) {
ensureInit();
try{return boost::any_cast<valType>(defaults[key]);}
Modified: trunk/extra/clump/Shop.hpp
===================================================================
--- trunk/extra/clump/Shop.hpp 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/extra/clump/Shop.hpp 2008-12-09 11:43:29 UTC (rev 1596)
@@ -101,6 +101,7 @@
//! Get unbalanced force of the whole simulation
static Real unbalancedForce(bool useMaxForce=false, MetaBody* _rb=NULL);
+ static Real kineticEnergy(MetaBody* _rb=NULL);
//! create transientInteraction between 2 bodies, using existing MetaEngine in Omega
static shared_ptr<Interaction> createExplicitInteraction(body_id_t id1, body_id_t id2);
Modified: trunk/gui/py/_eudoxos.cpp
===================================================================
--- trunk/gui/py/_eudoxos.cpp 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/gui/py/_eudoxos.cpp 2008-12-09 11:43:29 UTC (rev 1596)
@@ -1,6 +1,6 @@
#include<yade/extra/Brefcom.hpp>
#include<boost/python.hpp>
-using namespace boost;
+using namespace boost::python;
using namespace std;
# if 0
Real elasticEnergyDensityInAABB(python::tuple AABB){
@@ -27,3 +27,34 @@
}
#endif
+// copied from _utils.cpp
+Vector3r tuple2vec(const python::tuple& t){return Vector3r(extract<double>(t[0])(),extract<double>(t[1])(),extract<double>(t[2])());}
+
+/*! Set velocity of all dynamic particles so that if their motion were unconstrained,
+ * axis given by axisPoint and axisDirection would be reached in timeToAxis
+ * (or, point at distance subtractDist from axis would be reached).
+ *
+ * The code is analogous to AxialGravityEngine and is intended to give initial motion
+ * to particles subject to axial compaction to speed up the process. */
+void velocityTowardsAxis(python::tuple _axisPoint, python::tuple _axisDirection, Real timeToAxis, Real subtractDist=0.){
+ Vector3r axisPoint=tuple2vec(_axisPoint), axisDirection=tuple2vec(_axisDirection);
+ FOREACH(const shared_ptr<Body>&b, *(Omega::instance().getRootBody()->bodies)){
+ if(!b->isDynamic) continue;
+ ParticleParameters* pp=YADE_CAST<ParticleParameters*>(b->physicalParameters.get());
+ const Vector3r& x0=pp->se3.position;
+ const Vector3r& x1=axisPoint;
+ const Vector3r x2=axisPoint+axisDirection;
+ Vector3r closestAxisPoint=(x2-x1) * /* t */ (-(x1-x0).Dot(x2-x1))/((x2-x1).SquaredLength());
+ Vector3r toAxis=closestAxisPoint-x0;
+ if(subtractDist>0) toAxis*=(toAxis.Length()-subtractDist)/toAxis.Length();
+ pp->velocity=toAxis/timeToAxis;
+ }
+}
+BOOST_PYTHON_FUNCTION_OVERLOADS(velocityTowardsAxis_overloads,velocityTowardsAxis,3,4);
+
+
+
+BOOST_PYTHON_MODULE(_eudoxos){
+ def("velocityTowardsAxis",velocityTowardsAxis,velocityTowardsAxis_overloads(args("axisPoint","axisDirection","timeToAxis","subtractDist")));
+}
+
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/gui/py/_utils.cpp 2008-12-09 11:43:29 UTC (rev 1596)
@@ -171,7 +171,7 @@
* is position relative to axisPt; moment from moment is m; such moment per body is
* projected onto axis.
*/
-Real sumBexMoments(int mask, python::tuple _axis, python::tuple _axisPt){
+Real sumBexTorques(int mask, python::tuple _axis, python::tuple _axisPt){
Shop::Bex::initCache();
shared_ptr<MetaBody> rb=Omega::instance().getRootBody();
Real ret;
@@ -264,6 +264,7 @@
void Shop__createExplicitInteraction(body_id_t id1, body_id_t id2){ (void) Shop::createExplicitInteraction(id1,id2);}
BOOST_PYTHON_FUNCTION_OVERLOADS(unbalancedForce_overloads,Shop::unbalancedForce,0,1);
+Real Shop__kineticEnergy(){return Shop::kineticEnergy();}
BOOST_PYTHON_MODULE(_utils){
// http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
@@ -279,8 +280,9 @@
def("elasticEnergy",elasticEnergyInAABB);
def("inscribedCircleCenter",inscribedCircleCenter);
def("unbalancedForce",&Shop::unbalancedForce,unbalancedForce_overloads(args("useMaxForce")));
+ def("kineticEnergy",Shop__kineticEnergy);
def("sumBexForces",sumBexForces);
- def("sumBexMoments",sumBexMoments);
+ def("sumBexTorques",sumBexTorques);
def("createInteraction",Shop__createExplicitInteraction);
def("spiralProject",spiralProject,spiralProject_overloads(args("axis","theta0")));
def("pointInsidePolygon",pointInsidePolygon);
Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/gui/py/eudoxos.py 2008-12-09 11:43:29 UTC (rev 1596)
@@ -5,6 +5,7 @@
#
from yade.wrapper import *
from math import *
+from yade._eudoxos import * ## c++ implementations
def estimateStress(strain,cutoff=0.):
"""Use summed stored energy in contacts to compute macroscopic stress over the same volume, provided known strain."""
Modified: trunk/gui/py/utils.py
===================================================================
--- trunk/gui/py/utils.py 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/gui/py/utils.py 2008-12-09 11:43:29 UTC (rev 1596)
@@ -16,7 +16,7 @@
# c++ implementations for performance reasons
from yade._utils import *
-def saveVars(mark='',**kw):
+def saveVars(mark='',loadNow=False,**kw):
"""Save passed variables into the simulation so that it can be recovered when the simulation is loaded again.
For example, variables a=5, b=66 and c=7.5e-4 are defined. To save those, use
@@ -26,12 +26,15 @@
those variables will be save in the .xml file, when the simulation itself is saved. To recover those variables once
the .xml is loaded again, use
- utils.loadVars()
+ utils.loadVars(mark)
and they will be defined in the __builtin__ namespace (i.e. available from anywhere in the python code).
+
+ If loadParam==True, variables will be loaded immediately after saving. That effectively makes **kw available in builtin namespace.
"""
import cPickle
Omega().tags['pickledPythonVariablesDictionary'+mark]=cPickle.dumps(kw)
+ if loadNow: loadVars(mark)
def loadVars(mark=''):
import cPickle
@@ -40,10 +43,21 @@
for k in d: __builtin__.__dict__[k]=d[k]
+def SpherePWaveTimeStep(radius,density,young):
+ """Compute P-wave critical timestep for a single sphere.
+ If you want to compute minimum critical timestep for all spheres in the simulation, use utils.PWaveTimeStep() instead."""
+ from math import sqrt
+ return radius/sqrt(young/density)
+
def randomColor(): return [random.random(),random.random(),random.random()]
def typedEngine(name): return [e for e in Omega().engines if e.name==name][0]
+def downCast(obj,newClassName):
+ """Cast given object to class deriving from the same yade root class and copy all parameters from given object.
+ Obj should be up in the inheritance tree, otherwise some attributes may not be defined in the new class."""
+ return obj.__class__(newClassName,dict([ (key,obj[key]) for key in obj.keys() ]))
+
def sphere(center,radius,density=1,young=30e9,poisson=.3,frictionAngle=0.5236,dynamic=True,wire=False,color=None,physParamsClass='BodyMacroParameters'):
"""Create default sphere, with given parameters. Physical properties such as mass and inertia are calculated automatically."""
s=Body()
Modified: trunk/gui/qt3/GLViewer.cpp
===================================================================
--- trunk/gui/qt3/GLViewer.cpp 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/gui/qt3/GLViewer.cpp 2008-12-09 11:43:29 UTC (rev 1596)
@@ -149,7 +149,7 @@
string GLViewer::getState(){
QString origStateFileName=stateFileName();
- char tmpnam_str [L_tmpnam]; tmpnam(tmpnam_str);
+ char tmpnam_str [L_tmpnam]; char* result=tmpnam(tmpnam_str); if(!result) throw runtime_error("Failure getting temporary filename.");
setStateFileName(tmpnam_str); saveStateToFile(); setStateFileName(origStateFileName);
LOG_DEBUG("State saved to temp file "<<tmpnam_str);
// read tmp file contents and return it as string
@@ -160,7 +160,7 @@
}
void GLViewer::setState(string state){
- char tmpnam_str [L_tmpnam]; tmpnam(tmpnam_str);
+ char tmpnam_str [L_tmpnam]; char* result=tmpnam(tmpnam_str); if(!result) throw runtime_error("Failure getting temporary filename.");
std::ofstream out(tmpnam_str);
if(!out.good()){ LOG_ERROR("Error opening temp file `"<<tmpnam_str<<"', loading aborted."); return; }
out<<state; out.close();
@@ -492,7 +492,8 @@
drawArrow(wholeDiameter/6);
glPopMatrix();
}
-
+
+ MetaBody* rb=Omega::instance().getRootBody().get();
#define _W3 setw(3)<<setfill('0')
#define _W2 setw(2)<<setfill('0')
if(timeDispMask!=0){
@@ -519,6 +520,7 @@
if(timeDispMask & GLViewer::TIME_ITER){
ostringstream oss;
oss<<"#"<<Omega::instance().getCurrentIteration();
+ if(rb->stopAtIteration>rb->currentIteration) oss<<" ("<<setiosflags(ios::fixed)<<setw(3)<<setprecision(1)<<setfill('0')<<(100.*rb->currentIteration)/rb->stopAtIteration<<"%)";
QGLViewer::drawText(x,y,oss.str());
y-=lineHt;
}
@@ -530,7 +532,7 @@
ostringstream oss;
time_duration t=Omega::instance().getComputationDuration();
unsigned d=t.hours()/24,h=t.hours()%24,m=t.minutes(),s=t.seconds();
- oss<<"wall ";
+ oss<<"clock ";
if(d>0) oss<<d<<"days "<<_W2<<h<<":"<<_W2<<m<<":"<<_W2<<s;
else if(h>0) oss<<_W2<<h<<":"<<_W2<<m<<":"<<_W2<<s;
else oss<<_W2<<m<<":"<<_W2<<s;
Modified: trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp 2008-12-08 09:27:08 UTC (rev 1595)
+++ trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp 2008-12-09 11:43:29 UTC (rev 1596)
@@ -40,9 +40,7 @@
}
void AxialGravityEngine::applyCondition(MetaBody* rootBody){
- BodyContainer::iterator Iend=rootBody->bodies->end();
- for(BodyContainer::iterator I=rootBody->bodies->begin(); I!=Iend; ++I){
- const shared_ptr<Body>& b(*I);
+ FOREACH(const shared_ptr<Body>&b, *rootBody->bodies){
if(b->isClumpMember()) continue;
Real myMass=YADE_PTR_CAST<ParticleParameters>(b->physicalParameters)->mass;
/* http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html */
