yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-03-23 20:35:09 +0100 (Mon, 23 Mar 2009)
New Revision: 1728

Modified:
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
   trunk/scripts/test/triax-identical-results.py
Log:
1. Remove some garbage from SpheresContactGeometry
2. Verify that SCG_SHEAR doesn't alter behavior if ElasticContactLaw::useShear is false
3. Implement SCG_SHEAR for sphere-box interactions
4. sphere-box interactions no longer call goReverse, but swap interaction order instead, as facets do.
5. Fix triax-idnetical-results.py to reload generated initial config to avoid rounding issues of sphere coords in text file.



Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2009-03-23 19:35:09 UTC (rev 1728)
@@ -10,12 +10,12 @@
 SpheresContactGeometry::~SpheresContactGeometry(){};
 
 #ifdef SCG_SHEAR
-void SpheresContactGeometry::updateShear(const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
+Vector3r SpheresContactGeometry::updateShear(const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
 
 	Vector3r axis;
 	Real angle;
 
-	shearIncrement=Vector3r::ZERO;
+	Vector3r shearIncrement(Vector3r::ZERO);
 
 	// approximated rotations
 		axis = prevNormal.Cross(normal); 
@@ -62,7 +62,7 @@
 	shearIncrement -= shearDisplacement;
 
 	shear+=shearIncrement;
-	shearUpdateIter=Omega::instance().getCurrentIteration();
+	return shearIncrement;
 }
 #endif
 

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-03-23 19:35:09 UTC (rev 1728)
@@ -24,20 +24,13 @@
 			Vector3r shear;
 			//! Normal of the contact in the previous step
 			Vector3r prevNormal;
-			//! Increment of shear displacement in last step (is already added to shear); perhaps useful for viscosity or something similar
-			Vector3r shearIncrement;
-			long shearUpdateIter; // debugging only, to check when shear was updated last time
-			//! update shear on this contact given dynamic parameters of bodies. Should be called from constitutive law, exactly once per iteration
-			void updateShear(const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
-			//const Vector3r& getShear(){ if(Omega::instance().getCurrentIteration()>shearUpdateIter) throw runtime_error("SpheresContactGeometry::updateShear must be called prior to getShear()."); return shear; }
+			//! update shear on this contact given dynamic parameters of bodies. Should be called from constitutive law, exactly once per iteration. Returns shear increment in this update, which is already added to shear.
+			Vector3r updateShear(const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
 		#endif
 		
 		// all the rest here will hopefully disappear at some point...
 
 		// begin abusive storage
-		bool hasNormalViscosity;
-		Real NormalViscisity;
-		Real NormalRelativeVelocity;
 		//! Whether the original contact was on the positive (1) or negative (-1) facet side; this is to permit repulsion to the right side even if the sphere passes, under extreme pressure, geometrically to the other facet's side. This is used only in InteractingFacet2IteractingSphere4SpheresContactGeometry. If at any point the contact is with the edge or vertex instead of face, this attribute is reset so that contact with the other face is possible.
 		int facetContactFace;
 		// end abusive storage
@@ -99,7 +92,7 @@
 
 		SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();
 		#ifdef SCG_SHEAR
-			shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO /*initialized to aproper value by geom functor*/; shearIncrement=Vector3r::ZERO; shearUpdateIter=-1 /* proper value from geom functor */;
+			shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO /*initialized to proper value by geom functor*/;
 		#endif	
 		}
 		virtual ~SpheresContactGeometry();
@@ -114,8 +107,6 @@
 			#ifdef SCG_SHEAR
 				(shear)
 				(prevNormal)
-				(shearIncrement)
-				(shearUpdateIter)
 			#endif
 			(facetContactFace)
 			// hasShear

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp	2009-03-23 19:35:09 UTC (rev 1728)
@@ -87,6 +87,11 @@
 		shared_ptr<SpheresContactGeometry> scm;
 		if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
 		else scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
+
+		#ifdef SCG_SHEAR
+			if(c->isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); }
+			else {scm->prevNormal=scm->normal;}
+		#endif
 			
 		// contact point is in the middle of overlapping volumes
 		//(in the direction of penetration, which is normal to the box surface closest to sphere center) of overlapping volumes
@@ -129,6 +134,10 @@
 		shared_ptr<SpheresContactGeometry> scm;
 		if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
 		else scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);	
