yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00627
[svn] r1513 - in trunk: core extra gui gui/py pkg/common/DataClass/InteractionPhysics pkg/common/DataClass/PhysicalAction
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
