yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2008-09-24 10:42:04 +0200 (Wed, 24 Sep 2008)
New Revision: 1525

Modified:
   trunk/core/yade.cpp
   trunk/extra/Brefcom.cpp
   trunk/extra/Brefcom.hpp
   trunk/extra/clump/Shop.cpp
   trunk/extra/clump/Shop.hpp
   trunk/gui/py/PythonUI_rc.py
   trunk/gui/py/_utils.cpp
   trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
   trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp
   trunk/pkg/dem/DataClass/InteractionPhysics/ElasticContactInteraction.hpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp
Log:
1. Add NormalInteraction::normalForce and NormalShearInteraction::normalForce (moved up the hierarchy from ElasticContactInteraction) so that unbalancedForce can be calculated on any subclassed interaction type.
2. Adapt Brefcom for this.
3. Add Shop::unbalancedForce and utils.unbalancedForce
4. Add workaround for saving ipython history as suggested by http://lists.ipython.scipy.org/pipermail/ipython-user/2008-September/005839.html
5. Fix AxialGravityEngine to register parent class attributes.



Modified: trunk/core/yade.cpp
===================================================================
--- trunk/core/yade.cpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/core/yade.cpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -253,8 +253,14 @@
 	int ok = frontEnd->run(argc,argv);
 
 	#ifdef EMBED_PYTHON
+		/* pyFinalize crashes with boost:python<=1.35
+		 * http://www.boost.org/libs/python/todo.html#pyfinalize-safety has explanation 
+		 * once this is fixed, you should remove workaround that saves history from ipython session in gui/py/PythonUI_rc.py:63
+		 *   import IPython.ipapi
+		 *   IPython.ipapi.get().IP.atexit_operations()
+		 */
 		// LOG_DEBUG("Finalizing Python...");
-		// Py_Finalize(); // FIXME: http://www.boost.org/libs/python/todo.html#pyfinalize-safety says this is unsafe with boost::python
+		// Py_Finalize();
 	#endif
 	#ifdef YADE_DEBUG
 		signal(SIGABRT,SIG_DFL); signal(SIGHUP,SIG_DFL); // default handlers

Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/extra/Brefcom.cpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -173,8 +173,8 @@
 		// store Fn (and Fs?), for use with GlobalStiffnessCounter?
 		NNAN(sigmaN); NNANV(sigmaT);
 		NNAN(crossSection);
-		Fn=sigmaN*crossSection;
-		Fs=sigmaT*crossSection;
+		Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
+		Fs=sigmaT*crossSection; BC->shearForce=Fs;
 		applyForce(crossSection*(contGeom->normal*sigmaN + sigmaT),I->getId1(),I->getId2()); /* this is the force applied on the _first_ body; inverted applied to the second */
 	}
 }

Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/extra/Brefcom.hpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -83,6 +83,7 @@
 		bool neverDamage;
 
 		/*! auxiliary variable for visualization, recalculated in BrefcomLaw at every iteration */
+		// FIXME: Fn and Fs are stored as Vector3r normalForce, shearForce in NormalShearInteraction 
 		Real omega, Fn, sigmaN, epsN; Vector3r epsT, sigmaT, Fs;
 
 		BrefcomContact(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), tau(0), expDmgRate(0) { createIndex(); epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; }

Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/extra/clump/Shop.cpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -105,6 +105,31 @@
 #undef __BEX_ACCESS
 
 
