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 {
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
_______________________________________________
yade-dev mailing list
[email protected]
https://lists.berlios.de/mailman/listinfo/yade-dev