+		#ifdef SCG_SHEAR
+			if(c->isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
+			else {scm->prevNormal=scm->normal;}
+		#endif
 		scm->contactPoint = 0.5*(pt1+pt2);
 		//scm->normal = pt1-pt2; scm->normal.Normalize();
 		//scm->penetrationDepth = (pt1-pt2).Length();
@@ -150,19 +159,9 @@
 						const Se3r& se32,
 						const shared_ptr<Interaction>& c)
 {
-	bool isInteracting = go(cm2,cm1,se32,se31,c);
-	if (isInteracting)
-	{
-		SpheresContactGeometry* scm = static_cast<SpheresContactGeometry*>(c->interactionGeometry.get());
-		//Vector3r tmp = scm->closestsPoints[0].first;		
-		//scm->closestsPoints[0].first = scm->closestsPoints[0].second;
-		//scm->closestsPoints[0].second = tmp;
-		scm->normal = -scm->normal;
-		Real tmpR  = scm->radius1;
-		scm->radius1 = scm->radius2;
-		scm->radius2 = tmpR;
-	}
-	return isInteracting;
+	assert(c->isNew);
+	c->swapOrder();
+	return go(cm2,cm1,se32,se31,c);
 }
 
 YADE_PLUGIN();

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp	2009-03-23 19:35:09 UTC (rev 1728)
@@ -165,23 +165,10 @@
 								const Se3r& se32,
 								const shared_ptr<Interaction>& c)
 {
+	assert(c->isNew);
 	c->swapOrder();
 	//LOG_WARN("Swapped interaction order for "<<c->getId2()<<"&"<<c->getId1());
 	return go(cm2,cm1,se32,se31,c);
-#if 0	
-	bool isInteracting = go(cm2,cm1,se32,se31,c);
-	if (isInteracting)
-	{
-	    SpheresContactGeometry* scm = static_cast<SpheresContactGeometry*>(c->interactionGeometry.get());
-		 scm->normal*=-1;
-		 std::swap(scm->radius1,scm->radius2);
-		 if(hasShear){
-			 swap(scm->pos1,scm->pos2); swap(scm->ori1,scm->ori2);
-			 if(c->isNew){ swap(scm->cp1rel,scm->cp2rel); swap(scm->d1,scm->d2); }
-		 }
-	}
-	return isInteracting;
-#endif
 }
 
 YADE_PLUGIN();

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp	2009-03-23 19:35:09 UTC (rev 1728)
@@ -38,14 +38,8 @@
 		else { scm=shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry()); c->interactionGeometry=scm; }
 
 		#ifdef SCG_SHEAR
-			if(c->isNew){
-				scm->prevNormal=normal;
-				scm->shearUpdateIter=Omega::instance().getCurrentIteration(); /* no shear at the very beginning; shear initialized to zero vector in SCG ctor */
-			} else {
-				scm->prevNormal=scm->normal;
-				// make sure updateShear was properly called at last iteration; debugging only
-				//assert(scm->shearUpdateIter==Omega::instance().getCurrentIteration()-1);
-			}
+			if(c->isNew) scm->prevNormal=normal; 
+			else scm->prevNormal=scm->normal;
 		#endif
 
 		Real penetrationDepth=s1->radius+s2->radius-normal.Normalize(); /* Normalize() works in-place and returns length before normalization; from here, normal is unit vector */

Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp	2009-03-23 19:35:09 UTC (rev 1728)
@@ -123,10 +123,10 @@
 					currentContactGeometry->updateShear(de1,de2,dt);
 					shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
 				} else {
+			#endif
 					currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+			#ifdef SCG_SHEAR
 				}
-			#else
-				currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
 			#endif
 			
 			// PFC3d SlipModel, is using friction angle. CoulombCriterion

Modified: trunk/scripts/test/triax-identical-results.py
===================================================================
--- trunk/scripts/test/triax-identical-results.py	2009-03-22 18:08:17 UTC (rev 1727)
+++ trunk/scripts/test/triax-identical-results.py	2009-03-23 19:35:09 UTC (rev 1728)
@@ -13,18 +13,18 @@
 	if not exists(outSph): break
 	i+=1
 inSph='%s-in.spheres'%sph
-if exists(inSph):
-	print "Using existing initial configuration",inSph
-Preprocessor('TriaxialTest',{'importFilename':(inSph if exists(inSph) else ''),'numberOfGrains':400}).load()
-if not exists(inSph):
-	print "Saving new initial configuration to",inSph
+if exists(inSph): print "Using existing initial configuration",inSph
+else:
+	Preprocessor('TriaxialTest').load()
+	print "Using new initial configuration in",inSph
 	utils.spheresToFile(inSph)
+Preprocessor('TriaxialTest',{'importFilename':inSph}).load()
 O.usesTimeStepper=False
 O.dt=utils.PWaveTimeStep()
 #
 # uncomment this line to enable shear computation in SpheresContactGeometry and then compare results with this line commented
 #
-#[e for e in O.engines if e.name=='ElasticContactLaw'][0]['useShear']=True
+[e for e in O.engines if e.name=='ElasticContactLaw'][0]['useShear']=True
 if 1:
 	O.run(2000,True)
 	utils.spheresToFile(outSph)