Author: eudoxos Date: 2009-03-13 10:08:31 +0100 (Fri, 13 Mar 2009) New Revision: 1719
Added: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp Removed: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofContactGeometry.hpp Log: 1. Initial checkout for the DemXDofGeom classes (will contain the hasShear code from SpheresContactGeometry and more eventually). Deleted: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofContactGeometry.hpp =================================================================== --- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofContactGeometry.hpp 2009-03-12 17:27:35 UTC (rev 1718) +++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofContactGeometry.hpp 2009-03-13 09:08:31 UTC (rev 1719) @@ -1,41 +0,0 @@ -// 2009 © Václav Šmilauer <[email protected]> - -/*! Abstract class for geometry providing normal and shear strains (or displacements) */ -class Dem3DofContactGeometry: public InteractionGeometry { - public: - //! reference length of the contact (usually initial length) - Real refLength; - //! unit vector in the interaction direction (normal to the contact plane), always from id1 towards id2 - Vector3r normal; - //! contact point - Vector3r contactPoint; - //! Absolute displacement in the normal direction - virtual Real displacementN()=0; - //! Strain in the normal direction - virtual Real epsN(){ return displacementN()/refLength; } - //! Absolute displacement in the shear direction (in global coordinates) - virtual Vector3r displacementT()=0; - //! Strain in the shear direction (in global coordinates) - virtual Vector3r epsT(){ return displacementT/refLength; } - REGISTER_CLASS_AND_BASE(Dem3DofContactGeometry,InteractionGeometry); - REGISTER_ATTRIBUTES(InteractionGeometry,(refLength)(normal)(contactPoint)); -}; -REGISTER_SERIALIZABLE(Dem3DofContactGeometry); - -#if 0 -/*! Abstract class for providing torsion and bending, in addition to inherited normal and shear strains. */ -class Dem6DofContactGeometry: public Dem3DofContactGeometry { - public: - //! Absolute value of rotation along the interaction normal - virtual Real rotationN()=0; - //! Relative rotation along the contact normal - virtual Real torsion(){ return rotationN()/refLength; } - //! Absolute value of bending (in global coordinates) - virtual Vector3r rotationT()=0; - //! Relative rotation perpendicular to contact normal (in global coordinates) - virtual Vector3r bending(){ return rotationT()/refLength; } - REGISTER_CLASS_AND_BASE(Dem6DofContactGeometry,Dem3DofContactGeometry); - REGISTER_ATTRIBUTES(Dem3DofContactGeometry,); -}; -REGISTER_SERIALIZABLE(Dem6DofContactGeometry); -#endif Added: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp =================================================================== --- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-03-12 17:27:35 UTC (rev 1718) +++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-03-13 09:08:31 UTC (rev 1719) @@ -0,0 +1,86 @@ +#include "Dem3DofGeom_SphereSphere.hpp" + +Dem3DofGeom_SphereSphere::~Dem3DofGeom_SphereSphere(){} + +/*! Project point from sphere surface to tangent plane, + * such that the angle of shortest arc from (1,0,0) pt on the sphere to the point itself is the same + * as the angle of segment of the same length on the tangent plane. + * + * This function is (or should be) inverse of SpheresContactGeometry::rollPlanePtToSphere. + * + * @param fromXtoPtOri gives orientation of the vector from sphere center to the sphere point from the global +x axis. + * @param radius the distance from sphere center to the contact plane + * @param planeNormal unit vector pointing away from the sphere center, determining plane orientation on which the projected point lies. + * @returns The projected point coordinates (with origin at the contact point). + */ +Vector3r Dem3DofGeom_SphereSphere::unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& planeNormal){ + Quaternionr normal2pt; normal2pt.Align(planeNormal,fromXtoPtOri*Vector3r::UNIT_X); + Vector3r axis; Real angle; normal2pt.ToAxisAngle(axis,angle); + return (angle*radius) /* length */ *(axis.Cross(planeNormal)) /* direction: both are unit vectors */; +} + +/*! Project point from tangent plane to the sphere. + * + * This function is (or should be) inverse of SpheresContactGeometry::unrollSpherePtToPlane. + * + * @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 + * @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; +} + + + +/* Set contact points on both spheres such that their projection is the one given + * (should be on the plane passing through origin and oriented with normal; not checked!) + */ +void Dem3DofGeom_SphereSphere::setTgPlanePts(Vector3r p1new, Vector3r p2new){ + cp1rel=ori1.Conjugate()*rollPlanePtToSphere(p1new,effR1,normal); + cp2rel=ori2.Conjugate()*rollPlanePtToSphere(p2new,effR2,-normal); +} + + + +/*! Perform slip of the projected contact points so that their distance becomes equal (or remains smaller) than the given one. + * The slipped distance is returned. + */ +Real Dem3DofGeom_SphereSphere::slipToDisplacementTMax(Real displacementTMax){ + // very close, 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); + return 2*diff.Length(); +} + + +/* Move contact point on both spheres in such way that their relative position (displacementT) is the same; + * this should be done regularly to ensure that the angle doesn't go over π, since then quaternion would + * flip axis and the point would project on other side of the tangent plane piece. */ +void Dem3DofGeom_SphereSphere::relocateContactPoints(){ + relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2()); +} + +/*! Like Dem3DofGeom_SphereSphere::relocateContactPoints(), but use already computed tangent plane points. */ +void Dem3DofGeom_SphereSphere::relocateContactPoints(const Vector3r& p1, const Vector3r& p2){ + Vector3r midPt=(effR1/(effR1+effR2))*(p1+p2); // proportionally to radii, so that angle would be the same + if((p1.SquaredLength()>pow(effR1,2) || p2.SquaredLength()>pow(effR2,2)) && midPt.SquaredLength()>pow(min(effR1,effR2),2)){ + //cerr<<"RELOCATION with displacementT="<<displacementT(); // should be the same before and after relocation + setTgPlanePts(p1-midPt,p2-midPt); + //cerr<<" → "<<displacementT()<<endl; + } +} + + Added: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp =================================================================== --- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2009-03-12 17:27:35 UTC (rev 1718) +++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2009-03-13 09:08:31 UTC (rev 1719) @@ -0,0 +1,41 @@ +#pragma once +#include<yade/pkg-dem/DemXDofGeom.hpp> + +class Dem3DofGeom_SphereSphere: public Dem3DofGeom{ + public: + // auxiliary functions; the quaternion magic is all in here + static Vector3r unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& normal); + static Quaternionr rollPlanePtToSphere(const Vector3r& planePt, const Real& radius, const Vector3r& normal); + private: + Vector3r contPtInTgPlane1() const { return unrollSpherePtToPlane(ori1*cp1rel,effR1,normal); } + Vector3r contPtInTgPlane2() const { return unrollSpherePtToPlane(ori2*cp2rel,effR2,-normal);} + void setTgPlanePts(Vector3r p1new, Vector3r p2new); + 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; + const Vector3r& pos1, pos2; const Quaternionr& ori1, ori2; + Quaternionr cp1rel, cp2rel; + Dem3DofGeom_SphereSphere(): pos1(se31.position), pos2(se32.position), ori1(se31.orientation), ori2(se32.orientation){} + ~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; } + Vector3r displacementT() { + // enabling automatic relocation decreases overall simulation speed by about 3%; perhaps: bool largeStrains ... ? + #if 1 + Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2(); relocateContactPoints(p1,p2); return p2-p1; // shear before relocation, but that is OK + #else + return contPtInTgPlane2()-contPtInTgPlane1(); + #endif + } + Real slipToDisplacementTMax(Real displacementTMax); + /********* end API ***********/ + + REGISTER_ATTRIBUTES(Dem3DofGeom,(se31)(se32)(effR1)(effR2)); +}; +REGISTER_SERIALIZABLE(Dem3DofGeom_SphereSphere); + Added: trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp =================================================================== --- trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2009-03-12 17:27:35 UTC (rev 1718) +++ trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2009-03-13 09:08:31 UTC (rev 1719) @@ -0,0 +1,44 @@ +// 2009 © Václav Šmilauer <[email protected]> +#pragma once +#include<yade/core/InteractionGeometry.hpp> + +/*! Abstract base class for representing contact geometry of 2 elements that has 3 degrees of freedom: + * normal (1 component) and shear (Vector3r, but in plane perpendicular to the normal) + */ +class Dem3DofGeom: public InteractionGeometry { + public: + //! some length used to convert displacements to strains + Real refLength; + //! some unit reference vector (in the direction of the interaction) + Vector3r normal; + //! some reference point for the interaction + Vector3r contactPoint; + + // API that needs to be implemented in derived classes + virtual Real displacementN()=0; + virtual Vector3r displacementT()=0; + virtual Real slipToDisplacementTMax(Real displacementTMax)=0; // plasticity + // end API + + Real strainN(){return displacementN()/refLength;} + 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); + +#if 0 +/*! Abstract class for providing torsion and bending, in addition to inherited normal and shear strains. */ +class Dem6DofGeom: public Dem3DofGeom { + public: + //! rotations perpendicular to the normal (bending; in global coords) and parallel with the normal (torsion) + void bendingTorsionAbs(Vector3r& bend, Real& tors)=0; + void bendingTorsionRel(Vector3r& bend, Real& tors){ bendingTorsionAbs(bend,tors); bend/=refLength; tors/=refLength;} + REGISTER_CLASS_AND_BASE(Dem6DofGeom,Dem3DofGeom); + REGISTER_ATTRIBUTES(Dem3DofGeom, /*nothing*/); +}; +REGISTER_SERIALIZABLE(Dem6DofGeom); +#endif + _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
