Author: eudoxos
Date: 2008-09-09 12:07:09 +0200 (Tue, 09 Sep 2008)
New Revision: 1513
Modified:
trunk/core/yade.cpp
trunk/extra/Brefcom.hpp
trunk/gui/SConscript
trunk/gui/py/_utils.cpp
trunk/gui/py/eudoxos.py
trunk/gui/py/utils.py
trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
trunk/pkg/common/DataClass/PhysicalAction/Force.cpp
trunk/pkg/common/DataClass/PhysicalAction/Force.hpp
trunk/pkg/common/DataClass/PhysicalAction/Momentum.cpp
trunk/pkg/common/DataClass/PhysicalAction/Momentum.hpp
Log:
1. add utility function to compute elastic energy within a volume (the
dynamic_cast to NormalShearInteraction still fails; what is going on?)
2. add utils.fractionalBox, an AABB reduced to its fraction
3. fix fixmes in Force and Momentum about unused functions about to be deleted
4. plotting interaction direction histogram now works on contrained volume (to
study boundary influence on the distribution)
Modified: trunk/core/yade.cpp
===================================================================
--- trunk/core/yade.cpp 2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/core/yade.cpp 2008-09-09 10:07:09 UTC (rev 1513)
@@ -171,7 +171,7 @@
std::string logConf=configPath+"/logging.conf";
if(filesystem::exists(logConf)){
log4cxx::PropertyConfigurator::configure(logConf);
- LOG_INFO("Loading "<<logConf);
+ LOG_INFO("Loading "<<logConf);
} else { // otherwise use simple console-directed logging
log4cxx::BasicConfigurator::configure();
logger->setLevel(log4cxx::Level::WARN);
Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp 2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/extra/Brefcom.hpp 2008-09-09 10:07:09 UTC (rev 1513)
@@ -95,9 +95,10 @@
BrefcomContact(): NormalShearInteraction(),E(0), G(0),
tanFrictionAngle(0), undamagedCohesion(0), equilibriumDist(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; }
// BrefcomContact(Real _E, Real _G, Real
_tanFrictionAngle, Real _undamagedCohesion, Real _equilibriumDist, Real
_crossSection, Real _epsCrackOnset, Real _epsFracture, Real _expBending, Real
_xiShear, Real _tau=0, Real _expDmgRate=1): InteractionPhysics(), E(_E), G(_G),
tanFrictionAngle(_tanFrictionAngle), undamagedCohesion(_undamagedCohesion),
equilibriumDist(_equilibriumDist), crossSection(_crossSection),
epsCrackOnset(_epsCrackOnset), epsFracture(_epsFracture),
expBending(_expBending), xiShear(_xiShear), tau(_tau), expDmgRate(_expDmgRate)
{ epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0;
Fn=0; Fs=Vector3r::ZERO;
/*TRVAR5(epsCrackOnset,epsFracture,Kn,crossSection,equilibriumDist); */ }
+ virtual ~BrefcomContact(){}
- void registerAttributes(){
+ virtual void registerAttributes(){
InteractionPhysics::registerAttributes();
REGISTER_ATTRIBUTE(E);
REGISTER_ATTRIBUTE(G);
Modified: trunk/gui/SConscript
===================================================================
--- trunk/gui/SConscript 2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/gui/SConscript 2008-09-09 10:07:09 UTC (rev 1513)
@@ -65,7 +65,7 @@
],
),
env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',
- LIBS=env['LIBS']+['Shop','libboost_python']),
+ LIBS=env['LIBS']+['Shop','libboost_python','Brefcom']),
env.SharedLibrary('log',['py/log.cpp'],SHLIBPREFIX=''),
env.File('__init__.py','py'),
env.File('utils.py','py'),
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/gui/py/_utils.cpp 2008-09-09 10:07:09 UTC (rev 1513)
@@ -4,9 +4,14 @@
#include<yade/core/Omega.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-dem/SpheresContactGeometry.hpp>
+#include<yade/extra/Brefcom.hpp>
#include<cmath>
using namespace boost::python;
+#ifdef LOG4CXX
+ log4cxx::LoggerPtr logger=log4cxx::Logger::getLogger("yade.utils");
+#endif
+
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(extract<double>(t[0])(),extract<double>(t[1])(),extract<double>(t[2])());}
bool isInBB(Vector3r p, Vector3r bbMin, Vector3r bbMax){return p[0]>bbMin[0]
&& p[0]<bbMax[0] && p[1]>bbMin[1] && p[1]<bbMax[1] && p[2]>bbMin[2] &&
p[2]<bbMax[2];}
@@ -59,7 +64,37 @@
Real PWaveTimeStep(){return Shop::PWaveTimeStep();};
+Real elasticEnergyInAABB(python::tuple AABB){
+ Vector3r bbMin=tuple2vec(extract<python::tuple>(AABB[0])()),
bbMax=tuple2vec(extract<python::tuple>(AABB[1])());
+ shared_ptr<MetaBody> rb=Omega::instance().getRootBody();
+ Real E=0;
+ FOREACH(const shared_ptr<Interaction>&i, *rb->transientInteractions){
+ if(!i->interactionPhysics) continue;
+ //FIXME: dynamic_pointer_cast instead of static should be used
here, but the downcast fails for reasons that are not know :-(
+ shared_ptr<BrefcomContact>
bc=static_pointer_cast<BrefcomContact>(i->interactionPhysics); if(!bc) continue;
+ const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb),
b2=Body::byId(i->getId2(),rb);
+ bool
isIn1=isInBB(b1->physicalParameters->se3.position,bbMin,bbMax),
isIn2=isInBB(b2->physicalParameters->se3.position,bbMin,bbMax);
+ if(!isIn1 && !isIn2) continue;
+ LOG_DEBUG("Interaction #"<<i->getId1()<<"--#"<<i->getId2());
+ //cerr<<"Interaction #"<<i->getId1()<<"--#"<<i->getId2()<<endl;
+ Real weight=1.;
+ if((!isIn1 && isIn2) || (isIn1 && !isIn2)){
+ //shared_ptr<Body> bIn=isIn1?b1:b2, bOut=isIn2?b2:b1;
+ Vector3r
vIn=(isIn1?b1:b2)->physicalParameters->se3.position,
vOut=(isIn2?b1:b2)->physicalParameters->se3.position;
+ #define _WEIGHT_COMPONENT(axis)
if(vOut[axis]<bbMin[axis])
weight=min(weight,abs((vOut[axis]-bbMin[axis])/(vOut[axis]-vIn[axis]))); else
if(vOut[axis]>bbMax[axis])
weight=min(weight,abs((vOut[axis]-bbMax[axis])/(vOut[axis]-vIn[axis])));
+ _WEIGHT_COMPONENT(0); _WEIGHT_COMPONENT(1);
_WEIGHT_COMPONENT(2);
+ assert(weight>=0 && weight<=1);
+ //cerr<<"Interaction crosses AABB boundary, weight is
"<<weight<<endl;
+ //LOG_DEBUG("Interaction crosses AABB boundary, weight
is "<<weight);
+ } else {assert(isIn1 && isIn2); /* cerr<<"Interaction inside,
weight is "<<weight<<endl;*/ /*LOG_DEBUG("Interaction inside, weight is
"<<weight);*/}
+
E+=weight*(.5*bc->E*bc->crossSection*pow(bc->epsN,2)+.5*bc->G*bc->crossSection*pow(bc->epsT.Length(),2));
+
//TRVAR5(bc->crossSection,bc->E,bc->G,bc->epsN,bc->epsT.Length());;
+
//cerr<<bc->crossSection<<","<<bc->E<<","<<bc->G<<","<<bc->epsN<<","<<bc->epsT.Length()<<endl;
+ }
+ return E;
+}
+
/* return histogram ([bin1Min,bin2Min,…],[value1,value2,…]) from projections
of interactions
* to the plane perpendicular to axis∈[0,1,2]; the number of bins can be
specified and they cover
* the range (0,π), since interactions are bidirecional, hence periodically
the same on (π,2π).
@@ -68,16 +103,20 @@
* Both bodies must be in the mask (except if mask is 0, when all bodies are
considered)
* If the projection is shorter than minProjLen, it is skipped.
*
- * @todo implement groupMask
+ * If both bodies are _outside_ the aabb (if specified), the interaction is
skipped.
+ *
*/
-python::tuple interactionAnglesHistogram(int axis, int mask=0, size_t bins=20,
Real minProjLen=1e-6){
+python::tuple interactionAnglesHistogram(int axis, int mask=0, size_t bins=20,
python::tuple aabb=python::tuple(), Real minProjLen=1e-6){
if(axis<0||axis>2) throw invalid_argument("Axis must be from
{0,1,2}=x,y,z.");
+ Vector3r bbMin(Vector3r::ZERO), bbMax(Vector3r::ZERO); bool
useBB=python::len(aabb)>0;
if(useBB){bbMin=tuple2vec(extract<python::tuple>(aabb[0])());bbMax=tuple2vec(extract<python::tuple>(aabb[1])());}
Real binStep=Mathr::PI/bins; int axis2=(axis+1)%3, axis3=(axis+2)%3;
vector<Real> cummProj(bins,0.);
shared_ptr<MetaBody> rb=Omega::instance().getRootBody();
FOREACH(const shared_ptr<Interaction>& i, *rb->transientInteractions){
if(!i->isReal) continue;
- if(!(mask==0 || (mask>0 &&
(Body::byId(i->getId1(),rb)->getGroupMask() & mask) &&
(Body::byId(i->getId2(),rb)->getGroupMask() & mask)))) continue;
+ const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb),
b2=Body::byId(i->getId2(),rb);
+ if(!b1->maskOk(mask) || !b2->maskOk(mask)) continue;
+ if(useBB &&
!isInBB(b1->physicalParameters->se3.position,bbMin,bbMax) &&
!isInBB(b2->physicalParameters->se3.position,bbMin,bbMax)) continue;
shared_ptr<SpheresContactGeometry>
scg=dynamic_pointer_cast<SpheresContactGeometry>(i->interactionGeometry);
if(!scg) continue;
Vector3r n(scg->normal); n[axis]=0.; Real nLen=n.Length();
if(nLen<minProjLen) continue; // this interaction is (almost)
exactly parallel to our axis; skip that one
@@ -89,7 +128,7 @@
for(size_t i=0; i<(size_t)bins; i++){ val.append(cummProj[i]);
binMid.append(i*binStep);}
return python::make_tuple(binMid,val);
}
-BOOST_PYTHON_FUNCTION_OVERLOADS(interactionAnglesHistogram_overloads,interactionAnglesHistogram,1,3);
+BOOST_PYTHON_FUNCTION_OVERLOADS(interactionAnglesHistogram_overloads,interactionAnglesHistogram,1,4);
BOOST_PYTHON_MODULE(_utils){
@@ -98,7 +137,8 @@
def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor")));
def("coordsAndDisplacements",coordsAndDisplacements,coordsAndDisplacements_overloads(args("AABB")));
def("setRefSe3",setRefSe3);
-
def("interactionAnglesHistogram",interactionAnglesHistogram,interactionAnglesHistogram_overloads(args("axis","mask","bins")));
+
def("interactionAnglesHistogram",interactionAnglesHistogram,interactionAnglesHistogram_overloads(args("axis","mask","bins","aabb")));
+ def("elasticEnergy",elasticEnergyInAABB);
}
Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py 2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/gui/py/eudoxos.py 2008-09-09 10:07:09 UTC (rev 1513)
@@ -5,15 +5,15 @@
#
-def plotDirections(mask=0,bins=20):
+def plotDirections(mask=0,bins=20,aabb=()):
"Plot 3 histograms for distribution of interaction directions, in yz,xz
and xy planes."
- import pylab
+ import pylab,math
from yade import utils
for axis in [0,1,2]:
- d=utils.interactionAnglesHistogram(axis,mask=mask,bins=bins)
+
d=utils.interactionAnglesHistogram(axis,mask=mask,bins=bins,aabb=aabb)
fc=[0,0,0]; fc[axis]=1.
pylab.subplot(220+axis+1,polar=True);
-
bar(d[0],d[1],width=math.pi/(1.2*bins),fc=fc,alpha=.7,label=['yz','xz','xy'][axis])
+
pylab.bar(d[0],d[1],width=math.pi/(1.2*bins),fc=fc,alpha=.7,label=['yz','xz','xy'][axis])
pylab.show()
def estimatePoissonYoung(principalAxis,stress=0,plot=False,cutoff=0.):
@@ -32,9 +32,7 @@
dd=[] # storage for linear regression parameters
import pylab,numpy,stats
from yade import utils
- if cutoff>0:
- aabb=utils.aabbExtrema(); half=[.5*(aabb[1][i]-aabb[0][i]) for
i in [0,1,2]]
- cut=(tuple([aabb[0][i]+cutoff*half[i] for i in
[0,1,2]]),tuple([aabb[1][i]-cutoff*half[i] for i in [0,1,2]]))
+ if cutoff>0: cut=utils.fractionalBox(fraction=1-cutoff)
for axis in [0,1,2]:
if cutoff>0:
#print cutoff,utils.aabbExtrema(),cut
Modified: trunk/gui/py/utils.py
===================================================================
--- trunk/gui/py/utils.py 2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/gui/py/utils.py 2008-09-09 10:07:09 UTC (rev 1513)
@@ -72,6 +72,14 @@
other=((axis+1)%3,(axis+2)%3)
return
(ext[1][other[0]]-ext[0][other[0]])*(ext[1][other[1]]-ext[0][other[1]])
+def fractionalBox(fraction=1.,minMax=None):
+ """retrurn (min,max) that is the original minMax box (or aabb of the
whole simulation if not specified)
+ linearly scaled around its center to the fraction factor"""
+ if not minMax: minMax=aabbExtrema()
+ half=[.5*(minMax[1][i]-minMax[0][i]) for i in [0,1,2]]
+ return (tuple([minMax[0][i]+(1-fraction)*half[i] for i in
[0,1,2]]),tuple([minMax[1][i]-(1-fraction)*half[i] for i in [0,1,2]]))
+
+
def randomizeColors(onShapes=True,onMolds=False,onlyDynamic=False):
"""Assign random colors to shape's (GeometricalModel) and/or mold's
(InteractingGeometry) diffuseColor.
Modified:
trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
2008-09-08 08:56:42 UTC (rev 1512)
+++ trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
2008-09-09 10:07:09 UTC (rev 1513)
@@ -28,7 +28,7 @@
public:
//! shear stiffness
Real ks;
- NormalShearInteraction(){createIndex(); }
+ NormalShearInteraction(){ createIndex(); }
virtual ~NormalShearInteraction(){};
protected:
virtual void registerAttributes(){
Modified: trunk/pkg/common/DataClass/PhysicalAction/Force.cpp
===================================================================
--- trunk/pkg/common/DataClass/PhysicalAction/Force.cpp 2008-09-08 08:56:42 UTC
(rev 1512)
+++ trunk/pkg/common/DataClass/PhysicalAction/Force.cpp 2008-09-09 10:07:09 UTC
(rev 1513)
@@ -10,41 +10,12 @@
#include "Force.hpp"
+Force::Force() : PhysicalAction(){ createIndex(); }
-Force::Force() : PhysicalAction()
-{
- createIndex();
-}
+Force::~Force(){}
-Force::~Force()
-{
-}
+void Force::reset(){ force = Vector3r::ZERO; }
+shared_ptr<PhysicalAction> Force::clone(){ return shared_ptr<Force>(new
Force(*this));}
-/* FIXME - not used
-void Force::add(const shared_ptr<PhysicalAction>& a)
-{
- Force * f = static_cast<Force*>(a.get());
- force += f->force;
-}
-
-
-void Force::sub(const shared_ptr<PhysicalAction>& a)
-{
- Force * f = static_cast<Force*>(a.get());
- force -= f->force;
-}
-*/
-
-void Force::reset()
-{
- force = Vector3r::ZERO;
-}
-
-
-shared_ptr<PhysicalAction> Force::clone()
-{
- return shared_ptr<Force>(new Force(*this));
-}
-
YADE_PLUGIN();
Modified: trunk/pkg/common/DataClass/PhysicalAction/Force.hpp
===================================================================
--- trunk/pkg/common/DataClass/PhysicalAction/Force.hpp 2008-09-08 08:56:42 UTC
(rev 1512)
+++ trunk/pkg/common/DataClass/PhysicalAction/Force.hpp 2008-09-09 10:07:09 UTC
(rev 1513)
@@ -21,20 +21,14 @@
Force();
virtual ~Force();
-/// Methods
-// virtual void add(const shared_ptr<PhysicalAction>& a); // FIXME
- not used
-// virtual void sub(const shared_ptr<PhysicalAction>& a); // FIXME
- not used
-
virtual void reset();
virtual shared_ptr<PhysicalAction> clone();
-/// Serialization
REGISTER_CLASS_NAME(Force);
REGISTER_BASE_CLASS_NAME(PhysicalAction);
+ virtual void registerAttributes(){REGISTER_ATTRIBUTE(force);}
-/// Indexable
REGISTER_CLASS_INDEX(Force,PhysicalAction);
-
};
REGISTER_SERIALIZABLE(Force,false);
Modified: trunk/pkg/common/DataClass/PhysicalAction/Momentum.cpp
===================================================================
--- trunk/pkg/common/DataClass/PhysicalAction/Momentum.cpp 2008-09-08
08:56:42 UTC (rev 1512)
+++ trunk/pkg/common/DataClass/PhysicalAction/Momentum.cpp 2008-09-09
10:07:09 UTC (rev 1513)
@@ -10,40 +10,12 @@
#include "Momentum.hpp"
+Momentum::Momentum() : PhysicalAction(){createIndex();}
-Momentum::Momentum() : PhysicalAction()
-{
- createIndex();
-}
+Momentum::~Momentum(){}
-Momentum::~Momentum()
-{
-}
+void Momentum::reset(){ momentum = Vector3r::ZERO; }
-/* FIXME - not used
-void Momentum::add(const shared_ptr<PhysicalAction>& a)
-{
- Momentum * m = static_cast<Momentum*>(a.get());
- momentum += m->momentum;
-}
+shared_ptr<PhysicalAction> Momentum::clone(){return shared_ptr<Momentum>(new
Momentum(*this));}
-
-void Momentum::sub(const shared_ptr<PhysicalAction>& a)
-{
- Momentum * m = static_cast<Momentum*>(a.get());
- momentum -= m->momentum;
-}
-*/
-
-void Momentum::reset()
-{
- momentum = Vector3r::ZERO;
-}
-
-
-shared_ptr<PhysicalAction> Momentum::clone()
-{
- return shared_ptr<Momentum>(new Momentum(*this));
-}
-
YADE_PLUGIN();
Modified: trunk/pkg/common/DataClass/PhysicalAction/Momentum.hpp
===================================================================
--- trunk/pkg/common/DataClass/PhysicalAction/Momentum.hpp 2008-09-08
08:56:42 UTC (rev 1512)
+++ trunk/pkg/common/DataClass/PhysicalAction/Momentum.hpp 2008-09-09
10:07:09 UTC (rev 1513)
@@ -22,20 +22,15 @@
Momentum();
virtual ~Momentum();
-// virtual void add(const shared_ptr<PhysicalAction>& a); // FIXME
- not used
-// virtual void sub(const shared_ptr<PhysicalAction>& a); // FIXME
- not used
-
virtual void reset();
virtual shared_ptr<PhysicalAction> clone();
-/// Serialization
REGISTER_CLASS_NAME(Momentum);
REGISTER_BASE_CLASS_NAME(PhysicalAction);
+ virtual void registerAttributes(){REGISTER_ATTRIBUTE(momentum);}
-/// Indexable
REGISTER_CLASS_INDEX(Momentum,PhysicalAction);
};
-
REGISTER_SERIALIZABLE(Momentum,false);
#endif // MOMENTUM_HPP
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
