yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01632
[svn] r1933 - in trunk: extra lib pkg/common/Engine/DeusExMachina pkg/dem/DataClass pkg/dem/meta py py/yadeWrapper scripts/test
Author: eudoxos
Date: 2009-08-10 13:32:23 +0200 (Mon, 10 Aug 2009)
New Revision: 1933
Added:
trunk/pkg/dem/DataClass/SpherePack.cpp
trunk/pkg/dem/DataClass/SpherePack.hpp
trunk/scripts/test/periodic-compress.py
Modified:
trunk/extra/PeriodicInsertionSortCollider.cpp
trunk/extra/PeriodicInsertionSortCollider.hpp
trunk/extra/SConscript
trunk/lib/SConscript
trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp
trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp
trunk/pkg/dem/meta/Shop.cpp
trunk/pkg/dem/meta/Shop.hpp
trunk/py/SConscript
trunk/py/_packSpheres.cpp
trunk/py/_utils.cpp
trunk/py/yadeWrapper/yadeWrapper.cpp
trunk/scripts/test/periodic-grow.py
Log:
1. Add PeriIsoCompressor for periodic isotropic compression; see scripts/test/periodic-compress.py. Not so much slower than regular triax, nice surprise.
2. Make SpherePack standalone class (previously just in python's _packSpheres), add the function to generate cloud of spheres, both periodic and non-periodic. _packSpheres wrapper is update, but otherwise not changed.
3. Fix python-less compilation for se3 interpolator engine
4. Add Omega().runEngine()
5. Fix some compilation/linking issues.
Modified: trunk/extra/PeriodicInsertionSortCollider.cpp
===================================================================
--- trunk/extra/PeriodicInsertionSortCollider.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/extra/PeriodicInsertionSortCollider.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -16,7 +16,7 @@
using namespace std;
-YADE_PLUGIN((PeriodicInsertionSortCollider))
+YADE_PLUGIN((PeriodicInsertionSortCollider)(PeriIsoCompressor))
CREATE_LOGGER(PeriodicInsertionSortCollider);
// return floating value wrapped between x0 and x1 and saving period number to period
@@ -348,3 +348,70 @@
}
}
}
+
+
+
+
+
+/***********************************************************************************************************
+******************************* PeriIsoCompressor *********************************************************/
+#include<yade/pkg-dem/Shop.hpp>
+
+
+CREATE_LOGGER(PeriIsoCompressor);
+void PeriIsoCompressor::action(MetaBody* rb){
+ if(!rb->isPeriodic){ LOG_FATAL("Being used on non-periodic simulation!"); throw; }
+ if(state>=stresses.size()) return;
+ // initialize values
+ if(charLen<0){
+ BoundingVolume* bv=Body::byId(0,rb)->boundingVolume.get();
+ if(!bv){ LOG_FATAL("No charLen defined and body #0 has no boundingVolume"); throw; }
+ const Vector3r sz=bv->max-bv->min;
+ charLen=(sz[0]+sz[1]+sz[2])/3.;
+ LOG_INFO("No charLen defined, taking avg bbox size of body #0 = "<<charLen);
+ }
+ if(maxDisplPerStep<0) maxDisplPerStep=1e-2*charLen; // this should be tuned somehow…
+ const long& step=rb->currentIteration;
+ Vector3r cellSize=rb->cellMax-rb->cellMin; //unused: Real cellVolume=cellSize[0]*cellSize[1]*cellSize[2];
+ Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
+ if(((step%globalUpdateInt)==0) || avgStiffness<0 || sigma[0]<0 || sigma[1]<0 || sigma[2]<0){
+ Vector3r sumForces=Shop::totalForceInVolume(avgStiffness,rb);
+ sigma=Vector3r(sumForces[0]/cellArea[0],sumForces[1]/cellArea[1],sumForces[2]/cellArea[2]);
+ LOG_TRACE("Updated sigma="<<sigma<<", avgStiffness="<<avgStiffness);
+ }
+ Real sigmaGoal=stresses[state]; assert(sigmaGoal>0);
+ // expansion of cell in this step (absolute length)
+ Vector3r cellGrow(Vector3r::ZERO);
+ // is the stress condition satisfied in all directions?
+ bool allStressesOK=true;
+ Vector3r cg_;
+ for(int axis=0; axis<3; axis++){
+ // Δσ=ΔεE=(Δl/l)×(l×K/A) ↔ Δl=Δσ×A/K
+ // FIXME: either NormalShearInteraction::{kn,ks} is computed wrong or we have dimensionality problem here
+ // FIXME: that is why the fixup 1e-2 is needed here
+ // FIXME: or perhaps maxDisplaPerStep=1e-2*charLen is too big??
+ cellGrow[axis]=1e-4*(sigma[axis]-sigmaGoal)*cellArea[axis]/(avgStiffness>0?avgStiffness:1);
+ if(abs(cellGrow[axis])>maxDisplPerStep) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxDisplPerStep;
+ // crude way of predicting sigma, for steps when it is not computed from intrs
+ if(avgStiffness>0) sigma[axis]-=cellGrow[axis]*avgStiffness;
+ if(abs((sigma[axis]-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
+ }
+ TRVAR4(cellGrow,sigma,sigmaGoal,avgStiffness);
+ rb->cellMin-=.5*cellGrow; rb->cellMax+=.5*cellGrow;
+ // handle state transitions
+ if(allStressesOK){
+ if((step%globalUpdateInt)==0) currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,rb);
+ if(currUnbalanced<maxUnbalanced){
+ state+=1;
+ // sigmaGoal reached and packing stable
+ if(state==stresses.size()){ // no next stress to go for
+ LOG_INFO("Finished");
+ #ifdef YADE_PYTHON
+ if(!doneHook.empty()){ LOG_DEBUG("Running doneHook: "<<doneHook); PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(doneHook.c_str()); PyGILState_Release(gstate); }
+ #endif
+ } else { LOG_INFO("Loading to "<<sigmaGoal<<" done, starting going to "<<stresses[state]<<" now"); }
+ } else {
+ if((step%globalUpdateInt)==0) LOG_DEBUG("Stress="<<sigma<<", goal="<<sigmaGoal<<", unbalanced="<<currUnbalanced);
+ }
+ }
+}
Modified: trunk/extra/PeriodicInsertionSortCollider.hpp
===================================================================
--- trunk/extra/PeriodicInsertionSortCollider.hpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/extra/PeriodicInsertionSortCollider.hpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -138,3 +138,29 @@
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(PeriodicInsertionSortCollider);
+
+/* Simple engine for compressing random cloud of spheres until
+it reaches given average stress and then goes back to given residual stress
+*/
+class PeriIsoCompressor: public StandAloneEngine{
+ Real avgStiffness; Real maxDisplPerStep; Vector3r sumForces, sigma;
+ public:
+ //! Stresses that should be reached, one after another
+ vector<Real> stresses;
+ //! Characteristic length, should be something like mean particle diameter (default -1=invalid value))
+ Real charLen;
+ //! if actual unbalanced force is smaller than this number, the packing is considered stable (default 1e-4)
+ Real maxUnbalanced, currUnbalanced;
+ //! how often to recompute average stress, stiffness and unbalanced force (default 100)
+ int globalUpdateInt;
+ //! Where are we at in the process
+ size_t state;
+ //! python command to be run when reaching the last specified stress
+ string doneHook;
+ void action(MetaBody*);
+ PeriIsoCompressor(): avgStiffness(-1), maxDisplPerStep(-1), sumForces(Vector3r::ZERO), sigma(Vector3r::ZERO), charLen(-1), maxUnbalanced(1e-4), currUnbalanced(-1), globalUpdateInt(100), state(0){}
+ REGISTER_ATTRIBUTES(StandAloneEngine,(stresses)(charLen)(maxUnbalanced)(globalUpdateInt)(state)(doneHook));
+ DECLARE_LOGGER;
+ REGISTER_CLASS_AND_BASE(PeriIsoCompressor,StandAloneEngine);
+};
+REGISTER_SERIALIZABLE(PeriIsoCompressor);
Modified: trunk/extra/SConscript
===================================================================
--- trunk/extra/SConscript 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/extra/SConscript 2009-08-10 11:32:23 UTC (rev 1933)
@@ -6,7 +6,7 @@
import os.path, os
env.Install('$PREFIX/lib/yade$SUFFIX/extra',[
- env.SharedLibrary('PeriodicInsertionSortCollider',['PeriodicInsertionSortCollider.cpp'],LIBS=env['LIBS']+linkPlugins(['InsertionSortCollider'])),
+ env.SharedLibrary('PeriodicInsertionSortCollider',['PeriodicInsertionSortCollider.cpp'],LIBS=env['LIBS']+linkPlugins(['InsertionSortCollider','Shop'])),
])
Modified: trunk/lib/SConscript
===================================================================
--- trunk/lib/SConscript 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/lib/SConscript 2009-08-10 11:32:23 UTC (rev 1933)
@@ -59,7 +59,7 @@
'factory/ClassFactory.cpp','factory/DynLibManager.cpp','factory/FactoryExceptions.cpp',
'multimethods/Indexable.cpp','multimethods/MultiMethodsExceptions.cpp',
'serialization-xml/XMLFormatManager.cpp','serialization-xml/XMLSaxParser.cpp','serialization/Archive.cpp',
- 'serialization/IOFormatManager.cpp','serialization/IOManagerExceptions.cpp','serialization/FormatChecker.cpp','serialization/Serializable.cpp','serialization/SerializableSingleton.cpp','serialization/SerializationExceptions.cpp']))
+ 'serialization/IOFormatManager.cpp','serialization/IOManagerExceptions.cpp','serialization/FormatChecker.cpp','serialization/Serializable.cpp','serialization/SerializableSingleton.cpp','serialization/SerializationExceptions.cpp']),LIBS=['dl'])
#######################
###### 3rd party libs
Modified: trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -31,6 +31,8 @@
if(t>=1.){
done=true;
LOG_DEBUG("Goal reached.");
- if(!goalHook.empty()){ PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(goalHook.c_str()); PyGILState_Release(gstate); }
+ #ifdef YADE_PYTHON
+ if(!goalHook.empty()){ PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(goalHook.c_str()); PyGILState_Release(gstate); }
+ #endif
}
}
Modified: trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -4,7 +4,9 @@
#include<yade/core/DeusExMachina.hpp>
#include<yade/core/PhysicalParameters.hpp>
#include<yade/pkg-common/RigidBodyParameters.hpp>
-#include<Python.h>
+#ifdef YADE_PYTHON
+ #include<Python.h>
+#endif
/* Engine interpolating between starting (current) and goal (given) se3, both position and orientation.
*
Added: trunk/pkg/dem/DataClass/SpherePack.cpp
===================================================================
--- trunk/pkg/dem/DataClass/SpherePack.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/pkg/dem/DataClass/SpherePack.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -0,0 +1,113 @@
+// © 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
+
+#include<yade/pkg-dem/SpherePack.hpp>
+
+#include<yade/core/Omega.hpp>
+#include<yade/core/MetaBody.hpp>
+#include<yade/pkg-common/InteractingSphere.hpp>
+#include<yade/pkg-dem/Shop.hpp>
+
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/variate_generator.hpp>
+
+#include<yade/core/Timing.hpp>
+
+// not a serializable in the sense of YADE_PLUGIN
+
+CREATE_LOGGER(SpherePack);
+
+using namespace std;
+using namespace boost;
+
+#ifdef YADE_PYTHON
+ // if we need explicit conversions at a few places
+ python::tuple vec2tuple(const Vector3r& v){return boost::python::make_tuple(v[0],v[1],v[2]);}
+ Vector3r tuple2vec(const python::tuple& t){return Vector3r(python::extract<double>(t[0])(),python::extract<double>(t[1])(),python::extract<double>(t[2])());}
+ Vector3r tuple2vec(python::tuple& t){return Vector3r(python::extract<double>(t[0])(),python::extract<double>(t[1])(),python::extract<double>(t[2])());}
+
+ void SpherePack::fromList(const python::list& l){
+ pack.clear();
+ size_t len=python::len(l);
+ for(size_t i=0; i<len; i++){
+ const python::tuple& t=python::extract<python::tuple>(l[i]);
+ //python::extract<python::tuple> tup(t[0]);
+ //if(tup.check()) { pack.push_back(Sph(tuple2vec(tup()),python::extract<double>(t[1]))); continue;}
+ python::extract<Vector3r> vec(t[0]);
+ if(vec.check()) { pack.push_back(Sph(vec(),python::extract<double>(t[1]))); continue; }
+ PyErr_SetString(PyExc_TypeError, "List elements must be (tuple or Vector3, float)!");
+ python::throw_error_already_set();
+ }
+ };
+
+ python::list SpherePack::toList() const {
+ python::list ret;
+ FOREACH(const Sph& s, pack) ret.append(python::make_tuple(s.c,s.r));
+ return ret;
+ };
+
+ python::list SpherePack::toList_pointsAsTuples() const {
+ python::list ret;
+ FOREACH(const Sph& s, pack) ret.append(python::make_tuple(vec2tuple(s.c),s.r));
+ return ret;
+ };
+#endif
+
+void SpherePack::fromFile(string file) {
+ typedef pair<Vector3r,Real> pairVector3rReal;
+ vector<pairVector3rReal> ss;
+ Vector3r mn,mx;
+ ss=Shop::loadSpheresFromFile(file,mn,mx);
+ pack.clear();
+ FOREACH(const pairVector3rReal& s, ss) pack.push_back(Sph(s.first,s.second));
+}
+
+void SpherePack::toFile(const string fname) const {
+ ofstream f(fname.c_str());
+ if(!f.good()) throw runtime_error("Unable to open file `"+fname+"'");
+ FOREACH(const Sph& s, pack) f<<s.c[0]<<" "<<s.c[1]<<" "<<s.c[2]<<" "<<s.r<<endl;
+ f.close();
+};
+
+void SpherePack::fromSimulation() {
+ pack.clear();
+ MetaBody* rootBody=Omega::instance().getRootBody().get();
+ FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
+ shared_ptr<InteractingSphere> intSph=dynamic_pointer_cast<InteractingSphere>(b->interactingGeometry);
+ if(!intSph) continue;
+ pack.push_back(Sph(b->physicalParameters->se3.position,intSph->radius));
+ }
+}
+
+long SpherePack::makeCloud(Vector3r mn, Vector3r mx, Real rMean, Real rRelFuzz, size_t num, bool periodic){
+ static boost::minstd_rand randGen(TimingInfo::getNow(/* get the number even if timing is disabled globally */ true));
+ static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > rnd(randGen, boost::uniform_real<>(0,1));
+ const size_t maxTry=1000;
+ Vector3r size=mx-mn;
+ if(periodic)(cellSize=size);
+ for(size_t i=0; i<num; i++) {
+ size_t t;
+ for(t=0; t<maxTry; ++t){
+ Real r=(rnd()-.5)*rRelFuzz+rMean; Vector3r c;
+ if(!periodic) { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+r+(size[axis]-2*r)*rnd(); }
+ else { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+size[axis]*rnd(); }
+ size_t packSize=pack.size(); bool overlap=false;
+ if(!periodic){
+ for(size_t j=0; j<packSize; j++){ if(pow(pack[j].r+r,2) >= (pack[j].c-c).SquaredLength()) { overlap=true; break; } }
+ } else {
+ for(size_t j=0; j<packSize; j++){
+ Vector3r dr;
+ for(int axis=0; axis<3; axis++) dr[axis]=min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis]));
+ if(pow(pack[j].r+r,2)>= dr.SquaredLength()){ overlap=true; break; }
+ }
+ }
+ if(!overlap) { pack.push_back(Sph(c,r)); break; }
+ }
+ if (t==maxTry) {
+ LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num);
+ return i;
+ }
+ }
+ return pack.size();
+}
+
Added: trunk/pkg/dem/DataClass/SpherePack.hpp
===================================================================
--- trunk/pkg/dem/DataClass/SpherePack.hpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/pkg/dem/DataClass/SpherePack.hpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -0,0 +1,100 @@
+// © 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
+
+#include<vector>
+#include<string>
+using namespace std; // sorry
+
+#ifdef YADE_PYTHON
+ #include<boost/python.hpp>
+ #include<yade/extra/boost_python_len.hpp>
+ using namespace boost;
+#endif
+
+#include<boost/foreach.hpp>
+#ifndef FOREACH
+ #define FOREACH BOOST_FOREACH
+#endif
+
+#include<yade/lib-base/Logging.hpp>
+#include<yade/lib-base/yadeWm3Extra.hpp>
+
+/*! Class representing geometry of spherical packing, with some utility functions.
+
+*/
+class SpherePack{
+ // return coordinate wrapped to x0…x1, relative to x0; don't care about period
+ // copied from PeriodicInsertionSortCollider
+ Real cellWrapRel(const Real x, const Real x0, const Real x1){
+ Real xNorm=(x-x0)/(x1-x0);
+ return (xNorm-floor(xNorm))*(x1-x0);
+ }
+public:
+ struct Sph{
+ Vector3r c; Real r;
+ Sph(const Vector3r& _c, Real _r): c(_c), r(_r){};
+ #ifdef YADE_PYTHON
+ python::tuple asTuple() const { return python::make_tuple(c,r); }
+ #endif
+ };
+ std::vector<Sph> pack;
+ Vector3r cellSize;
+ SpherePack(): cellSize(Vector3r::ZERO){};
+ #ifdef YADE_PYTHON
+ SpherePack(const python::list& l):cellSize(Vector3r::ZERO){ fromList(l); }
+ #endif
+ // add single sphere
+ void add(const Vector3r& c, Real r){ pack.push_back(Sph(c,r)); }
+
+ // I/O
+ #ifdef YADE_PYTHON
+ void fromList(const python::list& l);
+ python::list toList() const;
+ python::list toList_pointsAsTuples() const;
+ #endif
+ void fromFile(const string file);
+ void toFile(const string file) const;
+ void fromSimulation();
+
+ // random generation
+ long makeCloud(Vector3r min, Vector3r max, Real rMean, Real rFuzz, size_t num, bool periodic=false);
+
+ // spatial characteristics
+ Vector3r dim() const {Vector3r mn,mx; aabb(mn,mx); return mx-mn;}
+ #ifdef YADE_PYTHON
+ python::tuple aabb_py() const { Vector3r mn,mx; aabb(mn,mx); return python::make_tuple(mn,mx); }
+ #endif
+ void aabb(Vector3r& mn, Vector3r& mx) const {
+ Real inf=std::numeric_limits<Real>::infinity(); mn=Vector3r(inf,inf,inf); mx=Vector3r(-inf,-inf,-inf);
+ FOREACH(const Sph& s, pack){ Vector3r r(s.r,s.r,s.r); mn=componentMinVector(mn,s.c-r); mx=componentMaxVector(mx,s.c+r);}
+ }
+ Vector3r midPt() const {Vector3r mn,mx; aabb(mn,mx); return .5*(mn+mx);}
+ Real relDensity() const {
+ Real sphVol=0; Vector3r dd=dim();
+ FOREACH(const Sph& s, pack) sphVol+=pow(s.r,3);
+ sphVol*=(4/3.)*Mathr::PI;
+ return sphVol/(dd[0]*dd[1]*dd[2]);
+ }
+
+ // transformations
+ void translate(const Vector3r& shift){ FOREACH(Sph& s, pack) s.c+=shift; }
+ void rotate(const Vector3r& axis, Real angle){ Vector3r mid=midPt(); Quaternionr q(axis,angle); FOREACH(Sph& s, pack) s.c=q*(s.c-mid)+mid; }
+ void scale(Real scale){ Vector3r mid=midPt(); FOREACH(Sph& s, pack) {s.c=scale*(s.c-mid)+mid; s.r*=abs(scale); } }
+
+ // iteration
+ #ifdef YADE_PYTHON
+ size_t len() const{ return pack.size(); }
+ python::tuple getitem(size_t idx){ if(idx<0 || idx>=pack.size()) throw runtime_error("Index "+lexical_cast<string>(idx)+" out of range 0.."+lexical_cast<string>(pack.size()-1)); return pack[idx].asTuple(); }
+ struct iterator{
+ const SpherePack& sPack; size_t pos;
+ iterator(const SpherePack& _sPack): sPack(_sPack), pos(0){}
+ iterator iter(){ return *this;}
+ python::tuple next(){
+ if(pos==sPack.pack.size()){ PyErr_SetNone(PyExc_StopIteration); python::throw_error_already_set(); }
+ return sPack.pack[pos++].asTuple();
+ }
+ };
+ SpherePack::iterator getIterator() const{ return SpherePack::iterator(*this);};
+ #endif
+ DECLARE_LOGGER;
+};
+
Modified: trunk/pkg/dem/meta/Shop.cpp
===================================================================
--- trunk/pkg/dem/meta/Shop.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/pkg/dem/meta/Shop.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -85,22 +85,25 @@
}
-/*! Compute sum of forces in the whole simulation.
+/*! Compute sum of forces in the whole simulation and averages stiffness.
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){
+Vector3r Shop::totalForceInVolume(Real& avgIsoStiffness, MetaBody* _rb){
MetaBody* rb=_rb ? _rb : Omega::instance().getRootBody().get();
- Vector3r ret(Vector3r::ZERO);
+ Vector3r force(Vector3r::ZERO); Real stiff=0; long n=0;
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]));
+ force+=Vector3r(abs(nsi->normalForce[0]+nsi->shearForce[0]),abs(nsi->normalForce[1]+nsi->shearForce[1]),abs(nsi->normalForce[2]+nsi->shearForce[2]));
+ stiff+=(1/3.)*nsi->kn+(2/3.)*nsi->ks; // count kn in one direction and ks in the other two
+ n++;
}
- return ret;
+ avgIsoStiffness= n>0 ? (1./n)*stiff : -1;
+ return force;
}
Real Shop::unbalancedForce(bool useMaxForce, MetaBody* _rb){
@@ -1167,3 +1170,4 @@
return x0+xxNorm*(x1-x0);
}
+
Modified: trunk/pkg/dem/meta/Shop.hpp
===================================================================
--- trunk/pkg/dem/meta/Shop.hpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/pkg/dem/meta/Shop.hpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -94,7 +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);
+ static Vector3r totalForceInVolume(Real& avgIsoStiffness, 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/SConscript
===================================================================
--- trunk/py/SConscript 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/py/SConscript 2009-08-10 11:32:23 UTC (rev 1933)
@@ -18,7 +18,7 @@
LIBS=env['LIBS']+(['_gts__python-module'] if 'GTS' in env['features'] else []),
),
env.SharedLibrary('_packSpheres',['_packSpheres.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+[
- linkPlugins(['Shop']),
+ linkPlugins(['Shop','SpherePack']),
]),
env.SharedLibrary('_packObb',['_packObb.cpp'],SHLIBPREFIX=''),
env.File('utils.py'),
Modified: trunk/py/_packSpheres.cpp
===================================================================
--- trunk/py/_packSpheres.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/py/_packSpheres.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -1,132 +1,7 @@
// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
-#include<boost/python.hpp>
-#include<yade/extra/boost_python_len.hpp>
-#include<yade/lib-base/Logging.hpp>
-#include<yade/lib-base/yadeWm3.hpp>
-#include<yade/lib-base/yadeWm3Extra.hpp>
-// #include<yade/gui-py/_utils.hpp> // will be: yade/lib-py/_utils.hpp> at some point
-#include<Wm3Vector3.h>
-#include<yade/core/Omega.hpp>
-#include<yade/core/MetaBody.hpp>
-#include<yade/pkg-common/InteractingSphere.hpp>
+#include<yade/pkg-dem/SpherePack.hpp>
-#include<yade/pkg-dem/Shop.hpp>
-
-using namespace boost;
-using namespace std;
-#ifdef LOG4CXX
- log4cxx::LoggerPtr logger=log4cxx::Logger::getLogger("yade.pack.spheres");
-#endif
-
-// copied from _packPredicates.cpp :-(
-python::tuple vec2tuple(const Vector3r& v){return boost::python::make_tuple(v[0],v[1],v[2]);}
-Vector3r tuple2vec(const python::tuple& t){return Vector3r(python::extract<double>(t[0])(),python::extract<double>(t[1])(),python::extract<double>(t[2])());}
-Vector3r tuple2vec(python::tuple& t){return Vector3r(python::extract<double>(t[0])(),python::extract<double>(t[1])(),python::extract<double>(t[2])());}
-
-struct SpherePack{
- struct Sph{ Vector3r c; Real r; Sph(const Vector3r& _c, Real _r): c(_c), r(_r){};
- python::tuple asTuple() const { return python::make_tuple(c,r); }
- };
- vector<Sph> pack;
- SpherePack(){};
- SpherePack(const python::list& l){ fromList(l); }
- // add single sphere
- void add(const Vector3r& c, Real r){ pack.push_back(Sph(c,r)); }
- // I/O
- void fromList(const python::list& l);
- python::list toList() const;
- python::list toList_pointsAsTuples() const;
- void fromFile(const string file);
- void toFile(const string file) const;
- void fromSimulation();
- // spatial characteristics
- Vector3r dim() const {Vector3r mn,mx; aabb(mn,mx); return mx-mn;}
- python::tuple aabb_py() const { Vector3r mn,mx; aabb(mn,mx); return python::make_tuple(mn,mx); }
- void aabb(Vector3r& mn, Vector3r& mx) const {
- Real inf=std::numeric_limits<Real>::infinity(); mn=Vector3r(inf,inf,inf); mx=Vector3r(-inf,-inf,-inf);
- FOREACH(const Sph& s, pack){ Vector3r r(s.r,s.r,s.r); mn=componentMinVector(mn,s.c-r); mx=componentMaxVector(mx,s.c+r);}
- }
- Vector3r midPt() const {Vector3r mn,mx; aabb(mn,mx); return .5*(mn+mx);}
- Real relDensity() const {
- Real sphVol=0; Vector3r dd=dim();
- FOREACH(const Sph& s, pack) sphVol+=pow(s.r,3);
- sphVol*=(4/3.)*Mathr::PI;
- return sphVol/(dd[0]*dd[1]*dd[2]);
- }
- // transformations
- void translate(const Vector3r& shift){ FOREACH(Sph& s, pack) s.c+=shift; }
- void rotate(const Vector3r& axis, Real angle){ Vector3r mid=midPt(); Quaternionr q(axis,angle); FOREACH(Sph& s, pack) s.c=q*(s.c-mid)+mid; }
- void scale(Real scale){ Vector3r mid=midPt(); FOREACH(Sph& s, pack) {s.c=scale*(s.c-mid)+mid; s.r*=abs(scale); } }
- // iteration
- size_t len() const{ return pack.size(); }
- python::tuple getitem(size_t idx){ if(idx<0 || idx>=pack.size()) throw runtime_error("Index "+lexical_cast<string>(idx)+" out of range 0.."+lexical_cast<string>(pack.size()-1)); return pack[idx].asTuple(); }
- struct iterator{
- const SpherePack& sPack; size_t pos;
- iterator(const SpherePack& _sPack): sPack(_sPack), pos(0){}
- iterator iter(){ return *this;}
- python::tuple next(){
- if(pos==sPack.pack.size()-1){ PyErr_SetNone(PyExc_StopIteration); python::throw_error_already_set(); }
- return sPack.pack[pos++].asTuple();
- }
- };
- SpherePack::iterator getIterator() const{ return SpherePack::iterator(*this);};
-};
-
-void SpherePack::fromList(const python::list& l){
- pack.clear();
- size_t len=python::len(l);
- for(size_t i=0; i<len; i++){
- const python::tuple& t=python::extract<python::tuple>(l[i]);
- python::extract<python::tuple> tup(t[0]);
- if(tup.check()) { pack.push_back(Sph(tuple2vec(tup()),python::extract<double>(t[1]))); continue;}
- python::extract<Vector3r> vec(t[0]);
- if(vec.check()) { pack.push_back(Sph(vec(),python::extract<double>(t[1]))); continue; }
- PyErr_SetString(PyExc_TypeError, "List elements must be (tuple or Vector3, float)!");
- python::throw_error_already_set();
- }
-};
-
-python::list SpherePack::toList() const {
- python::list ret;
- FOREACH(const Sph& s, pack) ret.append(python::make_tuple(s.c,s.r));
- return ret;
-};
-
-python::list SpherePack::toList_pointsAsTuples() const {
- python::list ret;
- FOREACH(const Sph& s, pack) ret.append(python::make_tuple(vec2tuple(s.c),s.r));
- return ret;
-};
-
-void SpherePack::fromFile(string file) {
- typedef pair<Vector3r,Real> pairVector3rReal;
- vector<pairVector3rReal> ss;
- Vector3r mn,mx;
- ss=Shop::loadSpheresFromFile(file,mn,mx);
- pack.clear();
- FOREACH(const pairVector3rReal& s, ss) pack.push_back(Sph(s.first,s.second));
-}
-
-void SpherePack::toFile(const string fname) const {
- ofstream f(fname.c_str());
- if(!f.good()) throw runtime_error("Unable to open file `"+fname+"'");
- FOREACH(const Sph& s, pack) f<<s.c[0]<<" "<<s.c[1]<<" "<<s.c[2]<<" "<<s.r<<endl;
- f.close();
-};
-
-void SpherePack::fromSimulation() {
- pack.clear();
- MetaBody* rootBody=Omega::instance().getRootBody().get();
- FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
- shared_ptr<InteractingSphere> intSph=dynamic_pointer_cast<InteractingSphere>(b->interactingGeometry);
- if(!intSph) continue;
- pack.push_back(Sph(b->physicalParameters->se3.position,intSph->radius));
- }
-}
-
-
-
BOOST_PYTHON_MODULE(_packSpheres){
python::class_<SpherePack>("SpherePack","Set of spheres as centers and radii",python::init<python::optional<python::list> >(python::args("list"),"Empty constructor, optionally taking list [ ((cx,cy,cz),r), … ] for initial data." ))
.def("add",&SpherePack::add,"Add single sphere to packing, given center as 3-tuple and radius")
@@ -136,9 +11,11 @@
.def("load",&SpherePack::fromFile,"Load packing from external text file (current data will be discarded).")
.def("save",&SpherePack::toFile,"Save packing to external text file (will be overwritten).")
.def("fromSimulation",&SpherePack::fromSimulation,"Make packing corresponding to the current simulation. Discards current data.")
+ .def("makeCloud",&SpherePack::makeCloud,python::args("minCorner","maxCorner","rMean","rRelFuzz","num","periodic"),"Create random packing encosed in box given by minCorner and maxCorner, containing num spheres. Returns number of created spheres, which can be < num if the packing is too tight.")
.def("aabb",&SpherePack::aabb_py,"Get axis-aligned bounding box coordinates, as 2 3-tuples.")
.def("dim",&SpherePack::dim,"Return dimensions of the packing in terms of aabb(), as a 3-tuple.")
.def("center",&SpherePack::midPt,"Return coordinates of the bounding box center.")
+ .def_readonly("cellSize",&SpherePack::cellSize,"Size of periodic cell; is Vector3(0,0,0) if not periodic.")
.def("relDensity",&SpherePack::relDensity,"Relative packing density, measured as sum of spheres' volumes / aabb volume.\n(Sphere overlaps are ignored.)")
.def("translate",&SpherePack::translate,"Translate all spheres by given vector.")
.def("rotate",&SpherePack::rotate,"Rotate all spheres around packing center (in terms of aabb()), given axis and angle of the rotation.")
Modified: trunk/py/_utils.cpp
===================================================================
--- trunk/py/_utils.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/py/_utils.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -396,7 +396,7 @@
BOOST_PYTHON_FUNCTION_OVERLOADS(unbalancedForce_overloads,Shop::unbalancedForce,0,1);
Real Shop__kineticEnergy(){return Shop::kineticEnergy();}
-Vector3r Shop__totalForceInVolume(){return Shop::totalForceInVolume();}
+python::tuple Shop__totalForceInVolume(){Real stiff; Vector3r ret=Shop::totalForceInVolume(stiff); return python::make_tuple(ret,stiff); }
BOOST_PYTHON_MODULE(_utils){
// http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
@@ -420,7 +420,7 @@
def("sumBexTorques",sumBexTorques);
def("forcesOnPlane",forcesOnPlane);
def("forcesOnCoordPlane",forcesOnCoordPlane);
- def("totalForceInVolume",Shop__totalForceInVolume);
+ def("totalForceInVolume",Shop__totalForceInVolume,"Return summed forces on all interactions and average isotropic stiffness, as tuple (Vector3,float)");
def("createInteraction",Shop__createExplicitInteraction);
def("spiralProject",spiralProject,spiralProject_overloads(args("axis","periodStart","theta0")));
def("pointInsidePolygon",pointInsidePolygon);
Modified: trunk/py/yadeWrapper/yadeWrapper.cpp
===================================================================
--- trunk/py/yadeWrapper/yadeWrapper.cpp 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/py/yadeWrapper/yadeWrapper.cpp 2009-08-10 11:32:23 UTC (rev 1933)
@@ -462,6 +462,7 @@
signal(SIGSEGV,termHandler); /* unset the handler that runs gdb and prints backtrace */
exit(status);
}
+ void runEngine(const shared_ptr<Engine>& e){ e->action(OMEGA.getRootBody().get()); }
};
BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(omega_run_overloads,run,0,2);
BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(omega_saveTmp_overloads,saveTmp,0,1);
@@ -648,6 +649,7 @@
#ifdef YADE_BOOST_SERIALIZATION
.def("saveXML",&pyOmega::saveXML,"[EXPERIMENTAL] function saving to XML file using boost::serialization.")
#endif
+ .def("runEngine",&pyOmega::runEngine,"Run given engine exactly once; simulation time, step number etc. will not be incremented (use only if you know what you do).")
;
python::class_<pyTags>("TagsWrapper",python::init<pyTags&>())
.def("__getitem__",&pyTags::getItem)
Added: trunk/scripts/test/periodic-compress.py
===================================================================
--- trunk/scripts/test/periodic-compress.py 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/scripts/test/periodic-compress.py 2009-08-10 11:32:23 UTC (rev 1933)
@@ -0,0 +1,27 @@
+O.periodicCell=(0,0,0),(20,20,20)
+from yade import pack,log,timing
+p=pack.SpherePack()
+p.makeCloud(O.periodicCell[0],O.periodicCell[1],1,.5,700,True)
+for sph in p:
+ O.bodies.append(utils.sphere(sph[0],sph[1],density=1000,dynamic=True))
+
+#log.setLevel("PeriIsoCompressor",log.TRACE)
+O.timingEnabled=True
+O.engines=[
+ BexResetter(),
+ BoundingVolumeMetaEngine([InteractingSphere2AABB(),MetaInteractingGeometry2AABB()]),
+ PeriodicInsertionSortCollider(),
+ InteractionDispatchers(
+ [ef2_Sphere_Sphere_Dem3DofGeom()],
+ [SimpleElasticRelationships()],
+ [Law2_Dem3Dof_Elastic_Elastic()],
+ ),
+ PeriIsoCompressor(charLen=.1,stresses=[50e9,1e8],doneHook="print 'FINISHED'; O.pause() "),
+ NewtonsDampedLaw(damping=.4)
+]
+O.dt=utils.PWaveTimeStep()
+O.saveTmp()
+from yade import qt; qt.Controller(); qt.View()
+O.run()
+O.wait()
+timing.stats()
Modified: trunk/scripts/test/periodic-grow.py
===================================================================
--- trunk/scripts/test/periodic-grow.py 2009-08-09 17:17:53 UTC (rev 1932)
+++ trunk/scripts/test/periodic-grow.py 2009-08-10 11:32:23 UTC (rev 1933)
@@ -32,7 +32,7 @@
step=(mx-mn); step=Vector3(.002*step[0],.002*step[1],.002*step[2])
O.periodicCell=mn+step,mx-step
if (i%10==0):
- F=utils.totalForceInVolume()
+ F,stiff=utils.totalForceInVolume()
dim=mx-mn; A=Vector3(dim[1]*dim[2],dim[0]*dim[2],dim[0]*dim[1])
avgStress=sum([F[i]/A[i] for i in 0,1,2])/3.
print 'strain',(cubeSize-dim[0])/cubeSize,'avg. stress ',avgStress,'unbalanced ',utils.unbalancedForce()
