← Back to team overview

yade-dev team mailing list archive

  1. yade-dev team
  2. Mailing list archive
  3. Message #04587

[Branch ~yade-dev/yade/trunk] Rev 2255: This commit complete the implementation of periodicity for incremental formulation, with updated ...

  Thread Previous • Date Previous • • Thread Next 
------------------------------------------------------------
revno: 2255
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Wed 2010-05-26 10:18:36 +0200
message:
  This commit complete the implementation of periodicity for incremental formulation, with updated relative velocity across periods.
  It has been tested carefully in periodic triaxial tests.
  1. updateShearForce is renamed rotateAndGetShear and DOES NOT update shear force any more : this is Law2's job. Two versions are available, for resp. periodic and
  non-periodic case so that people don't have to worry about the additional "shift" in function parameters.
  2. put the line fs+=ks*dus in laws using the former ScGeom::updateShearForce.
  3. Put back Hsize in Cell (I need that each time I write a line in periodicity, really, please don't remove it).
  4. The rest is details (fix wm3 incompatibilies, formating, documentation, etc.).
modified:
  core/Cell.cpp
  core/Cell.hpp
  pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.cpp
  pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
  pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_WallSphere.cpp
  pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
  pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp
  pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp
  pkg/dem/Engine/GlobalEngine/HertzMindlin.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 'core/Cell.cpp'
--- core/Cell.cpp	2010-05-03 12:17:44 +0000
+++ core/Cell.cpp	2010-05-26 08:18:36 +0000
@@ -6,8 +6,8 @@
 	_trsfInc=dt*velGrad;
 	// total transformation; M = (Id+G).M = F.M
 	trsf+=_trsfInc*trsf;
-	// Hsize will contain colums with transformed base vectors
-	Matrix3r Hsize(Matrix3r::Zero()); Hsize(0,0)=refSize[0]; Hsize(1,1)=refSize[1]; Hsize(2,2)=refSize[2]; // later with eigen: Hsize=Matrix::Zero(); Hsize.diagonal=refSize;
+	// Hsize will contain colums with transformed base vectors FIXME : no need to initialize all the time, even if it doesn't hurt, it doesn't make sense.
+	Hsize=Matrix3r::Zero(); Hsize(0,0)=refSize[0]; Hsize(1,1)=refSize[1]; Hsize(2,2)=refSize[2]; // later with eigen: Hsize=Matrix::Zero(); Hsize.diagonal=refSize;
 	Hsize=trsf*Hsize;
 	// lengths of transformed cell vectors, skew cosines
 	Matrix3r Hnorm; // normalized transformed base vectors

=== modified file 'core/Cell.hpp'
--- core/Cell.hpp	2010-05-04 13:56:05 +0000
+++ core/Cell.hpp	2010-05-26 08:18:36 +0000
@@ -95,8 +95,8 @@
 		(Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
 		((Vector3r,refSize,Vector3r(1,1,1),"[will be overridden below]"))
 		((Matrix3r,trsf,Matrix3r::Identity(),"[will be overridden below]"))
-		((Matrix3r,velGrad,Matrix3r::Zero(),"Velocity gradient of the transformation; used in NewtonIntegrator.")),
-
+		((Matrix3r,velGrad,Matrix3r::Zero(),"Velocity gradient of the transformation; used in NewtonIntegrator."))
+		((Matrix3r,Hsize,Matrix3r::Zero(),"The current period size (one column per box edge) |yupdate|")),
 		/*ctor*/ integrateAndUpdate(0),
 
 		/*py*/

=== modified file 'pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.cpp'
--- pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.cpp	2010-04-18 17:40:36 +0000
+++ pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.cpp	2010-05-26 08:18:36 +0000
@@ -1,5 +1,5 @@
 /*************************************************************************
-*  Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxx>         *
+*  Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>         *
 *                                                                        *
 *  This program is free software; it is licensed under the terms of the  *
 *  GNU General Public License v2 or later. See file LICENSE for details. *

=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp	2010-05-20 18:39:42 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp	2010-05-26 08:18:36 +0000
@@ -86,7 +86,7 @@
 	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
 	Vector3r diff=.5*(multiplier-1)*(p2-p1);
 	setTgPlanePts(p1-diff,p2+diff);
-	return 2*diff;
+	return diff*2;
 }
 
 

=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_WallSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_WallSphere.cpp	2010-05-24 15:42:48 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_WallSphere.cpp	2010-05-26 08:18:36 +0000
@@ -43,7 +43,7 @@
 	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
 	Vector3r diff=.5*(multiplier-1)*(p2-p1);
 	setTgPlanePts(p1-diff,p2+diff);
-	return 2*diff;
+	return diff*2;
 }
 
 CREATE_LOGGER(Ig2_Wall_Sphere_Dem3DofGeom);

=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-05-20 18:39:42 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-05-26 08:18:36 +0000
@@ -1,7 +1,7 @@
 // © 2004 Olivier Galizzi <olivier.galizzi@xxxxxxx>
 // © 2004 Janek Kozicki <cosurgi@xxxxxxxxxx>
 // © 2008 Václav Šmilauer <eudoxos@xxxxxxxx>
-// © 2008 Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
+// © 2006 Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
 
 #include "ScGeom.hpp"
 #include<yade/core/Omega.hpp>
@@ -64,7 +64,12 @@
 	return shearIncrement;
 }
 
-Vector3r ScGeom::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
+Vector3r ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
+	//Just pass null shift to the periodic version
+	return rotateAndGetShear(shearForce,prevNormal,rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
+}
+
+Vector3r ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
 	Vector3r axis;
 	Real angle;
 	// approximated rotations
@@ -89,28 +94,42 @@
 	Vector3r& x = contactPoint;
 	Vector3r c1x, c2x;
 	if(avoidGranularRatcheting){
-		/* The following definition of c1x and c2x is to avoid "granular ratcheting" 
-		 *  (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN, 
-		 *  "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
-		 *  ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors) */
-
+		//FIXME : put the long comment on the wiki and keep only a small abstract and link here.
+		/* B.C. Comment : The following definition of c1x and c2x is to avoid "granular ratcheting". This phenomenon has been introduced to me by S. McNamara in a seminar help in Paris, 2004 (GDR MiDi). The concept is also mentionned in many McNamara, Hermann, Lüding, and co-workers papers (see online : "Discrete element methods for the micro-mechanical investigation of granular ratcheting", R. García-Rojo, S. McNamara, H. J. Herrmann, Proceedings ECCOMAS 2004, @ http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it refers to an accumulation of strain in cyclic loadings.
+         
+	        Unfortunately, published papers tend to focus on the "good" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while the talk of McNamara in Paris clearly mentioned a possible "bad" ratcheting, purely linked with the formulation of the contact laws in what he called "molecular dynamics" (i.e. Cundall's model, as opposed to "contact dynamics" from Moreau and Jean).
+		
+		Giving a short explanation :
+		The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
+		1. translation "dx" in the normal direction
+		2. rotation "a"
+		3. translation "-dx" (back to initial position)
+		4. rotation "-a" (back to initial orientation)
+		
+		If the branch vector used to define the relative shear in rotation*branch is not constant (typically if it is defined from the vector center->contactPoint like in the "else" below), then the shear displacement at the end of this cycle is not null : rotations a and -a are multiplied by branches of different lengths.
+		It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an inconsistent energy creation or loss. It is BAD! And given the fact that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.
+		The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
+		 */
 		// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
 		c1x =   radius1*normal; 
 		c2x =  -radius2*normal;
-	}
-	else {
+	} else {
 		// FIXME: It is correct for sphere-sphere and sphere-facet contact
 		c1x = (x - rbp1->pos);
 		c2x = (x - rbp2->pos);
 	}
 	Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
-	Vector3r shearVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
-	Vector3r shearDisplacement = shearVelocity*dt;
-	shearForce -= ks*shearDisplacement;
+	
+	//update relative velocity for interactions across periods
+	//if (scene->cell->isPeriodic) shearDisplacement+=scene->cell->velGrad*scene->cell->Hsize*cellDist; FIXME : we need pointer to scene. For now, everything will be computed outside this function, which doens't make much sense.
+	relativeVelocity+=shiftVel;
+	//keep the shear part only
+	relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
+	Vector3r shearDisplacement = relativeVelocity*dt;
+	//shearForce -= ks*shearDisplacement;//this is constitutive behaviour -> moved to Law2 functors
 	return shearDisplacement;
 }
 
