yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05682
[Branch ~yade-dev/yade/trunk] Rev 2439: 1. Resurrect SnapshotEngine
------------------------------------------------------------
revno: 2439
committer: Chiara Modenese <chia@engs-018373>
branch nick: trunk
timestamp: Sat 2010-09-18 15:10:30 +0100
message:
1. Resurrect SnapshotEngine
2. Add Gl1_NormPhys for displaying interaction network (http://beta.arcig.cz/~eudoxos/temp/Gl1_NormPhys.avi)
3. Make VTKRecroder work nicely with periodic boundary conditions (repeating interaction on both sides if it crosses the cell boundary using a few hacks)
4. Add Shop::normalShearStressTensors for computing (once more) stress tensor, but normal and shear contributions separately
5. Allow specify particle counts for SpherePack.particleSD
6. utils.makeVideo now uses mencoder (and works, unlike the older gstreamer-based thing, which would not work anymore)
7. Add Matrix3.diagonal()
8. Integrate regression tests suite in the main program (yade --test)
9. Set plot.scientific==True only if the version of pylab installed supports it (thanks to Anton for reporting)
added:
pkg/common/RenderingEngine/Gl1_NormPhys.cpp
pkg/common/RenderingEngine/Gl1_NormPhys.hpp
modified:
core/Cell.hpp
core/main/main.py.in
core/main/yade-batch.in
debian/rules
doc/references.bib
gui/SConscript
gui/qt4/GLViewer.cpp
gui/qt4/GLViewer.hpp
gui/qt4/_GLViewer.cpp
pkg/common/RenderingEngine/OpenGLRenderer.cpp
pkg/dem/DataClass/SpherePack.cpp
pkg/dem/DataClass/SpherePack.hpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp
pkg/dem/meta/Shop.cpp
pkg/dem/meta/Shop.hpp
py/_utils.cpp
py/mathWrap/miniEigen.cpp
py/pack/_packSpheres.cpp
py/plot.py
py/tests/omega.py
py/utils.py
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'core/Cell.hpp'
--- core/Cell.hpp 2010-09-02 09:11:04 +0000
+++ core/Cell.hpp 2010-09-18 14:10:30 +0000
@@ -61,43 +61,44 @@
//! "integrate" velGrad, update cached values used by public getter
void integrateAndUpdate(Real dt);
/*! Return point inside periodic cell, even if shear is applied */
- Vector3r wrapShearedPt(const Vector3r& pt){ return shearPt(wrapPt(unshearPt(pt))); }
+ Vector3r wrapShearedPt(const Vector3r& pt) const { return shearPt(wrapPt(unshearPt(pt))); }
/*! Return point inside periodic cell, even if shear is applied; store cell coordinates in period. */
- Vector3r wrapShearedPt(const Vector3r& pt, Vector3i& period){ return shearPt(wrapPt(unshearPt(pt),period)); }
+ Vector3r wrapShearedPt(const Vector3r& pt, Vector3i& period) const { return shearPt(wrapPt(unshearPt(pt),period)); }
/*! Apply inverse shear on point; to put it inside (unsheared) periodic cell, apply wrapPt on the returned value. */
- Vector3r unshearPt(const Vector3r& pt){ return _unshearTrsf*pt; }
+ Vector3r unshearPt(const Vector3r& pt) const { return _unshearTrsf*pt; }
//! Apply shear on point.
- Vector3r shearPt(const Vector3r& pt){ return _shearTrsf*pt; }
+ Vector3r shearPt(const Vector3r& pt) const { return _shearTrsf*pt; }
/*! Wrap point to inside the periodic cell; don't compute number of periods wrapped */
- Vector3r wrapPt(const Vector3r& pt)const{
+ Vector3r wrapPt(const Vector3r& pt) const {
Vector3r ret; for(int i=0;i<3;i++) ret[i]=wrapNum(pt[i],_size[i]); return ret;
}
/*! Wrap point to inside the periodic cell; period will contain by how many cells it was wrapped. */
- Vector3r wrapPt(const Vector3r& pt, Vector3i& period)const{
+ Vector3r wrapPt(const Vector3r& pt, Vector3i& period) const {
Vector3r ret; for(int i=0; i<3; i++){ ret[i]=wrapNum(pt[i],_size[i],period[i]); } return ret;
}
/*! Wrap number to interval 0…sz */
- static Real wrapNum(const Real& x, const Real& sz){
+ static Real wrapNum(const Real& x, const Real& sz) {
Real norm=x/sz; return (norm-floor(norm))*sz;
}
/*! Wrap number to interval 0…sz; store how many intervals were wrapped in period */
- static Real wrapNum(const Real& x, const Real& sz, int& period){
+ static Real wrapNum(const Real& x, const Real& sz, int& period) {
Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz;
}
- Vector3r getRefSize(){ return refSize; }
+ Vector3r getRefSize() const { return refSize; }
void setRefSize(const Vector3r& s){ refSize=s; integrateAndUpdate(0); }
- Matrix3r getTrsf(){ return trsf; }
+ Matrix3r getTrsf() const { return trsf; }
void setTrsf(const Matrix3r& m){ trsf=m; integrateAndUpdate(0); }
// return current cell volume
- Real getVolume(){ return Hsize.determinant(); }
+ Real getVolume() const { return Hsize.determinant(); }
void postLoad(Cell&){ integrateAndUpdate(0); }
- Vector3r wrapShearedPt_py(const Vector3r& pt){ return wrapShearedPt(pt);}
- Vector3r wrapPt_py(const Vector3r& pt){ return wrapPt(pt);}
+ Vector3r wrapShearedPt_py(const Vector3r& pt) const { return wrapShearedPt(pt);}
+ Vector3r wrapPt_py(const Vector3r& pt) const { return wrapPt(pt);}
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
+ /* TODO: remove the overrides, add Attr::pyCallPostLoad */
((Vector3r,refSize,Vector3r(1,1,1),,"[will be overridden below]"))
((Matrix3r,trsf,Matrix3r::Identity(),,"[will be overridden below]"))
((Matrix3r,velGrad,Matrix3r::Zero(),,"Velocity gradient of the transformation; used in NewtonIntegrator."))
=== modified file 'core/main/main.py.in'
--- core/main/main.py.in 2010-09-02 09:11:04 +0000
+++ core/main/main.py.in 2010-09-18 14:10:30 +0000
@@ -29,6 +29,7 @@
par.add_option('-n',help="Run without graphical interface (equivalent to unsetting the DISPLAY environment variable)",dest='nogui',action='store_true')
par.add_option('--generate-manpage',help="Generate man page documenting this program and exit",dest='manpage',metavar='FILE')
par.add_option('--rebuild',help="Re-run build in the source directory, then run the updated yade with the same command line except --rebuild. The build profile for this build (${profile}) and its stored parameters will be used.",dest='rebuild',action='store_true')
+par.add_option('--test',help="Run regression test suite and exit; the exists status is 0 if all tests pass, 1 if a test fails and 2 for an unspecified exception.",dest="test",action='store_true')
par.disable_interspersed_args()
opts,args=par.parse_args()
@@ -46,6 +47,26 @@
print 'Running yade using',' '.join(argv)
sys.exit(subprocess.call(argv))
+# run regression test suite and exit
+if opts.test:
+ import yade.tests
+ try:
+ result=yade.tests.testAll()
+ except:
+ print 20*'*'+' UNEXPECTED EXCEPTION WHILE RUNNING TESTS '+20*'*'
+ print 20*'*'+' '+str(sys.exc_info()[0])
+ print 20*'*'+" Please report bug at http://bugs.launchpad.net/yade providing the following traceback:"
+ import traceback; traceback.print_exc()
+ print 20*'*'+' Thank you '+20*'*'
+ sys.exit(2)
+ if result.wasSuccessful():
+ print "*** ALL TESTS PASSED ***"
+ sys.exit(0)
+ else:
+ print 20*'*'+' SOME TESTS FAILED '+20*'*'
+ sys.exit(1)
+
+
# c++ boot code checks for YADE_DEBUG at some places; debug verbosity is equivalent
# do this early, to have debug messages in the boot code (plugin registration etc)
if opts.verbosity>1: os.environ['YADE_DEBUG']='1'
=== modified file 'core/main/yade-batch.in'
--- core/main/yade-batch.in 2010-09-15 10:06:42 +0000
+++ core/main/yade-batch.in 2010-09-18 14:10:30 +0000
@@ -109,6 +109,7 @@
else:
failed=len([j for j in jobs if j.exitStatus!=0])
lastFinished=max([j.finished for j in jobs])
+ # FIXME: do not report sum of runnign time of all jobs, only the real timespan
ret='<p><span style="background-color:%s">Finished</span>, idle for %s, running time %s since %s.</p>'%('red' if failed else 'lime',t2hhmmss(time.time()-lastFinished),t2hhmmss(sum([j.finished-j.started for j in jobs if j.started is not None])),time.ctime(t0))
ret+='<p>Pid %d</p>'%(os.getpid())
ret+='<p>%d slots available, %d used, %d free.</p>'%(maxJobs,usedSlots,maxJobs-usedSlots)
@@ -125,10 +126,10 @@
else:
jobMatch=re.match('/jobs/([0-9]+)/(.*)',self.path)
if not jobMatch:
- self.send_error(404); return
+ self.send_error(404,self.path); return
jobId=int(jobMatch.group(1))
if jobId>=len(jobs):
- self.send_error(404); return
+ self.send_error(404,self.path); return
job=jobs[jobId]
rest=jobMatch.group(2)
if rest=='plots':
@@ -136,7 +137,7 @@
self.sendFile(job.plotsFile,contentType=yade.remote.plotImgMimetype,refresh=(0 if job.status=='DONE' else 5))
elif rest=='log':
if not os.path.exists(job.log):
- self.send_error(404); return
+ self.send_error(404,self.path); return
## once we authenticate properly, send the whole file
## self.sendTextFile(jobs[jobId].log,refresh=5)
## now we have to filter away the cookie
@@ -150,7 +151,7 @@
self.sendPygmentizedFile(job.script,linenostep=5)
elif rest=='table':
self.sendPygmentizedFile(job.table,hl_lines=[job.lineNo],linenostep=1)
- else: self.send_error(404)
+ else: self.send_error(404,self.path)
return
def log_request(self,req): pass
def sendGlobal(self):
@@ -330,8 +331,8 @@
else:
logging.warning('WARNING: job #%d will use %d slots but only %d are available'%(i,nSlots,maxJobs))
# compose command-line: build the hyper-linked variant, then strip HTML tags (ugly, but ensures consistency)
- env='PARAM_TABLE=<a href="/jobs/%d/table">%s:%d</a> DISPLAY= %s '%(i,table,l,' '.join(envVars))
- cmd='%s --threads=%d %s -x <a href="/jobs/%d/script">%s</a>'%(executable,int(nSlots),'--nice=%s'%nice if nice!=None else '',i,simul)
+ env='PARAM_TABLE=<a href="jobs/%d/table">%s:%d</a> DISPLAY= %s '%(i,table,l,' '.join(envVars))
+ cmd='%s --threads=%d %s -x <a href="jobs/%d/script">%s</a>'%(executable,int(nSlots),'--nice=%s'%nice if nice!=None else '',i,simul)
log='> <a href="jobs/%d/log">%s</a> 2>&1'%(i,pipes.quote(logFile))
hrefCmd=env+cmd+log
fullCmd=re.sub('(<a href="[^">]+">|</a>)','',hrefCmd)
=== modified file 'debian/rules'
--- debian/rules 2010-07-24 18:10:24 +0000
+++ debian/rules 2010-09-18 14:10:30 +0000
@@ -69,8 +69,8 @@
check: install
dh_testdir
dh_testroot
- YADE_PREFIX=debian/yade${_VERSION}-dbg/usr debian/yade${_VERSION}-dbg/usr/bin/yade${_VERSION}-dbg -x scripts/regression-tests.py
- YADE_PREFIX=debian/yade${_VERSION}/usr debian/yade${_VERSION}/usr/bin/yade${_VERSION} -x scripts/regression-tests.py
+ YADE_PREFIX=debian/yade${_VERSION}-dbg/usr debian/yade${_VERSION}-dbg/usr/bin/yade${_VERSION}-dbg --test
+ YADE_PREFIX=debian/yade${_VERSION}/usr debian/yade${_VERSION}/usr/bin/yade${_VERSION} --test
# Build architecture-independent files here.
binary-indep: build install
=== modified file 'doc/references.bib'
--- doc/references.bib 2010-09-09 14:02:49 +0000
+++ doc/references.bib 2010-09-18 14:10:30 +0000
@@ -329,3 +329,14 @@
}
+@article{Thornton2000,
+ title={Numerical simulations of deviatoric shear deformation of granular media},
+ author={Colin Thornton},
+ journal={Géotechnique},
+ pages={43--53},
+ volume={50},
+ year={2000},
+ number={1},
+ doi={10.1680/geot.2000.50.1.43}
+}
+
=== modified file 'gui/SConscript'
--- gui/SConscript 2010-08-29 10:39:02 +0000
+++ gui/SConscript 2010-09-18 14:10:30 +0000
@@ -9,10 +9,7 @@
env.File('qt4/SerializableEditor.py'),
env.File('qt4/Inspector.py'),
env.File('qt4/__init__.py'),
- env.SharedLibrary('_GLViewer',['qt4/GLViewer.cpp','qt4/_GLViewer.cpp','qt4/OpenGLManager.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+[env['QGLVIEWER_LIB']],RPATH=env['RPATH']+[env.Literal('\\$$ORIGIN/../../../gui')])
+ env.SharedLibrary('_GLViewer',['qt4/GLViewer.cpp','qt4/_GLViewer.cpp','qt4/OpenGLManager.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+[env['QGLVIEWER_LIB']]+linkPlugins(['PeriodicEngines']),RPATH=env['RPATH']+[env.Literal('\\$$ORIGIN/../../../gui')])
])
env.Command('qt4/img_rc.py','qt4/img.qrc','pyrcc4 -o $buildDir/gui/qt4/img_rc.py gui/qt4/img.qrc')
env.Command('qt4/ui_controller.py','qt4/controller.ui','pyuic4 -o $buildDir/gui/qt4/ui_controller.py gui/qt4/controller.ui')
- #env.Install('$PREFIX/lib/yade$SUFFIX/gui',[
- # env.SharedLibrary('GLViewer',['qt4/GLViewer.cpp'],LIBS=env['LIBS']+[env['QGLVIEWER_LIB']]+linkPlugins(['OpenGLRenderer']))
- #])
=== modified file 'gui/qt4/GLViewer.cpp'
--- gui/qt4/GLViewer.cpp 2010-08-29 10:39:02 +0000
+++ gui/qt4/GLViewer.cpp 2010-09-18 14:10:30 +0000
@@ -20,6 +20,9 @@
#include<boost/algorithm/string.hpp>
#include<boost/version.hpp>
#include<boost/python.hpp>
+#include<sstream>
+#include<iomanip>
+#include<boost/algorithm/string/case_conv.hpp>
#include<yade/lib-serialization/ObjectIO.hpp>
#include<QtGui/qevent.h>
@@ -30,6 +33,36 @@
#include<gl2ps.h>
#endif
+YADE_PLUGIN((SnapshotEngine));
+
+/*****************************************************************************
+*********************************** SnapshotEngine ***************************
+*****************************************************************************/
+
+CREATE_LOGGER(SnapshotEngine);
+void SnapshotEngine::action(){
+ shared_ptr<GLViewer> glv;
+ const int viewNo=0;
+ if(!OpenGLManager::self || ((size_t)viewNo>=OpenGLManager::self->views.size()) || !(glv=OpenGLManager::self->views[viewNo])){
+ if(!ignoreErrors) throw invalid_argument("View #"+lexical_cast<string>(viewNo)+" doesn't exist.");
+ return;
+ }
+ ostringstream fss; fss<<fileBase<<setw(5)<<setfill('0')<<counter++<<"."<<boost::algorithm::to_lower_copy(format);
+ LOG_DEBUG("GL view #"<<viewNo<<" → "<<fss.str())
+ glv->setSnapshotFormat(QString(format.c_str()));
+ glv->nextFrameSnapshotFilename=fss.str();
+ // wait for the renderer to save the frame (will happen at next postDraw)
+ timespec t1,t2; t1.tv_sec=0; t1.tv_nsec=10000000; /* 10 ms */
+ long waiting=0;
+ while(!glv->nextFrameSnapshotFilename.empty()){
+ nanosleep(&t1,&t2);
+ if(((++waiting) % 1000)==0) LOG_WARN("Already waiting "<<waiting/100<<"s for snapshot to be saved. Something went wrong?");
+ }
+ savedSnapshots.push_back(fss.str());
+ usleep((long)(msecSleep*1000));
+}
+
+
CREATE_LOGGER(GLViewer);
GLLock::GLLock(GLViewer* _glv): boost::try_mutex::scoped_lock(Omega::instance().renderMutex), glv(_glv){
@@ -38,29 +71,6 @@
GLLock::~GLLock(){ glv->doneCurrent(); }
-
-// this KILLS performace with qt4!!
-
-#if 0
-// void GLViewer::updateGL(void){/*GLLock lock(this); */QGLViewer::updateGL();}
-/* additionally try: doneCurrent; glFlush; glSwapBuffers after paintGL */
-void GLViewer::paintGL(void){
- /* paintGL encapsulated preDraw, draw and postDraw within QGLViewer. If the mutex cannot be locked,
- * we just return without repainting */
- boost::try_mutex::scoped_try_lock lock(Omega::instance().renderMutex);
- #if BOOST_VERSION<103500
- if(lock.locked())
- #else
- if(lock.owns_lock())
- #endif
- {
- this->makeCurrent();
- QGLViewer::paintGL();
- }
- this->doneCurrent();
-}
-#endif
-
GLViewer::~GLViewer(){ /* get the GL mutex when closing */ GLLock lock(this); /* cerr<<"Destructing view #"<<viewId<<endl;*/ }
void GLViewer::closeEvent(QCloseEvent *e){
=== modified file 'gui/qt4/GLViewer.hpp'
--- gui/qt4/GLViewer.hpp 2010-07-24 18:10:24 +0000
+++ gui/qt4/GLViewer.hpp 2010-09-18 14:10:30 +0000
@@ -11,6 +11,29 @@
#include<boost/date_time/posix_time/posix_time.hpp>
#include<set>
+#include<yade/pkg-common/PeriodicEngines.hpp>
+
+/*****************************************************************************
+*********************************** SnapshotEngine ***************************
+*****************************************************************************/
+/* NOTE: this engine bypasses usual registration of plugins;
+ pyRegisterClass must be called in the module import stanza in gui/qt4/_GLViewer.cpp
+*/
+class SnapshotEngine: public PeriodicEngine{
+ public:
+ virtual void action();
+ YADE_CLASS_BASE_DOC_ATTRS(SnapshotEngine,PeriodicEngine,"Periodically save snapshots of GLView(s) as .png files. Files are named *fileBase*+*counter*+'.png' (counter is left-padded by 0s, i.e. snap00004.png).",
+ ((string,format,"PNG",,"Format of snapshots (one of JPEG, PNG, EPS, PS, PPM, BMP) `QGLViewer documentation <http://www.libqglviewer.com/refManual/classQGLViewer.html#abbb1add55632dced395e2f1b78ef491c>`_. File extension will be lowercased *format*. Validity of format is not checked."))
+ ((string,fileBase,"",,"Basename for snapshots"))
+ ((int,counter,0,,"Number that will be appended to fileBase when the next snapshot is saved (incremented at every save). |yupdate|"))
+ // ((int,viewNo,((void)"primary view",0),"The GLView number that we save."))
+ ((bool,ignoreErrors,true,,"Silently return if selected view doesn't exist"))
+ ((vector<string>,savedSnapshots,,,"Files that have been created so far"))
+ ((int,msecSleep,0,,"number of msec to sleep after snapshot (to prevent 3d hw problems) [ms]"))
+ );
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(SnapshotEngine);
/*! Class handling user interaction with the openGL rendering of simulation.
*
@@ -148,3 +171,6 @@
virtual float zNear() const;
virtual void setCuttingDistance(float s){cuttingDistance=s;};
};
+
+
+
=== modified file 'gui/qt4/_GLViewer.cpp'
--- gui/qt4/_GLViewer.cpp 2010-08-29 10:39:02 +0000
+++ gui/qt4/_GLViewer.cpp 2010-09-18 14:10:30 +0000
@@ -86,7 +86,9 @@
OpenGLManager* glm=new OpenGLManager(); // keep this singleton object forever
glm->emitStartTimer();
-
+
+ // HACK: register SnapshotEngine here
+ SnapshotEngine se; se.pyRegisterClass(py::scope());
py::def("View",createView,"Create a new 3d view.");
py::def("center",centerViews,"Center all views.");
=== added file 'pkg/common/RenderingEngine/Gl1_NormPhys.cpp'
--- pkg/common/RenderingEngine/Gl1_NormPhys.cpp 1970-01-01 00:00:00 +0000
+++ pkg/common/RenderingEngine/Gl1_NormPhys.cpp 2010-09-18 14:10:30 +0000
@@ -0,0 +1,60 @@
+
+#include<yade/core/Scene.hpp>
+#include<yade/pkg-common/Gl1_NormPhys.hpp>
+#include<yade/pkg-common/NormShearPhys.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
+#include<yade/pkg-dem/Shop.hpp>
+
+YADE_PLUGIN((Gl1_NormPhys));
+YADE_REQUIRE_FEATURE(OPENGL);
+
+GLUquadric* Gl1_NormPhys::gluQuadric=NULL;
+Real Gl1_NormPhys::maxFn;
+Real Gl1_NormPhys::refRadius;
+Real Gl1_NormPhys::maxRadius;
+int Gl1_NormPhys::signFilter;
+int Gl1_NormPhys::slices;
+int Gl1_NormPhys::stacks;
+
+void Gl1_NormPhys::go(const shared_ptr<InteractionPhysics>& ip, const shared_ptr<Interaction>& i, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
+ if(!gluQuadric){ gluQuadric=gluNewQuadric(); if(!gluQuadric) throw runtime_error("Gl1_NormPhys::go unable to allocate new GLUquadric object (out of memory?)."); }
+ NormPhys* np=static_cast<NormPhys*>(ip.get());
+ shared_ptr<InteractionGeometry> ig(i->interactionGeometry); if(!ig) return; // changed meanwhile?
+ GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*>(ig.get());
+ //if(!gsc) cerr<<"Gl1_NormPhys: InteractionGeometry is not a GenericSpheresContact, but a "<<ig->getClassName()<<endl;
+ Real fnNorm=np->normalForce.dot(gsc->normal);
+ if((signFilter>0 && fnNorm<0) || (signFilter<0 && fnNorm>0)) return;
+ int fnSign=fnNorm>0?1:-1;
+ fnNorm=abs(fnNorm);
+ maxFn=max(fnNorm,maxFn);
+ Real realMaxRadius;
+ if(maxRadius<0){
+ if(gsc->refR1>0) refRadius=min(gsc->refR1,refRadius);
+ if(gsc->refR2>0) refRadius=min(gsc->refR2,refRadius);
+ realMaxRadius=refRadius;
+ }
+ else realMaxRadius=maxRadius;
+ Real radius=realMaxRadius*(fnNorm/maxFn); // use logarithmic scale here?
+ Vector3r color=Shop::scalarOnColorScale(fnNorm*fnSign,-maxFn,maxFn);
+ Vector3r p1=b1->state->pos, p2=b2->state->pos;
+ Vector3r relPos;
+ if(scene->isPeriodic){
+ relPos=p2+scene->cell->Hsize*i->cellDist.cast<Real>()-p1;
+ p1=scene->cell->wrapShearedPt(p1);
+ p2=p1+relPos;
+ } else {
+ relPos=p2-p1;
+ }
+ Real dist=relPos.norm();
+
+ glDisable(GL_CULL_FACE);
+ glPushMatrix();
+ glTranslatef(p1[0],p1[1],p1[2]);
+ Quaternionr q(Quaternionr().setFromTwoVectors(Vector3r(0,0,1),relPos/dist /* normalized */));
+ // using Transform with OpenGL: http://eigen.tuxfamily.org/dox/TutorialGeometry.html
+ glMultMatrixd(Eigen::Transform3d(q).data());
+ glColor3v(color);
+ gluCylinder(gluQuadric,radius,radius,dist,slices,stacks);
+ glPopMatrix();
+}
+
=== added file 'pkg/common/RenderingEngine/Gl1_NormPhys.hpp'
--- pkg/common/RenderingEngine/Gl1_NormPhys.hpp 1970-01-01 00:00:00 +0000
+++ pkg/common/RenderingEngine/Gl1_NormPhys.hpp 2010-09-18 14:10:30 +0000
@@ -0,0 +1,24 @@
+#pragma once
+
+#include<yade/lib-opengl/OpenGLWrapper.hpp>
+#include<yade/lib-opengl/GLUtils.hpp>
+#include<yade/pkg-common/GLDrawFunctors.hpp>
+#include<GL/glu.h>
+
+
+class Gl1_NormPhys: public GlInteractionPhysicsFunctor{
+ static GLUquadric* gluQuadric; // needed for gluCylinder, initialized by ::go if no initialized yet
+ public:
+ virtual void go(const shared_ptr<InteractionPhysics>&,const shared_ptr<Interaction>&,const shared_ptr<Body>&,const shared_ptr<Body>&,bool wireFrame);
+ YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_NormPhys,GlInteractionPhysicsFunctor,"Renders :yref:`NormPhys` objects as cylinders of which diameter and color depends on :yref:`NormPhys:normForce` magnitude.",
+ ((Real,maxFn,0,,"Value of :yref:`NormPhys.normalForce` corresponding to :yref:`maxDiameter<Gl1_NormPhys.maxDiameter>`. This value will be increased (but *not decreased*) automatically."))
+ ((int,signFilter,0,,"If non-zero, only display contacts with negative (-1) or positive (+1) normal forces; if zero, all contacts will be displayed."))
+ ((Real,refRadius,std::numeric_limits<Real>::infinity(),,"Reference (minimum) particle radius; used only if :yref:`maxRadius<Gl1_NormPhys.maxRadius>` is negative. This value will be decreased (but *not increased*) automatically. |yupdate|"))
+ ((Real,maxRadius,-1,,"Cylinder radius corresponding to the maximum normal force. If negative, auto-updated :yref:`refRadius<Gl1_NormPhys.refRadius>` will be used instead."))
+ ((int,slices,6,,"Number of sphere slices; (see `glutCylinder reference <http://www.opengl.org/sdk/docs/man/xhtml/gluCylinder.xml>`__)"))
+ ((int,stacks,1,,"Number of sphere stacks; (see `glutCylinder reference <http://www.opengl.org/sdk/docs/man/xhtml/gluCylinder.xml>`__)"))
+ );
+ RENDERS(NormPhys);
+};
+REGISTER_SERIALIZABLE(Gl1_NormPhys);
+
=== modified file 'pkg/common/RenderingEngine/OpenGLRenderer.cpp'
--- pkg/common/RenderingEngine/OpenGLRenderer.cpp 2010-08-15 05:27:53 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderer.cpp 2010-09-18 14:10:30 +0000
@@ -155,6 +155,11 @@
scene=_scene;
// assign scene inside functors
+ boundDispatcher.updateScenePtr();
+ interactionGeometryDispatcher.updateScenePtr();
+ interactionPhysicsDispatcher.updateScenePtr();
+ shapeDispatcher.updateScenePtr();
+ stateDispatcher.updateScenePtr();
// just to make sure, since it is not initialized by default
if(!scene->bound) scene->bound=shared_ptr<Aabb>(new Aabb);
@@ -203,9 +208,7 @@
if (dof || id) renderDOF_ID();
if (bound) renderBoundingVolume();
- if (shape){
- renderShape();
- }
+ if (shape) renderShape();
if (intrAllWire) renderAllInteractionsWire();
if (intrGeom) renderInteractionGeometry();
if (intrPhys) renderInteractionPhysics();
=== modified file 'pkg/dem/DataClass/SpherePack.cpp'
--- pkg/dem/DataClass/SpherePack.cpp 2010-09-10 11:33:21 +0000
+++ pkg/dem/DataClass/SpherePack.cpp 2010-09-18 14:10:30 +0000
@@ -216,22 +216,26 @@
// New code to include the psd giving few points of it
const float pi = 3.1415926;
-long SpherePack::particleSD(Vector3r mn, Vector3r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing){
- Real Vtot=numSph*4./3.*pi*pow(rMean,3.); // total volume of the packing (computed with rMean)
-
- // calculate number of spheres necessary per each radius to match the wanted psd
- // passing has to contain increasing values
+long SpherePack::particleSD(Vector3r mn, Vector3r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount){
vector<Real> numbers;
- for (size_t i=0; i<radii.size(); i++){
- Real volS=4./3.*pi*pow(radii[i],3.);
- if (i==0) {numbers.push_back(passing[i]/100.*Vtot/volS);}
- else {numbers.push_back((passing[i]-passing[i-1])/100.*Vtot/volS);} //
+ if(!passingIsNotPercentageButCount){
+ Real Vtot=numSph*4./3.*pi*pow(rMean,3.); // total volume of the packing (computed with rMean)
- cout<<"numbers["<<i<<"] = "<<numbers[i]<<endl;
- cout<<"radii["<<i<<"] = "<<radii[i]<<endl;
- cout<<"vol tot = "<<Vtot<<endl;
- cout<<"v_sphere = "<<volS<<endl;
- cout<<"passing["<<i<<"] = "<<passing[i]<<endl;
+ // calculate number of spheres necessary per each radius to match the wanted psd
+ // passing has to contain increasing values
+ for (size_t i=0; i<radii.size(); i++){
+ Real volS=4./3.*pi*pow(radii[i],3.);
+ if (i==0) {numbers.push_back(passing[i]/100.*Vtot/volS);}
+ else {numbers.push_back((passing[i]-passing[i-1])/100.*Vtot/volS);} //
+
+ cout<<"numbers["<<i<<"] = "<<numbers[i]<<endl;
+ cout<<"radii["<<i<<"] = "<<radii[i]<<endl;
+ cout<<"vol tot = "<<Vtot<<endl;
+ cout<<"v_sphere = "<<volS<<endl;
+ cout<<"passing["<<i<<"] = "<<passing[i]<<endl;
+ }
+ } else {
+ FOREACH(Real p, passing) numbers.push_back(p);
}
static boost::minstd_rand randGen(TimingInfo::getNow(true));
=== modified file 'pkg/dem/DataClass/SpherePack.hpp'
--- pkg/dem/DataClass/SpherePack.hpp 2010-09-09 14:02:49 +0000
+++ pkg/dem/DataClass/SpherePack.hpp 2010-09-18 14:10:30 +0000
@@ -61,7 +61,7 @@
int psdGetPiece(Real x, const vector<Real>& cumm, Real& norm);
// function to generate the packing based on a given psd
- long particleSD(Vector3r min, Vector3r max, Real rMean, bool periodic=false, string name="", int numSph=400, const vector<Real>& radii=vector<Real>(), const vector<Real>& passing=vector<Real>());
+ long particleSD(Vector3r min, Vector3r max, Real rMean, bool periodic=false, string name="", int numSph=400, const vector<Real>& radii=vector<Real>(), const vector<Real>& passing=vector<Real>(),bool passingIsNotPercentageButCount=false);
// interpolate a variable with power distribution (exponent -3) between two margin values, given uniformly distributed x∈(0,1)
Real pow3Interp(Real x,Real a,Real b){ return pow(x*(pow(b,-2)-pow(a,-2))+pow(a,-2),-1./2); }
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-09-10 11:41:37 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-09-18 14:10:30 +0000
@@ -166,7 +166,7 @@
Real Fn = phys->kno*std::pow(uN,1.5); // normal Force (scalar)
if (includeAdhesion) {
Fn -= phys->adhesionForce; // include adhesion force to account for the effect of Van der Waals interactions
- if (Fn < 0.0) {phys->isAdhesive = true;} // set true the bool to count the number of adhesive contacts
+ phys->isAdhesive = (Fn<0); // set true the bool to count the number of adhesive contacts
}
phys->normalForce = Fn*scg->normal; // normal Force (vector)
=== modified file 'pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp 2010-09-15 10:06:42 +0000
+++ pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp 2010-09-18 14:10:30 +0000
@@ -181,11 +181,21 @@
rpmSpecDiam->SetName("rpmSpecDiam");
if(recActive[REC_INTR]){
+ // holds information about cell distance between spatial and displayed position of each particle
+ vector<Vector3i> wrapCellDist; if (scene->isPeriodic){ wrapCellDist.resize(scene->bodies->size()); }
// save body positions, referenced by ids by vtkLine
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if (!b) { /* we must keep ids contiguous */ intrBodyPos->InsertNextPoint(NaN,NaN,NaN); continue; }
- const Vector3r& pos=b->state->pos;
- intrBodyPos->InsertNextPoint(pos[0],pos[1],pos[2]);
+ if (!b) {
+ /* must keep ids contiguous, so that position in the array corresponds to Body::id */
+ intrBodyPos->InsertNextPoint(NaN,NaN,NaN);
+ continue;
+ }
+ if(!scene->isPeriodic){ intrBodyPos->InsertNextPoint(b->state->pos[0],b->state->pos[1],b->state->pos[2]); }
+ else {
+ Vector3r pos=scene->cell->wrapShearedPt(b->state->pos,wrapCellDist[b->id]);
+ intrBodyPos->InsertNextPoint(pos[0],pos[1],pos[2]);
+ }
+ assert(intrBodyPos->GetNumberOfPoints()==b->id+1);
}
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
@@ -193,15 +203,47 @@
if(!(dynamic_cast<Sphere*>(Body::byId(I->getId1())->shape.get()))) continue;
if(!(dynamic_cast<Sphere*>(Body::byId(I->getId2())->shape.get()))) continue;
}
- vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
- line->GetPointIds()->SetId(0,I->getId1());
- line->GetPointIds()->SetId(1,I->getId2());
- intrCells->InsertNextCell(line);
+ /* For the periodic boundary conditions,
+ find out whether the interaction crosses the boundary of the periodic cell;
+ if it does, display the interaction on both sides of the cell, with one of the
+ points sticking out in each case.
+ Since vtkLines must connect points with an ID assigned, we will create a new additional
+ point for each point outside the cell. It might create some data redundancy, but
+ let us suppose that the number of interactions crossing the cell boundary is low compared
+ to total numer of interactions
+ */
+ // how many times to add values defined on interactions, depending on how many times the interaction is saved
+ int numAddValues=1;
+ // aperiodic boundary, or interaction is inside the cell
+ if(!scene->isPeriodic || (scene->isPeriodic && (I->cellDist==wrapCellDist[I->getId2()]-wrapCellDist[I->getId1()]))){
+ vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
+ line->GetPointIds()->SetId(0,I->getId1());
+ line->GetPointIds()->SetId(1,I->getId2());
+ intrCells->InsertNextCell(line);
+ } else {
+ assert(scene->isPeriodic);
+ // spatial positions of particles
+ const Vector3r& p01(Body::byId(I->getId1())->state->pos); const Vector3r& p02(Body::byId(I->getId2())->state->pos);
+ // create two line objects; each of them has one endpoint inside the cell and the other one sticks outside
+ // A,B are the "fake" bodies outside the cell for id1 and id2 respectively, p1,p2 are the displayed points
+ // distance in cell units for shifting A away from p1; negated value is shift of B away from p2
+ Vector3r ptA(p01+scene->cell->Hsize*(wrapCellDist[I->getId2()]-I->cellDist).cast<Real>());
+ Vector3r ptB(p02+scene->cell->Hsize*(wrapCellDist[I->getId1()]-I->cellDist).cast<Real>());
+ vtkIdType idPtA=intrBodyPos->InsertNextPoint(ptA[0],ptA[1],ptA[2]), idPtB=intrBodyPos->InsertNextPoint(ptB[0],ptB[1],ptB[2]);
+ vtkSmartPointer<vtkLine> line1B(vtkSmartPointer<vtkLine>::New()); line1B->GetPointIds()->SetId(0,I->getId1()); line1B->GetPointIds()->SetId(1,idPtB);
+ vtkSmartPointer<vtkLine> lineA2(vtkSmartPointer<vtkLine>::New()); lineA2->GetPointIds()->SetId(0,idPtA); line1B->GetPointIds()->SetId(1,I->getId2());
+ numAddValues=2;
+ }
const NormShearPhys* phys = YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
- float fn[3]={abs(phys->normalForce[0]),abs(phys->normalForce[1]),abs(phys->normalForce[2])};
+ const GenericSpheresContact* geom = YADE_CAST<GenericSpheresContact*>(I->interactionGeometry.get());
+ // gives _signed_ scalar of normal force, following the convention used in the respective constitutive law
+ float fn=phys->normalForce.dot(geom->normal);
float fs[3]={abs(phys->shearForce[0]),abs(phys->shearForce[1]),abs(phys->shearForce[2])};
- intrForceN->InsertNextTupleValue(fn);
- intrAbsForceT->InsertNextTupleValue(fs);
+ // add the value once for each interaction object that we created (might be 2 for the periodic boundary)
+ for(int i=0; i<numAddValues; i++){
+ intrForceN->InsertNextValue(fn);
+ intrAbsForceT->InsertNextTupleValue(fs);
+ }
}
}
@@ -217,7 +259,7 @@
if (sphere){
if(skipNondynamic && !b->isDynamic()) continue;
vtkIdType pid[1];
- const Vector3r& pos = b->state->pos;
+ Vector3r pos(scene->isPeriodic ? scene->cell->wrapShearedPt(b->state->pos) : b->state->pos);
pid[0] = spheresPos->InsertNextPoint(pos[0], pos[1], pos[2]);
spheresCells->InsertNextCell(1,pid);
radii->InsertNextValue(sphere->radius);
@@ -273,13 +315,13 @@
if (recActive[REC_FACETS]){
const Facet* facet = dynamic_cast<Facet*>(b->shape.get());
if (facet){
- const Se3r& O = b->state->se3;
+ Vector3r pos(scene->isPeriodic ? scene->cell->wrapShearedPt(b->state->pos) : b->state->pos);
const vector<Vector3r>& localPos = facet->vertices;
- Matrix3r facetAxisT=O.orientation.toRotationMatrix();
+ Matrix3r facetAxisT=b->state->ori.toRotationMatrix();
vtkSmartPointer<vtkTriangle> tri = vtkSmartPointer<vtkTriangle>::New();
vtkIdType nbPoints=facetsPos->GetNumberOfPoints();
for (int i=0;i<3;++i){
- Vector3r globalPos = O.position + facetAxisT * localPos[i];
+ Vector3r globalPos = pos + facetAxisT * localPos[i];
facetsPos->InsertNextPoint(globalPos[0], globalPos[1], globalPos[2]);
tri->GetPointIds()->SetId(i,nbPoints+i);
}
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2010-09-12 08:44:00 +0000
+++ pkg/dem/meta/Shop.cpp 2010-09-18 14:10:30 +0000
@@ -1204,6 +1204,42 @@
}
}
+/* Return the stress tensor decomposed in 2 contributions, from normal and shear forces.
+The formulation follows the [Thornton2000]_ article
+"Numerical simulations of deviatoric shear deformation of granular media", eq (3) and (4)
+ */
+py::tuple Shop::normalShearStressTensors(bool compressionPositive){
+ Scene* scene=Omega::instance().getScene().get();
+ if (!scene->isPeriodic){ throw runtime_error("Can't compute stress of periodic cell in aperiodic simulation."); }
+ Matrix3r sigN(Matrix3r::Zero()),sigT(Matrix3r::Zero());
+ const Matrix3r& cellHsize(scene->cell->Hsize);
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->interactionGeometry.get());
+ NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
+ const Vector3r& n=geom->normal;
+ // if compression has positive sign, we need to change sign for both normal and shear force
+ // make copy to Fs since shearForce is used at multiple places (less efficient, but avoids confusion)
+ Vector3r Fs=(compressionPositive?-1:1)*phys->shearForce;
+ Real N=(compressionPositive?-1:1)*phys->normalForce.dot(n), T=Fs.norm();
+ bool hasShear=(T>0);
+ Vector3r t; if(hasShear) t=Fs/T;
+ // Real R=(Body::byId(I->getId2(),scene)->state->pos+cellHsize*I->cellDist.cast<Real>()-Body::byId(I->getId1(),scene)->state->pos).norm();
+ Real R=.5*(geom->refR1+geom->refR2);
+ for(int i=0; i<3; i++) for(int j=i; j<3; j++){
+ sigN(i,j)+=R*N*n[i]*n[j];
+ if(hasShear) sigT(i,j)+=R*T*n[i]*t[j];
+ }
+ }
+ Real vol=scene->cell->getVolume();
+ sigN*=2/vol; sigT*=2/vol;
+ // fill terms under the diagonal
+ sigN(1,0)=sigN(0,1); sigN(2,0)=sigN(0,2); sigN(2,1)=sigN(1,2);
+ sigT(1,0)=sigT(0,1); sigT(2,0)=sigT(0,2); sigT(2,1)=sigT(1,2);
+ return py::make_tuple(sigN,sigT);
+}
+
+
Matrix3r Shop::stressTensorOfPeriodicCell(bool smallStrains){
Scene* scene=Omega::instance().getScene().get();
if (!scene->isPeriodic){ throw runtime_error("Can't compute stress of periodic cell in aperiodic simulation."); }
=== modified file 'pkg/dem/meta/Shop.hpp'
--- pkg/dem/meta/Shop.hpp 2010-09-12 08:44:00 +0000
+++ pkg/dem/meta/Shop.hpp 2010-09-18 14:10:30 +0000
@@ -24,6 +24,7 @@
using namespace std;
using boost::shared_ptr;
+namespace py = boost::python;
/*! Miscillaneous utility functions which are believed to be generally useful.
*
@@ -132,4 +133,7 @@
//! Function to compute overall ("macroscopic") stress of periodic cell
static Matrix3r stressTensorOfPeriodicCell(bool smallStrains=true);
+ //! Compute overall ("macroscopic") stress of periodic cell, returning 2 tensors
+ //! (contribution of normal and shear forces)
+ static py::tuple normalShearStressTensors(bool compressionPositive=false);
};
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2010-09-15 10:06:42 +0000
+++ py/_utils.cpp 2010-09-18 14:10:30 +0000
@@ -423,6 +423,7 @@
Real Shop__getPorosity(Real volume=-1){ return Shop::getPorosity(Omega::instance().getScene(),volume); }
Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
+py::tuple Shop__normalShearStressTensors(bool compressionPositive=false){return Shop::normalShearStressTensors(compressionPositive);}
BOOST_PYTHON_MODULE(_utils){
// http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
@@ -466,5 +467,6 @@
py::def("flipCell",&Shop::flipCell,(py::arg("flip")=Matrix3r(Matrix3r::Zero())),"Flip periodic cell so that angles between $R^3$ axes and transformed axes are as small as possible. This function relies on the fact that periodic cell defines by repetition or its corners regular grid of points in $R^3$; however, all cells generating identical grid are equivalent and can be flipped one over another. This necessiatates adjustment of :yref:`Interaction.cellDist` for interactions that cross boundary and didn't before (or vice versa), and re-initialization of collider. The *flip* argument can be used to specify desired flip: integers, each column for one axis; if zero matrix, best fit (minimizing the angles) is computed automatically.\n\nIn c++, this function is accessible as ``Shop::flipCell``.\n\n.. warning::\n\t This function is currently broken and should not be used.");
py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("m"),py::arg("tc"),py::arg("en"),py::arg("es")),"Get viscoelastic interaction parameters from analytical solution of a pair spheres collision problem. \n\n:Parameters:\n\t`m` : float\n\t\tsphere mass\n\t`tc` : float\n\t\tcollision time\n\t`en` : float\n\t\tnormal restitution coefficient\n\t`es` : float\n\t\ttangential restitution coefficient.\n\n\n:return: \n\tdict with keys:\n\n\tkn : float\n\t\tnormal elastic coefficient computed as:\n\n.. math:: k_n=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\n\n\tcn : float\n\t\tnormal viscous coefficient computed as:\n\n.. math:: c_n=-\\frac{2m}{t_c}\\ln e_n\n\n\n\tkt : float\n\t\ttangential elastic coefficient computed as:\n\n.. math:: k_t=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)\n\n\tct : float\n\t\ttangential viscous coefficient computed as:\n\n.. math:: c_t=-\\frac27\\frac{m}{t_c}\\ln e_t.\n\n\nFor details see [Pournin2001]_. ");
py::def("stressTensorOfPeriodicCell",Shop__stressTensorOfPeriodicCell,(py::args("smallStrains")=false),"Compute overall (macroscopic) stress of periodic cell using equation published in [Kuhl2001]_:\n\n.. math:: \\vec{\\sigma}=\\frac{1}{V}\\sum_cl^c[\\vec{N}^cf_N^c+\\vec{T}^{cT}\\cdot\\vec{f}^c_T],\n\nwhere $V$ is volume of the cell, $l^c$ length of interaction $c$, $f^c_N$ normal force and $\\vec{f}^c_T$ shear force. Sumed are values over all interactions $c$. $\\vec{N}^c$ and $\\vec{T}^{cT}$ are projection tensors (see the original publication for more details):\n\n.. math:: \\vec{N}=\\vec{n}\\otimes\\vec{n}\\rightarrow N_{ij}=n_in_j\n\n.. math:: \\vec{T}^T=\\vec{I}_{sym}\\cdot\\vec{n}-\\vec{n}\\otimes\\vec{n}\\otimes\\vec{n}\\rightarrow T^T_{ijk}=\\frac{1}{2}(\\delta_{ik}\\delta_{jl}+\\delta_{il}\\delta_{jk})n_l-n_in_jn_k\n\n.. math:: \\vec{T}^T\\cdot\\vec{f}_T\\equiv T^T_{ijk}f_k=(\\delta_{ik}n_j/2+\\delta_{jk}n_i/2-n_in_jn_k)f_k=n_jf_i/2+n_if_j/2-n_in_jn_kf_k,\n\nwhere $n$ is unit vector oriented along the interaction (:yref:`normal<GenericSpheresContact::normal>`) and $\\delta$ is Kronecker's delta. As $\\vec{n}$ and $\\vec{f}_T$ are perpendicular (therfore $n_if_i=0$) we can write\n\n.. math:: \\sigma_{ij}=\\frac{1}{V}\\sum l[n_in_jf_N+n_jf^T_i/2+n_if^T_j/2]\n\n:param bool smallStrains: if false (large strains), real values of volume and interaction lengths are computed. If true, only :yref:`refLength<Dem3DofGeom::refLength>` of interactions and initial volume are computed (can save some time).\n\n:return: macroscopic stress tensor as Matrix3");
+ py::def("normalShearStressTensors",Shop__normalShearStressTensors,(py::args("compressionPositive")=false),"Compute overall stress tensor of the periodic cell decomposed in 2 parts, one contributed by normal forces, the other by shear forces. The formulation can be found in [Thornton2000]_, eq. (3):\n\n.. math:: \\tensor{sigma}_{ij}=\\frac{2}{V}\\sum R N \\vec{n}_i \\vec{n}_j+\frac{2}{V}\\sum R T \\vec{n}_i\\vec{t}_j\n\nwhere $V$ is the cell volume, $R$ is \"contact radius\" (in our implementation, current distance between particle centroids), $\\vec{n}$ is the normal vector, $\\vec{t}$ is a vector perpendicular to $\\vec{n}$, $N$ and $T$ are norms of normal and shear forces.");
py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap ration in interactions (with :yref:`ScGeom`) of two :yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1 r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$, $r_2$ are radii of the two spheres in contact.");
}
=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp 2010-07-28 06:50:31 +0000
+++ py/mathWrap/miniEigen.cpp 2010-09-18 14:10:30 +0000
@@ -97,6 +97,7 @@
};
static Matrix3r* Matrix3r_fromElements(Real m00, Real m01, Real m02, Real m10, Real m11, Real m12, Real m20, Real m21, Real m22){ Matrix3r* m(new Matrix3r); (*m)<<m00,m01,m02,m10,m11,m12,m20,m21,m22; return m; }
+static Vector3r Matrix3r_diagonal(const Matrix3r& m){ return Vector3r(m.diagonal()); }
static Quaternionr Quaternionr_setFromTwoVectors(Quaternionr& q, const Vector3r& u, const Vector3r& v){ return q.setFromTwoVectors(u,v); }
static Vector3r Quaternionr_Rotate(Quaternionr& q, const Vector3r& u){ return q*u; }
// supposed to return raw pointer (or auto_ptr), boost::python takes care of the lifetime management
@@ -188,6 +189,7 @@
.def("inverse",&Matrix3r_inverse)
.def("transpose",&Matrix3r_transpose)
.def("polarDecomposition",&Matrix3r_polarDecomposition)
+ .def("diagonal",&Matrix3r_diagonal)
//
.def("__neg__",&Matrix3r__neg__)
=== modified file 'py/pack/_packSpheres.cpp'
--- py/pack/_packSpheres.cpp 2010-09-09 14:02:49 +0000
+++ py/pack/_packSpheres.cpp 2010-09-18 14:10:30 +0000
@@ -21,7 +21,7 @@
"\n\n:param Vector3 minCorner: lower corner of the box\n:param Vector3 maxCorner: upper corner of the box\n:param float rMean: mean radius or spheres\n:param float rRelFuzz: dispersion of radius relative to rMean\n:param int num: number of spheres to be generated (if negavite, generate as many as possible, ending after a fixed number of tries to place the sphere in space)\n:param bool periodic: whether the packing to be generated should be periodic\n:param float porosity: if specified, estimate mean radius $r_m$ (*rMean* must not be given) using *rRelFuzz* ($z$) and *num* ($N$) so that the porosity given ($\\rho$) is approximately achieved at the end of generation, $r_m=\\sqrt[3]{\\frac{V(1-\\rho)}{\\frac{4}{3}\\pi(1+z^2)N}}$. The value of $\\rho$=0.65 is recommended.\n:param psdSizes: sieve sizes (particle diameters) when particle size distribution (PSD) is specified\n:param psdCumm: cummulative fractions of particle sizes given by *psdSizes*; must be the same length as *psdSizes* and should be non-decreasing\n:param bool distributeMass: if ``True``, given distribution will be used to distribute sphere's mass rather than radius of them.\n:returns: number of created spheres, which can be lower than *num* is the packing is too tight.\n")
.def("psd",&SpherePack::psd,(python::arg("bins")=10,python::arg("mass")=false),"Return `particle size distribution <http://en.wikipedia.org/wiki/Particle_size_distribution>`__ of the packing.\n:param int bins: number of bins between minimum and maximum diameter\n:param mass: Compute relative mass rather than relative particle count for each bin. Corresponds to :yref:`distributeMass parameter for makeCloud<yade.pack.SpherePack.makeCloud>`.\n:returns: tuple of ``(cumm,edges)``, where ``cumm`` are cummulative fractions for respective diameters and ``edges`` are those diameter values. Dimension of both arrays is equal to ``bins+1``.")
// new psd
- .def("particleSD",&SpherePack::particleSD,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("rMean"),python::arg("periodic")=false,python::arg("name"),python::arg("numSph"),python::arg("radii")=vector<Real>(),python::arg("passing")=vector<Real>()),"Create random packing enclosed in box given by minCorner and maxCorner, containing numSph spheres. Returns number of created spheres, which can be < num if the packing is too tight. The computation is done according to the given psd.")
+ .def("particleSD",&SpherePack::particleSD,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("rMean"),python::arg("periodic")=false,python::arg("name"),python::arg("numSph"),python::arg("radii")=vector<Real>(),python::arg("passing")=vector<Real>(),python::arg("passingIsNotPercentageButCount")=false),"Create random packing enclosed in box given by minCorner and maxCorner, containing numSph spheres. Returns number of created spheres, which can be < num if the packing is too tight. The computation is done according to the given psd.")
//
.def("aabb",&SpherePack::aabb_py,"Get axis-aligned bounding box coordinates, as 2 3-tuples.")
.def("dim",&SpherePack::dim,"Return dimensions of the packing in terms of aabb(), as a 3-tuple.")
=== modified file 'py/plot.py'
--- py/plot.py 2010-09-15 10:06:42 +0000
+++ py/plot.py 2010-09-18 14:10:30 +0000
@@ -52,7 +52,7 @@
"Interval for the live plot updating, in seconds."
autozoom=True
"Enable/disable automatic plot rezooming after data update."
-scientific=True
+scientific=True if hasattr(pylab,'ticklabel_format') else False ## safe default for older matplotlib versions
"Use scientific notation for axes ticks."
axesWd=0
"Linewidth (in points) to make *x* and *y* axes better visible; not activated if non-positive."
@@ -228,6 +228,7 @@
You can use
>>> from yade import plot
+ >>> plot.plots={'foo':('bar',)}
>>> plot.plot(noShow=True).savefig('someFile.pdf')
>>> import os
>>> os.path.exists('someFile.pdf')
=== modified file 'py/tests/omega.py'
--- py/tests/omega.py 2010-08-24 12:54:14 +0000
+++ py/tests/omega.py 2010-09-18 14:10:30 +0000
@@ -19,6 +19,15 @@
class TestIO(unittest.TestCase): pass
class TestTags(unittest.TestCase): pass
+
+class TestCell(unittest.TestCase):
+ def setUp(self):
+ O.reset(); O.periodic=True
+ def testAttributesAreCrossUpdated(self):
+ "updates Hsize automatically when refSize is updated"
+ O.cell.refSize=(2.55,11,45)
+ self.assert_(O.cell.Hsize==Matrix3(2.55,0,0, 0,11,0, 0,0,45));
+
class TestMaterialStateAssociativity(unittest.TestCase):
def setUp(self): O.reset()
def testThrowsAtBadCombination(self):
=== modified file 'py/utils.py'
--- py/utils.py 2010-09-15 10:06:42 +0000
+++ py/utils.py 2010-09-18 14:10:30 +0000
@@ -466,13 +466,31 @@
color=(random.random(),random.random(),random.random())
if b.dynamic or not onlyDynamic: b.shape.color=color
-def avgNumInteractions(cutoff=0.):
- """Return average numer of interactions per particle, also known as *coordination number*.
+def avgNumInteractions(cutoff=0.,skipFree=False):
+ r"""Return average numer of interactions per particle, also known as *coordination number* $Z$. This number is defined as
+
+ .. math:: Z=2C/N
+
+ where $C$ is number of contacts and $N$ is number of particles.
+
+ With *skipFree*, particles not contributing to stable state of the packing are skipped, following equation (8) given in [Thornton2000]_:
+
+ .. math:: Z_m=\frac{2C-N_1}{N-N_0-N_1}
:param cutoff: cut some relative part of the sample's bounding box away.
- """
- nums,counts=bodyNumInteractionsHistogram(aabbExtrema(cutoff))
- return sum([nums[i]*counts[i] for i in range(len(nums))])/(1.*sum(counts))
+ :param skipFree: """
+ if cutoff==0 and not skipFree: return 2*O.interactions.countReal()*1./len(O.bodies)
+ else:
+ nums,counts=bodyNumInteractionsHistogram(aabbExtrema(cutoff))
+ ## CC is 2*C
+ CC=sum([nums[i]*counts[i] for i in range(len(nums))]); N=sum(counts)
+ if not skipFree: return CC*1./N
+ ## find bins with 0 and 1 spheres
+ N0=0 if (0 not in nums) else counts[nums.index(0)]
+ N1=0 if (1 not in nums) else counts[nums.index(1)]
+ NN=N-N0-N1
+ return (CC-N1)*1./NN if NN>0 else float('nan')
+
def plotNumInteractionsHistogram(cutoff=0.):
"Plot histogram with number of interactions per body, optionally cutting away *cutoff* relative axis-aligned box from specimen margin."
@@ -512,51 +530,31 @@
else: pylab.show()
-
-
-
-def encodeVideoFromFrames(frameSpec,out,renameNotOverwrite=True,fps=24):
- """Create .ogg video from external image files.
-
- :parameters:
- `frameSpec`: wildcard | sequence of filenames
- If string, wildcard in format understood by GStreamer's multifilesrc plugin (e.g. '/tmp/frame-%04d.png'). If list or tuple, filenames to be encoded in given order.
-
- .. warning::
- GStreamer is picky about the wildcard; if you pass a wrong one, if will not complain, but silently stall.
- `out`: filename
- file to save video into
- `renameNotOverwrite`: bool
- if True, existing same-named video file will have ~[number] appended; will be overwritten otherwise.
- `fps`:
- Frames per second.
+def encodeVideoFromFrames(*args,**kw):
+ "|ydeprecated|"
+ _deprecatedUtilsFunction('utils.encodeVideoFromFrames','utils.makeVideo')
+ return makeVideo(*args,**kw)
+
+def makeVideo(frameSpec,out,renameNotOverwrite=True,fps=24,bps=2400):
+ """Create a video from external image files using `mencoder <http://www.mplayerhq.hu>`__. Two-pass encoding using the default mencoder codec (mpeg4) is performed, running multi-threaded with number of threads equal to number of OpenMP threads allocated for Yade.
+
+ :param frameSpec: wildcard | sequence of filenames. If list or tuple, filenames to be encoded in given order; otherwise wildcard understood by mencoder's mf:// URI option (shell wildcards such as ``/tmp/snap-*.png`` or and printf-style pattern like ``/tmp/snap-%05d.png``)
+ :param str out: file to save video into
+ :param bool renameNotOverwrite: if True, existing same-named video file will have -*number* appended; will be overwritten otherwise.
+ :param int fps: Frames per second (``-mf fps=…``)
+ :param int bps: Bitrate (``-lavcopts vbitrate=…``)
"""
- import pygst,sys,gobject,os,tempfile,shutil,os.path
- pygst.require("0.10")
- import gst
+ import os,os.path,subprocess
if renameNotOverwrite and os.path.exists(out):
- i=0;
+ i=0
while(os.path.exists(out+"~%d"%i)): i+=1
os.rename(out,out+"~%d"%i); print "Output file `%s' already existed, old file renamed to `%s'"%(out,out+"~%d"%i)
- if frameSpec.__class__==list or frameSpec.__class__==tuple:
- # create a new common prefix, symlink given files to prefix-%05d.png in their order, adjust wildcard, go ahead.
- tmpDir=tempfile.mkdtemp()
- for no,frame in enumerate(frameSpec):
- os.symlink(os.path.abspath(frame),os.path.join(tmpDir,'%07d'%no))
- wildcard=os.path.join(tmpDir,'%07d')
- else:
- tmpDir=None; wildcard=frameSpec
- print "Encoding video from %s to %s"%(wildcard,out)
- pipeline=gst.parse_launch('multifilesrc location="%s" index=0 caps="image/png,framerate=\(fraction\)%d/1" ! pngdec ! ffmpegcolorspace ! theoraenc sharpness=2 quality=63 ! oggmux ! filesink location="%s"'%(wildcard,fps,out))
- bus=pipeline.get_bus()
- bus.add_signal_watch()
- mainloop=gobject.MainLoop();
- bus.connect("message::eos",lambda bus,msg: mainloop.quit())
- pipeline.set_state(gst.STATE_PLAYING)
- mainloop.run()
- pipeline.set_state(gst.STATE_NULL); pipeline.get_state()
- # remove symlinked frames, if any
- if tmpDir: shutil.rmtree(tmpDir)
+ if isinstance(frameSpec,list) or isinstance(frameSpec,tuple): frameSpec=','.join(frameSpec)
+ for passNo in (1,2):
+ cmd=['mencoder','mf://%s'%frameSpec,'-mf','fps=%d'%int(fps),'-ovc','lavc','-lavcopts','vbitrate=%d:vpass=%d:threads=%d:%s'%(int(bps),passNo,O.numThreads,'turbo' if passNo==1 else ''),'-o',('/dev/null' if passNo==1 else out)]
+ print 'Pass %d:'%passNo,' '.join(cmd)
+ ret=subprocess.call(cmd)
+ if ret!=0: raise RuntimeError("Error when running mencoder.")
def replaceCollider(colliderEngine):
Follow ups
