yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07268
[Branch ~yade-dev/yade/trunk] Rev 2784: Added yade.utils.voxelPorosity and yade.utils.voxelPorosityTriaxial. It allows calculating porosi...
------------------------------------------------------------
revno: 2784
committer: Janek Kozicki <cosurgi@xxxxxxxxxx>
branch nick: yade
timestamp: Wed 2011-03-16 22:02:40 +0100
message:
Added yade.utils.voxelPorosity and yade.utils.voxelPorosityTriaxial. It allows calculating porosity at any arbitrary sub-volume of a whole sample. See included documentation.
modified:
SConstruct
pkg/dem/Shop.cpp
pkg/dem/Shop.hpp
py/_utils.cpp
py/utils.py
--
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 'SConstruct'
--- SConstruct 2011-02-20 19:33:26 +0000
+++ SConstruct 2011-03-16 21:02:40 +0000
@@ -193,7 +193,7 @@
print cmd; os.system(cmd)
sys.argv.remove('tags')
if 'doc' in sys.argv:
- raise RuntimeError("'doc' argument not supported by scons now")
+ raise RuntimeError("'doc' argument not supported by scons now. See doc/README")
# cmd="cd doc; doxygen Doxyfile"
# print cmd; os.system(cmd)
# sys.argv.remove('doc')
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2011-02-28 11:32:03 +0000
+++ pkg/dem/Shop.cpp 2011-03-16 21:02:40 +0000
@@ -302,6 +302,71 @@
return (V-Vs)/V;
}
+Real Shop::getVoxelPorosity(const shared_ptr<Scene>& _scene, int _resolution, Vector3r _start,Vector3r _end){
+ const shared_ptr<Scene> scene=(_scene?_scene:Omega::instance().getScene());
+ if(_start==_end) throw std::invalid_argument("utils.voxelPorosity: cannot calculate porosity when start==end of the volume box.");
+ if(_resolution<50) throw std::invalid_argument("utils.voxelPorosity: it doesn't make sense to calculate porosity with voxel resolution below 50.");
+
+ // prepare the gird, it eats a lot of memory.
+ // I am not optimizing for using bits. A single byte for each cell is used.
+ std::vector<std::vector<std::vector<unsigned char> > > grid;
+ int S(_resolution);
+ grid.resize(S);
+ for(int i=0 ; i<S ; ++i) {
+ grid[i].resize(S);
+ for(int j=0 ; j<S ; ++j) {
+ grid[i][j].resize(S,0);
+ }
+ }
+
+ Vector3r start(_start), size(_end - _start);
+
+ FOREACH(shared_ptr<Body> bi, *scene->bodies){
+ if((bi)->isClump()) continue;
+ const shared_ptr<Body>& b = bi;
+ if ( b->isDynamic() || b->isClumpMember() ) {
+ const shared_ptr<Sphere>& sphere = YADE_PTR_CAST<Sphere> ( b->shape );
+ Real r = sphere->radius;
+ Real rr = r*r;
+ Vector3r pos = b->state->se3.position;
+ // we got sphere with radius r, at position pos.
+ // and a box of size S, scaled to 'size'
+ // mark cells that are iniside a sphere
+ int ii(0),II(S),jj(0),JJ(S),kk(0),KK(S);
+ // make sure to loop only in AABB of that sphere. No need to waste cycles outside of it.
+ ii = std::max((int)((Real)(pos[0]-start[0]-r)*(Real)(S)/(Real)(size[0]))-1,0); II = std::min(ii+(int)((Real)(2.0*r)*(Real)(S)/(Real)(size[0]))+3,S);
+ jj = std::max((int)((Real)(pos[1]-start[1]-r)*(Real)(S)/(Real)(size[1]))-1,0); JJ = std::min(jj+(int)((Real)(2.0*r)*(Real)(S)/(Real)(size[1]))+3,S);
+ kk = std::max((int)((Real)(pos[2]-start[2]-r)*(Real)(S)/(Real)(size[2]))-1,0); KK = std::min(kk+(int)((Real)(2.0*r)*(Real)(S)/(Real)(size[2]))+3,S);
+ for(int i=ii ; i<II ; ++i) {
+ for(int j=jj ; j<JJ ; ++j) {
+ for(int k=kk ; k<KK ; ++k) {
+ Vector3r a(i,j,k);
+ a = a/(Real)(S);
+ Vector3r b( a[0]*size[0] , a[1]*size[1] , a[2]*size[2] );
+ b = b+start;
+ Vector3r c(0,0,0);
+ c = pos-b;
+ Real x = c[0];
+ Real y = c[1];
+ Real z = c[2];
+ if(x*x+y*y+z*z < rr)
+ grid[i][j][k]=1;
+ }}}
+ }
+ }
+
+ Real Vv=0;
+ for(int i=0 ; i<S ; ++i ) {
+ for(int j=0 ; j<S ; ++j ) {
+ for(int k=0 ; k<S ; ++k ) {
+ if(grid[i][j][k]==1)
+ Vv += 1.0;
+ }
+ }
+ }
+
+ return ( std::pow(S,3) - Vv ) / std::pow(S,3);
+};
vector<pair<Vector3r,Real> > Shop::loadSpheresFromFile(const string& fname, Vector3r& minXYZ, Vector3r& maxXYZ, Vector3r* cellSize){
if(!boost::filesystem::exists(fname)) {
=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp 2011-02-28 11:32:03 +0000
+++ pkg/dem/Shop.hpp 2011-03-16 21:02:40 +0000
@@ -56,6 +56,9 @@
//! Compute porosity; volume must be given for aperiodic simulations
static Real getPorosity(const shared_ptr<Scene>& rb=shared_ptr<Scene>(),Real volume=-1);
+ //! Compute porosity by dividing given volume into a grid of voxels;
+ static Real getVoxelPorosity(const shared_ptr<Scene>& rb=shared_ptr<Scene>(),int resolution=500,Vector3r start=Vector3r(0,0,0),Vector3r end=Vector3r(0,0,0));
+
//! Estimate timestep based on P-wave propagation speed
static Real PWaveTimeStep(const shared_ptr<Scene> rb=shared_ptr<Scene>());
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2011-03-04 15:51:30 +0000
+++ py/_utils.cpp 2011-03-16 21:02:40 +0000
@@ -428,6 +428,7 @@
}
Real Shop__getPorosity(Real volume=-1){ return Shop::getPorosity(Omega::instance().getScene(),volume); }
+Real Shop__getVoxelPorosity(int resolution=200, Vector3r start=Vector3r(0,0,0),Vector3r end=Vector3r(0,0,0)){ return Shop::getVoxelPorosity(Omega::instance().getScene(),resolution,start,end); }
Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
py::tuple Shop__fabricTensor(bool splitTensor=false, bool revertSign=false, Real thresholdForce=NaN){return Shop::fabricTensor(splitTensor,revertSign,thresholdForce);}
@@ -456,6 +457,7 @@
py::def("RayleighWaveTimeStep",RayleighWaveTimeStep,"Determination of time step according to Rayleigh wave speed of force propagation.");
py::def("getSpheresVolume",Shop__getSpheresVolume,"Compute the total volume of spheres in the simulation (might crash for now if dynamic bodies are not spheres)");
py::def("porosity",Shop__getPorosity,(py::arg("volume")=-1),"Compute packing poro sity $\\frac{V-V_s}{V}$ where $V$ is overall volume and $V_s$ is volume of spheres.\n\n:param float volume: overall volume which must be specified for aperiodic simulations. For periodic simulations, current volume of the :yref:`Cell` is used.\n");
+ py::def("voxelPorosity",Shop__getVoxelPorosity,(py::arg("resolution")=200,py::arg("start")=Vector3r(0,0,0),py::arg("end")=Vector3r(0,0,0)),"Compute packing porosity $\\frac{V-V_v}{V}$ where $V$ is a specified volume (from start to end) and $V_v$ is volume of voxels that fall inside any sphere. The calculation method is to divide whole volume into a dense grid of voxels (at given resolution), and count the voxels that fall inside any of the spheres. This method allows to calculate porosity in any given sub-volume of a whole sample. It is properly excluding part of a sphere that does not fall inside a specified volume.\n\n:param int resolution: voxel grid resolution, values bigger than resolution=1600 require a 64 bit operating system, because more than 4GB of RAM is used, a resolution=800 will use 500MB of RAM.\n:param Vector3 start: start corner of the volume.\n:param Vector3 end: end corner of the volume.\n");
py::def("aabbExtrema",aabbExtrema,(py::arg("cutoff")=0.0,py::arg("centers")=false),"Return coordinates of box enclosing all bodies\n\n:param bool centers: do not take sphere radii in account, only their centroids\n:param float∈〈0…1〉 cutoff: relative 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",isInBB,"Return True/False whether the point p is within box given by its min and max 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."));
=== modified file 'py/utils.py'
--- py/utils.py 2011-02-17 11:32:00 +0000
+++ py/utils.py 2011-03-16 21:02:40 +0000
@@ -493,6 +493,57 @@
negIds,posIds=negPosExtremeIds(axis=axis,distFactor=2.2)
return {'negIds':negIds,'posIds':posIds,'axis':axis,'area':min(areas)}
+def voxelPorosityTriaxial(triax,resolution=200,offset=0):
+ """
+ Calculate the porosity of a sample, given the TriaxialCompressionEngine.
+
+ A function :yref:`yade.utils.voxelPorosity` is invoked, with the volume of a box enclosed by TriaxialCompressionEngine walls.
+ The additional parameter offset allows using a smaller volume inside the box, where each side of the volume is at offset distance
+ from the walls. By this way it is possible to find a more precise porosity of the sample, since at walls' contact the porosity is usually reduced.
+
+ A recommended value of offset is bigger or equal to the average radius of spheres inside.
+
+ The value of resolution depends on size of spheres used. It can be calibrated by invoking voxelPorosityTriaxial with offset=0 and
+ comparing the result with TriaxialCompressionEngine.porosity. After calibration, the offset can be set to radius, or a bigger value, to get
+ the result.
+
+ :param triax: the TriaxialCompressionEngine handle
+ :param offset: offset distance
+ :param resolution: voxel grid resolution
+ :return: the porosity of the sample inside given volume
+
+ Example invocation:
+ >>> O.engines[5].porosity
+ 0.3901808
+ >>> utils.voxelPorosityTriaxial(O.engines[5],200,0)
+ 0.3908145
+ >>> utils.voxelPorosityTriaxial(O.engines[5],800,0)
+ 0.3902078
+ >>> utils.voxelPorosityTriaxial(O.engines[5],200,avg_radius)
+ 0.3604238
+ >>> utils.voxelPorosityTriaxial(O.engines[5],800,avg_radius)
+ 0.3604175
+
+ """
+ p_bottom = O.bodies[triax.wall_bottom_id].state.se3[0]
+ p_top = O.bodies[triax.wall_top_id ].state.se3[0]
+ p_left = O.bodies[triax.wall_left_id ].state.se3[0]
+ p_right = O.bodies[triax.wall_right_id ].state.se3[0]
+ p_front = O.bodies[triax.wall_front_id ].state.se3[0]
+ p_back = O.bodies[triax.wall_back_id ].state.se3[0]
+ th = (triax.thickness)*0.5+offset
+ x_0 = p_left [0] + th
+ x_1 = p_right [0] - th
+ y_0 = p_bottom[1] + th
+ y_1 = p_top [1] - th
+ z_0 = p_back [2] + th
+ z_1 = p_front [2] - th
+ a=Vector3(x_0,y_0,z_0)
+ b=Vector3(x_1,y_1,z_1)
+ return voxelPorosity(resolution,a,b)
+
+
+
def NormalRestitution2DampingRate(en):
r"""Compute the normal damping rate as a function of the normal coefficient of restitution $e_n$. For $e_n\in\langle0,1\rangle$ damping rate equals
Follow ups