+
+Real Shop::unbalancedForce(bool useMaxForce, MetaBody* _rb){
+	MetaBody* rb=_rb ? _rb : Omega::instance().getRootBody().get();
+
+	// get maximum force on a body and sum of all forces (for averaging)
+	Real sumF=0,maxF=0,currF;
+	FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+		if(!b->isDynamic) continue;
+		currF=Shop::Bex::force(b->id,rb).Length(); maxF=max(currF,maxF); sumF+=currF;
+	}
+	Real meanF=sumF/rb->bodies->size(); 
+	// get max force on contacts
+	Real maxContactF=0;
+	FOREACH(const shared_ptr<Interaction>& I, *rb->transientInteractions){
+		if(!I->isReal) continue;
+		shared_ptr<NormalShearInteraction> nsi=YADE_PTR_CAST<NormalShearInteraction>(I->interactionPhysics); assert(nsi);
+		maxContactF=max(maxContactF,(nsi->normalForce+nsi->shearForce).Length());
+	}
+	return (useMaxForce?maxF:meanF)/maxContactF;
+}
+
+
+
+
+
 template <typename valType> valType Shop::getDefault(const string& key) {
 	ensureInit();
 	try{return boost::any_cast<valType>(defaults[key]);}

Modified: trunk/extra/clump/Shop.hpp
===================================================================
--- trunk/extra/clump/Shop.hpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/extra/clump/Shop.hpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -99,4 +99,7 @@
 
 		//! Calculate inscribed circle center of trianlge
 		static Vector3r inscribedCircleCenter(const Vector3r& v0, const Vector3r& v1, const Vector3r& v2);
+
+		//! Get unbalanced force of the whole simulation
+		static Real unbalancedForce(bool useMaxForce=false, MetaBody* _rb=NULL);
 };

Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/gui/py/PythonUI_rc.py	2008-09-24 08:42:04 UTC (rev 1525)
@@ -58,6 +58,10 @@
 	,rc_override={'execfile':[runtime.prefix+'/lib/yade'+runtime.suffix+'/gui/yade/ipython.py']})
 
 	ipshell()
+	# save history -- a workaround for atexit handlers not being run (why?)
+	# http://lists.ipython.scipy.org/pipermail/ipython-user/2008-September/005839.html
+	import IPython.ipapi
+	IPython.ipapi.get().IP.atexit_operations()
 	try:
 		import yade.qt
 		yade.qt.close()

Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/gui/py/_utils.cpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -159,6 +159,7 @@
 	return vec2tuple(Shop::inscribedCircleCenter(Vector3r(python::extract<double>(v0[0]),python::extract<double>(v0[1]),python::extract<double>(v0[2])), Vector3r(python::extract<double>(v1[0]),python::extract<double>(v1[1]),python::extract<double>(v1[2])), Vector3r(python::extract<double>(v2[0]),python::extract<double>(v2[1]),python::extract<double>(v2[2]))));
 }
 
+BOOST_PYTHON_FUNCTION_OVERLOADS(unbalancedForce_overloads,Shop::unbalancedForce,0,1);
 
 BOOST_PYTHON_MODULE(_utils){
 	def("PWaveTimeStep",PWaveTimeStep);
@@ -170,6 +171,7 @@
 	def("bodyNumInteractionsHistogram",bodyNumInteractionsHistogram,bodyNumInteractionsHistogram_overloads(args("aabb")));
 	def("elasticEnergy",elasticEnergyInAABB);
 	def("inscribedCircleCenter",inscribedCircleCenter);
+	def("unbalancedForce",&Shop::unbalancedForce,unbalancedForce_overloads(args("useMaxForce")));
 }
 
 

Modified: trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -9,11 +9,14 @@
 	public:
 		//! normal stiffness
 		Real kn;
+		//! normal force
+		Vector3r normalForce;
 		NormalInteraction(){createIndex(); }
 		virtual ~NormalInteraction();
 	protected:
 		virtual void registerAttributes(){
 			REGISTER_ATTRIBUTE(kn);
+			REGISTER_ATTRIBUTE(normalForce);
 		}
 	REGISTER_CLASS_NAME(NormalInteraction);
 	REGISTER_BASE_CLASS_NAME(InteractionPhysics);
@@ -28,12 +31,15 @@
 	public:
 		//! shear stiffness
 		Real ks;
+		//! shear force
+		Vector3r shearForce;
 		NormalShearInteraction(){ createIndex(); }
 		virtual ~NormalShearInteraction();
 	protected:
 		virtual void registerAttributes(){	
 			NormalInteraction::registerAttributes();
 			REGISTER_ATTRIBUTE(ks);
+			REGISTER_ATTRIBUTE(shearForce);
 		}
 	REGISTER_CLASS_NAME(NormalShearInteraction);
 	REGISTER_BASE_CLASS_NAME(NormalInteraction);

