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 <[email protected]>
+
+#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 <[email protected]>
+
+#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 <[email protected]>
-#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()
_______________________________________________
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
