yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2818
committer: Janek Kozicki <cosurgi@xxxxxxxxxx>
branch nick: yade
timestamp: Sat 2011-04-16 13:35:08 +0200
message:
  last correction for simple scene energy tracking example - I forgot to commit changes in CohesiveFrictionalContactLaw
  but in fact maybe better if those methods for calculating elastic normal & shear energy would be moved somewhere else, maybe to shop? I'm not sure.
modified:
  examples/simple-scene/simple-scene-energy-tracking.py
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/CohesiveFrictionalContactLaw.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 'examples/simple-scene/simple-scene-energy-tracking.py'
--- examples/simple-scene/simple-scene-energy-tracking.py	2011-04-15 19:57:48 +0000
+++ examples/simple-scene/simple-scene-energy-tracking.py	2011-04-16 11:35:08 +0000
@@ -80,13 +80,12 @@
 	'E_kin_rotation':'Rotation energy: E_kin=I*$\omega$^2/2',
 	'E_pot':'Gravitational potential: E_pot=m*g*h',
 	'total':'total',
-	'total_plus_damp':'total + daping'}
+	'total_plus_damp':'total + damping'}
 
 plot.plots={'t':('normal_Work',
 		'shear_Work',
 		'E_kin_translation',
 		'E_kin_rotation',
-		#'E_kin_r','E_kin_tr','E_pot_',  ## those are from energy tracker
 		'E_pot',
 		'total',
 		'total_plus_damp')}
@@ -95,9 +94,8 @@
 ## it should add data with the labels that we will plot
 ## if a datum is not specified (but exists), it will be NaN and will not be plotted
 def myAddPlotData():
-	sph=O.bodies[1]
-	normal_Work	  = law.normElastEnergyForce()
-	shear_Work	  = law.shearElastEnergyForce()
+	normal_Work	  = law.normElastEnergy()
+	shear_Work	  = law.shearElastEnergy()
 	E_kin_translation = 0
 	E_kin_rotation    = 0
 	E_pot		  = 0
@@ -108,6 +106,7 @@
 		E_pot		  = dict(O.energy.items())['gravWork'] 
 
 	else: ## for one sphere we can just calculate, and it will be correct
+		sph=O.bodies[1]
 		h=sph.state.pos[2]
 		V=sph.state.vel.norm()
 		w=sph.state.angVel.norm()

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp	2011-01-26 18:19:52 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2011-04-16 11:35:08 +0000
@@ -18,6 +18,33 @@
 
 Vector3r translation_vect_ ( 0.10,0,0 );
 
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy()
+{
+	Real normEnergy=0;
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		ScGeom6D* scg = YADE_CAST<ScGeom6D*>(I->geom.get());
+		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+		if (phys) {
+			normEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
+		}
+	}
+	return normEnergy;
+}
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::shearElastEnergy()
+{
+	Real shearEnergy=0;
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		ScGeom6D* scg = YADE_CAST<ScGeom6D*>(I->geom.get());
+		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+		if (phys) {
+			shearEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
+		}
+	}
+	return shearEnergy;
+}
+
 void CohesiveFrictionalContactLaw::action()
 {
 	if(!functor) functor=shared_ptr<Law2_ScGeom6D_CohFrictPhys_CohesionMoment>(new Law2_ScGeom6D_CohFrictPhys_CohesionMoment);

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp	2011-01-26 18:19:52 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp	2011-04-16 11:35:08 +0000
@@ -17,14 +17,19 @@
 
 class Law2_ScGeom6D_CohFrictPhys_CohesionMoment: public LawFunctor{
 	public:
+		Real normElastEnergy();
+		Real shearElastEnergy();
 	virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
-	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. There is no maximum value of moments in the current implementation, though they could be added in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.\n\n.. note::\n  Periodicity is not handled yet in this law.",
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. There is no maximum value of moments in the current implementation, though they could be added in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.\n\n.. note::\n  Periodicity is not handled yet in this law.",
 		((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
 		((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
 		((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
 		((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
 		((bool,useIncrementalForm,false,,"use the incremental formulation to compute bending and twisting moments. Creep on the twisting moment is not included in such a case."))
 		((Real,creep_viscosity,1,,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_CohFrictMat_CohFrictMat_CohFrictPhys..."))
+		,,
+		.def("normElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
+		.def("shearElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
 	);
 	FUNCTOR2D(ScGeom6D,CohFrictPhys);
 	DECLARE_LOGGER;