yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05933
[Branch ~yade-dev/yade/trunk] Rev 2499: 1. undefine pi in MindlinPhys
------------------------------------------------------------
revno: 2499
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Mon 2010-10-18 21:27:07 +0200
message:
1. undefine pi in MindlinPhys
2. Add isAspherical() flag to Body
3. Adjust kineticEnergy with PBC (not yet fully tested)
4. Make createInteraction work with PBC
5. Fix getIncidentVelocity with PBC and !avoidGranularRatcheting
6. Add some regression tests related to Cell (not yet fully done)
added:
py/tests/cell.py
modified:
.kdev4/yade.kdev4
core/Body.hpp
core/main/main.py.in
core/main/yade-batch.in
pkg/common/Dispatching.cpp
pkg/dem/Clump.cpp
pkg/dem/HertzMindlin.cpp
pkg/dem/L3Geom.hpp
pkg/dem/NewtonIntegrator.cpp
pkg/dem/NewtonIntegrator.hpp
pkg/dem/Shop.cpp
py/SConscript
py/remote.py
py/tests/__init__.py
py/tests/omega.py
py/tests/wrapper.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 '.kdev4/yade.kdev4'
--- .kdev4/yade.kdev4 2010-07-16 10:52:03 +0000
+++ .kdev4/yade.kdev4 2010-10-18 19:27:07 +0000
@@ -1,4 +1,5 @@
[Buildset]
+BuildItems=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x00\x01\x00\x00\x00\x08\x00y\x00a\x00d\x00e)
[MakeBuilder]
Make Binary=scons
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2010-09-23 09:17:40 +0000
+++ core/Body.hpp 2010-10-18 19:27:07 +0000
@@ -25,7 +25,7 @@
public:
typedef int id_t;
// bits for Body::flags
- enum { FLAG_DYNAMIC=1, FLAG_BOUNDED=2 }; /* add powers of 2 as needed */
+ enum { FLAG_DYNAMIC=1, FLAG_BOUNDED=2, FLAG_ASPHERICAL=4 }; /* add powers of 2 as needed */
//! symbolic constant for body that doesn't exist.
static const Body::id_t ID_NONE;
//! get Body pointer given its id.
@@ -48,6 +48,8 @@
void setDynamic(bool d){ if(d) flags|=FLAG_DYNAMIC; else flags&=~(FLAG_DYNAMIC); }
bool isBounded() const {return flags & FLAG_BOUNDED; }
void setBounded(bool d){ if(d) flags|=FLAG_BOUNDED; else flags&=~(FLAG_BOUNDED); }
+ bool isAspherical() const {return flags & FLAG_ASPHERICAL; }
+ void setAspherical(bool d){ if(d) flags|=FLAG_ASPHERICAL; else flags&=~(FLAG_ASPHERICAL); }
/*! Hook for clump to update position of members when user-forced reposition and redraw (through GUI) occurs.
* This is useful only in cases when engines that do that in every iteration are not active - i.e. when the simulation is paused.
@@ -81,6 +83,7 @@
.add_property("isDynamic",&Body::isDynamic,&Body::setDynamic,"Deprecated synonym for :yref:`Body::dynamic` |ydeprecated|")
.add_property("dynamic",&Body::isDynamic,&Body::setDynamic,"Whether this body will be moved by forces. (In c++, use ``Body::isDynamic``/``Body::setDynamic``) :ydefault:`true`")
.add_property("bounded",&Body::isBounded,&Body::setBounded,"Whether this body should have :yref:`Body.bound` created. Note that bodies without a :yref:`bound <Body.bound>` do not participate in collision detection. (In c++, use ``Body::isBounded``/``Body::setBounded``) :ydefault:`true`")
+ .add_property("aspherical",&Body::isAspherical,&Body::setAspherical,"Whether this body has different inertia along principal axes; :yref:`NewtonIntegrator` makes use of this flag to call rotation integration routine for aspherical bodies, which is more expensive. :ydefault:`false`")
.def_readwrite("mask",&Body::groupMask,"Shorthand for :yref:`Body::groupMask`")
.add_property("isStandalone",&Body::isStandalone,"True if this body is neither clump, nor clump member; false otherwise.")
.add_property("isClumpMember",&Body::isClumpMember,"True if this body is clump member, false otherwise.")
=== modified file 'core/main/main.py.in'
--- core/main/main.py.in 2010-10-15 16:49:26 +0000
+++ core/main/main.py.in 2010-10-18 19:27:07 +0000
@@ -165,6 +165,7 @@
# print 'worker started'
#else:
runScript(arg0)
+ if yade.runtime.opts.exitAfter: sys.exit(0)
# show python console
if 1:
from IPython.Shell import IPShellEmbed
=== modified file 'core/main/yade-batch.in'
--- core/main/yade-batch.in 2010-10-11 16:28:49 +0000
+++ core/main/yade-batch.in 2010-10-18 19:27:07 +0000
@@ -269,6 +269,7 @@
parser.add_option('-l','--lines',dest='lineList',help='Lines of TABLE to use, in the format 2,3-5,8,11-13 (default: all available lines in TABLE)',metavar='LIST')
parser.add_option('--nice',dest='nice',type='int',help='Nice value of spawned jobs (default: 10)',default=10)
parser.add_option('--executable',dest='executable',help='Name of the program to run (default: %s)'%executable,default=executable,metavar='FILE')
+parser.add_option('--rebuild',dest='rebuild',help='Run executable(s) with --rebuild prior to running any jobs.',default=False,action='store_true')
parser.add_option('--debug',dest='debug',action='store_true',help='Run the executable with --debug.',default=False)
parser.add_option('--gnuplot',dest='gnuplotOut',help='Gnuplot file where gnuplot from all jobs should be put together',default=None,metavar='FILE')
parser.add_option('--dry-run',action='store_true',dest='dryRun',help='Do not actually run (useful for getting gnuplot only, for instance)',default=False)
@@ -323,6 +324,7 @@
print "Will use lines ",', '.join([str(i)+' (%s)'%params[i]['description'] for i in useLines])+'.'
jobs=[]
+executables=set()
for i,l in enumerate(useLines):
logFile=logFormat.replace('%',str(l))
logFile=logFile.replace('@',params[l]['description'])
@@ -340,6 +342,7 @@
nSlots=maxJobs
else:
logging.warning('WARNING: job #%d will use %d slots but only %d are available'%(i,nSlots,maxJobs))
+ executables.add(executable)
# 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%s --threads=%d %s -x <a href="jobs/%d/script">%s</a>'%(executable,' --debug' if opts.debug else '',int(nSlots),'--nice=%s'%nice if nice!=None else '',i,simul)
@@ -350,11 +353,21 @@
print "Master process pid",os.getpid()
+if opts.rebuild:
+ print "Rebuilding all active executables, since --rebuild was specified"
+ for e in executables:
+ import subprocess
+ if subprocess.call([e,'--rebuild','-x']+(['--debug'] if opts.debug else [])):
+ raise RuntimeError('Error rebuilding %s (--rebuild).'%e)
+ print "Rebuilding done."
+
+
print "Job summary:"
for job in jobs:
print ' #%d (%s%s):'%(job.num,job.id,'' if job.nSlots==1 else '/%d'%job.nSlots),job.command
sys.stdout.flush()
+
httpLastServe=0
runHttpStatsServer()
if opts.plotAlwaysUpdateTime>0:
=== modified file 'pkg/common/Dispatching.cpp'
--- pkg/common/Dispatching.cpp 2010-10-13 16:23:08 +0000
+++ pkg/common/Dispatching.cpp 2010-10-18 19:27:07 +0000
@@ -60,7 +60,12 @@
shared_ptr<Interaction> IGeomDispatcher::explicitAction(const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool force){
scene=Omega::instance().getScene().get(); // to make sure if called from outside of the loop
- if(scene->isPeriodic) throw logic_error("IGeomDispatcher::explicitAction does not support periodic boundary conditions (O.periodic==True)");
+ Vector3i cellDist=Vector3i::Zero();
+ if(scene->isPeriodic) {
+ //throw logic_error("IGeomDispatcher::explicitAction does not support periodic boundary conditions (O.periodic==True)");
+ for(int i=0; i<3; i++) cellDist[i]=-(int)((b2->state->pos[i]-b1->state->pos[i])/scene->cell->getSize()[i]+.5);
+ }
+ Vector3r shift2=scene->cell->Hsize*cellDist.cast<Real>();
updateScenePtr();
if(force){
#ifdef YADE_DEVIRT_FUNCTORS
@@ -68,6 +73,7 @@
#endif
assert(b1->shape && b2->shape);
shared_ptr<Interaction> I(new Interaction(b1->getId(),b2->getId()));
+ I->cellDist=cellDist;
// FIXME: this code is more or less duplicated from InteractionLoop :-(
bool swap=false;
I->functorCache.geom=getFunctor2D(b1->shape,b2->shape,swap);
@@ -75,12 +81,13 @@
if(swap){I->swapOrder();}
const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
- bool succ=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,/*shift2*/Vector3r::Zero(),/*force*/true,I);
+ bool succ=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/true,I);
if(!succ) throw logic_error("Functor "+I->functorCache.geom->getClassName()+"::go returned false, even if asked to force IGeom creation. Please report bug.");
return I;
} else {
shared_ptr<Interaction> I(new Interaction(b1->getId(),b2->getId()));
- b1->shape && b2->shape && operator()( b1->shape , b2->shape , *b1->state , *b2->state , Vector3r::Zero(), /*force*/ false, I);
+ I->cellDist=cellDist;
+ b1->shape && b2->shape && operator()(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/ false,I);
return I;
}
}
=== modified file 'pkg/dem/Clump.cpp'
--- pkg/dem/Clump.cpp 2010-10-17 12:55:29 +0000
+++ pkg/dem/Clump.cpp 2010-10-18 19:27:07 +0000
@@ -248,6 +248,8 @@
// TODO: these might be calculated from members... but complicated... - someone needs that?!
state->vel=state->angVel=Vector3r::Zero();
+ if(state->inertia[0]!=state->inertia[1] || state->inertia[0]!=state->inertia[2]) this->setAspherical(true);
+
// update subBodySe3s; subtract clump orientation (=apply its inverse first) to subBody's orientation
// Conjugate is equivalent to Inverse for normalized quaternions
for(memberMap::iterator I=members.begin(); I!=members.end(); I++){
=== modified file 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp 2010-10-16 18:31:17 +0000
+++ pkg/dem/HertzMindlin.cpp 2010-10-18 19:27:07 +0000
@@ -6,8 +6,6 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-#define pi 3.14159265
-
YADE_PLUGIN(
(MindlinPhys)
(Ip2_FrictMat_FrictMat_MindlinPhys)
@@ -66,7 +64,7 @@
Real Kso = 2*sqrt(4*R)*G/(2-V); // coefficient for shear stiffness
Real frictionAngle = std::min(fa,fb);
- Real Adhesion = 4.*pi*R*gamma; // calculate adhesion force as predicted by DMT theory
+ Real Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory
/* pass values calculated from above to MindlinPhys */
mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle);
=== modified file 'pkg/dem/L3Geom.hpp'
--- pkg/dem/L3Geom.hpp 2010-10-16 18:31:17 +0000
+++ pkg/dem/L3Geom.hpp 2010-10-18 19:27:07 +0000
@@ -1,4 +1,5 @@
+#pragma once
#include<yade/core/IGeom.hpp>
#include<yade/core/IPhys.hpp>
#include<yade/core/Shape.hpp>
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2010-10-18 19:27:07 +0000
@@ -121,7 +121,7 @@
if(isnan(m[0])||isnan(m[1])||isnan(m[2])) throw runtime_error(("NewtonIntegrator: NaN torque acting on #"+lexical_cast<string>(id)+".").c_str());
#endif
// exactAsphericalRot is disabled or the body is spherical
- if (!exactAsphericalRot || (state->inertia[0]==state->inertia[1] && state->inertia[1]==state->inertia[2])){
+ if (!exactAsphericalRot || !b->isAspherical()){
state->angAccel=m.cwise()/state->inertia;
cundallDamp(dt,m,state->angVel,state->angAccel);
leapfrogSphericalRotate(scene,state,id,dt);
@@ -139,7 +139,7 @@
Vector3r M(m);
// sum force on clump memebrs
// exactAsphericalRot enabled and clump is aspherical
- if (exactAsphericalRot && ((state->inertia[0]!=state->inertia[1] || state->inertia[1]!=state->inertia[2]))){
+ if (exactAsphericalRot && b->isAspherical()){
FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
const Body::id_t& memberId=mm.first;
const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
@@ -197,7 +197,7 @@
if (homotheticCellResize>1) state->vel+=prevVelGrad*state->vel*dt;
//In all cases, reflect macroscopic (periodic cell) acceleration in the velocity. This is the dominant term in the update in most cases
- Vector3r dVel=(scene->cell->velGrad-prevVelGrad)*/*scene->cell->wrapShearedPt(*/state->pos/*)*/;
+ Vector3r dVel=(scene->cell->velGrad-prevVelGrad)*state->pos;
state->vel+=dVel;
}
state->vel+=dt*state->accel;
=== modified file 'pkg/dem/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NewtonIntegrator.hpp 2010-10-18 19:27:07 +0000
@@ -49,7 +49,7 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR(NewtonIntegrator,GlobalEngine,"Engine integrating newtonian motion equations.",
((Real,damping,0.2,,"damping coefficient for Cundall's non viscous damping (see [Chareyre2005]_) [-]"))
((Real,maxVelocitySq,NaN,,"store square of max. velocity, for informative purposes; computed again at every step. |yupdate|"))
- ((bool,exactAsphericalRot,true,,"Enable more exact body rotation integrator for aspherical bodies *only*, using formulation from [Allen1989]_, pg. 89."))
+ ((bool,exactAsphericalRot,true,,"Enable more exact body rotation integrator for :yref:`aspherical bodies<Body.aspherical>` *only*, using formulation from [Allen1989]_, pg. 89."))
((unsigned short,homotheticCellResize,0,,"Enable velocity updates insuring that all bodies move homogeneously while the periodic cell deforms. 0: No update, 1: Reflect on each body the changes in macroscopic velocity field :yref:`Cell::velGrad`, using Δv_i=Δ(grad(v_macro))*x_i. 2: Emulate the Cundall-style equation Δx_i=(grad(v_macro))*x_i, by adding a convective term in the velocity update."))
((Matrix3r,prevVelGrad,Matrix3r::Zero(),,"Store previous velocity gradient (:yref:`Cell::velGrad`) to track acceleration. |yupdate|"))
((vector<shared_ptr<BodyCallback> >,callbacks,,,"List (std::vector in c++) of :yref:`BodyCallbacks<BodyCallback>` which will be called for each body as it is being processed."))
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/Shop.cpp 2010-10-18 19:27:07 +0000
@@ -180,15 +180,62 @@
return (useMaxForce?maxF:meanF)/(sumF);
}
-Real Shop::kineticEnergy(Scene* _rb, Body::id_t* maxId){
- Scene* rb=_rb ? _rb : Omega::instance().getScene().get();
+Real Shop::kineticEnergy(Scene* _scene, Body::id_t* maxId){
+ Scene* scene=_scene ? _scene : Omega::instance().getScene().get();
Real ret=0.;
Real maxE=0; if(maxId) *maxId=Body::ID_NONE;
- FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ const bool isPeriodic=scene->isPeriodic;
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
if(!b || !b->isDynamic()) continue;
+ const State* state(b->state.get());
// ½(mv²+ωIω)
- Real E=.5*(b->state->mass*b->state->vel.squaredNorm()+b->state->angVel.dot(b->state->inertia.cwise()*b->state->angVel));
- if(maxId) { if(E>maxE) *maxId=b->getId(); maxE=max(E,maxE); }
+ Vector3r vel;
+ if(!isPeriodic) vel=b->state->vel;
+ else{
+ /* Kinetic energy is defined with regards to absolute space (or off by a constant for system in steady motion).
+ However, in the periodic case with non-zero velocity gradient, there is no such inertial system,
+ as velocity depends on absolute coordinate in space (NewtonIntegrator::homotheticCellResize).
+ There are several inconsistent options (and no consistent one):
+
+ 1. Subtract the contribution of homotheticCellResize; such kinetic energy however materializes
+ when there is an interaction between particles, and it would come apparently from nowhere.
+
+ 2. Compute the kinetic energy for particles translated inside the periodic cell. That way,
+ we avoid the dependency on absolute space position, while still keeping the contribution
+ of homothetic resize. This comes at the cost of discontinuity which is apparent for 1 particle.
+
+ Suppose that particle A is inside the base cell (i.e. between the origin and O.cell.size) and that
+ it has some velocity away from the origin. Applying velocity gradient, velocity of the particle will
+ be augmented by the amount corrseponding to its distance from the origin. At some moment, A will
+ leave the base cell. Kinetic energy computed from coordinates wrapped inside the base cell will suddenly
+ drop, since the wrapped coordinate will jump from O.cell.size[i] to 0.
+
+ With increasing number of particles, this effect will become less and less apparent,
+ eventually approaching zero. Moreoved, for dense packings, two directions of motion accross the cell
+ boundary will compensate each other.
+
+ 3. Propose some special form of potential energy related to the periodic cell. Not sure if that is doable,
+ since the amount of energy depends also on what interaction (with which particles) will be established
+ in the future.
+
+ */
+ Vector3i period; Vector3r basePos=scene->cell->wrapShearedPt(state->pos,period);
+ vel=state->vel-scene->cell->velGrad*(state->pos-basePos);
+ // TODO: move NewtonIntegrator.homotheticCellResize to Cell.homotheticDeformation so that
+ // we have access to how is the velocity adjusted
+ // in addition, create function in Cell that will compute velocity compensations etc for us
+ // since it might be used in more places than just here (code audit needed)
+ }
+ Real E=.5*(state->mass*vel.squaredNorm());
+ if(b->isAspherical()){
+ Matrix3r T(state->ori);
+ // the tensorial expression http://en.wikipedia.org/wiki/Moment_of_inertia#Moment_of_inertia_tensor
+ // inertia tensor rotation from http://www.kwon3d.com/theory/moi/triten.html
+ Matrix3r mI; mI.diagonal()=state->inertia;
+ E+=.5*state->angVel.transpose().dot((T.transpose()*mI*T)*state->angVel);
+ }
+ else { E+=state->angVel.dot(state->inertia.cwise()*state->angVel);}
+ if(maxId && E>maxE) { *maxId=b->getId(); maxE=E; }
ret+=E;
}
return ret;
=== modified file 'py/SConscript'
--- py/SConscript 2010-10-14 12:36:04 +0000
+++ py/SConscript 2010-10-18 19:27:07 +0000
@@ -45,7 +45,8 @@
env.Install('$LIBDIR/py/yade/tests',[
env.File('__init__.py','tests'),
env.File('wrapper.py','tests'),
- env.File('omega.py','tests')
+ env.File('omega.py','tests'),
+ env.File('cell.py','tests'),
])
# 3rd party modules:
=== modified file 'py/remote.py'
--- py/remote.py 2010-09-10 11:33:21 +0000
+++ py/remote.py 2010-10-18 19:27:07 +0000
@@ -26,16 +26,22 @@
sys.stdout.flush(); sys.stderr.flush()
return ret
def plot(self):
- from yade import plot
- if len(plot.plots)==0: return None
- fig=plot.plot(subPlots=True,noShow=True)
- img=O.tmpFilename()+'.'+plotImgFormat
- sqrtFigs=math.sqrt(len(plot.plots))
- fig.set_size_inches(5*sqrtFigs,7*sqrtFigs)
- fig.savefig(img)
- f=open(img,'rb'); data=f.read(); f.close(); os.remove(img)
- #print 'returning '+plotImgFormat
- return xmlrpclib.Binary(data)
+ try
+ from yade import plot
+ if len(plot.plots)==0: return None
+ fig=plot.plot(subPlots=True,noShow=True)
+ img=O.tmpFilename()+'.'+plotImgFormat
+ sqrtFigs=math.sqrt(len(plot.plots))
+ fig.set_size_inches(5*sqrtFigs,7*sqrtFigs)
+ fig.savefig(img)
+ f=open(img,'rb'); data=f.read(); f.close(); os.remove(img)
+ #print 'returning '+plotImgFormat
+ return xmlrpclib.Binary(data)
+ except:
+ print 'Error updating plots:'
+ import traceback
+ traceback.print_exc()
+ return None
class PythonConsoleSocketEmulator(SocketServer.BaseRequestHandler):
=== modified file 'py/tests/__init__.py'
--- py/tests/__init__.py 2010-08-24 12:54:14 +0000
+++ py/tests/__init__.py 2010-10-18 19:27:07 +0000
@@ -4,7 +4,7 @@
import unittest,inspect
# add any new test suites to the list here, so that they are picked up by testAll
-allTests=['wrapper','omega']
+allTests=['wrapper','omega','cell']
# all yade modules (ugly...)
import yade.eudoxos,yade.export,yade.linterpolation,yade.log,yade.pack,yade.plot,yade.post2d,yade.timing,yade.utils,yade.ymport
=== added file 'py/tests/cell.py'
--- py/tests/cell.py 1970-01-01 00:00:00 +0000
+++ py/tests/cell.py 2010-10-18 19:27:07 +0000
@@ -0,0 +1,57 @@
+# encoding: utf-8
+# 2010 © Václav Šmilauer <eudoxos@xxxxxxxx>
+
+'''
+Various computations affected by the periodic boundary conditions.
+'''
+
+import unittest
+import random
+from yade.wrapper import *
+from miniEigen import *
+from yade._customConverters import *
+from yade import utils
+from yade import *
+
+
+class TestCell(unittest.TestCase):
+ def setUp(self):
+ O.reset(); O.periodic=True;
+ O.cell.refSize=Vector3(2.5,2.5,3)
+ self.cellDist=Vector3i(0,0,10) # how many cells away we go
+ self.relDist=Vector3(0,.999,0) # rel position of the 2nd ball within the cell
+ self.initVel=Vector3(0,0,5)
+ O.bodies.append(utils.sphere((1,1,1),.5))
+ self.initPos=Vector3([O.bodies[0].state.pos[i]+self.relDist[i]+self.cellDist[i]*O.cell.refSize[i] for i in (0,1,2)])
+ O.bodies.append(utils.sphere(self.initPos,.5))
+ print O.bodies[1].state.pos
+ O.bodies[1].state.vel=self.initVel
+ O.dt=1e-5
+ O.engines=[NewtonIntegrator(homotheticCellResize=2)]
+ O.cell.velGrad=Matrix3(0,0,0, 0,0,0, 0,0,-1)
+ def testHomotheticResize(self):
+ "Homothetic cell resize adjusts particle velocity"
+ O.step()
+ s1=O.bodies[1].state
+ self.assertAlmostEqual(s1.vel[2],self.initVel[2]+self.initPos[2]*O.cell.velGrad[2,2])
+ def testScGeomIncidentVelocity(self):
+ O.step()
+ O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
+ i=utils.createInteraction(0,1)
+ if 0:
+ print i.geom.incidentVel(i,avoidGranularRatcheting=True)
+ print i.geom.incidentVel(i,avoidGranularRatcheting=False)
+ print i.geom.penetrationDepth
+ print i.geom.dict()
+ print '=========='
+ print 'v1',O.bodies[1].state.vel
+ print 'shiftVel',O.cell.velGrad*O.cell.Hsize*(0,0,-10)
+
+ def testKineticEnergy(self):
+ O.step() # updates velocity with homotheticCellResize
+ # ½(mv²+ωIω)
+ # #0 is still, no need to add it; #1 has zero angular velocity
+ # we must take self.initVel since O.bodies[1].state.vel now contains the homothetic resize which utils.kineticEnergy is supposed to compensate back
+ Ek=.5*O.bodies[1].state.mass*self.initVel.squaredNorm()
+
+
=== modified file 'py/tests/omega.py'
--- py/tests/omega.py 2010-09-18 14:10:30 +0000
+++ py/tests/omega.py 2010-10-18 19:27:07 +0000
@@ -14,7 +14,7 @@
## TODO tests
class TestInteractions(unittest.TestCase): pass
-class TestBex(unittest.TestCase): pass
+class TestForce(unittest.TestCase): pass
class TestEngines(unittest.TestCase): pass
class TestIO(unittest.TestCase): pass
class TestTags(unittest.TestCase): pass
=== modified file 'py/tests/wrapper.py'
--- py/tests/wrapper.py 2010-10-15 16:49:26 +0000
+++ py/tests/wrapper.py 2010-10-18 19:27:07 +0000
@@ -22,29 +22,34 @@
def setUp(self):
pass # no setup needed for tests here
def testClassCtors(self):
+ "Correct types are instantiated"
# correct instances created with Foo() syntax
for r in allClasses:
obj=eval(r)();
self.assert_(obj.__class__.__name__==r,'Failed for '+r)
def testRootDerivedCtors_attrs_few(self):
+ "Class ctor's attributes"
# attributes passed when using the Foo(attr1=value1,attr2=value2) syntax
gm=Shape(wire=True); self.assert_(gm.wire==True)
def testDispatcherCtor(self):
+ "Dispatcher ctors with functors"
# dispatchers take list of their functors in the ctor
# same functors are collapsed in one
cld1=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_CundallStrack(),Law2_Dem3DofGeom_FrictPhys_CundallStrack()]); self.assert_(len(cld1.functors)==1)
# two different make two different, right?
cld2=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_CundallStrack(),Law2_Dem3DofGeom_CpmPhys_Cpm()]); self.assert_(len(cld2.functors)==2)
def testInteractionLoopCtor(self):
+ "InteractionLoop special ctor"
# InteractionLoop takes 3 lists
id=InteractionLoop([Ig2_Facet_Sphere_Dem3DofGeom(),Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],)
self.assert_(len(id.geomDispatcher.functors)==2)
- self.assert_(id.geomDispatcher.__class__==IGeomDispatcher)
- self.assert_(id.physDispatcher.functors[0].__class__==Ip2_FrictMat_FrictMat_FrictPhys)
- self.assert_(id.lawDispatcher.functors[0].__class__==Law2_Dem3DofGeom_FrictPhys_CundallStrack)
+ self.assert_(id.geomDispatcher.__class__==IGeomDispatcher().__class__)
+ self.assert_(id.physDispatcher.functors[0].__class__==Ip2_FrictMat_FrictMat_FrictPhys().__class__)
+ self.assert_(id.lawDispatcher.functors[0].__class__==Law2_Dem3DofGeom_FrictPhys_CundallStrack().__class__)
def testParallelEngineCtor(self):
+ "ParallelEngine special ctor"
pe=ParallelEngine([InsertionSortCollider(),[BoundDispatcher(),ForceResetter()]])
- self.assert_(pe.slaves[0].__class__==InsertionSortCollider)
+ self.assert_(pe.slaves[0].__class__==InsertionSortCollider().__class__)
self.assert_(len(pe.slaves[1])==2)
pe.slaves=[]
self.assert_(len(pe.slaves)==0)
@@ -65,11 +70,13 @@
self.assertEqual(len(v1),len(v2));
for i in range(len(v1)): self.assertAlmostEqual(v1[i],v2[i],msg='Component '+str(i)+' of '+str(v1)+' and '+str(v2))
def testVector2(self):
+ "Vector2 operations"
v=Vector2(1,2); v2=Vector2(3,4)
self.assert_(v+v2==Vector2(4,6))
self.assert_(Vector2().UnitX.dot(Vector2().UnitY)==0)
self.assert_(Vector2().Zero.norm()==0)
def testVector3(self):
+ "Vector3 operations"
v=Vector3(3,4,5); v2=Vector3(3,4,5)
self.assert_(v[0]==3 and v[1]==4 and v[2]==5)
self.assert_(v.squaredNorm()==50)
@@ -81,6 +88,7 @@
self.assert_(x.dot(y)==0)
self.assert_(x.cross(y)==z)
def testQuaternion(self):
+ "Quaternion operations"
# construction
q1=Quaternion((0,0,1),pi/2)
q2=Quaternion(Vector3(0,0,1),pi/2)
@@ -92,6 +100,7 @@
self.assertSeqAlmostEqual(q1.toAxisAngle()[0],(0,0,1))
self.assertAlmostEqual(q1.toAxisAngle()[1],pi/2)
def testMatrix3(self):
+ "Matrix3 operations"
#construction
m1=Matrix3(1,0,0,0,1,0,0,0,1)
# comparison
=== modified file 'py/utils.py'
--- py/utils.py 2010-09-18 14:10:30 +0000
+++ py/utils.py 2010-10-18 19:27:07 +0000
@@ -196,6 +196,7 @@
_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material)
b.state.pos=b.state.refPos=center
b.dynamic=dynamic
+ b.aspherical=False
b.mask=mask
return b
@@ -213,6 +214,7 @@
b.state.pos=b.state.refPos=center
b.dynamic=dynamic
b.mask=mask
+ b.aspherical=True
return b
@@ -261,6 +263,7 @@
else: pos2=position
b.state.pos=b.state.refPos=pos2
b.dynamic=False
+ b.aspherical=True # never used, since the wall is not dynamic
b.mask=mask
return b
@@ -281,6 +284,7 @@
_commonBodySetup(b,0,Vector3(0,0,0),material,noBound=noBound)
b.state.pos=b.state.refPos=center
b.dynamic=dynamic
+ b.aspherical=True
b.mask=mask
return b
