yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-08-07 21:11:33 +0200 (Fri, 07 Aug 2009)
New Revision: 1931

Added:
   trunk/scripts/test/periodic-grow.py
Modified:
   trunk/pkg/dem/meta/Shop.cpp
   trunk/pkg/dem/meta/Shop.hpp
   trunk/py/_utils.cpp
Log:
1. Script that shrinks the periodic cell progressively (reveals bugs in the collider, though ;-) )
2. Add function to compute total force in volume from interactions to Shop, wrapped in _utils as well.


Modified: trunk/pkg/dem/meta/Shop.cpp
===================================================================
--- trunk/pkg/dem/meta/Shop.cpp	2009-08-07 13:22:16 UTC (rev 1930)
+++ trunk/pkg/dem/meta/Shop.cpp	2009-08-07 19:11:33 UTC (rev 1931)
@@ -60,7 +60,12 @@
 
 #include<yade/pkg-dem/Tetra.hpp>
 
+#include<boost/foreach.hpp>
+#ifndef FOREACH
+	#define FOREACH BOOST_FOREACH
+#endif
 
+
 #define _SPEC_CAST(orig,cast) template<> void Shop::setDefault<orig>(string key, orig val){setDefault(key,cast(val));}
 _SPEC_CAST(const char*,string);
 _SPEC_CAST(char*,string);
@@ -80,6 +85,24 @@
 }
 
 
+/*! Compute sum of forces in the whole simulation.
+
+Designed for being used with periodic cell, where diving the resulting components by
+areas of the cell will give average stress in that direction.
+
+Requires all .isReal() interaction to have interactionPhysics deriving from NormalShearInteraction.
+*/
+Vector3r Shop::totalForceInVolume(MetaBody* _rb){
+	MetaBody* rb=_rb ? _rb : Omega::instance().getRootBody().get();
+	Vector3r ret(Vector3r::ZERO);
+	FOREACH(const shared_ptr<Interaction>&I, *rb->interactions){
+		if(!I->isReal()) continue;
+		NormalShearInteraction* nsi=YADE_CAST<NormalShearInteraction*>(I->interactionPhysics.get());
+		ret+=Vector3r(abs(nsi->normalForce[0]+nsi->shearForce[0]),abs(nsi->normalForce[1]+nsi->shearForce[1]),abs(nsi->normalForce[2]+nsi->shearForce[2]));
+	}
+	return ret;
+}
+
 Real Shop::unbalancedForce(bool useMaxForce, MetaBody* _rb){
 	MetaBody* rb=_rb ? _rb : Omega::instance().getRootBody().get();
 	rb->bex.sync();

Modified: trunk/pkg/dem/meta/Shop.hpp
===================================================================
--- trunk/pkg/dem/meta/Shop.hpp	2009-08-07 13:22:16 UTC (rev 1930)
+++ trunk/pkg/dem/meta/Shop.hpp	2009-08-07 19:11:33 UTC (rev 1931)
@@ -94,6 +94,7 @@
 		//! Get unbalanced force of the whole simulation
 		static Real unbalancedForce(bool useMaxForce=false, MetaBody* _rb=NULL);
 		static Real kineticEnergy(MetaBody* _rb=NULL);
+		static Vector3r totalForceInVolume(MetaBody *_rb=NULL);
 
 		//! create transientInteraction between 2 bodies, using existing MetaEngine in Omega
 		static shared_ptr<Interaction> createExplicitInteraction(body_id_t id1, body_id_t id2);

Modified: trunk/py/_utils.cpp
===================================================================
--- trunk/py/_utils.cpp	2009-08-07 13:22:16 UTC (rev 1930)
+++ trunk/py/_utils.cpp	2009-08-07 19:11:33 UTC (rev 1931)
@@ -396,6 +396,7 @@
 
 BOOST_PYTHON_FUNCTION_OVERLOADS(unbalancedForce_overloads,Shop::unbalancedForce,0,1);
 Real Shop__kineticEnergy(){return Shop::kineticEnergy();}
+Vector3r Shop__totalForceInVolume(){return Shop::totalForceInVolume();}
 
 BOOST_PYTHON_MODULE(_utils){
 	// http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
@@ -419,6 +420,7 @@
 	def("sumBexTorques",sumBexTorques);
 	def("forcesOnPlane",forcesOnPlane);
 	def("forcesOnCoordPlane",forcesOnCoordPlane);
+	def("totalForceInVolume",Shop__totalForceInVolume);
 	def("createInteraction",Shop__createExplicitInteraction);
 	def("spiralProject",spiralProject,spiralProject_overloads(args("axis","periodStart","theta0")));
 	def("pointInsidePolygon",pointInsidePolygon);

Added: trunk/scripts/test/periodic-grow.py
===================================================================
--- trunk/scripts/test/periodic-grow.py	2009-08-07 13:22:16 UTC (rev 1930)
+++ trunk/scripts/test/periodic-grow.py	2009-08-07 19:11:33 UTC (rev 1931)
@@ -0,0 +1,47 @@
+"""Script that either grows spheres inside the cell or shrinks
+the cell progressively. It prints the total volume force once in a while.
+This script also shows that the collider misses some interactions as spheres
+are getting one over another. That is not acceptable, of course. """
+from yade import log,timing
+
+O.engines=[
+	BexResetter(),
+	BoundingVolumeMetaEngine([InteractingSphere2AABB(),MetaInteractingGeometry2AABB()]),
+	PeriodicInsertionSortCollider(label='collider'),
+	InteractionDispatchers(
+		[ef2_Sphere_Sphere_Dem3DofGeom()],
+		[SimpleElasticRelationships()],
+		[Law2_Dem3Dof_Elastic_Elastic()],
+	),
+	NewtonsDampedLaw(damping=.6)
+]
+import random
+O.bodies.append(utils.sphere((0,0,0),.5,dynamic=False,density=1000)) # stable point
+for i in xrange(150):
+	O.bodies.append(utils.sphere(Vector3(10*random.random(),10*random.random(),10*random.random()),.2+.2*random.random(),density=1000))
+O.periodicCell=((-5,-5,0),(5,5,10))
+O.dt=.8*utils.PWaveTimeStep()
+O.saveTmp()
+from yade import qt
+qt.Controller(); qt.View()
+step=.01
+O.run(200,True)
+if 0:
+	for i in range(0,500):
+		O.run(500,True)
+		for b in O.bodies:
+			b.shape['radius']=b.shape['radius']+step
+			b.mold['radius']=b.mold['radius']+step
+		for i in O.interactions:
+			if not i.isReal: continue
+			i.geom['effR1']=i.geom['effR1']+step
+			i.geom['effR2']=i.geom['effR2']+step
+		print O.iter,utils.totalForceInVolume()
+else:
+	for i in range(0,500):
+		O.run(100,True)
+		mn,mx=O.periodicCell
+		step=(mx-mn); step=Vector3(.002*step[0],.002*step[1],.002*step[2])
+		O.periodicCell=mn+step,mx-step
+		if (i%10==0): print O.iter,utils.totalForceInVolume()
+#O.timingEnabled=True; timing.reset(); O.run(200000,True); timing.stats()