yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02521
[Branch ~yade-dev/yade/trunk] Rev 1861: 1. Enhance PeriTriaxController for controlling strain/stress in periodic simulations; stress and ...
------------------------------------------------------------
revno: 1861
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Sun 2009-12-06 18:18:16 +0100
message:
1. Enhance PeriTriaxController for controlling strain/stress in periodic simulations; stress and strains appear to be physically correct now (however, see https://bugs.launchpad.net/yade/+bug/493102)
2. Add examples/concrete/periodic.py, showing uniaxial/biaxial periodic loading. Results are awesome.
3. Fix stiffness computation for Cpm.
4. Add (unpolished) pack.randomPeriPack for generating periodic packings.
5. Add some sphinx doc files (no contents, just skeleton).
added:
core/GLConfig.hpp
doc/sphinx/modules.rst
doc/sphinx/yade.pack.rst
doc/sphinx/yade.utils.rst
examples/concrete/periodic.py
modified:
core/main/main.py.in
examples/concrete/README
pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp
pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp
pkg/dem/meta/ConcretePM.cpp
pkg/dem/meta/ConcretePM.hpp
py/pack.py
py/utils.py
scripts/default-test.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.
=== added file 'core/GLConfig.hpp'
--- core/GLConfig.hpp 1970-01-01 00:00:00 +0000
+++ core/GLConfig.hpp 2009-12-06 17:18:16 +0000
@@ -0,0 +1,39 @@
+// code not yet for use (6/12/2009); if long here uselessly, delete.
+//
+
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#include<yade/core/Serializable.hpp>
+/*! Storage for general 3d view settings.
+
+Is saved along with simulation and passed to every call to render(...).
+Contains more or less what used to be inside OpenGLRenderingEngine.
+
+*/
+class GLConfig: public Serializable{
+
+ Vector3r lightPos,bgColor;
+ body_id_t currSel;
+ bool dof,id,bbox,geom,wire,intrGeom,intrPhys;
+ int mask;
+ bool scaleDisplacements,scaleRotations;
+ Vector3r displacementScale; Real rotationScale;
+ vector<Se3r> clipPlaneSe3;
+ vector<int> clipPlaneActive; // should be bool, but serialization doesn't handle vector<bool>
+ const int clipPlaneNum;
+
+ // not saved
+ Vector3r highlightEmission0;
+ Vector3r highlightEmission1;
+
+ // normalized saw signal with given periodicity, with values ∈ 〈0,1〉 */
+ Real normSaw(Real t, Real period){ Real xi=(t-period*((int)(t/period)))/period; /* normalized value, (0-1〉 */ return (xi<.5?2*xi:2-2*xi); }
+ Real normSquare(Real t, Real period){ Real xi=(t-period*((int)(t/period)))/period; /* normalized value, (0-1〉 */ return (xi<.5?0:1); }
+
+ //! wrap number to interval x0…x1
+ Real wrapCell(const Real x, const Real x0, const Real x1);
+ //! wrap point to inside Scene's cell (identity if !Scene::isPeriodic)
+ Vector3r wrapCellPt(const Vector3r& pt, Scene* rb);
+ void drawPeriodicCell(Scene*);
+
+ REGISTER_ATTRIBUTES(Serializable,(dof)(id)(bbox)(geom)(wire)(intrGeom)(intrPhys)(mask)(scaleDisplacements)(scaleRotations)(displacementScale)(rotationScale)(clipPlaneSe3)(clipPlaneActive));
+};
=== modified file 'core/main/main.py.in'
--- core/main/main.py.in 2009-12-05 23:12:46 +0000
+++ core/main/main.py.in 2009-12-06 17:18:16 +0000
@@ -32,6 +32,7 @@
import yade.wrapper
import yade.log
import yade.system
+import yade.plot
# continue option processing
@@ -82,7 +83,7 @@
ipshell=IPShellEmbed(
argv=[],
#exit_msg='Bye.',
- banner='[[ ^L clears screen, ^U kills line.'+(' F12 controller, F11 3d view, F10 both, F9 generator.' if qtEnabled else '')+' ]]',
+ banner='[[ ^L clears screen, ^U kills line.'+(' F12 controller, F11 3d view, F10 both, F9 generator, F8 plot.' if qtEnabled else '')+' ]]',
rc_override={
'execfile':[prefix+'/lib/yade'+suffix+'/py/yade/ipython.py'],
'readline_parse_and_bind':[
@@ -93,6 +94,7 @@
'"\e[23~": "\C-Uyade.qt.View();\C-M"', # F11
'"\e[21~": "\C-Uyade.qt.Controller(), yade.qt.View();\C-M"', # F10
'"\e[20~": "\C-Uyade.qt.Generator();\C-M"', #F9
+ '"\e[19~": "\C-Uyade.plot.plot();\C-M"', #F9
] if qtEnabled else [])
})
ipshell()
=== added file 'doc/sphinx/modules.rst'
--- doc/sphinx/modules.rst 1970-01-01 00:00:00 +0000
+++ doc/sphinx/modules.rst 2009-12-06 17:18:16 +0000
@@ -0,0 +1,16 @@
+Yade modules
+=============
+
+.. toctree::
+
+ yade.pack.rst
+
+ yade.utils.rst
+
+
+.. * :ref:`yade.utils module <yade.utils>`
+ * :ref:`yade.pack modue <yade.pack>`
+ * :ref:`yade.log`
+ * :ref:`yade.timing`
+ * :ref:`yade.post2d`
+ * :ref:`yade.plot`
=== added file 'doc/sphinx/yade.pack.rst'
--- doc/sphinx/yade.pack.rst 1970-01-01 00:00:00 +0000
+++ doc/sphinx/yade.pack.rst 2009-12-06 17:18:16 +0000
@@ -0,0 +1,10 @@
+.. _yade.pack
+
+yade.pack
+=========
+
+.. automodule:: yade.pack
+ :members:
+ :undoc-members:
+
+
=== added file 'doc/sphinx/yade.utils.rst'
--- doc/sphinx/yade.utils.rst 1970-01-01 00:00:00 +0000
+++ doc/sphinx/yade.utils.rst 2009-12-06 17:18:16 +0000
@@ -0,0 +1,11 @@
+.. _yade.utils
+
+yade.utils
+==========
+
+.. automodule:: yade.utils
+
+
+.. :members:
+ #:undoc-members:
+
=== modified file 'examples/concrete/README'
--- examples/concrete/README 2009-08-22 18:21:25 +0000
+++ examples/concrete/README 2009-12-06 17:18:16 +0000
@@ -8,3 +8,6 @@
Type yade.plot.plot() to see strain/stress plot during the simulation.
See comments at the beginning of uniax.py for more information.
+
+The periodic.py is uniax.py, but running on a periodic sample.
+It is only experimental at this moment.
=== added file 'examples/concrete/periodic.py'
--- examples/concrete/periodic.py 1970-01-01 00:00:00 +0000
+++ examples/concrete/periodic.py 2009-12-06 17:18:16 +0000
@@ -0,0 +1,195 @@
+# -*- encoding=utf-8 -*-
+from __future__ import division
+
+from yade import utils,plot,pack
+import time, sys, os, copy
+
+"""
+A fairly complex script performing uniaxial tension-compression test on hyperboloid-shaped specimen.
+
+Most parameters of the model (and of the setup) can be read from table using yade-multi.
+
+After the simulation setup, tension loading is run and stresses are periodically saved for plotting
+as well as checked for getting below the maximum value so far. This indicates failure (see stopIfDamaged
+function). After failure in tension, the original setup is loaded anew and the sense of loading reversed.
+After failure in compression, strain-stress curves are saved via plot.saveGnuplot and we exit,
+giving some useful information like peak stresses in tension/compression.
+
+Running this script for the first time can take long time, as the specimen is prepared using triaxial
+compression. Next time, however, an attempt is made to load previously-generated packing
+(from /tmp/triaxPackCache.sqlite) and this expensive procedure is avoided.
+
+The specimen length can be specified, its diameter is half of the length and skirt of the hyperboloid is
+4/5 of the width.
+
+The particle size is constant and can be specified using the sphereRadius parameter.
+
+The 3d display has displacement scaling applied, so that the fracture looks more spectacular. The scale
+is 1000 for tension and 100 for compression.
+
+"""
+
+
+
+# default parameters or from table
+utils.readParamsFromTable(noTableOk=True, # unknownOk=True,
+ young=24e9,
+ poisson=.2,
+ G_over_E=.20,
+
+ sigmaT=3.5e6,
+ frictionAngle=atan(0.8),
+ epsCrackOnset=1e-4,
+ relDuctility=30,
+
+ intRadius=1.5,
+ dtSafety=.8,
+ damping=0.4,
+ strainRateTension=10,
+ strainRateCompression=50,
+ # 1=tension, 2=compression (ANDed; 3=both)
+ doModes=3,
+ biaxial=True,
+
+ # isotropic confinement (should be negative)
+ isoPrestress=0,
+)
+
+if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']
+
+packingFile='periCube.pickle'
+# got periodic packing? Memoization not (yet) supported, so just generate it if there is not the right file
+# Save and reuse next time.
+if not os.path.exists(packingFile):
+ sp=pack.randomPeriPack(radius=.05e-3,rRelFuzz=0.,initSize=Vector3().ONE*1.5e-3)
+ dd=dict(cell=(sp.cellSize[0],sp.cellSize[1],sp.cellSize[2]),spheres=sp.toList_pointsAsTuples())
+ import cPickle as pickle
+ pickle.dump(dd,open(packingFile,'w'))
+
+#
+# load the packing (again);
+#
+import cPickle as pickle
+concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800))
+sphDict=pickle.load(open(packingFile))
+from yade import pack
+sp=pack.SpherePack()
+sp.fromList(sphDict['spheres'])
+sp.cellSize=sphDict['cell']
+
+import numpy
+avgRadius=numpy.average([r for c,r in sp])
+O.bodies.append([utils.sphere(c,r,color=utils.randomColor()) for c,r in sp])
+O.periodicCell=(Vector3(0,0,0),sp.cellSize)
+axis=2
+ax1=(axis+1)%3
+ax2=(axis+2)%3
+O.dt=dtSafety*utils.PWaveTimeStep()
+
+import yade.log
+#yade.log.setLevel('PeriTriaxController',yade.log.TRACE)
+import yade.plot as yp
+
+O.engines=[
+ BexResetter(),
+ BoundDispatcher([InteractingSphere2AABB(aabbEnlargeFactor=intRadius,label='is2aabb'),MetaInteractingGeometry2AABB()]),
+ InsertionSortCollider(sweepLength=.05*avgRadius,nBins=5,binCoeff=5),
+ InteractionDispatchers(
+ [Ig2_Sphere_Sphere_Dem3DofGeom(distFactor=intRadius,label='ss2d3dg')],
+ [Ip2_CpmMat_CpmMat_CpmPhys(sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E,isoPrestress=isoPrestress)],
+ [Law2_Dem3DofGeom_CpmPhys_Cpm()],
+ ),
+ NewtonsDampedLaw(damping=damping,label='newton'),
+ CpmStateUpdater(realPeriod=1,label='updater'),
+ #
+ #UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
+ #
+ PeriTriaxController(goal=[1,1,1],stressMask=( (7^(1<<axis | 1<<ax1)) if biaxial else (7^(1<<axis)) ),maxStrainRate=Vector3(1,1,1),label='strainer',reversedForces=False,globUpdate=2),
+ PeriodicPythonRunner(virtPeriod=1e-5/strainRateTension,realLim=2,command='addPlotData()',label='plotDataCollector'),
+ PeriodicPythonRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
+]
+#O.miscParams=[GLDrawCpmPhys(dmgLabel=False,colorStrain=False,epsNLabel=False,epsT=False,epsTAxes=False,normal=False,contactLine=True)]
+
+# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that includes confinement, though
+plot.plots={'eps':('sigma','sig1','sig2','|||','eps','eps1','eps2'),'eps_':('sigma','|||','relResid',),} #'sigma.25','sigma.50','sigma.75')}
+plot.maxDataLen=4000
+
+O.saveTmp('initial');
+
+global mode
+mode='tension' if doModes & 1 else 'compression'
+
+def initTest():
+ global mode
+ print "init"
+ if O.iter>0:
+ O.wait();
+ O.loadTmp('initial')
+ print "Reversing plot data"; plot.reverseData()
+ maxStrainRate=Vector3(1,1,1)
+ goal=Vector3(1,1,1);
+ if not biaxial: # uniaxial
+ maxStrainRate[axis]=abs(strainRateTension) if mode=='tension' else abs(strainRateCompression)
+ maxStrainRate[ax1]=maxStrainRate[ax2]=100*maxStrainRate[axis]
+ goal[axis]=1 if mode=='tension' else -1;
+ else:
+ maxStrainRate[axis]=abs(strainRateTension) if mode=='tension' else abs(strainRateCompression)
+ maxStrainRate[ax1]=maxStrainRate[axis]
+ maxStrainRate[ax2]=100*maxStrainRate[axis]
+ goal[axis]=1 if mode=='tension' else -1;
+ goal[ax1]=goal[axis]
+ strainer['maxStrainRate']=maxStrainRate
+ strainer['goal']=goal
+ try:
+ pass
+ #from yade import qt
+ #renderer=qt.Renderer()
+ #renderer['scaleDisplacements']=True
+ #renderer['displacementScale']=(1000,1000,1000) if mode=='tension' else (100,100,100)
+ except ImportError: pass
+ print "init done, will now run."
+ O.step(); O.step(); # to create initial contacts
+ # now reset the interaction radius and go ahead
+ ss2d3dg['distFactor']=-1.
+ is2aabb['aabbEnlargeFactor']=-1.
+ O.run()
+
+def stopIfDamaged():
+ global mode
+ if O.iter<2 or not plot.data.has_key('sigma'): return # do nothing at the very beginning
+ sigma,eps=plot.data['sigma'],plot.data['eps']
+ extremum=max(sigma) if (strainer['maxStrainRate']>0) else min(sigma)
+ # FIXME: only temporary, should be .5
+ minMaxRatio=0.5 if mode=='tension' else 0.5
+ if extremum==0: return
+ cellSize=O.periodicCell[1]-O.periodicCell[0]
+ print O.tags['id'],mode,strainer['strain'][axis],sigma[-1]
+ #print 'strain',strainer['strain'],'stress',strainer['stress']
+ import sys; sys.stdout.flush()
+ if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer['strain'][axis])>6e-3:
+ if mode=='tension' and doModes & 2: # only if compression is enabled
+ mode='compression'
+ #O.save('/tmp/uniax-tension.xml.bz2')
+ print "Damaged, switching to compression... "; O.pause()
+ # important! initTest must be launched in a separate thread;
+ # otherwise O.load would wait for the iteration to finish,
+ # but it would wait for initTest to return and deadlock would result
+ import thread; thread.start_new_thread(initTest,())
+ return
+ else:
+ print "Damaged, stopping."
+ ft,fc=max(sigma),min(sigma)
+ print 'Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft))
+ title=O.tags['description'] if 'description' in O.tags.keys() else O.tags['params']
+ print'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
+ print 'Bye.'
+ # O.pause()
+ sys.exit(0)
+
+def addPlotData():
+ yade.plot.addData(t=O.time,i=O.iter,eps=strainer['strain'][axis],eps_=strainer['strain'][axis],sigma=strainer['stress'][axis]+isoPrestress,eps1=strainer['strain'][ax1],eps2=strainer['strain'][ax2],sig1=strainer['stress'][ax1],sig2=strainer['stress'][ax2],relResid=updater['avgRelResidual'])
+
+initTest()
+# sleep forever if run by yade-multi, exit is called from stopIfDamaged
+if os.environ.has_key('PARAM_TABLE'): time.sleep(1e12)
+
=== modified file 'pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp 2009-12-05 23:12:46 +0000
+++ pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.cpp 2009-12-06 17:18:16 +0000
@@ -103,10 +103,10 @@
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];
+ Real f=(reversedForces?-1.:1.)*(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;
+ sumForce[i]+=Mathr::Sign(f*gsc->normal[i])*abs(f);
sumLength[i]+=absNormI*(gsc->refR1+gsc->refR2);
sumStiff[i]+=absNormI*nsi->kn+(1-absNormI)*nsi->ks;
}
@@ -143,8 +143,8 @@
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
+ Real maxGrow=refSize[axis]*maxStrainRate[axis]*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.);
@@ -160,11 +160,9 @@
// 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){
+ 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 {
@@ -175,9 +173,18 @@
if(doUpdate) LOG_DEBUG("Strain not OK; "<<abs(curr-goal[axis])<<">1e-6");
}
}
- // change cell size now
- scene->cellMax+=cellGrow;
- }
+ }
+ // update stress and strain
+ for(int axis=0; axis<3; axis++){
+ strain[axis]+=cellGrow[axis]/refSize[axis];
+ // take in account something like poisson's effect here…
+ //Real bogusPoisson=0.25; int ax1=(axis+1)%3,ax2=(axis+2)%3;
+ if(stiff[axis]>0) stress[axis]+=(cellGrow[axis]/refSize[axis])*(stiff[axis]/cellArea[axis]); //-bogusPoisson*(cellGrow[ax1]/refSize[ax1])*(stiff[ax1]/cellArea[ax1])-bogusPoisson*(cellGrow[ax2]/refSize[ax2])*(stiff[ax2]/cellArea[ax2]);
+ }
+ // change cell size now
+ scene->cellMax+=cellGrow;
+ strainRate=cellGrow/scene->dt;
+
if(allOk){
if(doUpdate || currUnbalanced<0) {
currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,scene);
=== modified file 'pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp 2009-12-05 23:12:46 +0000
+++ pkg/dem/Engine/StandAloneEngine/PeriIsoCompressor.hpp 2009-12-06 17:18:16 +0000
@@ -41,12 +41,15 @@
*/
class PeriTriaxController: public StandAloneEngine{
public:
+ //! For broken constitutive laws, normalForce and shearForce on interactions are in the reverse sense
+ //! see https://bugs.launchpad.net/yade/+bug/493102
+ bool reversedForces;
//! 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;
+ Vector3r maxStrainRate;
//! maximum unbalanced force (defaults to 1e-4)
Real maxUnbalanced;
//! Absolute stress tolerance (1e3)
@@ -69,6 +72,8 @@
Vector3r stress;
//! cell strain, updated at every step
Vector3r strain;
+ //! cell strain rate, updated at every step
+ Vector3r strainRate;
//! average stiffness, updated at every step (only every globUpdate steps recomputed from interactions)
Vector3r stiff;
//! current unbalanced force (updated every globUpdate)
@@ -78,8 +83,8 @@
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));
+ PeriTriaxController(): reversedForces(false),goal(Vector3r::ZERO),stressMask(0),maxStrainRate(Vector3r(1,1,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),strainRate(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO){}
+ REGISTER_ATTRIBUTES(StandAloneEngine,(reversedForces)(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(refSize)(stress)(strain)(strainRate)(stiff));
DECLARE_LOGGER;
REGISTER_CLASS_AND_BASE(PeriTriaxController,StandAloneEngine);
};
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2009-12-04 23:46:45 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2009-12-06 17:18:16 +0000
@@ -43,8 +43,8 @@
contPhys->epsCrackOnset=epsCrackOnset;
contPhys->epsFracture=relDuctility*epsCrackOnset;
// inherited from NormalShearInteracion, used in the timestepper
- contPhys->kn=contPhys->E*contPhys->crossSection;
- contPhys->ks=contPhys->G*contPhys->crossSection;
+
+ // contPhys->kn, contPhys->ks assigned in the constitutive law, as they depend on area of the contact as well
if(neverDamage) contPhys->neverDamage=true;
if(cohesiveThresholdIter<0 || (Omega::instance().getCurrentIteration()<cohesiveThresholdIter)) contPhys->isCohesive=true;
@@ -134,10 +134,16 @@
#endif
-void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody){
+void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* scene){
Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
+ // just the first time
+ if(I->isFresh(scene)){
+ BC->kn=BC->crossSection*BC->E/contGeom->refLength;
+ BC->ks=BC->crossSection*BC->G/contGeom->refLength;
+ }
+
// shorthands
Real& epsN(BC->epsN);
Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
@@ -198,20 +204,20 @@
// handle broken contacts
if(epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)){
if(isCohesive){
- const shared_ptr<Body>& body1=Body::byId(I->getId1(),rootBody), body2=Body::byId(I->getId2(),rootBody); assert(body1); assert(body2);
+ const shared_ptr<Body>& body1=Body::byId(I->getId1(),scene), body2=Body::byId(I->getId2(),scene); assert(body1); assert(body2);
const shared_ptr<CpmState>& st1=YADE_PTR_CAST<CpmState>(body1->state), st2=YADE_PTR_CAST<CpmState>(body2->state);
// nice article about openMP::critical vs. scoped locks: http://www.thinkingparallel.com/2006/08/21/scoped-locking-vs-critical-in-openmp-a-personal-shootout/
{ boost::mutex::scoped_lock lock(st1->updateMutex); st1->numBrokenCohesive+=1; st1->epsPlBroken+=epsPlSum; }
{ boost::mutex::scoped_lock lock(st2->updateMutex); st2->numBrokenCohesive+=1; st2->epsPlBroken+=epsPlSum; }
}
- rootBody->interactions->requestErase(I->getId1(),I->getId2());
+ scene->interactions->requestErase(I->getId1(),I->getId2());
return;
}
Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
Fs=sigmaT*crossSection; BC->shearForce=Fs;
- applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position, rootBody);
+ applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position, scene);
}
#ifdef YADE_OPENGL
@@ -345,12 +351,14 @@
/********************** CpmStateUpdater ****************************/
CREATE_LOGGER(CpmStateUpdater);
Real CpmStateUpdater::maxOmega=0.;
+Real CpmStateUpdater::avgRelResidual=0.;
-void CpmStateUpdater::update(Scene* _rootBody){
- Scene *rootBody=_rootBody?_rootBody:Omega::instance().getScene().get();
- vector<BodyStats> bodyStats; bodyStats.resize(rootBody->bodies->size());
+void CpmStateUpdater::update(Scene* _scene){
+ Scene *scene=_scene?_scene:Omega::instance().getScene().get();
+ vector<BodyStats> bodyStats; bodyStats.resize(scene->bodies->size());
assert(bodyStats[0].nCohLinks==0); // should be initialized by dfault ctor
- FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
+ avgRelResidual=0; Real nAvgRelResidual=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
shared_ptr<CpmPhys> phys=dynamic_pointer_cast<CpmPhys>(I->interactionPhysics);
if(!phys) continue;
@@ -372,8 +380,11 @@
bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum+=(1-phys->relResidualStrength); bodyStats[id1].epsPlSum+=phys->epsPlSum;
bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum+=(1-phys->relResidualStrength); bodyStats[id2].epsPlSum+=phys->epsPlSum;
maxOmega=max(maxOmega,phys->omega);
+
+ avgRelResidual+=phys->relResidualStrength;
+ nAvgRelResidual+=1;
}
- FOREACH(shared_ptr<Body> B, *rootBody->bodies){
+ FOREACH(shared_ptr<Body> B, *scene->bodies){
const body_id_t& id=B->getId();
// add damaged contacts that have already been deleted
CpmState* state=dynamic_cast<CpmState*>(B->state.get());
@@ -390,5 +401,7 @@
}
else { state->normDmg=0; state->normEpsPl=0;}
B->shape->diffuseColor=Vector3r(state->normDmg,1-state->normDmg,B->isDynamic?0:1);
+ nAvgRelResidual+=0.5*state->numBrokenCohesive; // add half or broken interactions, other body has the other half
}
+ avgRelResidual/=nAvgRelResidual;
}
=== modified file 'pkg/dem/meta/ConcretePM.hpp'
--- pkg/dem/meta/ConcretePM.hpp 2009-12-04 23:07:34 +0000
+++ pkg/dem/meta/ConcretePM.hpp 2009-12-06 17:18:16 +0000
@@ -334,13 +334,15 @@
class CpmStateUpdater: public PeriodicEngine {
struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Vector3r sigma, tau; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), sigma(Vector3r::ZERO), tau(Vector3r::ZERO) {} };
public:
+ //! average residual strength
+ static Real avgRelResidual;
//! maximum damage over all contacts
static Real maxOmega;
CpmStateUpdater(){maxOmega=0; /* run at the very beginning */ initRun=true;}
virtual void action(Scene* rb){ update(rb); }
static void update(Scene* rb=NULL);
DECLARE_LOGGER;
- REGISTER_ATTRIBUTES(PeriodicEngine,(maxOmega));
+ REGISTER_ATTRIBUTES(PeriodicEngine,(maxOmega)(avgRelResidual));
REGISTER_CLASS_AND_BASE(CpmStateUpdater,PeriodicEngine);
};
REGISTER_SERIALIZABLE(CpmStateUpdater);
=== modified file 'py/pack.py'
--- py/pack.py 2009-12-05 23:12:46 +0000
+++ py/pack.py 2009-12-06 17:18:16 +0000
@@ -354,4 +354,35 @@
import warnings; warnings.warn("pack.triaxialPack was renamed to pack.randomDensePack, update your code!",DeprecationWarning,stacklevel=2);
return randomDensePack(*args,**kw)
+def randomPeriPack(radius,rRelFuzz,initSize):
+ """Generate periodic dense packing. EXPERIMENTAL, you at your own risk.
+
+ A cell of initSize is stuffed with as many spheres as possible (ignore the warning from SpherePack::makeCloud about
+ not being able to add any more spheres), then we run periodic compression with PeriIsoCompressor, just like with
+ randomDensePack.
+
+ :param radius: mean sphere radius
+ :param rRelFuzz: relative fuzz of sphere radius (equal distribution); see the same param for randomDensePack.
+ :param initSize: initial size of the periodic cell.
+
+ :return: SpherePack object, which also contains periodicity information.
+
+ :todo: memoization in db; what criteria??
+
+ """
+ from math import pi
+ O.switchScene(); O.resetThisScene()
+ sp=SpherePack()
+ O.periodicCell=((0,0,0),Vector3(initSize))
+ sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,rRelFuzz,int(initSize[0]*initSize[1]*initSize[2]/((4/3.)*pi*radius**3)),True)
+ O.engines=[BexResetter(),BoundDispatcher([InteractingSphere2AABB()]),InsertionSortCollider(nBins=2,sweepLength=.05*radius),InteractionDispatchers([Ig2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=20,keepProportions=True),NewtonsDampedLaw(damping=.8)]
+ O.materials.append(GranularMat(young=30e9,frictionAngle=.1,poisson=.3,density=1e3))
+ for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
+ O.dt=utils.PWaveTimeStep()
+ O.timingEnabled=True
+ O.run(); O.wait()
+ ret=SpherePack()
+ ret.fromSimulation()
+ O.switchScene()
+ return ret
=== modified file 'py/utils.py'
--- py/utils.py 2009-12-04 23:27:26 +0000
+++ py/utils.py 2009-12-06 17:18:16 +0000
@@ -62,7 +62,7 @@
from math import sqrt
return radius/sqrt(young/density)
-def randomColor(): return [random.random(),random.random(),random.random()]
+def randomColor(): return Vector3(random.random(),random.random(),random.random())
def typedEngine(name): return [e for e in Omega().engines if e.name==name][0]
@@ -162,7 +162,6 @@
See utils.sphere's documentation for meaning of other parameters."""
b=Body()
- if not color: color=randomColor()
center=inscribedCircleCenter(vertices[0],vertices[1],vertices[2])
vertices=Vector3(vertices[0])-center,Vector3(vertices[1])-center,Vector3(vertices[2])-center
b.shape=Shape('InteractingFacet',diffuseColor=color if color else randomColor(),wire=wire,highlight=highlight)
=== modified file 'scripts/default-test.py'
--- scripts/default-test.py 2009-08-04 16:25:28 +0000
+++ scripts/default-test.py 2009-12-06 17:18:16 +0000
@@ -6,29 +6,36 @@
# using the default SMTP settings (sendmail?) on your system
#
import os,time,sys
-import yade.runtime
+from yade import *
+import yade.runtime,yade.system,yade.config
simulFile='/tmp/yade-test-%d.xml'%(os.getpid()) # generated simulations here
pyCmdFile='/tmp/yade-test-%d.py'%(os.getpid()) # generated script here
msgFile='/tmp/yade-test-%d.msg'%(os.getpid()) # write message here
runSimul="""
# generated file
+from yade import *
simulFile='%s'; msgFile='%s'; nIter=%d;
import time
try:
- o=Omega(); o.load(simulFile)
- o.run(10); o.wait() # run first 10 iterations
- start=time.time(); o.run(nIter); o.wait(); finish=time.time() # run nIter iterations, wait to finish, measure elapsed time
+ O.load(simulFile)
+ O.run(10); O.wait() # run first 10 iterations
+ start=time.time(); O.run(nIter); O.wait(); finish=time.time() # run nIter iterations, wait to finish, measure elapsed time
speed=nIter/(finish-start); open(msgFile,'w').write('%%g iter/sec'%%speed)
except:
import sys, traceback
traceback.print_exc()
sys.exit(1)
+print 'main: Yade: normal exit.'
+O.exitNoBacktrace()
quit()
"""%(simulFile,msgFile,100)
runGenerator="""
#generated file
-FileGenerator('%%s'%%s).generate('%s')
+from yade import *
+%%s(%%s).generate('%s')
+print 'main: Yade: normal exit.'
+O.exitNoBacktrace()
quit()
"""%(simulFile)
@@ -37,7 +44,7 @@
import subprocess,os,os.path,yade.runtime
f=open(pyCmdFile,'w'); f.write(cmd); f.close();
if os.path.exists(msgFile): os.remove(msgFile)
- p=subprocess.Popen([yade.runtime.executable,'-N','PythonUI','--','-n','-s',pyCmdFile],stdout=subprocess.PIPE,stderr=subprocess.STDOUT)
+ p=subprocess.Popen([sys.executable,pyCmdFile],stdout=subprocess.PIPE,stderr=subprocess.STDOUT,env=dict(os.environ,**{'PYTHONPATH':os.path.join(yade.config.prefix,'lib','yade'+yade.config.suffix,'py'),'DISPLAY':''}))
pout=p.communicate()[0]
retval=p.wait()
if not quiet: print pout
@@ -56,9 +63,8 @@
genParams={
#'USCTGen':{'spheresFile':'examples/small.sdec.xyz'}
}
-o=Omega()
-for pp in o.childClasses('FileGenerator'):
+for pp in yade.system.childClasses('FileGenerator'):
if pp in broken:
summary.append(pp,'skipped (broken)','');
params='' if pp not in genParams else (","+",".join(["%s=%s"%(k,repr(genParams[pp][k])) for k in genParams[pp]]))
