yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01138
[svn] r1721 - in trunk: core lib pkg/dem pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/EngineUnit pkg/snow/Engine scripts/test
Author: eudoxos
Date: 2009-03-17 16:14:49 +0100 (Tue, 17 Mar 2009)
New Revision: 1721
Added:
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp
trunk/scripts/test/Dem3DofGeom.py
Modified:
trunk/core/Collider.cpp
trunk/lib/SConscript
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/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp
trunk/pkg/dem/SConscript
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp
Log:
1. Preliminary version for sphere-facet collisions with total formulation (moderately tested).
2. Same for sphere-sphere collision.
3. test script for facet-sphere collision geometry.
Modified: trunk/core/Collider.cpp
===================================================================
--- trunk/core/Collider.cpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/core/Collider.cpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -33,6 +33,7 @@
if(!I->isReal && !I->isNew) return false; // should be deleted
assert(false); // unreachable
+ return false;
}
bool Collider::mayCollide(const Body* b1, const Body* b2){
Modified: trunk/lib/SConscript
===================================================================
--- trunk/lib/SConscript 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/lib/SConscript 2009-03-17 15:14:49 UTC (rev 1721)
@@ -113,8 +113,9 @@
['opengl/FpsTracker.cpp',
'opengl/GLTextLabel.cpp',
'opengl/GLWindow.cpp',
- 'opengl/GLWindowsManager.cpp'],
- LIBS=env['LIBS']+['glut']),
+ 'opengl/GLWindowsManager.cpp',
+ 'opengl/GLUtils.cpp'],
+ LIBS=env['LIBS']+['glut','$QGLVIEWER_LIB']),
env.SharedLibrary('XMLFormatManager',
Added: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -0,0 +1,135 @@
+#include "Dem3DofGeom_FacetSphere.hpp"
+#include<yade/pkg-common/InteractingSphere.hpp>
+#include<yade/pkg-common/InteractingFacet.hpp>
+YADE_PLUGIN("Dem3DofGeom_FacetSphere","GLDraw_Dem3DofGeom_FacetSphere","ef2_Facet_Sphere_Dem3DofGeom");
+
+CREATE_LOGGER(Dem3DofGeom_FacetSphere);
+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);
+ TRVAR3(cp1pt,cp2rel,contPtInTgPlane2()-contPtInTgPlane1());
+}
+
+void Dem3DofGeom_FacetSphere::relocateContactPoints(const Vector3r& p1, const Vector3r& p2){
+ //TRVAR2(p2.Length(),effR2);
+ if(p2.SquaredLength()>pow(effR2,2)){
+ setTgPlanePts(Vector3r::ZERO,p2-p1);
+ }
+}
+
+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());
+ Real sphereRadius=static_cast<InteractingSphere*>(cm2.get())->radius;
+ // begin facet-local coordinates
+ Vector3r contactLine=se31.orientation.Conjugate()*(se32.position-se31.position);
+ Vector3r normal=facet->nf;
+ Real L=normal.Dot(contactLine); // height/depth of sphere's center from facet's plane
+ if(L>sphereRadius && !c->isReal) return false; // sphere too far away from the plane
+
+ Vector3r contactPt=contactLine-L*normal; // projection of sphere's center to facet's plane (preliminary contact point)
+ const Vector3r* edgeNormals=facet->ne; // array[3] of edge normals (in facet plane)
+ int edgeMax=0; Real distMax=edgeNormals[0].Dot(contactPt);
+ for(int i=1; i<3; i++){
+ Real dist=edgeNormals[i].Dot(contactPt);
+ if(distMax<dist){edgeMax=i; distMax=dist;}
+ }
+ //TRVAR2(distMax,edgeMax);
+ // OK, what's the logic here? Copying from IF2IS4SCG…
+ Real sphereRReduced=shrinkFactor*sphereRadius;
+ Real inCircleR=facet->icr-sphereRReduced;
+ Real penetrationDepth;
+ if(inCircleR<0){inCircleR=facet->icr; sphereRReduced=0;}
+ if(distMax<inCircleR){// contact with facet's surface
+ penetrationDepth=sphereRadius-L;
+ normal.Normalize();
+ } else { // contact with the edge
+ contactPt+=edgeNormals[edgeMax]*(inCircleR-distMax);
+ //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);
+ // contact with vertex no. edgeMax+1
+ else if(contactPt.Dot(edgeNormals[edgeMax=(edgeMax+1)%3])>inCircleR) contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
+ normal=contactLine-contactPt;
+ //TRVAR4(contactLine,contactPt,normal,normal.Length());
+ //TRVAR3(se31.orientation*contactLine,se31.position+se31.orientation*contactPt,se31.orientation*normal);
+ penetrationDepth=sphereRadius-normal.Normalize();
+ //TRVAR1(penetrationDepth);
+ }
+ // end facet-local coordinates
+
+ if(penetrationDepth<0 && !c->isReal) return false;
+
+ shared_ptr<Dem3DofGeom_FacetSphere> fs;
+ Vector3r normalGlob=se31.orientation*normal;
+ if(c->interactionGeometry) fs=YADE_PTR_CAST<Dem3DofGeom_FacetSphere>(c->interactionGeometry);
+ else {
+ fs=shared_ptr<Dem3DofGeom_FacetSphere>(new Dem3DofGeom_FacetSphere());
+ c->interactionGeometry=fs;
+ fs->effR2=sphereRadius-penetrationDepth;
+ fs->refLength=fs->effR2;
+ fs->cp1pt=contactPt; // facet-local intial contact point
+ fs->localFacetNormal=normal;
+ fs->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normalGlob)); // initial sphere-local center-contactPt orientation WRT +x
+ fs->cp2rel.Normalize();
+ }
+ fs->se31=se31; fs->se32=se32;
+ fs->normal=normalGlob;
+ fs->contactPoint=se32.position+(-normalGlob)*(sphereRadius-penetrationDepth);
+ if(c->isNew){
+ TRVAR1(penetrationDepth);
+ TRVAR3(fs->refLength,fs->cp1pt,fs->localFacetNormal);
+ TRVAR3(fs->effR2,fs->cp2rel,fs->normal);
+ TRVAR2(fs->se31.orientation,fs->se32.orientation);
+ TRVAR2(fs->contPtInTgPlane1(),fs->contPtInTgPlane2());
+ }
+ return true;
+}
+
+#include<yade/lib-opengl/OpenGLWrapper.hpp>
+#include<yade/lib-opengl/GLUtils.hpp>
+
+bool GLDraw_Dem3DofGeom_FacetSphere::normal=false;
+bool GLDraw_Dem3DofGeom_FacetSphere::rolledPoints=false;
+bool GLDraw_Dem3DofGeom_FacetSphere::unrolledPoints=false;
+bool GLDraw_Dem3DofGeom_FacetSphere::shear=false;
+bool GLDraw_Dem3DofGeom_FacetSphere::shearLabel=false;
+
+void GLDraw_Dem3DofGeom_FacetSphere::go(const shared_ptr<InteractionGeometry>& ig, const shared_ptr<Interaction>& ip, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
+ Dem3DofGeom_FacetSphere* fs = static_cast<Dem3DofGeom_FacetSphere*>(ig.get());
+ //const Se3r& se31=b1->physicalParameters->dispSe3,se32=b2->physicalParameters->dispSe3;
+ const Se3r& se31=b1->physicalParameters->se3,se32=b2->physicalParameters->se3;
+ const Vector3r& pos1=se31.position; const Vector3r& pos2=se32.position;
+ const Quaternionr& ori1=se31.orientation; const Quaternionr& ori2=se32.orientation;
+ const Vector3r& contPt=fs->contactPoint;
+
+ if(normal){
+ GLUtils::GLDrawArrow(contPt,contPt+fs->refLength*fs->normal); // normal of the contact
+ }
+ // sphere center to point on the sphere
+ if(rolledPoints){
+ //cerr<<pos1<<" "<<pos1+ori1*fs->cp1pt<<" "<<contPt<<endl;
+ GLUtils::GLDrawLine(pos1+ori1*fs->cp1pt,contPt,Vector3r(0,.5,1));
+ GLUtils::GLDrawLine(pos2,pos2+(ori2*fs->cp2rel*Vector3r::UNIT_X*fs->effR2),Vector3r(0,1,.5));
+ }
+ // contact point to projected points
+ if(unrolledPoints||shear){
+ Vector3r ptTg1=fs->contPtInTgPlane1(), ptTg2=fs->contPtInTgPlane2();
+ if(unrolledPoints){
+ //TRVAR3(ptTg1,ptTg2,ss->normal)
+ GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1));
+ GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5)); GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
+ }
+ if(shear){
+ GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
+ if(shearLabel) GLUtils::GLDrawNum(fs->displacementT().Length(),contPt,Vector3r(1,1,1));
+ }
+ }
+}
+
+
+
Added: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.hpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -0,0 +1,72 @@
+#pragma once
+#include<yade/pkg-dem/DemXDofGeom.hpp>
+// for static roll/unroll functions in Dem3DofGeom_SphereSphere
+#include<yade/pkg-dem/Dem3DofGeom_SphereSphere.hpp>
+
+class Dem3DofGeom_FacetSphere: public Dem3DofGeom{
+ //! turn planePt in plane with normal0 around line passing through origin to plane with normal1
+ static Vector3r turnPlanePt(const Vector3r& planePt, const Vector3r& normal0, const Vector3r& normal1){ Quaternionr q; q.Align(normal0,normal1); return q*planePt; }
+
+ Vector3r contPtInTgPlane1() const { return turnPlanePt(se31.position+se31.orientation*cp1pt-contactPoint,se31.orientation*localFacetNormal,normal); }
+ Vector3r contPtInTgPlane2() const { return Dem3DofGeom_SphereSphere::unrollSpherePtToPlane(se32.orientation*cp2rel,effR2,-normal);}
+
+ public:
+ Dem3DofGeom_FacetSphere(){createIndex();}
+ virtual ~Dem3DofGeom_FacetSphere();
+ /******* 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; }
+ /***** end API ******/
+
+ void setTgPlanePts(const Vector3r&, const Vector3r&);
+ void relocateContactPoints(){ relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2()); }
+ void relocateContactPoints(const Vector3r& p1, const Vector3r& p2);
+ //! reference contact point on the facet in facet-local coords
+ 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_CLASS_AND_BASE(Dem3DofGeom_FacetSphere,Dem3DofGeom);
+ DECLARE_LOGGER;
+ friend class GLDraw_Dem3DofGeom_FacetSphere;
+ friend class ef2_Facet_Sphere_Dem3DofGeom;
+};
+REGISTER_SERIALIZABLE(Dem3DofGeom_FacetSphere);
+
+#include<yade/pkg-common/GLDrawFunctors.hpp>
+class GLDraw_Dem3DofGeom_FacetSphere:public GLDrawInteractionGeometryFunctor{
+ public:
+ virtual void go(const shared_ptr<InteractionGeometry>&,const shared_ptr<Interaction>&,const shared_ptr<Body>&,const shared_ptr<Body>&,bool wireFrame);
+ static bool normal,rolledPoints,unrolledPoints,shear,shearLabel;
+ RENDERS(Dem3DofGeom_FacetSphere);
+ REGISTER_CLASS_AND_BASE(GLDraw_Dem3DofGeom_FacetSphere,GLDrawInteractionGeometryFunctor);
+ REGISTER_ATTRIBUTES(GLDrawInteractionGeometryFunctor, (normal)(rolledPoints)(unrolledPoints)(shear)(shearLabel) );
+};
+REGISTER_SERIALIZABLE(GLDraw_Dem3DofGeom_FacetSphere);
+
+#include<yade/pkg-common/InteractionGeometryEngineUnit.hpp>
+class ef2_Facet_Sphere_Dem3DofGeom:public InteractionGeometryEngineUnit{
+ public:
+ virtual bool go(const shared_ptr<InteractingGeometry>& cm1, const shared_ptr<InteractingGeometry>& cm2, const Se3r& se31, const Se3r& se32, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<InteractingGeometry>& cm1, const shared_ptr<InteractingGeometry>& cm2, const Se3r& se31, const Se3r& se32, const shared_ptr<Interaction>& c){
+ c->swapOrder(); return go(cm2,cm1,se32,se31,c);
+ LOG_ERROR("!! goReverse maybe doesn't work in ef2_Facet_Sphere_Dem3DofGeom. InteractionGeometryMetaEngine should swap interaction members first and call go(...) afterwards.");
+ }
+ //! Reduce the facet's size, probably to avoid singularities at common facets' edges (?)
+ Real shrinkFactor;
+ ef2_Facet_Sphere_Dem3DofGeom(): shrinkFactor(0.) {}
+ FUNCTOR2D(InteractingFacet,InteractingSphere);
+ DEFINE_FUNCTOR_ORDER_2D(InteractingFacet,InteractingSphere);
+ REGISTER_CLASS_AND_BASE(ef2_Facet_Sphere_Dem3DofGeom,InteractionGeometryEngineUnit);
+ REGISTER_ATTRIBUTES(InteractionGeometryEngineUnit,(shrinkFactor));
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(ef2_Facet_Sphere_Dem3DofGeom);
+
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -1,5 +1,9 @@
#include "Dem3DofGeom_SphereSphere.hpp"
+#include<yade/pkg-common/InteractingSphere.hpp>
+YADE_PLUGIN("Dem3DofGeom_SphereSphere","GLDraw_Dem3DofGeom_SphereSphere","ef2_Sphere_Sphere_Dem3DofGeom");
+
+
Dem3DofGeom_SphereSphere::~Dem3DofGeom_SphereSphere(){}
/*! Project point from sphere surface to tangent plane,
@@ -25,17 +29,17 @@
*
* @param planePt point on the tangent plane, with origin at the contact point (i.e. at sphere center + normal*radius)
* @param radius sphere radius
- * @param planeNorma _unit_ vector pointing away from sphere center
+ * @param planeNormal _unit_ vector pointing away from sphere center
* @returns orientation that transforms +x axis to the vector between sphere center and point on the sphere that corresponds to planePt.
*
* @note It is not checked whether planePt relly lies on the tangent plane. If not, result will be incorrect.
*/
Quaternionr Dem3DofGeom_SphereSphere::rollPlanePtToSphere(const Vector3r& planePt, const Real& radius, const Vector3r& planeNormal){
- Vector3r axis=planeNormal.Cross(planePt); axis.Normalize();
- Real angle=planePt.Length()/radius;
- Quaternionr normal2pt(axis,angle);
- Quaternionr ret; ret.Align(Vector3r::UNIT_X,normal2pt*planeNormal);
- return ret;
+ Vector3r axis=planeNormal.Cross(planePt); axis.Normalize();
+ Real angle=planePt.Length()/radius;
+ Quaternionr normal2pt(axis,angle);
+ Quaternionr ret; ret.Align(Vector3r::UNIT_X,normal2pt*planeNormal);
+ return ret;
}
@@ -84,3 +88,90 @@
}
+#include<yade/lib-opengl/OpenGLWrapper.hpp>
+#include<yade/lib-opengl/GLUtils.hpp>
+
+
+
+bool GLDraw_Dem3DofGeom_SphereSphere::normal=false;
+bool GLDraw_Dem3DofGeom_SphereSphere::rolledPoints=false;
+bool GLDraw_Dem3DofGeom_SphereSphere::unrolledPoints=false;
+bool GLDraw_Dem3DofGeom_SphereSphere::shear=false;
+bool GLDraw_Dem3DofGeom_SphereSphere::shearLabel=false;
+
+void GLDraw_Dem3DofGeom_SphereSphere::go(const shared_ptr<InteractionGeometry>& ig, const shared_ptr<Interaction>& ip, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
+ Dem3DofGeom_SphereSphere* ss = static_cast<Dem3DofGeom_SphereSphere*>(ig.get());
+ //const Se3r& se31=b1->physicalParameters->dispSe3,se32=b2->physicalParameters->dispSe3;
+ const Se3r& se31=b1->physicalParameters->se3,se32=b2->physicalParameters->se3;
+ const Vector3r& pos1=se31.position,pos2=se32.position;
+ Vector3r& contPt=ss->contactPoint;
+
+ if(normal){
+ GLUtils::GLDrawArrow(contPt,contPt+ss->normal*.5*ss->refLength); // normal of the contact
+ }
+ #if 0
+ // never used, since bending/torsion not used
+ //Vector3r contPt=se31.position+(ss->effR1/ss->refLength)*(se32.position-se31.position); // must be recalculated to not be unscaled if scaling displacements ...
+ GLUtils::GLDrawLine(pos1,pos2,Vector3r(.5,.5,.5));
+ Vector3r bend; Real tors;
+ ss->bendingTorsionRel(bend,tors);
+ GLUtils::GLDrawLine(contPt,contPt+10*ss->radius1*(bend+ss->normal*tors),Vector3r(1,0,0));
+ #if 0
+ GLUtils::GLDrawNum(bend[0],contPt-.2*ss->normal*ss->radius1,Vector3r(1,0,0));
+ GLUtils::GLDrawNum(bend[1],contPt,Vector3r(0,1,0));
+ GLUtils::GLDrawNum(bend[2],contPt+.2*ss->normal*ss->radius1,Vector3r(0,0,1));
+ GLUtils::GLDrawNum(tors,contPt+.5*ss->normal*ss->radius2,Vector3r(1,1,0));
+ #endif
+ #endif
+ // sphere center to point on the sphere
+ if(rolledPoints){
+ GLUtils::GLDrawLine(pos1,pos1+(ss->ori1*ss->cp1rel*Vector3r::UNIT_X*ss->effR1),Vector3r(0,.5,1));
+ GLUtils::GLDrawLine(pos2,pos2+(ss->ori2*ss->cp2rel*Vector3r::UNIT_X*ss->effR2),Vector3r(0,1,.5));
+ }
+ //TRVAR4(pos1,ss->ori1,pos2,ss->ori2);
+ //TRVAR2(ss->cp2rel,pos2+(ss->ori2*ss->cp2rel*Vector3r::UNIT_X*ss->effR2));
+ // contact point to projected points
+ if(unrolledPoints||shear){
+ Vector3r ptTg1=ss->contPtInTgPlane1(), ptTg2=ss->contPtInTgPlane2();
+ if(unrolledPoints){
+ //TRVAR3(ptTg1,ptTg2,ss->normal)
+ GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1)); GLUtils::GLDrawLine(pos1,contPt+ptTg1,Vector3r(0,.5,1));
+ GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5)); GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
+ }
+ if(shear){
+ GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
+ if(shearLabel) GLUtils::GLDrawNum(ss->displacementT().Length(),contPt,Vector3r(1,1,1));
+ }
+ }
+}
+
+CREATE_LOGGER(ef2_Sphere_Sphere_Dem3DofGeom);
+
+bool ef2_Sphere_Sphere_Dem3DofGeom::go(const shared_ptr<InteractingGeometry>& cm1, const shared_ptr<InteractingGeometry>& cm2, const Se3r& se31, const Se3r& se32, const shared_ptr<Interaction>& c){
+ InteractingSphere *s1=static_cast<InteractingSphere*>(cm1.get()), *s2=static_cast<InteractingSphere*>(cm2.get());
+ Vector3r normal=se32.position-se31.position;
+ Real penetrationDepthSq=pow(distanceFactor*(s1->radius+s2->radius),2)-normal.SquaredLength();
+ if (penetrationDepthSq<0 && !c->isReal) return false;
+
+ Real dist=normal.Normalize(); /* Normalize() works in-place and returns length before normalization; from here, normal is unit vector */
+ shared_ptr<Dem3DofGeom_SphereSphere> ss;
+ if(c->interactionGeometry) ss=YADE_PTR_CAST<Dem3DofGeom_SphereSphere>(c->interactionGeometry);
+ else {
+ ss=shared_ptr<Dem3DofGeom_SphereSphere>(new Dem3DofGeom_SphereSphere());
+ c->interactionGeometry=ss;
+ // constants
+ ss->refLength=dist;
+ 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;
+ // quasi-constants
+ ss->cp1rel.Align(Vector3r::UNIT_X,se31.orientation.Conjugate()*normal);
+ ss->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normal));
+ ss->cp1rel.Normalize(); ss->cp2rel.Normalize();
+ }
+ ss->normal=normal;
+ ss->contactPoint=se31.position+(ss->effR1-.5*(ss->refLength-dist))*ss->normal;
+ ss->se31=se31; ss->se32=se32;
+ return true;
+}
+
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -15,15 +15,17 @@
public:
//! Positions and orientations of both spheres; must be updated at every iteration by the geom functor
Se3r se31, se32;
- const Vector3r& pos1, pos2; const Quaternionr& ori1, ori2;
+ //! relative orientation of the contact point with regards to sphere-local +x axis (quasi-constant)
Quaternionr cp1rel, cp2rel;
- Dem3DofGeom_SphereSphere(): pos1(se31.position), pos2(se32.position), ori1(se31.orientation), ori2(se32.orientation){}
+ //! shorthands
+ const Vector3r &pos1; const Quaternionr &ori1; const Vector3r &pos2; const Quaternionr &ori2;
+ Dem3DofGeom_SphereSphere(): pos1(se31.position), ori1(se31.orientation), pos2(se32.position), ori2(se32.orientation){createIndex();}
~Dem3DofGeom_SphereSphere();
//! effective radii of spheres for this interaction; can be smaller/larger than actual radii, but quasi-constant throughout the interaction life
Real effR1, effR2;
/********* API **********/
- Real displacementN(){ return (pos2-pos2).Length()-refLength; }
+ Real displacementN(){ return (pos2-pos1).Length()-refLength; }
Vector3r displacementT() {
// enabling automatic relocation decreases overall simulation speed by about 3%; perhaps: bool largeStrains ... ?
#if 1
@@ -35,7 +37,37 @@
Real slipToDisplacementTMax(Real displacementTMax);
/********* end API ***********/
- REGISTER_ATTRIBUTES(Dem3DofGeom,(se31)(se32)(effR1)(effR2));
+ REGISTER_ATTRIBUTES(Dem3DofGeom,(se31)(se32)(effR1)(effR2)(cp1rel)(cp2rel));
+ REGISTER_CLASS_AND_BASE(Dem3DofGeom_SphereSphere,Dem3DofGeom);
+ friend class GLDraw_Dem3DofGeom_SphereSphere;
+ friend class ef2_Sphere_Sphere_Dem3DofGeom;
};
REGISTER_SERIALIZABLE(Dem3DofGeom_SphereSphere);
+#include<yade/pkg-common/GLDrawFunctors.hpp>
+class GLDraw_Dem3DofGeom_SphereSphere:public GLDrawInteractionGeometryFunctor{
+ public:
+ virtual void go(const shared_ptr<InteractionGeometry>&,const shared_ptr<Interaction>&,const shared_ptr<Body>&,const shared_ptr<Body>&,bool wireFrame);
+ static bool normal,rolledPoints,unrolledPoints,shear,shearLabel;
+ //RENDERS(Dem3DofGeom_SphereSphere);
+ //REGISTER_CLASS_AND_BASE(GLDraw_Dem3DofGeom_SphereSphere,GLDrawInteractionGeometryFunctor);
+ REGISTER_ATTRIBUTES(GLDrawInteractionGeometryFunctor,(normal)(rolledPoints)(unrolledPoints)(shear)(shearLabel));
+};
+REGISTER_SERIALIZABLE(GLDraw_Dem3DofGeom_SphereSphere);
+
+#include<yade/pkg-common/InteractionGeometryEngineUnit.hpp>
+class ef2_Sphere_Sphere_Dem3DofGeom:public InteractionGeometryEngineUnit{
+ public:
+ virtual bool go(const shared_ptr<InteractingGeometry>& cm1, const shared_ptr<InteractingGeometry>& cm2, const Se3r& se31, const Se3r& se32, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<InteractingGeometry>&, const shared_ptr<InteractingGeometry>&, const Se3r&, const Se3r&, const shared_ptr<Interaction>&){throw runtime_error("goReverse on symmetric functor should never be called!");}
+ //! Factor of sphere radius such that sphere "touch" if their centers are not further than distanceFactor*(r1+r2); defaults to 1.
+ Real distanceFactor;
+ ef2_Sphere_Sphere_Dem3DofGeom(): distanceFactor(1.) {}
+ FUNCTOR2D(InteractingSphere,InteractingSphere);
+ DEFINE_FUNCTOR_ORDER_2D(InteractingSphere,InteractingSphere);
+ REGISTER_CLASS_AND_BASE(ef2_Sphere_Sphere_Dem3DofGeom,InteractionGeometryEngineUnit);
+ REGISTER_ATTRIBUTES(InteractionGeometryEngineUnit,(distanceFactor));
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(ef2_Sphere_Sphere_Dem3DofGeom);
+
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -13,21 +13,27 @@
Vector3r normal;
//! some reference point for the interaction
Vector3r contactPoint;
+ //! make strain go to -∞ for length going to zero
+ bool logCompression;
// API that needs to be implemented in derived classes
- virtual Real displacementN()=0;
- virtual Vector3r displacementT()=0;
- virtual Real slipToDisplacementTMax(Real displacementTMax)=0; // plasticity
+ virtual Real displacementN(){throw;}
+ virtual Vector3r displacementT(){throw;}
+ virtual Real slipToDisplacementTMax(Real displacementTMax){throw;}; // plasticity
// end API
- Real strainN(){return displacementN()/refLength;}
+ Real strainN(){
+ //if(!logCompression)
+ return displacementN()/refLength;
+ //else{Real dn=displacementN(); }
+ }
Vector3r strainT(){return displacementT()/refLength;}
Real slipToStrainTMax(Real strainTMax){return slipToDisplacementTMax(strainTMax*refLength)/refLength;}
REGISTER_CLASS_AND_BASE(Dem3DofGeom,InteractionGeometry);
REGISTER_ATTRIBUTES(InteractionGeometry,(refLength)(normal)(contactPoint));
};
-//REGISTER_SERIALIZABLE(Dem3DofGeom);
+REGISTER_SERIALIZABLE(Dem3DofGeom);
#if 0
/*! Abstract class for providing torsion and bending, in addition to inherited normal and shear strains. */
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -36,7 +36,7 @@
{
InteractingFacet* facet = static_cast<InteractingFacet*>(cm1.get());
/* could be written as (needs to be tested):
- * Vector3r cl=se32.orientation.Conjugate()*(se32.position-se32.position);
+ * Vector3r cl=se31.orientation.Conjugate()*(se32.position-se31.position);
*/
Matrix3r facetAxisT; se31.orientation.ToRotationMatrix(facetAxisT);
Matrix3r facetAxis = facetAxisT.Transpose();
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -33,7 +33,7 @@
FUNCTOR2D(InteractingSphere,InteractingSphere);
- //FIXME: what is this good for?!
+ // needed for the dispatcher, even if it is symmetric
DEFINE_FUNCTOR_ORDER_2D(InteractingSphere,InteractingSphere);
protected :
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/dem/SConscript 2009-03-17 15:14:49 UTC (rev 1721)
@@ -29,6 +29,10 @@
env.Install('$PREFIX/lib/yade$SUFFIX/pkg-dem',[
+ 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('SQLiteRecorder',
['Engine/StandAloneEngine/SQLiteRecorder.cpp'],
LIBS=env['LIBS']+['sqlite3x']),
Modified: trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp
===================================================================
--- trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp 2009-03-17 15:14:49 UTC (rev 1721)
@@ -231,7 +231,8 @@
//return true;
#else
- std::cerr << "ERROR: full version of wm3 library is needed\n";
+ LOG_FATAL("Using miniWm3; recompile with full Wm3 support to make snow folly functional.");
+ throw runtime_error("full wm3 required (message above).");
#endif
}
// if(! m1->depths[id2].empty()) m1->depths[id2].clear();
@@ -435,7 +436,8 @@
//return true;
#else
- std::cerr << "ERROR: full version of wm3 library is needed\n";
+ LOG_FATAL("Using miniWm3; recompile with full Wm3 support to make snow folly functional.");
+ throw runtime_error("full wm3 required (message above).");
#endif
}
if(! m1->depths[id2].empty()) m1->depths[id2].clear();
Added: trunk/scripts/test/Dem3DofGeom.py
===================================================================
--- trunk/scripts/test/Dem3DofGeom.py 2009-03-16 00:20:35 UTC (rev 1720)
+++ trunk/scripts/test/Dem3DofGeom.py 2009-03-17 15:14:49 UTC (rev 1721)
@@ -0,0 +1,39 @@
+
+O.bodies.append([
+ #utils.sphere([0,0,0],1,dynamic=False,color=(0,1,0),wire=True),
+ utils.facet(([2,2,1],[-2,0,1],[2,-2,1]),dynamic=False,color=(0,1,0),wire=False),
+ utils.sphere([-1,0,2],1,dynamic=True,color=(1,0,0),wire=True),
+])
+O.engines=[
+ MetaEngine('BoundingVolumeMetaEngine',[
+ EngineUnit('InteractingSphere2AABB'),
+ EngineUnit('InteractingFacet2AABB'),
+ ]),
+ StandAloneEngine('PersistentSAPCollider'),
+ MetaEngine('InteractionGeometryMetaEngine',[
+ EngineUnit('ef2_Sphere_Sphere_Dem3DofGeom'),
+ EngineUnit('ef2_Facet_Sphere_Dem3DofGeom')
+ ]),
+ #DeusExMachina('GravityEngine',{'gravity':[0,0,-10]}),
+ DeusExMachina('RotationEngine',{'rotationAxis':[0,1,0],'angularVelocity':10,'subscribedBodies':[1]}),
+ DeusExMachina('TranslationEngine',{'translationAxis':[1,0,0],'velocity':10,'subscribedBodies':[1]}),
+ #DeusExMachina('NewtonsDampedLaw')
+]
+O.miscParams=[
+ Generic('GLDraw_Dem3DofGeom_SphereSphere',{'normal':True,'rolledPoints':True,'unrolledPoints':True,'shear':True,'shearLabel':True}),
+ Generic('GLDraw_Dem3DofGeom_FacetSphere',{'normal':False,'rolledPoints':True,'unrolledPoints':True,'shear':True,'shearLabel':True}),
+ Generic('GLDrawSphere',{'glutUse':True})
+]
+
+try:
+ from yade import qt
+ renderer=qt.Renderer()
+ renderer['Body_wire']=True
+ renderer['Interaction_geometry']=True
+ qt.Controller()
+ qt.View()
+except ImportError: pass
+
+O.dt=1e-6
+O.saveTmp('init')
+O.run(50)
