yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05968
[Branch ~yade-dev/yade/trunk] Rev 2511: 1. Deprecate NewtonIntegrator::homotheticCellResize in favor of Cell::homoDeform (compatibility i...
------------------------------------------------------------
revno: 2511
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Tue 2010-10-26 15:41:30 +0200
message:
1. Deprecate NewtonIntegrator::homotheticCellResize in favor of Cell::homoDeform (compatibility interface with warning)
2. Add the possibility of homothetically changing positions rather than velocities; avoid the meanfield/fluctuation velocity jazz at the expense of vel not being time derivative of pos
3. Add PBC tests for the previous variant
4. Add defThreads arg to scons, which, if specified, gives default number of threads for simulations, if not overridden with -j
5. Fix a few example scripts, more work ahead come (please help!)
6. Add regression tests for saving/loading yade objects; on maverick, this makes the unregistered class to re-surface (!!)
BTW Bruno, can you run --test and fix chained cylinders attributes you changed?
modified:
SConstruct
core/Cell.cpp
core/Cell.hpp
core/main/main.py.in
debian/control-template
examples/baraban/baraban.py
examples/bulldozer/bulldozer.py
pkg/common/Callbacks.cpp
pkg/common/Callbacks.hpp
pkg/common/Dispatching.cpp
pkg/dem/NewtonIntegrator.cpp
pkg/dem/NewtonIntegrator.hpp
pkg/dem/ScGeom.cpp
pkg/dem/Shop.cpp
py/system.py
py/tests/omega.py
py/tests/pbc.py
py/ymport.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 'SConstruct'
--- SConstruct 2010-10-14 08:18:42 +0000
+++ SConstruct 2010-10-26 13:41:30 +0000
@@ -134,6 +134,7 @@
('march','Architecture to use with -march=... when optimizing','native',None,None),
BoolVariable('mono','[experimental] Build only one shared library and make all other files (python objects, for instance) only be symlinks.',0),
('execCheck','Name of the main script that should be installed; if the current one differs, an erro is raised (do not use directly, only intended for --rebuild',None),
+ ('defThreads','Default number of OpenMP threads to run with, unless overridden with -j. Keep unset (negative) to default to using all cores.',-1),
#('SHLINK','Linker for shared objects','g++'),
('SHCCFLAGS','Additional compiler flags for linking (for plugins).',None,None,Split),
BoolVariable('QUAD_PRECISION','typedef Real as long double (=quad)',0),
=== modified file 'core/Cell.cpp'
--- core/Cell.cpp 2010-05-26 08:18:36 +0000
+++ core/Cell.cpp 2010-10-26 13:41:30 +0000
@@ -35,6 +35,8 @@
_hasShear=(trsf(0,1)!=0 || trsf(0,2)!=0 || trsf(1,0)!=0 || trsf(1,2)!=0 || trsf(2,0)!=0 || trsf(2,1)!=0);
// OpenGL shear matrix (used frequently)
fillGlShearTrsfMatrix(_glShearTrsfMatrix);
+
+ if(!(homoDeform==HOMO_NONE || homoDeform==HOMO_POS || homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND)) throw std::invalid_argument("Cell.homoDeform must be in {0,1,2,3}.");
}
void Cell::fillGlShearTrsfMatrix(double m[16]){
=== modified file 'core/Cell.hpp'
--- core/Cell.hpp 2010-09-23 09:17:40 +0000
+++ core/Cell.hpp 2010-10-26 13:41:30 +0000
@@ -85,6 +85,12 @@
Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz;
}
+ // relative position and velocity for interaction accross multiple cells
+ Vector3r intrShiftPos(const Vector3i& cellDist) const { return Hsize*cellDist.cast<Real>(); }
+ Vector3r intrShiftVel(const Vector3i& cellDist) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return velGrad*Hsize*cellDist.cast<Real>(); return Vector3r::Zero(); }
+ // return body velocity while taking away mean field velocity (coming from velGrad) if the mean field velocity is applied on velocity
+ Vector3r bodyFluctuationVel(const Vector3r& pos, const Vector3r& vel) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return (vel-velGrad*pos); return vel; }
+
Vector3r getRefSize() const { return refSize; }
void setRefSize(const Vector3r& s){ refSize=s; integrateAndUpdate(0); }
Matrix3r getTrsf() const { return trsf; }
@@ -97,12 +103,16 @@
// to resolve overloads
Vector3r wrapShearedPt_py(const Vector3r& pt) const { return wrapShearedPt(pt);}
Vector3r wrapPt_py(const Vector3r& pt) const { return wrapPt(pt);}
+
+ enum { HOMO_NONE=0, HOMO_POS=1, HOMO_VEL=2, HOMO_VEL_2ND=3 };
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
((Vector3r,refSize,Vector3r(1,1,1),Attr::triggerPostLoad,"Reference size of the cell."))
((Matrix3r,trsf,Matrix3r::Identity(),Attr::triggerPostLoad,"Current transformation matrix of the cell."))
((Matrix3r,velGrad,Matrix3r::Zero(),,"Velocity gradient of the transformation; used in NewtonIntegrator."))
- ((Matrix3r,Hsize,Matrix3r::Zero(),Attr::readonly,"The current cell size (one column per box edge), computed from *refSize* and *trsf* |yupdate|")),
+ ((Matrix3r,prevVelGrad,Matrix3r::Zero(),Attr::readonly,"Velocity gradient in the previous step."))
+ ((Matrix3r,Hsize,Matrix3r::Zero(),Attr::readonly,"The current cell size (one column per box edge), computed from *refSize* and *trsf* |yupdate|"))
+ ((int,homoDeform,3,Attr::triggerPostLoad,"Deform (:yref:`velGrad<Cell.velGrad>`) the cell homothetically, by adjusting positions or velocities of particles. The values have the following meaning: 0: no homothetic deformation, 1: set absolute particle positions directly (when ``velGrad`` is non-zero), but without changing their velocity, 2: adjust particle velocity (only when ``velGrad`` changed) with Δv_i=Δ ∇v x_i. 3: as 2, but include a 2nd order term in addition -- the derivative of 1 (convective term in the velocity update).")),
/*ctor*/ integrateAndUpdate(0),
/*py*/
.def_readonly("size",&Cell::getSize_copy,"Current size of the cell, i.e. lengths of 3 cell lateral vectors after applying current trsf. Update automatically at every step.")
=== modified file 'core/main/main.py.in'
--- core/main/main.py.in 2010-10-18 19:27:07 +0000
+++ core/main/main.py.in 2010-10-26 13:41:30 +0000
@@ -98,6 +98,8 @@
# this must be done before loading yade libs ("import yade" below)
# has no effeect after libgomp initializes
if opts.threads: os.environ['OMP_NUM_THREADS']=str(opts.threads)
+elif int('$defThreads')>0:
+ os.environ['OMP_NUM_THREADS']='$defThreads'
sys.stderr.write('Welcome to Yade '+version+'%s\n'%(' (debug build)' if opts.debug else ''))
=== modified file 'debian/control-template'
--- debian/control-template 2010-10-11 16:28:49 +0000
+++ debian/control-template 2010-10-26 13:41:30 +0000
@@ -7,7 +7,7 @@
Package: yade@_VERSION@
Architecture: any
-Depends: ${shlibs:Depends}, ${misc:Depends}, python-numpy, ipython, python-matplotlib, python-tk, python-qt4, python-xlib, mencoder
+Depends: ${shlibs:Depends}, ${misc:Depends}, python-numpy, ipython, python-matplotlib, python-tk, python-qt4, python-xlib, mencoder, gnuplot
Description: Platform for dynamical modeling.
Yet Another Dynamic Engine. etc.
.
=== modified file 'examples/baraban/baraban.py'
--- examples/baraban/baraban.py 2010-09-27 17:47:59 +0000
+++ examples/baraban/baraban.py 2010-10-26 13:41:30 +0000
@@ -10,7 +10,7 @@
es = 0.3
## Import wall's geometry
-params=utils.getViscoelasticFromSpheresInteraction(10e3,tc,en,es)
+params=utils.getViscoelasticFromSpheresInteraction(tc,en,es)
facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) # **params sets kn, cn, ks, cs
sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,**params))
from yade import ymport
@@ -26,7 +26,7 @@
y = (j*2 - nbSpheres[1])*sphereRadius*1.1
z = (k*2 - nbSpheres[2])*sphereRadius*1.1
s=utils.sphere([x,y,z],sphereRadius,material=sphereMat)
- p=utils.getViscoelasticFromSpheresInteraction(s.state.mass,tc,en,es)
+ p=utils.getViscoelasticFromSpheresInteraction(tc,en,es)
s.mat.kn,s.mat.cn,s.mat.ks,s.mat.cs=p['kn'],p['cn'],p['ks'],p['cs']
O.bodies.append(s)
@@ -55,7 +55,7 @@
## Saving results
#VTKRecorder(virtPeriod=0.04,fileName='/tmp/stlimp-',recorders=['spheres','facets']),
## Apply kinematics to walls
- RotationEngine(subscribedBodies=fctIds,rotationAxis=[0,0,1],rotateAroundZero=True,angularVelocity=0.5)
+ RotationEngine(ids=fctIds,rotationAxis=[0,0,1],rotateAroundZero=True,angularVelocity=0.5)
]
from yade import qt
=== modified file 'examples/bulldozer/bulldozer.py'
--- examples/bulldozer/bulldozer.py 2010-09-27 17:47:59 +0000
+++ examples/bulldozer/bulldozer.py 2010-10-26 13:41:30 +0000
@@ -29,6 +29,7 @@
KnifeP=[Knife,[p+Vector3(0,lengthKnife,0) for p in Knife]]
KnifePoly=pack.sweptPolylines2gtsSurface(KnifeP,threshold=1e-4)
+KnifeIDs=[]
KnifeIDs=O.bodies.append(pack.gtsSurface2Facets(KnifePoly,color=(1,0,0),wire=False))
@@ -49,14 +50,17 @@
colorsph1.normalize();
colorsph2.normalize();
colorSph=colorsph1
+#O.bodies.append(utils.sphere([0,0,0],1))
for xyz in itertools.product(arange(0,numBoxes[0]),arange(0,numBoxes[1]),arange(0,numBoxes[2])):
- ids_spheres=O.bodies.appendClumped(pack.regularHexa(pack.inEllipsoid((xyz[0]*(sizeBox+gapBetweenBoxes),xyz[1]*(sizeBox+gapBetweenBoxes)+sizeBox*0.5,xyz[2]*(sizeBox+gapBetweenBoxes)-radiusKnife+sizeBox*0.6),(sizeBox/2,sizeBox/2,sizeBox/2)),radius=radiusSph,gap=0,color=colorSph))
+ continue
+ #ids_spheres=O.bodies.appendClumped(pack.regularHexa(pack.inEllipsoid((xyz[0]*(sizeBox+gapBetweenBoxes),xyz[1]*(sizeBox+gapBetweenBoxes)+sizeBox*0.5,xyz[2]*(sizeBox+gapBetweenBoxes)-radiusKnife+sizeBox*0.6),(sizeBox/2,sizeBox/2,sizeBox/2)),radius=radiusSph,gap=0,color=colorSph))
if (colorSph==colorsph1):
colorSph=colorsph2
else:
colorSph=colorsph1
+from yade import qt
O.dt=2*utils.PWaveTimeStep() # We do not need now a high accuracy
O.engines=[
@@ -68,16 +72,17 @@
[Law2_Dem3DofGeom_FrictPhys_Basic()],
),
GravityEngine(gravity=(0,0,-9.8)),
- TranslationEngine(translationAxis=[1,0,0],velocity=5,subscribedBodies=KnifeIDs), # Buldozer motion
+ TranslationEngine(translationAxis=[1,0,0],velocity=5,ids=KnifeIDs), # Buldozer motion
NewtonIntegrator(damping=.3),
- SnapshotEngine(iterPeriod=100,fileBase='/tmp/bulldozer-',viewNo=0,label='snapshooter'),
+ #qt.SnapshotEngine(iterPeriod=100,fileBase='/tmp/bulldozer-',label='snapshooter'),
]
+O.engines=[]
+print len(O.bodies)
O.saveTmp()
-from yade import qt
qt.Controller()
qt.View()
r=qt.Renderer()
r.lightPos=Vector3(0,0,50)
O.run(2000); O.wait()
-utils.encodeVideoFromFrames(snapshooter.savedSnapshots,out='/tmp/bulldozer.ogg',fps=2)
+#utils.encodeVideoFromFrames(snapshooter.savedSnapshots,out='/tmp/bulldozer.ogg',fps=2)
=== modified file 'pkg/common/Callbacks.cpp'
--- pkg/common/Callbacks.cpp 2010-10-13 16:23:08 +0000
+++ pkg/common/Callbacks.cpp 2010-10-26 13:41:30 +0000
@@ -1,3 +1,4 @@
#include<yade/pkg-common/Callbacks.hpp>
-
+BodyCallback::~BodyCallback(){};
+IntrCallback::~IntrCallback(){};
YADE_PLUGIN((IntrCallback)(BodyCallback));
=== modified file 'pkg/common/Callbacks.hpp'
--- pkg/common/Callbacks.hpp 2010-10-13 16:23:08 +0000
+++ pkg/common/Callbacks.hpp 2010-10-26 13:41:30 +0000
@@ -10,6 +10,7 @@
class IntrCallback: public Serializable{
public:
+ virtual ~IntrCallback(); // vtable
typedef void(*FuncPtr)(IntrCallback*,Interaction*);
// should be set at every step by InteractionLoop
Scene* scene;
@@ -25,6 +26,7 @@
class BodyCallback: public Serializable{
public:
+ virtual ~BodyCallback(); // vtable
typedef void(*FuncPtr)(BodyCallback*,Body*);
// set at every step, before stepInit() is called
Scene* scene;
=== modified file 'pkg/common/Dispatching.cpp'
--- pkg/common/Dispatching.cpp 2010-10-18 19:27:07 +0000
+++ pkg/common/Dispatching.cpp 2010-10-26 13:41:30 +0000
@@ -101,7 +101,8 @@
updateScenePtr();
shared_ptr<BodyContainer>& bodies = scene->bodies;
- Matrix3r cellHsize; if(scene->isPeriodic) cellHsize=scene->cell->Hsize;
+ const bool isPeriodic(scene->isPeriodic);
+ Matrix3r cellHsize; if(isPeriodic) cellHsize=scene->cell->Hsize;
bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->iter);
#ifdef YADE_OPENMP
const long size=scene->interactions->size();
@@ -123,7 +124,7 @@
bool wasReal=I->isReal();
if (!b1->shape || !b2->shape) { assert(!wasReal); continue; } // some bodies do not have shape
bool geomCreated;
- if(!scene->isPeriodic){
+ if(!isPeriodic){
geomCreated=operator()(b1->shape, b2->shape, *b1->state, *b2->state, Vector3r::Zero(), /*force*/ false, I);
} else{
Vector3r shift2=cellHsize*I->cellDist.cast<Real>();
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2010-10-18 19:27:07 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2010-10-26 13:41:30 +0000
@@ -63,10 +63,17 @@
{
scene->forces.sync();
Real dt=scene->dt;
+ if(scene->isPeriodic && homotheticCellResize>=0){
+ LOG_WARN("Newton.homotheticCellResize is deprecated, use Cell.homoDeform instead. Setting Cell.homoDeform for you now, but the compatibility interface will be removed in the future.");
+ scene->cell->homoDeform=(homotheticCellResize==0?Cell::HOMO_NONE:(homotheticCellResize==1?Cell::HOMO_VEL:Cell::HOMO_VEL_2ND));
+ homotheticCellResize=-1;
+ }
+ homoDeform=(scene->isPeriodic ? scene->cell->homoDeform : -1); // -1 for aperiodic simulations
+ dVelGrad=scene->cell->velGrad-prevVelGrad;
// account for motion of the periodic boundary, if we remember its last position
// its velocity will count as max velocity of bodies
// otherwise the collider might not run if only the cell were changing without any particle motion
- if(scene->isPeriodic && (prevCellSize!=scene->cell->getSize())){ cellChanged=true; maxVelocitySq=(prevCellSize-scene->cell->getSize()).squaredNorm()/pow(dt,2); }
+ if(scene->isPeriodic && ((prevCellSize!=scene->cell->getSize())) && /* initial value */!isnan(prevCellSize[0]) ){ cellChanged=true; maxVelocitySq=(prevCellSize-scene->cell->getSize()).squaredNorm()/pow(dt,2); }
else { maxVelocitySq=0; cellChanged=false; }
haveBins=(bool)velocityBins;
if(haveBins) velocityBins->binVelSqInitialize(maxVelocitySq);
@@ -108,7 +115,7 @@
if(isnan(f[0])||isnan(f[1])||isnan(f[2])) throw runtime_error(("NewtonIntegrator: NaN force acting on #"+lexical_cast<string>(id)+".").c_str());
#endif
state->accel=f/state->mass;
- if (!scene->isPeriodic || homotheticCellResize==0)
+ if (!scene->isPeriodic || homoDeform==Cell::HOMO_NONE)
cundallDamp(dt,f,state->vel,state->accel);
else {//Apply damping on velocity fluctuations only rather than true velocity meanfield+fluctuation.
Vector3r velFluctuation(state->vel-prevVelGrad*state->pos);
@@ -181,7 +188,7 @@
FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
#endif
if(haveBins) velocityBins->binVelSqFinalize();
- if(scene->isPeriodic) { prevCellSize=scene->cell->getSize();prevVelGrad=scene->cell->velGrad; }
+ if(scene->isPeriodic) { prevCellSize=scene->cell->getSize(); prevVelGrad=scene->cell->velGrad; }
}
inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const Body::id_t& id, const Real& dt)
@@ -189,16 +196,18 @@
blockTranslateDOFs(state->blockedDOFs, state->accel);
if (scene->forces.getMoveRotUsed()) state->pos+=scene->forces.getMove(id);
- if (homotheticCellResize>0) {
+ if (homoDeform==Cell::HOMO_VEL || homoDeform==Cell::HOMO_VEL_2ND) {
// update velocity reflecting changes in the macroscopic velocity field, making the problem homothetic.
//NOTE : if the velocity is updated before moving the body, it means the current velGrad (i.e. before integration in cell->integrateAndUpdate) will be effective for the current time-step. Is it correct? If not, this velocity update can be moved just after "state->pos += state->vel*dt", meaning the current velocity impulse will be applied at next iteration, after the contact law. (All this assuming the ordering is resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw. Any other might fool us.)
//NOTE : dVel defined without wraping the coordinates means bodies out of the (0,0,0) period can move realy fast. It has to be compensated properly in the definition of relative velocities (see Ig2 functors and contact laws).
//This is the convective term, appearing in the time derivation of Cundall/Thornton expression (dx/dt=velGrad*pos -> d²x/dt²=dvelGrad/dt+velGrad*vel), negligible in many cases but not for high speed large deformations (gaz or turbulent flow). Emulating Cundall is an option, I don't especially recommend it. I know homothetic 1 and 2 expressions tend to identical values in the limit of dense quasi-static situations. They can give slitghly different results in other cases, and I'm still not sure which one should be considered better, if any (Cundall formula is in contradiction with molecular dynamics litterature).
- if (homotheticCellResize>1) state->vel+=prevVelGrad*state->vel*dt;
+ if (homoDeform==Cell::HOMO_VEL_2ND) state->vel+=scene->cell->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)*state->pos;
+ Vector3r dVel=dVelGrad*state->pos;
state->vel+=dVel;
+ } else if (homoDeform==Cell::HOMO_POS){
+ state->pos+=scene->cell->velGrad*state->pos*dt;
}
state->vel+=dt*state->accel;
state->pos += state->vel*dt;
=== modified file 'pkg/dem/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp 2010-10-18 19:27:07 +0000
+++ pkg/dem/NewtonIntegrator.hpp 2010-10-26 13:41:30 +0000
@@ -34,10 +34,10 @@
Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q);
inline void blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v);
inline void blockRotateDOFs(unsigned blockedDOFs, Vector3r& v);
- // cell size from previous step, used to detect change, find max velocity and update positions if linearCellResize enabled
- Vector3r prevCellSize;
// whether the cell has changed from the previous step
bool cellChanged;
+ int homoDeform; // updated from scene at every call; -1 for aperiodic simulations, otherwise equal to scene->cell->homoDeform
+ Matrix3r dVelGrad; // dtto
public:
#ifdef YADE_OPENMP
@@ -50,12 +50,13 @@
((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 :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."))
+ ((int,homotheticCellResize,-1,Attr::hidden,"[This attribute is deprecated, use Cell::homoDeform instead.]"))
+ // 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."))
+ ((Vector3r,prevCellSize,Vector3r(NaN,NaN,NaN),Attr::hidden,"cell size from previous step, used to detect change and find max velocity"))
,
/*ctor*/
- prevCellSize=Vector3r::Zero();
#ifdef YADE_OPENMP
threadMaxVelocitySq.resize(omp_get_max_threads());
#endif
=== modified file 'pkg/dem/ScGeom.cpp'
--- pkg/dem/ScGeom.cpp 2010-10-20 11:31:49 +0000
+++ pkg/dem/ScGeom.cpp 2010-10-26 13:41:30 +0000
@@ -30,7 +30,7 @@
//Update contact normal
normal=currentNormal;
//Precompute shear increment
- Vector3r relativeVelocity=getIncidentVel(&rbp1,&rbp2,scene->dt,shift2,scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>()) : Vector3r::Zero(),avoidGranularRatcheting);
+ Vector3r relativeVelocity=getIncidentVel(&rbp1,&rbp2,scene->dt,shift2,scene->isPeriodic ? scene->cell->intrShiftVel(c->cellDist) : Vector3r::Zero(),avoidGranularRatcheting);
//keep the shear part only
relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
shearInc = relativeVelocity*scene->dt;
@@ -80,8 +80,8 @@
Scene* scene=Omega::instance().getScene().get();
return getIncidentVel(Body::byId(i->getId1(),scene)->state.get(),Body::byId(i->getId2(),scene)->state.get(),
scene->dt,
- scene->isPeriodic ? Vector3r( scene->cell->Hsize*i->cellDist.cast<Real>()): Vector3r::Zero(), // shift2
- scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*i->cellDist.cast<Real>()) : Vector3r::Zero(), // shiftVel
+ scene->isPeriodic ? scene->cell->intrShiftPos(i->cellDist): Vector3r::Zero(), // shift2
+ scene->isPeriodic ? scene->cell->intrShiftVel(i->cellDist): Vector3r::Zero(), // shiftVel
avoidGranularRatcheting);
}
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2010-10-22 14:30:53 +0000
+++ pkg/dem/Shop.cpp 2010-10-26 13:41:30 +0000
@@ -192,10 +192,8 @@
Vector3r vel=b->state->vel;
if(isPeriodic){
/* Only take in account the fluctuation velocity, not the mean velocity of homothetic resize. */
- vel-=scene->cell->velGrad*state->pos;
- // 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
+ vel=scene->cell->bodyFluctuationVel(state->pos,vel);
+ // 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());
=== modified file 'py/system.py'
--- py/system.py 2010-10-13 16:00:36 +0000
+++ py/system.py 2010-10-26 13:41:30 +0000
@@ -21,15 +21,6 @@
return ret | ret2
_allSerializables=childClasses('Serializable')
-# classes that cannot be instantiated in python directly, and will have no properties generated for them
-_noPropsClasses=set(['InteractionContainer','BodyContainer','Functor','Engine','Dispatcher'])
-# classes that have special wrappers; only the most-bottom ones, with their names as it is in c++
-_pyRootClasses=set([
- 'GlobalEngine','PartialEngine','Shape','Bound','InteractionGeometry','InteractionPhysics','FileGenerator',
- 'BoundFunctor','IGeomFunctor','InteractionPhysicsFunctor','LawFunctor','Material','State']
- # childless classes
- +['BoundDispatcher','IGeomDispatcher','IPhysDispatcher','LawDispatcher','InteractionDispatchers','ParallelEngine']
-)
## set of classes for which the proxies were created
_proxiedClasses=set()
@@ -127,15 +118,11 @@
'Ip2_BMP_BMP_CSPhys':'Ip2_2xFrictMat_CSPhys', # Wed Mar 10 15:08:56 2010, eudoxos@frigo
'CinemDTCEngine':'KinemCTDEngine', # Tue Mar 16 13:54:21 2010, jduriez@c1solimara-l
'NormalInelasticityLaw':'Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity', # Wed Mar 17 15:50:59 2010, jduriez@c1solimara-l
- 'OldName':'NewName', # Tue Mar 30 14:11:13 2010, sch50p@fluent-ph
'CapillaryCohesiveLaw':'CapillaryLaw', # Tue Mar 30 14:11:36 2010, sch50p@fluent-ph
- 'Simplecd':'/home/sch50p/YADE-bzr/yade/pkg/dem/PreProcessor', # Tue Mar 30 14:13:03 2010, sch50p@fluent-ph
'SimpleElasticRelationshipsWater':'Ip2_Frictmat_FrictMat_CapillaryLawPhys', # Tue Mar 30 14:20:36 2010, sch50p@fluent-ph
- 'OldName':'NewName', # Wed Mar 31 09:22:48 2010, sch50p@fluent-ph
'CapillaryLaw':'Law2_ScGeom_CapillaryPhys_Capillarity', # Wed Mar 31 09:23:36 2010, sch50p@fluent-ph
'CapillaryParameters':'CapillaryPhys', # Wed Mar 31 09:25:03 2010, sch50p@fluent-ph
'Ip2_FrictMat_FrictMat_CapillaryLawPhys':'Ip2_FrictMat_FrictMat_CapillaryPhys', # Wed Mar 31 09:26:04 2010, sch50p@fluent-ph
- 'Ip2_Frictmat_FrictMat_CapillaryLawPhys':'Ip2_FrictMat_FrictMat_CapillaryPhys', # Wed Mar 31 09:26:56 2010, sch50p@fluent-ph
'SimpleViscoelasticMat':'ViscElMat', # Fri Apr 9 19:25:38 2010, vaclav@flux
'SimpleViscoelasticPhys':'ViscElPhys', # Fri Apr 9 19:26:34 2010, vaclav@flux
'Law2_Spheres_Viscoelastic_SimpleViscoelastic':'Law2_ScGeom_ViscElPhys_Basic', # Fri Apr 9 19:28:02 2010, vaclav@flux
=== modified file 'py/tests/omega.py'
--- py/tests/omega.py 2010-10-19 14:57:27 +0000
+++ py/tests/omega.py 2010-10-26 13:41:30 +0000
@@ -16,9 +16,24 @@
class TestInteractions(unittest.TestCase): pass
class TestForce(unittest.TestCase): pass
class TestEngines(unittest.TestCase): pass
-class TestIO(unittest.TestCase): pass
class TestTags(unittest.TestCase): pass
+class TestIO(unittest.TestCase):
+ def testSaveAllClasses(self):
+ 'All classes can be saved and loaded with boost::serialization'
+ import yade.system
+ failed=set()
+ for c in yade.system.childClasses('Serializable'):
+ O.reset()
+ try:
+ O.miscParams=[eval(c)()]
+ O.saveTmp()
+ O.loadTmp()
+ except (RuntimeError,ValueError):
+ failed.add(c)
+ self.assert_(len(failed)==0,'Failed classes were: '+' '.join(failed))
+
+
class TestCell(unittest.TestCase):
def setUp(self):
=== modified file 'py/tests/pbc.py'
--- py/tests/pbc.py 2010-10-19 14:57:27 +0000
+++ py/tests/pbc.py 2010-10-26 13:41:30 +0000
@@ -27,17 +27,34 @@
O.bodies.append(utils.sphere(self.initPos,.5))
#print O.bodies[1].state.pos
O.bodies[1].state.vel=self.initVel
+ O.engines=[NewtonIntegrator()]
+ O.cell.velGrad=Matrix3(0,0,0, 0,0,0, 0,0,-1)
+ O.cell.homoDeform=3
+ O.dt=0 # do not change positions with dt=0 in NewtonIntegrator, but still update velocities from velGrad
+ def testHomotheticResizeVel(self):
+ "PBC: (homoDeform==3) homothetic cell deformation adjusts particle velocity"
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):
- "PBC: 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 testHomotheticResizePos(self):
+ "PBC: (homoDeform==1) homothetic cell deformation adjusts particle position"
+ O.cell.homoDeform=1
+ O.step()
+ s1=O.bodies[1].state
+ self.assertAlmostEqual(s1.vel[2],self.initVel[2])
+ self.assertAlmostEqual(s1.pos[2],self.initPos[2]+self.initPos[2]*O.cell.velGrad[2,2]*O.dt)
def testScGeomIncidentVelocity(self):
- "PBC: ScGeom computes incident velocity correctly"
- O.dt=0 # do not change positions with dt=0 in NewtonIntegrator, but still update velocities from velGrad
+ "PBC: (homoDeform==3) ScGeom computes incident velocity correctly"
+ O.step()
+ O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
+ i=utils.createInteraction(0,1)
+ self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=True))
+ self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
+ self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
+ def testScGeomIncidentVelocity_homoPos(self):
+ "PBC: (homoDeform==1) ScGeom computes incident velocity correctly"
+ O.cell.homoDeform=1
O.step()
O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
i=utils.createInteraction(0,1)
@@ -45,13 +62,16 @@
self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
def testKineticEnergy(self):
- "PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
- O.dt=1e-50 # with timestep almost zero, just update the velocity in the integrator
+ "PBC: (homoDeform==3) utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
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()
self.assertAlmostEqual(Ek,utils.kineticEnergy())
+ def testKineticEnergy_homoPos(self):
+ "PBC: (homoDeform==1) utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
+ O.cell.homoDeform=1; O.step()
+ self.assertAlmostEqual(.5*O.bodies[1].state.mass*self.initVel.squaredNorm(),utils.kineticEnergy())
=== modified file 'py/ymport.py'
--- py/ymport.py 2010-06-29 22:46:24 +0000
+++ py/ymport.py 2010-10-26 13:41:30 +0000
@@ -69,8 +69,7 @@
facets=imp.ymport(file)
for b in facets:
b.dynamic=dynamic
- b.shape.postProcessAttributes(True)
- b.shape.Color=color if color else utils.randomColor()
+ b.shape.color=color if color else utils.randomColor()
b.shape.wire=wire
b.shape.highlight=highlight
pos,ori=b.state.pos,b.state.ori
Follow ups
