yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2008-09-18 11:50:19 +0200 (Thu, 18 Sep 2008)
New Revision: 1517

Modified:
   trunk/extra/Brefcom.cpp
   trunk/extra/clump/Shop.cpp
   trunk/extra/clump/Shop.hpp
   trunk/lib/opengl/GLUtils.hpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp
   trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
Log:
1. Remove GL things from Shop, moved to lib/opengl/GLUtils.hpp
2. Cleanup of SpheresContactGeometry code, rename exactRot to hasShear



Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/extra/Brefcom.cpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -4,6 +4,7 @@
 #include<yade/pkg-dem/BodyMacroParameters.hpp>
 #include<yade/pkg-common/Sphere.hpp>
 #include<yade/lib-QGLViewer/qglviewer.h>
+#include<yade/lib-opengl/GLUtils.hpp>
 
 
 YADE_PLUGIN("BrefcomMakeContact","BrefcomContact","BrefcomLaw","GLDrawBrefcomContact","BrefcomDamageColorizer", "BrefcomPhysParams", "BrefcomGlobalCharacteristics" /* ,"BrefcomStiffnessComputer"*/ );
@@ -314,9 +315,9 @@
 		min(1.,max(0.,BC->epsTrans/BC->epsCrackOnset)),
 		min(1.,max(0.,abs(BC->epsTrans)/BC->epsCrackOnset-1)));
 