-
 /* keep this for reference; declarations in the header */
 #if 0
 	Vector3r ScGeom::relRotVector() const{

=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-05-20 18:39:42 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-05-26 08:18:36 +0000
@@ -24,8 +24,10 @@
 		Vector3r updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
 		
 		virtual ~ScGeom();
-		//! update the shear force on this contact given dynamic parameters of bodies (rotate the force and add the force increment). Returns shear displacement increment, don't update shear.
-		Vector3r updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
+		//! Update the shear force orientation (or any other vector, think "shear") on this contact given dynamic parameters of bodies. Compute and return shear displacement increment.
+		Vector3r rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt,bool avoidGranularRatcheting=true);
+		//!  Periodic variant. Needs the velocity shift across periods for periodic BCs (else it is safe to pass Vector3r::Zero()). Typically obtained as scene->cell->velGrad*scene->cell->Hsize*cellDist. It would be better to define the shift transparently inside the function, but it needs scene and interaction pointers, which we don't have here.
+		Vector3r rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, 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/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp	2010-05-19 21:05:54 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp	2010-05-26 08:18:36 +0000
@@ -85,7 +85,9 @@
 		///////////////////////// CREEP START ///////////
 		if (shear_creep) shearForce -= currentContactPhysics->ks*(shearForce*dt/creep_viscosity);
 		///////////////////////// CREEP END ////////////
