yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02157
[Branch ~yade-dev/yade/trunk] Rev 1771: 1. Enhance VTKRecorder to save more CPM data
------------------------------------------------------------
revno: 1771
committer: Václav Šmilauer <vaclav@flux>
branch nick: trunk
timestamp: Sun 2009-10-11 13:22:25 +0200
message:
1. Enhance VTKRecorder to save more CPM data
2. Add confinement storage to CpmMat; CpmStateUpdater now computes that.
3. Add method to get aabb of Gauss-smoothed domain.
modified:
lib/smoothing/WeightedAverage2d.hpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp
pkg/dem/PreProcessor/UniaxialStrainerGen.cpp
pkg/dem/meta/ConcretePM.cpp
pkg/dem/meta/ConcretePM.hpp
py/WeightedAverage2d.cpp
py/_eudoxos.cpp
py/_utils.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'lib/smoothing/WeightedAverage2d.hpp'
--- lib/smoothing/WeightedAverage2d.hpp 2009-10-04 14:37:41 +0000
+++ lib/smoothing/WeightedAverage2d.hpp 2009-10-11 11:22:25 +0000
@@ -32,6 +32,8 @@
Vector2r cellSizes;
Vector2i nCells;
public:
+ Vector2r getLo(){return lo;}
+ Vector2r getHi(){return hi;}
vector<vector<vector<T> > > grid;
/* construct grid from lower left corner, upper right corner and number of cells in each direction */
GridContainer(Vector2r _lo, Vector2r _hi, Vector2i _nCells): lo(_lo), hi(_hi), nCells(_nCells){
@@ -155,13 +157,13 @@
//struct Scalar2d{Vector2r pos; Real val;};
Vector2r tuple2vec2r(const python::tuple& t){return Vector2r(python::extract<Real>(t[0])(),python::extract<Real>(t[1])());}
Vector2i tuple2vec2i(const python::tuple& t){return Vector2i(python::extract<int>(t[0])(),python::extract<int>(t[1])());}
- shared_ptr<SGDA_Scalar2d> sgka;
+ shared_ptr<SGDA_Scalar2d> sgda;
struct Poly2d{vector<Vector2r> vertices; bool inclusive;};
vector<Poly2d> clips;
public:
pyGaussAverage(python::tuple lo, python::tuple hi, python::tuple nCells, Real stDev){
shared_ptr<GridContainer<Scalar2d> > g(new GridContainer<Scalar2d>(tuple2vec2r(lo),tuple2vec2r(hi),tuple2vec2i(nCells)));
- sgka=shared_ptr<SGDA_Scalar2d>(new SGDA_Scalar2d(g,stDev));
+ sgda=shared_ptr<SGDA_Scalar2d>(new SGDA_Scalar2d(g,stDev));
}
bool pointInsidePolygon(const Vector2r&,const vector<Vector2r>&);
bool ptIsClipped(const Vector2r& pt){
@@ -171,10 +173,11 @@
}
return false;
}
- bool addPt(Real val, python::tuple pos){Scalar2d d; d.pos=tuple2vec2r(pos); if(ptIsClipped(d.pos)) return false; d.val=val; sgka->grid->add(d,d.pos); return true; }
- Real avg(python::tuple _pt){Vector2r pt=tuple2vec2r(_pt); if(ptIsClipped(pt)) return std::numeric_limits<Real>::quiet_NaN(); return sgka->computeAverage(pt);}
- Real stDev_get(){return sgka->stDev;} void stDev_set(Real s){sgka->stDev=s;}
- Real relThreshold_get(){return sgka->relThreshold;} void relThreshold_set(Real rt){sgka->relThreshold=rt;}
+ bool addPt(Real val, python::tuple pos){Scalar2d d; d.pos=tuple2vec2r(pos); if(ptIsClipped(d.pos)) return false; d.val=val; sgda->grid->add(d,d.pos); return true; }
+ Real avg(python::tuple _pt){Vector2r pt=tuple2vec2r(_pt); if(ptIsClipped(pt)) return std::numeric_limits<Real>::quiet_NaN(); return sgda->computeAverage(pt);}
+ Real stDev_get(){return sgda->stDev;} void stDev_set(Real s){sgda->stDev=s;}
+ Real relThreshold_get(){return sgda->relThreshold;} void relThreshold_set(Real rt){sgda->relThreshold=rt;}
+ python::tuple aabb_get(){return python::make_tuple(sgda->grid->getLo(),sgda->grid->getHi());}
python::list clips_get(){
python::list ret;
FOREACH(const Poly2d& poly, clips){
@@ -199,10 +202,10 @@
}
python::tuple data_get(){
python::list x,y,val;
- const Vector2i& dim=sgka->grid->getSize();
+ const Vector2i& dim=sgda->grid->getSize();
for(int i=0; i<dim[0]; i++){
for(int j=0; j<dim[1]; j++){
- FOREACH(const Scalar2d& element, sgka->grid->grid[i][j]){
+ FOREACH(const Scalar2d& element, sgda->grid->grid[i][j]){
x.append(element.pos[0]); y.append(element.pos[1]); val.append(element.val);
}
}
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-10-04 14:37:41 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-10-11 11:22:25 +0000
@@ -44,7 +44,7 @@
if(rec=="spheres") recActive[REC_SPHERES]=true;
else if(rec=="facets") recActive[REC_FACETS]=true;
else if(rec=="colors") recActive[REC_COLORS]=true;
- else if(rec=="cpmDamage") recActive[REC_CPM_DAMAGE]=true;
+ else if(rec=="cpm") recActive[REC_CPM]=true;
else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: spheres, facets, colors, cpmDamage). Ignored.");
}
@@ -56,9 +56,14 @@
vtkSmartPointer<vtkFloatArray> spheresColors = vtkSmartPointer<vtkFloatArray>::New();
spheresColors->SetNumberOfComponents(3);
spheresColors->SetName("Colors");
- vtkSmartPointer<vtkFloatArray> damage = vtkSmartPointer<vtkFloatArray>::New();
- damage->SetNumberOfComponents(1);
- damage->SetName("cpmDamage");
+
+ if(recActive[REC_CPM]) CpmStateUpdater::update(rootBody);
+ vtkSmartPointer<vtkFloatArray> cpmDamage = vtkSmartPointer<vtkFloatArray>::New();
+ cpmDamage->SetNumberOfComponents(1);
+ cpmDamage->SetName("cpmDamage");
+ vtkSmartPointer<vtkFloatArray> cpmStress = vtkSmartPointer<vtkFloatArray>::New();
+ cpmStress->SetNumberOfComponents(1);
+ cpmStress->SetName("cpmStress");
vtkSmartPointer<vtkPoints> facetsPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> facetsCells = vtkSmartPointer<vtkCellArray>::New();
@@ -84,8 +89,9 @@
float c[3] = {color[0],color[1],color[2]};
spheresColors->InsertNextTupleValue(c);
}
- if (recActive[REC_CPM_DAMAGE]) {
- damage->InsertNextValue(YADE_PTR_CAST<CpmMat>(b->physicalParameters)->normDmg);
+ if (recActive[REC_CPM]) {
+ cpmDamage->InsertNextValue(YADE_PTR_CAST<CpmMat>(b->physicalParameters)->normDmg);
+ cpmStress->InsertNextValue(YADE_PTR_CAST<CpmMat>(b->physicalParameters)->avgStress);
}
continue;
}
@@ -127,7 +133,10 @@
spheresUg->SetCells(VTK_VERTEX, spheresCells);
spheresUg->GetPointData()->AddArray(radii);
if (recActive[REC_COLORS]) spheresUg->GetPointData()->AddArray(spheresColors);
- if (recActive[REC_CPM_DAMAGE]) spheresUg->GetPointData()->AddArray(damage);
+ if (recActive[REC_CPM]) {
+ spheresUg->GetPointData()->AddArray(cpmDamage);
+ spheresUg->GetPointData()->AddArray(cpmStress);
+ }
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
if(compress) writer->SetCompressor(compressor);
string fn=fileName+"spheres."+lexical_cast<string>(rootBody->currentIteration)+".vtu";
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-10-04 14:37:41 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-10-11 11:22:25 +0000
@@ -4,7 +4,7 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM_DAMAGE,REC_SENTINEL};
+ enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM,REC_SENTINEL};
//! A stuff to record: spheres,facets,colors
vector<string> recorders;
string fileName;
=== modified file 'pkg/dem/PreProcessor/UniaxialStrainerGen.cpp'
--- pkg/dem/PreProcessor/UniaxialStrainerGen.cpp 2009-08-05 07:32:18 +0000
+++ pkg/dem/PreProcessor/UniaxialStrainerGen.cpp 2009-10-11 11:22:25 +0000
@@ -117,7 +117,7 @@
newton->damping=damping;
rootBody->engines.push_back(newton);
- rootBody->engines.push_back(shared_ptr<CpmPhysDamageColorizer>(new CpmPhysDamageColorizer));
+ rootBody->engines.push_back(shared_ptr<CpmStateUpdater>(new CpmStateUpdater));
}
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2009-08-03 10:02:11 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2009-10-11 11:22:25 +0000
@@ -10,7 +10,7 @@
#ifdef YADE_OPENGL
(GLDrawCpmPhys)
#endif
- (CpmPhysDamageColorizer));
+ (CpmStateUpdater));
/********************** Ip2_CpmMat_CpmMat_CpmPhys ****************************/
@@ -135,6 +135,7 @@
#include"../../../../brefcom-mm.hh"
#endif
+
void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
@@ -289,7 +290,8 @@
#endif
/********************** CpmGlobalCharacteristics ****************************/
-
+/*** DEPRECATED ***/
+#if 0
CREATE_LOGGER(CpmGlobalCharacteristics);
void CpmGlobalCharacteristics::compute(MetaBody* rb, bool useMaxForce){
rb->bex.sync();
@@ -336,28 +338,42 @@
}
#endif
}
-
-
-/********************** CpmPhysDamageColorizer ****************************/
-CREATE_LOGGER(CpmPhysDamageColorizer);
-void CpmPhysDamageColorizer::action(MetaBody* rootBody){
- //vector<pair<short,Real> > bodyDamage; /* number of cohesive interactions per body; cummulative damage of interactions */
- //vector<pair<short,
+#endif
+
+/********************** CpmStateUpdater ****************************/
+CREATE_LOGGER(CpmStateUpdater);
+Real CpmStateUpdater::maxOmega=0.;
+
+void CpmStateUpdater::update(MetaBody* _rootBody){
+ MetaBody *rootBody=_rootBody?_rootBody:Omega::instance().getRootBody().get();
vector<BodyStats> bodyStats; bodyStats.resize(rootBody->bodies->size());
assert(bodyStats[0].nCohLinks==0); // should be initialized by dfault ctor
FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
- shared_ptr<CpmPhys> BC=dynamic_pointer_cast<CpmPhys>(I->interactionPhysics);
- if(!BC || !BC->isCohesive) continue;
+ if(!I->isReal()) continue;
+ shared_ptr<CpmPhys> phys=dynamic_pointer_cast<CpmPhys>(I->interactionPhysics);
+ if(!phys) continue;
const body_id_t id1=I->getId1(), id2=I->getId2();
- bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum+=(1-BC->relResidualStrength); bodyStats[id1].epsPlSum+=BC->epsPlSum;
- bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum+=(1-BC->relResidualStrength); bodyStats[id2].epsPlSum+=BC->epsPlSum;
- maxOmega=max(maxOmega,BC->omega);
+ Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
+
+ Vector3r f[]={phys->normalForce,-phys->normalForce}; const Vector3r pos[]={geom->se31.position,geom->se32.position}; const body_id_t ids[]={id1,id2};
+ Real stress=phys->normalForce.Length()/phys->crossSection;
+ for(int i=0; i<2; i++){
+ int sgn=(geom->contactPoint-pos[i]).Dot(f[i])>0?1:-1;
+ bodyStats[ids[i]].avgStress+=sgn*stress;
+ bodyStats[ids[i]].nLinks++;
+ }
+
+ if(!phys->isCohesive) continue;
+ bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum+=(1-phys->relResidualStrength); bodyStats[id1].epsPlSum+=phys->epsPlSum;
+ bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum+=(1-phys->relResidualStrength); bodyStats[id2].epsPlSum+=phys->epsPlSum;
+ maxOmega=max(maxOmega,phys->omega);
}
FOREACH(shared_ptr<Body> B, *rootBody->bodies){
const body_id_t& id=B->getId();
// add damaged contacts that have already been deleted
CpmMat* bpp=dynamic_cast<CpmMat*>(B->physicalParameters.get());
if(!bpp) continue;
+ bpp->avgStress=bodyStats[id].avgStress/bodyStats[id].nLinks;
int cohLinksWhenever=bodyStats[id].nCohLinks+bpp->numBrokenCohesive;
if(cohLinksWhenever>0){
bpp->normDmg=(bodyStats[id].dmgSum+bpp->numBrokenCohesive)/cohLinksWhenever;
=== modified file 'pkg/dem/meta/ConcretePM.hpp'
--- pkg/dem/meta/ConcretePM.hpp 2009-08-13 08:01:41 +0000
+++ pkg/dem/meta/ConcretePM.hpp 2009-10-11 11:22:25 +0000
@@ -37,7 +37,7 @@
* GLDrawCpmPhys draws interaction physics (color for damage and a few other); rarely used, though.
- * CpmPhysDamageColorizer changes bodies' colors depending on average damage of their interactions
+ * CpmStateUpdater changes bodies' colors depending on average damage of their interactions
and number of interactions that were already fully broken and have disappeared. This engine
contains its own loop (2 loops, more precisely) over all bodies and is run periodically
to update colors.
@@ -69,8 +69,10 @@
Real epsPlBroken;
//! sum of plastic strains normalized by number of contacts
Real normEpsPl;
- CpmMat(): epsVolumetric(0.), numBrokenCohesive(0), numContacts(0), normDmg(0.), epsPlBroken(0.), normEpsPl(0.) {createIndex();};
- REGISTER_ATTRIBUTES(BodyMacroParameters, (epsVolumetric) (numBrokenCohesive) (numContacts) (normDmg) (epsPlBroken) (normEpsPl));
+ //! averaged stress on the particle
+ Real avgStress;
+ CpmMat(): epsVolumetric(0.), numBrokenCohesive(0), numContacts(0), normDmg(0.), epsPlBroken(0.), normEpsPl(0.), avgStress(0.) {createIndex();};
+ REGISTER_ATTRIBUTES(BodyMacroParameters, (epsVolumetric) (numBrokenCohesive) (numContacts) (normDmg) (epsPlBroken) (normEpsPl) (avgStress));
REGISTER_CLASS_AND_BASE(CpmMat,BodyMacroParameters);
REGISTER_CLASS_INDEX(CpmMat,BodyMacroParameters);
};
@@ -263,6 +265,10 @@
static Real relKnSoft;
Law2_Dem3DofGeom_CpmPhys_Cpm() { }
void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
+ // utility functions
+ //! Update avgStress on all bodies (called from VTKRecorder and yade.eudoxos.particleConfinement)
+ //! Might be anywhere else as well (static method)
+ static void updateBodiesState(MetaBody*);
FUNCTOR2D(Dem3DofGeom,CpmPhys);
REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_CpmPhys_Cpm,ConstitutiveLaw);
REGISTER_ATTRIBUTES(ConstitutiveLaw,(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(omegaThreshold)(epsSoft)(relKnSoft));
@@ -270,6 +276,8 @@
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CpmPhys_Cpm);
+//deprecated code
+#if 0
/* Engine encompassing several computations looping over all bodies/interactions
*
* * Compute and store unbalanced force over the whole simulation.
@@ -292,6 +300,7 @@
REGISTER_CLASS_AND_BASE(CpmGlobalCharacteristics,PeriodicEngine);
};
REGISTER_SERIALIZABLE(CpmGlobalCharacteristics);
+#endif
#ifdef YADE_OPENGL
#include<yade/pkg-common/GLDrawFunctors.hpp>
@@ -307,16 +316,17 @@
REGISTER_SERIALIZABLE(GLDrawCpmPhys);
#endif
-class CpmPhysDamageColorizer: public PeriodicEngine {
- struct BodyStats{ int nCohLinks; Real dmgSum; Real epsPlSum; BodyStats(): nCohLinks(0), dmgSum(0.), epsPlSum(0.){} };
+class CpmStateUpdater: public PeriodicEngine {
+ struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum, avgStress; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), avgStress(0.) {} };
public:
//! maximum damage over all contacts
- Real maxOmega;
- CpmPhysDamageColorizer(){maxOmega=0; /* run at the very beginning */ initRun=true;}
- virtual void action(MetaBody*);
+ static Real maxOmega;
+ CpmStateUpdater(){maxOmega=0; /* run at the very beginning */ initRun=true;}
+ virtual void action(MetaBody* rb){ update(rb); }
+ static void update(MetaBody* rb=NULL);
DECLARE_LOGGER;
REGISTER_ATTRIBUTES(PeriodicEngine,(maxOmega));
- REGISTER_CLASS_AND_BASE(CpmPhysDamageColorizer,PeriodicEngine);
+ REGISTER_CLASS_AND_BASE(CpmStateUpdater,PeriodicEngine);
};
-REGISTER_SERIALIZABLE(CpmPhysDamageColorizer);
+REGISTER_SERIALIZABLE(CpmStateUpdater);
=== modified file 'py/WeightedAverage2d.cpp'
--- py/WeightedAverage2d.cpp 2009-07-27 17:08:23 +0000
+++ py/WeightedAverage2d.cpp 2009-10-11 11:22:25 +0000
@@ -23,6 +23,7 @@
.add_property("relThreshold",&pyGaussAverage::relThreshold_get,&pyGaussAverage::relThreshold_set)
.add_property("clips",&pyGaussAverage::clips_get,&pyGaussAverage::clips_set)
.add_property("data",&pyGaussAverage::data_get)
+ .add_property("aabb",&pyGaussAverage::aabb_get)
;
};
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2009-10-04 14:37:41 +0000
+++ py/_eudoxos.cpp 2009-10-11 11:22:25 +0000
@@ -125,21 +125,8 @@
return ret;
}
-std::vector<Real> particleConfinement(){
- MetaBody* rb=Omega::instance().getRootBody().get();
- vector<Real> ret(rb->bodies->size(),0.);
- FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
- if(!I->isReal()) continue;
- CpmPhys* phys=YADE_CAST<CpmPhys*>(I->interactionPhysics.get());
- Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
- Vector3r f[]={phys->normalForce,-phys->normalForce}; Vector3r pos[]={geom->se31.position,geom->se32.position}; body_id_t ids[]={I->getId1(),I->getId2()};
- Real stress=phys->normalForce.Length()/phys->crossSection;
- for(int i=0; i<2; i++){
- int sgn=(geom->contactPoint-pos[i]).Dot(f[i])>0?1:-1;
- ret[ids[i]]+=sgn*stress;
- }
- }
- return ret;
+void particleConfinement(){
+ CpmStateUpdater::update();
}
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2009-10-04 14:37:41 +0000
+++ py/_utils.cpp 2009-10-11 11:22:25 +0000
@@ -234,13 +234,19 @@
return ret;
}
-/* Sum force acting on facets given by their ids in the sense of their respective normals. */
-Real sumFacetNormalForces(vector<body_id_t> ids){
+/* Sum force acting on facets given by their ids in the sense of their respective normals.
+ If axis is given, it will sum forces perpendicular to given axis only (not the the facet normals).
+*/
+Real sumFacetNormalForces(vector<body_id_t> ids, int axis=-1){
shared_ptr<MetaBody> rb=Omega::instance().getRootBody(); rb->bex.sync();
Real ret=0;
FOREACH(const body_id_t id, ids){
InteractingFacet* f=YADE_CAST<InteractingFacet*>(Body::byId(id,rb)->interactingGeometry.get());
- ret+=rb->bex.getForce(id).Dot(f->nf);
+ if(axis<0) ret+=rb->bex.getForce(id).Dot(f->nf);
+ else {
+ Vector3r ff=rb->bex.getForce(id); ff[axis]=0;
+ ret+=ff.Dot(f->nf);
+ }
}
return ret;
}
@@ -402,7 +408,7 @@
def("kineticEnergy",Shop__kineticEnergy);
def("sumBexForces",sumBexForces);
def("sumBexTorques",sumBexTorques);
- def("sumFacetNormalForces",sumFacetNormalForces);
+ def("sumFacetNormalForces",sumFacetNormalForces,(python::arg("axis")=-1));
def("forcesOnPlane",forcesOnPlane);
def("forcesOnCoordPlane",forcesOnCoordPlane);
def("totalForceInVolume",Shop__totalForceInVolume,"Return summed forces on all interactions and average isotropic stiffness, as tuple (Vector3,float)");