yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10831
[Branch ~yade-pkg/yade/git-trunk] Rev 3957: Polyhedra implementation improvement (Contributed by Jan Elias)
------------------------------------------------------------
revno: 3957
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Fri 2014-05-16 13:58:59 +0200
message:
Polyhedra implementation improvement (Contributed by Jan Elias)
modified:
pkg/dem/Polyhedra.cpp
pkg/dem/Polyhedra.hpp
pkg/dem/Polyhedra_Ig2.cpp
pkg/dem/Polyhedra_splitter.cpp
pkg/dem/Polyhedra_support.cpp
py/_polyhedra_utils.cpp
py/polyhedra_utils.py
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp 2013-10-17 11:59:10 +0000
+++ pkg/dem/Polyhedra.cpp 2014-05-16 11:58:59 +0000
@@ -16,7 +16,7 @@
YADE_PLUGIN(/* self-contained in hpp: */ (Polyhedra) (PolyhedraGeom) (Bo1_Polyhedra_Aabb) (PolyhedraPhys) (PolyhedraMat) (Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys) (PolyhedraVolumetricLaw)
/* some code in cpp (this file): */
#ifdef YADE_OPENGL
- (Gl1_Polyhedra)
+ (Gl1_Polyhedra) (Gl1_PolyhedraGeom) (Gl1_PolyhedraPhys)
#endif
);
@@ -197,9 +197,9 @@
Triangulation dt(nuclei.begin(), nuclei.end());
Triangulation::Vertex_handle zero_point = dt.insert(CGALpoint(5.,5.,5.));
v.clear();
- std::vector<Triangulation::Cell_handle> ch_cells;
+ std::vector<typename Triangulation::Cell_handle> ch_cells;
dt.incident_cells(zero_point,std::back_inserter(ch_cells));
- for(std::vector<Triangulation::Cell_handle>::iterator ci = ch_cells.begin(); ci !=ch_cells.end(); ++ci){
+ for(std::vector<typename Triangulation::Cell_handle>::iterator ci = ch_cells.begin(); ci !=ch_cells.end(); ++ci){
v.push_back(FromCGALPoint(dt.dual(*ci))-Vector3r(5.,5.,5.));
}
@@ -256,14 +256,15 @@
#ifdef YADE_OPENGL
#include<yade/lib/opengl/OpenGLWrapper.hpp>
- void Gl1_Polyhedra::go(const shared_ptr<Shape>& cm, const shared_ptr<State>&,bool,const GLViewInfo&)
+ bool Gl1_Polyhedra::wire;
+ void Gl1_Polyhedra::go(const shared_ptr<Shape>& cm, const shared_ptr<State>&,bool wire2,const GLViewInfo&)
{
glMaterialv(GL_FRONT,GL_AMBIENT_AND_DIFFUSE,Vector3r(cm->color[0],cm->color[1],cm->color[2]));
glColor3v(cm->color);
Polyhedra* t=static_cast<Polyhedra*>(cm.get());
vector<int> faceTri = t->GetSurfaceTriangulation();
- if (0) {
+ if (wire || wire2) {
glDisable(GL_LIGHTING);
glBegin(GL_LINES);
#define __ONEWIRE(a,b) glVertex3v(t->v[a]);glVertex3v(t->v[b])
@@ -283,6 +284,58 @@
glEnd();
}
}
+
+
+ void Gl1_PolyhedraGeom::go(const shared_ptr<IGeom>& ig, const shared_ptr<Interaction>&, const shared_ptr<Body>&, const shared_ptr<Body>&, bool){ draw(ig);}
+
+ void Gl1_PolyhedraGeom::draw(const shared_ptr<IGeom>& ig){
+ };
+
+ GLUquadric* Gl1_PolyhedraPhys::gluQuadric=NULL;
+ Real Gl1_PolyhedraPhys::maxFn;
+ Real Gl1_PolyhedraPhys::refRadius;
+ Real Gl1_PolyhedraPhys::maxRadius;
+ int Gl1_PolyhedraPhys::signFilter;
+ int Gl1_PolyhedraPhys::slices;
+ int Gl1_PolyhedraPhys::stacks;
+
+ void Gl1_PolyhedraPhys::go(const shared_ptr<IPhys>& ip, const shared_ptr<Interaction>& i, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
+ if(!gluQuadric){ gluQuadric=gluNewQuadric(); if(!gluQuadric) throw runtime_error("Gl1_PolyhedraPhys::go unable to allocate new GLUquadric object (out of memory?)."); }
+ PolyhedraPhys* np=static_cast<PolyhedraPhys*>(ip.get());
+ shared_ptr<IGeom> ig(i->geom); if(!ig) return; // changed meanwhile?
+ PolyhedraGeom* geom=YADE_CAST<PolyhedraGeom*>(ig.get());
+ Real fnNorm=np->normalForce.dot(geom->normal);
+ if((signFilter>0 && fnNorm<0) || (signFilter<0 && fnNorm>0)) return;
+ int fnSign=fnNorm>0?1:-1;
+ fnNorm=abs(fnNorm);
+ Real radiusScale=1.;
+ maxFn=max(fnNorm,maxFn);
+ Real realMaxRadius;
+ if(maxRadius<0){
+ refRadius=min(0.03,refRadius);
+ realMaxRadius=refRadius;
+ }
+ else realMaxRadius=maxRadius;
+ Real radius=radiusScale*realMaxRadius*(fnNorm/maxFn);
+ if (radius<=0.) radius = 1E-8;
+ Vector3r color=Shop::scalarOnColorScale(fnNorm*fnSign,-maxFn,maxFn);
+
+ Vector3r p1=b1->state->pos, p2=b2->state->pos;
+ Vector3r relPos;
+ relPos=p2-p1;
+ Real dist=relPos.norm();
+
+ glDisable(GL_CULL_FACE);
+ glPushMatrix();
+ glTranslatef(p1[0],p1[1],p1[2]);
+ Quaternionr q(Quaternionr().setFromTwoVectors(Vector3r(0,0,1),relPos/dist /* normalized */));
+ // using Transform with OpenGL: http://eigen.tuxfamily.org/dox/TutorialGeometry.html
+ //glMultMatrixd(Eigen::Affine3d(q).data());
+ glMultMatrix(Eigen::Transform<Real,3,Eigen::Affine>(q).data());
+ glColor3v(color);
+ gluCylinder(gluQuadric,radius,radius,dist,slices,stacks);
+ glPopMatrix();
+ }
#endif
//**********************************************************************************
@@ -362,9 +415,9 @@
// Apply forces on polyhedrons in collision based on geometric configuration
void PolyhedraVolumetricLaw::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
- if (!I->geom) {I->phys=shared_ptr<IPhys>(); return;}
+ if (!I->geom) {return;}
const shared_ptr<PolyhedraGeom>& contactGeom(dynamic_pointer_cast<PolyhedraGeom>(I->geom));
- if(!contactGeom) {I->phys=shared_ptr<IPhys>(); return;}
+ if(!contactGeom) {return;}
const Body::id_t idA=I->getId1(), idB=I->getId2();
const shared_ptr<Body>& A=Body::byId(idA), B=Body::byId(idB);
@@ -373,11 +426,12 @@
//erase the interaction when aAbB shows separation, otherwise keep it to be able to store previous separating plane for fast detection of separation
if (A->bound->min[0] >= B->bound->max[0] || B->bound->min[0] >= A->bound->max[0] || A->bound->min[1] >= B->bound->max[1] || B->bound->min[1] >= A->bound->max[1] || A->bound->min[2] >= B->bound->max[2] || B->bound->min[2] >= A->bound->max[2]) {
scene->interactions->requestErase(I);
+ phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.);
return;
}
//zero penetration depth means no interaction force
- if(!(contactGeom->penetrationVolume > 0) ) {I->phys=shared_ptr<IPhys>(); return;}
+ if(!(contactGeom->equivalentPenetrationDepth > 1E-18) || !(contactGeom->penetrationVolume > 0)) {phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.); return;}
Vector3r normalForce=contactGeom->normal*contactGeom->penetrationVolume*phys->kn;
//shear force: in case the polyhdras are separated and come to contact again, one should not use the previous shear force
@@ -418,7 +472,22 @@
scene->forces.addForce (idB, -F);
scene->forces.addTorque(idA, -(A->state->pos-contactGeom->contactPoint).cross(F));
scene->forces.addTorque(idB, (B->state->pos-contactGeom->contactPoint).cross(F));
-
+
+ /*
+ FILE * fin = fopen("Forces.dat","a");
+ fprintf(fin,"************** IDS %d %d **************\n",idA, idB);
+ Vector3r T = (B->state->pos-contactGeom->contactPoint).cross(F);
+ fprintf(fin,"volume\t%e\n",contactGeom->penetrationVolume);
+ fprintf(fin,"normal_force\t%e\t%e\t%e\n",normalForce[0],normalForce[1],normalForce[2]);
+ fprintf(fin,"shear_force\t%e\t%e\t%e\n",shearForce[0],shearForce[1],shearForce[2]);
+ fprintf(fin,"total_force\t%e\t%e\t%e\n",F[0],F[1],F[2]);
+ fprintf(fin,"torsion\t%e\t%e\t%e\n",T[0],T[1],T[2]);
+ fprintf(fin,"A\t%e\t%e\t%e\n",A->state->pos[0],A->state->pos[1],A->state->pos[2]);
+ fprintf(fin,"B\t%e\t%e\t%e\n",B->state->pos[0],B->state->pos[1],B->state->pos[2]);
+ fprintf(fin,"centroid\t%e\t%e\t%e\n",contactGeom->contactPoint[0],contactGeom->contactPoint[1],contactGeom->contactPoint[2]);
+ fclose(fin);
+ */
+
//needed to be able to acces interaction forces in other parts of yade
phys->normalForce = normalForce;
phys->shearForce = shearForce;
=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp 2014-01-28 08:43:35 +0000
+++ pkg/dem/Polyhedra.hpp 2014-05-16 11:58:59 +0000
@@ -145,19 +145,6 @@
REGISTER_SERIALIZABLE(Bo1_Polyhedra_Aabb);
//***************************************************************************
-#ifdef YADE_OPENGL
- #include<yade/pkg/common/GLDrawFunctors.hpp>
- /*! Draw Polyhedra using OpenGL */
- class Gl1_Polyhedra: public GlShapeFunctor{
- public:
- virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>&,bool,const GLViewInfo&);
- YADE_CLASS_BASE_DOC(Gl1_Polyhedra,GlShapeFunctor,"Renders :yref:`Polyhedra` object");
- RENDERS(Polyhedra);
- };
- REGISTER_SERIALIZABLE(Gl1_Polyhedra);
-#endif
-
-//***************************************************************************
/*! Elastic material */
class PolyhedraMat: public Material{
public:
@@ -193,6 +180,53 @@
REGISTER_SERIALIZABLE(PolyhedraPhys);
//***************************************************************************
+#ifdef YADE_OPENGL
+ #include<yade/pkg/common/GLDrawFunctors.hpp>
+ #include<yade/lib/opengl/OpenGLWrapper.hpp>
+ #include<yade/lib/opengl/GLUtils.hpp>
+ #include<GL/glu.h>
+ #include<yade/pkg/dem/Shop.hpp>
+
+ /*! Draw Polyhedra using OpenGL */
+ class Gl1_Polyhedra: public GlShapeFunctor{
+ public:
+ virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>&,bool,const GLViewInfo&);
+ YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_Polyhedra,GlShapeFunctor,"Renders :yref:`Polyhedra` object",
+ ((bool,wire,false,,"Only show wireframe"))
+ );
+ RENDERS(Polyhedra);
+ };
+ REGISTER_SERIALIZABLE(Gl1_Polyhedra);
+
+ struct Gl1_PolyhedraGeom: public GlIGeomFunctor{
+ RENDERS(PolyhedraGeom);
+ void go(const shared_ptr<IGeom>&, const shared_ptr<Interaction>&, const shared_ptr<Body>&, const shared_ptr<Body>&, bool);
+ void draw(const shared_ptr<IGeom>&);
+ YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_PolyhedraGeom,GlIGeomFunctor,"Render :yref:`PolyhedraGeom` geometry.",
+ );
+ };
+ REGISTER_SERIALIZABLE(Gl1_PolyhedraGeom);
+
+ class Gl1_PolyhedraPhys: public GlIPhysFunctor{
+ static GLUquadric* gluQuadric; // needed for gluCylinder, initialized by ::go if no initialized yet
+ public:
+ virtual void go(const shared_ptr<IPhys>&,const shared_ptr<Interaction>&,const shared_ptr<Body>&,const shared_ptr<Body>&,bool wireFrame);
+ YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_PolyhedraPhys,GlIPhysFunctor,"Renders :yref:`PolyhedraPhys` objects as cylinders of which diameter and color depends on :yref:`PolyhedraPhys::normForce` magnitude.",
+ ((Real,maxFn,0,,"Value of :yref:`NormPhys.normalForce` corresponding to :yref:`maxDiameter<Gl1_NormPhys.maxDiameter>`. This value will be increased (but *not decreased* ) automatically."))
+ ((Real,refRadius,std::numeric_limits<Real>::infinity(),,"Reference (minimum) particle radius"))
+ ((int,signFilter,0,,"If non-zero, only display contacts with negative (-1) or positive (+1) normal forces; if zero, all contacts will be displayed."))
+ ((Real,maxRadius,-1,,"Cylinder radius corresponding to the maximum normal force."))
+ ((int,slices,6,,"Number of sphere slices; (see `glutCylinder reference <http://www.opengl.org/sdk/docs/man/xhtml/gluCylinder.xml>`__)"))
+ ( (int,stacks,1,,"Number of sphere stacks; (see `glutCylinder reference <http://www.opengl.org/sdk/docs/man/xhtml/gluCylinder.xml>`__)"))
+ );
+ RENDERS(PolyhedraPhys);
+ };
+ REGISTER_SERIALIZABLE(Gl1_PolyhedraPhys);
+
+#endif
+
+
+//***************************************************************************
class Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys: public IPhysFunctor{
public:
virtual void go(const shared_ptr<Material>& b1,
@@ -247,6 +281,8 @@
Polyhedron Polyhedron_Polyhedron_intersection(Polyhedron A, Polyhedron B, CGALpoint X, CGALpoint centroidA, CGALpoint centroidB, std::vector<int> &code);
//return intersection of plane & polyhedron
Polyhedron Polyhedron_Plane_intersection(Polyhedron A, Plane B, CGALpoint centroid, CGALpoint X);
+//Test if point is inside Polyhedron
+bool Is_inside_Polyhedron(Polyhedron P, CGALpoint inside);
//return approximate intersection of sphere & polyhedron
bool Sphere_Polyhedron_intersection(Polyhedron A, double r, CGALpoint C, CGALpoint centroid, double volume, CGALvector normal, double area);
//return volume and centroid of polyhedra
@@ -263,12 +299,13 @@
Polyhedron Simplify(Polyhedron P, double lim);
//list of facets and edges
void PrintPolyhedron(Polyhedron P);
+void PrintPolyhedron2File(Polyhedron P,FILE* X);
//normal by least square fitting of separating segments
Vector3r FindNormal(Polyhedron Int, Polyhedron PA, Polyhedron PB);
//calculate area of projection of polyhedron into the plane
double CalculateProjectionArea(Polyhedron Int, CGALvector CGALnormal);
//split polyhedron
-void SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction);
+shared_ptr<Body> SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction, Vector3r point);
//new polyhedra
shared_ptr<Body> NewPolyhedra(vector<Vector3r> v, shared_ptr<Material> mat);
=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp 2013-10-16 15:24:49 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp 2014-05-16 11:58:59 +0000
@@ -29,6 +29,7 @@
Transformation t_rot_trans(rot_mat(0,0),rot_mat(0,1),rot_mat(0,2), trans_vec[0],rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),trans_vec[1],rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),trans_vec[2],1.);
Polyhedron PA = A->GetPolyhedron();
std::transform( PA.points_begin(), PA.points_end(), PA.points_begin(), t_rot_trans);
+ std::transform( PA.facets_begin(), PA.facets_end(), PA.planes_begin(),Plane_equation());
//move and rotate 2nd the CGAL structure Polyhedron
@@ -37,6 +38,7 @@
t_rot_trans = Transformation(rot_mat(0,0),rot_mat(0,1),rot_mat(0,2), trans_vec[0],rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),trans_vec[1],rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),trans_vec[2],1.);
Polyhedron PB = B->GetPolyhedron();
std::transform( PB.points_begin(), PB.points_end(), PB.points_begin(), t_rot_trans);
+ std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());
shared_ptr<PolyhedraGeom> bang;
if (isNew) {
@@ -62,6 +64,7 @@
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > min(A->GetVolume(),B->GetVolume())) {bang->equivalentPenetrationDepth=0; return true;}
+ if ( (!Is_inside_Polyhedron(PA, ToCGALPoint(centroid))) or (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid)))) {bang->equivalentPenetrationDepth=0; return true;}
//find normal direction
Vector3r normal = FindNormal(Int, PA, PB);
@@ -80,6 +83,17 @@
bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
bang->normal=normal;
+ /*
+ FILE * fin = fopen("Interactions.dat","a");
+ fprintf(fin,"************** IDS %d %d **************\n",interaction->id1, interaction->id2);
+ fprintf(fin,"volume\t%e\n",volume);
+ fprintf(fin,"centroid\t%e\t%e\t%e\n",centroid[0],centroid[1],centroid[2]);
+ PrintPolyhedron2File(PA,fin);
+ PrintPolyhedron2File(PB,fin);
+ PrintPolyhedron2File(Int,fin);
+ fclose(fin);
+ */
+
return true;
}
@@ -131,6 +145,7 @@
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
+ if (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid))) {bang->equivalentPenetrationDepth=0; return true;}
//calculate area of projection of Intersection into the normal plane
double area = volume/1E-8;
@@ -167,6 +182,7 @@
Transformation t_rot_trans(rot_mat(0,0),rot_mat(0,1),rot_mat(0,2), trans_vec[0],rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),trans_vec[1],rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),trans_vec[2],1.);
Polyhedron PB = B->GetPolyhedron();
std::transform( PB.points_begin(), PB.points_end(), PB.points_begin(), t_rot_trans);
+ std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());
//move and rotate facet
vector<CGALpoint> v;
@@ -224,6 +240,7 @@
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
+ if (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid))) {bang->equivalentPenetrationDepth=0; return true;}
//find normal direction
Vector3r normal = FindNormal(Int, PA, PB);
=== modified file 'pkg/dem/Polyhedra_splitter.cpp'
--- pkg/dem/Polyhedra_splitter.cpp 2013-10-16 15:24:49 +0000
+++ pkg/dem/Polyhedra_splitter.cpp 2014-05-16 11:58:59 +0000
@@ -22,11 +22,9 @@
PolyhedraPhys* phys=YADE_CAST<PolyhedraPhys*>(I->phys.get());
if(!geom || !phys) continue;
Vector3r f=phys->normalForce+phys->shearForce;
- //Sum f_i*l_j for each contact of each particle
+ //Sum f_i*l_j for each contact of each particle
bStresses[I->getId1()] -=f*((geom->contactPoint-Body::byId(I->getId1(),scene)->state->pos).transpose());
bStresses[I->getId2()] +=f*((geom->contactPoint-Body::byId(I->getId2(),scene)->state->pos).transpose());
-
-
}
}
@@ -52,6 +50,20 @@
bStress(2,1) = bStress(1,2);
}
+//**********************************************************************************
+//split polyhedra
+void SplitPolyhedraDouble(const shared_ptr<Body>& body, Vector3r direction1, Vector3r direction2){
+
+ const Se3r& se3=body->state->se3;
+ Vector3r point = se3.position;
+
+ shared_ptr<Body> B2 = SplitPolyhedra(body, direction1, point);
+ shared_ptr<Body> B3 = SplitPolyhedra(B2, direction2, point);
+ shared_ptr<Body> B4 = SplitPolyhedra(body, direction2, point);
+
+}
+
+
//*********************************************************************************
/* Split if stress exceed strength */
@@ -61,7 +73,7 @@
shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
vector<shared_ptr<Body> > bodies;
- vector<Vector3r > directions;
+ vector<Vector3r > directions1, directions2;
vector<double > sigmas;
@@ -69,16 +81,14 @@
vector<Matrix3r> bStresses;
getStressForEachBody(bStresses);
- int i = -1;
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
- i++;
if(!b || !b->material || !b->shape) continue;
shared_ptr<Polyhedra> p=dynamic_pointer_cast<Polyhedra>(b->shape);
shared_ptr<PolyhedraMat> m=dynamic_pointer_cast<PolyhedraMat>(b->material);
if(p && m->IsSplitable){
//not real strees, to get real one, it has to be divided by body volume
- Matrix3r stress = bStresses[i];
+ Matrix3r stress = bStresses[b->id];
//get eigenstresses
Symmetrize(stress);
@@ -93,16 +103,18 @@
//division of stress by volume
double comp_stress = I_valu(min_i,min_i)/p->GetVolume();
- double tens_stress = I_valu(max_i,max_i)/p->GetVolume();
+ double tens_stress = I_valu(max_i,max_i)/p->GetVolume();
Vector3r dirC = I_vect.col(max_i);
Vector3r dirT = I_vect.col(min_i);
- Vector3r dir = dirC.normalized() + dirT.normalized();;
- double sigma_t = -comp_stress/2.+ tens_stress;
- if (sigma_t > getStrength(p->GetVolume(),m->GetStrength())) {bodies.push_back(b); directions.push_back(dir); sigmas.push_back(sigma_t);};
+ Vector3r dir1 = dirC.normalized() + dirT.normalized();
+ Vector3r dir2 = dirC.normalized() - dirT.normalized();
+ //double sigma_t = -comp_stress/2.+ tens_stress;
+ double sigma_t = pow((pow(I_valu(0,0)-I_valu(1,1),2)+pow(I_valu(0,0)-I_valu(2,2),2)+pow(I_valu(1,1)-I_valu(2,2),2))/2.,0.5)/p->GetVolume();
+ if (sigma_t > getStrength(p->GetVolume(),m->GetStrength())) {bodies.push_back(b); directions1.push_back(dir1.normalized()); directions2.push_back(dir2.normalized()); sigmas.push_back(sigma_t);};
}
}
for(int i=0; i<int(bodies.size()); i++){
- SplitPolyhedra(bodies[i], directions[i]);
+ SplitPolyhedraDouble(bodies[i], directions1[i], directions2[i]);
}
}
=== modified file 'pkg/dem/Polyhedra_support.cpp'
--- pkg/dem/Polyhedra_support.cpp 2013-10-16 15:24:49 +0000
+++ pkg/dem/Polyhedra_support.cpp 2014-05-16 11:58:59 +0000
@@ -10,13 +10,15 @@
//EMPRIRICAL CONSTANTS - ADJUST IF SEGMENTATION FAULT OCCUR, IT IS A PROBLEM OF CGAL. THESE ARE USED TO CHECK CGAL IMPUTS
//DISTANCE_LIMIT controls numerical issues in calculating intersection. It should be small enough to neglect only extremely small overlaps, but large enough to prevent errors during computation of convex hull
-#define DISTANCE_LIMIT 1E-11 //-11
+#define DISTANCE_LIMIT 2E-11 //-11
//MERGE_PLANES_LIMIT - if two facets of two intersecting polyhedron differ less, then they are treated ose one only
#define MERGE_PLANES_LIMIT 1E-18 //18
//SIMPLIFY_LIMIT - if two facets of one polyhedron differ less, then they are joint into one facet
-#define SIMPLIFY_LIMIT 1E-20 //19
+#define SIMPLIFY_LIMIT 1E-19 //19
//FIND_NORMAL_LIMIT - to determine which facet of intersection belongs to which polyhedron
#define FIND_NORMAL_LIMIT 1E-40
+//SPLITTER_GAP - creates gap between splitted polyhedrons
+#define SPLITTER_GAP 1E-8
//**********************************************************************************
@@ -288,6 +290,30 @@
//**********************************************************************************
//list of facets + edges
+void PrintPolyhedron2File(Polyhedron P,FILE* X){
+ Vector3r A,B,C;
+ fprintf(X,"*** faces ***\n");
+ for (Polyhedron::Facet_iterator fIter = P.facets_begin(); fIter!=P.facets_end(); ++fIter){
+ Polyhedron::Halfedge_around_facet_circulator hfc0;
+ hfc0 = fIter->facet_begin();
+ int n = fIter->facet_degree();
+ A = FromCGALPoint(hfc0->vertex()->point());
+ C = FromCGALPoint(hfc0->next()->vertex()->point());
+ for (int i=2; i<n; i++){
+ ++hfc0;
+ B = C;
+ C = FromCGALPoint(hfc0->next()->vertex()->point());
+ fprintf(X,"%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",A[0],A[1],A[2],B[0],B[1],B[2],C[0],C[1],C[2]);
+ }
+ }
+ fprintf(X,"*** edges ***\n");
+ for (Polyhedron::Edge_iterator hei = P.edges_begin(); hei!=P.edges_end(); ++hei){
+ fprintf(X,"%e\t%e\t%e\t%e\t%e\t%e\n",hei->vertex()->point()[0],hei->vertex()->point()[1],hei->vertex()->point()[2], hei->opposite()->vertex()->point()[0], hei->opposite()->vertex()->point()[1], hei->opposite()->vertex()->point()[2]);
+ }
+}
+
+//**********************************************************************************
+//list of facets + edges
void PrintPolyhedron(Polyhedron P){
Vector3r A,B,C;
cout << "*** faces ***" << endl;
@@ -478,6 +504,7 @@
// test if do intersect, find some intersecting point
bool intersection_found = false;
double lim = pow(DISTANCE_LIMIT,2);
+ std::transform( A.facets_begin(), A.facets_end(), A.planes_begin(),Plane_equation());
// test centroid of previous intersection
if(Is_inside_Polyhedron(A,X,DISTANCE_LIMIT) && Oriented_squared_distance(B,X)<=-lim) {
intersection_found = true;
@@ -486,20 +513,24 @@
} else {
for(Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end() && !intersection_found ; vIter++){
if(Oriented_squared_distance(B,vIter->point())<=-lim) {
- intersection_found = true;
- CGAL::Object result = CGAL::intersection(Line(vIter->point(),centroid), B);
- if (const CGALpoint *ipoint = CGAL::object_cast<CGALpoint>(&result)) {
- inside = vIter->point() + 0.5*CGALvector(vIter->point(),*ipoint);
+ if (Oriented_squared_distance(B,centroid)<lim/10.){
+ inside = vIter->point() + 0.5*CGALvector(vIter->point(),centroid);
}else{
- cerr << "Error in line-plane intersection" << endl;
+ CGAL::Object result = CGAL::intersection(Line(vIter->point(),centroid), B);
+ if (const CGALpoint *ipoint = CGAL::object_cast<CGALpoint>(&result)) {
+ inside = vIter->point() + 0.5*CGALvector(vIter->point(),*ipoint);
+ }else{
+ cerr << "Error in line-plane intersection" << endl;
+ }
}
+ if(Is_inside_Polyhedron(A,inside,DISTANCE_LIMIT) && Oriented_squared_distance(B,inside)<=-lim) intersection_found = true;
}
}
}
//no intersectiong point => no intersection polyhedron
if(!intersection_found) return Intersection;
-
+
//set the intersection point to origin
Transformation transl_back(CGAL::TRANSLATION, inside - CGALpoint(0.,0.,0.));
Transformation transl(CGAL::TRANSLATION, CGALpoint(0.,0.,0.)-inside);
@@ -520,7 +551,7 @@
//compute convex hull of it
Intersection = ConvexHull(dual_planes);
- if (Intersection.empty()) {cout << "0" << endl; cout.flush(); return Intersection;}
+ if (Intersection.empty()) return Intersection;
//simplify
std::transform( Intersection.facets_begin(), Intersection.facets_end(), Intersection.planes_begin(),Plane_equation());
@@ -542,7 +573,7 @@
}
//**********************************************************************************
-//returnes intersecting polyhedron of polyhedron & plane defined by diriction and point
+//returnes intersecting polyhedron of polyhedron & plane defined by direction and point
Polyhedron Polyhedron_Plane_intersection(Polyhedron A, Vector3r direction, Vector3r plane_point){
Plane B(ToCGALPoint(plane_point), ToCGALVector(direction));
CGALpoint X = ToCGALPoint(plane_point) - 1E-8*ToCGALVector(direction);
@@ -633,7 +664,7 @@
//compute convex hull of it
Intersection = ConvexHull(dual_planes);
- if (Intersection.empty()) {cout << "0" << endl; cout.flush(); return Intersection;};
+ if (Intersection.empty()) return Intersection;
//simplify
std::transform( Intersection.facets_begin(), Intersection.facets_end(), Intersection.planes_begin(),Plane_equation());
@@ -759,10 +790,11 @@
return body;
}
-
//**********************************************************************************
//split polyhedra
-void SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction){
+shared_ptr<Body> SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction, Vector3r point){
+
+ direction.normalize();
const Se3r& se3=body->state->se3;
Polyhedra* A = static_cast<Polyhedra*>(body->shape.get());
@@ -779,10 +811,14 @@
Polyhedron PA = A->GetPolyhedron();
std::transform( PA.points_begin(), PA.points_end(), PA.points_begin(), t_rot_trans);
- //calculate splitted polyhedrons
- Polyhedron S1 = Polyhedron_Plane_intersection(PA, direction, trans_vec);
- Polyhedron S2 = Polyhedron_Plane_intersection(PA, (-1)*direction, trans_vec);
-
+ //calculate first splitted polyhedrons
+ Plane B(ToCGALPoint(point-direction*SPLITTER_GAP), ToCGALVector(direction));
+ Polyhedron S1 = Polyhedron_Plane_intersection(PA, B, ToCGALPoint(se3.position), B.projection(ToCGALPoint(OrigPos)) - 1E-6*ToCGALVector(direction));
+ B = Plane(ToCGALPoint(point+direction*SPLITTER_GAP), ToCGALVector((-1.)*direction));
+ Polyhedron S2 = Polyhedron_Plane_intersection(PA, B, ToCGALPoint(se3.position), B.projection(ToCGALPoint(OrigPos)) + 1E-6*ToCGALVector(direction));
+ Scene* scene=Omega::instance().getScene().get();
+ //scene->bodies->erase(body->id);
+
//replace original polyhedron
A->Clear();
@@ -793,6 +829,8 @@
X->mass=body->material->density*A->GetVolume();
X->refPos[0] = X->refPos[0]+1.;
X->inertia =A->GetInertia()*body->material->density;
+ X->vel = OrigVel + OrigAngVel.cross(X->pos-OrigPos);
+ X->angVel = OrigAngVel;
//second polyhedron
vector<Vector3r> v2;
@@ -800,18 +838,12 @@
shared_ptr<Body> BP = NewPolyhedra(v2, body->material);
//BP->shape->color = body->shape->color;
BP->shape->color = Vector3r(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
-
- //append to the rest
- Scene* scene=Omega::instance().getScene().get();
- //scene->bodies->erase(body->id);
scene->bodies->insert(BP);
-
//set proper state variables
- X->vel = OrigVel + OrigAngVel.cross(X->pos-OrigPos);
BP->state->vel = OrigVel + OrigAngVel.cross(BP->state->pos-OrigPos);
- X->angVel = OrigAngVel;
BP->state->angVel = OrigAngVel;
+ return BP;
}
#endif // YADE_CGAL
=== modified file 'py/_polyhedra_utils.cpp'
--- py/_polyhedra_utils.cpp 2014-04-02 15:19:41 +0000
+++ py/_polyhedra_utils.cpp 2014-05-16 11:58:59 +0000
@@ -103,8 +103,7 @@
Real density=b->state->mass/p->GetVolume();
//get equivalent radius and use same equation as for sphere
Real equi_radius=pow(p->GetVolume()/((4./3.)*Mathr::PI),1./3.);
- //dt=min(dt,equi_radius/sqrt(ebp->Kn*equi_radius/density));
- dt=min(dt,equi_radius/sqrt(ebp->Kn/equi_radius/density));
+ dt=min(dt,equi_radius/sqrt(ebp->Kn*equi_radius/density));
}
}
if (dt==std::numeric_limits<Real>::infinity()) {
@@ -326,10 +325,182 @@
return v;
}
+//**************************************************************************
+/* Generate truncated icosahedron*/
+vector<Vector3r> TruncIcosaHedPoints(Vector3r radii){
+ vector<Vector3r> v;
+
+ double p = (1.+sqrt(5.))/2.;
+ Vector3r f,c,b;
+ f = radii/sqrt(9.*p + 1.);
+ vector<Vector3r> A,B;
+ A.push_back(Vector3r(0.,1.,3.*p)); A.push_back(Vector3r(2.,1.+2.*p,p)); A.push_back(Vector3r(1.,2.+p,2.*p));
+ for(int i=0;i<(int) A.size() ;i++){
+ B.clear();
+ c = Vector3r(A[i][0]*f[0],A[i][1]*f[1],A[i][2]*f[2]);
+ B.push_back(c); B.push_back(Vector3r(c[1],c[2],c[0])); B.push_back(Vector3r(c[2],c[0],c[1]));
+ for(int j=0;j<(int) B.size() ;j++){
+ b = B[j];
+ v.push_back(b);
+ if (b[0] != 0.){
+ v.push_back(Vector3r(-b[0], b[1], b[2]));
+ if (b[1] != 0.){
+ v.push_back(Vector3r(-b[0],-b[1], b[2]));
+ if (b[2] != 0.) v.push_back(Vector3r(-b[0],-b[1],-b[2]));
+ }
+ if (b[2] != 0.) v.push_back(Vector3r(-b[0], b[1],-b[2]));
+ }
+ if (b[1] != 0.){
+ v.push_back(Vector3r( b[0],-b[1], b[2]));
+ if (b[2] != 0.) v.push_back(Vector3r( b[0],-b[1],-b[2]));
+ }
+ if (b[2] != 0.) v.push_back(Vector3r( b[0], b[1],-b[2]));
+ }
+ }
+ return v;
+}
+
+//**************************************************************************
+/* Generate SnubCube*/
+vector<Vector3r> SnubCubePoints(Vector3r radii){
+ vector<Vector3r> v;
+
+ double c1 = 0.337754;
+ double c2 = 1.14261;
+ double c3 = 0.621226;
+ Vector3r f,b;
+ f = radii/1.3437133737446;
+ vector<Vector3r> A;
+ A.push_back(Vector3r(c2,c1,c3)); A.push_back(Vector3r(c1,c3,c2)); A.push_back(Vector3r(c3,c2,c1));
+ A.push_back(Vector3r(-c1,-c2,-c3)); A.push_back(Vector3r(-c2,-c3,-c1)); A.push_back(Vector3r(-c3,-c1,-c2));
+ for(int i=0;i<(int) A.size();i++){
+ b = Vector3r(A[i][0]*f[0],A[i][1]*f[1],A[i][2]*f[2]);
+ v.push_back(b);
+ v.push_back(Vector3r(-b[0],-b[1], b[2]));
+ v.push_back(Vector3r(-b[0], b[1],-b[2]));
+ v.push_back(Vector3r( b[0],-b[1],-b[2]));
+ }
+ return v;
+}
+
+//**************************************************************************
+/* Generate ball*/
+vector<Vector3r> BallPoints(Vector3r radii, int NumFacets,int seed){
+ vector<Vector3r> v;
+ if (NumFacets == 60) v = TruncIcosaHedPoints(radii);
+ if (NumFacets == 24) v = SnubCubePoints(radii);
+ else{
+ double inc = Mathr::PI * (3. - pow(5.,0.5));
+ double off = 2. / double(NumFacets);
+ double y,r,phi;
+ for(int k=0; k<NumFacets; k++){
+ y = double(k) * off - 1. + (off / 2.);
+ r = pow(1. - y*y,0.5);
+ phi = double(k) * inc;
+ v.push_back(Vector3r(cos(phi)*r*radii[0], y*radii[1], sin(phi)*r*radii[2]));
+ }
+ }
+
+ // randomly rotate
+ srand(seed);
+ Quaternionr Rot(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
+ Rot.normalize();
+ for(int i=0; i< (int) v.size();i++) {
+ v[i] = Rot*(Vector3r(v[i][0],v[i][1],v[i][2]));
+ }
+ return v;
+}
+
+//**********************************************************************************
+//generate "packing" of non-overlapping balls
+vector<Vector3r> fillBoxByBalls_cpp(Vector3r minCoord, Vector3r maxCoord, Vector3r sizemin, Vector3r sizemax, Vector3r ratio, int seed, shared_ptr<Material> mat, int NumPoints){
+ vector<Vector3r> v;
+ Polyhedra trialP;
+ Polyhedron trial, trial_moved;
+ srand(seed);
+ int it = 0;
+ vector<Polyhedron> polyhedrons;
+ vector<vector<Vector3r> > vv;
+ Vector3r position;
+ bool intersection;
+ int count = 0;
+ Vector3r radii;
+
+
+ bool fixed_ratio = 0;
+ if (ratio[0] > 0 && ratio[1] > 0 && ratio[2]>0){
+ fixed_ratio = 1;
+ sizemax[0] = min(min(sizemax[0]/ratio[0], sizemax[1]/ratio[1]), sizemax[2]/ratio[2]);
+ sizemin[0] = max(max(sizemin[0]/ratio[0], sizemin[1]/ratio[1]), sizemin[2]/ratio[2]);
+ }
+
+ fixed_ratio = 1; //force spherical
+
+ //it - number of trials to make packing possibly more/less dense
+ Vector3r random_size;
+ while (it<1000){
+ it = it+1;
+ if (it == 1){
+ if (fixed_ratio) {
+ double rrr = (rand()*(sizemax[0]-sizemin[0])/RAND_MAX + sizemin[0])/2.;
+ radii = Vector3r(rrr,rrr,rrr);
+ }else {
+ radii = Vector3r(rand()*(sizemax[0]-sizemin[0])/2.,rand()*(sizemax[1]-sizemin[1])/2.,rand()*(sizemax[2]-sizemin[2])/2.)/RAND_MAX + sizemin/2.;
+ }
+ trialP.v = BallPoints(radii,NumPoints,rand());
+ trialP.Initialize();
+ trial = trialP.GetPolyhedron();
+ Matrix3r rot_mat = (trialP.GetOri()).toRotationMatrix();
+ Transformation t_rot(rot_mat(0,0),rot_mat(0,1),rot_mat(0,2),rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),1.);
+ std::transform( trial.points_begin(), trial.points_end(), trial.points_begin(), t_rot);
+ }
+ position = Vector3r(rand()*(maxCoord[0]-minCoord[0]),rand()*(maxCoord[1]-minCoord[1]),rand()*(maxCoord[2]-minCoord[2]))/RAND_MAX + minCoord;
+
+ //move CGAL structure Polyhedron
+ Transformation transl(CGAL::TRANSLATION, ToCGALVector(position));
+ trial_moved = trial;
+ std::transform( trial_moved.points_begin(), trial_moved.points_end(), trial_moved.points_begin(), transl);
+ //calculate plane equations
+ std::transform( trial_moved.facets_begin(), trial_moved.facets_end(), trial_moved.planes_begin(),Plane_equation());
+
+ intersection = false;
+ //call test with boundary
+ for(Polyhedron::Vertex_iterator vi = trial_moved.vertices_begin(); (vi != trial_moved.vertices_end()) && (!intersection); vi++){
+ intersection = (vi->point().x()<minCoord[0]) || (vi->point().x()>maxCoord[0]) || (vi->point().y()<minCoord[1]) || (vi->point().y()>maxCoord[1]) || (vi->point().z()<minCoord[2]) || (vi->point().z()>maxCoord[2]);
+ }
+ //call test with other polyhedrons
+ for(vector<Polyhedron>::iterator a = polyhedrons.begin(); (a != polyhedrons.end()) && (!intersection); a++){
+ intersection = do_intersect(*a,trial_moved);
+ if (intersection) break;
+ }
+ if (!intersection){
+ polyhedrons.push_back(trial_moved);
+ v.clear();
+ for(Polyhedron::Vertex_iterator vi = trial_moved.vertices_begin(); vi != trial_moved.vertices_end(); vi++){
+ v.push_back(FromCGALPoint(vi->point()));
+ }
+ vv.push_back(v);
+ it = 0;
+ count ++;
+
+ }
+ }
+ cout << "generated " << count << " polyhedrons"<< endl;
+
+ //can't be used - no information about material
+ Scene* scene=Omega::instance().getScene().get();
+ for(vector<vector<Vector3r> >::iterator p=vv.begin(); p!=vv.end(); ++p){
+ shared_ptr<Body> BP = NewPolyhedra(*p, mat);
+ BP->shape->color = Vector3r(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
+ scene->bodies->insert(BP);
+ }
+ return v;
+}
+
//**********************************************************************************
//split polyhedra
-void Split(const shared_ptr<Body> body, Vector3r direction){
- SplitPolyhedra(body, direction);
+void Split(const shared_ptr<Body> body, Vector3r direction, Vector3r point){
+ SplitPolyhedra(body, direction, point);
}
//**********************************************************************************
@@ -361,14 +532,15 @@
py::def("PWaveTimeStep",PWaveTimeStep,"Get timestep accoring to the velocity of P-Wave propagation; computed from sphere radii, rigidities and masses.");
py::def("do_Polyhedras_Intersect",do_Polyhedras_Intersect,"check polyhedras intersection");
py::def("fillBox_cpp",fillBox_cpp,"Generate non-overlaping polyhedrons in box");
+ py::def("fillBoxByBalls_cpp",fillBoxByBalls_cpp,"Generate non-overlaping 'spherical' polyhedrons in box");
py::def("MinCoord",MinCoord,"returns min coordinates");
py::def("MaxCoord",MaxCoord,"returns max coordinates");
py::def("SieveSize",SieveSize,"returns approximate sieve size of polyhedron");
py::def("SieveCurve",SieveCurve,"save sieve curve coordinates into file");
py::def("SizeOfPolyhedra",SizeOfPolyhedra,"returns max, middle an min size in perpendicular directions");
py::def("SizeRatio",SizeRatio,"save sizes of polyhedra into file");
- py::def("convexHull",convexHull,"");
- py::def("Split",Split,"split polyhedron perpendicularly to given direction direction");
+ py::def("convexHull",convexHull,"....");
+ py::def("Split",Split,"split polyhedron perpendicularly to given direction through given point");
}
#endif // YADE_CGAL
=== modified file 'py/polyhedra_utils.py'
--- py/polyhedra_utils.py 2013-10-15 16:14:46 +0000
+++ py/polyhedra_utils.py 2014-05-16 11:58:59 +0000
@@ -79,8 +79,54 @@
ball = polyhedra(material,v=pts)
ball.state.pos = center
return ball
+
+#**********************************************************************************
+def polyhedraTruncIcosaHed(radius, material, centre,mask=1):
+ pts = []
+
+ p = (1.+math.sqrt(5.))/2.
+ f = radius/math.sqrt(9.*p + 1.)
+ A = [[0.,1.,3.*p],[2.,1.+2.*p,p],[1.,2.+p,2.*p]]
+ for a in A:
+ a = [a[0]*f,a[1]*f,a[2]*f]
+ B = [a,[a[1],a[2],a[0]],[a[2],a[0],a[1]]]
+ for b in B:
+ pts.append(b)
+ if not b[0] == 0:
+ pts.append([-b[0], b[1], b[2]])
+ if not b[1] == 0:
+ pts.append([-b[0],-b[1], b[2]])
+ if not b[2] == 0: pts.append([-b[0],-b[1],-b[2]])
+ if not b[2] == 0:
+ pts.append([-b[0], b[1],-b[2]])
+ if not b[1] == 0:
+ pts.append([ b[0],-b[1], b[2]])
+ if not b[2] == 0:
+ pts.append([ b[0],-b[1],-b[2]])
+ if not b[2] == 0: pts.append([ b[0], b[1],-b[2]])
+ ball = polyhedra(material,v=pts)
+ ball.state.pos = centre
+ return ball
#**********************************************************************************
+def polyhedraSnubCube(radius, material, centre, mask=1):
+ pts = []
+
+ f = radius/1.3437133737446
+ c1 = 0.337754
+ c2 = 1.14261
+ c3 = 0.621226
+ A = [[c2,c1,c3],[c1,c3,c2],[c3,c2,c1],[-c1,-c2,-c3],[-c2,-c3,-c1],[-c3,-c1,-c2]]
+ for a in A:
+ a = [a[0]*f,a[1]*f,a[2]*f]
+ pts.append([-a[0],-a[1], a[2]])
+ pts.append([ a[0],-a[1],-a[2]])
+ pts.append([-a[0], a[1],-a[2]])
+ pts.append([ a[0], a[1], a[2]])
+ ball = polyhedra(material,v=pts)
+ ball.state.pos = centre
+ return ball
+#**********************************************************************************
#fill box [mincoord, maxcoord] by non-overlaping polyhedrons with random geometry and sizes within the range (uniformly distributed)
def fillBox(mincoord, maxcoord,material,sizemin=[1,1,1],sizemax=[1,1,1],ratio=[0,0,0],seed=None,mask=1):
"""fill box [mincoord, maxcoord] by non-overlaping polyhedrons with random geometry and sizes within the range (uniformly distributed)
@@ -98,5 +144,16 @@
# if(math.isnan(v[i][0])):
# O.bodies.append(polyhedra(material,seed=random.randint(0,1E6),v=v[lastnan+1:i],mask=1,fixed=False))
# lastnan = i
+
+#**********************************************************************************
+#fill box [mincoord, maxcoord] by non-overlaping polyhedrons with random geometry and sizes within the range (uniformly distributed)
+def fillBoxByBalls(mincoord, maxcoord,material,sizemin=[1,1,1],sizemax=[1,1,1],ratio=[0,0,0],seed=None,mask=1,numpoints=60):
+ random.seed(seed);
+ v = fillBoxByBalls_cpp(mincoord, maxcoord, sizemin,sizemax, ratio, random.randint(0,1E6), material,numpoints)
+ #lastnan = -1
+ #for i in range(0,len(v)):
+ # if(math.isnan(v[i][0])):
+ # O.bodies.append(polyhedra(material,seed=random.randint(0,1E6),v=v[lastnan+1:i],mask=1,fixed=False))
+ # lastnan = i
