yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10646
[Branch ~yade-pkg/yade/git-trunk] Rev 3875: Use MatchMaker to set tc, et and en in ViscEl model
------------------------------------------------------------
revno: 3875
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-02 17:33:41 +0200
message:
Use MatchMaker to set tc, et and en in ViscEl model
added:
examples/ViscElMatchMaker.py
modified:
pkg/dem/ViscoelasticCapillarPM.cpp
pkg/dem/ViscoelasticCapillarPM.hpp
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== added file 'examples/ViscElMatchMaker.py'
--- examples/ViscElMatchMaker.py 1970-01-01 00:00:00 +0000
+++ examples/ViscElMatchMaker.py 2014-04-02 15:33:41 +0000
@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+# This example shows, how matchmaker can be used to
+# set the parameters of ViscoElastic model.
+
+from yade import utils, plot
+o = Omega()
+fr = 0.5;rho=2000
+tc = 0.001; en = 0.5; et = 0.5; o.dt = 0.002*tc
+
+
+r1 = 0.002381
+r2 = 0.002381
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,tc=tc,en=en,et=et,density=rho))
+mat2 = O.materials.append(ViscElMat(frictionAngle=fr,tc=tc,en=en,et=et,density=rho))
+mat3 = O.materials.append(ViscElMat(frictionAngle=fr,tc=tc,en=en,et=et,density=rho))
+
+
+id11 = O.bodies.append(sphere(center=[0,0,0],radius=r1,material=mat1,fixed=True,color=[0,0,1]))
+id12 = O.bodies.append(sphere(center=[0,0,(r1+r2+0.005*r2)],radius=r2,material=mat2,fixed=False,color=[0,0,1]))
+
+id21 = O.bodies.append(sphere(center=[3*r1,0,0],radius=r1,material=mat1,fixed=True,color=[0,1,0]))
+id22 = O.bodies.append(sphere(center=[3*r1,0,(r1+r2+0.005*r2)],radius=r2,material=mat3,fixed=False,color=[0,1,0]))
+
+id31 = O.bodies.append(sphere(center=[6*r1,0,0],radius=r1,material=mat2,fixed=True,color=[1,0,0]))
+id32 = O.bodies.append(sphere(center=[6*r1,0,(r1+r2+0.005*r2)],radius=r2,material=mat3,fixed=False,color=[1,0,0]))
+
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=(r1+r2)*5.0),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom()],
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys(
+ en=MatchMaker(matches=((mat1,mat2,.9),(mat1,mat3,.5),(mat2,mat3,.1))), # Set parameters
+ et=MatchMaker(matches=((mat1,mat2,.9),(mat1,mat3,.5),(mat2,mat3,.1)))
+ )],
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
+ PyRunner(command='addPlotData()',iterPeriod=100),
+]
+
+vel=-0.50
+O.bodies[id12].state.vel=[0,0,vel]
+O.bodies[id22].state.vel=[0,0,vel]
+O.bodies[id32].state.vel=[0,0,vel]
+
+def addPlotData():
+ s1 = (O.bodies[id12].state.pos[2]-O.bodies[id11].state.pos[2])-(r1+r2)
+ s2 = (O.bodies[id22].state.pos[2]-O.bodies[id11].state.pos[2])-(r1+r2)
+ s3 = (O.bodies[id32].state.pos[2]-O.bodies[id11].state.pos[2])-(r1+r2)
+ plot.addData(mat1mat2=s1,mat1mat3=s2,mat2mat3=s3,it=O.iter)
+
+
+
+plot.plots={'it':('mat1mat2','mat1mat3','mat2mat3')}; plot.plot()
+
+O.step()
+from yade import qt
+qt.View()
+
+#O.run(100000, True)
+#plot.saveGnuplot('sim-data_')
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-04-02 15:33:41 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-04-02 15:33:41 +0000
@@ -18,16 +18,8 @@
// no updates of an existing contact
if(interaction->phys) return;
- ViscElPhys* physT = Calculate_ViscElMat_ViscElMat_ViscElPhys(b1, b2, interaction);
-
shared_ptr<ViscElCapPhys> phys (new ViscElCapPhys());
- phys->kn = physT->kn;
- phys->ks = physT->ks;
- phys->cn = physT->cn;
- phys->cs = physT->cs;
- phys->tangensOfFrictionAngle = physT->tangensOfFrictionAngle;
- phys->shearForce = physT->shearForce;
- phys->mRtype = physT->mRtype;
+ Calculate_ViscElMat_ViscElMat_ViscElPhys(b1, b2, interaction, phys);
ViscElCapMat* mat1 = static_cast<ViscElCapMat*>(b1.get());
ViscElCapMat* mat2 = static_cast<ViscElCapMat*>(b2.get());
=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp 2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp 2014-04-02 15:33:41 +0000
@@ -36,12 +36,12 @@
REGISTER_SERIALIZABLE(ViscElCapPhys);
/// Convert material to interaction physics.
-class Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys: public IPhysFunctor {
+class Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys: public Ip2_ViscElMat_ViscElMat_ViscElPhys {
public :
virtual void go(const shared_ptr<Material>& b1,
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
- YADE_CLASS_BASE_DOC(Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElCapMat` to :yref:`ViscElCapPhys` using the rule of consecutive connection.");
+ YADE_CLASS_BASE_DOC(Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys,Ip2_ViscElMat_ViscElMat_ViscElPhys,"Convert 2 instances of :yref:`ViscElCapMat` to :yref:`ViscElCapPhys` using the rule of consecutive connection.");
FUNCTOR2D(ViscElCapMat,ViscElCapMat);
};
REGISTER_SERIALIZABLE(Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys);
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2014-04-02 15:33:41 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2014-04-02 15:33:41 +0000
@@ -18,8 +18,9 @@
void Ip2_ViscElMat_ViscElMat_ViscElPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
// no updates of an existing contact
if(interaction->phys) return;
- ViscElPhys * phys = Calculate_ViscElMat_ViscElMat_ViscElPhys(b1, b2, interaction);
- interaction->phys = shared_ptr<ViscElPhys>(phys);
+ shared_ptr<ViscElPhys> phys (new ViscElPhys());
+ Calculate_ViscElMat_ViscElMat_ViscElPhys(b1, b2, interaction, phys);
+ interaction->phys = phys;
}
/* Law2_ScGeom_ViscElPhys_Basic */
@@ -116,9 +117,7 @@
torque2 = c2x.cross(force)-momentResistance;
}
-ViscElPhys* Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
- ViscElPhys* phys = new ViscElPhys();
-
+void Ip2_ViscElMat_ViscElMat_ViscElPhys::Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction, shared_ptr<ViscElPhys> phys) {
ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
ViscElMat* mat2 = static_cast<ViscElMat*>(b2.get());
Real mass1 = 1.0;
@@ -160,17 +159,17 @@
Real ks1 = 0.0; Real ks2 = 0.0;
Real cs1 = 0.0; Real cs2 = 0.0;
- if ((isnormal(mat1->tc)) and (isnormal(mat1->en)) and (isnormal(mat1->et))) {
+ if (((isnormal(mat1->tc)) and (isnormal(mat1->en)) and (isnormal(mat1->et))) or ((tc) and (en) and (et))) {
//Set parameters according to [Pournin2001]
- const Real tc = (mat1->tc+mat2->tc)/2.0;
- const Real en = (mat1->en+mat2->en)/2.0;
- const Real et = (mat1->et+mat2->et)/2.0;
+ const Real Tc = (tc) ? (*tc)(mat1->id,mat2->id) : (mat1->tc+mat2->tc)/2.0;
+ const Real En = (en) ? (*en)(mat1->id,mat2->id) : (mat1->en+mat2->en)/2.0;
+ const Real Et = (et) ? (*et)(mat1->id,mat2->id) : (mat1->et+mat2->et)/2.0;
- kn1 = kn2 = 1/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(en),2) )*massR;
- cn1 = cn2 = -2.0 /tc * log(en)*massR;
- ks1 = ks2 = 2.0/7.0 /tc/tc * ( Mathr::PI*Mathr::PI + pow(log(et),2) )*massR;
- cs1 = cs2 = -2.0/7.0 /tc * log(et)*massR;
+ kn1 = kn2 = 1/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR;
+ cn1 = cn2 = -2.0 /Tc * log(En)*massR;
+ ks1 = ks2 = 2.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR;
+ cs1 = cs2 = -2.0/7.0 /Tc * log(Et)*massR;
if (abs(cn1) <= Mathr::ZERO_TOLERANCE ) cn1=0;
if (abs(cn2) <= Mathr::ZERO_TOLERANCE ) cn2=0;
@@ -223,8 +222,6 @@
} else {
phys->mRtype = mRtype1;
}
-
- return phys;
}
/* Contact parameter calculation function */
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-04-02 15:33:41 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-04-02 15:33:41 +0000
@@ -9,6 +9,7 @@
#include<yade/pkg/common/Dispatching.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/pkg/dem/DemXDofGeom.hpp>
+#include<yade/pkg/common/MatchMaker.hpp>
/* Simple viscoelastic model */
@@ -61,9 +62,12 @@
virtual void go(const shared_ptr<Material>& b1,
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
- YADE_CLASS_BASE_DOC(Ip2_ViscElMat_ViscElMat_ViscElPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElMat` to :yref:`ViscElPhys` using the rule of consecutive connection.");
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_ViscElMat_ViscElMat_ViscElPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElMat` to :yref:`ViscElPhys` using the rule of consecutive connection.",
+ ((shared_ptr<MatchMaker>,tc,,,"Contact time"))
+ ((shared_ptr<MatchMaker>,en,,,"Restitution coefficient in normal direction"))
+ ((shared_ptr<MatchMaker>,et,,,"Restitution coefficient in tangential direction")));
+ virtual void Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction, shared_ptr<ViscElPhys> phys);
FUNCTOR2D(ViscElMat,ViscElMat);
-
};
REGISTER_SERIALIZABLE(Ip2_ViscElMat_ViscElMat_ViscElPhys);
@@ -80,5 +84,4 @@
REGISTER_SERIALIZABLE(Law2_ScGeom_ViscElPhys_Basic);
Real contactParameterCalculation(const Real& l1,const Real& l2, const bool& massMultiply);
-ViscElPhys* Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction);
void computeForceTorqueViscEl(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2);
