yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

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 <eudoxos@xxxxxxxx>
-
-/*! 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 <eudoxos@xxxxxxxx>
+#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
+