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) _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp