yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04030
[Branch ~yade-dev/yade/trunk] Rev 2153: 1. Completely remove HydraulicForceEngine which didn't exist in sources (removed by Bruno in prev...
------------------------------------------------------------
revno: 2153
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Sun 2010-04-18 22:46:54 +0200
message:
1. Completely remove HydraulicForceEngine which didn't exist in sources (removed by Bruno in previous commit)
2. Experimental binary io using boost::serialization
3. Some function in yade.eudoxos
modified:
core/Omega.cpp
lib/factory/ClassFactory.hpp
pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp
py/_eudoxos.cpp
py/_utils.cpp
py/utils.py
py/yadeWrapper/yadeWrapper.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 'core/Omega.cpp'
--- core/Omega.cpp 2010-03-15 14:35:44 +0000
+++ core/Omega.cpp 2010-04-18 20:46:54 +0000
@@ -283,11 +283,6 @@
resetScene();
RenderMutexLock lock; IOFormatManager::loadFromFile("XMLFormatManager",simulationFileName,"scene",scene);
}
- else if(algorithm::ends_with(simulationFileName,".yade")){
- joinSimulationLoop();
- resetScene();
- RenderMutexLock lock; IOFormatManager::loadFromFile("BINFormatManager",simulationFileName,"scene",scene);
- }
else if(algorithm::starts_with(simulationFileName,":memory:")){
if(memSavedSimulations.count(simulationFileName)==0) throw runtime_error("Cannot load nonexistent memory-saved simulation "+simulationFileName);
istringstream iss(memSavedSimulations[simulationFileName]);
@@ -316,12 +311,6 @@
FormatChecker::format=FormatChecker::XML;
IOFormatManager::saveToFile("XMLFormatManager",name,"scene",scene);
}
- #if 0
- else if(algorithm::ends_with(name,".yade")){
- FormatChecker::format=FormatChecker::BIN;
- IOFormatManager::saveToFile("BINFormatManager",name,"scene",scene);
- }
- #endif
else if(algorithm::starts_with(name,":memory:")){
if(memSavedSimulations.count(simulationFileName)>0) LOG_INFO("Overwriting in-memory saved simulation "<<name);
ostringstream oss;
=== modified file 'lib/factory/ClassFactory.hpp'
--- lib/factory/ClassFactory.hpp 2010-02-09 16:50:30 +0000
+++ lib/factory/ClassFactory.hpp 2010-04-18 20:46:54 +0000
@@ -42,8 +42,8 @@
//#define YADE_BOOST_SERIALIZATION
#ifdef YADE_BOOST_SERIALIZATION
- //#include<boost/archive/binary_oarchive.hpp>
- //#include<boost/archive/binary_iarchive.hpp>
+ #include<boost/archive/binary_oarchive.hpp>
+ #include<boost/archive/binary_iarchive.hpp>
#include<boost/archive/xml_oarchive.hpp>
#include<boost/archive/xml_iarchive.hpp>
=== modified file 'pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp'
--- pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp 2010-03-27 22:18:10 +0000
+++ pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp 2010-04-18 20:46:54 +0000
@@ -26,7 +26,7 @@
#include<yade/pkg-common/BoundDispatcher.hpp>
#include<yade/pkg-common/GravityEngines.hpp>
-#include<yade/pkg-dem/HydraulicForceEngine.hpp>
+// #include<yade/pkg-dem/HydraulicForceEngine.hpp>
#include<yade/pkg-dem/NewtonIntegrator.hpp>
#include<yade/pkg-dem/Ig2_Sphere_Sphere_ScGeom.hpp>
#include<yade/pkg-dem/Ig2_Box_Sphere_ScGeom.hpp>
@@ -416,8 +416,8 @@
triaxialcompressionEngine->internalCompaction = internalCompaction;
triaxialcompressionEngine->maxMultiplier = maxMultiplier;
- shared_ptr<HydraulicForceEngine> hydraulicForceEngine = shared_ptr<HydraulicForceEngine> (new HydraulicForceEngine);
- hydraulicForceEngine->dummyParameter = true;
+ // shared_ptr<HydraulicForceEngine> hydraulicForceEngine = shared_ptr<HydraulicForceEngine> (new HydraulicForceEngine);
+ // hydraulicForceEngine->dummyParameter = true;
//cerr << "fin de section triaxialcompressionEngine = shared_ptr<TriaxialCompressionEngine> (new TriaxialCompressionEngine);" << std::endl;
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2010-04-01 06:53:09 +0000
+++ py/_eudoxos.cpp 2010-04-18 20:46:54 +0000
@@ -92,10 +92,124 @@
ret["vel"]=vel;
return ret;
}
+
+#if 0
+/* Compute stress tensor on each particle */
+void particleMacroStress(void){
+ Scene* scene=Omega::instance().getScene().get();
+ // find interactions of each particle
+ std::vector<std::list<body_id_t> > bIntr(scene->bodies->size());
+ FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
+ if(!i->isReal) continue;
+ // only contacts between 2 spheres
+ Sphere* s1=dynamic_cast<Sphere*>(Body::byId(i->getId1(),scene)->shape.get())
+ Sphere* s2=dynamic_cast<Sphere*>(Body::byId(i->getId2(),scene)->shape.get())
+ if(!s1 || !s2) continue;
+ bIntr[i.getId1()].push_back(i.getId2());
+ bIntr[i.getId2()].push_back(i.getId1());
+ }
+ for(body_id_t id1=0; id1<(body_id_t)bIntr.size(); id1++){
+ if(bIntr[id1].size()==0) continue;
+ Matrix3r ss(Matrix3r::ZERO); // stress tensor on current particle
+ FOREACH(body_id_t id2, bIntr[id1]){
+ const shared_ptr<Interaction> i=scene->interactions->find(id1,id2);
+ assert(i);
+ Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(i->interactionGeometry);
+ CpmPhys* phys=YADE_CAST<CpmPhys*>(i->interactionPhysics);
+ Real d=(geom->se31->pos-geom->se32->pos).Length(); // current contact length
+ const Vector3r& n=geom->normal;
+ const Real& A=phys->crossSection;
+ const Vector3r& sigmaT=phys->sigmaT;
+ const Real& sigmaN=phys->sigmaN;
+ for(int i=0; i<3; i++) for(int j=0;j<3; j++){
+ ss[i][j]=d*A*(sigmaN*n[i]*n[j]+.5*(sigmaT[i]*n[j]+sigmaT[j]*n[i]));
+ }
+ }
+ // divide ss by V of the particle
+ // FIXME: for now, just use 2*(4/3)*π*r³ (.5 porosity)
+ ss*=1/(2*(4/3)*Mathr::PI*);
+ }
+}
+#endif
+#include<yade/lib-smoothing/WeightedAverage2d.hpp>
+/* Fastly locate interactions within given distance from a point in 2d (projected to plane) */
+struct SpiralInteractionLocator2d{
+ struct FlatInteraction{ Real r,h,theta; shared_ptr<Interaction> i; FlatInteraction(Real _r, Real _h, Real _theta, const shared_ptr<Interaction>& _i): r(_r), h(_h), theta(_theta), i(_i){}; };
+ shared_ptr<GridContainer<FlatInteraction> > grid;
+ Real thetaSpan;
+ int axis;
+ SpiralInteractionLocator2d(Real dH_dTheta, int _axis, Real theta0): axis(_axis){
+ Scene* scene=Omega::instance().getScene().get();
+ Real inf=std::numeric_limits<Real>::infinity();
+ Vector2r lo=Vector2r(inf,inf), hi(-inf,-inf);
+ Real minD0(inf),maxD0(-inf);
+ Real minTheta(inf), maxTheta(-inf);
+ std::list<FlatInteraction> intrs;
+ // first pass: find extrema for positions and interaction lengths, build interaction list
+ FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
+ Dem3DofGeom* ge=dynamic_cast<Dem3DofGeom*>(i->interactionGeometry.get());
+ CpmPhys* ph=dynamic_cast<CpmPhys*>(i->interactionPhysics.get());
+ if(!ge || !ph) continue;
+ Real r,h,theta;
+ boost::tie(r,h,theta)=Shop::spiralProject(ge->contactPoint,dH_dTheta,axis,NaN,theta0);
+ lo=componentMinVector(lo,Vector2r(r,h)); hi=componentMaxVector(hi,Vector2r(r,h));
+ minD0=min(minD0,ge->refLength); maxD0=max(maxD0,ge->refLength);
+ minTheta=min(minTheta,theta); maxTheta=max(maxTheta,theta);
+ intrs.push_back(FlatInteraction(r,h,theta,i));
+ }
+ // create grid, put interactions inside
+ Vector2i nCells=Vector2i(max(10,(int)((hi[0]-lo[0])/(2*minD0))),max(10,(int)((hi[1]-lo[1])/(2*minD0))));
+ Vector2r hair=1e-2*Vector2r((hi[0]-lo[0])/nCells[0],(hi[1]-lo[1])/nCells[1]); // avoid rounding issue on boundary, enlarge by 1/100 cell size on each side
+ grid=shared_ptr<GridContainer<FlatInteraction> >(new GridContainer<FlatInteraction>(lo-hair,hi+hair,nCells));
+ FOREACH(const FlatInteraction& fi, intrs){
+ grid->add(fi,Vector2r(fi.r,fi.h));
+ }
+ thetaSpan=maxTheta-minTheta;
+ }
+ py::list intrsAroundPt(const Vector2r& pt, Real radius){
+ py::list ret;
+ FOREACH(const Vector2i& v, grid->circleFilter(pt,radius)){
+ FOREACH(const FlatInteraction& fi, grid->grid[v[0]][v[1]]){
+ if((pow(fi.r-pt[0],2)+pow(fi.h-pt[1],2))>radius*radius) continue; // too far
+ ret.append(fi.i);
+ }
+ }
+ return ret;
+ }
+ // return macroscopic stress around interactions that project around given point
+ // stresses are rotated around axis back by theta angle
+ Matrix3r macroStressAroundPt(const Vector2r& pt, Real radius){
+ Matrix3r ss(Matrix3r::ZERO);
+ FOREACH(const Vector2i& v, grid->circleFilter(pt,radius)){
+ FOREACH(const FlatInteraction& fi, grid->grid[v[0]][v[1]]){
+ if((pow(fi.r-pt[0],2)+pow(fi.h-pt[1],2))>radius*radius) continue; // too far
+ Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(fi.i->interactionGeometry.get());
+ CpmPhys* phys=YADE_CAST<CpmPhys*>(fi.i->interactionPhysics.get());
+ // transformation matrix, to rotate to the plane
+ Vector3r ax(Vector3r::ZERO); ax[axis]=1.;
+ Quaternionr q; q.FromAxisAngle(ax,fi.theta); q=q.Conjugate();
+ Matrix3r TT; q.ToRotationMatrix(TT);
+ //
+ Real d=(geom->se31.position-geom->se32.position).Length(); // current contact length
+ const Vector3r& n=TT*geom->normal;
+ const Real& A=phys->crossSection;
+ const Vector3r& sigmaT=TT*phys->sigmaT;
+ const Real& sigmaN=phys->sigmaN;
+ for(int i=0; i<3; i++) for(int j=0;j<3; j++){
+ ss[i][j]+=d*A*(sigmaN*n[i]*n[j]+.5*(sigmaT[i]*n[j]+sigmaT[j]*n[i]));
+ }
+ }
+ }
+ // divide by approx spatial volume over which we averaged:
+ // spiral cylinder with two half-spherical caps at ends
+ ss*=1/((4/3.)*Mathr::PI*pow(radius,3)+Mathr::PI*pow(radius,2)*(thetaSpan*pt[0]-2*radius));
+ return ss;
+ }
+};
+
#ifdef YADE_VTK
-
-/* Fastly locate interactions withing given distance from a given point. See python docs for details. */
+/* Fastly locate interactions within given distance from a given point. See python docs for details. */
class InteractionLocator{
// object doing the work, see http://www.vtk.org/doc/release/5.2/html/a01048.html
vtkPointLocator *locator;
@@ -168,4 +282,10 @@
.add_property("count",&InteractionLocator::getCnt,"Number of interactions held")
;
#endif
+ py::class_<SpiralInteractionLocator2d>("SpiralInteractionLocator2d",
+ "Locate all real interactions in 2d plane (reduced by spiral projection from 3d, using ``Shop::spiralProject``, which is the same as :yref:`yade.utils.spiralProject`) using their :yref:`contact points<Dem3DofGeom::contactPoint>`. \n\n.. note::\n\tDo not run simulation while using this object.",
+ python::init<Real,int,Real>((python::arg("dH_dTheta"),python::arg("axis")=0,python::arg("theta0")=0),":Parameters:\n\n\tdH_dTheta: float\n\t\tSpiral inclination, i.e. height increase per 1 radian turn;\n\taxis: int\n\t\tAxis of rotation (0=x,1=y,2=z)\n\ttheta: float\n\t\tSpiral angle at zero height (theta intercept)\n\n")
+ )
+ .def("intrsAroundPt",&SpiralInteractionLocator2d::intrsAroundPt,(python::arg("pt2d"),python::arg("radius")),"Return list of interaction objects that are not further from *pt2d* than *radius* in the projection plane")
+ .def("macroStressAroundPt",&SpiralInteractionLocator2d::macroStressAroundPt,(python::arg("pt2d"),python::arg("radius")),"Compute macroscopic stress around given point, rotating the interaction to the projection plane first. The formula used is\n\n.. math::\n\n \\sigma_{ij}=\\frac{1}{V}\\sum_{IJ}d^{IJ}A^{IJ}\\left[\\sigma^{N,IJ}n_i^{IJ}n_j^{IJ}+\\frac{1}{2}\\left(\\sigma_i^{T,IJ}n_j^{IJ}+\\sigma_j^{T,IJ}n_i^{IJ}\\right)\\right]\n\nwhere the sum is taken over volume $V$ containing interactions $IJ$ between spheres $I$ and $J$;\n* $i$, $j$ indices denote Cartesian components of vectors and tensors,\n* $d^{IJ}$ is current distance between spheres $I$ and $J$,\n* $A^{IJ}$ is area of contact $IJ$,\n* $n$ is interaction normal (unit vector pointing from center of $I$ to the center of $J$)\n* $\\sigma^{N,IJ}$ is normal stress (as scalar) in contact $IJ$,\n* $\\sigma^{T,IJ}$ is shear stress in contact $IJ$ in global coordinates.\n\n$\\sigma^{T}$ and $n$ are transformed by angle $\\theta$ as given by :yref:`utils.spiralProject`.");
}
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2010-04-10 15:11:48 +0000
+++ py/_utils.cpp 2010-04-18 20:46:54 +0000
@@ -413,7 +413,7 @@
py::def("PWaveTimeStep",PWaveTimeStep,"Get timestep accoring to the velocity of P-Wave propagation; computed from sphere radii, rigidities and masses.");
- py::def("aabbExtrema",aabbExtrema,(py::arg("cutoff")=0.0,py::arg("centers")=false),"Return coordinates of box enclosing all bodies\n\n:Parameters:\n\tcenters : bool\n\t\tdo not take sphere radii in account, only their centroids\n\tcutoff : float ∈〈0…1〉\n\t\trelative dimension by which the box will be cut away at its boundaries.\n\n:return: (lower corner, upper corner) as (Vector3,Vector3)\n");
+ py::def("aabbExtrema",aabbExtrema,(py::arg("cutoff")=0.0,py::arg("centers")=false),"Return coordinates of box enclosing all bodies\n\n:Parameters:\n\tcenters : bool\n\t\tdo not take sphere radii in account, only their centroids\n\tcutoff : float ∈〈0…1〉\n\t\trelative dimension by which the box will be cut away at its boundaries.\n\n\n:return: (lower corner, upper corner) as (Vector3,Vector3)\n\n");
py::def("ptInAABB",ptInAABB,"Return True/False whether the point (3-tuple) p is within box given by its min (3-tuple) and max (3-tuple) corners");
py::def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(py::args("axis","distFactor"),"Return list of ids for spheres (only) that are on extremal ends of the specimen along given axis; distFactor multiplies their radius so that sphere that do not touch the boundary coordinate can also be returned."));
py::def("approxSectionArea",approxSectionArea,"Compute area of convex hull when when taking (swept) spheres crossing the plane at coord, perpendicular to axis.");
=== modified file 'py/utils.py'
--- py/utils.py 2010-04-10 15:11:48 +0000
+++ py/utils.py 2010-04-18 20:46:54 +0000
@@ -552,6 +552,10 @@
if given, force strained axis, rather than computing it from predominant length
:return: dictionary with keys 'negIds', 'posIds', 'axis', 'area'.
+
+.. warning::
+ The function :yref:`yade.utils.approxSectionArea` uses convex hull algorithm to find the area, but the implementation is reported to be *buggy* (bot works in some cases). Always check this number, or fix the convex hull algorithm (it is documented in the source, see :ysrc:`py/_utils.cpp`).
+
"""
if filename: ids=spheresFromFile(filename,**kw)
else: ids=[b.id for b in O.bodies]
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2010-04-10 15:11:48 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2010-04-18 20:46:54 +0000
@@ -464,15 +464,28 @@
std::ofstream ofs(filename.c_str());
boost::archive::xml_oarchive oa(ofs);
const shared_ptr<Scene>& scene=OMEGA.getScene();
- oa << BOOST_SERIALIZATION_NVP(scene);
+ oa << boost::serialization::make_nvp("scene",*scene);
}
void loadXML(string filename){
std::ifstream ifs(filename.c_str());
boost::archive::xml_iarchive ia(ifs);
shared_ptr<Scene> scene(new Scene);
- ia >> BOOST_SERIALIZATION_NVP(scene); // this might work…
- OMEGA.setScene(scene);
- //throw runtime_error("Not yet implemented");
+ ia >> boost::serialization::make_nvp("scene",*scene);
+ OMEGA.setScene(scene);
+ }
+ // binary format
+ void saveBin(string filename){
+ std::ofstream ofs(filename.c_str());
+ boost::archive::binary_oarchive oa(ofs);
+ const shared_ptr<Scene>& scene=OMEGA.getScene();
+ oa << boost::serialization::make_nvp("scene",*scene);
+ }
+ void loadBin(string filename){
+ std::ifstream ifs(filename.c_str());
+ boost::archive::binary_iarchive ia(ifs);
+ shared_ptr<Scene> scene(new Scene);
+ ia >> boost::serialization::make_nvp("scene",*scene);
+ OMEGA.setScene(scene);
}
#endif
@@ -556,7 +569,9 @@
.def("disableGdb",&pyOmega::disableGdb,"Revert SEGV and ABRT handlers to system defaults.")
#ifdef YADE_BOOST_SERIALIZATION
.def("saveXML",&pyOmega::saveXML,"[EXPERIMENTAL] save to XML using boost::serialization.")
- .def("loadXML",&pyOmega::saveXML,"[EXPERIMENTAL] load from XML using boost::serialization.")
+ .def("loadXML",&pyOmega::loadXML,"[EXPERIMENTAL] load from XML using boost::serialization.")
+ .def("saveBin",&pyOmega::saveBin,"[EXPERIMENTAL] save to binary file using boost::serialization.")
+ .def("loadBin",&pyOmega::loadBin,"[EXPERIMENTAL] load from binary 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).")
.def("tmpFilename",&pyOmega::tmpFilename,"Return unique name of file in temporary directory which will be deleted when yade exits.")
