yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01202
[svn] r1757 - in trunk: extra pkg/common/DataClass/InteractingGeometry pkg/dem pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/EngineUnit pkg/dem/Engine/StandAloneEngine scripts/test
Author: eudoxos
Date: 2009-04-22 08:57:45 +0200 (Wed, 22 Apr 2009)
New Revision: 1757
Added:
trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp
Modified:
trunk/extra/Brefcom.cpp
trunk/extra/Brefcom.hpp
trunk/extra/SConscript
trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp
trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp
trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp
trunk/pkg/dem/SConscript
trunk/scripts/test/facet-topo.py
Log:
1. Finish implementation of Dof3DofGeom_FacetSphere (with plastic slip as well)
2. Adapt Brefcom optionally to Dem3DofGeom
3. Add d0fixup to SpheresContactGeometry to get correct normal strain in sphere-facet contact with fictive sphere
(zero for sphere-sphere)
4. Finish FacetTopologyAnalyzer; angle usage in Dem4DofGeom_FacetSphere not yet tested.
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/extra/Brefcom.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -5,6 +5,7 @@
#include<yade/pkg-common/Sphere.hpp>
#include<yade/lib-QGLViewer/qglviewer.h>
#include<yade/lib-opengl/GLUtils.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
YADE_PLUGIN("BrefcomMakeContact","BrefcomContact","BrefcomLaw","GLDrawBrefcomContact","BrefcomDamageColorizer", "BrefcomPhysParams", "BrefcomGlobalCharacteristics", "ef2_Spheres_Brefcom_BrefcomLaw" /* ,"BrefcomStiffnessComputer"*/ );
@@ -66,7 +67,11 @@
//! @todo Formulas in the following should be verified
void BrefcomMakeContact::go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction){
- const shared_ptr<SpheresContactGeometry>& contGeom=YADE_PTR_CAST<SpheresContactGeometry>(interaction->interactionGeometry);
+ #ifdef BREFCOM_DEM3DOF
+ const Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
+ #else
+ const SpheresContactGeometry* contGeom=YADE_CAST<SpheresContactGeometry*>(interaction->interactionGeometry.get());
+ #endif
assert(contGeom); // for now, don't handle anything other than SpheresContactGeometry
if(!interaction->isNew && interaction->interactionPhysics){ /* relax */ }
@@ -78,7 +83,12 @@
Real E12=2*elast1->young*elast2->young/(elast1->young+elast2->young); // harmonic Young's modulus average
//Real nu12=2*elast1->poisson*elast2->poisson/(elast1->poisson+elast2->poisson); // dtto for Poisson ratio
- Real S12=Mathr::PI*pow(min(contGeom->radius1,contGeom->radius2),2); // "surface" of interaction
+ #ifdef BREFCOM_DEM3DOF
+ Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
+ Real S12=Mathr::PI*pow(minRad,2); // "surface" of interaction
+ #else
+ Real S12=Mathr::PI*pow(min(contGeom->radius1,contGeom->radius2),2); // "surface" of interaction
+ #endif
//Real E=(E12 /* was here for Kn: *S12/d0 */)*((1+alpha)/(beta*(1+nu12)+gamma*(1-alpha*nu12)));
//Real E=E12; // apply alpha, beta, gamma: garbage values of E !?
@@ -189,7 +199,12 @@
void ef2_Spheres_Brefcom_BrefcomLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
//timingDeltas->start();
- SpheresContactGeometry* contGeom=static_cast<SpheresContactGeometry*>(_geom.get());
+ #ifdef BREFCOM_DEM3DOF
+ Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
+ #else
+ SpheresContactGeometry* contGeom=static_cast<SpheresContactGeometry*>(_geom.get());
+ assert(contGeom->hasShear);
+ #endif
BrefcomContact* BC=static_cast<BrefcomContact*>(_phys.get());
/* kept fully damaged contacts; note that normally the contact is deleted _after_ the BREFCOM_MATERIAL_MODEL,
@@ -206,11 +221,17 @@
#define NNAN(a) assert(!isnan(a));
#define NNANV(v) assert(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
- assert(contGeom->hasShear);
//timingDeltas->checkpoint("setup");
- epsN=contGeom->epsN(); epsT=contGeom->epsT();
+ #ifdef BREFCOM_DEM3DOF
+ epsN=contGeom->strainN(); epsT=contGeom->strainT();
+ #else
+ epsN=contGeom->epsN(); epsT=contGeom->epsT();
+ #endif
if(isnan(epsN)){
- LOG_ERROR("d0="<<contGeom->d0<<", d1,d2="<<contGeom->d1<<","<<contGeom->d2<<"; pos1,pos2="<<contGeom->pos1<<","<<contGeom->pos2);
+ #ifndef BREFCOM_DEM3DOF
+ LOG_ERROR("d0="<<contGeom->d0<<", d1,d2="<<contGeom->d1<<","<<contGeom->d2<<"; pos1,pos2="<<contGeom->pos1<<","<<contGeom->pos2);
+ #endif
+ throw;
}
NNAN(epsN); NNANV(epsT);
// already in SpheresContactGeometry:
@@ -237,7 +258,11 @@
Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
Fs=sigmaT*crossSection; BC->shearForce=Fs;
- applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->pos1, I->getId2(), contGeom->pos2, rootBody);
+ #ifdef BREFCOM_DEM3DOF
+ applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position, rootBody);
+ #else
+ applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->pos1, I->getId2(), contGeom->pos2, rootBody);
+ #endif
//timingDeltas->checkpoint("rest");
}
Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/extra/Brefcom.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -174,6 +174,8 @@
};
REGISTER_SERIALIZABLE(BrefcomPhysParams);
+//#define BREFCOM_DEM3DOF
+
class ef2_Spheres_Brefcom_BrefcomLaw: public ConstitutiveLaw{
public:
/*! Damage evolution law */
@@ -184,7 +186,11 @@
bool logStrain;
ef2_Spheres_Brefcom_BrefcomLaw(): logStrain(false){ /*timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/ }
void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
- FUNCTOR2D(SpheresContactGeometry,BrefcomContact);
+ #ifdef BREFCOM_DEM3DOF
+ FUNCTOR2D(Dem3DofGeom,BrefcomContact);
+ #else
+ FUNCTOR2D(SpheresContactGeometry,BrefcomContact);
+ #endif
REGISTER_CLASS_AND_BASE(ef2_Spheres_Brefcom_BrefcomLaw,ConstitutiveLaw);
REGISTER_ATTRIBUTES(ConstitutiveLaw,(logStrain));
DECLARE_LOGGER;
Modified: trunk/extra/SConscript
===================================================================
--- trunk/extra/SConscript 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/extra/SConscript 2009-04-22 06:57:45 UTC (rev 1757)
@@ -42,9 +42,9 @@
env.SharedLibrary('UniaxialStrainControlledTest',['usct/UniaxialStrainControlledTest.cpp'],LIBS=env['LIBS']+['Shop','GlobalStiffnessTimeStepper','Brefcom','PositionOrientationRecorder']),
- env.SharedLibrary('Brefcom',['Brefcom.cpp'],CXXFLAGS=env['CXXFLAGS']+brefcomInclude,LIBS=env['LIBS']+['Shop','InteractingSphere2InteractingSphere4SpheresContactGeometry']),
+ env.SharedLibrary('Brefcom',['Brefcom.cpp'],CXXFLAGS=env['CXXFLAGS']+brefcomInclude,LIBS=env['LIBS']+['Shop','InteractingSphere2InteractingSphere4SpheresContactGeometry','DemXDofGeom']),
- env.SharedLibrary('BrefcomTestGen',['BrefcomTestGen.cpp'],LIBS=env['LIBS']+['Shop','UniaxialStrainControlledTest','PositionOrientationRecorder']),
+ env.SharedLibrary('BrefcomTestGen',['BrefcomTestGen.cpp'],LIBS=env['LIBS']+['Shop','UniaxialStrainControlledTest','PositionOrientationRecorder','Brefcom']),
env.SharedLibrary('SimpleScene',['SimpleScene.cpp'],LIBS=env['LIBS']+['Shop','SimpleElasticRelationships']),
Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -12,7 +12,7 @@
createIndex();
#ifdef FACET_TOPO
edgeAdjIds.resize(3,Body::ID_NONE);
- edgeAdjAngle.resize(3,0);
+ edgeAdjHalfAngle.resize(3,0);
#endif
}
@@ -26,7 +26,7 @@
REGISTER_ATTRIBUTE(vertices);
#ifdef FACET_TOPO
REGISTER_ATTRIBUTE(edgeAdjIds);
- REGISTER_ATTRIBUTE(edgeAdjAngle);
+ REGISTER_ATTRIBUTE(edgeAdjHalfAngle);
#endif
}
Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -41,8 +41,8 @@
#ifdef FACET_TOPO
//! facet id's that are adjacent to respective edges
vector<body_id_t> edgeAdjIds;
- //! angle between normals of this facet and the adjacent facet
- vector<body_id_t> edgeAdjAngle;
+ //! half angle between normals of this facet and the adjacent facet
+ vector<Real> edgeAdjHalfAngle;
#endif
protected:
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -7,7 +7,6 @@
Dem3DofGeom_FacetSphere::~Dem3DofGeom_FacetSphere(){}
void Dem3DofGeom_FacetSphere::setTgPlanePts(const Vector3r& p1new, const Vector3r& p2new){
- // FIXME: not tested yet
TRVAR3(cp1pt,cp2rel,contPtInTgPlane2()-contPtInTgPlane1());
cp1pt=se31.orientation.Conjugate()*(turnPlanePt(p1new,normal,se31.orientation*localFacetNormal)+contactPoint-se31.position);
cp2rel=se32.orientation.Conjugate()*Dem3DofGeom_SphereSphere::rollPlanePtToSphere(p2new,effR2,-normal);
@@ -21,6 +20,19 @@
}
}
+Real Dem3DofGeom_FacetSphere::slipToDisplacementTMax(Real displacementTMax){
+ //FIXME: not yet tested
+ // negative or zero: reset shear
+ if(displacementTMax<=0.){ setTgPlanePts(Vector3r(0,0,0),Vector3r(0,0,0)); return displacementTMax;}
+ // otherwise
+ Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
+ Real currDistSq=(p2-p1).SquaredLength();
+ if(currDistSq<pow(displacementTMax,2)) return 0; // close enough, no slip needed
+ //Vector3r diff=.5*(sqrt(currDistSq)/displacementTMax-1)*(p2-p1); setTgPlanePts(p1+diff,p2-diff);
+ Real scale=displacementTMax/sqrt(currDistSq); setTgPlanePts(scale*p1,scale*p2);
+ return (displacementTMax/scale)*(1-scale);
+}
+
CREATE_LOGGER(ef2_Facet_Sphere_Dem3DofGeom);
bool ef2_Facet_Sphere_Dem3DofGeom::go(const shared_ptr<InteractingGeometry>& cm1, const shared_ptr<InteractingGeometry>& cm2, const Se3r& se31, const Se3r& se32, const shared_ptr<Interaction>& c){
InteractingFacet* facet=static_cast<InteractingFacet*>(cm1.get());
@@ -49,12 +61,29 @@
normal.Normalize();
} else { // contact with the edge
contactPt+=edgeNormals[edgeMax]*(inCircleR-distMax);
+ bool noVertexContact=false;
//TRVAR3(edgeNormals[edgeMax],inCircleR,distMax);
// contact with vertex no. edgeMax
- if (contactPt.Dot(edgeNormals[(edgeMax-1)%3])>inCircleR) contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
+ // FIXME: this is the original version, but why (edgeMax-1)%3? IN that case, edgeNormal to edgeMax would never be tried
+ // if (contactPt.Dot(edgeNormals[ (edgeMax-1)%3])>inCircleR) contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
+ if (contactPt.Dot(edgeNormals[ edgeMax ])>inCircleR) contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
// contact with vertex no. edgeMax+1
else if(contactPt.Dot(edgeNormals[edgeMax=(edgeMax+1)%3])>inCircleR) contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
+ // contact with edge no. edgeMax
+ else noVertexContact=true;
normal=contactLine-contactPt;
+ #ifdef FACET_TOPO
+ if(noVertexContact && facet->edgeAdjIds[edgeMax]!=Body::ID_NONE){
+ // find angle between our normal and the facet's normal (still local coords)
+ Quaternionr q; q.Align(facet->nf,normal); Vector3r axis; Real angle; q.ToAxisAngle(axis,angle);
+ assert(angle>=0 && angle<=Mathr::PI);
+ if(edgeNormals[edgeMax].Dot(axis)<0) angle*=-1.;
+ bool negFace=normal.Dot(facet->nf)<0; // contact in on the negative facet's face
+ Real halfAngle=(negFace?-1.:1.)*facet->edgeAdjHalfAngle[edgeMax];
+ if(halfAngle<0 && angle>halfAngle) return false; // on concave boundary, and if in the other facet's sector, no contact
+ // otherwise the contact will be created
+ }
+ #endif
//TRVAR4(contactLine,contactPt,normal,normal.Length());
//TRVAR3(se31.orientation*contactLine,se31.position+se31.orientation*contactPt,se31.orientation*normal);
penetrationDepth=sphereRadius-normal.Normalize();
@@ -71,6 +100,7 @@
fs=shared_ptr<Dem3DofGeom_FacetSphere>(new Dem3DofGeom_FacetSphere());
c->interactionGeometry=fs;
fs->effR2=sphereRadius-penetrationDepth;
+ fs->refR1=-1; fs->refR2=sphereRadius;
fs->refLength=fs->effR2;
fs->cp1pt=contactPt; // facet-local intial contact point
fs->localFacetNormal=normal;
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -16,7 +16,7 @@
/******* API ********/
virtual Real displacementN(){ return (se32.position-contactPoint).Length()-refLength;}
virtual Vector3r displacementT(){ relocateContactPoints(); return contPtInTgPlane2()-contPtInTgPlane1(); }
- virtual Real slipToDisplacementTMax(Real displacementTMax){ LOG_FATAL("Not implemented yet."); throw; }
+ virtual Real slipToDisplacementTMax(Real displacementTMax);
/***** end API ******/
void setTgPlanePts(const Vector3r&, const Vector3r&);
@@ -26,13 +26,11 @@
Vector3r cp1pt;
//! orientation between +x and the reference contact point (on the sphere) in sphere-local coords
Quaternionr cp2rel;
- //! positions and orientations of both bodies; updated at every iteration
- Se3r se31, se32;
//! unit normal of the facet plane in facet-local coordinates
Vector3r localFacetNormal;
// effective radius of sphere
Real effR2;
- REGISTER_ATTRIBUTES(Dem3DofGeom,(cp1pt)(cp2rel)(se31)(se32)(localFacetNormal)(effR2) );
+ REGISTER_ATTRIBUTES(Dem3DofGeom,(cp1pt)(cp2rel)(localFacetNormal)(effR2) );
REGISTER_CLASS_AND_BASE(Dem3DofGeom_FacetSphere,Dem3DofGeom);
DECLARE_LOGGER;
friend class GLDraw_Dem3DofGeom_FacetSphere;
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -161,6 +161,7 @@
c->interactionGeometry=ss;
// constants
ss->refLength=dist;
+ ss->refR1=s1->radius; ss->refR2=s2->radius;
Real penetrationDepth=s1->radius+s2->radius-ss->refLength;
ss->effR1=s1->radius-.5*penetrationDepth; ss->effR2=s2->radius-.5*penetrationDepth;
// for bending only: ss->initRelOri12=se31.orientation.Conjugate()*se32.orientation;
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -13,8 +13,6 @@
void relocateContactPoints();
void relocateContactPoints(const Vector3r& tgPlanePt1, const Vector3r& tgPlanePt2);
public:
- //! Positions and orientations of both spheres; must be updated at every iteration by the geom functor
- Se3r se31, se32;
//! relative orientation of the contact point with regards to sphere-local +x axis (quasi-constant)
Quaternionr cp1rel, cp2rel;
//! shorthands
@@ -37,7 +35,7 @@
Real slipToDisplacementTMax(Real displacementTMax);
/********* end API ***********/
- REGISTER_ATTRIBUTES(Dem3DofGeom,(se31)(se32)(effR1)(effR2)(cp1rel)(cp2rel));
+ REGISTER_ATTRIBUTES(Dem3DofGeom,(effR1)(effR2)(cp1rel)(cp2rel));
REGISTER_CLASS_AND_BASE(Dem3DofGeom_SphereSphere,Dem3DofGeom);
friend class GLDraw_Dem3DofGeom_SphereSphere;
friend class ef2_Sphere_Sphere_Dem3DofGeom;
Added: trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -0,0 +1,3 @@
+#include"DemXDofGeom.hpp"
+YADE_PLUGIN("Dem3DofGeom");
+Real Dem3DofGeom::displacementN(){throw;}
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -15,11 +15,16 @@
Vector3r contactPoint;
//! make strain go to -∞ for length going to zero
bool logCompression;
+ //! se3 of both bodies (needed to compute torque from the contact, strains etc). Must be updated at every step.
+ Se3r se31, se32;
+ //! reference radii of particles (for initial contact stiffness computation)
+ Real refR1, refR2;
// API that needs to be implemented in derived classes
- virtual Real displacementN(){throw;}
+ virtual Real displacementN();
virtual Vector3r displacementT(){throw;}
virtual Real slipToDisplacementTMax(Real displacementTMax){throw;}; // plasticity
+ // reference radii, for contact stiffness computation; set to negative for nonsense values
// end API
Real strainN(){
@@ -31,7 +36,7 @@
Real slipToStrainTMax(Real strainTMax){return slipToDisplacementTMax(strainTMax*refLength)/refLength;}
REGISTER_CLASS_AND_BASE(Dem3DofGeom,InteractionGeometry);
- REGISTER_ATTRIBUTES(InteractionGeometry,(refLength)(normal)(contactPoint));
+ REGISTER_ATTRIBUTES(InteractionGeometry,(refLength)(normal)(contactPoint)(se31)(se32)(refR1)(refR2));
};
REGISTER_SERIALIZABLE(Dem3DofGeom);
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -49,7 +49,8 @@
// interaction "radii" and total length; this is _really_ constant throughout the interaction life
// d1 is really distance from the sphere1 center to the contact plane, it may not correspond to the sphere radius!
// therefore, d1+d2=d0 (distance at which the contact was created)
- Real d1, d2, d0;
+ // d0fixup is added to d0 when computing normal strain; it should fix problems with sphere-facet interactions never getting enough compressed
+ Real d1, d2, d0, d0fixup;
// initial relative orientation, used for bending and torsion computation
Quaternionr initRelOri12;
@@ -64,7 +65,7 @@
Vector3r contPt() const {return contactPoint; /*pos1+(d1/d0)*(pos2-pos1)*/}
Real displacementN() const {assert(hasShear); return (pos2-pos1).Length()-d0;}
- Real epsN() const {return displacementN()*(1./d0);}
+ Real epsN() const {return displacementN()*(1./(d0+d0fixup));}
Vector3r displacementT() { assert(hasShear);
// enabling automatic relocation decreases overall simulation speed by about 3%
// perhaps: bool largeStrains ... ?
@@ -76,11 +77,11 @@
return contPtInTgPlane2()-contPtInTgPlane1();
#endif
}
- Vector3r epsT() {return displacementT()*(1./d0);}
+ Vector3r epsT() {return displacementT()*(1./(d0+d0fixup));}
Real slipToDisplacementTMax(Real displacementTMax);
//! slip to epsTMax if current epsT is greater; return the relative slip magnitude
- Real slipToEpsTMax(Real epsTMax){ return (1/d0)*slipToDisplacementTMax(epsTMax*d0);}
+ Real slipToStrainTMax(Real epsTMax){ return (1/d0)*slipToDisplacementTMax(epsTMax*d0);}
void relocateContactPoints();
void relocateContactPoints(const Vector3r& tgPlanePt1, const Vector3r& tgPlanePt2);
@@ -90,7 +91,7 @@
Vector3r relRotVector() const;
- 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();
+ 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),d0fixup(0),initRelOri12(Quaternionr::IDENTITY){createIndex();
#ifdef SCG_SHEAR
shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO /*initialized to proper value by geom functor*/;
#endif
@@ -120,6 +121,7 @@
(d1)
(d2)
(d0)
+ (d0fixup)
(initRelOri12));
REGISTER_CLASS_AND_BASE(SpheresContactGeometry,InteractionGeometry);
REGISTER_CLASS_INDEX(SpheresContactGeometry,InteractionGeometry);
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -95,7 +95,6 @@
if (bm<icr) // contact with facet's surface
{
- assert(contactFace!=0);
penetrationDepth = sphereRadius - L;
normal.Normalize();
}
@@ -145,6 +144,7 @@
if(c->isNew){
scm->d0=scm->radius1+scm->radius2-penetrationDepth;
scm->d1=scm->radius1-.5*penetrationDepth; scm->d2=scm->radius2-.5*penetrationDepth;
+ scm->d0fixup=-scm->radius1-.5*penetrationDepth;
// quasi-constants
scm->cp1rel.Align(Vector3r::UNIT_X,se31.orientation.Conjugate()*normal);
scm->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normal));
Modified: trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -40,7 +40,7 @@
abs(vi->pos[2]-vj->pos[2])<=tolerance &&
(vi->pos-vj->pos).SquaredLength()<=tolerance*tolerance){
// OK, these two vertices touch
- LOG_TRACE("Vertices "<<vi->id<<"/"<<vi->vertexNo<<" and "<<vj->id<<"/"<<vj->vertexNo<<" close enough.");
+ // LOG_TRACE("Vertices "<<vi->id<<"/"<<vi->vertexNo<<" and "<<vj->id<<"/"<<vj->vertexNo<<" close enough.");
// add vertex to the nextIndetical of the one that has lower index; the one that is added will have isLowestIndex=false
if(vi->index<vj->index){ vi->nextIdentical.push_back(vj); vj->isLowestIndex=false; }
else{ vj->nextIdentical.push_back(vi); vi->isLowestIndex=false; }
@@ -81,11 +81,18 @@
}
topo1.push_back(t);
}
+ std::sort(topo1.begin(),topo1.end(),FacetTopology::MinVertexComparator());
size_t nTopo=topo1.size();
for(size_t i=0; i<nTopo; i++){
size_t j=i;
+ const shared_ptr<FacetTopology>& ti(topo1[i]);
+ long tiMaxVertex=ti->maxVertex();
while(++j<nTopo){
- const shared_ptr<FacetTopology>& ti(topo1[i]), &tj(topo1[j]);
+ const shared_ptr<FacetTopology> &tj(topo1[j]);
+ /* since facets are sorted by their min vertex id,
+ we know that it is safe to skip all the rest
+ as soon as max vertex of ti one is smaller than min vertex of tj, as i<=j */
+ if(tj->minVertex()>tiMaxVertex) break;
vector<size_t> vvv; // array of common vertices
for(size_t k=0; k<3; k++){
if (ti->vertices[k]==tj->vertices[0]) vvv.push_back(ti->vertices[k]);
@@ -93,9 +100,9 @@
else if(ti->vertices[k]==tj->vertices[2]) vvv.push_back(ti->vertices[k]);
}
if(vvv.size()<2) continue;
- assert(vvv.size()!=3); // same coords? nonsense
- assert(vvv.size()==2);
- vector<int> edge(2,0);
+ assert(vvv.size()!=3 /* same coords? nonsense*/ ); assert(vvv.size()==2);
+ // from here, we know ti and tj are adjacent
+ vector<int> edge(2,0); int &ei(edge[0]),&ej(edge[1]);
// identify what edge are we at, for both facets
for(int k=0; k<2; k++){
for(edge[k]=0; edge[k]<3; edge[k]++){
@@ -109,10 +116,25 @@
}
}
// add adjacency information to the facet itself
- YADE_PTR_CAST<InteractingFacet>((*rb->bodies)[ti->id]->interactingGeometry)->edgeAdjIds[edge[0]]=tj->id;
- YADE_PTR_CAST<InteractingFacet>((*rb->bodies)[tj->id]->interactingGeometry)->edgeAdjIds[edge[1]]=ti->id;
+ InteractingFacet *f1=YADE_CAST<InteractingFacet*>((*rb->bodies)[ti->id]->interactingGeometry.get()), *f2=YADE_CAST<InteractingFacet*>((*rb->bodies)[tj->id]->interactingGeometry.get());
+ f1->edgeAdjIds[ei]=ti->id; f2->edgeAdjIds[ej]=tj->id;
+ // normals are in the sense of vertex rotation (right-hand rule); therefore, if vertices of the adjacent edge are opposite on each facet, normals are in the same direction
+ bool invNormals=(ti->vertices[ei]==tj->vertices[ej]);
+ assert(
+ ( invNormals && (ti->vertices[ ei ]==tj->vertices[ej]) && (ti->vertices[(ei+1)%3]==tj->vertices[(ej+1)%3]) ) ||
+ (!invNormals && (ti->vertices[(ei+1)%3]==tj->vertices[ej]) && (ti->vertices[ ei ]==tj->vertices[(ej+1)%3]) ));
+ // angle between normals
+ PhysicalParameters *pp1=(*rb->bodies)[ti->id]->physicalParameters.get(), *pp2=(*rb->bodies)[tj->id]->physicalParameters.get();
+ Vector3r n1g=pp1->se3.orientation*f1->nf, n2g=pp2->se3.orientation*f2->nf;
+ //TRVAR2(n1g,n2g);
+ Vector3r contEdge1g=pp1->se3.orientation*(f1->vertices[(ei+1)%3]-f1->vertices[ei]); // vector of the edge of contact in global coords
+ Quaternionr q12; q12.Align(n1g,(invNormals?-1.:1.)*n2g); Real halfAngle; Vector3r axis; q12.ToAxisAngle(axis,halfAngle); halfAngle*=.5;
+ assert(halfAngle>=0 && halfAngle<=Mathr::HALF_PI);
+ if(axis.Dot(contEdge1g)<0 /* convex contact from the side of +n1 */ ) halfAngle*=-1.;
+ f1->edgeAdjHalfAngle[ei]=halfAngle;
+ f2->edgeAdjHalfAngle[ej]=(invNormals ? -halfAngle : halfAngle);
commonEdgesFound++;
- LOG_TRACE("Added adjacency information for #"<<ti->id<<"+#"<<tj->id<<" (common edges "<<edge[0]<<"+"<<edge[1]<<")");
+ LOG_TRACE("Found adjacent #"<<ti->id<<"+#"<<tj->id<<"; common edges "<<ei<<"+"<<ej<<", normals "<<(invNormals?"inversed":"congruent")<<", halfAngle "<<halfAngle<<")");
}
}
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp 2009-04-22 06:57:45 UTC (rev 1757)
@@ -39,6 +39,11 @@
long vertices[3];
//! facet id, for back reference
body_id_t id;
+ long minVertex(){return min(vertices[0],min(vertices[1],vertices[2]));}
+ long maxVertex(){return max(vertices[0],max(vertices[1],vertices[2]));}
+ struct MinVertexComparator{
+ bool operator()(const shared_ptr<FacetTopology>& t1, const shared_ptr<FacetTopology>& t2){ return t1->minVertex()<t2->minVertex();}
+ };
};
public:
//! Axis along which to do the initial vertex sort
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/pkg/dem/SConscript 2009-04-22 06:57:45 UTC (rev 1757)
@@ -26,10 +26,11 @@
])
env.Install('$PREFIX/lib/yade$SUFFIX/pkg-dem',[
+
+ env.SharedLibrary('DemXDofGeom',['DataClass/InteractionGeometry/DemXDofGeom.cpp']),
+ env.SharedLibrary('Dem3DofGeom_SphereSphere',['DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'],LIBS=env['LIBS']+['RigidBodyParameters','InteractingSphere','yade-opengl','DemXDofGeom']),
+ env.SharedLibrary('Dem3DofGeom_FacetSphere',['DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp'],LIBS=env['LIBS']+['RigidBodyParameters','InteractingSphere','InteractingFacet','yade-opengl','Dem3DofGeom_SphereSphere','DemXDofGeom']),
- env.SharedLibrary('Dem3DofGeom_SphereSphere',['DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'],LIBS=env['LIBS']+['RigidBodyParameters','InteractingSphere','yade-opengl']),
- env.SharedLibrary('Dem3DofGeom_FacetSphere',['DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp'],LIBS=env['LIBS']+['RigidBodyParameters','InteractingSphere','InteractingFacet','yade-opengl','Dem3DofGeom_SphereSphere']),
-
env.SharedLibrary('FacetTopologyAnalyzer',['Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp'],LIBS=env['LIBS']+['InteractingFacet']),
Modified: trunk/scripts/test/facet-topo.py
===================================================================
--- trunk/scripts/test/facet-topo.py 2009-04-17 15:11:23 UTC (rev 1756)
+++ trunk/scripts/test/facet-topo.py 2009-04-22 06:57:45 UTC (rev 1757)
@@ -1,5 +1,5 @@
import yade.log
-#yade.log.setLevel('FacetTopologyAnalyzer',yade.log.TRACE)
+yade.log.setLevel('FacetTopologyAnalyzer',yade.log.TRACE)
# Note: FacetTopologyAnalyzer is normally run as an initializer;
# it is only for testing sake that it is in O.engines here.
@@ -24,6 +24,17 @@
assert(topo['commonEdgesFound']==1)
if 1:
O.bodies.clear()
+ O.bodies.append([
+ utils.facet([(0,0,0),(1,0,0),(0,1,0)]),
+ utils.facet([(1,1,1),(1,0,0),(0,1,0)]),
+ ])
+ O.step()
+ assert(O.bodies[0].phys['edgeAdjIds'][1]==1 and O.bodies[1].phys['edgeAdjIds'][0]==1)
+ assert(topo['commonEdgesFound']==1)
+ assert(abs(O.bodies[0].mold['edgeAdjHalfAngle'][1]-(-5*atan(2/sqrt(2))))<1e-6)
+
+if 1:
+ O.bodies.clear()
r=.5 # radius of the sphere
nPoly=12; # try 128, it is still quite fast
def sphPt(i,j):
