yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2237
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Tue 2010-05-18 18:29:56 +0200
message:
  - Simplify equations inside plastic condition of Dem3Dof (1 sqrt instead of 3, less norm()), add a new function that uses a multiplier instead of maxDisp to take 
  advantage of this new formulation (Vaclav, could you review and tell if we should merge maxDisp/multiplier in one single function with a bool? I didn't want to break any 
  other part so I didn't touch the previous function).
  - Fix a wrong equation in shop::unbalanced force (max is for body force, not contact force).
  - hardcode tensor product in PeriIso, as this is instanciating a matrix and it is done for each contact at each step.
modified:
  pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
  pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp
  pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp
  pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
  pkg/dem/meta/Shop.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 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp	2010-05-11 12:25:19 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp	2010-05-18 16:29:56 +0000
@@ -78,6 +78,17 @@
 	return 2*diff.norm();
 }
 
+/*! As above : perform slip of the projected contact points. Here, we directly give the multiplier applied on the distance for faster results.
+ * The slipped distance is returned.
+ */
+Real Dem3DofGeom_SphereSphere::scaleToDisplacementTMax(Real multiplier){
+	assert(multiplier>=0 && multiplier<=1);
+	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
+	Vector3r diff=.5*(multiplier-1)*(p2-p1);
+	setTgPlanePts(p1-diff,p2+diff);
+	return 2*diff.norm();
+}
+
 
 /* Move contact point on both spheres in such way that their relative position (displacementT) is the same;
  * this should be done regularly to ensure that the angle doesn't go over π, since then quaternion would

=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp	2010-04-25 15:46:26 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp	2010-05-18 16:29:56 +0000
@@ -28,6 +28,7 @@
 			#endif
 		}
 		Real slipToDisplacementTMax(Real displacementTMax);
+		Real scaleToDisplacementTMax(Real multiplier);
 		/********* end API ***********/
 
 	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(Dem3DofGeom_SphereSphere,Dem3DofGeom,"Class representing 2 spheres in contact which computes 3 degrees of freedom (normal and shear deformation).",

=== modified file 'pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp	2010-03-10 20:30:46 +0000
+++ pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp	2010-05-18 16:29:56 +0000
@@ -30,6 +30,7 @@
 		virtual Real displacementN();
 		virtual Vector3r displacementT(){throw;}
 		virtual Real slipToDisplacementTMax(Real displacementTMax){throw;}; // plasticity
+		virtual Real scaleToDisplacementTMax(Real multiplier){throw;}; // plasticity (variant using multiplier dispMax/disp)
 		// reference radii, for contact stiffness computation; set to negative for nonsense values
 		// end API
 

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-05-18 10:28:02 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-05-18 16:29:56 +0000
@@ -112,7 +112,12 @@
 	phys->normalForce=phys->kn*displN*geom->normal;
 	Real maxFsSq=phys->normalForce.squaredNorm()*pow(phys->tangensOfFrictionAngle,2);
 	Vector3r trialFs=phys->ks*geom->displacementT();
-	if(trialFs.squaredNorm()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)/phys->ks); trialFs*=sqrt(maxFsSq/(trialFs.squaredNorm()));}
+	Real multiplier=maxFsSq/trialFs.squaredNorm();
+// 	if(trialFs.squaredNorm()>maxFsSq){
+// 		geom->slipToDisplacementTMax(sqrt(maxFsSq)/phys->ks); trialFs*=sqrt(maxFsSq/(trialFs.squaredNorm()));}
+	if(multiplier<1){
+		multiplier = sqrt(multiplier);
+		geom->scaleToDisplacementTMax(multiplier); trialFs*=multiplier;}
 	phys->shearForce=trialFs;
 	applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,scene);
 }

=== modified file 'pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp'
--- pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp	2010-05-18 10:28:02 +0000
+++ pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp	2010-05-18 16:29:56 +0000
@@ -32,7 +32,6 @@
 		void findTimeStepFromBody(const shared_ptr<Body>& body, Scene * ncb);
 	
 	public :
-		//! used as a multiplier on the theoretical critical timestep (compensate some approximations in the computation)
 		virtual ~GlobalStiffnessTimeStepper();
 	
 		virtual void computeTimeStep(Scene*);

=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2010-05-11 12:38:36 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2010-05-18 16:29:56 +0000
@@ -126,8 +126,11 @@
 		//    Vector3r branch=
 		//   (reversedForces?-1.:1.)*(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
 		Vector3r branch= ( reversedForces?-1.:1. ) *gsc->normal* ( gsc->refR1+gsc->refR2 );
-		// tensorial product f*branch
-		stressTensor+=makeTensorProduct( f, branch );
+		// tensorial product f*branch (hand-write the tensor product to prevent matrix instanciation inside the loop by makeTensorProduct)
+		stressTensor(0,0)+=f(0)*branch(0); stressTensor(1,0)+=f(1)*branch(0); stressTensor(2,0)+=f(2)*branch(0);
+		stressTensor(0,1)+=f(0)*branch(1); stressTensor(1,1)+=f(1)*branch(1); stressTensor(2,1)+=f(2)*branch(1);
+		stressTensor(0,2)+=f(0)*branch(2); stressTensor(1,2)+=f(1)*branch(2); stressTensor(2,2)+=f(2)*branch(2);
+		//stressTensor+=makeTensorProduct( f, branch );
 		if( !dynCell )
 		{
 			for ( int i=0; i<3; i++ ) sumStiff[i]+=abs ( gsc->normal[i] ) *nsi->kn+ ( 1-abs ( gsc->normal[i] ) ) *nsi->ks;

=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp	2010-05-18 10:28:02 +0000
+++ pkg/dem/meta/Shop.cpp	2010-05-18 16:29:56 +0000
@@ -168,20 +168,21 @@
 	Scene* rb=_rb ? _rb : Omega::instance().getScene().get();
 	rb->forces.sync();
 	// get maximum force on a body and sum of all forces (for averaging)
-	Real sumF=0,maxF=0,currF;
+	Real sumF=0,maxF=0,currF; int nb=0;
 	FOREACH(const shared_ptr<Body>& b, *rb->bodies){
 		if(!b->isDynamic) continue;
-		currF=rb->forces.getForce(b->id).norm(); maxF=max(currF,maxF); sumF+=currF;
+		currF=rb->forces.getForce(b->id).norm(); maxF=max(currF,maxF); sumF+=currF; nb++;
 	}
-	Real meanF=sumF/rb->bodies->size(); 
+	Real meanF=sumF/nb;
 	// get max force on contacts
-	Real maxContactF=0;
+	Real maxContactF=0; sumF=0; nb=0;
 	FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
 		if(!I->isReal()) continue;
 		shared_ptr<NormShearPhys> nsi=YADE_PTR_CAST<NormShearPhys>(I->interactionPhysics); assert(nsi);
-		maxContactF=max(maxContactF,(nsi->normalForce+nsi->shearForce).norm());
+		sumF+=(nsi->normalForce+nsi->shearForce).norm(); nb++;
 	}
-	return (useMaxForce?maxF:meanF)/maxContactF;
+	sumF/=nb;
+	return (useMaxForce?maxF:meanF)/(sumF);
 }
 
 Real Shop::kineticEnergy(Scene* _rb){