yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01183
[svn] r1743 - in trunk: extra extra/usct gui/py
Author: eudoxos
Date: 2009-04-02 11:28:06 +0200 (Thu, 02 Apr 2009)
New Revision: 1743
Modified:
trunk/extra/Brefcom.cpp
trunk/extra/Brefcom.hpp
trunk/extra/usct/UniaxialStrainControlledTest.cpp
trunk/extra/usct/UniaxialStrainControlledTest.hpp
trunk/gui/py/plot.py
trunk/gui/py/yadeControl.cpp
Log:
1. cleanup UniaxialStrainer code, add vars to control length of the acceleration phase and the ability to set absolute speed
2. EngineUnits within InteractionDispatchers can be labeled and accessed from python now.
3. A few minor things in the plot module
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-03-31 23:10:43 UTC (rev 1742)
+++ trunk/extra/Brefcom.cpp 2009-04-02 09:28:06 UTC (rev 1743)
@@ -144,7 +144,9 @@
const Real maxError=1e-12;
Real f, ret=0.;
for(int i=0; i<maxIter; i++){
- cummBetaIter++;
+ #ifdef YADE_DEBUG
+ cummBetaIter++;
+ #endif
Real aux=c*exp(N*ret)+exp(ret);
f=log(aux);
if(fabs(f)<maxError) return ret;
@@ -166,12 +168,13 @@
Real deltaDmgStrain=(epsN*omega-dmgStrain)*exp(beta);
dmgStrain+=deltaDmgStrain;
LOG_DEBUG("deltaDmgStrain="<<deltaDmgStrain<<", viscous overstress "<<(epsN*omega-dmgStrain)*E);
+ /* σN=Kn(εN-εd); dmgOverstress=σN-(1-ω)*Kn*εN=…=Kn(ω*εN-εd) */
return (epsN*omega-dmgStrain)*E;
}
Real BrefcomContact::computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt){
if(sigmaTNorm<sigmaTYield) return 1.;
- Real c=/* should this be sigmaT0?? */ sigmaTNorm*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-1.);
+ Real c=sigmaTNorm*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-1.);
Real beta=solveBeta(c,plRateExp);
return 1.-exp(beta)*(1-sigmaTYield/sigmaTNorm);
}
Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp 2009-03-31 23:10:43 UTC (rev 1742)
+++ trunk/extra/Brefcom.hpp 2009-04-02 09:28:06 UTC (rev 1743)
@@ -172,8 +172,6 @@
class ef2_Spheres_Brefcom_BrefcomLaw: public ConstitutiveLaw{
public:
- /*! Cohesion evolution law (it is 1-funcH here) */
- Real funcH(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{ return 1-funcG(kappaD,epsCrackOnset,epsFracture,neverDamage); }
/*! Damage evolution law */
static Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) {
if(kappaD<epsCrackOnset || neverDamage) return 0;
@@ -198,8 +196,6 @@
MetaBody* rootBody;
//! aplly calculated force on both particles (applied in the inverse sense on B)
void applyForce(const Vector3r&, const body_id_t&, const body_id_t&);
- /*! Cohesion evolution law (it is 1-funcH here) */
- Real funcH(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{ return 1-funcG(kappaD,epsCrackOnset,epsFracture,neverDamage); }
/*! Damage evolution law */
Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{ return ef2_Spheres_Brefcom_BrefcomLaw::funcG(kappaD,epsCrackOnset,epsFracture,neverDamage); }
Modified: trunk/extra/usct/UniaxialStrainControlledTest.cpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.cpp 2009-03-31 23:10:43 UTC (rev 1742)
+++ trunk/extra/usct/UniaxialStrainControlledTest.cpp 2009-04-02 09:28:06 UTC (rev 1743)
@@ -48,6 +48,11 @@
if(isnan(originalLength)) LOG_FATAL("Initial length is NaN!");
assert(originalLength>0 && !isnan(originalLength));
+ assert(!isnan(strainRate) || !isnan(absSpeed));
+ if(isnan(strainRate)){ strainRate=absSpeed/originalLength; }
+
+ initAccelTime_s=initAccelTime>=0 ? initAccelTime : Omega::instance().getTimeStep()*(-initAccelTime);
+
/* if we have default (<0) crossSectionArea, try to get it from root's AABB;
* this will not work if there are foreign bodies in the simulation,
* in which case you must give the value yourself as engine attribute.
@@ -82,7 +87,10 @@
//nothing to do
if(posIds.size()==0 || negIds.size()==0) return;
// linearly increase strain to the desired value
- if(abs(currentStrainRate)<abs(strainRate))currentStrainRate+=strainRate*.005; else currentStrainRate=strainRate;
+ if(abs(currentStrainRate)<abs(strainRate)){
+ Real t=Omega::instance().getSimulationTime();
+ currentStrainRate=(t/initAccelTime_s)*strainRate;
+ } else currentStrainRate=strainRate;
// how much do we move (in total, symmetry handled below)
Real dAX=currentStrainRate*originalLength*Omega::instance().getTimeStep();
if(!isnan(stopStrain)){
@@ -114,12 +122,6 @@
if(Omega::instance().getCurrentIteration()%10==0) {
computeAxialForce(rootBody);
- #if 0
- vector<Real> widths;
- pushTransStrainSensors(rootBody,widths);
- assert(widths.size()==originalWidths.size());
- for(size_t i=0; i<widths.size(); i++) avgTransStrain+=(widths[i]/originalWidths[i]-1); avgTransStrain/=widths.size();
- #endif
avgStress=(sumPosForces+sumNegForces)/(2*crossSectionArea); // average nominal stress
if(!recordFile.empty() && recStream.good()) recStream<<Omega::instance().getCurrentIteration()<<" "<<strain<<" "<<avgStress<<endl; // <<" "<<avgTransStrain<<endl;
}
@@ -285,108 +287,3 @@
-
-
-
-
-
-/******************* DEPRECATED *********************/
-
-
-
-
-
-#if 0
-/* Initialize UniaxialStrainSensorPusher so that subscribed bodies and forces are consistent.
- * Apply small forces on those bodies, cap their velocity and reset orientation.
- * Return vector of widths in the direction of sensor pairs.
- */
-void UniaxialStrainer::pushTransStrainSensors(MetaBody* rb, vector<Real>& widths){
- if(transStrainSensors.size()==0) return;
- if(!sensorsPusher){
- int count=0;
- FOREACH(const shared_ptr<Engine>& e,rb->engines){
- if(e->getClassName()=="UniaxialStrainSensorPusher"){ count++; sensorsPusher=static_pointer_cast<UniaxialStrainSensorPusher>(e); }
- }
- if(count>1) LOG_ERROR("Multiple UniaxialStrainSensorPusher's found, using the last one for transversal strain sensors!");
- if(count<1) { LOG_ERROR("No UniaxialStrainSensorPusher found, transversal strain sensors will not work!"); return; }
- sensorsPusher->subscribedBodies.clear();
- FOREACH(body_id_t id, transStrainSensors) sensorsPusher->subscribedBodies.push_back(id);
- sensorsPusher->forces.resize(transStrainSensors.size());
- //TRVAR3(transStrainSensors.size(),sensorsPusher->subscribedBodies.size(),sensorsPusher->forces.size());
- }
- assert((sensorsPusher->subscribedBodies.size()==transStrainSensors.size()) && (sensorsPusher->subscribedBodies.size()==sensorsPusher->forces.size()));
- Real forceMagnitude=.001*abs(avgStress)*transStrainSensorArea;
- Real maxVelocity=2*abs(strainRate)*originalLength; // move at max 5 × faster than strained ends
- /* reset orientation to identity and limit velocity */
- FOREACH(body_id_t id, transStrainSensors){
- const shared_ptr<Body>& b=Body::byId(id); const shared_ptr<ParticleParameters>& pp=YADE_PTR_CAST<ParticleParameters>(b->physicalParameters);
- pp->se3.orientation=Quaternionr::IDENTITY;
- if(pp->velocity.SquaredLength()>pow(maxVelocity,2)){ pp->velocity.Normalize(); pp->velocity*=maxVelocity; }
- }
- /* calcuate and store force that will be applied in UniaxialStrainSensorPusher */
- widths.clear();
- for(int i=1; i<=(transStrainSensors.size()==2?1:2); i++){
- int transAxis=(axis+i)%3, perpTransAxis=(i==1?(axis+2)%3:/* i==2 */ (axis+1)%3);
- Vector3r F; F[axis]=0; F[perpTransAxis]=0; F[transAxis]=forceMagnitude;
- body_id_t n1=2*(i-1), n2=2*(i-1)+1;
- sensorsPusher->forces[n1]=+F; sensorsPusher->forces[n2]=-F;
- const shared_ptr<Body>& lo=Body::byId(transStrainSensors[n1]), hi=Body::byId(transStrainSensors[n2]);
- Real wd=hi->physicalParameters->se3.position[transAxis]-lo->physicalParameters->se3.position[transAxis]-static_pointer_cast<Box>(hi->geometricalModel)->extents[transAxis]-static_pointer_cast<Box>(lo->geometricalModel)->extents[transAxis];
- // negative width? Apply no more force, reset velocity to 0
- if(wd<=0) {
- /* doesn't work... Why? */
- //LOG_ERROR("Width is negative, resetting forces and velocities");
- sensorsPusher->forces[n1]=sensorsPusher->forces[n2]=Vector3r::ZERO;
- YADE_PTR_CAST<ParticleParameters>(lo->physicalParameters)->velocity=Vector3r::ZERO;
- YADE_PTR_CAST<ParticleParameters>(hi->physicalParameters)->velocity=Vector3r::ZERO;
- wd=0;
- }
- widths.push_back(wd);
- }
-}
-
-void UniaxialStrainer::setupTransStrainSensors(){
- assert(transStrainSensors.size()==0 || transStrainSensors.size()==2 || transStrainSensors.size()==4);
- if(transStrainSensors.size()==0) return;
- assert(Omega::instance().getRootBody()->boundingVolume);
- shared_ptr<AABB> rbAABB=dynamic_pointer_cast<AABB>(Omega::instance().getRootBody()->boundingVolume);
- assert(rbAABB);
- TRVAR2(rbAABB->center,rbAABB->halfSize);
- transStrainSensorArea=0;
- int numCouples=(transStrainSensors.size()==2?1:2);
- originalWidths.clear();
- for(int i=1; i<=numCouples; i++){
- int transAxis=(axis+i)%3, perpTransAxis=(i==1?(axis+2)%3:/* i==2 */ (axis+1)%3);
- TRVAR3(axis,transAxis,perpTransAxis);
- originalWidths.push_back(2*rbAABB->halfSize[transAxis]);
- body_id_t sensId[]={transStrainSensors[2*(i-1)],transStrainSensors[2*(i-1)+1]};
- // do the same on either side, only positioning will be different
- FOREACH(body_id_t id, sensId){
- shared_ptr<Body> b=Body::byId(id);
- shared_ptr<RigidBodyParameters> rbp=dynamic_pointer_cast<RigidBodyParameters>(b->physicalParameters);
- shared_ptr<Box> box=dynamic_pointer_cast<Box>(b->geometricalModel);
- shared_ptr<InteractingBox> iBox=dynamic_pointer_cast<InteractingBox>(b->interactingGeometry);
- assert(rbp && box && iBox);
- LOG_DEBUG("Setting up transversal strain sensor, #"<<id);
- // change box size: axis,transAxis: length_specimen_along_axis)/10; perpTransAxis: width_specimen*1.5
- Vector3r ext;
- ext[axis]=ext[transAxis]=.1*rbAABB->halfSize[axis]; ext[perpTransAxis]=1.5*rbAABB->halfSize[perpTransAxis];
- box->extents=iBox->extents=ext;
- // reset orientation
- rbp->se3.orientation=Quaternionr::IDENTITY;
- // set isDynamic==True, GeometricalMode::wire=true;
- b->isDynamic=true; box->wire=true;
- // set position so that it touches AABB of rootBody
- int sign=(id==sensId[0]?-1:1); // first id goes underneath (pushed in the dir of +transAxis), the other one goes up
- rbp->se3.position[axis]=rbAABB->center[axis];
- rbp->se3.position[perpTransAxis]=rbAABB->center[perpTransAxis];
- rbp->se3.position[transAxis]=rbAABB->center[transAxis]+sign*rbAABB->halfSize[transAxis]+sign*iBox->extents[transAxis];
- TRVAR2(box->extents,rbp->se3.position);
- transStrainSensorArea+=(box->extents[axis]*rbAABB->halfSize[perpTransAxis]);
- }
- }
- transStrainSensorArea/=transStrainSensors.size();
-}
-#endif
-
Modified: trunk/extra/usct/UniaxialStrainControlledTest.hpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.hpp 2009-03-31 23:10:43 UTC (rev 1742)
+++ trunk/extra/usct/UniaxialStrainControlledTest.hpp 2009-04-02 09:28:06 UTC (rev 1743)
@@ -37,30 +37,7 @@
};
REGISTER_SERIALIZABLE(USCTGen);
-#if 0
-/* This class applies force on transversal strain sensors from UniaxialStrainer.
- */
-class UniaxialStrainSensorPusher: public DeusExMachina{
- public:
- UniaxialStrainSensorPusher(){ Shop::Bex::initCache(); }
- ~UniaxialStrainSensorPusher(){}
- vector<Vector3r> forces;
- virtual void applyCondition(MetaBody* mb){
- assert(subscribedBodies.size()==forces.size());
- for(size_t i=0; i<subscribedBodies.size(); i++) Shop::Bex::force(subscribedBodies[i])+=forces[i];
- }
- void registerAttributes(){
- DeusExMachina::registerAttributes();
- REGISTER_ATTRIBUTE(forces);
- }
- REGISTER_CLASS_NAME(UniaxialStrainSensorPusher);
- REGISTER_BASE_CLASS_NAME(DeusExMachina);
- //DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(UniaxialStrainSensorPusher);
-#endif
-
/*! Axial displacing two groups of bodies in the opposite direction with given strain rate.
*
* Takes two groups of body IDs (in posIds and negIds) and displaces them at each timestep in the direction given by axis∈{0,1,2} (for axes x,y,z respectively). These bodies automatically have Body::isDynamic==false.
@@ -79,12 +56,14 @@
void init();
public:
virtual bool isActivated(){return active;}
+ //! strain rate, starting at 0, linearly raising to strainRate
Real strainRate,currentStrainRate;
+ //! alternatively, absolute speed of boundary motion can be specified; this is effective only at the beginning and if strainRate is not set; changing absSpeed directly during simulation wil have no effect.
+ Real absSpeed;
//! strain at which we will pause simulation; inactive (nan) by default; must be reached from below (in absolute value)
Real stopStrain;
//! distance of reference bodies in the direction of axis before straining started
Real originalLength;
- vector<Real> originalWidths;
//! invert the sense of straining (sharply, without transition) one this value of strain is reached. Not effective if 0.
Real limitStrain;
//! Flag whether the sense of straining has already been reversed
@@ -103,6 +82,8 @@
bool active;
//! Number of iterations that will pass without straining activity after stopStrain has been reached (default: 0)
long idleIterations;
+ //! Time for strain reaching the requested value (linear interpolation). If negative, the time is dt*(-initAccelTime), where dt is the timestep at the first iteration.
+ Real initAccelTime, initAccelTime_s /* value always in s, computed from initAccelTime */;
/** bodies on which straining will be applied (on the positive and negative side of axis) */
vector<body_id_t> posIds, negIds;
@@ -129,10 +110,12 @@
Real strain, avgStress;
virtual void applyCondition(MetaBody* rootBody);
- UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0; originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1; needsInit=true; /* sensorsPusher=shared_ptr<UniaxialStrainSensorPusher>(); */ recordFile=""; strain=avgStress=/*avgTransStrain=*/0; blockRotations=false; blockDisplacements=false; stopStrain=numeric_limits<Real>::quiet_NaN(); active=true; idleIterations=0; };
+ UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0; originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1; needsInit=true; /* sensorsPusher=shared_ptr<UniaxialStrainSensorPusher>(); */ recordFile=""; strain=avgStress=/*avgTransStrain=*/0; blockRotations=false; blockDisplacements=false; absSpeed=strainRate=stopStrain=numeric_limits<Real>::quiet_NaN(); active=true; idleIterations=0; };
virtual ~UniaxialStrainer(){};
REGISTER_ATTRIBUTES(DeusExMachina,
(strainRate)
+ (absSpeed)
+ (initAccelTime)
(stopStrain)
(active)
(idleIterations)
@@ -142,7 +125,6 @@
(posIds)
(negIds)
(originalLength)
- (originalWidths)
(limitStrain)
(notYetReversed)
(crossSectionArea)
Modified: trunk/gui/py/plot.py
===================================================================
--- trunk/gui/py/plot.py 2009-03-31 23:10:43 UTC (rev 1742)
+++ trunk/gui/py/plot.py 2009-04-02 09:28:06 UTC (rev 1743)
@@ -47,6 +47,11 @@
for attr in ['virtPeriod','realPeriod','iterPeriod']:
if(plotDataCollector[attr]>0): plotDataCollector[attr]=2*plotDataCollector[attr]
+def splitData():
+ "Make all plots discontinuous at this point (adds nan's to all data fields)"
+ addData({})
+
+
def reverseData():
for k in data: data[k].reverse()
@@ -84,6 +89,7 @@
plotLines[p]=pylab.plot(*sum([[data[p],data[d[0]],d[1]] for d in plots_p],[]))
pylab.legend([_p[0] for _p in plots_p])
pylab.xlabel(p)
+ if 'title' in O.tags.keys(): pylab.title(O.tags['title'])
pylab.show()
updatePeriod=0
def periodicUpdate(period):
@@ -141,6 +147,7 @@
if term in ['wxt','x11']: fPlot.write("set term %s %d persist\n"%(term,i))
else: fPlot.write("set term %s; set output '%s.%d.%s'\n"%(term,baseNameNoPath,i,extension))
fPlot.write("set xlabel '%s'\n"%p)
+ fPlot.write("set datafile missing 'nan'\n"%p)
if title: fPlot.write("set title '%s'\n"%title)
fPlot.write("plot "+",".join([" %s using %d:%d title '%s(%s)' with lines"%(dataFile,vars.index(p)+1,vars.index(pp[0])+1,pp[0],p) for pp in plots_p])+"\n")
i+=1
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-03-31 23:10:43 UTC (rev 1742)
+++ trunk/gui/py/yadeControl.cpp 2009-04-02 09:28:06 UTC (rev 1743)
@@ -445,12 +445,22 @@
shared_ptr<MetaEngine> ee=dynamic_pointer_cast<MetaEngine>(e);
FOREACH(const shared_ptr<EngineUnit>& f, ee->functorArguments){
if(!f->label.empty()){
- PyGILState_STATE gstate; gstate = PyGILState_Ensure();
- PyRun_SimpleString(("__builtins__."+f->label+"=Omega().labeledEngine('"+f->label+"')").c_str());
- PyGILState_Release(gstate);
+ PyGILState_STATE gstate; gstate = PyGILState_Ensure(); PyRun_SimpleString(("__builtins__."+f->label+"=Omega().labeledEngine('"+f->label+"')").c_str()); PyGILState_Release(gstate);
}
}
}
+ if(isChildClassOf(e->getClassName(),"InteractionDispatchers")){
+ shared_ptr<InteractionDispatchers> ee=dynamic_pointer_cast<InteractionDispatchers>(e);
+ list<shared_ptr<EngineUnit> > eus;
+ FOREACH(const shared_ptr<EngineUnit>& eu,ee->geomDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<EngineUnit>& eu,ee->physDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<EngineUnit>& eu,ee->constLawDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<EngineUnit>& eu,eus){
+ if(!eu->label.empty()){
+ PyGILState_STATE gstate; gstate = PyGILState_Ensure(); PyRun_SimpleString(("__builtins__."+eu->label+"=Omega().labeledEngine('"+eu->label+"')").c_str()); PyGILState_Release(gstate);
+ }
+ }
+ }
}
}
@@ -557,6 +567,16 @@
if(eu->label==label) return python::object(pyEngineUnit(eu));
}
}
+ shared_ptr<InteractionDispatchers> ee=dynamic_pointer_cast<InteractionDispatchers>(eng);
+ if(ee){
+ list<shared_ptr<EngineUnit> > eus;
+ FOREACH(const shared_ptr<EngineUnit>& eu,ee->geomDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<EngineUnit>& eu,ee->physDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<EngineUnit>& eu,ee->constLawDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<EngineUnit>& eu,eus){
+ if(eu->label==label) return python::object(pyEngineUnit(eu));
+ }
+ }
}
throw std::invalid_argument(string("No engine labeled `")+label+"'");
}
@@ -607,7 +627,7 @@
BASIC_PY_PROXY_HEAD(pyFileGenerator,FileGenerator)
void generate(string outFile){ensureAcc(); proxee->setFileName(outFile); proxee->setSerializationLibrary("XMLFormatManager"); bool ret=proxee->generateAndSave(); LOG_INFO((ret?"SUCCESS:\n":"FAILURE:\n")<<proxee->message); if(ret==false) throw runtime_error("Generator reported error: "+proxee->message); };
- void load(){ ensureAcc(); char tmpnam_str [L_tmpnam]; tmpnam(tmpnam_str); string xml(tmpnam_str+string(".xml.bz2")); LOG_DEBUG("Using temp file "<<xml); this->generate(xml); pyOmega().load(xml); }
+ void load(){ ensureAcc(); char tmpnam_str [L_tmpnam]; char* result=tmpnam(tmpnam_str); if(result!=tmpnam_str) throw runtime_error(__FILE__ ": tmpnam(char*) failed!"); string xml(tmpnam_str+string(".xml.bz2")); LOG_DEBUG("Using temp file "<<xml); this->generate(xml); pyOmega().load(xml); }
BASIC_PY_PROXY_TAIL;
class pySTLImporter : public STLImporter {
