------------------------------------------------------------
revno: 2249
committer: Bruno Chareyre <bchare...@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];
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
