yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2828
committer: Janek Kozicki <cosurgi@xxxxxxxxxx>
branch nick: yade
timestamp: Thu 2011-04-21 18:03:07 +0200
message:
  fix warning,
  add plastic dissipation tracking
modified:
  examples/simple-scene/simple-scene-energy-tracking.py
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/CohesiveFrictionalContactLaw.hpp
  pkg/dem/ElasticContactLaw.cpp


--
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-16 11:35:08 +0000
+++ examples/simple-scene/simple-scene-energy-tracking.py	2011-04-21 16:03:07 +0000
@@ -10,9 +10,9 @@
 # initial angular velocity
 angVel = 3.0
 # use two spheres?
-two_spheres = False
+two_spheres =True
 # sphere rotating more?
-rotate_in_two_directions = False
+rotate_in_two_directions = True
 
 ############################################
 ##### material                         #####
@@ -79,16 +79,20 @@
 	'E_kin_translation':'Translation energy: E_kin=m*V^2/2',
 	'E_kin_rotation':'Rotation energy: E_kin=I*$\omega$^2/2',
 	'E_pot':'Gravitational potential: E_pot=m*g*h',
+	'E_plastic':'Plastic dissipation on shearing: E_pl=F*$\Delta$F/k',
 	'total':'total',
 	'total_plus_damp':'total + damping'}
 
-plot.plots={'t':('normal_Work',
-		'shear_Work',
-		'E_kin_translation',
-		'E_kin_rotation',
-		'E_pot',
-		'total',
-		'total_plus_damp')}
+plot.plots={'t':(
+		('normal_Work','b-'),
+		('shear_Work','r-'),
+		('E_kin_translation','b-.'),
+		('E_kin_rotation','r-.'),
+		('E_plastic','c-'),
+		('E_pot','y-'),
+		('total','k:'),
+		('total_plus_damp','k-')
+		)}
 
 ## this function is called by plotDataCollector
 ## it should add data with the labels that we will plot
@@ -99,11 +103,13 @@
 	E_kin_translation = 0
 	E_kin_rotation    = 0
 	E_pot		  = 0
+	E_plastic	  = 0
+	E_tracker	  = dict(O.energy.items())
 
 	if(two_spheres):## for more bodies we better use the energy tracker, because it's tracking all bodies
-		E_kin_translation = dict(O.energy.items())['kinTrans']
-		E_kin_rotation    = dict(O.energy.items())['kinRot']
-		E_pot		  = dict(O.energy.items())['gravWork'] 
+		E_kin_translation = E_tracker['kinTrans']
+		E_kin_rotation    = E_tracker['kinRot']
+		E_pot		  = E_tracker['gravWork'] 
 
 	else: ## for one sphere we can just calculate, and it will be correct
 		sph=O.bodies[1]
@@ -116,10 +122,13 @@
 		E_kin_rotation    = I*w**2.0/2.0
 		E_pot		  = m*g*h
 
-	total = normal_Work + shear_Work + E_kin_translation + E_kin_rotation + E_pot
+	if('plastDissip' in E_tracker):
+		E_plastic	  = E_tracker['plastDissip']
+
+	total = normal_Work + shear_Work + E_plastic + E_kin_translation + E_kin_rotation + E_pot
 	total_plus_damp	  = 0
 	if(damping!=0):
-		total_plus_damp	  = total + dict(O.energy.items())['nonviscDamp']
+		total_plus_damp	  = total + E_tracker['nonviscDamp']
 	else:	
 		total_plus_damp	  = total
 	plot.addData(
@@ -129,6 +138,7 @@
 		E_kin_translation = E_kin_translation,
 		E_kin_rotation    = E_kin_rotation   ,
 		E_pot		  = E_pot		 ,
+		E_plastic	  = E_plastic ,
 		total		  = total		 ,
 		total_plus_damp	  = total_plus_damp	 ,
 	)

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp	2011-04-21 13:59:13 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2011-04-21 16:03:07 +0000
@@ -16,14 +16,11 @@
 YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom6D_CohFrictPhys_CohesionMoment));
 CREATE_LOGGER(Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
 
-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());  //Commented due to warning
 		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
 		if (phys) {
 			normEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
@@ -36,7 +33,6 @@
 	Real shearEnergy=0;
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
 		if(!I->isReal()) continue;
-		//ScGeom6D* scg = YADE_CAST<ScGeom6D*>(I->geom.get());  //Commented due to warning
 		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
 		if (phys) {
 			shearEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
@@ -105,7 +101,10 @@
 				maxFs = max((Real) 0, Fn*currentContactPhysics->tangensOfFrictionAngle);
 			}
 			maxFs = maxFs / Fs;
+			Vector3r trialForce=shearForce;
 			shearForce *= maxFs;
+			Real dissip=((1/currentContactPhysics->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
+			if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);
 			if (Fn<0)  currentContactPhysics->normalForce = Vector3r::Zero();//Vector3r::Zero()
 		}
 		applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position);

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp	2011-04-16 11:35:08 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp	2011-04-21 16:03:07 +0000
@@ -26,6 +26,7 @@
 		((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."))
+		((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
 		((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.")

=== modified file 'pkg/dem/ElasticContactLaw.cpp'
--- pkg/dem/ElasticContactLaw.cpp	2010-12-31 14:35:21 +0000
+++ pkg/dem/ElasticContactLaw.cpp	2011-04-21 16:03:07 +0000
@@ -71,7 +71,8 @@
 			Real ratio = sqrt(maxFs) / shearForce.norm();
 			shearForce *= ratio;}
 	} else {
-		//almost the same with additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
+		//almost the same with additional Vector3r instatinated for energy tracing, 
+		//duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
 		if(shearForce.squaredNorm() > maxFs){
 			Real ratio = sqrt(maxFs) / shearForce.norm();
 			Vector3r trialForce=shearForce;//store prev force for definition of plastic slip