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()
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
_______________________________________________
yade-dev mailing list
[email protected]
https://lists.berlios.de/mailman/listinfo/yade-dev