Modified: trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -43,7 +43,7 @@
 		virtual ~CentralGravityEngine(){};
 		virtual void applyCondition(MetaBody*);
 	protected:
-		virtual void registerAttributes(){REGISTER_ATTRIBUTE(centralBody); REGISTER_ATTRIBUTE(accel); REGISTER_ATTRIBUTE(reciprocal); }
+		virtual void registerAttributes(){DeusExMachina::registerAttributes(); REGISTER_ATTRIBUTE(centralBody); REGISTER_ATTRIBUTE(accel); REGISTER_ATTRIBUTE(reciprocal); }
 		NEEDS_BEX("Force");
 		REGISTER_CLASS_NAME(CentralGravityEngine);
 		REGISTER_BASE_CLASS_NAME(DeusExMachina);
@@ -68,7 +68,7 @@
 		virtual ~AxialGravityEngine(){};
 		virtual void applyCondition(MetaBody*);
 	protected:
-		virtual void registerAttributes(){REGISTER_ATTRIBUTE(axisPoint); REGISTER_ATTRIBUTE(axisDirection); REGISTER_ATTRIBUTE(acceleration); }
+		virtual void registerAttributes(){DeusExMachina::registerAttributes(); REGISTER_ATTRIBUTE(axisPoint); REGISTER_ATTRIBUTE(axisDirection); REGISTER_ATTRIBUTE(acceleration); }
 		NEEDS_BEX("Force");
 		REGISTER_CLASS_NAME(AxialGravityEngine);
 		REGISTER_BASE_CLASS_NAME(DeusExMachina);

Modified: trunk/pkg/dem/DataClass/InteractionPhysics/ElasticContactInteraction.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionPhysics/ElasticContactInteraction.hpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/pkg/dem/DataClass/InteractionPhysics/ElasticContactInteraction.hpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -14,9 +14,7 @@
 				,frictionAngle 			// angle of friction, according to Coulumb criterion
 				,tangensOfFrictionAngle;
 	
-		Vector3r	prevNormal			// unit normal of the contact plane.
-				,normalForce			// normal force applied on a DE
-				,shearForce;			// shear force applied on a DE
+		Vector3r	prevNormal;			// unit normal of the contact plane.
 
 		ElasticContactInteraction(){createIndex();};
 		virtual ~ElasticContactInteraction(){};
@@ -24,7 +22,6 @@
 		virtual void registerAttributes(){
 			NormalShearInteraction::registerAttributes();
 			REGISTER_ATTRIBUTE(prevNormal);
-			REGISTER_ATTRIBUTE(shearForce);
 			REGISTER_ATTRIBUTE(initialKn);
 			REGISTER_ATTRIBUTE(initialKs);
 			REGISTER_ATTRIBUTE(tangensOfFrictionAngle);

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -13,6 +13,7 @@
 
 #include<yade/lib-base/yadeWm3Extra.hpp>
 
+CREATE_LOGGER(InteractingFacet2InteractingSphere4SpheresContactGeometry);
 
 InteractingFacet2InteractingSphere4SpheresContactGeometry::InteractingFacet2InteractingSphere4SpheresContactGeometry() 
 {

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp	2008-09-23 13:44:55 UTC (rev 1524)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp	2008-09-24 08:42:04 UTC (rev 1525)
@@ -31,6 +31,8 @@
 	REGISTER_CLASS_NAME(InteractingFacet2InteractingSphere4SpheresContactGeometry);
 	REGISTER_BASE_CLASS_NAME(InteractionGeometryEngineUnit);
 
+	DECLARE_LOGGER;
+
 	FUNCTOR2D(InteractingFacet,InteractingSphere);
 
 	DEFINE_FUNCTOR_ORDER_2D(InteractingFacet,InteractingSphere);