yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02515
[Branch ~yade-dev/yade/trunk] Rev 1858: 1. Add PeriTriaxController, with sample script in scripts/test/periodic-triax.py; works perfectly
------------------------------------------------------------
revno: 1858
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Sun 2009-12-06 00:12:46 +0100
message:
1. Add PeriTriaxController, with sample script in scripts/test/periodic-triax.py; works perfectly
2. Fix stress sign in PeriIsoCompressor (summing abs stresses still broken)
3. Handle periodic boundary moves with velocityBins
4. Rename ncb to scene in NewtonsDampedLaw (gradual change)
added:
scripts/test/periodic-triax.py
modified:
core/main/main.py.in
pkg/common/DataClass/VelocityBins.cpp
pkg/common/DataClass/VelocityBins.hpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp
pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp
py/pack.py
scripts/test/periodic-compress.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/main/main.py.in'
--- core/main/main.py.in 2009-12-05 10:45:14 +0000
+++ core/main/main.py.in 2009-12-05 23:12:46 +0000
@@ -66,7 +66,13 @@
if args[0].endswith('.py'):
print "Running script "+args[0]
yade.runtime.argv=args[1:]
- execfile(args[0])
+ try:
+ execfile(args[0])
+ except SystemExit: raise
+ except: # all other exceptions
+ import traceback
+ traceback.print_exc()
+ if opts.exitAfter: sys.exit(1)
if opts.exitAfter: sys.exit(0)
# show python command-line
=== modified file 'pkg/common/DataClass/VelocityBins.cpp'
--- pkg/common/DataClass/VelocityBins.cpp 2009-12-04 23:07:34 +0000
+++ pkg/common/DataClass/VelocityBins.cpp 2009-12-05 23:12:46 +0000
@@ -96,7 +96,7 @@
/* non-parallel implementations */
#ifdef YADE_OPENMP
- void VelocityBins::binVelSqInitialize(){ FOREACH(Bin& bin, bins){ FOREACH(Real& vSq, bin.threadMaxVelSq) vSq=0.; }}
+ void VelocityBins::binVelSqInitialize(Real vSqInit){ FOREACH(Bin& bin, bins){ FOREACH(Real& vSq, bin.threadMaxVelSq) vSq=vSqInit; }}
void VelocityBins::binVelSqUse(body_id_t id, Real velSq){
Real& maxVelSq(bins[bodyBins[id]].threadMaxVelSq[omp_get_thread_num()]);
maxVelSq=max(maxVelSq,velSq);
@@ -108,7 +108,7 @@
}
}
#else
- void VelocityBins::binVelSqInitialize(){ FOREACH(Bin& bin, bins) bin.currMaxVelSq=0; }
+ void VelocityBins::binVelSqInitialize(Real vSqInit=0.){ FOREACH(Bin& bin, bins) bin.currMaxVelSq=vSqInit; }
void VelocityBins::binVelSqUse(body_id_t id, Real velSq){
Real& maxVelSq(bins[bodyBins[id]].currMaxVelSq);
maxVelSq=max(maxVelSq,velSq);
=== modified file 'pkg/common/DataClass/VelocityBins.hpp'
--- pkg/common/DataClass/VelocityBins.hpp 2009-12-04 23:07:34 +0000
+++ pkg/common/DataClass/VelocityBins.hpp 2009-12-05 23:12:46 +0000
@@ -69,7 +69,7 @@
body in that bin and binVelSqFinalize() will not do nothing.
*/
// reset per-bin max velocities
- void binVelSqInitialize();
+ void binVelSqInitialize(Real velSqInit=0.);
// use body max velocity -- called for every body at every step (from NewtonsDampedLaw, normally)
void binVelSqUse(body_id_t id, Real velSq);
// actually assign max velocities to their respective bins
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-12-04 23:07:34 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-12-05 23:12:46 +0000
@@ -31,26 +31,26 @@
if((blockedDOFs & State::DOF_RY)!=0) v[1]=0;
if((blockedDOFs & State::DOF_RZ)!=0) v[2]=0;
}
-void NewtonsDampedLaw::handleClumpMemberAccel(Scene* ncb, const body_id_t& memberId, State* memberState, State* clumpState){
- const Vector3r& f=ncb->bex.getForce(memberId);
+void NewtonsDampedLaw::handleClumpMemberAccel(Scene* scene, const body_id_t& memberId, State* memberState, State* clumpState){
+ const Vector3r& f=scene->bex.getForce(memberId);
Vector3r diffClumpAccel=f/clumpState->mass;
// damp increment of accels on the clump, using velocities of the clump MEMBER
- cundallDamp(ncb->dt,f,memberState->vel,diffClumpAccel);
+ cundallDamp(scene->dt,f,memberState->vel,diffClumpAccel);
clumpState->accel+=diffClumpAccel;
}
-void NewtonsDampedLaw::handleClumpMemberAngAccel(Scene* ncb, const body_id_t& memberId, State* memberState, State* clumpState){
+void NewtonsDampedLaw::handleClumpMemberAngAccel(Scene* scene, const body_id_t& memberId, State* memberState, State* clumpState){
// angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
- const Vector3r& m=ncb->bex.getTorque(memberId); const Vector3r& f=ncb->bex.getForce(memberId);
+ const Vector3r& m=scene->bex.getTorque(memberId); const Vector3r& f=scene->bex.getForce(memberId);
Vector3r diffClumpAngularAccel=diagDiv(m,clumpState->inertia)+diagDiv((memberState->pos-clumpState->pos).Cross(f),clumpState->inertia);
// damp increment of accels on the clump, using velocities of the clump MEMBER
- cundallDamp(ncb->dt,m,memberState->angVel,diffClumpAngularAccel);
+ cundallDamp(scene->dt,m,memberState->angVel,diffClumpAngularAccel);
clumpState->angAccel+=diffClumpAngularAccel;
}
-void NewtonsDampedLaw::handleClumpMemberTorque(Scene* ncb, const body_id_t& memberId, State* memberState, State* clumpState, Vector3r& M){
- const Vector3r& m=ncb->bex.getTorque(memberId); const Vector3r& f=ncb->bex.getForce(memberId);
+void NewtonsDampedLaw::handleClumpMemberTorque(Scene* scene, const body_id_t& memberId, State* memberState, State* clumpState, Vector3r& M){
+ const Vector3r& m=scene->bex.getTorque(memberId); const Vector3r& f=scene->bex.getForce(memberId);
M+=(memberState->pos-clumpState->pos).Cross(f)+m;
}
-void NewtonsDampedLaw::saveMaximaVelocity(Scene* ncb, const body_id_t& id, State* state){
+void NewtonsDampedLaw::saveMaximaVelocity(Scene* scene, const body_id_t& id, State* state){
if(haveBins) velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(state));
#ifdef YADE_OPENMP
Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,state->vel.SquaredLength());
@@ -59,89 +59,94 @@
#endif
}
-void NewtonsDampedLaw::action(Scene * ncb)
+void NewtonsDampedLaw::action(Scene* scene)
{
- ncb->bex.sync();
- Real dt=Omega::instance().getTimeStep();
- maxVelocitySq=-1;
+ scene->bex.sync();
+ Real dt=scene->dt;
+ // 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!=Vector3r::ZERO){ maxVelocitySq=(prevCellSize-(scene->cellMax-scene->cellMin)).SquaredLength()/pow(dt,2); }
+ else maxVelocitySq=0;
haveBins=(bool)velocityBins;
- if(haveBins) velocityBins->binVelSqInitialize();
+ if(haveBins) velocityBins->binVelSqInitialize(maxVelocitySq);
+ if(scene->isPeriodic) prevCellSize=scene->cellMax-scene->cellMin;
#ifdef YADE_OPENMP
FOREACH(Real& thrMaxVSq, threadMaxVelocitySq) { thrMaxVSq=0; }
- const BodyContainer& bodies=*(ncb->bodies.get());
+ const BodyContainer& bodies=*(scene->bodies.get());
const long size=(long)bodies.size();
#pragma omp parallel for schedule(static)
for(long _id=0; _id<size; _id++){
const shared_ptr<Body>& b(bodies[_id]);
#else
- FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
#endif
if(!b) continue;
State* state=b->state.get();
const body_id_t& id=b->getId();
// clump members are non-dynamic; we only get their velocities here
if (!b->isDynamic || b->isClumpMember()){
- saveMaximaVelocity(ncb,id,state);
+ saveMaximaVelocity(scene,id,state);
continue;
}
if (b->isStandalone()){
// translate equation
- const Vector3r& f=ncb->bex.getForce(id);
+ const Vector3r& f=scene->bex.getForce(id);
state->accel=f/state->mass;
cundallDamp(dt,f,state->vel,state->accel);
- leapfrogTranslate(ncb,state,id,dt);
+ leapfrogTranslate(scene,state,id,dt);
// rotate equation
if (/*b->isSpheral || */ !exactAsphericalRot){ // spheral body or exactAsphericalRot disabled
- const Vector3r& m=ncb->bex.getTorque(id);
+ const Vector3r& m=scene->bex.getTorque(id);
state->angAccel=diagDiv(m,state->inertia);
cundallDamp(dt,m,state->angVel,state->angAccel);
- leapfrogSphericalRotate(ncb,state,id,dt);
+ leapfrogSphericalRotate(scene,state,id,dt);
} else { // non spheral body and exactAsphericalRot enabled
- const Vector3r& m=ncb->bex.getTorque(id);
- leapfrogAsphericalRotate(ncb,state,id,dt,m);
+ const Vector3r& m=scene->bex.getTorque(id);
+ leapfrogAsphericalRotate(scene,state,id,dt,m);
}
} else if (b->isClump()){
// clump mass forces
- const Vector3r& f=ncb->bex.getForce(id);
+ const Vector3r& f=scene->bex.getForce(id);
Vector3r dLinAccel=f/state->mass;
cundallDamp(dt,f,state->vel,dLinAccel);
state->accel+=dLinAccel;
- const Vector3r& m=ncb->bex.getTorque(id);
+ const Vector3r& m=scene->bex.getTorque(id);
Vector3r M(m);
// sum force on clump memebrs
if (exactAsphericalRot){
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,ncb); assert(b->isClumpMember());
+ const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(b->isClumpMember());
State* memberState=member->state.get();
- handleClumpMemberAccel(ncb,memberId,memberState,state);
- handleClumpMemberTorque(ncb,memberId,memberState,state,M);
- saveMaximaVelocity(ncb,memberId,memberState);
+ handleClumpMemberAccel(scene,memberId,memberState,state);
+ handleClumpMemberTorque(scene,memberId,memberState,state,M);
+ saveMaximaVelocity(scene,memberId,memberState);
}
// motion
- leapfrogTranslate(ncb,state,id,dt);
- leapfrogAsphericalRotate(ncb,state,id,dt,M);
+ leapfrogTranslate(scene,state,id,dt);
+ leapfrogAsphericalRotate(scene,state,id,dt,M);
} else { // exactAsphericalRot disabled
Vector3r dAngAccel=diagDiv(M,state->inertia);
cundallDamp(dt,M,state->angVel,dAngAccel);
state->angAccel+=dAngAccel;
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,ncb); assert(member->isClumpMember());
+ const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
State* memberState=member->state.get();
- handleClumpMemberAccel(ncb,memberId,memberState,state);
- handleClumpMemberAngAccel(ncb,memberId,memberState,state);
- saveMaximaVelocity(ncb,memberId,memberState);
+ handleClumpMemberAccel(scene,memberId,memberState,state);
+ handleClumpMemberAngAccel(scene,memberId,memberState,state);
+ saveMaximaVelocity(scene,memberId,memberState);
}
// motion
- leapfrogTranslate(ncb,state,id,dt);
- leapfrogSphericalRotate(ncb,state,id,dt);
+ leapfrogTranslate(scene,state,id,dt);
+ leapfrogSphericalRotate(scene,state,id,dt);
}
static_cast<Clump*>(b.get())->moveMembers();
}
- saveMaximaVelocity(ncb,id,state);
+ saveMaximaVelocity(scene,id,state);
}
#ifdef YADE_OPENMP
FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
@@ -149,23 +154,23 @@
if(haveBins) velocityBins->binVelSqFinalize();
}
-inline void NewtonsDampedLaw::leapfrogTranslate(Scene* ncb, State* state, const body_id_t& id, const Real& dt )
+inline void NewtonsDampedLaw::leapfrogTranslate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
{
blockTranslateDOFs(state->blockedDOFs, state->accel);
state->vel+=dt*state->accel;
- state->pos += state->vel*dt + ncb->bex.getMove(id);
+ state->pos += state->vel*dt + scene->bex.getMove(id);
}
-inline void NewtonsDampedLaw::leapfrogSphericalRotate(Scene* ncb, State* state, const body_id_t& id, const Real& dt )
+inline void NewtonsDampedLaw::leapfrogSphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
{
blockRotateDOFs(state->blockedDOFs, state->angAccel);
state->angVel+=dt*state->angAccel;
Vector3r axis = state->angVel; Real angle = axis.Normalize();
Quaternionr q; q.FromAxisAngle ( axis,angle*dt );
state->ori = q*state->ori;
- if(ncb->bex.getMoveRotUsed() && ncb->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(ncb->bex.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
+ if(scene->bex.getMoveRotUsed() && scene->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(scene->bex.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
state->ori.Normalize();
}
-void NewtonsDampedLaw::leapfrogAsphericalRotate(Scene* ncb, State* state, const body_id_t& id, const Real& dt, const Vector3r& M){
+void NewtonsDampedLaw::leapfrogAsphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt, const Vector3r& M){
Matrix3r A; state->ori.Conjugate().ToRotationMatrix(A); // rotation matrix from global to local r.f.
const Vector3r l_n = state->angMom + dt/2 * M; // global angular momentum at time n
const Vector3r l_b_n = A*l_n; // local angular momentum at time n
@@ -179,7 +184,7 @@
const Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
state->ori+=dt*dotQ_half; // Q at time n+1
state->angVel=state->ori.Rotate(angVel_b_half); // global angular velocity at time n+1/2
- if(ncb->bex.getMoveRotUsed() && ncb->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(ncb->bex.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
+ if(scene->bex.getMoveRotUsed() && scene->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(scene->bex.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
state->ori.Normalize();
}
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-12-04 23:07:34 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-12-05 23:12:46 +0000
@@ -50,6 +50,9 @@
Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q);
inline void blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v);
inline void blockRotateDOFs(unsigned blockedDOFs, Vector3r& v);
+ Vector3r prevCellSize;
+
+
public:
///damping coefficient for Cundall's non viscous damping
@@ -58,13 +61,14 @@
Real maxVelocitySq;
/// Enable of the exact aspherical body rotation integrator
bool exactAsphericalRot;
+
#ifdef YADE_OPENMP
vector<Real> threadMaxVelocitySq;
#endif
/// velocity bins (not used if not created)
shared_ptr<VelocityBins> velocityBins;
virtual void action(Scene *);
- NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false){
+ NewtonsDampedLaw(): prevCellSize(Vector3r::ZERO),damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false){
#ifdef YADE_OPENMP
threadMaxVelocitySq.resize(omp_get_max_threads());
#endif
=== modified file 'pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp 2009-12-04 23:46:45 +0000
+++ pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp 2009-12-05 23:12:46 +0000
@@ -4,10 +4,12 @@
#include<yade/pkg-dem/PeriIsoCompressor.hpp>
#include<yade/pkg-dem/Shop.hpp>
#include<yade/core/Scene.hpp>
+#include<yade/pkg-common/NormalShearInteractions.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
using namespace std;
-YADE_PLUGIN((PeriIsoCompressor))
+YADE_PLUGIN((PeriIsoCompressor)(PeriTriaxController))
CREATE_LOGGER(PeriIsoCompressor);
@@ -24,7 +26,7 @@
}
if(maxSpan<=0){
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
- if(!b->bound) continue;
+ if(!b || !b->bound) continue;
for(int i=0; i<3; i++) maxSpan=max(maxSpan,b->bound->max[i]-b->bound->min[i]);
}
@@ -37,7 +39,7 @@
if(minSize<2.1*maxSpan){ throw runtime_error("Minimum cell size is smaller than 2.1*span_of_the_biggest_body! (periodic collider requirement)"); }
if(((step%globalUpdateInt)==0) || avgStiffness<0 || sigma[0]<0 || sigma[1]<0 || sigma[2]<0){
Vector3r sumForces=Shop::totalForceInVolume(avgStiffness,rb);
- sigma=Vector3r(sumForces[0]/cellArea[0],sumForces[1]/cellArea[1],sumForces[2]/cellArea[2]);
+ sigma=-Vector3r(sumForces[0]/cellArea[0],sumForces[1]/cellArea[1],sumForces[2]/cellArea[2]);
LOG_TRACE("Updated sigma="<<sigma<<", avgStiffness="<<avgStiffness);
}
Real sigmaGoal=stresses[state]; assert(sigmaGoal>0);
@@ -47,11 +49,11 @@
bool allStressesOK=true;
if(keepProportions){ // the same algo as below, but operating on quantitites averaged over all dimensions
Real sigAvg=(sigma[0]+sigma[1]+sigma[2])/3., avgArea=(cellArea[0]+cellArea[1]+cellArea[2])/3., avgSize=(cellSize[0]+cellSize[1]+cellSize[2])/3.;
- Real avgGrow=1e-4*(sigAvg-sigmaGoal)*avgArea/(avgStiffness>0?avgStiffness:1);
+ Real avgGrow=1e-4*(sigmaGoal-sigAvg)*avgArea/(avgStiffness>0?avgStiffness:1);
Real maxToAvg=maxSize/avgSize;
if(abs(maxToAvg*avgGrow)>maxDisplPerStep) avgGrow=Mathr::Sign(avgGrow)*maxDisplPerStep/maxToAvg;
avgGrow=max(avgGrow,-(minSize-2.1*maxSpan)/maxToAvg);
- if(avgStiffness>0) { sigma-=(avgGrow*avgStiffness)*Vector3r::ONE; sigAvg-=avgGrow*avgStiffness; }
+ if(avgStiffness>0) { sigma+=(avgGrow*avgStiffness)*Vector3r::ONE; sigAvg+=avgGrow*avgStiffness; }
if(abs((sigAvg-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
cellGrow=(avgGrow/avgSize)*cellSize;
}
@@ -61,11 +63,11 @@
// FIXME: either NormalShearInteraction::{kn,ks} is computed wrong or we have dimensionality problem here
// FIXME: that is why the fixup 1e-4 is needed here
// FIXME: or perhaps maxDisplaPerStep=1e-2*charLen is too big??
- cellGrow[axis]=1e-4*(sigma[axis]-sigmaGoal)*cellArea[axis]/(avgStiffness>0?avgStiffness:1);
+ cellGrow[axis]=1e-4*(sigmaGoal-sigma[axis])*cellArea[axis]/(avgStiffness>0?avgStiffness:1); // FIXME: wrong dimensions? See PeriTriaxController
if(abs(cellGrow[axis])>maxDisplPerStep) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxDisplPerStep;
cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxSpan));
// crude way of predicting sigma, for steps when it is not computed from intrs
- if(avgStiffness>0) sigma[axis]-=cellGrow[axis]*avgStiffness;
+ if(avgStiffness>0) sigma[axis]+=cellGrow[axis]*avgStiffness; // FIXME: dimensions
if(abs((sigma[axis]-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
}
}
@@ -87,4 +89,106 @@
}
}
-
+void PeriTriaxController::strainStressStiffUpdate(){
+ // update strain first
+ const Vector3r cellSize(scene->cellMax-scene->cellMin);
+ for(int i=0; i<3; i++) strain[i]=(cellSize[i]-refSize[i])/refSize[i];
+ // stress and stiffness
+ Vector3r sumForce(Vector3r::ZERO), sumStiff(Vector3r::ZERO), sumLength(Vector3r::ZERO);
+ int n=0;
+ FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
+ if(!I->isReal()) continue;
+ n++;
+ NormalShearInteraction* nsi=YADE_CAST<NormalShearInteraction*>(I->interactionPhysics.get());
+ GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*>(I->interactionGeometry.get());
+ for(int i=0; i<3; i++){
+ const Real absNormI=abs(gsc->normal[i]);
+ Real f=nsi->normalForce[i]+nsi->shearForce[i];
+ // force is as applied to id1, and normal is from id1 towards id2;
+ // therefore if they have the same sense, it is tensile force (positive)
+ sumForce[i]+=Mathr::Sign(f*gsc->normal[i])*f;
+ sumLength[i]+=absNormI*(gsc->refR1+gsc->refR2);
+ sumStiff[i]+=absNormI*nsi->kn+(1-absNormI)*nsi->ks;
+ }
+ }
+ if(n>0){
+ // sumForce/area is sum of stresses on fictious slices; the number of slices is (average length)/length in the respective sense
+ for(int i=0; i<3; i++) stress[i]=((1./n)*sumLength[i]/cellSize[i])*sumForce[i]/(cellSize[(i+1)%3]*cellSize[(i+2)%3]);
+ stiff=(1./n)*sumStiff;
+ } else { stiff=Vector3r::ZERO; stress=Vector3r::ZERO; }
+}
+
+CREATE_LOGGER(PeriTriaxController);
+
+void PeriTriaxController::action(Scene* scene){
+ if(!scene->isPeriodic){ throw runtime_error("PeriTriaxController run on aperiodic simulation."); }
+ Vector3r cellSize=scene->cellMax-scene->cellMin;
+ Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
+ // initial updates
+ if(refSize[0]<0) refSize=cellSize;
+ if(maxBodySpan[0]<=0){
+ FOREACH(const shared_ptr<Body>& b,*scene->bodies){
+ if(!b || !b->bound) continue;
+ for(int i=0; i<3; i++) maxBodySpan[i]=max(maxBodySpan[i],b->bound->max[i]-b->bound->min[i]);
+ }
+ }
+ // check current size
+ if(2.1*maxBodySpan[0]>cellSize[0] || 2.1*maxBodySpan[1]>cellSize[1] || 2.1*maxBodySpan[2]>cellSize[2]){
+ LOG_DEBUG("cellSize="<<cellSize<<", maxBodySpan="<<maxBodySpan);
+ throw runtime_error("Minimum cell size is smaller than 2.1*maxBodySpan (periodic collider requirement)");
+ }
+ bool doUpdate((scene->currentIteration%globUpdate)==0);
+ if(doUpdate || min(stiff[0],min(stiff[1],stiff[2]))<=0){ strainStressStiffUpdate(); }
+
+ bool allOk=true; Vector3r cellGrow(Vector3r::ZERO);
+ // apply condition along each axis separately (stress or strain)
+ for(int axis=0; axis<3; axis++){
+ Real maxGrow=refSize[axis]*maxStrainRate*scene->dt;
+ if(stressMask & 1<<axis){ // control stress
+ // stiffness K=EA; σâ‚=goal stress; Δσ wanted stress difference to apply
+ // ΔεE=(Δl/lâ‚€)(K/A) = Δσ=σâ‚-σ ↔ Δl=(σâ‚-σ)lâ‚€(A/K)
+ cellGrow[axis]=(goal[axis]-stress[axis])*refSize[axis]*cellArea[axis]/(stiff[axis]>0?stiff[axis]:1.);
+ LOG_TRACE(axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis]);
+ } else { // control strain
+ cellGrow[axis]=(goal[axis]-strain[axis])*refSize[axis];
+ LOG_TRACE(axis<<": strain="<<strain[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis]);
+ }
+ // limit maximum strain rate
+ if(abs(cellGrow[axis])>maxGrow) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxGrow;
+ // do not shrink below minimum cell size (periodic collider condition), although it is suboptimal WRT resulting stress
+ cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxBodySpan[axis]));
+ // steady evolution with fluctuations; see TriaxialStressController
+ cellGrow[axis]=(1-growDamping)*cellGrow[axis]+.8*prevGrow[axis];
+ // crude way of predicting stress, for steps when it is not computed from intrs
+ if(stiff[axis]>0) stress[axis]+=(cellGrow[axis]/refSize[axis])*(stiff[axis]/cellArea[axis]);
+ strain[axis]+=cellGrow[axis]/refSize[axis];
+ if(doUpdate) LOG_DEBUG(axis<<": cellGrow="<<cellGrow[axis]<<", new stress="<<stress[axis]<<", new strain="<<strain[axis]);
+ // signal if condition not satisfied
+ if(stressMask & 1<<axis){
+ Real curr=stress[axis];
+ if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>relStressTol) || abs(curr-goal[axis])>absStressTol) allOk=false;
+ } else {
+ Real curr=strain[axis];
+ // since strain is prescribed exactly, tolerances need just to accomodate rounding issues
+ if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>1e-6) || abs(curr-goal[axis])>1e-6) {
+ allOk=false;
+ if(doUpdate) LOG_DEBUG("Strain not OK; "<<abs(curr-goal[axis])<<">1e-6");
+ }
+ }
+ // change cell size now
+ scene->cellMax+=cellGrow;
+ }
+ if(allOk){
+ if(doUpdate || currUnbalanced<0) {
+ currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,scene);
+ LOG_DEBUG("Stress/strain="<<(stressMask&1?stress[0]:strain[0])<<","<<(stressMask&2?stress[1]:strain[1])<<","<<(stressMask&4?stress[2]:strain[2])<<", goal="<<goal<<", unbalanced="<<currUnbalanced);
+ }
+ if(currUnbalanced<maxUnbalanced){
+ // LOG_INFO("Goal reached, packing stable.");
+ if(!doneHook.empty()){
+ LOG_DEBUG("Running doneHook: "<<doneHook);
+ PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(doneHook.c_str()); PyGILState_Release(gstate);
+ } else { Omega::instance().stopSimulationLoop(); }
+ }
+ }
+}
=== modified file 'pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp 2009-12-04 23:07:34 +0000
+++ pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp 2009-12-05 23:12:46 +0000
@@ -33,4 +33,54 @@
REGISTER_CLASS_AND_BASE(PeriIsoCompressor,StandAloneEngine);
};
REGISTER_SERIALIZABLE(PeriIsoCompressor);
-
+/* Engine for independently controlling stress or strain in periodic simulations.
+
+strainStress contains absolute values for the controlled quantity, and stressMask determines
+meaning of those values (0 for strain, 1 for stress): e.g. ( 1<<0 | 1<<2 ) = 1 | 4 = 5 means that
+strainStress[0] and strainStress[2] are stress values, and strainStress[1] is strain.
+*/
+class PeriTriaxController: public StandAloneEngine{
+ public:
+ //! Desired stress or strain values (depending on stressMask)
+ Vector3r goal;
+ //! mask determining strain/stress (0/1) meaning for goal components
+ int stressMask;
+ //! Maximum strain rate of the periodic cell
+ Real maxStrainRate;
+ //! maximum unbalanced force (defaults to 1e-4)
+ Real maxUnbalanced;
+ //! Absolute stress tolerance (1e3)
+ Real absStressTol;
+ //! Relative stress tolerance (3e-5)
+ Real relStressTol;
+ //! Damping of cell resizing (0=perfect control, 1=no control at all); see also TriaxialStressController::wallDamping.
+ Real growDamping;
+ //! how often to recompute average stress, stiffness and unbalaced force (defaults to 100)
+ int globUpdate;
+ //! python command to be run when the desired state is reached
+ string doneHook;
+ //! reference cell size (set automatically; if refSize[0]<0 (initial default), the current size is referenced)
+ Vector3r refSize;
+ //! maximum body dimension (set automatically)
+ Vector3r maxBodySpan;
+
+
+ //! average stresses, updated at every step (only every globUpdate steps recomputed from interactions)
+ Vector3r stress;
+ //! cell strain, updated at every step
+ Vector3r strain;
+ //! average stiffness, updated at every step (only every globUpdate steps recomputed from interactions)
+ Vector3r stiff;
+ //! current unbalanced force (updated every globUpdate)
+ Real currUnbalanced;
+ //! previous cell grow
+ Vector3r prevGrow;
+
+ void action(Scene*);
+ void strainStressStiffUpdate();
+ PeriTriaxController(): goal(Vector3r::ZERO),stressMask(0),maxStrainRate(1.),maxUnbalanced(1e-4),absStressTol(1e3),relStressTol(3e-5),growDamping(.25),globUpdate(100),refSize(Vector3r(-1,-1,-1)),maxBodySpan(Vector3r(-1,-1,-1)),stress(Vector3r::ZERO),strain(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO){}
+ REGISTER_ATTRIBUTES(StandAloneEngine,(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(refSize)(stress)(strain)(stiff));
+ DECLARE_LOGGER;
+ REGISTER_CLASS_AND_BASE(PeriTriaxController,StandAloneEngine);
+};
+REGISTER_SERIALIZABLE(PeriTriaxController);
=== modified file 'py/pack.py'
--- py/pack.py 2009-12-04 23:27:26 +0000
+++ py/pack.py 2009-12-05 23:12:46 +0000
@@ -308,7 +308,7 @@
#print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.periodicCell
#print x1,y1,z1,radius,rRelFuzz
num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,rRelFuzz,spheresInCell,True)
- O.engines=[BexResetter(),BoundDispatcher([InteractingSphere2AABB()]),InsertionSortCollider(nBins=5,sweepLength=.05*radius),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[100e9,1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonsDampedLaw(damping=.6)]
+ O.engines=[BexResetter(),BoundDispatcher([InteractingSphere2AABB()]),InsertionSortCollider(nBins=5,sweepLength=.05*radius),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonsDampedLaw(damping=.6)]
O.materials.append(GranularMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3))
for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
O.dt=utils.PWaveTimeStep()
=== modified file 'scripts/test/periodic-compress.py'
--- scripts/test/periodic-compress.py 2009-12-04 23:27:26 +0000
+++ scripts/test/periodic-compress.py 2009-12-05 23:12:46 +0000
@@ -16,7 +16,7 @@
[SimpleElasticRelationships()],
[Law2_Dem3Dof_Elastic_Elastic()],
),
- PeriIsoCompressor(charLen=.5,stresses=[50e9,1e8],doneHook="print 'FINISHED'; O.pause() ",keepProportions=True),
+ PeriIsoCompressor(charLen=.5,stresses=[-50e9,-1e8],doneHook="print 'FINISHED'; O.pause() ",keepProportions=True),
NewtonsDampedLaw(damping=.4)
]
O.dt=utils.PWaveTimeStep()
=== added file 'scripts/test/periodic-triax.py'
--- scripts/test/periodic-triax.py 1970-01-01 00:00:00 +0000
+++ scripts/test/periodic-triax.py 2009-12-05 23:12:46 +0000
@@ -0,0 +1,44 @@
+# coding: utf-8
+# 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+"Test and demonstrate use of PeriTriaxController."
+from yade import *
+from yade import pack,log,qt
+log.setLevel('PeriTriaxController',log.DEBUG)
+O.periodicCell=(0,0,0),(.1,.1,.1)
+sp=pack.SpherePack()
+radius=5e-3
+num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
+O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
+
+
+O.engines=[
+ BexResetter(),
+ BoundDispatcher([InteractingSphere2AABB()]),
+ InsertionSortCollider(nBins=5,sweepLength=.05*radius),
+ InteractionDispatchers(
+ [Ig2_Sphere_Sphere_Dem3DofGeom()],
+ [SimpleElasticRelationships()],
+ [Law2_Dem3Dof_Elastic_Elastic()]
+ ),
+ PeriTriaxController(goal=[-1e5,-1e5,0],stressMask=3,globUpdate=5,maxStrainRate=1.,doneHook='triaxDone()',label='triax'),
+ NewtonsDampedLaw(damping=.6),
+]
+O.dt=utils.PWaveTimeStep()
+O.run();
+qt.View()
+
+phase=0
+def triaxDone():
+ global phase
+ if phase==0:
+ print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']
+ print 'Now εz will go from 0 to .2 while σx and σy will be kept the same.'
+ triax['goal']=[-1e5,-1e5,.2]
+ phase+=1
+ elif phase==1:
+ print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']
+ print 'Done, pausing now.'
+ O.pause()
+
+
+
