yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2249
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Mon 2010-05-24 21:56:22 +0200
message:
  1 - stress definition in PeriIsoCompressor was wrong for ScGeom (sign mistake due to double-reverse);
  
  2 - applyForceAtContactPoint cannot be used in periodic + ScGeom, since this function needs the wrapped position of each body (in Dem3Dof, wrapped position is an 
  attribute, in ScGeom there is no such attribute and we don't want to define it as it costs cpu). The function call is replaced by separate applyForce/applyTorque using a 
  "branch" computed using radius, normal, and penetration depth.
  
  3 - Implement homothetic=2 properly. The acceleration of the period (i.e. change in velGrad) is applied on particles as an equivalent impulse. If velGrad is constant, 
  there is no effect on particles. If velGrad changes, each particle is accelerated with (delta_velGrad*position). This is like giving initial velocities to particles so 
  that they move homothetically.
  
  TODO :
  4 - Define relative velocity across period correctly. It should be done adding a velGrad*shift2 term, but shift2 is unknown in Law2... I'll maybe exploit the "shear" 
  attribute of ScGeom to precompute a shifted shear disp. increment in Ig2, and use it later in Law2.
  (btw : isn't that an interesting option to merge ScGeom and Dem3? Letting shear be the total OR incremental displacement?)
modified:
  pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.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/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp	2010-04-25 15:46:26 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp	2010-05-24 19:56:22 +0000
@@ -29,10 +29,8 @@
 		Real norm=normal.norm(); normal/=norm; // normal is unit vector now
 		Real penetrationDepth=s1->radius+s2->radius-norm;
 		scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
-		#ifdef SCG_SHEAR
-			if(isNew) scm->prevNormal=normal; 
-			else scm->prevNormal=scm->normal;
-		#endif
+		if(isNew) scm->prevNormal=normal; 
+		else scm->prevNormal=scm->normal;
 		scm->normal=normal;
 		scm->penetrationDepth=penetrationDepth;
 		scm->radius1=s1->radius;

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-05-24 15:42:48 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-05-24 19:56:22 +0000
@@ -98,8 +98,16 @@
 			((1/currentContactPhysics->ks)*(trialForce-shearForce))//plastic disp.
 			.dot(shearForce);//active force
 		}
-	}	
+	}
+	if (!scene->isPeriodic)
 	applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
+	else {//we need to use correct branches in the periodic case, the following apply for spheres only
+		Vector3r force = -currentContactPhysics->normalForce-shearForce;
+		ncb->forces.addForce(id1,force);
+		ncb->forces.addForce(id2,-force);
+		ncb->forces.addTorque(id1,(currentContactGeometry->radius1-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
+		ncb->forces.addTorque(id2,(currentContactGeometry->radius2-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
+	}
 	currentContactPhysics->prevNormal = currentContactGeometry->normal;
 }
 

=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2010-05-18 21:22:14 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2010-05-24 19:56:22 +0000
@@ -123,9 +123,8 @@
 		//Contact force
 		Vector3r f= ( reversedForces?-1.:1. ) * ( nsi->normalForce+nsi->shearForce );
 		//branch vector, FIXME : the first definition generalizes to non-spherical bodies but needs wrapped coords.
-		//    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 );
+		//    Vector3r branch=(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
+		Vector3r branch= gsc->normal* ( gsc->refR1+gsc->refR2 );
 		// 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];