yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2264
committer: Chiara Modenese <chia@engs-018373>
branch nick: trunk
timestamp: Tue 2010-06-01 13:04:04 +0100
message:
  - Small changes in HM to handle PBC if scene is periodic (same logic as in ElasticContactLaw.cpp)
modified:
  pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
  pkg/dem/Engine/GlobalEngine/HertzMindlin.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/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-05-26 08:18:36 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-06-01 12:04:04 +0000
@@ -93,7 +93,9 @@
 	/*** SHEAR FORCE ***/
 	phys->ks = phys->kso*std::pow(uN,0.5); // get tangential stiffness (this is a tangent value, so we can pass it to the GSTS)
 	Vector3r& trialFs = phys->shearForce;
-  	Vector3r dus =scg->rotateAndGetShear(trialFs, phys->prevNormal, de1, de2, dt, /*define this for periodicity, see :ElasticContactLaw.cpp:71*/ Vector3r::Zero(), preventGranularRatcheting);
+	// define shift to handle periodicity
+	Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
+  	Vector3r dus =scg->rotateAndGetShear(trialFs, phys->prevNormal, de1, de2, dt, shiftVel, preventGranularRatcheting);
 	//Linear elasticity giving "trial" shear force
 	trialFs -= phys->ks*dus;
 
@@ -105,7 +107,15 @@
 
 
 	/*** APPLY FORCES ***/
+	if (!scene->isPeriodic)
 	applyForceAtContactPoint(-phys->normalForce-trialFs , scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position, ncb);
+	else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used
+		Vector3r force = -phys->normalForce-trialFs;
+		ncb->forces.addForce(id1,force);
+		ncb->forces.addForce(id2,-force);
+		ncb->forces.addTorque(id1,(scg->radius1-0.5*scg->penetrationDepth)* scg->normal.cross(force));
+		ncb->forces.addTorque(id2,(scg->radius2-0.5*scg->penetrationDepth)* scg->normal.cross(force));
+	}
 	phys->prevNormal = scg->normal;
 }
 

=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp	2010-02-23 14:22:35 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp	2010-06-01 12:04:04 +0000
@@ -49,7 +49,7 @@
 	virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
 	FUNCTOR2D(ScGeom,MindlinPhys);
 	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Mindlin's formulation.",
-			((bool,preventGranularRatcheting,false,"bool to avoid granular ratcheting"))
+			((bool,preventGranularRatcheting,true,"bool to avoid granular ratcheting"))
 
 	);
 	DECLARE_LOGGER;