-	if(contactLine) Shop::GLDrawLine(b1->physicalParameters->dispSe3.position,b2->physicalParameters->dispSe3.position,lineColor);
-	if(dmgLabel){ Shop::GLDrawNum(BC->omega,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
-	else if(epsNLabel){ Shop::GLDrawNum(BC->epsN,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
+	if(contactLine) GLUtils::GLDrawLine(b1->physicalParameters->dispSe3.position,b2->physicalParameters->dispSe3.position,lineColor);
+	if(dmgLabel){ GLUtils::GLDrawNum(BC->omega,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
+	else if(epsNLabel){ GLUtils::GLDrawNum(BC->epsN,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
 
 	const Vector3r& cp=static_pointer_cast<SpheresContactGeometry>(i->interactionGeometry)->contactPoint;
 	if(epsT){
@@ -325,16 +326,16 @@
 		Real scale=.5*BC->equilibriumDist;
 		Vector3r dirShear=BC->epsT; dirShear.Normalize();
 		if(epsTAxes){
-			Shop::GLDrawLine(cp-Vector3r(scale,0,0),cp+Vector3r(scale,0,0));
-			Shop::GLDrawLine(cp-Vector3r(0,scale,0),cp+Vector3r(0,scale,0));
-			Shop::GLDrawLine(cp-Vector3r(0,0,scale),cp+Vector3r(0,0,scale));
+			GLUtils::GLDrawLine(cp-Vector3r(scale,0,0),cp+Vector3r(scale,0,0));
+			GLUtils::GLDrawLine(cp-Vector3r(0,scale,0),cp+Vector3r(0,scale,0));
+			GLUtils::GLDrawLine(cp-Vector3r(0,0,scale),cp+Vector3r(0,0,scale));
 		}
-		Shop::GLDrawArrow(cp,cp+dirShear*relShear*scale,Vector3r(1.,0.,0.));
-		Shop::GLDrawLine(cp+dirShear*relShear*scale,cp+dirShear*scale,Vector3r(.3,.3,.3));
+		GLUtils::GLDrawArrow(cp,cp+dirShear*relShear*scale,Vector3r(1.,0.,0.));
+		GLUtils::GLDrawLine(cp+dirShear*relShear*scale,cp+dirShear*scale,Vector3r(.3,.3,.3));
 
-		/* normal strain */ Shop::GLDrawArrow(cp,cp+BC->prevNormal*(BC->epsN/maxShear),Vector3r(0.,1.,0.));
+		/* normal strain */ GLUtils::GLDrawArrow(cp,cp+BC->prevNormal*(BC->epsN/maxShear),Vector3r(0.,1.,0.));
 	}
-	//if(normal) Shop::GLDrawArrow(cp,cp+BC->prevNormal*.5*BC->equilibriumDist,Vector3r(0.,1.,0.));
+	//if(normal) GLUtils::GLDrawArrow(cp,cp+BC->prevNormal*.5*BC->equilibriumDist,Vector3r(0.,1.,0.));
 }
 
 /********************** BrefcomDamageColorizer ****************************/

Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/extra/clump/Shop.cpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -1071,56 +1071,3 @@
 	{0.370512,0.415453,0.055970,0.010546},
 	{-1.,-1.,-1.,-1. } /* sentinel: non-positive radius */
 };
-#if 0
-Real Shop::ElasticWaveTimestepEstimate(shared_ptr<MetaBody> rootBody){
-	Real minDt=std::numeric_limits<Real>::infinity();
-	FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
-		shared_ptr<Sphere> sphere=dynamic_pointer_cast<Sphere>(b->geometricalModel);
-		shared_ptr<ElasticBodyParameters> elast=dynamic_pointer_cast<ElasticBodyParameters>(b->physicalParameters);
-		if(!sphere || !elast) continue;
-		Real density=elast->mass/((4./3.)*pow(Mathr::PI*sphere->radius,3));
-		Real thisDt=2*sphere->radius/sqrt(elast->young/density);
-		minDt=std::min(minDt,thisDt);
-	}
-	return minDt;
-}
-#endif
-
-void Shop::GLDrawLine(Vector3r from, Vector3r to, Vector3r color){
-	glEnable(GL_LIGHTING); glColor3v(color);
-	glBegin(GL_LINES); glVertex3v(from); glVertex3v(to); glEnd();
-}
-
-void Shop::GLDrawArrow(Vector3r from, Vector3r to, Vector3r color){
-	glEnable(GL_LIGHTING); glColor3v(color); qglviewer::Vec a(from[0],from[1],from[2]),b(to[0],to[1],to[2]); QGLViewer::drawArrow(a,b);	
-}
-
-void Shop::GLDrawNum(Real n, Vector3r pos, Vector3r color, unsigned precision){
-	ostringstream oss; oss<<setprecision(precision)<< /* "w="<< */ (double)n;
-	Shop::GLDrawText(oss.str(),pos,color);
-}
-
-void Shop::GLDrawText(std::string txt, Vector3r pos, Vector3r color){
-	glPushMatrix();
-	glTranslatev(pos);
-	glColor3(color[0],color[1],color[2]);
-	glRasterPos2i(0,0);
-	for(unsigned int i=0;i<txt.length();i++)
-		glutBitmapCharacter(GLUT_BITMAP_HELVETICA_12, txt[i]);
-	glPopMatrix();
-
-}
-Real Shop::PWaveTimeStep(shared_ptr<MetaBody> _rb){
-	shared_ptr<MetaBody> rb=_rb;
-	if(!rb)rb=Omega::instance().getRootBody();
-	Real dt=std::numeric_limits<Real>::infinity();
-	FOREACH(const shared_ptr<Body>& b, *rb->bodies){
-		if(!b->physicalParameters || !b->geometricalModel) continue;
-		shared_ptr<ElasticBodyParameters> ebp=dynamic_pointer_cast<ElasticBodyParameters>(b->physicalParameters);
-		shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->geometricalModel);
-		if(!ebp || !s) continue;
-		Real density=ebp->mass/((4/3.)*Mathr::PI*pow(s->radius,3));
-		dt=min(dt,s->radius/sqrt(ebp->young/density));
-	}
-	return dt;
-}

Modified: trunk/extra/clump/Shop.hpp
===================================================================
--- trunk/extra/clump/Shop.hpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/extra/clump/Shop.hpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -79,11 +79,6 @@
 		//! Save spheres in the current simulation into a text file
 		static void saveSpheresToFile(string fileName);
 
-		static void GLDrawLine(Vector3r from, Vector3r to, Vector3r color=Vector3r(1,1,1));
-		static void GLDrawArrow(Vector3r from, Vector3r to, Vector3r color=Vector3r(1,1,1));
-		static void GLDrawText(std::string text, Vector3r pos, Vector3r color=Vector3r(1,1,1));
-		static void GLDrawNum(Real n, Vector3r pos, Vector3r color=Vector3r(1,1,1), unsigned precision=3);
-
 		/*! Cache for class indices for physical actions (body external variables, Bex)
 		 *
 		 * It is necessary to populate the cache by calling initCache(); then supported

Modified: trunk/lib/opengl/GLUtils.hpp
===================================================================
--- trunk/lib/opengl/GLUtils.hpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/lib/opengl/GLUtils.hpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -12,20 +12,20 @@
 
 struct GLUtils{
 
-	static void GLDrawArrow(Vector3r from, Vector3r to, Vector3r color){
+	static void GLDrawArrow(const Vector3r& from, const Vector3r& to, const Vector3r& color=Vector3r(1,1,1)){
 		glEnable(GL_LIGHTING); glColor3v(color); qglviewer::Vec a(from[0],from[1],from[2]),b(to[0],to[1],to[2]); QGLViewer::drawArrow(a,b);	
 	}
-	static void GLDrawLine(Vector3r from, Vector3r to, Vector3r color){
+	static void GLDrawLine(const Vector3r& from, const Vector3r& to, const Vector3r& color=Vector3r(1,1,1)){
 		glEnable(GL_LIGHTING); glColor3v(color);
 		glBegin(GL_LINES); glVertex3v(from); glVertex3v(to); glEnd();
 	}
 
-	static void GLDrawNum(Real n, Vector3r pos, Vector3r color, unsigned precision){
+	static void GLDrawNum(const Real& n, const Vector3r& pos, const Vector3r& color=Vector3r(1,1,1), unsigned precision=3){
 		std::ostringstream oss; oss<<std::setprecision(precision)<< /* "w="<< */ (double)n;
 		GLUtils::GLDrawText(oss.str(),pos,color);
 	}
 
-	static void GLDrawText(std::string txt, Vector3r pos, Vector3r color){
+	static void GLDrawText(const std::string& txt, const Vector3r& pos, const Vector3r& color=Vector3r(1,1,1)){
 		glPushMatrix();
 		glTranslatev(pos);
 		glColor3(color[0],color[1],color[2]);

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -9,6 +9,7 @@
  * (should be on the plane passing through origin and oriented with normal; not checked!)
  */
 void SpheresContactGeometry::setTgPlanePts(Vector3r p1new, Vector3r p2new){
+	assert(hasShear);
 	cp1rel=ori1.Conjugate()*rollPlanePtToSphere(p1new,d1,normal);
 	cp2rel=ori2.Conjugate()*rollPlanePtToSphere(p2new,d2,-normal);
 }
@@ -19,6 +20,7 @@
  * flip axis and the point would project on other side of the tangent plane piece.
  */
 void SpheresContactGeometry::relocateContactPoints(){
+	assert(hasShear);
 	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
 	Vector3r midPt=.5*(p1+p2);
 	if((p1.SquaredLength()>4*d1 || p2.SquaredLength()>4*d2) && midPt.SquaredLength()>.5*min(d1,d2)){
@@ -29,15 +31,17 @@
 }
 
 /*! Perform slip of the projected contact points so that their distance becomes equal (or remains smaller) than the given one.
- * TODO: not yet tested
+ * The slipped distance is returned.
  */
-void SpheresContactGeometry::slipToDistIfNeeded(Real dist){
+Real SpheresContactGeometry::slipToDisplacementTMax(Real displacementTMax){
+	assert(hasShear);
+	assert(displacementTMax>Mathr::ZERO_TOLERANCE);
 	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
 	Real currDist=(p2-p1).Length();
-	if(currDist<dist) return; // close enough, no slip needed
-	assert(dist>Mathr::ZERO_TOLERANCE);
-	Vector3r diff=.5*(currDist/dist-1)*(p2-p1);
+	if(currDist<displacementTMax) return 0; // close enough, no slip needed
+	Vector3r diff=.5*(currDist/displacementTMax-1)*(p2-p1);
 	setTgPlanePts(p1+diff,p2-diff);
+	return 2*diff.Length();
 }
 
 /*! Project point from sphere surface to tangent plane,

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -30,7 +30,7 @@
 			contactPoint;
 		Real radius1,radius2,penetrationDepth;
 
-		bool exactRot; // whether the exact rotation code is being used -- following fields are needed for that
+		bool hasShear; // whether the exact rotation code is being used -- following fields are needed for that
 		//! positions and orientations of both spheres -- must be updated at every iteration
 		Vector3r pos1, pos2; Quaternionr ori1, ori2;
 		/*! Orientation of the contact point relative to each sphere-local coordinates.
@@ -61,18 +61,21 @@
 		Vector3r displacementT(bool relocate=true){ return contPtInTgPlane2()-contPtInTgPlane1();}
 		Vector3r epsT(){return displacementT()*(1./d0);}
 	
-		void slipToDistIfNeeded(Real slip);
+		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);}
+
 		void relocateContactPoints();
 
-		SpheresContactGeometry():InteractionGeometry(),contactPoint(Vector3r::ZERO),radius1(0),radius2(0),exactRot(false){createIndex();}
+		SpheresContactGeometry():InteractionGeometry(),contactPoint(Vector3r::ZERO),radius1(0),radius2(0),hasShear(false){createIndex();}
 		virtual ~SpheresContactGeometry(){};
 	protected :
 		virtual void registerAttributes(){
 			REGISTER_ATTRIBUTE(radius1);
 			REGISTER_ATTRIBUTE(radius2);
 			REGISTER_ATTRIBUTE(contactPoint); // to allow access from python
-			// exactRot
-			REGISTER_ATTRIBUTE(exactRot);
+			// hasShear
+			REGISTER_ATTRIBUTE(hasShear);
 			REGISTER_ATTRIBUTE(pos1);
 			REGISTER_ATTRIBUTE(pos2);
 			REGISTER_ATTRIBUTE(ori1);

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -24,7 +24,7 @@
 void InteractingSphere2InteractingSphere4SpheresContactGeometry::registerAttributes()
 {	
 	REGISTER_ATTRIBUTE(interactionDetectionFactor);
-	REGISTER_ATTRIBUTE(exactRot);
+	REGISTER_ATTRIBUTE(hasShear);
 }
 
 bool InteractingSphere2InteractingSphere4SpheresContactGeometry::go(	const shared_ptr<InteractingGeometry>& cm1,
@@ -51,8 +51,8 @@
 		scm->radius1 = s1->radius;
 		scm->radius2 = s2->radius;
 		if (!c->interactionGeometry) c->interactionGeometry = scm;
-		if(exactRot){
-			scm->exactRot=true;
+		if(hasShear){
+			scm->hasShear=true;
 			scm->pos1=se31.position; scm->pos2=se32.position;
 			scm->ori1=se31.orientation; scm->ori2=se32.orientation;
 			if(c->isNew){
@@ -64,8 +64,13 @@
 				scm->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normal));
 				scm->cp1rel.Normalize(); scm->cp2rel.Normalize();
 			}
-			// FIXME: temporary, testing only!
-			if((Omega::instance().getCurrentIteration()%500)==0) scm->relocateContactPoints();
+			// testing only
+			#if 0
+			if((Omega::instance().getCurrentIteration()%10000)==0) scm->relocateContactPoints();
+			if((Omega::instance().getCurrentIteration()%10000)==0) {
+				Real slip=scm->slipToEpsNMax(1.); if(slip>0) cerr<<"SLIP by Δε_N "<<slip<<" → ε_N="<<scm->epsN()<<endl;
+			}
+			#endif
 		}
 		return true;
 	} else return false;

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -26,8 +26,8 @@
 		 * @note This parameter is functionally coupled with InteractinSphere2AABB::aabbEnlargeFactor,
 		 * which will create larger bounding boxes and should be of the same value. */
 		double interactionDetectionFactor;
-		/*! Whether we create SpheresContactGeometry with data necessary for exact rotation computation */
-		bool exactRot;
+		/*! Whether we create SpheresContactGeometry with data necessary for (exact) shear computation */
+		bool hasShear;
 
 	REGISTER_CLASS_NAME(InteractingSphere2InteractingSphere4SpheresContactGeometry);
 	REGISTER_BASE_CLASS_NAME(InteractionGeometryEngineUnit);

Modified: trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp	2008-09-17 21:48:17 UTC (rev 1516)
+++ trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp	2008-09-18 09:50:19 UTC (rev 1517)
@@ -74,7 +74,7 @@
 		#endif
 	}
 
-	if(sc->exactRot){
+	if(sc->hasShear){
 		GLUtils::GLDrawLine(sc->pos1,sc->pos2,Vector3r(.5,.5,.5));
 		// sphere center to point on the sphere
 		//GLUtils::GLDrawText("[1]",sc->pos1,Vector3r(1,1,1)); GLUtils::GLDrawText("[2]",sc->pos2,Vector3r(1,1,1));