Author: eudoxos
Date: 2009-03-23 20:35:09 +0100 (Mon, 23 Mar 2009)
New Revision: 1728
Modified:
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/scripts/test/triax-identical-results.py
Log:
1. Remove some garbage from SpheresContactGeometry
2. Verify that SCG_SHEAR doesn't alter behavior if ElasticContactLaw::useShear
is false
3. Implement SCG_SHEAR for sphere-box interactions
4. sphere-box interactions no longer call goReverse, but swap interaction order
instead, as facets do.
5. Fix triax-idnetical-results.py to reload generated initial config to avoid
rounding issues of sphere coords in text file.
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
2009-03-23 19:35:09 UTC (rev 1728)
@@ -10,12 +10,12 @@
SpheresContactGeometry::~SpheresContactGeometry(){};
#ifdef SCG_SHEAR
-void SpheresContactGeometry::updateShear(const RigidBodyParameters* rbp1,
const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
+Vector3r SpheresContactGeometry::updateShear(const RigidBodyParameters* rbp1,
const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
Vector3r axis;
Real angle;
- shearIncrement=Vector3r::ZERO;
+ Vector3r shearIncrement(Vector3r::ZERO);
// approximated rotations
axis = prevNormal.Cross(normal);
@@ -62,7 +62,7 @@
shearIncrement -= shearDisplacement;
shear+=shearIncrement;
- shearUpdateIter=Omega::instance().getCurrentIteration();
+ return shearIncrement;
}
#endif
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
2009-03-23 19:35:09 UTC (rev 1728)
@@ -24,20 +24,13 @@
Vector3r shear;
//! Normal of the contact in the previous step
Vector3r prevNormal;
- //! Increment of shear displacement in last step (is
already added to shear); perhaps useful for viscosity or something similar
- Vector3r shearIncrement;
- long shearUpdateIter; // debugging only, to check when
shear was updated last time
- //! update shear on this contact given dynamic
parameters of bodies. Should be called from constitutive law, exactly once per
iteration
- void updateShear(const RigidBodyParameters* rbp1, const
RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
- //const Vector3r& getShear(){
if(Omega::instance().getCurrentIteration()>shearUpdateIter) throw
runtime_error("SpheresContactGeometry::updateShear must be called prior to
getShear()."); return shear; }
+ //! update shear on this contact given dynamic
parameters of bodies. Should be called from constitutive law, exactly once per
iteration. Returns shear increment in this update, which is already added to
shear.
+ Vector3r updateShear(const RigidBodyParameters* rbp1,
const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
#endif
// all the rest here will hopefully disappear at some point...
// begin abusive storage
- bool hasNormalViscosity;
- Real NormalViscisity;
- Real NormalRelativeVelocity;
//! Whether the original contact was on the positive (1) or
negative (-1) facet side; this is to permit repulsion to the right side even if
the sphere passes, under extreme pressure, geometrically to the other facet's
side. This is used only in
InteractingFacet2IteractingSphere4SpheresContactGeometry. If at any point the
contact is with the edge or vertex instead of face, this attribute is reset so
that contact with the other face is possible.
int facetContactFace;
// end abusive storage
@@ -99,7 +92,7 @@
SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();
#ifdef SCG_SHEAR
- shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO
/*initialized to aproper value by geom functor*/;
shearIncrement=Vector3r::ZERO; shearUpdateIter=-1 /* proper value from geom
functor */;
+ shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO
/*initialized to proper value by geom functor*/;
#endif
}
virtual ~SpheresContactGeometry();
@@ -114,8 +107,6 @@
#ifdef SCG_SHEAR
(shear)
(prevNormal)
- (shearIncrement)
- (shearUpdateIter)
#endif
(facetContactFace)
// hasShear
Modified:
trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
---
trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
2009-03-22 18:08:17 UTC (rev 1727)
+++
trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
2009-03-23 19:35:09 UTC (rev 1728)
@@ -87,6 +87,11 @@
shared_ptr<SpheresContactGeometry> scm;
if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new
SpheresContactGeometry());
else scm =
YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
+
+ #ifdef SCG_SHEAR
+ if(c->isNew) { /* same as below */
scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); }
+ else {scm->prevNormal=scm->normal;}
+ #endif
// contact point is in the middle of overlapping volumes
//(in the direction of penetration, which is normal to the box
surface closest to sphere center) of overlapping volumes
@@ -129,6 +134,10 @@
shared_ptr<SpheresContactGeometry> scm;
if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new
SpheresContactGeometry());
else scm =
YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
+ #ifdef SCG_SHEAR
+ if(c->isNew) { /* same as below */
scm->prevNormal=-cOnBox_sphere; }
+ else {scm->prevNormal=scm->normal;}
+ #endif
scm->contactPoint = 0.5*(pt1+pt2);
//scm->normal = pt1-pt2; scm->normal.Normalize();
//scm->penetrationDepth = (pt1-pt2).Length();
@@ -150,19 +159,9 @@
const Se3r& se32,
const shared_ptr<Interaction>&
c)
{
- bool isInteracting = go(cm2,cm1,se32,se31,c);
- if (isInteracting)
- {
- SpheresContactGeometry* scm =
static_cast<SpheresContactGeometry*>(c->interactionGeometry.get());
- //Vector3r tmp = scm->closestsPoints[0].first;
- //scm->closestsPoints[0].first = scm->closestsPoints[0].second;
- //scm->closestsPoints[0].second = tmp;
- scm->normal = -scm->normal;
- Real tmpR = scm->radius1;
- scm->radius1 = scm->radius2;
- scm->radius2 = tmpR;
- }
- return isInteracting;
+ assert(c->isNew);
+ c->swapOrder();
+ return go(cm2,cm1,se32,se31,c);
}
YADE_PLUGIN();
Modified:
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
---
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
2009-03-22 18:08:17 UTC (rev 1727)
+++
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
2009-03-23 19:35:09 UTC (rev 1728)
@@ -165,23 +165,10 @@
const Se3r&
se32,
const
shared_ptr<Interaction>& c)
{
+ assert(c->isNew);
c->swapOrder();
//LOG_WARN("Swapped interaction order for
"<<c->getId2()<<"&"<<c->getId1());
return go(cm2,cm1,se32,se31,c);
-#if 0
- bool isInteracting = go(cm2,cm1,se32,se31,c);
- if (isInteracting)
- {
- SpheresContactGeometry* scm =
static_cast<SpheresContactGeometry*>(c->interactionGeometry.get());
- scm->normal*=-1;
- std::swap(scm->radius1,scm->radius2);
- if(hasShear){
- swap(scm->pos1,scm->pos2); swap(scm->ori1,scm->ori2);
- if(c->isNew){ swap(scm->cp1rel,scm->cp2rel);
swap(scm->d1,scm->d2); }
- }
- }
- return isInteracting;
-#endif
}
YADE_PLUGIN();
Modified:
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
---
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
2009-03-22 18:08:17 UTC (rev 1727)
+++
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
2009-03-23 19:35:09 UTC (rev 1728)
@@ -38,14 +38,8 @@
else { scm=shared_ptr<SpheresContactGeometry>(new
SpheresContactGeometry()); c->interactionGeometry=scm; }
#ifdef SCG_SHEAR
- if(c->isNew){
- scm->prevNormal=normal;
-
scm->shearUpdateIter=Omega::instance().getCurrentIteration(); /* no shear at
the very beginning; shear initialized to zero vector in SCG ctor */
- } else {
- scm->prevNormal=scm->normal;
- // make sure updateShear was properly called at
last iteration; debugging only
-
//assert(scm->shearUpdateIter==Omega::instance().getCurrentIteration()-1);
- }
+ if(c->isNew) scm->prevNormal=normal;
+ else scm->prevNormal=scm->normal;
#endif
Real penetrationDepth=s1->radius+s2->radius-normal.Normalize();
/* Normalize() works in-place and returns length before normalization; from
here, normal is unit vector */
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-22
18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-23
19:35:09 UTC (rev 1728)
@@ -123,10 +123,10 @@
currentContactGeometry->updateShear(de1,de2,dt);
shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
} else {
+ #endif
currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+ #ifdef SCG_SHEAR
}
- #else
-
currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
#endif
// PFC3d SlipModel, is using friction angle.
CoulombCriterion
Modified: trunk/scripts/test/triax-identical-results.py
===================================================================
--- trunk/scripts/test/triax-identical-results.py 2009-03-22 18:08:17 UTC
(rev 1727)
+++ trunk/scripts/test/triax-identical-results.py 2009-03-23 19:35:09 UTC
(rev 1728)
@@ -13,18 +13,18 @@
if not exists(outSph): break
i+=1
inSph='%s-in.spheres'%sph
-if exists(inSph):
- print "Using existing initial configuration",inSph
-Preprocessor('TriaxialTest',{'importFilename':(inSph if exists(inSph) else
''),'numberOfGrains':400}).load()
-if not exists(inSph):
- print "Saving new initial configuration to",inSph
+if exists(inSph): print "Using existing initial configuration",inSph
+else:
+ Preprocessor('TriaxialTest').load()
+ print "Using new initial configuration in",inSph
utils.spheresToFile(inSph)
+Preprocessor('TriaxialTest',{'importFilename':inSph}).load()
O.usesTimeStepper=False
O.dt=utils.PWaveTimeStep()
#
# uncomment this line to enable shear computation in SpheresContactGeometry
and then compare results with this line commented
#
-#[e for e in O.engines if e.name=='ElasticContactLaw'][0]['useShear']=True
+[e for e in O.engines if e.name=='ElasticContactLaw'][0]['useShear']=True
if 1:
O.run(2000,True)
utils.spheresToFile(outSph)
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
_______________________________________________
yade-dev mailing list
[email protected]
https://lists.berlios.de/mailman/listinfo/yade-dev
