yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06495
[Branch ~yade-dev/yade/trunk] Rev 2603: 1. Fixed problem in postLoad (regression tests should work now)
------------------------------------------------------------
revno: 2603
committer: Klaus Thoeni <klaus.thoeni@xxxxxxxxx>
branch nick: yade
timestamp: Fri 2010-12-10 15:42:10 +1300
message:
1. Fixed problem in postLoad (regression tests should work now)
2. Add some examples for testing the new WirePM with 2 particles
added:
scripts/test/WIreMatPM/
scripts/test/WIreMatPM/net-2part-displ-unloading.py
scripts/test/WIreMatPM/net-2part-displ.py
scripts/test/WIreMatPM/net-2part-strain.py
modified:
pkg/dem/WirePM.cpp
pkg/dem/WirePM.hpp
--
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 'pkg/dem/WirePM.cpp'
--- pkg/dem/WirePM.cpp 2010-12-09 05:46:51 +0000
+++ pkg/dem/WirePM.cpp 2010-12-10 02:42:10 +0000
@@ -12,26 +12,18 @@
CREATE_LOGGER(WireMat);
void WireMat::postLoad(WireMat&){
- if(StrainStressValues.size() < 2)
- throw invalid_argument("WireMat.StrainStressValues: at least two points must be given.");
- if(StrainStressValues[0](0) == 0. && StrainStressValues[0](1) == 0.)
- throw invalid_argument("WireMat.StrainStressValues: Definition must start with values greather then zero (strain>0,stress>0)");
+
+ //BUG: ????? postLoad is called twice,
+ LOG_TRACE( "WireMat::postLoad - update material parameters" );
+
+ if(strainStressValues.empty()) return; // uninitialized object, don't do nothing at all
+ if(strainStressValues.size() < 2)
+ throw invalid_argument("WireMat.strainStressValues: at least two points must be given.");
+ if(strainStressValues[0](0) == 0. && strainStressValues[0](1) == 0.)
+ throw invalid_argument("WireMat.strainStressValues: Definition must start with values greather then zero (strain>0,stress>0)");
// compute cross-sectin area
- As = pow(diameter*0.5,2)*Mathr::PI;
-
- // debug printout
-// Vector2r vbegin = *StrainStressValues.begin();
-// // Vector2r vbegin = StrainStressValues.front();
-// std::cerr << vbegin(0) << endl;
-// std::cerr << vbegin(1) << endl;
-// std::cerr << endl;
-//
-// Real v0 = (*(StrainStressValues.begin()))(0);
-// Real v1 = (*StrainStressValues.begin())(1);
-// std::cerr << v0 << endl;
-// std::cerr << v1 << endl;
-
- //BUG ????? postLoad is called twice,
+ as = pow(diameter*0.5,2)*Mathr::PI;
+
}
@@ -39,6 +31,9 @@
CREATE_LOGGER(Law2_ScGeom_WirePhys_WirePM);
void Law2_ScGeom_WirePhys_WirePM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+
+ LOG_TRACE( "Law2_ScGeom_WirePhys_WirePM::go - contact law" );
+
ScGeom* geom = static_cast<ScGeom*>(ig.get());
WirePhys* phys = static_cast<WirePhys*>(ip.get());
const int &id1 = contact->getId1();
@@ -49,8 +44,8 @@
Real displN = geom->penetrationDepth; // NOTE: ScGeom->penetrationDepth>0 when spheres interpenetrate, and therefore, for wire always negative
/* get reference to values since values are updated for unloading */
- vector<Vector2r> &DFValues = phys->DisplForceValues;
- vector<Real> &kValues = phys->StiffnessValues;
+ vector<Vector2r> &DFValues = phys->displForceValues;
+ vector<Real> &kValues = phys->stiffnessValues;
Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance
@@ -73,7 +68,7 @@
/* compute normal force Fn */
Real Fn = 0.;
if ( D > DFValues[0](0) )
- Fn = kValues[0] * (D-phys->Dplast);
+ Fn = kValues[0] * (D-phys->plastD);
else {
bool isDone = false;
unsigned int i = 0;
@@ -81,7 +76,7 @@
i++;
if ( D > DFValues[i](0) ) {
Fn = DFValues[i-1](1) + (D-DFValues[i-1](0))*kValues[i];
- phys->Dplast = D - Fn/kValues[0];
+ phys->plastD = D - Fn/kValues[0];
isDone = true;
// update values for unloading
DFValues[0](0) = D;
@@ -119,6 +114,8 @@
/* avoid any updates if interactions which already exist */
if(interaction->phys) return;
+ LOG_TRACE( "Ip2_WireMat_WireMat_WirePhys::go - create interaction physics" );
+
ScGeom* geom=dynamic_cast<ScGeom*>(interaction->geom.get());
assert(geom);
@@ -136,8 +133,8 @@
/* ckeck properties of interaction */
if ( mat1->id == mat2->id ) { // interaction of two bodies of the same material
- crossSection = mat1->As;
- SSValues = mat1->StrainStressValues;
+ crossSection = mat1->as;
+ SSValues = mat1->strainStressValues;
if ( (mat1->isDoubleTwist) && (abs(interaction->getId1()-interaction->getId2())==1) ) // bodies which id differs by 1 are double twisted
contactPhysics->isDoubleTwist = true;
else
@@ -146,21 +143,23 @@
else { // interaction of two bodies of two different materials, take weaker material and no double-twist
contactPhysics->isDoubleTwist = false;
if ( mat1->diameter <= mat2->diameter){
- crossSection = mat1->As;
- SSValues = mat1->StrainStressValues;
+ crossSection = mat1->as;
+ SSValues = mat1->strainStressValues;
}
else {
- crossSection = mat2->As;
- SSValues = mat2->StrainStressValues;
+ crossSection = mat2->as;
+ SSValues = mat2->strainStressValues;
}
}
+ if(SSValues.empty()) throw std::invalid_argument("WireMat.strainStressValue is empty!");
+
Real R1 = geom->radius1;
Real R2 = geom->radius2;
Real l0 = R1 + R2 - contactPhysics->initD; // initial lenght of the wire (can be single or double twisted)
- /* compute displscement-force values (tension negative since ScGem is used!) */
+ /* compute thresholddisplscement-force values (tension negative since ScGem is used!) */
vector<Vector2r> DFValues;
for ( vector<Vector2r>::iterator it = SSValues.begin(); it != SSValues.end(); it++ ) {
Vector2r values = Vector2r::Zero();
@@ -188,18 +187,18 @@
}
/* store displscement-force values in physics */
- contactPhysics->DisplForceValues = DFValues;
+ contactPhysics->displForceValues = DFValues;
/* compute stiffness-values of wire */
contactPhysics->kn = k;
kValues.push_back(k);
- for( unsigned int i = 1 ; i < DFValues.size()-1; i++ ) {
+ for( unsigned int i = 1 ; i < DFValues.size(); i++ ) {
Real deltau = -DFValues[i](0) + DFValues[i-1](0);
Real deltaF = -DFValues[i](1) + DFValues[i-1](1);
k = deltaF/deltau;
kValues.push_back(k);
}
- contactPhysics->StiffnessValues = kValues;
+ contactPhysics->stiffnessValues = kValues;
/* set particles as linked */
if ( (scene->iter < linkThresholdIteration))
=== modified file 'pkg/dem/WirePM.hpp'
--- pkg/dem/WirePM.hpp 2010-12-09 05:46:51 +0000
+++ pkg/dem/WirePM.hpp 2010-12-10 02:42:10 +0000
@@ -44,11 +44,11 @@
DECLARE_LOGGER;
YADE_CLASS_BASE_DOC_ATTRS_CTOR(WireMat,ElastMat,"Material for use with the Wire classes",
((Real,diameter,0.0027,," (Diameter of the single wire in [m] (the diameter is used to compute the cross-section area of the wire)."))
- ((vector<Vector2r>,StrainStressValues,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for one single wire. Tension only is considered and the point (0,0) is not needed!"))
+ ((vector<Vector2r>,strainStressValues,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for one single wire. Tension only is considered and the point (0,0) is not needed!"))
((bool,isDoubleTwist,false,,"Type of the mesh. If true two particles of the same material which body ids differ by one will be considered as double-twisted interaction."))
((Real,lambdaEps,0.4,,"Parameter between 0 and 1 to reduce the failure strain of the double-twisted wire (as used by [Bertrand2008]_). [-]"))
- ((Real,lambdak,0.21,,"Parameter between 0 and 1 to compute the elastic stiffness of the double-twisted wire (as used by [Bertrand2008]_): $\\k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S). [-]"))
- ((Real,As,0,Attr::readonly,"Cross-section area of a single wire used for the computation of the limit normal contact forces. [m²]"))
+ ((Real,lambdak,0.21,,"Parameter between 0 and 1 to compute the elastic stiffness of the double-twisted wire (as used by [Bertrand2008]_): $\\k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S)$. [-]"))
+ ((Real,as,0,Attr::readonly,"Cross-section area of a single wire used for the computation of the limit normal contact forces. [m²]"))
,
createIndex();
);
@@ -65,9 +65,9 @@
((Real,initD,0,,"Equilibrium distance for particles. Computed as the initial inter-particular distance when particle are linked."))
((bool,isLinked,false,,"If true particles are linked and will interact. Interactions are linked automatically by the definition of the corresponding interaction radius. The value is false if the wire breaks (no more interaction)."))
((bool,isDoubleTwist,false,,"If true the properties of the interaction will be defined as a double-twisted wire."))
- ((vector<Vector2r>,DisplForceValues,,(Attr::readonly|Attr::noSave),"Defines the values for force-displacement curve."))
- ((vector<Real>,StiffnessValues,,(Attr::readonly|Attr::noSave),"Defines the values for the different stiffness (first value corresponds to elastic stiffness kn)."))
- ((Real,Dplast,0,Attr::readonly,"Plastic part of the inter-particular distance of the previous step. \n\n.. note::\n\t Only elastic displacements are reversible (the elastic stiffness is used for unloading) and compressive forces are inadmissible. The compressive stiffness is assumed to be equal to zero (see [Bertrand2005]_).\n\n.."))
+ ((vector<Vector2r>,displForceValues,,(Attr::readonly|Attr::noSave),"Defines the values for force-displacement curve."))
+ ((vector<Real>,stiffnessValues,,(Attr::readonly|Attr::noSave),"Defines the values for the different stiffness (first value corresponds to elastic stiffness kn)."))
+ ((Real,plastD,0,Attr::readonly,"Plastic part of the inter-particular distance of the previous step. \n\n.. note::\n\t Only elastic displacements are reversible (the elastic stiffness is used for unloading) and compressive forces are inadmissible. The compressive stiffness is assumed to be equal to zero (see [Bertrand2005]_).\n\n.."))
,
createIndex();
,
@@ -85,7 +85,7 @@
FUNCTOR2D(WireMat,WireMat);
DECLARE_LOGGER;
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_WireMat_WireMat_WirePhys,IPhysFunctor,"Converts 2 :yref:`path<WireMat.path>` instances to :yref:`path<WirePhys.path>` with corresponding parameters.",
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_WireMat_WireMat_WirePhys,IPhysFunctor,"Converts 2 :yref:`WireMat` instances to :yref:`WirePhys` with corresponding parameters.",
((int,linkThresholdIteration,1,,"Iteration to create the link."))
);
};
=== added directory 'scripts/test/WIreMatPM'
=== added file 'scripts/test/WIreMatPM/net-2part-displ-unloading.py'
--- scripts/test/WIreMatPM/net-2part-displ-unloading.py 1970-01-01 00:00:00 +0000
+++ scripts/test/WIreMatPM/net-2part-displ-unloading.py 2010-12-10 02:42:10 +0000
@@ -0,0 +1,140 @@
+# -*- coding: utf-8 -*-
+# encoding: utf-8
+from yade import utils, ymport, qt, plot
+
+from yade import log
+log.setLevel('Law2_ScGeom_WirePhys_WirePM',log.TRACE) # must compile with debug option to get logs
+
+## definition of some colors for colored text output in terminal
+BLUE = '\033[94m'
+GREEN = '\033[92m'
+YELLOW = '\033[93m'
+RED = '\033[91m'
+BLACK = '\033[0m'
+
+#### short description of script
+print BLUE+'''
+Simple test for two particles to test the contact law for the WireMat
+by unsing the '''+RED+'''StepDisplacer'''+BLUE+''' with loading and unloading.
+'''+BLACK
+
+#### define parameters for the net
+# mesh opening size
+mos = 80./1000.
+a = mos/sqrt(3)
+# wire diameter
+d = 2.7/1000.
+# particle radius
+radius = d*5.
+# define piecewise lineare stress-strain curve
+strainStressValues=[(0.0019230769,2.5e8),(0.0192,3.2195e8),(0.05,3.8292e8),(0.15,5.1219e8),(0.25,5.5854e8),(0.3,5.6585e8),(0.35,5.6585e8)]
+# elastic material properties
+particleVolume = 4./3.*pow(radius,3)*pi
+particleMass = 3.9/1000.
+density = particleMass/particleVolume
+young = strainStressValues[0][1] / strainStressValues[0][0]
+poisson = 0.3
+
+
+#### material definition
+netMat = O.materials.append(WireMat(young=young,poisson=poisson,density=density,isDoubleTwist=False,diameter=d,strainStressValues=strainStressValues,lambdaEps=0.4,lambdak=0.21))
+
+
+#### create boddies, default: dynamic=True
+O.bodies.append( utils.sphere([0,0,0], radius, wire=False, color=[1,0,0], highlight=False, material=netMat) )
+O.bodies.append( utils.sphere([0,a,0], radius, wire=False, color=[0,1,0], highlight=False, material=netMat) )
+
+FixedSphere=O.bodies[0]
+MovingSphere=O.bodies[1]
+
+FixedSphere.dynamic=False
+MovingSphere.dynamic=False
+
+def addPlotData():
+ if O.iter < 1:
+ plot.addData( Fn=0., un=0. )
+ plot.saveGnuplot('net-2part-displ-unloading')
+ else:
+ try:
+ i=O.interactions[FixedSphere.id,MovingSphere.id]
+ plot.addData( Fn=i.phys.normalForce.norm(), un=(O.bodies[1].state.pos[1]-O.bodies[0].state.pos[1])-a )
+ plot.saveGnuplot('net-2part-displ-unloading')
+ except:
+ print "No interaction!"
+ O.pause()
+
+#### define simulation to create link
+interactionRadius=2.
+O.engines = [
+ ForceResetter(),
+ InsertionSortCollider( [Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='aabb')] ),
+
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='Ig2ssGeom')],
+ [Ip2_WireMat_WireMat_WirePhys(linkThresholdIteration=1,label='interactionPhys')],
+ [Law2_ScGeom_WirePhys_WirePM(linkThresholdIteration=1,label='interactionLaw')]
+ ),
+ PyRunner(initRun=True,iterPeriod=1,command='addPlotData()')
+]
+
+
+#### plot some results
+plot.plots={'un':('Fn',)}
+plot.plot()
+
+
+#### create link (no time step needed since loading is involved in this step)
+O.step() # create cohesive link (cohesiveTresholdIteration=1)
+
+
+#### initializes now the interaction detection factor
+aabb.aabbEnlargeFactor=-1.
+Ig2ssGeom.interactionDetectionFactor=-1.
+
+## time step definition
+## no time step definition is required since setVelocities=False in StepDisplacer
+
+
+#### define simulation loading
+O.engines = [StepDisplacer( ids=[1],mov=Vector3(0,+1e-5,0),rot=Quaternion().Identity,setVelocities=False )] + O.engines
+
+print 'Loading (press enter)'
+raw_input()
+O.run(100,True)
+
+#### define simulation unloading
+O.engines = [StepDisplacer( ids=[1],mov=Vector3(0,-1.3e-5,0),rot=Quaternion().Identity,setVelocities=False )] + O.engines[1:]
+
+print 'Unloading (press enter)'
+raw_input()
+O.run(50,True)
+
+#### define simulation reloading
+O.engines = [StepDisplacer( ids=[1],mov=Vector3(0,+1.6e-5,0),rot=Quaternion().Identity,setVelocities=False )] + O.engines[1:]
+
+print 'Reloading (press enter)'
+raw_input()
+O.run(500,True)
+
+
+#### define simulation unloading
+O.engines = [StepDisplacer( ids=[1],mov=Vector3(0,-1.45e-5,0),rot=Quaternion().Identity,setVelocities=False )] + O.engines[1:]
+
+print 'Reunloading (press enter)'
+raw_input()
+O.run(10,True)
+
+
+#### define simulation reloading
+O.engines = [StepDisplacer( ids=[1],mov=Vector3(0,+1.6e-5,0),rot=Quaternion().Identity,setVelocities=False )] + O.engines[1:]
+
+print 'Reloading (press enter)'
+raw_input()
+O.run(500,True)
+
+
+#### to see it
+v=qt.Controller()
+v=qt.View()
+rr=qt.Renderer()
+rr.intrAllWire=True
=== added file 'scripts/test/WIreMatPM/net-2part-displ.py'
--- scripts/test/WIreMatPM/net-2part-displ.py 1970-01-01 00:00:00 +0000
+++ scripts/test/WIreMatPM/net-2part-displ.py 2010-12-10 02:42:10 +0000
@@ -0,0 +1,104 @@
+# -*- coding: utf-8 -*-
+# encoding: utf-8
+from yade import utils, ymport, qt, plot
+
+## definition of some colors for colored text output in terminal
+BLUE = '\033[94m'
+GREEN = '\033[92m'
+YELLOW = '\033[93m'
+RED = '\033[91m'
+BLACK = '\033[0m'
+
+#### short description of script
+print BLUE+'Simple test for two particles to test contact law with '+RED+'StepDisplacer'+BLUE+'.'+BLACK
+
+#### define parameters for the net
+# mesh opening size
+mos = 80./1000.
+a = mos/sqrt(3)
+# wire diameter
+d = 2.7/1000.
+# particle radius
+radius = d*5.
+# define piecewise lineare stress-strain curve
+strainStressValues=[(0.0019230769,2.5e8),(0.0192,3.2195e8),(0.05,3.8292e8),(0.15,5.1219e8),(0.25,5.5854e8),(0.3,5.6585e8),(0.35,5.6585e8)]
+# elastic material properties
+particleVolume = 4./3.*pow(radius,3)*pi
+particleMass = 3.9/1000.
+density = particleMass/particleVolume
+young = strainStressValues[0][1] / strainStressValues[0][0]
+poisson = 0.3
+
+
+#### material definition
+netMat = O.materials.append(WireMat(young=young,poisson=poisson,density=density,isDoubleTwist=False,diameter=d,strainStressValues=strainStressValues,lambdaEps=0.4,lambdak=0.21))
+
+
+#### create boddies, default: dynamic=True
+O.bodies.append( utils.sphere([0,0,0], radius, wire=False, color=[1,0,0], highlight=False, material=netMat) )
+O.bodies.append( utils.sphere([0,a,0], radius, wire=False, color=[0,1,0], highlight=False, material=netMat) )
+
+FixedSphere=O.bodies[0]
+MovingSphere=O.bodies[1]
+
+FixedSphere.dynamic=False
+MovingSphere.dynamic=False
+
+
+#### define addPlotData
+def addPlotData():
+ if O.iter < 1:
+ plot.addData( Fn=0., un=0. )
+ plot.saveGnuplot('net-2part-displ')
+ else:
+ try:
+ i=O.interactions[FixedSphere.id,MovingSphere.id]
+ plot.addData( Fn=i.phys.normalForce.norm(), un=(O.bodies[1].state.pos[1]-O.bodies[0].state.pos[1])-a )
+ plot.saveGnuplot('net-2part-displ')
+ except:
+ print "No interaction!"
+ O.pause()
+
+
+#### define simulation to create cohesive link
+interactionRadius=2.
+O.engines = [
+ ForceResetter(),
+ InsertionSortCollider( [Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='aabb')] ),
+
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='Ig2ssGeom')],
+ [Ip2_WireMat_WireMat_WirePhys(linkThresholdIteration=1,label='interactionPhys')],
+ [Law2_ScGeom_WirePhys_WirePM(linkThresholdIteration=1,label='interactionLaw')]
+ ),
+ PyRunner(initRun=True,iterPeriod=1,command='addPlotData()')
+]
+
+
+#### plot some results
+plot.plots={'un':('Fn',)}
+plot.plot()
+
+
+#### create link (no time step needed since loading is involved in this step)
+O.step() # create cohesive link (cohesiveTresholdIteration=1)
+
+
+#### initializes now the interaction detection factor
+aabb.aabbEnlargeFactor=-1.
+Ig2ssGeom.interactionDetectionFactor=-1.
+
+
+#### time step definition
+## no time step definition is required since setVelocities=False in StepDisplacer
+
+
+#### define simulation loading
+O.engines = [StepDisplacer( ids=[1],mov=Vector3(0,+1e-5,0),rot=Quaternion().Identity,setVelocities=False )] + O.engines
+
+
+#### to see it
+v=qt.Controller()
+v=qt.View()
+rr=qt.Renderer()
+rr.intrAllWire=True
=== added file 'scripts/test/WIreMatPM/net-2part-strain.py'
--- scripts/test/WIreMatPM/net-2part-strain.py 1970-01-01 00:00:00 +0000
+++ scripts/test/WIreMatPM/net-2part-strain.py 2010-12-10 02:42:10 +0000
@@ -0,0 +1,113 @@
+# -*- coding: utf-8 -*-
+# encoding: utf-8
+from yade import utils, ymport, qt
+
+#### logging
+from yade import log
+log.setLevel('Law2_ScGeom_WirePhys_WirePM',log.TRACE) # must compile with debug option to get logs
+#log.setLevel('Law2_ScGeom_WirePhys_WirePM',log.DEBUG)
+#log.setLevel('',log.WARN)
+
+## definition of some colors for colored text output in terminal
+BLUE = '\033[94m'
+GREEN = '\033[92m'
+YELLOW = '\033[93m'
+RED = '\033[91m'
+BLACK = '\033[0m'
+
+#### short description of script
+print BLUE+'Simple test for two particles to test contact law with '+RED+'UniaxialStrainer'+BLUE+'.'+BLACK
+
+#### define parameters for the net
+# mesh opening size
+mos = 80./1000.
+a = mos/sqrt(3)
+# wire diameter
+d = 2.7/1000.
+# particle radius
+radius = d*5.
+# define piecewise lineare stress-strain curve
+strainStressValues=[(0.0019230769,2.5e8),(0.0192,3.2195e8),(0.05,3.8292e8),(0.15,5.1219e8),(0.25,5.5854e8),(0.3,5.6585e8),(0.35,5.6585e8)]
+# elastic material properties
+particleVolume = 4./3.*pow(radius,3)*pi
+particleMass = 3.9/1000.
+density = particleMass/particleVolume
+young = strainStressValues[0][1] / strainStressValues[0][0]
+poisson = 0.3
+
+
+#### material definition
+netMat = O.materials.append(WireMat(young=young,poisson=poisson,density=density,isDoubleTwist=False,diameter=d,strainStressValues=strainStressValues,lambdaEps=0.4,lambdak=0.21))
+
+
+#### create boddies, default: dynamic=True
+O.bodies.append( utils.sphere([0,0,0], radius, wire=False, color=[1,0,0], highlight=False, material=netMat) )
+O.bodies.append( utils.sphere([0,a,0], radius, wire=False, color=[0,1,0], highlight=False, material=netMat) )
+
+FixedSphere=O.bodies[0]
+MovingSphere=O.bodies[1]
+
+FixedSphere.dynamic=True
+MovingSphere.dynamic=True
+
+
+#### initialize values for UniaxialStrainer
+bb = utils.uniaxialTestFeatures(axis=1)
+negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
+strainRateTension = 1./a
+setSpeeds = True
+
+
+#### define simulation to create cohesive link
+interactionRadius=2.
+O.engines = [
+ ForceResetter(),
+ InsertionSortCollider( [Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='aabb')] ),
+
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='Ig2ssGeom')],
+ [Ip2_WireMat_WireMat_WirePhys(linkThresholdIteration=1,label='interactionPhys')],
+ [Law2_ScGeom_WirePhys_WirePM(linkThresholdIteration=1,label='interactionLaw')]
+ ),
+]
+
+
+#### create link (no time step needed since loading is involved in this step)
+O.step() # create cohesive link (cohesiveTresholdIteration=1)
+
+
+#### initializes now the interaction detection factor
+aabb.aabbEnlargeFactor=-1.
+Ig2ssGeom.interactionDetectionFactor=-1.
+
+## time step definition
+O.dt = 1e-5
+## critical time step proposed by Bertrand
+#O.dt = 0.2*sqrt(particleMass/(2.*O.interactions[0,1].phys.kn))
+
+#### plot some results
+from math import *
+from yade import plot
+
+plot.plots={'un':('Fn',)}
+plot.plot()
+
+def addPlotData():
+ try:
+ i=O.interactions[FixedSphere.id,MovingSphere.id]
+ plot.addData( Fn=i.phys.normalForce.norm(), un=(O.bodies[1].state.pos[1]-O.bodies[0].state.pos[1])-a )
+ plot.saveGnuplot('net-2part-strain')
+ except:
+ print "No interaction!"
+ O.pause()
+
+
+#### define simulation
+O.engines += [UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=1,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=True,blockRotations=True,setSpeeds=setSpeeds,label='strainer')] + [NewtonIntegrator(damping=0.0)] + [PyRunner(initRun=True,iterPeriod=1,command='addPlotData()')]
+
+
+#### to see it
+v=qt.Controller()
+v=qt.View()
+rr=qt.Renderer()
+rr.intrAllWire=True
