yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2170
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Wed 2010-04-21 19:17:54 +0200
message:
  - Make the radius of fictious sphere equal to the one of real sphere in box-sphere geometry.
  - Implement energy tracing in ElasticContactLaw
  - Remove some "#ifdef SCG_SHEAR"  
modified:
  pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
  pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp
  pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.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/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2009-12-13 20:11:31 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-04-21 17:17:54 +0000
@@ -65,7 +65,7 @@
 }
 #endif
 
-void ScGeom::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
+Vector3r ScGeom::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
 
 	Vector3r axis;
 	Real angle;
@@ -113,6 +113,7 @@
 	Vector3r shearVelocity = relativeVelocity-normal.Dot(relativeVelocity)*normal;
 	Vector3r shearDisplacement = shearVelocity*dt;
 	shearForce -= ks*shearDisplacement;
+	return shearDisplacement;
 }
 
 

=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-04-20 14:57:37 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-04-21 17:17:54 +0000
@@ -27,7 +27,7 @@
 		
 		virtual ~ScGeom();
 
-		void updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
+		Vector3r updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
 
 	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`spheres<Sphere>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
 		((Vector3r,contactPoint,Vector3r::ZERO,"Reference point of the contact. |ycomp|"))

=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp	2010-02-09 20:22:04 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp	2010-04-21 17:17:54 +0000
@@ -93,17 +93,15 @@
 		if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
 		else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
 
-		#ifdef SCG_SHEAR
-			if(isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); }
-			else {scm->prevNormal=scm->normal;}
-		#endif
+		if(isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); }
+		else {scm->prevNormal=scm->normal;}
 			
 		// contact point is in the middle of overlapping volumes
 		//(in the direction of penetration, which is normal to the box surface closest to sphere center) of overlapping volumes
 		scm->contactPoint = 0.5*(pt1+pt2);
 		scm->normal = pt1-pt2; scm->normal.Normalize();
 		scm->penetrationDepth = (pt1-pt2).Length();
-		scm->radius1 = s->radius*2;
+		scm->radius1 = s->radius;
 		scm->radius2 = s->radius;
 		c->interactionGeometry = scm;
 	} else { // outside
@@ -140,10 +138,9 @@
 		bool isNew=!c->interactionGeometry;
 		if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
 		else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);	
-		#ifdef SCG_SHEAR
-			if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
-			else {scm->prevNormal=scm->normal;}
-		#endif
+		if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
+		else {scm->prevNormal=scm->normal;}
+		
 		scm->contactPoint = 0.5*(pt1+pt2);
 		//scm->normal = pt1-pt2; scm->normal.Normalize();
 		//scm->penetrationDepth = (pt1-pt2).Length();
@@ -151,7 +148,7 @@
 		scm->penetrationDepth = depth;
 		
 		
-		scm->radius1 = s->radius*2;
+		scm->radius1 = s->radius;
 		scm->radius2 = s->radius;
 		c->interactionGeometry = scm;
 	}

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-04-19 13:33:45 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-04-21 17:17:54 +0000
@@ -33,13 +33,23 @@
 	}
 }
 
+Real Law2_ScGeom_FrictPhys_Basic::elasticEnergy()
+{
+	Real energy=0;
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		FrictPhys* phys = dynamic_cast<FrictPhys*>(I->interactionPhysics.get());
+		if(phys) {
+			energy += 0.5*(phys->normalForce.SquaredLength()/phys->kn + phys->shearForce.SquaredLength()/phys->ks);}
+	}
+	return energy;
+}
 
 CREATE_LOGGER(Law2_ScGeom_FrictPhys_Basic);
 void Law2_ScGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb){
 	Real dt = Omega::instance().getTimeStep();
 	int id1 = contact->getId1(), id2 = contact->getId2();
-// 			// FIXME: mask handling should move to LawFunctor itself, outside the functors
-// 			if( !(Body::byId(id1,ncb)->getGroupMask() & Body::byId(id2,ncb)->getGroupMask() & sdecGroupMask) ) return;
+
 	ScGeom*    currentContactGeometry= static_cast<ScGeom*>(ig.get());
 	FrictPhys* currentContactPhysics = static_cast<FrictPhys*>(ip.get());
 	if(currentContactGeometry->penetrationDepth <0){
@@ -54,18 +64,42 @@
 	Real un=currentContactGeometry->penetrationDepth;
 	TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
 	currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
-	if(useShear){
-		currentContactGeometry->updateShear(de1,de2,dt);
-		shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
-	} else {
-		currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);}
-	// PFC3d SlipModel, is using friction angle. CoulombCriterion
-	Real maxFs = currentContactPhysics->normalForce.SquaredLength()*
-			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
-	if( shearForce.SquaredLength() > maxFs ){
-		Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length();
-		shearForce *= ratio;
-		if(useShear) currentContactGeometry->shear*=ratio;}
+	
+	if (!traceEnergy){//Update force but don't compute energy terms
+		if(useShear){
+			currentContactGeometry->updateShear(de1,de2,dt);
+			shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
+		} else {
+			currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);}
+		// PFC3d SlipModel, is using friction angle. CoulombCriterion
+		Real maxFs = currentContactPhysics->normalForce.SquaredLength()*
+			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
+		if( shearForce.SquaredLength() > maxFs ){
+			Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length();
+			shearForce *= ratio;
+			if(useShear) currentContactGeometry->shear*=ratio;}
+	} else {//almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
+		Vector3r prevForce=shearForce;//store prev force for definition of plastic slip
+		if(useShear) throw ("energy tracing not defined with useShear==true");
+		/*{
+			currentContactGeometry->updateShear(de1,de2,dt);
+			shearForce=currentContactPhysics->ks*currentContactGeometry->shear;//FIXME : energy terms if useShear?
+		} else {*/
+		Vector3r shearDisp = currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+ 		//}
+		// PFC3d SlipModel, is using friction angle. CoulombCriterion
+		Real maxFs = currentContactPhysics->normalForce.SquaredLength()*
+			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
+		if( shearForce.SquaredLength() > maxFs ){
+			Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length();
+			//define the plastic work input and increment the total plastic energy dissipated
+			plasticDissipation +=
+			(shearDisp+(1/currentContactPhysics->ks)*(shearForce-prevForce))//plastic disp.
+			.Dot(shearForce);//active force
+			shearForce *= ratio;
+			//if(useShear) currentContactGeometry->shear*=ratio;
+		}
+	}	
 	applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
 	currentContactPhysics->prevNormal = currentContactGeometry->normal;
 }

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp	2010-04-19 13:33:45 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp	2010-04-21 17:17:54 +0000
@@ -21,10 +21,15 @@
 class Law2_ScGeom_FrictPhys_Basic: public LawFunctor{
 	public:
 		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
-	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom`.",
+		Real elasticEnergy ();
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom`.",
 		((int,sdecGroupMask,1,"Bitmask for allowing collision between particles :yref:`Body::groupMask`"))
 		((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,useShear,false,"Use ScGeom::updateShear rather than ScGeom::updateShearForce for shear force computation."))
+		((bool,traceEnergy,false,"Define the total energy dissipated in plastic slips at all contacts."))
+		((Real,plasticDissipation,0,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_Basic::traceEnergy` is true. |yupdate|"))
+		,,
+		.def("elasticEnergy",&Law2_ScGeom_FrictPhys_Basic::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
 	);
 	FUNCTOR2D(ScGeom,FrictPhys);
 	DECLARE_LOGGER;