-		currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+		Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt);
+		//Linear elasticity giving "trial" shear force
+		shearForce -= currentContactPhysics->ks*dus;
 
 		Real Fs = currentContactPhysics->shearForce.norm();
 		Real maxFs = currentContactPhysics->shearAdhesion;

=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp	2010-05-13 20:19:38 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp	2010-05-26 08:18:36 +0000
@@ -63,8 +63,9 @@
 	// using scGeom function updateShear(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting)	
 	State* st1 = Body::byId(id1,rootBody)->state.get();
 	State* st2 = Body::byId(id2,rootBody)->state.get();
-	geom->updateShearForce(shearForce, phys->ks, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
-	
+	Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
+	//Linear elasticity giving "trial" shear force
+	shearForce -= phys->ks*dus;
 	// needed for the next timestep
 	phys->prevNormal = geom->normal;
 		

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-05-24 19:56:22 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-05-26 08:18:36 +0000
@@ -67,14 +67,18 @@
 	Real un=currentContactGeometry->penetrationDepth;
 	TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
 	currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
-	
+	//FIXME : it would make more sense to compute the shift in some Ig2->getShear, using precomputed velGrad*Hsize in Cell and I->cellDist, but it is not possible with current design (Ig doesn't know scene nor I).
+	Vector3r shiftVel = scene->isPeriodic ? (scene->cell->velGrad*scene->cell->Hsize)*Vector3r(contact->cellDist[0],contact->cellDist[1],contact->cellDist[2]) : Vector3r::Zero();
 	if (!traceEnergy){
 		//Update force but don't compute energy terms (see below))
 		if(useShear){
 			currentContactGeometry->updateShear(de1,de2,dt);
 			shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
 		} else {
-			currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);}
+			Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt,shiftVel,/*avoid ratcheting*/true);
+			//Linear elasticity giving "trial" shear force
+			shearForce -= currentContactPhysics->ks*dus;
+		}
 		// PFC3d SlipModel, is using friction angle. CoulombCriterion
 		Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
 			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
@@ -85,8 +89,9 @@
 	} 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
 		if(useShear) throw ("energy tracing not defined with useShear==true");
+		Vector3r shearDisp = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt,shiftVel,true);
 		Vector3r prevForce=shearForce;//store prev force for definition of plastic slip
-		Vector3r shearDisp = currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+		shearForce -= currentContactPhysics->ks*shearDisp;
 		Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
 			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
 		if( shearForce.squaredNorm() > maxFs ){

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp	2010-05-20 18:39:42 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp	2010-05-26 08:18:36 +0000
@@ -25,7 +25,7 @@
 		void initPlasticDissipation(Real initVal=0);
 		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` (sphere-box interactions are not implemented for the latest).",
 		((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,useShear,false,"Use ScGeom::updateShear rather than ScGeom::rotateAndGetShear for shear force computation."))
 		((bool,traceEnergy,false,"Define the total energy dissipated in plastic slips at all contacts."))
 		,,
 		.def("elasticEnergy",&Law2_ScGeom_FrictPhys_Basic::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
@@ -66,7 +66,7 @@
 		void action();
 	YADE_CLASS_BASE_DOC_ATTRS(ElasticContactLaw,GlobalEngine,"[DEPRECATED] Loop over interactions applying :yref:`Law2_ScGeom_FrictPhys_Basic` on all interactions.\n\n.. note::\n  Use :yref:`InteractionDispatchers` and :yref:`Law2_ScGeom_FrictPhys_Basic` instead of this class for performance reasons.",
 		((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 :yref:`ScGeom`::updateShear rather than :yref:`ScGeom`::updateShearForce for shear force computation."))
+		((bool,useShear,false,"Use :yref:`ScGeom`::updateShear rather than :yref:`ScGeom`::rotateAndGetShear for shear force computation."))
 	);
 };
 

=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-05-25 13:55:40 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-05-26 08:18:36 +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;
-  	scg->updateShearForce(trialFs, phys->ks, phys->prevNormal, de1, de2, dt, preventGranularRatcheting);
+  	Vector3r dus =scg->rotateAndGetShear(trialFs, phys->prevNormal, de1, de2, dt, /*define this for periodicity, see :ElasticContactLaw.cpp:71*/ Vector3r::Zero(), preventGranularRatcheting);
+	//Linear elasticity giving "trial" shear force
+	trialFs -= phys->ks*dus;
 
  
 	/*** MOHR-COULOMB LAW ***/
  Thread Previous • Date Previous • • Thread Next