yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10138
[Branch ~yade-pkg/yade/git-trunk] Rev 3715: Polyhedra implementation + examples (Contributed by Jan Elias).
------------------------------------------------------------
revno: 3715
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Tue 2013-10-15 18:14:46 +0200
message:
Polyhedra implementation + examples (Contributed by Jan Elias).
added:
examples/polyhedra/
examples/polyhedra/ball.py
examples/polyhedra/free-fall.py
examples/polyhedra/irregular.py
examples/polyhedra/wall.py
pkg/dem/Polyhedra.cpp
pkg/dem/Polyhedra.hpp
pkg/dem/Polyhedra_Ig2.cpp
pkg/dem/Polyhedra_Ig2.hpp
pkg/dem/Polyhedra_splitter.cpp
pkg/dem/Polyhedra_splitter.hpp
pkg/dem/Polyhedra_support.cpp
py/_polyhedra_utils.cpp
py/polyhedra_utils.py
modified:
.gitignore
doc/sphinx/yadeSphinx.py
py/CMakeLists.txt
--
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 '.gitignore'
--- .gitignore 2013-10-01 07:37:13 +0000
+++ .gitignore 2013-10-15 16:14:46 +0000
@@ -25,6 +25,7 @@
doc/sphinx/yade.utils.rst
doc/sphinx/yade.wrapper.rst
doc/sphinx/yade.ymport.rst
+doc/sphinx/yade.polyhedra_utils.rst
.kdev4/
*.kdev4
*.pyc
=== modified file 'doc/sphinx/yadeSphinx.py'
--- doc/sphinx/yadeSphinx.py 2013-09-11 10:55:59 +0000
+++ doc/sphinx/yadeSphinx.py 2013-10-15 16:14:46 +0000
@@ -35,7 +35,7 @@
#
# don't forget to put the module in index.rst as well!
#
-mods={'export':[],'post2d':[],'pack':['_packSpheres','_packPredicates','_packObb'],'plot':[],'timing':[],'utils':['_utils'],'ymport':[],'geom':[],'bodiesHandling':[],'qt':['qt._GLViewer'],'linterpolation':[]}
+mods={'export':[],'post2d':[],'pack':['_packSpheres','_packPredicates','_packObb'],'plot':[],'timing':[],'utils':['_utils'],'polyhedra_utils':['_polyhedra_utils'],'ymport':[],'geom':[],'bodiesHandling':[],'qt':['qt._GLViewer'],'linterpolation':[]}
#
# generate documentation, in alphabetical order
mm=mods.keys(); mm.sort()
=== added directory 'examples/polyhedra'
=== added file 'examples/polyhedra/ball.py'
--- examples/polyhedra/ball.py 1970-01-01 00:00:00 +0000
+++ examples/polyhedra/ball.py 2013-10-15 16:14:46 +0000
@@ -0,0 +1,47 @@
+from yade import plot, polyhedra_utils
+from yade import qt
+
+m = PolyhedraMat()
+m.density = 2600 #kg/m^3
+m.Ks = 20000
+m.Kn = 1E6 #Pa
+m.frictionAngle = 0.6 #rad
+
+maxLoad = 3E6
+minLoad = 1E3
+
+O.bodies.append(utils.wall(0,axis=2,sense=1, material = m))
+
+t = polyhedra_utils.polyhedralBall(0.025, 100, m, (0,0,0))
+t.state.pos = (0,0,0.5)
+O.bodies.append(t)
+
+def checkUnbalanced():
+ print "unbalanced forces = %.5f, position %f, %f, %f"%(utils.unbalancedForce(), t.state.pos[0], t.state.pos[1], t.state.pos[2])
+
+
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
+ InteractionLoop(
+ [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
+ [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
+ [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ ),
+ #GravityEngine(gravity=(0,0,-9.81)),
+ NewtonIntegrator(damping=0.5,gravity=(0,0,-9.81)),
+ PyRunner(command='checkUnbalanced()',realPeriod=3,label='checker')
+]
+
+
+#O.dt=0.025*polyhedra_utils.PWaveTimeStep()
+O.dt=0.00025
+
+
+qt.Controller()
+V = qt.View()
+
+
+O.saveTmp()
+#O.run()
=== added file 'examples/polyhedra/free-fall.py'
--- examples/polyhedra/free-fall.py 1970-01-01 00:00:00 +0000
+++ examples/polyhedra/free-fall.py 2013-10-15 16:14:46 +0000
@@ -0,0 +1,58 @@
+from yade import plot, polyhedra_utils
+
+
+gravel = PolyhedraMat()
+gravel.density = 2600 #kg/m^3
+gravel.Ks = 20000
+gravel.Kn = 1E7 #Pa
+gravel.frictionAngle = 0.5 #rad
+
+steel = PolyhedraMat()
+steel.density = 7850 #kg/m^3
+steel.Ks = 10*gravel.Ks
+steel.Kn = 10*gravel.Kn
+steel.frictionAngle = 0.4 #rad
+
+rubber = PolyhedraMat()
+rubber.density = 1000 #kg/m^3
+rubber.Ks = gravel.Ks/10
+rubber.Kn = gravel.Kn/10
+rubber.frictionAngle = 0.7 #rad
+
+O.bodies.append(polyhedra_utils.polyhedra(gravel,v=((0,0,-0.05),(0.3,0,-0.05),(0.3,0.3,-0.05),(0,0.3,-0.05),(0,0,0),(0.3,0,0),(0.3,0.3,0),(0,0.3,0)),fixed=True, color=(0.35,0.35,0.35)))
+#O.bodies.append(utils.wall(0,axis=1,sense=1, material = gravel))
+#O.bodies.append(utils.wall(0,axis=0,sense=1, material = gravel))
+#O.bodies.append(utils.wall(0.3,axis=1,sense=-1, material = gravel))
+#O.bodies.append(utils.wall(0.3,axis=0,sense=-1, material = gravel))
+
+polyhedra_utils.fillBox((0,0,0), (0.3,0.3,0.3),gravel,sizemin=[0.025,0.025,0.025],sizemax=[0.05,0.05,0.05],seed=4)
+
+def checkUnbalancedI():
+ print "iter %d, time elapsed %f, time step %.5e, unbalanced forces = %.5f"%(O.iter, O.realtime, O.dt, utils.unbalancedForce())
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
+ InteractionLoop(
+ [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
+ [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
+ [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ ),
+ #GravityEngine(gravity=(0,0,-9.81)),
+ NewtonIntegrator(damping=0.3,gravity=(0,0,-9.81)),
+ PyRunner(command='checkUnbalancedI()',realPeriod=5,label='checker')
+]
+
+
+#O.dt=0.25*polyhedra_utils.PWaveTimeStep()
+O.dt=0.0025*polyhedra_utils.PWaveTimeStep()
+
+
+from yade import qt
+qt.Controller()
+V = qt.View()
+V.screenSize = (550,450)
+V.sceneRadius = 1
+V.eyePosition = (0.7,0.5,0.1)
+V.upVector = (0,0,1)
+V.lookAt = (0.15,0.15,0.1)
=== added file 'examples/polyhedra/irregular.py'
--- examples/polyhedra/irregular.py 1970-01-01 00:00:00 +0000
+++ examples/polyhedra/irregular.py 2013-10-15 16:14:46 +0000
@@ -0,0 +1,63 @@
+# gravity deposition, continuing with oedometric test after stabilization
+# shows also how to run parametric studies with yade-batch
+
+# The components of the batch are:
+# 1. table with parameters, one set of parameters per line (ccc.table)
+# 2. utils.readParamsFromTable which reads respective line from the parameter file
+# 3. the simulation muse be run using yade-batch, not yade
+#
+# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
+#
+
+
+# create box with free top, and ceate loose packing inside the box
+from yade import plot, polyhedra_utils
+from yade import qt
+
+m = PolyhedraMat()
+m.density = 2600 #kg/m^3
+m.Ks = 20000
+m.Kn = 1E6 #Pa
+m.frictionAngle = 0.6 #rad
+
+O.bodies.append(utils.wall(0,axis=2,sense=1, material = m))
+
+t = polyhedra_utils.polyhedra(m,size = (0.06,0.06,0.06),seed = 5)
+t.state.pos = (0.,0.,0.5)
+O.bodies.append(t)
+
+def checkUnbalanced():
+ # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
+ print "unbalanced forces = %.5f, position %f, %f, %f"%(utils.unbalancedForce(), t.state.pos[0], t.state.pos[1], t.state.pos[2])
+
+
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
+ InteractionLoop(
+ [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
+ [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
+ [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ ),
+ #GravityEngine(gravity=(0,0,-9.81)),
+ NewtonIntegrator(damping=0.5,gravity=(0,0,-9.81)),
+ PyRunner(command='checkUnbalanced()',realPeriod=3,label='checker')
+
+]
+
+
+#O.dt=0.025*polyhedra_utils.PWaveTimeStep()
+O.dt=0.00025
+
+
+
+
+qt.Controller()
+V = qt.View()
+
+
+O.saveTmp()
+#O.run()
+#O.save('./done')
+utils.waitIfBatch()
=== added file 'examples/polyhedra/wall.py'
--- examples/polyhedra/wall.py 1970-01-01 00:00:00 +0000
+++ examples/polyhedra/wall.py 2013-10-15 16:14:46 +0000
@@ -0,0 +1,72 @@
+from yade import qt
+import numpy as np
+import random
+from yade import polyhedra_utils
+
+m = PolyhedraMat()
+m.density = 1000
+m.Ks = 5E6
+m.Kn = 5E8
+m.frictionAngle = 0.7
+
+size = 0.1;
+vertices = [[0,0,0],[size,0,0],[size,size,0],[size,size,size],[0,size,0],[0,size,size],[0,0,size],[size,0,size]]
+
+for i in range(0,10):
+ for j in range(0,10):
+ t = polyhedra_utils.polyhedra(m,v=vertices)
+ t.state.pos = (0,(i+0.5)*size-5*size,(j+0.5)*size*1)
+ O.bodies.append(t)
+
+
+
+qt.Controller()
+V = qt.View()
+V.screenSize = (750,550)
+V.sceneRadius = 1
+V.eyePosition = (-1.5,1.5,0.5)
+V.lookAt = (0,-0.5,0)
+V.upVector = (0,0,1)
+
+
+O.bodies.append(utils.wall(0,axis=2,sense=1,material=m))
+
+def checkUnbalanced():
+ print ('iter %d, time elapsed %f, unbalanced forces = %f'%(O.iter, O.realtime, utils.unbalancedForce()))
+ if O.iter<200: return
+ if utils.unbalancedForce()>0.01: return
+ AddBall()
+
+def AddBall():
+ ball = polyhedra_utils.polyhedralBall(size*1, 40, m, (-size*2.5,0,size*7),mask=1)
+ ball.shape.color = (0,0,0.9)
+ ball.state.vel = (15,0,0)
+ O.bodies.append(ball)
+ checker.dead = True
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb()]),
+ InteractionLoop(
+ [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
+ [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
+ [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ ),
+ #GravityEngine(gravity=(0,0,-9.81)),
+ NewtonIntegrator(damping=0.5,gravity=(0,0,-9.81)),
+ PyRunner(command='checkUnbalanced()',realPeriod=3,label='checker')#,
+ # wideo_recording
+ #qt.SnapshotEngine(fileBase='W',iterPeriod=50,label='snapshot')
+]
+
+#O.dt=.5*polyhedra_utils.PWaveTimeStep()
+O.dt = 0.000025
+O.saveTmp()
+
+#O.run()
+
+
+#comment ball to see stability
+
+
+
=== added file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra.cpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,427 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+#ifdef YADE_CGAL
+
+#include<yade/core/Interaction.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/core/State.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/common/Aabb.hpp>
+#include"Polyhedra.hpp"
+
+#define _USE_MATH_DEFINES
+
+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)
+ #endif
+ );
+
+//*********************************************************************************
+/* Polyhedra Constructor */
+
+void Polyhedra::Initialize(){
+
+ if (init) return;
+
+ bool isRandom = false;
+
+ //get vertices
+ int N = (int) v.size();
+ if (N==0) {
+ //generate randomly
+ while ((int) v.size()<4) GenerateRandomGeometry();
+ N = (int) v.size();
+ isRandom = true;
+ }
+
+ //compute convex hull of vertices
+ std::vector<CGALpoint> points;
+ points.resize(v.size());
+ for(int i=0;i<N;i++) {
+ points[i] = CGALpoint(v[i][0],v[i][1],v[i][2]);
+ }
+
+ CGAL::convex_hull_3(points.begin(), points.end(), P);
+
+ //connect triagular facets if possible
+ std::transform(P.facets_begin(), P.facets_end(), P.planes_begin(),Plane_equation());
+ P = Simplify(P, 1E-9);
+
+ //modify order of v according to CGAl polyhedron
+ int i = 0;
+ v.clear();
+ for (Polyhedron::Vertex_iterator vIter = P.vertices_begin(); vIter != P.vertices_end(); ++vIter, i++){
+ v.push_back(Vector3r(vIter->point().x(),vIter->point().y(),vIter->point().z()));
+ }
+
+ //list surface triangles for plotting
+ faceTri.clear();
+ std::transform(P.facets_begin(), P.facets_end(), P.planes_begin(),Plane_equation());
+ for (Polyhedron::Facet_iterator fIter = P.facets_begin(); fIter != P.facets_end(); fIter++){
+ Polyhedron::Halfedge_around_facet_circulator hfc0;
+ int n = fIter->facet_degree();
+ hfc0 = fIter->facet_begin();
+ int a = std::distance(P.vertices_begin(), hfc0->vertex());
+ for (int i=2; i<n; i++){
+ ++hfc0;
+ faceTri.push_back(a);
+ faceTri.push_back(std::distance(P.vertices_begin(), hfc0->vertex()));
+ faceTri.push_back(std::distance(P.vertices_begin(), hfc0->next()->vertex()));
+ }
+ }
+
+ //compute centroid and volume
+ P_volume_centroid(P, &volume, ¢roid);
+ //check vierd behavior of CGAL in tessalation
+ if(isRandom && volume*1.75<4./3.*3.14*size[0]/2.*size[1]/2.*size[2]/2.) {
+ v.clear();
+ seed = rand();
+ Initialize();
+ }
+ Vector3r translation((-1)*centroid);
+
+ //set centroid to be [0,0,0]
+ for(int i=0;i<N;i++) {
+ v[i] = v[i]-centroid;
+ }
+ if(isRandom) centroid = Vector3r::Zero();
+
+ Vector3r origin(0,0,0);
+
+ //move and rotate also the CGAL structure Polyhedron
+ Transformation t_trans(1.,0.,0.,translation[0],0.,1.,0.,translation[1],0.,0.,1.,translation[2],1.);
+ std::transform( P.points_begin(), P.points_end(), P.points_begin(), t_trans);
+
+ //compute inertia
+ double vtet;
+ Vector3r ctet;
+ Matrix3r Itet1, Itet2;
+ Matrix3r inertia_tensor(Matrix3r::Zero());
+ for(int i=0; i<(int) faceTri.size(); i+=3){
+ vtet = std::abs((origin-v[faceTri[i+2]]).dot((v[faceTri[i]]-v[faceTri[i+2]]).cross(v[faceTri[i+1]]-v[faceTri[i+2]]))/6.);
+ ctet = (origin+v[faceTri[i]]+v[faceTri[i+1]]+v[faceTri[i+2]]) / 4.;
+ Itet1 = TetraInertiaTensor(origin-ctet, v[faceTri[i]]-ctet, v[faceTri[i+1]]-ctet, v[faceTri[i+2]]-ctet);
+ ctet = ctet-origin;
+ Itet2<<
+ ctet[1]*ctet[1]+ctet[2]*ctet[2], -ctet[0]*ctet[1], -ctet[0]*ctet[2],
+ -ctet[0]*ctet[1], ctet[0]*ctet[0]+ctet[2]*ctet[2], -ctet[2]*ctet[1],
+ -ctet[0]*ctet[2], -ctet[2]*ctet[1], ctet[1]*ctet[1]+ctet[0]*ctet[0];
+ inertia_tensor = inertia_tensor + Itet1 + Itet2*vtet;
+ }
+
+ if(abs(inertia_tensor(0,1))+abs(inertia_tensor(0,2))+abs(inertia_tensor(1,2)) < 1E-13){
+ // no need to rotate, inertia already diagonal
+ orientation = Quaternionr::Identity();
+ inertia = Vector3r(inertia_tensor(0,0),inertia_tensor(1,1),inertia_tensor(2,2));
+ }else{
+ // calculate eigenvectors of I
+ Vector3r rot;
+ Matrix3r I_rot(Matrix3r::Zero()), I_new(Matrix3r::Zero());
+ matrixEigenDecomposition(inertia_tensor,I_rot,I_new);
+ // I_rot = eigenvectors of inertia_tensors in columns
+ // I_new = eigenvalues on diagonal
+ // set positove direction of vectors - otherwise reloading does not work
+ Matrix3r sign(Matrix3r::Zero());
+ double max_v_signed = I_rot(0,0);
+ double max_v = abs(I_rot(0,0));
+ if (max_v < abs(I_rot(1,0))) {max_v_signed = I_rot(1,0); max_v = abs(I_rot(1,0));}
+ if (max_v < abs(I_rot(2,0))) {max_v_signed = I_rot(2,0); max_v = abs(I_rot(2,0));}
+ sign(0,0) = max_v_signed/max_v;
+ max_v_signed = I_rot(0,1);
+ max_v = abs(I_rot(0,1));
+ if (max_v < abs(I_rot(1,1))) {max_v_signed = I_rot(1,1); max_v = abs(I_rot(1,1));}
+ if (max_v < abs(I_rot(2,1))) {max_v_signed = I_rot(2,1); max_v = abs(I_rot(2,1));}
+ sign(1,1) = max_v_signed/max_v;
+ sign(2,2) = 1.;
+ I_rot = I_rot*sign;
+ // force the eigenvectors to be right-hand oriented
+ Vector3r third = (I_rot.col(0)).cross(I_rot.col(1));
+ I_rot(0,2) = third[0];
+ I_rot(1,2) = third[1];
+ I_rot(2,2) = third[2];
+
+
+ inertia = Vector3r(I_new(0,0),I_new(1,1),I_new(2,2));
+ orientation = Quaternionr(I_rot);
+ //rotate the voronoi cell so that x - is maximal inertia axis and z - is minimal inertia axis
+ //orientation.normalize(); //not needed
+ for(int i=0; i< (int) v.size();i++) {
+ v[i] = orientation.conjugate()*v[i];
+ }
+
+ //rotate also the CGAL structure Polyhedron
+ Matrix3r rot_mat = (orientation.conjugate()).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( P.points_begin(), P.points_end(), P.points_begin(), t_rot);
+
+ }
+ //initialization done
+ init = 1;
+}
+
+//**************************************************************************
+/* Generator of randomly shaped polyhedron based on Voronoi tessellation*/
+
+void Polyhedra::GenerateRandomGeometry(){
+ srand(seed);
+
+ vector<CGALpoint> nuclei;
+ nuclei.push_back(CGALpoint(5.,5.,5.));
+ CGALpoint trial;
+ int iter = 0;
+ bool isOK;
+ //fill box 5x5x5 with randomly located nuclei with restricted minimal mutual distance 0.75 which gives approximate mean mutual distance 1;
+ double dist_min2 = 0.75*0.75;
+ while(iter<500){
+ isOK = true;
+ iter++;
+ trial = CGALpoint(double(rand())/RAND_MAX*5.+2.5,double(rand())/RAND_MAX*5.+2.5,double(rand())/RAND_MAX*5.+2.5);
+ for(int i=0;i< (int) nuclei.size();i++) {
+ isOK = pow(nuclei[i].x()-trial.x(),2)+pow(nuclei[i].y()-trial.y(),2)+pow(nuclei[i].z()-trial.z(),2) > dist_min2;
+ if (!isOK) break;
+ }
+
+ if(isOK){
+ iter = 0;
+ nuclei.push_back(trial);
+ }
+ }
+
+
+ //perform Voronoi tessellation
+ nuclei.erase(nuclei.begin());
+ Triangulation dt(nuclei.begin(), nuclei.end());
+ Triangulation::Vertex_handle zero_point = dt.insert(CGALpoint(5.,5.,5.));
+ v.clear();
+ std::vector<typename Triangulation::Cell_handle> ch_cells;
+ dt.incident_cells(zero_point,std::back_inserter(ch_cells));
+ 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.));
+ }
+
+
+ //resize and rotate the voronoi cell
+ 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]*size[0],v[i][1]*size[1],v[i][2]*size[2]));
+ }
+
+ //to avoid patological cases (that should not be present, but CGAL works somehow unpredicable)
+ if (v.size() < 8) {
+ cout << "wrong " << v.size() << endl;
+ v.clear();
+ seed = rand();
+ GenerateRandomGeometry();
+ }
+}
+
+//****************************************************************************************
+/* Destructor */
+
+Polyhedra::~Polyhedra(){}
+
+//****************************************************************************************
+/* Destructor */
+
+PolyhedraGeom::~PolyhedraGeom(){}
+
+//****************************************************************************************
+/* AaBb overlap checker */
+
+void Bo1_Polyhedra_Aabb::go(const shared_ptr<Shape>& ig, shared_ptr<Bound>& bv, const Se3r& se3, const Body*){
+
+ Polyhedra* t=static_cast<Polyhedra*>(ig.get());
+ if (!t->IsInitialized()) t->Initialize();
+ if(!bv){ bv=shared_ptr<Bound>(new Aabb); }
+ Aabb* aabb=static_cast<Aabb*>(bv.get());
+ //Quaternionr invRot=se3.orientation.conjugate();
+ int N = (int) t->v.size();
+ Vector3r v_g, mincoords(0.,0.,0.), maxcoords(0.,0.,0.);
+ for(int i=0; i<N; i++) {
+ v_g=se3.orientation*t->v[i]; // vertices in global coordinates
+ mincoords = Vector3r(min(mincoords[0],v_g[0]),min(mincoords[1],v_g[1]),min(mincoords[2],v_g[2]));
+ maxcoords = Vector3r(max(maxcoords[0],v_g[0]),max(maxcoords[1],v_g[1]),max(maxcoords[2],v_g[2]));
+ }
+ aabb->min=se3.position+mincoords;
+ aabb->max=se3.position+maxcoords;
+}
+
+//**********************************************************************************
+/* Plotting */
+
+#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&)
+ {
+ 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) {
+ glDisable(GL_LIGHTING);
+ glBegin(GL_LINES);
+ #define __ONEWIRE(a,b) glVertex3v(t->v[a]);glVertex3v(t->v[b])
+ for(int tri=0; tri < (int) faceTri.size(); tri+=3) {__ONEWIRE(faceTri[tri],faceTri[tri+1]); __ONEWIRE(faceTri[tri],faceTri[tri+2]); __ONEWIRE(faceTri[tri+1],faceTri[tri+2]);}
+ #undef __ONEWIRE
+ glEnd();
+ }
+ else
+ {
+ Vector3r centroid = t->GetCentroid();
+ Vector3r faceCenter, n;
+ glDisable(GL_CULL_FACE); glEnable(GL_LIGHTING);
+ glBegin(GL_TRIANGLES);
+ #define __ONEFACE(a,b,c) n=(t->v[b]-t->v[a]).cross(t->v[c]-t->v[a]); n.normalize(); faceCenter=(t->v[a]+t->v[b]+t->v[c])/3.; if((faceCenter-centroid).dot(n)<0)n=-n; glNormal3v(n); glVertex3v(t->v[a]); glVertex3v(t->v[b]); glVertex3v(t->v[c]);
+ for(int tri=0; tri < (int) faceTri.size(); tri+=3) {__ONEFACE(faceTri[tri],faceTri[tri+1],faceTri[tri+2]);}
+ #undef __ONEFACE
+ glEnd();
+ }
+ }
+#endif
+
+//**********************************************************************************
+//!Precompute data needed for rotating tangent vectors attached to the interaction
+
+void PolyhedraGeom::precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r&
+currentNormal, bool isNew, const Vector3r& shift2){
+
+ if(!isNew) {
+ orthonormal_axis = normal.cross(currentNormal);
+ Real angle = scene->dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
+ twist_axis = angle*normal;}
+ else twist_axis=orthonormal_axis=Vector3r::Zero();
+ //Update contact normal
+ normal=currentNormal;
+ //Precompute shear increment
+ Vector3r c1x = (contactPoint - rbp1.pos);
+ Vector3r c2x = (contactPoint - rbp2.pos + shift2);
+ Vector3r relativeVelocity = (rbp2.vel+rbp2.angVel.cross(c2x)) - (rbp1.vel+rbp1.angVel.cross(c1x));
+ //keep the shear part only
+ relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
+ shearInc = relativeVelocity*scene->dt;
+}
+
+//**********************************************************************************
+Vector3r& PolyhedraGeom::rotate(Vector3r& shearForce) const {
+ // approximated rotations
+ shearForce -= shearForce.cross(orthonormal_axis);
+ shearForce -= shearForce.cross(twist_axis);
+ //NOTE : make sure it is in the tangent plane? It's never been done before. Is it not adding rounding errors at the same time in fact?...
+ shearForce -= normal.dot(shearForce)*normal;
+ return shearForce;
+}
+
+
+//**********************************************************************************
+/* Material law, physics */
+
+void Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys::go( const shared_ptr<Material>& b1
+ , const shared_ptr<Material>& b2
+ , const shared_ptr<Interaction>& interaction)
+{
+ if(interaction->phys) return;
+ const shared_ptr<PolyhedraMat>& mat1 = YADE_PTR_CAST<PolyhedraMat>(b1);
+ const shared_ptr<PolyhedraMat>& mat2 = YADE_PTR_CAST<PolyhedraMat>(b2);
+ interaction->phys = shared_ptr<PolyhedraPhys>(new PolyhedraPhys());
+ const shared_ptr<PolyhedraPhys>& contactPhysics = YADE_PTR_CAST<PolyhedraPhys>(interaction->phys);
+ Real Kna = mat1->Kn;
+ Real Knb = mat2->Kn;
+ Real Ksa = mat1->Ks;
+ Real Ksb = mat2->Ks;
+ Real frictionAngle = std::min(mat1->frictionAngle,mat2->frictionAngle);
+ contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
+ contactPhysics->kn = Kna*Knb/(Kna+Knb);
+ contactPhysics->ks = Ksa*Ksb/(Ksa+Ksb);
+};
+
+//**************************************************************************************
+#if 1
+Real PolyhedraVolumetricLaw::getPlasticDissipation() {return (Real) plasticDissipation;}
+void PolyhedraVolumetricLaw::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
+Real PolyhedraVolumetricLaw::elasticEnergy()
+{
+ Real energy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ FrictPhys* phys = dynamic_cast<FrictPhys*>(I->phys.get());
+ if(phys) {
+ energy += 0.5*(phys->normalForce.squaredNorm()/phys->kn + phys->shearForce.squaredNorm()/phys->ks);}
+ }
+ return energy;
+}
+#endif
+
+
+//**************************************************************************************
+// 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;}
+ const shared_ptr<PolyhedraGeom>& contactGeom(dynamic_pointer_cast<PolyhedraGeom>(I->geom));
+ if(!contactGeom) {I->phys=shared_ptr<IPhys>(); return;}
+ const Body::id_t idA=I->getId1(), idB=I->getId2();
+ const shared_ptr<Body>& A=Body::byId(idA), B=Body::byId(idB);
+
+ PolyhedraPhys* phys = dynamic_cast<PolyhedraPhys*>(I->phys.get());
+
+ //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);
+ return;
+ }
+
+ //zero penetration depth means no interaction force
+ if(!(contactGeom->penetrationVolume > 0) ) {I->phys=shared_ptr<IPhys>(); 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
+ if (contactGeom->isShearNew)
+ shearForce = Vector3r::Zero();
+ else
+ shearForce = contactGeom->rotate(shearForce);
+
+ const Vector3r& shearDisp = contactGeom->shearInc;
+ shearForce -= phys->ks*shearDisp;
+
+ Real maxFs = normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
+
+ if (likely(!scene->trackEnergy && !traceEnergy)){//Update force but don't compute energy terms (see below))
+ // PFC3d SlipModel, is using friction angle. CoulombCriterion
+ if( shearForce.squaredNorm() > maxFs ){
+ Real ratio = sqrt(maxFs) / shearForce.norm();
+ shearForce *= ratio;}
+ } else {
+ //almost the same with additional Vector3r instatinated for energy tracing,
+ //duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
+ if(shearForce.squaredNorm() > maxFs){
+ Real ratio = sqrt(maxFs) / shearForce.norm();
+ Vector3r trialForce=shearForce;//store prev force for definition of plastic slip
+ //define the plastic work input and increment the total plastic energy dissipated
+ shearForce *= ratio;
+ Real dissip=((1/phys->ks)*(trialForce-shearForce)).dot(shearForce);
+ if (traceEnergy) plasticDissipation += dissip;
+ else if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,false);
+ }
+ // compute elastic energy as well
+ scene->energy->add(0.5*(normalForce.squaredNorm()/phys->kn+shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,true);
+ }
+
+ Vector3r F = -normalForce-shearForce;
+ if (contactGeom->equivalentPenetrationDepth != contactGeom->equivalentPenetrationDepth) exit(1);
+ scene->forces.addForce (idA,F);
+ 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));
+
+ //needed to be able to acces interaction forces in other parts of yade
+ phys->normalForce = normalForce;
+ phys->shearForce = shearForce;
+}
+
+#endif // YADE_CGAL
=== added file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra.hpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,274 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+
+#pragma once
+
+#ifdef YADE_CGAL
+
+#include<vector>
+#include<yade/core/Shape.hpp>
+#include<yade/core/IGeom.hpp>
+#include<yade/core/GlobalEngine.hpp>
+#include<yade/core/Material.hpp>
+#include<yade/pkg/common/Aabb.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
+#include<yade/pkg/common/Wall.hpp>
+#include<yade/pkg/common/Facet.hpp>
+#include<yade/lib/base/openmp-accu.hpp>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/Triangulation_data_structure_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Polyhedron_items_with_id_3.h>
+#include <CGAL/convex_hull_3.h>
+#include <CGAL/Tetrahedron_3.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include<time.h>
+
+#define likely(x) __builtin_expect((x),1)
+#define unlikely(x) __builtin_expect((x),0)
+
+//CGAL definitions - does not work with another kernel!! Why???
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Polyhedron_3<K> Polyhedron;
+typedef CGAL::Delaunay_triangulation_3<K> Triangulation;
+typedef K::Point_3 CGALpoint;
+typedef K::Vector_3 CGALvector;
+typedef CGAL::Aff_transformation_3<K> Transformation;
+typedef K::Segment_3 Segment;
+typedef CGAL::Triangle_3<K> Triangle;
+typedef CGAL::Plane_3<K> Plane;
+typedef CGAL::Line_3<K> Line;
+typedef CGAL::Origin CGAL_ORIGIN;
+
+//**********************************************************************************
+class Polyhedra: public Shape{
+ public:
+ //constructor from Vertices
+ Polyhedra(std::vector<Vector3r> V) { createIndex(); v.resize(V.size()); for(int i=0;i<(int) V.size();i++) v[i]=V[i]; Initialize();} //contructor of "random" polyhedra
+ Polyhedra(Vector3r xsize, int xseed) { createIndex(); seed=xseed; size=xsize; v.clear(); Initialize();}
+ virtual ~Polyhedra();
+ Vector3r GetCentroid(){Initialize(); return centroid;}
+ Vector3r GetInertia(){Initialize(); return inertia;}
+ vector<int> GetSurfaceTriangulation(){Initialize(); return faceTri;}
+ void Initialize();
+ bool IsInitialized(){return init;}
+ std::vector<Vector3r> GetOriginalVertices();
+ double GetVolume(){Initialize(); return volume;}
+ Quaternionr GetOri(){Initialize(); return orientation;}
+ Polyhedron GetPolyhedron(){return P;};
+ void Clear(){v.clear(); P.clear(); init = 0; size = Vector3r(1.,1.,1.); faceTri.clear();};
+
+ protected:
+ //triangulation of facets for plotting
+ vector<int> faceTri;
+ //centroid = (0,0,0) for random Polyhedra
+ Vector3r centroid;
+ //CGAL structure Polyhedron
+ Polyhedron P;
+ //sign of performed initialization
+ bool init;
+ //centroid Volume
+ double volume;
+ //centroid inerta - diagonal of the tensor
+ Vector3r inertia;
+ //orientation, that provides diagonal inertia tensor
+ Quaternionr orientation;
+ void GenerateRandomGeometry();
+
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(Polyhedra,Shape,"Polyhedral (convex) geometry.",
+ ((std::vector<Vector3r>,v,,,"Tetrahedron vertices in global coordinate system."))
+ ((int,seed, time(NULL),,"Seed for random generator."))
+ ((Vector3r, size, Vector3r(1.,1.,1.),,"Size of the grain in meters - x,y,z - before random rotation")),
+ /*init*/,
+ /*ctor*/
+ createIndex();
+ init = 0,
+ .def("Initialize",&Polyhedra::Initialize,"Initialization")
+ .def("GetVolume",&Polyhedra::GetVolume,"return polyhedra's volume")
+ .def("GetInertia",&Polyhedra::GetInertia,"return polyhedra's inertia tensor")
+ .def("GetOri",&Polyhedra::GetOri,"return polyhedra's orientation")
+ .def("GetCentroid",&Polyhedra::GetCentroid,"return polyhedra's centroid")
+ );
+ REGISTER_CLASS_INDEX(Polyhedra,Shape);
+};
+REGISTER_SERIALIZABLE(Polyhedra);
+
+
+//***************************************************************************
+/*! Collision configuration for Polyhedra and something.
+ * This is expressed as penetration volume properties: centroid, volume, depth ...
+ *
+ * Self-contained. */
+class PolyhedraGeom: public IGeom{
+ public:
+ virtual ~PolyhedraGeom();
+ //precompute data for shear evaluation
+ void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, const Vector3r& shift2);
+ Vector3r& rotate(Vector3r& shearForce) const;
+ //sep_plane is a code storing plane, that previously separated two polyhedras. It is used for faster detection of non-overlap.
+ std::vector<int> sep_plane;
+ bool isShearNew;
+ protected:
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraGeom,IGeom,"Geometry of interaction between 2 :yref:`vector<Polyhedra>`, including volumetric characteristics",
+ ((Real,penetrationVolume,NaN,,"Volume of overlap [m³]"))
+ ((Real,equivalentCrossSection,NaN,,"Cross-section area of the overlap (perpendicular to the normal) - not used"))
+ ((Real,equivalentPenetrationDepth,NaN,,"volume / equivalentCrossSection - not used"))
+ ((Vector3r,contactPoint,Vector3r::Zero(),,"Contact point (global coords), centriod of the overlapping polyhedron"))
+ ((Vector3r,shearInc,Vector3r::Zero(),,"Shear displacement increment in the last step"))
+ ((Vector3r,normal,Vector3r::Zero(),,"Normal direction of the interaction"))
+ ((Vector3r,twist_axis,Vector3r::Zero(),,""))
+ ((Vector3r,orthonormal_axis,Vector3r::Zero(),,"")),
+ createIndex();
+ sep_plane.assign(3,0);
+ );
+ //FUNCTOR2D(Tetra,Tetra);
+ REGISTER_CLASS_INDEX(PolyhedraGeom,IGeom);
+};
+REGISTER_SERIALIZABLE(PolyhedraGeom);
+
+//***************************************************************************
+/*! Creates Aabb from Polyhedra.
+ *
+ * Self-contained. */
+class Bo1_Polyhedra_Aabb: public BoundFunctor{
+ public:
+ void go(const shared_ptr<Shape>& ig, shared_ptr<Bound>& bv, const Se3r& se3, const Body*);
+ FUNCTOR1D(Polyhedra);
+ YADE_CLASS_BASE_DOC(Bo1_Polyhedra_Aabb,BoundFunctor,"Create/update :yref:`Aabb` of a :yref:`Polyhedra`");
+};
+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:
+ PolyhedraMat(double N, double S, double F){Kn=N; Ks=S; frictionAngle=F;};
+ double GetStrength(){return strength;};
+ virtual ~PolyhedraMat(){};
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraMat,Material,"Elastic material with Coulomb friction.",
+ ((Real,Kn,1e8,,"Normal volumetric 'stiffness' (N/m3)."))
+ ((Real,Ks,1e5,,"Shear stiffness (N/m)."))
+ ((Real,frictionAngle,.5,,"Contact friction angle (in radians)."))
+ ((bool,IsSplitable,0,,"To be splitted ... or not"))
+ ((double,strength,100,,"Stress at whis polyhedra of volume 4/3*pi [mm] breaks.")),
+ /*ctor*/ createIndex();
+ );
+ REGISTER_CLASS_INDEX(PolyhedraMat,Material);
+};
+REGISTER_SERIALIZABLE(PolyhedraMat);
+
+//***************************************************************************
+class PolyhedraPhys: public IPhys{
+ public:
+ virtual ~PolyhedraPhys(){};
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraPhys,IPhys,"Simple elastic material with friction for volumetric constitutive laws",
+ ((Real,kn,0,,"Normal stiffness"))
+ ((Vector3r,normalForce,Vector3r::Zero(),,"Normal force after previous step (in global coordinates)."))
+ ((Real,ks,0,,"Shear stiffness"))
+ ((Vector3r,shearForce,Vector3r::Zero(),,"Shear force after previous step (in global coordinates)."))
+ ((Real,tangensOfFrictionAngle,NaN,,"tangens of angle of internal friction")),
+ /*ctor*/ createIndex();
+ );
+ REGISTER_CLASS_INDEX(PolyhedraPhys,IPhys);
+};
+REGISTER_SERIALIZABLE(PolyhedraPhys);
+
+//***************************************************************************
+class Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys: public IPhysFunctor{
+ public:
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(PolyhedraMat,PolyhedraMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys,IPhysFunctor,"",
+ );
+};
+REGISTER_SERIALIZABLE(Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys);
+
+//***************************************************************************
+/*! Calculate physical response based on penetration configuration given by TTetraGeom. */
+
+class PolyhedraVolumetricLaw: public LawFunctor{
+ OpenMPAccumulator<Real> plasticDissipation;
+ virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ Real elasticEnergy ();
+ Real getPlasticDissipation();
+ void initPlasticDissipation(Real initVal=0);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(PolyhedraVolumetricLaw,LawFunctor,"Calculate physical response of 2 :yref:`vector<Polyhedra>` in interaction, based on penetration configuration given by :yref:`PolyhedraGeom`.",
+ ((Vector3r,shearForce,Vector3r::Zero(),,"Shear force from last step"))
+ ((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
+ ((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
+ ((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)"))
+ ,,
+ .def("elasticEnergy",&PolyhedraVolumetricLaw::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
+ .def("plasticDissipation",&PolyhedraVolumetricLaw::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
+ .def("initPlasticDissipation",&PolyhedraVolumetricLaw::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
+ );
+ FUNCTOR2D(PolyhedraGeom,PolyhedraPhys);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(PolyhedraVolumetricLaw);
+
+
+//***************************************************************************
+//compute plane equation from three points on the facet
+struct Plane_equation {
+ template <class Facet>
+ typename Facet::Plane_3 operator()( Facet& f) {
+ typename Facet::Halfedge_handle h = f.halfedge();
+ typedef typename Facet::Plane_3 Plane;
+ return Plane( h->vertex()->point(),
+ h->next()->vertex()->point(),
+ h->next()->next()->vertex()->point());
+ }
+};
+//get Tetrahedron inertia
+Matrix3r TetraInertiaTensor(Vector3r av,Vector3r bv,Vector3r cv,Vector3r dv);
+//return intersection of two polyhedrons
+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);
+//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
+bool P_volume_centroid(Polyhedron P, double * volume, Vector3r * centroid);
+//CGAL - miniEigen communication
+Vector3r FromCGALPoint(CGALpoint A);
+Vector3r FromCGALVector(CGALvector A);
+CGALpoint ToCGALPoint(Vector3r A);
+CGALvector ToCGALVector(Vector3r A);
+//determination of intersection of two polyhedras
+bool do_intersect(Polyhedron A, Polyhedron B);
+bool do_intersect(Polyhedron A, Polyhedron B, std::vector<int> &sep_plane);
+//connect triagular facets if possible
+Polyhedron Simplify(Polyhedron P, double lim);
+//list of facets and edges
+void PrintPolyhedron(Polyhedron P);
+//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);
+//new polyhedra
+shared_ptr<Body> NewPolyhedra(vector<Vector3r> v, shared_ptr<Material> mat);
+
+#endif // YADE_CGAL
=== added file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,305 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+#ifdef YADE_CGAL
+
+#include"Polyhedra.hpp"
+#include"Polyhedra_Ig2.hpp"
+
+#define _USE_MATH_DEFINES
+
+YADE_PLUGIN(/* self-contained in hpp: */ (Ig2_Polyhedra_Polyhedra_PolyhedraGeom) (Ig2_Wall_Polyhedra_PolyhedraGeom) (Ig2_Facet_Polyhedra_PolyhedraGeom) );
+
+//**********************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Polyhedras. */
+
+bool Ig2_Polyhedra_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& cm1,const shared_ptr<Shape>& cm2,const State& state1,const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& interaction){
+
+ //get polyhedras
+ const Se3r& se31=state1.se3;
+ const Se3r& se32=state2.se3;
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+ Polyhedra* B = static_cast<Polyhedra*>(cm2.get());
+
+ bool isNew = !interaction->geom;
+
+ //move and rotate 1st the CGAL structure Polyhedron
+ Matrix3r rot_mat = (se31.orientation).toRotationMatrix();
+ Vector3r trans_vec = se31.position;
+ 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);
+
+
+ //move and rotate 2nd the CGAL structure Polyhedron
+ rot_mat = (se32.orientation).toRotationMatrix();
+ trans_vec = se32.position;
+ 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);
+
+ shared_ptr<PolyhedraGeom> bang;
+ if (isNew) {
+ // new interaction
+ bang=shared_ptr<PolyhedraGeom>(new PolyhedraGeom());
+ bang->sep_plane.assign(3,0);
+ bang->contactPoint = Vector3r(0,0,0);
+ bang->isShearNew = true;
+ interaction->geom = bang;
+ }else{
+ // use data from old interaction
+ bang=YADE_PTR_CAST<PolyhedraGeom>(interaction->geom);
+ bang->isShearNew = bang->equivalentPenetrationDepth<=0;
+ }
+
+
+ //find intersection Polyhedra
+ Polyhedron Int;
+ Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position), bang->sep_plane);
+
+ //volume and centroid of intersection
+ double volume;
+ 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;}
+
+ //find normal direction
+ Vector3r normal = FindNormal(Int, PA, PB);
+ if((se32.position-centroid).dot(normal)<0) normal*=-1;
+
+ //calculate area of projection of Intersection into the normal plane
+ //double area = CalculateProjectionArea(Int, ToCGALVector(normal));
+ //if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
+ double area = volume/1E-8;
+
+ // store calculated stuff in bang; some is redundant
+ bang->equivalentCrossSection=area;
+ bang->contactPoint=centroid;
+ bang->penetrationVolume=volume;
+ bang->equivalentPenetrationDepth=volume/area;
+ bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
+ bang->normal=normal;
+
+ return true;
+}
+
+//**********************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Polyhedron and Wall. */
+
+bool Ig2_Wall_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& cm1,const shared_ptr<Shape>& cm2,const State& state1,const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& interaction){
+
+ const int& PA(cm1->cast<Wall>().axis);
+ const int& sense(cm1->cast<Wall>().sense);
+ const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
+ Polyhedra* B = static_cast<Polyhedra*>(cm2.get());
+
+ bool isNew = !interaction->geom;
+
+ //move and rotate also the CGAL structure Polyhedron
+ Matrix3r rot_mat = (se32.orientation).toRotationMatrix();
+ Vector3r trans_vec = se32.position;
+ 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 wall
+ Vector3r normal = Vector3r(0,0,0);
+ normal[PA] = sense;
+ CGALvector CGALnormal = CGALvector(normal[0],normal[1],normal[2]);
+ Plane A = Plane(CGALpoint(se31.position[0],se31.position[1],se31.position[2]),CGALvector(normal[0],normal[1],normal[2]));
+
+ shared_ptr<PolyhedraGeom> bang;
+ if (isNew) {
+ // new interaction
+ bang=shared_ptr<PolyhedraGeom>(new PolyhedraGeom());
+ bang->contactPoint = Vector3r(0,0,0);
+ bang->isShearNew = true;
+ interaction->geom = bang;
+ }else{
+ // use data from old interaction
+ bang=YADE_PTR_CAST<PolyhedraGeom>(interaction->geom);
+ bang->isShearNew = bang->equivalentPenetrationDepth<=0;
+ }
+
+ //find intersection Polyhedra
+ Polyhedron Int;
+ Int = Polyhedron_Plane_intersection(PB,A,ToCGALPoint(se32.position),ToCGALPoint(bang->contactPoint));
+
+ //volume and centroid of intersection
+ double volume;
+ Vector3r centroid;
+ P_volume_centroid(Int, &volume, ¢roid);
+ if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
+
+ //calculate area of projection of Intersection into the normal plane
+ double area = volume/1E-8;
+ //double area = CalculateProjectionArea(Int, CGALnormal);
+ //if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
+
+ // store calculated stuff in bang; some is redundant
+ bang->equivalentCrossSection=area;
+ bang->contactPoint=centroid;
+ bang->penetrationVolume=volume;
+ bang->equivalentPenetrationDepth=volume/area;
+ bang->precompute(state1,state2,scene,interaction,normal,isNew,shift2);
+ bang->normal=FromCGALVector(CGALnormal);
+
+ return true;
+}
+
+//**********************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Polyhedron and Facet. */
+
+bool Ig2_Facet_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& cm1,const shared_ptr<Shape>& cm2,const State& state1,const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& interaction){
+
+
+ const Se3r& se31=state1.se3;
+ const Se3r& se32=state2.se3;
+ Facet* A = static_cast<Facet*>(cm1.get());
+ Polyhedra* B = static_cast<Polyhedra*>(cm2.get());
+
+ bool isNew = !interaction->geom;
+
+ //move and rotate 1st the CGAL structure Polyhedron
+ Matrix3r rot_mat = (se32.orientation).toRotationMatrix();
+ Vector3r trans_vec = se32.position;
+ 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);
+
+ //move and rotate facet
+ vector<CGALpoint> v;
+ v.resize(6);
+ for (int i=0; i<3; i++) v[i] = ToCGALPoint(se31.orientation*A->vertices[i] + se31.position); // vertices in global coordinates
+
+ //determine
+ CGALpoint f_center((v[0].x()+v[1].x()+v[2].x())/3.,(v[0].y()+v[1].y()+v[2].y())/3.,(v[0].z()+v[1].z()+v[2].z())/3.);
+ CGALvector f_normal = CGAL::cross_product((v[1]-v[0]),(v[2]-v[0]));
+
+ //chage normal + change order of vertices to get outwarding facet normal for initial tetrahedron
+ if ((ToCGALPoint(se32.position)-f_center)*(f_normal) < 0) {
+ f_normal = (-1)*f_normal;
+ CGALpoint help;
+ help = v[0];
+ v[0]=v[1];
+ v[1]=help;
+ }
+
+ double f_area = sqrt(f_normal.squared_length());
+ for (int i=3; i<6; i++) v[i] = v[i-3]-f_normal/f_area*0.05*sqrt(f_area); // vertices in global coordinates
+
+ Polyhedron PA;
+ //use convex hull
+ //CGAL::convex_hull_3(v.begin(), v.end(), PA);
+
+ //construct polyhedron directly
+ Polyhedron::Halfedge_iterator hei = PA.make_tetrahedron(v[4],v[1],v[0],v[2]);
+ hei = PA.split_vertex(hei, hei->next()->opposite());
+ hei->vertex()->point() = v[3];
+ hei = PA.split_facet(hei->next()->next()->next(),hei->next());
+ hei = PA.split_vertex(hei, hei->next_on_vertex()->next_on_vertex());
+ hei->vertex()->point() = v[5];
+
+ shared_ptr<PolyhedraGeom> bang;
+ if (isNew) {
+ // new interaction
+ bang=shared_ptr<PolyhedraGeom>(new PolyhedraGeom());
+ bang->sep_plane.assign(3,0);
+ bang->contactPoint = Vector3r(0,0,0);
+ bang->isShearNew = true;
+ interaction->geom = bang;
+ }else{
+ // use data from old interaction
+ bang=YADE_PTR_CAST<PolyhedraGeom>(interaction->geom);
+ bang->isShearNew = bang->equivalentPenetrationDepth<=0;
+ }
+
+ //find intersection Polyhedra
+ Polyhedron Int;
+ Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position), bang->sep_plane);
+
+ //volume and centroid of intersection
+ double volume;
+ Vector3r centroid;
+ P_volume_centroid(Int, &volume, ¢roid);
+ if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
+
+ //find normal direction
+ Vector3r normal = FindNormal(Int, PA, PB);
+ if((se32.position-centroid).dot(normal)<0) normal*=-1;
+
+ //calculate area of projection of Intersection into the normal plane
+ double area = volume/1E-8;
+ //double area = CalculateProjectionArea(Int, ToCGALVector(normal));
+ //if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
+
+ // store calculated stuff in bang; some is redundant
+ bang->equivalentCrossSection=area;
+ bang->contactPoint=centroid;
+ bang->penetrationVolume=volume;
+ bang->equivalentPenetrationDepth=volume/area;
+ bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
+ bang->normal=normal;
+
+ return true;
+}
+
+//**********************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Polyhedron and Wall. */
+
+/*
+bool Ig2_Sphere_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& cm1,const shared_ptr<Shape>& cm2,const State& state1,const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& interaction){
+
+ const Se3r& se31=state1.se3;
+ const Se3r& se32=state2.se3;
+ Polyhedra* B = static_cast<Polyhedra*>(cm2.get());
+ const Real& r=cm1->cast<Sphere>().radius;
+
+ //cout << "Sphere x Polyhedra" << endl;
+
+ bool isNew = !interaction->geom;
+
+ //move and rotate 1st the CGAL structure Polyhedron
+ Matrix3r rot_mat = (se32.orientation).toRotationMatrix();
+ Vector3r trans_vec = se32.position;
+ 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());
+
+ shared_ptr<PolyhedraGeom> bang;
+ if (isNew) {
+ // new interaction
+ bang=shared_ptr<PolyhedraGeom>(new PolyhedraGeom());
+ bang->isShearNew = true;
+ interaction->geom = bang;
+ }else{
+ // use data from old interaction
+ bang=YADE_PTR_CAST<PolyhedraGeom>(interaction->geom);
+ }
+
+ //volume and centroid of intersection
+ double volume, area;
+ CGALvector normalCGAL;
+ CGALpoint centroidCGAL=ToCGALPoint(se32.position);
+ Sphere_Polyhedron_intersection(PB, r, ToCGALPoint(se31.position), centroidCGAL, volume, normalCGAL, area);
+ if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
+ Vector3r centroid = FromCGALPoint(centroidCGAL);
+ Vector3r normal = FromCGALVector(normalCGAL);
+
+
+
+ // store calculated stuff in bang; some is redundant
+ bang->equivalentCrossSection=area;
+ bang->contactPoint=centroid;
+ bang->penetrationVolume=volume;
+ bang->equivalentPenetrationDepth=volume/area;
+ bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
+ bang->normal=normal;
+
+ return true;
+}
+*/
+
+#endif // YADE_CGAL
=== added file 'pkg/dem/Polyhedra_Ig2.hpp'
--- pkg/dem/Polyhedra_Ig2.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra_Ig2.hpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,71 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+#ifdef YADE_CGAL
+
+#include"Polyhedra.hpp"
+#include<yade/pkg/common/Sphere.hpp>
+
+//***************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Polyhedras. */
+class Ig2_Polyhedra_Polyhedra_PolyhedraGeom: public IGeomFunctor
+{
+ public:
+ virtual ~Ig2_Polyhedra_Polyhedra_PolyhedraGeom(){};
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ FUNCTOR2D(Polyhedra,Polyhedra);
+ DEFINE_FUNCTOR_ORDER_2D(Polyhedra,Polyhedra);
+ YADE_CLASS_BASE_DOC(Ig2_Polyhedra_Polyhedra_PolyhedraGeom,IGeomFunctor,"Create/update geometry of collision between 2 Polyhedras");
+ DECLARE_LOGGER;
+ private:
+};
+REGISTER_SERIALIZABLE(Ig2_Polyhedra_Polyhedra_PolyhedraGeom);
+
+//***************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Wall & Polyhedra. */
+class Ig2_Wall_Polyhedra_PolyhedraGeom: public IGeomFunctor
+{
+ public:
+ virtual ~Ig2_Wall_Polyhedra_PolyhedraGeom(){};
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ FUNCTOR2D(Wall,Polyhedra);
+ DEFINE_FUNCTOR_ORDER_2D(Wall,Polyhedra);
+ YADE_CLASS_BASE_DOC(Ig2_Wall_Polyhedra_PolyhedraGeom,IGeomFunctor,"Create/update geometry of collision between Wall and Polyhedra");
+ DECLARE_LOGGER;
+ private:
+};
+REGISTER_SERIALIZABLE(Ig2_Wall_Polyhedra_PolyhedraGeom);
+
+//***************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Facet & Polyhedra. */
+class Ig2_Facet_Polyhedra_PolyhedraGeom: public IGeomFunctor
+{
+ public:
+ virtual ~Ig2_Facet_Polyhedra_PolyhedraGeom(){};
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ FUNCTOR2D(Facet,Polyhedra);
+ DEFINE_FUNCTOR_ORDER_2D(Facet,Polyhedra);
+ YADE_CLASS_BASE_DOC(Ig2_Facet_Polyhedra_PolyhedraGeom,IGeomFunctor,"Create/update geometry of collision between Facet and Polyhedra");
+ DECLARE_LOGGER;
+ private:
+};
+REGISTER_SERIALIZABLE(Ig2_Facet_Polyhedra_PolyhedraGeom);
+
+//***************************************************************************
+/*! Create Polyhedra (collision geometry) from colliding Sphere & Polyhedra. */
+/*
+class Ig2_Sphere_Polyhedra_PolyhedraGeom: public IGeomFunctor
+{
+ public:
+ virtual ~Ig2_Sphere_Polyhedra_PolyhedraGeom(){};
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ FUNCTOR2D(Sphere,Polyhedra);
+ DEFINE_FUNCTOR_ORDER_2D(Sphere,Polyhedra);
+ YADE_CLASS_BASE_DOC(Ig2_Sphere_Polyhedra_PolyhedraGeom,IGeomFunctor,"Create/update geometry of collision between Sphere and Polyhedra");
+ DECLARE_LOGGER;
+ private:
+};
+REGISTER_SERIALIZABLE(Ig2_Sphere_Polyhedra_PolyhedraGeom);
+*/
+
+#endif // YADE_CGAL
=== added file 'pkg/dem/Polyhedra_splitter.cpp'
--- pkg/dem/Polyhedra_splitter.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra_splitter.cpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,109 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+#ifdef YADE_CGAL
+
+#include<yade/pkg/dem/Polyhedra_splitter.hpp>
+
+
+YADE_PLUGIN((PolyhedraSplitter));
+CREATE_LOGGER(PolyhedraSplitter);
+
+//*********************************************************************************
+/* Evaluate tensorial stress estimation in polyhedras */
+
+void getStressForEachBody(vector<Matrix3r>& bStresses){
+ const shared_ptr<Scene>& scene=Omega::instance().getScene();
+ bStresses.resize(scene->bodies->size());
+ for (size_t k=0;k<scene->bodies->size();k++) bStresses[k]=Matrix3r::Zero();
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ PolyhedraGeom* geom=YADE_CAST<PolyhedraGeom*>(I->geom.get());
+ 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
+ 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());
+
+
+ }
+}
+
+//*********************************************************************************
+/* Size dependent strength */
+
+double PolyhedraSplitter::getStrength(double volume, double strength){
+ //equvalent radius
+ double r_eq = pow(volume*3./4./Mathr::PI,1./3.);
+ //r should be in milimeters
+ return strength/(r_eq/1000.);
+}
+
+//*********************************************************************************
+/* Symmetrization of stress tensor */
+
+void PolyhedraSplitter::Symmetrize(Matrix3r & bStress){
+ bStress(0,1) = (bStress(0,1) + bStress(1,0))/2.;
+ bStress(0,2) = (bStress(0,2) + bStress(2,0))/2.;
+ bStress(1,2) = (bStress(1,2) + bStress(2,1))/2.;
+ bStress(1,0) = bStress(0,1);
+ bStress(2,0) = bStress(0,2);
+ bStress(2,1) = bStress(1,2);
+}
+
+//*********************************************************************************
+/* Split if stress exceed strength */
+
+void PolyhedraSplitter::action()
+{
+ const shared_ptr<Scene> _rb=shared_ptr<Scene>();
+ shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
+
+ vector<shared_ptr<Body> > bodies;
+ vector<Vector3r > directions;
+ vector<double > sigmas;
+
+
+
+ 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];
+
+ //get eigenstresses
+ Symmetrize(stress);
+ Matrix3r I_vect(Matrix3r::Zero()), I_valu(Matrix3r::Zero());
+ matrixEigenDecomposition(stress,I_vect,I_valu);
+ int min_i = 0;
+ if (I_valu(min_i,min_i) > I_valu(1,1)) min_i = 1;
+ if (I_valu(min_i,min_i) > I_valu(2,2)) min_i = 2;
+ int max_i = 0;
+ if (I_valu(max_i,max_i) < I_valu(1,1)) max_i = 1;
+ if (I_valu(max_i,max_i) < I_valu(2,2)) max_i = 2;
+
+ //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();
+ 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);};
+ }
+ }
+ for(int i=0; i<int(bodies.size()); i++){
+ SplitPolyhedra(bodies[i], directions[i]);
+ }
+}
+
+#endif // YADE_CGAL
=== added file 'pkg/dem/Polyhedra_splitter.hpp'
--- pkg/dem/Polyhedra_splitter.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra_splitter.hpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,33 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+
+#pragma once
+
+#ifdef YADE_CGAL
+
+#include<yade/pkg/common/PeriodicEngines.hpp>
+#include<yade/core/Interaction.hpp>
+
+#include"Polyhedra.hpp"
+
+//*********************************************************************************
+/* Polyhedra Splitter */
+class PolyhedraSplitter : public PeriodicEngine{
+
+ public:
+ virtual void action();
+ double getStrength(double volume, double strength);
+ void Symmetrize(Matrix3r & bStress);
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(PolyhedraSplitter,PeriodicEngine,"Engine that splits polyhedras.",
+ ,
+ /*ctor*/
+ ,/*py*/
+ );
+ DECLARE_LOGGER;
+};
+
+REGISTER_SERIALIZABLE(PolyhedraSplitter);
+
+#endif // YADE_CGAL
=== added file 'pkg/dem/Polyhedra_support.cpp'
--- pkg/dem/Polyhedra_support.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Polyhedra_support.cpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,817 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+#ifdef YADE_CGAL
+
+#include"Polyhedra.hpp"
+
+#define _USE_MATH_DEFINES
+
+
+//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
+//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
+//FIND_NORMAL_LIMIT - to determine which facet of intersection belongs to which polyhedron
+#define FIND_NORMAL_LIMIT 1E-40
+
+
+//**********************************************************************************
+//return volume and centroid of polyhedron
+
+bool P_volume_centroid(Polyhedron P, double * volume, Vector3r * centroid){
+
+
+ Vector3r basepoint = FromCGALPoint(P.vertices_begin()->point());
+ Vector3r A,B,C,D;
+ (*volume) = 0;
+ double vtet;
+ (*centroid) = Vector3r(0.,0.,0.);
+
+ //compute centroid and volume
+ 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());
+ vtet = std::abs((basepoint-C).dot((A-C).cross(B-C)))/6.;
+ *volume += vtet;
+ *centroid = *centroid + (basepoint+A+B+C) / 4. * vtet;
+ }
+ }
+ *centroid = *centroid/(*volume);
+ return true;
+}
+
+//**********************************************************************************
+//STOLEN FROM TETRA body of Vaclav Smilauer
+
+/*! Calculates tetrahedron inertia relative to the origin (0,0,0), with unit density (scales linearly).
+
+See article F. Tonon, "Explicit Exact Formulas for the 3-D Tetrahedron Inertia Tensor in Terms of its Vertex Coordinates", http://www.scipub.org/fulltext/jms2/jms2118-11.pdf
+*/
+
+//Centroid MUST be [0,0,0]
+Matrix3r TetraInertiaTensor(Vector3r av,Vector3r bv,Vector3r cv,Vector3r dv){
+ #define x1 av[0]
+ #define y1 av[1]
+ #define z1 av[2]
+ #define x2 bv[0]
+ #define y2 bv[1]
+ #define z2 bv[2]
+ #define x3 cv[0]
+ #define y3 cv[1]
+ #define z3 cv[2]
+ #define x4 dv[0]
+ #define y4 dv[1]
+ #define z4 dv[2]
+
+ // Jacobian of transformation to the reference 4hedron
+ double detJ=(x2-x1)*(y3-y1)*(z4-z1)+(x3-x1)*(y4-y1)*(z2-z1)+(x4-x1)*(y2-y1)*(z3-z1)
+ -(x2-x1)*(y4-y1)*(z3-z1)-(x3-x1)*(y2-y1)*(z4-z1)-(x4-x1)*(y3-y1)*(z2-z1);
+ detJ=fabs(detJ);
+ double a=detJ*(y1*y1+y1*y2+y2*y2+y1*y3+y2*y3+
+ y3*y3+y1*y4+y2*y4+y3*y4+y4*y4+z1*z1+z1*z2+
+ z2*z2+z1*z3+z2*z3+z3*z3+z1*z4+z2*z4+z3*z4+z4*z4)/60.;
+ double b=detJ*(x1*x1+x1*x2+x2*x2+x1*x3+x2*x3+x3*x3+
+ x1*x4+x2*x4+x3*x4+x4*x4+z1*z1+z1*z2+z2*z2+z1*z3+
+ z2*z3+z3*z3+z1*z4+z2*z4+z3*z4+z4*z4)/60.;
+ double c=detJ*(x1*x1+x1*x2+x2*x2+x1*x3+x2*x3+x3*x3+x1*x4+
+ x2*x4+x3*x4+x4*x4+y1*y1+y1*y2+y2*y2+y1*y3+
+ y2*y3+y3*y3+y1*y4+y2*y4+y3*y4+y4*y4)/60.;
+ // a' in the article etc.
+ double a__=detJ*(2*y1*z1+y2*z1+y3*z1+y4*z1+y1*z2+
+ 2*y2*z2+y3*z2+y4*z2+y1*z3+y2*z3+2*y3*z3+
+ y4*z3+y1*z4+y2*z4+y3*z4+2*y4*z4)/120.;
+ double b__=detJ*(2*x1*z1+x2*z1+x3*z1+x4*z1+x1*z2+
+ 2*x2*z2+x3*z2+x4*z2+x1*z3+x2*z3+2*x3*z3+
+ x4*z3+x1*z4+x2*z4+x3*z4+2*x4*z4)/120.;
+ double c__=detJ*(2*x1*y1+x2*y1+x3*y1+x4*y1+x1*y2+
+ 2*x2*y2+x3*y2+x4*y2+x1*y3+x2*y3+2*x3*y3+
+ x4*y3+x1*y4+x2*y4+x3*y4+2*x4*y4)/120.;
+
+ Matrix3r ret; ret<<
+ a , -c__, -b__,
+ -c__, b , -a__,
+ -b__, -a__, c ;
+ return ret;
+
+ #undef x1
+ #undef y1
+ #undef z1
+ #undef x2
+ #undef y2
+ #undef z2
+ #undef x3
+ #undef y3
+ #undef z3
+ #undef x4
+ #undef y4
+ #undef z4
+}
+
+//**********************************************************************************
+//distace of point from a plane (squared) with sign
+double Oriented_squared_distance(Plane P, CGALpoint x){
+ double h = P.a()*x.x()+P.b()*x.y()+P.c()*x.z()+P.d();
+ return ((h>0.)-(h<0.))*pow(h,2)/(CGALvector(P.a(),P.b(),P.c())).squared_length();
+}
+
+//**********************************************************************************
+// test if point is inside polyhedra in strong sence, i.e. boundary location is not enough
+bool Is_inside_Polyhedron(Polyhedron P, CGALpoint inside){
+ Polyhedron::Plane_iterator pi;
+ for(pi=P.planes_begin(); pi!=P.planes_end(); ++pi){
+ if( ! pi->has_on_negative_side(inside)) return false;
+ }
+ return true;
+}
+
+//**********************************************************************************
+// test if point is inside polyhedra not closer than lim to its boundary
+bool Is_inside_Polyhedron(Polyhedron P, CGALpoint inside, double lim){
+ Polyhedron::Plane_iterator pi;
+ lim = pow(lim,2);
+ for(pi=P.planes_begin(); pi!=P.planes_end(); ++pi){
+ if(Oriented_squared_distance(*pi, inside)>-lim) return false;
+ }
+ return true;
+}
+
+//**********************************************************************************
+//test if two polyhedron intersect
+bool do_intersect(Polyhedron A, Polyhedron B){
+ std::vector<int> sep_plane;
+ sep_plane.assign(3,0);
+ return do_intersect(A, B, sep_plane);
+}
+
+//**********************************************************************************
+//test if two polyhedron intersect based on previous data
+bool do_intersect(Polyhedron A, Polyhedron B, std::vector<int> &sep_plane){
+
+
+ bool found;
+ //check previous separation plane
+ switch (sep_plane[0]){
+ case 1: //separation plane was previously determined as sep_plane[2]-th plane of A polyhedron
+ {
+ if(unlikely((unsigned) sep_plane[2]>=A.size_of_facets())) break;
+ Polyhedron::Facet_iterator fIter = A.facets_begin();
+ for (int i=0; i<sep_plane[2]; i++) ++fIter;
+ found = true;
+ for (Polyhedron::Vertex_iterator vIter = B.vertices_begin(); vIter != B.vertices_end(); ++vIter){
+ if (! fIter->plane().has_on_positive_side(vIter->point())) {found = false; break;};
+ }
+ if(found) return false;
+ }
+ break;
+ case 2: //separation plane was previously determined as sep_plane[2]-th plane of B polyhedron
+ {
+ if(unlikely((unsigned) sep_plane[2]>=B.size_of_facets())) break;
+ Polyhedron::Facet_iterator fIter = B.facets_begin();
+ for (int i=0; i<sep_plane[2]; i++) ++fIter;
+ found = true;
+ for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end(); ++vIter){
+ if (! fIter->plane().has_on_positive_side(vIter->point())) {found = false; break;};
+ }
+ if(found) return false;
+ }
+ break;
+ case 3: //separation plane was previously given by sep_plane[1]-th and sep_plane[2]-th edge of A & B polyhedrons
+ {
+ if(unlikely((unsigned) sep_plane[1]>=A.size_of_halfedges()/2)) break;
+ if(unlikely((unsigned) sep_plane[2]>=B.size_of_halfedges()/2)) break;
+ Polyhedron::Edge_iterator eIter1 = A.edges_begin();
+ Polyhedron::Edge_iterator eIter2 = B.edges_begin();
+ for (int i=0; i<sep_plane[1]; i++) ++eIter1;
+ for (int i=0; i<sep_plane[2]; i++) ++eIter2;
+ found = true;
+ Plane X(eIter1->vertex()->point(),CGAL::cross_product((eIter1->vertex()->point() - eIter1->opposite()->vertex()->point()),(eIter2->vertex()->point() - eIter2->opposite()->vertex()->point())));
+ if (!X.has_on_positive_side(B.vertices_begin()->point())) X = X.opposite();
+
+ double lim = pow(DISTANCE_LIMIT,2);
+ for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end(); ++vIter){
+ if (Oriented_squared_distance(X, vIter->point())>lim) {found = false; break;};
+ }
+ for (Polyhedron::Vertex_iterator vIter = B.vertices_begin(); vIter != B.vertices_end(); ++vIter){
+ if (! X.has_on_positive_side(vIter->point())) {found = false; break;};
+ }
+ if(found) return false;
+ }
+ break;
+ }
+
+ //regular test with no previous information about separating plane
+ //test all planes from A
+ int i = 0;
+ for (Polyhedron::Facet_iterator fIter = A.facets_begin(); fIter != A.facets_end(); fIter++, i++){
+ found = true;
+ for (Polyhedron::Vertex_iterator vIter = B.vertices_begin(); vIter != B.vertices_end(); ++vIter){
+ if (! fIter->plane().has_on_positive_side(vIter->point())) {found = false; break;};
+ }
+ if(found) {sep_plane[0] = 1; sep_plane[1] = 1; sep_plane[2] = i; return false;}
+ }
+ //test all planes from B
+ i = 0;
+ for (Polyhedron::Facet_iterator fIter = B.facets_begin(); fIter != B.facets_end(); fIter++, i++){
+ found = true;
+ for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end(); ++vIter){
+ if (! fIter->plane().has_on_positive_side(vIter->point())) {found = false; break;};
+ }
+ if(found) {sep_plane[0] = 2; sep_plane[1] = 2; sep_plane[2] = i; return false;}
+ }
+ //test all pairs of edges from A & B
+ Plane X;
+ CGALvector vA;
+ double lim = pow(DISTANCE_LIMIT,2);
+ i = 0;
+ for (Polyhedron::Edge_iterator eIter1 = A.edges_begin(); eIter1 != A.edges_end(); ++eIter1, i++){
+ vA = eIter1->vertex()->point() - eIter1->opposite()->vertex()->point();
+ int j = 0;
+ for (Polyhedron::Edge_iterator eIter2 = B.edges_begin(); eIter2 != B.edges_end(); ++eIter2, j++){
+ found = true;
+ X = Plane(eIter1->vertex()->point(),CGAL::cross_product(vA,(eIter2->vertex()->point() - eIter2->opposite()->vertex()->point()) ));
+ if (!X.has_on_positive_side(B.vertices_begin()->point())) X = X.opposite();
+ for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end(); ++vIter){
+ if (Oriented_squared_distance(X, vIter->point())>lim) {found = false; break;};
+ }
+ for (Polyhedron::Vertex_iterator vIter = B.vertices_begin(); vIter != B.vertices_end(); ++vIter){
+ if (! X.has_on_positive_side(vIter->point())) {found = false; break;};
+ }
+ if(found) {sep_plane[0] = 3; sep_plane[1] = i; sep_plane[2] = j; return false;}
+ }
+ }
+
+ sep_plane[0] = 0;
+ return true;
+}
+
+//**********************************************************************************
+//norm of difference between two planes
+double PlaneDifference(const Plane &a, const Plane &b){
+ double la = sqrt(pow(a.a(),2) + pow(a.b(),2) + pow(a.c(),2) + pow(a.d(),2));
+ double lb = sqrt(pow(b.a(),2) + pow(b.b(),2) + pow(b.c(),2) + pow(b.d(),2));
+ return pow(a.a()/la-b.a()/lb,2) + pow(a.b()/la-b.b()/lb,2) + pow(a.c()/la-b.c()/lb,2) + pow(a.d()/la-b.d()/lb,2);
+
+ //in case we do not care of the orientation
+ //return min(pow(a.a()/la-b.a()/lb,2) + pow(a.b()/la-b.b()/lb,2) + pow(a.c()/la-b.c()/lb,2) + pow(a.d()/la-b.d()/lb,2),pow(a.a()/la+b.a()/lb,2) + pow(a.b()/la+b.b()/lb,2) + pow(a.c()/la+b.c()/lb,2) + pow(a.d()/la+b.d()/lb,2));
+}
+
+//**********************************************************************************
+//connect triagular facets if possible
+Polyhedron Simplify(Polyhedron P, double limit){
+ bool elimination = true;
+ while(elimination){
+ elimination = false;
+ for (Polyhedron::Edge_iterator hei = P.edges_begin(); hei!=P.edges_end(); ++hei){
+ if (PlaneDifference(hei->facet()->plane(),hei->opposite()->facet()->plane()) < limit){
+ if (hei->vertex()->vertex_degree() < 3) hei = P.erase_center_vertex(hei);
+ else if(hei->opposite()->vertex()->vertex_degree() < 3) hei = P.erase_center_vertex(hei->opposite());
+ else hei = P.join_facet(hei);
+ elimination = true;
+ break;
+ }
+ }
+ }
+ if (P.size_of_facets() < 4) P.clear();
+ return P;
+}
+
+//**********************************************************************************
+//list of facets + edges
+void PrintPolyhedron(Polyhedron P){
+ Vector3r A,B,C;
+ cout << "*** faces ***" << endl;
+ 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());
+ cout << A << " "<< B << " "<< C << endl;
+ }
+ }
+ cout << "*** edges ***" << endl;
+ for (Polyhedron::Edge_iterator hei = P.edges_begin(); hei!=P.edges_end(); ++hei){
+ cout << hei->vertex()->point() << " " << hei->opposite()->vertex()->point() << endl;
+ }
+}
+
+//**********************************************************************************
+//list of facets
+void PrintPolyhedronFacets(Polyhedron P){
+ Vector3r A,B,C;
+ for (Polyhedron::Facet_iterator fIter = P.facets_begin(); fIter!=P.facets_end(); ++fIter){
+ cout << "***" << endl;
+ Polyhedron::Halfedge_around_facet_circulator hfc0;
+ hfc0 = fIter->facet_begin();
+ int n = fIter->facet_degree();
+ for (int i = 0; i<n; ++hfc0, i++){
+ cout << hfc0->vertex()->point() << endl;
+ }
+ }
+}
+
+//**********************************************************************************
+//calculate intersection of three planes (exceptional cases not checked) - not in use, can be deleted
+CGALpoint PlaneIntersection( Plane a, Plane b, Plane c){
+ Matrix3r A;
+ A<<
+ a.a(),a.b(),a.c(),
+ b.a(),b.b(),b.c(),
+ c.a(),c.b(),c.c();
+ Vector3r x = A.inverse()*Vector3r(-a.d(),-b.d(),-c.d());
+ return CGALpoint(x.x(),x.y(),x.z());
+}
+
+//**********************************************************************************
+//solve system of 3x3 by Cramers rule
+Vector3r SolveLinSys3x3(Matrix3r A, Vector3r y){
+ //only system 3x3 by Cramers rule
+ double det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1)-A(0,2)*A(1,1)*A(2,0)-A(0,1)*A(1,0)*A(2,2)-A(0,0)*A(1,2)*A(2,1);
+ if (det == 0) {cerr << "error in linear solver" << endl; return Vector3r(0,0,0);}
+ return Vector3r((y(0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*y(2)+A(0,2)*y(1)*A(2,1)-A(0,2)*A(1,1)*y(2)-A(0,1)*y(1)*A(2,2)-y(0)*A(1,2)*A(2,1))/det, (A(0,0)*y(1)*A(2,2)+y(0)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*y(2)-A(0,2)*y(1)*A(2,0)-y(0)*A(1,0)*A(2,2)-A(0,0)*A(1,2)*y(2))/det, (A(0,0)*A(1,1)*y(2)+A(0,1)*y(1)*A(2,0)+y(0)*A(1,0)*A(2,1)-y(0)*A(1,1)*A(2,0)-A(0,1)*A(1,0)*y(2)-A(0,0)*y(1)*A(2,1))/det);
+}
+
+//**********************************************************************************
+/*return convex hull of points
+critical point, because CGAL often returnes segmentation fault. The planes must be sufficiently "different". This is, however, checked elswhere by DISTANCE_LIMIT variable.
+*/
+Polyhedron ConvexHull(vector<CGALpoint> &planes){
+ Polyhedron Int;
+ //cout << "C"; cout.flush();
+ if (planes.size()>3) CGAL::convex_hull_3(planes.begin(), planes.end(), Int);
+ //cout << "-C"<< endl; cout.flush();
+ return Int;
+}
+
+//**********************************************************************************
+//determination of normal direction of intersection
+
+Vector3r FindNormal(Polyhedron Int, Polyhedron PA, Polyhedron PB){
+
+ //determine which plane is from which polyhedra
+ Polyhedron::Plane_iterator pi, pj;
+ std::transform( Int.facets_begin(), Int.facets_end(), Int.planes_begin(),Plane_equation());
+ std::transform( PA.facets_begin(), PA.facets_end(), PA.planes_begin(),Plane_equation());
+ std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());
+ vector<bool> from_A(Int.size_of_facets());
+ vector<double> minsA(Int.size_of_facets());
+ vector<double> minsB(Int.size_of_facets());
+ int i=0;
+ double minA, minB, k;
+ for(pi = Int.planes_begin(); pi!=Int.planes_end(); ++pi,i++){
+ minA = 1.;
+ minB = 1.;
+ for(pj=PA.planes_begin(); pj!=PA.planes_end(); ++pj){
+ k = PlaneDifference(*pi, *pj);
+ if (k<minA) {
+ minA = k;
+ minsA[i] = minA;
+ if (minA<FIND_NORMAL_LIMIT) {from_A[i] = true; break;} //already satisfactory
+ }
+
+ }
+ if (minA<FIND_NORMAL_LIMIT) continue;
+ for(pj=PB.planes_begin(); pj!=PB.planes_end(); ++pj){
+ k = PlaneDifference(*pi, *pj);
+
+ if (k<minB){
+ minB = k;
+ minsB[i] = minB;
+ if (minB<FIND_NORMAL_LIMIT || minB<minA) break; //already satisfactory
+ }
+ }
+ from_A[i] = ((minA < minB) ? true : false);
+ }
+ //check that not all belongs to A and not all belongs to B
+ if (*std::min_element(from_A.begin(),from_A.end())==1){
+ int loc = std::min_element(minsB.begin(),minsB.end()) - minsB.begin();
+ from_A[loc] = false;
+ }else if (*std::max_element(from_A.begin(),from_A.end())==0){
+ int loc = std::min_element(minsA.begin(),minsA.end()) - minsA.begin();
+ from_A[loc] = true;
+ }
+
+ //find intersecting segments
+ vector<Segment> segments;
+ int a,b;
+
+ for (Polyhedron::Edge_iterator hei = Int.edges_begin(); hei!=Int.edges_end(); ++hei){
+ a = std::distance(Int.facets_begin(), hei->facet());
+ b = std::distance(Int.facets_begin(), hei->opposite()->facet());
+ if ((from_A[a] && !from_A[b]) || (from_A[b] && !from_A[a])) segments.push_back(Segment(hei->vertex()->point(),hei->opposite()->vertex()->point()));
+ }
+
+ //find normal direction
+ Plane fit;
+ linear_least_squares_fitting_3(segments.begin(),segments.end(),fit,CGAL::Dimension_tag<1>());
+ CGALvector CGALnormal=fit.orthogonal_vector();
+ CGALnormal = CGALnormal/sqrt(CGALnormal.squared_length());
+ // reverse direction if projection of the (contact_point-centroid_of_B) vector onto the normal is negative (i.e. the normal points more towards B)
+
+ return FromCGALVector(CGALnormal);
+}
+
+//**********************************************************************************
+//not used
+double CalculateProjectionArea(Polyhedron Int, CGALvector CGALnormal){
+ //calculate area of projection of Intersection into the normal plane
+ double area = 0.;
+ double abs_size;
+ Polyhedron::Halfedge_around_facet_circulator hfc0;
+ CGALvector normal2;
+ double norm2;
+ std::transform( Int.facets_begin(), Int.facets_end(), Int.planes_begin(),Plane_equation());
+ for(Polyhedron::Facet_iterator fi = Int.facets_begin(); fi!= Int.facets_end(); ++fi){
+ assert(fi->is_triangle());
+ hfc0 = fi->facet_begin();
+ normal2 = fi->plane().orthogonal_vector();
+ norm2 = normal2.squared_length();
+ if(norm2 < 1E-20) continue;
+ abs_size = 0.5*sqrt((cross_product(CGALvector(hfc0->vertex()->point(),hfc0->next()->vertex()->point()),CGALvector(hfc0->vertex()->point(),hfc0->next()->next()->vertex()->point()))).squared_length());
+ // factor 0.5 because this procedure returnes doubles projected area
+ if (abs_size>0) area += 0.5*abs_size*abs(CGALnormal*normal2/sqrt(norm2));
+ }
+ return area;
+}
+
+//**********************************************************************************
+//prepare data for CGAL convex hull
+vector<Plane> MergePlanes(vector<Plane> planes1, vector<Plane> planes2, double limit){
+ vector<Plane> P = planes1;
+ bool add;
+ for(vector<Plane>::iterator i = planes2.begin(); i!=planes2.end(); ++i){
+ add = true;
+ for(vector<Plane>::iterator j = planes1.begin(); j!=planes1.end(); ++j){
+ if (PlaneDifference(*i,*j) < limit){
+ add = false;
+ break;
+ }
+ }
+ if (add) P.push_back(*i);
+ }
+ return P;
+}
+
+//**********************************************************************************
+//returnes intersecting polyhedron of polyhedron & plane (possibly empty)
+Polyhedron Polyhedron_Plane_intersection(Polyhedron A, Plane B, CGALpoint centroid, CGALpoint X){
+
+ Polyhedron Intersection;
+ CGALpoint inside;
+ vector<Plane> planes1, planes2;
+ vector<CGALpoint> dual_planes;
+ // test if do intersect, find some intersecting point
+ bool intersection_found = false;
+ double lim = pow(DISTANCE_LIMIT,2);
+ // test centroid of previous intersection
+ if(Is_inside_Polyhedron(A,X,DISTANCE_LIMIT) && Oriented_squared_distance(B,X)<=-lim) {
+ intersection_found = true;
+ inside = X;
+ // find new point by checking polyhedron vertices that lies on negative side of the plane
+ } 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);
+ }else{
+ cerr << "Error in line-plane intersection" << endl;
+ }
+ }
+ }
+ }
+ //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);
+
+ std::transform( A.points_begin(), A.points_end(), A.points_begin(), transl);
+ B = transl.transform(B);
+
+ //dualize plane
+ planes1.push_back(B);
+
+ //dualize polyhedron
+ std::transform( A.facets_begin(), A.facets_end(), A.planes_begin(),Plane_equation());
+ for(Polyhedron::Plane_iterator pi =A.planes_begin(); pi!=A.planes_end(); ++pi) planes2.push_back(*pi);;
+
+ //merge planes
+ planes1 = MergePlanes(planes1,planes2, MERGE_PLANES_LIMIT);//MERGE_PLANES_LIMIT);
+ for(vector<Plane>::iterator pi =planes1.begin(); pi!= planes1.end(); ++pi) dual_planes.push_back(CGALpoint(-pi->a()/pi->d(),-pi->b()/pi->d(),-pi->c()/pi->d()));
+
+ //compute convex hull of it
+ Intersection = ConvexHull(dual_planes);
+ if (Intersection.empty()) {cout << "0" << endl; cout.flush(); return Intersection;}
+
+ //simplify
+ std::transform( Intersection.facets_begin(), Intersection.facets_end(), Intersection.planes_begin(),Plane_equation());
+ Intersection = Simplify(Intersection, SIMPLIFY_LIMIT);
+ std::transform( Intersection.facets_begin(), Intersection.facets_end(), Intersection.planes_begin(),Plane_equation());
+
+ //dualize again
+ dual_planes.clear();
+ for(Polyhedron::Plane_iterator pi =Intersection.planes_begin(); pi!=Intersection.planes_end(); ++pi) dual_planes.push_back(CGALpoint(-pi->a()/pi->d(),-pi->b()/pi->d(),-pi->c()/pi->d()));
+
+ //compute convex hull of it
+ Intersection = ConvexHull(dual_planes);
+
+ //return to original position
+ std::transform( Intersection.points_begin(), Intersection.points_end(), Intersection.points_begin(), transl_back);
+
+ if (Intersection.size_of_facets() < 4) Intersection.clear();
+ return Intersection;
+}
+
+//**********************************************************************************
+//returnes intersecting polyhedron of polyhedron & plane defined by diriction 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);
+ return Polyhedron_Plane_intersection(A, B, ToCGALPoint(plane_point), X);
+}
+
+//**********************************************************************************
+//returnes intersecting polyhedron of two polyhedrons (possibly empty)
+Polyhedron Polyhedron_Polyhedron_intersection(Polyhedron A, Polyhedron B, CGALpoint X, CGALpoint centroidA, CGALpoint centroidB, std::vector<int> &sep_plane){
+
+ Polyhedron Intersection;
+
+ vector<Plane> planes1, planes2;
+ vector<CGALpoint> dual_planes;
+ Polyhedron::Plane_iterator pi;
+ CGALpoint inside(0,0,0);
+
+ bool intersection_found = false;
+ std::transform( A.facets_begin(), A.facets_end(), A.planes_begin(),Plane_equation());
+ std::transform( B.facets_begin(), B.facets_end(), B.planes_begin(),Plane_equation());
+ Matrix3r Amatrix;
+ Vector3r y;
+ // test that X is really inside
+ if(Is_inside_Polyhedron(A,X,DISTANCE_LIMIT) && Is_inside_Polyhedron(B,X,DISTANCE_LIMIT)) {
+ intersection_found = true;
+ inside = X;
+ } else {
+ if (!do_intersect(A, B, sep_plane)) return Intersection;
+ //some intersection point
+ double dist_S, dist_T;
+ double lim2 = pow(DISTANCE_LIMIT,2);
+ CGALvector d1;
+ double factor = sqrt(DISTANCE_LIMIT * 1.5);
+ //test vertices A - not needed, edges are enough
+ for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end() && !intersection_found; vIter++){
+ d1 = centroidA-vIter->point();
+ inside = vIter->point() + d1/sqrt(d1.squared_length())*DISTANCE_LIMIT*20.;
+ intersection_found = (Is_inside_Polyhedron(A,inside,DISTANCE_LIMIT) && Is_inside_Polyhedron(B,inside,DISTANCE_LIMIT));
+ }
+ //test vertices B - necessary
+ for (Polyhedron::Vertex_iterator vIter = B.vertices_begin(); vIter != B.vertices_end() && !intersection_found; vIter++){
+ d1 = centroidB-vIter->point();
+ inside = vIter->point() + d1/sqrt(d1.squared_length())*DISTANCE_LIMIT*20.;
+ intersection_found = (Is_inside_Polyhedron(A,inside,DISTANCE_LIMIT) && Is_inside_Polyhedron(B,inside,DISTANCE_LIMIT));
+ }
+
+ //test edges
+ for (Polyhedron::Edge_iterator eIter = A.edges_begin(); eIter != A.edges_end() && !intersection_found; eIter++){
+ for (Polyhedron::Facet_iterator fIter = B.facets_begin(); fIter != B.facets_end() && !intersection_found; fIter++){
+ dist_S = Oriented_squared_distance(fIter->plane(), eIter->vertex()->point());
+ dist_T = Oriented_squared_distance(fIter->plane(), eIter->opposite()->vertex()->point());
+ if (dist_S*dist_T >= 0 || abs(dist_S)<lim2 || abs(dist_T)<lim2) continue;
+ inside = eIter->vertex()->point() + (eIter->opposite()->vertex()->point()-eIter->vertex()->point())*sqrt(abs(dist_S))/(sqrt(abs(dist_S))+sqrt(abs(dist_T)));
+ // the fact that edge intersects the facet (not only its plane) is not explicitely checked, it sufices to check that the resulting point is inside both polyhedras
+ Plane p1 = fIter->plane();
+ Plane p2 = eIter->facet()->plane();
+ Plane p3 = eIter->opposite()->facet()->plane();
+ Amatrix << p1.a(),p1.b(),p1.c(),
+ p2.a(),p2.b(),p2.c(),
+ p3.a(),p3.b(),p3.c();
+ y = Vector3r(-p1.d()-factor*sqrt(pow(p1.a(),2)+pow(p1.b(),2)+pow(p1.c(),2)),-p2.d()-factor*sqrt(pow(p2.a(),2)+pow(p2.b(),2)+pow(p2.c(),2)),-p3.d()-factor*sqrt(pow(p3.a(),2)+pow(p3.b(),2)+pow(p3.c(),2)));
+ inside = ToCGALPoint(SolveLinSys3x3(Amatrix,y));
+ intersection_found = (Is_inside_Polyhedron(A,inside,DISTANCE_LIMIT) && Is_inside_Polyhedron (B,inside,DISTANCE_LIMIT));
+ }
+ }
+ }
+
+ //Polyhedrons do not intersect
+ 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);
+
+ std::transform( A.points_begin(), A.points_end(), A.points_begin(), transl);
+ std::transform( B.points_begin(), B.points_end(), B.points_begin(), transl);
+
+ //dualize polyhedrons
+ std::transform( A.facets_begin(), A.facets_end(), A.planes_begin(),Plane_equation());
+ std::transform( B.facets_begin(), B.facets_end(), B.planes_begin(),Plane_equation());
+ for(Polyhedron::Plane_iterator pi =A.planes_begin(); pi!=A.planes_end(); ++pi) planes1.push_back(*pi);
+ for(Polyhedron::Plane_iterator pi =B.planes_begin(); pi!=B.planes_end(); ++pi) planes2.push_back(*pi);
+
+
+ //merge planes
+ planes1 = MergePlanes(planes1,planes2, MERGE_PLANES_LIMIT);//MERGE_PLANES_LIMIT);
+ for(vector<Plane>::iterator pi =planes1.begin(); pi!= planes1.end(); ++pi) dual_planes.push_back(CGALpoint(-pi->a()/pi->d(),-pi->b()/pi->d(),-pi->c()/pi->d()));
+
+ //compute convex hull of it
+ Intersection = ConvexHull(dual_planes);
+ if (Intersection.empty()) {cout << "0" << endl; cout.flush(); return Intersection;};
+
+ //simplify
+ std::transform( Intersection.facets_begin(), Intersection.facets_end(), Intersection.planes_begin(),Plane_equation());
+ Intersection = Simplify(Intersection, SIMPLIFY_LIMIT);
+ std::transform( Intersection.facets_begin(), Intersection.facets_end(), Intersection.planes_begin(),Plane_equation());
+
+ //dualize again
+ dual_planes.clear();
+ for(Polyhedron::Plane_iterator pi =Intersection.planes_begin(); pi!=Intersection.planes_end(); ++pi) dual_planes.push_back(CGALpoint(-pi->a()/pi->d(),-pi->b()/pi->d(),-pi->c()/pi->d()));
+
+ //compute convex hull of it
+ Intersection = ConvexHull(dual_planes);
+ //return to original position
+ std::transform( Intersection.points_begin(), Intersection.points_end(), Intersection.points_begin(), transl_back);
+
+ if (Intersection.size_of_facets() < 4) Intersection.clear();
+ return Intersection;
+}
+
+//**********************************************************************************
+/*
+//returnes intersection of Sphere and polyhedron (possibly empty) - not finished
+bool Sphere_Polyhedron_intersection(Polyhedron A, double r, CGALpoint C, CGALpoint centroid, double volume, CGALvector normal, double area){
+ volume = 0;
+ area = 0;
+
+ double r2 = pow(r,2);
+ //check points
+ Polyhedron::Vertex_iterator closest_vert;
+ double dist2;
+ double min_dist2 = (C-centroid).squared_length();
+ //test vertices
+ for(Polyhedron::Vertex_iterator vi=A.vertices_begin(); vi!=A.vertices_end(); ++vi){
+ dist2 = (C-vi->point()).squared_length();
+ if(min_dist2 > dist2) {min_dist2 = dist2; closest_vert=vi;}
+ }
+ if(min_dist2 < r2){ //test if we are already inside the sphere
+
+ return true;
+ }
+
+ //test edges
+ CGALvector v1(C-closest_vert->point());
+ v1 = v1/sqrt(v1.squared_length());
+ Polyhedron::Halfedge_around_vertex_circulator closest_edge;
+ int i =0;
+ double cos_angle;
+ double max_cos_angle=-1;
+ CGALvector v2;
+ for(Polyhedron::Halfedge_around_vertex_circulator hfi=closest_vert->vertex_begin(); i<(int)closest_vert->degree(); ++hfi, i++){
+ v2 = hfi->opposite()->vertex()->point()-closest_vert->point();
+ v2 = v2/sqrt(v2.squared_length());
+ cos_angle = v1*v2;
+ if(cos_angle > max_cos_angle) closest_edge=hfi; max_cos_angle=cos_angle; //find the edge with minimum angle
+ }
+ if(max_cos_angle<=0) return false; //all edges goes with angle above or equal 90
+ double k = (C - closest_vert->point())*v2; //find closest point on the most descending edge (it must be on edge, otherwise opposite point would be already closer)
+ CGALpoint closest_point = closest_vert->point() + k*v2;
+ if((C-closest_point).squared_length()<r2){ //test if we are already inside the sphere
+
+ return true;
+ }
+
+
+ //check planes
+ v1 = (C-closest_point);
+ v1 = v1/sqrt(v1.squared_length());
+ v2 = closest_edge->vertex()->point()-closest_edge->opposite()->vertex()->point();
+ CGALvector v3 = closest_edge->facet()->plane().orthogonal_vector();
+ CGALvector v4 = closest_edge->opposite()->facet()->plane().orthogonal_vector();
+ v3 = CGAL::cross_product(v3,v2);
+ v4 = CGAL::cross_product(v4,v2*(-1));
+ Plane closest_plane;
+ if (v1*v3>=0 && v1*v4>=0) return false; // both planes running away
+ if (v1*v3<0) closest_plane = closest_edge->facet()->plane();
+ else closest_plane = closest_edge->opposite()->facet()->plane();
+ closest_point = closest_plane.projection(C);
+ if((C-closest_point).squared_length()<r2){ //test if we are already inside the sphere
+ double h = sqrt((C-closest_point).squared_length());
+ volume = 3.1415*pow(h,2)/3*(3*r-h);
+ return true;
+ }
+
+ return false;
+}
+*/
+
+//**********************************************************************************
+Vector3r FromCGALPoint(CGALpoint A){
+ return Vector3r(A.x(), A.y(), A.z());
+}
+
+//**********************************************************************************
+Vector3r FromCGALVector(CGALvector A){
+ return Vector3r(A.x(), A.y(), A.z());
+}
+
+//**********************************************************************************
+CGALpoint ToCGALPoint(Vector3r A){
+ return CGALpoint(A[0], A[1], A[2]);
+}
+
+//**********************************************************************************
+CGALvector ToCGALVector(Vector3r A){
+ return CGALvector(A[0], A[1], A[2]);
+}
+
+//**********************************************************************************
+//new polyhedra
+shared_ptr<Body> NewPolyhedra(vector<Vector3r> v, shared_ptr<Material> mat){
+ shared_ptr<Body> body(new Body);
+ body->material=mat;
+ body->shape=shared_ptr<Polyhedra>(new Polyhedra());
+ Polyhedra* A = static_cast<Polyhedra*>(body->shape.get());
+ A->v = v;
+ A->Initialize();
+ body->state->pos= A->GetCentroid();
+ body->state->mass=body->material->density*A->GetVolume();
+ body->state->inertia =A->GetInertia()*body->material->density;
+ body->state->ori = A->GetOri();
+ body->bound=shared_ptr<Aabb>(new Aabb);
+ body->setAspherical(true);
+ return body;
+}
+
+
+//**********************************************************************************
+//split polyhedra
+void SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction){
+
+ const Se3r& se3=body->state->se3;
+ Polyhedra* A = static_cast<Polyhedra*>(body->shape.get());
+ State* X = static_cast<State*>(body->state.get());
+
+ Vector3r OrigPos = X->pos;
+ Vector3r OrigVel = X->vel;
+ Vector3r OrigAngVel = X->angVel;
+
+ //move and rotate CGAL structure Polyhedron
+ Matrix3r rot_mat = (se3.orientation).toRotationMatrix();
+ Vector3r trans_vec = se3.position;
+ 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);
+
+ //calculate splitted polyhedrons
+ Polyhedron S1 = Polyhedron_Plane_intersection(PA, direction, trans_vec);
+ Polyhedron S2 = Polyhedron_Plane_intersection(PA, (-1)*direction, trans_vec);
+
+
+ //replace original polyhedron
+ A->Clear();
+ for(Polyhedron::Vertex_iterator vi = S1.vertices_begin(); vi != S1.vertices_end(); vi++){ A->v.push_back(FromCGALPoint(vi->point()));};
+ A->Initialize();
+ X->pos = A->GetCentroid();
+ X->ori = A->GetOri();
+ X->mass=body->material->density*A->GetVolume();
+ X->refPos[0] = X->refPos[0]+1.;
+ X->inertia =A->GetInertia()*body->material->density;
+
+ //second polyhedron
+ vector<Vector3r> v2;
+ for(Polyhedron::Vertex_iterator vi = S2.vertices_begin(); vi != S2.vertices_end(); vi++){ v2.push_back(FromCGALPoint(vi->point()));};
+ 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;
+
+}
+
+#endif // YADE_CGAL
=== modified file 'py/CMakeLists.txt'
--- py/CMakeLists.txt 2013-08-29 10:30:31 +0000
+++ py/CMakeLists.txt 2013-10-15 16:14:46 +0000
@@ -36,6 +36,12 @@
INSTALL(TARGETS _utils DESTINATION "${YADE_PY_PATH}/yade/")
+ADD_LIBRARY(_polyhedra_utils SHARED "${CMAKE_CURRENT_SOURCE_DIR}/_polyhedra_utils.cpp")
+SET_TARGET_PROPERTIES(_polyhedra_utils PROPERTIES PREFIX "")
+TARGET_LINK_LIBRARIES(_polyhedra_utils)
+INSTALL(TARGETS _polyhedra_utils DESTINATION "${YADE_PY_PATH}/yade/")
+
+
ADD_LIBRARY(_packPredicates SHARED "${CMAKE_CURRENT_SOURCE_DIR}/pack/_packPredicates.cpp")
SET_TARGET_PROPERTIES(_packPredicates PROPERTIES PREFIX "")
TARGET_LINK_LIBRARIES(_packPredicates yade)
=== added file 'py/_polyhedra_utils.cpp'
--- py/_polyhedra_utils.cpp 1970-01-01 00:00:00 +0000
+++ py/_polyhedra_utils.cpp 2013-10-15 16:14:46 +0000
@@ -0,0 +1,374 @@
+// © 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+#ifdef YADE_CGAL
+
+#include"yade/pkg/dem/Polyhedra.hpp"
+
+#include<boost/python.hpp>
+#include <boost/python/module.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/lib/pyutil/doc_opts.hpp>
+#include<cmath>
+
+#include<numpy/ndarrayobject.h>
+
+using namespace std;
+namespace py = boost::python;
+
+
+//**********************************************************************************
+//print polyhedron in basic position
+void PrintPolyhedra(const shared_ptr<Shape>& shape){
+ Polyhedra* A = static_cast<Polyhedra*>(shape.get());
+ Polyhedron PA = A->GetPolyhedron();
+ A->Initialize();
+ PrintPolyhedron(PA);
+}
+
+//**********************************************************************************
+//print polyhedron in actual position
+void PrintPolyhedraActualPos(const shared_ptr<Shape>& cm1,const State& state1){
+ const Se3r& se3=state1.se3;
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+ A->Initialize();
+
+ //move and rotate CGAL structure Polyhedron
+ Matrix3r rot_mat = (se3.orientation).toRotationMatrix();
+ Vector3r trans_vec = se3.position;
+ 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);
+
+ PrintPolyhedron(PA);
+}
+
+
+//**********************************************************************************
+//test of polyhedron intersection callable from python shell
+bool do_Polyhedras_Intersect(const shared_ptr<Shape>& cm1,const shared_ptr<Shape>& cm2,const State& state1,const State& state2){
+
+ const Se3r& se31=state1.se3;
+ const Se3r& se32=state2.se3;
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+ Polyhedra* B = static_cast<Polyhedra*>(cm2.get());
+
+ //move and rotate 1st the CGAL structure Polyhedron
+ Matrix3r rot_mat = (se31.orientation).toRotationMatrix();
+ Vector3r trans_vec = se31.position;
+ 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);
+
+ //move and rotate 2st the CGAL structure Polyhedron
+ rot_mat = (se32.orientation).toRotationMatrix();
+ trans_vec = se32.position;
+ 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);
+
+ //calculate plane equations
+ std::transform( PA.facets_begin(), PA.facets_end(), PA.planes_begin(),Plane_equation());
+ std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());
+
+
+ //call test
+ return do_intersect(PA,PB);
+}
+
+//**********************************************************************************
+//determination of critical time step for polyhedrons & spheres (just rough estimation)
+Real PWaveTimeStep(){
+ const shared_ptr<Scene> _rb=shared_ptr<Scene>();
+ shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
+ Real dt=std::numeric_limits<Real>::infinity();
+ FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ if(!b || !b->material || !b->shape) continue;
+ shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->shape);
+ shared_ptr<Polyhedra> p=dynamic_pointer_cast<Polyhedra>(b->shape);
+ if(!s && !p) continue;
+ if(!p){
+ //spheres
+ shared_ptr<ElastMat> ebp=dynamic_pointer_cast<ElastMat>(b->material);
+ if(!ebp) continue;
+ Real density=b->state->mass/((4./3.)*Mathr::PI*pow(s->radius,3));
+ dt=min(dt,s->radius/sqrt(ebp->young/density));
+ }else{
+ //polyhedrons
+ shared_ptr<PolyhedraMat> ebp=dynamic_pointer_cast<PolyhedraMat>(b->material);
+ if(!ebp) continue;
+ 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));
+ }
+ }
+ if (dt==std::numeric_limits<Real>::infinity()) {
+ dt = 1.0;
+ LOG_WARN("PWaveTimeStep has not found any suitable spherical or polyhedral body to calculate dt. dt is set to 1.0");
+ }
+ return dt;
+}
+
+//**********************************************************************************
+//returns approximate sieve size of polyhedron
+Real SieveSize(const shared_ptr<Shape>& cm1){
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+
+ double phi = M_PI/4.;
+ double x,y;
+ double minx = 0, maxx = 0, miny = 0, maxy = 0;
+
+ for (vector<Vector3r>::iterator i=A->v.begin(); i!=A->v.end(); ++i){
+ x = cos(phi)*(*i)[1] + sin(phi)*(*i)[2];
+ y = -sin(phi)*(*i)[1] + cos(phi)*(*i)[2];
+ minx = min(minx,x);
+ maxx = max(maxx,x);
+ miny = min(miny,y);
+ maxy = max(maxy,y);
+ }
+
+ return max(maxx-minx, maxy-miny);
+}
+
+
+//**********************************************************************************
+//returns approximate size of polyhedron
+Vector3r SizeOfPolyhedra(const shared_ptr<Shape>& cm1){
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+
+ double minx = 0, maxx = 0, miny = 0, maxy = 0, minz = 0, maxz = 0;
+
+ for (vector<Vector3r>::iterator i=A->v.begin(); i!=A->v.end(); ++i){
+ minx = min(minx,(*i)[0]);
+ maxx = max(maxx,(*i)[0]);
+ miny = min(miny,(*i)[1]);
+ maxy = max(maxy,(*i)[1]);
+ minz = min(minz,(*i)[2]);
+ maxz = max(maxz,(*i)[2]);
+ }
+
+ return Vector3r(maxx-minx, maxy-miny, maxz-minz);
+}
+
+//**********************************************************************************
+//save sieve curve points into a file
+void SieveCurve(){
+ const shared_ptr<Scene> _rb=shared_ptr<Scene>();
+ shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
+ std::vector< std::pair<double,double> > sieve_volume;
+ double total_volume = 0;
+ FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ if(!b || !b->shape) continue;
+ shared_ptr<Polyhedra> p=dynamic_pointer_cast<Polyhedra>(b->shape);
+ if(p){
+ sieve_volume.push_back(std::pair<double,double>(SieveSize(p),p->GetVolume()));
+ total_volume += p->GetVolume();
+ }
+ }
+
+ std::sort(sieve_volume.begin(), sieve_volume.end()) ;
+ double cumul_vol = 0;
+
+ ofstream myfile;
+ myfile.open ("sieve_curve.dat");
+ for(std::vector< std::pair<double,double> > :: iterator i = sieve_volume.begin(); i != sieve_volume.end(); ++i) {
+ cumul_vol += i->second/total_volume;
+ myfile << i->first << "\t" << cumul_vol << endl;
+ }
+ myfile.close();
+}
+
+//**********************************************************************************
+//save size of polyhedrons into a file
+void SizeRatio(){
+ const shared_ptr<Scene> _rb=shared_ptr<Scene>();
+ shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
+ ofstream myfile;
+ myfile.open ("sizes.dat");
+ FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ if(!b || !b->shape) continue;
+ shared_ptr<Polyhedra> p=dynamic_pointer_cast<Polyhedra>(b->shape);
+ if(p){
+ myfile << SizeOfPolyhedra(p) << endl;
+ }
+ }
+ myfile.close();
+}
+
+//**********************************************************************************
+//returns max coordinates
+Vector3r MaxCoord(const shared_ptr<Shape>& cm1,const State& state1){
+ const Se3r& se3=state1.se3;
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+
+ //move and rotate CGAL structure Polyhedron
+ Matrix3r rot_mat = (se3.orientation).toRotationMatrix();
+ Vector3r trans_vec = se3.position;
+ 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);
+
+ Vector3r maxccord = trans_vec;
+ for(Polyhedron::Vertex_iterator vi = PA.vertices_begin(); vi != PA.vertices_end(); ++vi){
+ if (vi->point()[0]>maxccord[0]) maxccord[0]=vi->point()[0];
+ if (vi->point()[1]>maxccord[1]) maxccord[1]=vi->point()[1];
+ if (vi->point()[2]>maxccord[2]) maxccord[2]=vi->point()[2];
+ }
+
+ return maxccord;
+}
+
+//**********************************************************************************
+//returns min coordinates
+Vector3r MinCoord(const shared_ptr<Shape>& cm1,const State& state1){
+ const Se3r& se3=state1.se3;
+ Polyhedra* A = static_cast<Polyhedra*>(cm1.get());
+
+ //move and rotate CGAL structure Polyhedron
+ Matrix3r rot_mat = (se3.orientation).toRotationMatrix();
+ Vector3r trans_vec = se3.position;
+ 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);
+
+ Vector3r minccord = trans_vec;
+ for(Polyhedron::Vertex_iterator vi = PA.vertices_begin(); vi != PA.vertices_end(); ++vi){
+ if (vi->point()[0]<minccord[0]) minccord[0]=vi->point()[0];
+ if (vi->point()[1]<minccord[1]) minccord[1]=vi->point()[1];
+ if (vi->point()[2]<minccord[2]) minccord[2]=vi->point()[2];
+ }
+
+ return minccord;
+}
+
+
+//**********************************************************************************
+//generate "packing" of non-overlapping polyhedrons
+vector<Vector3r> fillBox_cpp(Vector3r minCoord, Vector3r maxCoord, Vector3r sizemin, Vector3r sizemax, Vector3r ratio, int seed, shared_ptr<Material> mat){
+ 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;
+
+ 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]);
+ }
+
+ //it - number of trials to make packing possibly more/less dense
+ Vector3r random_size;
+ while (it<1000){
+ it = it+1;
+ if (it == 1){
+ trialP.Clear();
+ trialP.seed = rand();
+ if(fixed_ratio) trialP.size = (rand()*(sizemax[0]-sizemin[0])/RAND_MAX + sizemin[0])*ratio;
+ else trialP.size = Vector3r(rand()*(sizemax[0]-sizemin[0]),rand()*(sizemax[1]-sizemin[1]),rand()*(sizemax[2]-sizemin[2]))/RAND_MAX + sizemin;
+ 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);
+}
+
+//**********************************************************************************
+//distace of point from a plane (squared) with sign
+double Oriented_squared_distance2(Plane P, CGALpoint x){
+ double h = P.a()*x.x()+P.b()*x.y()+P.c()*x.z()+P.d();
+ return ((h>0.)-(h<0.))*pow(h,2)/(CGALvector(P.a(),P.b(),P.c())).squared_length();
+}
+
+//**********************************************************************************
+bool convexHull(vector<Vector3r> points){
+ vector<CGALpoint> pointsCGAL;
+ for(int i=0;i<(int) points.size() ;i++) {
+ pointsCGAL.push_back(ToCGALPoint(points[i]));
+ }
+ Polyhedron P;
+ CGAL::convex_hull_3(pointsCGAL.begin(), pointsCGAL.end(), P);
+ return true;
+}
+
+BOOST_PYTHON_MODULE(_polyhedra_utils){
+ // http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
+ import_array();
+
+ YADE_SET_DOCSTRING_OPTS;
+
+ py::def("PrintPolyhedra",PrintPolyhedra,"Print list of vertices sorted according to polyhedrons facets.");
+ py::def("PrintPolyhedraActualPos",PrintPolyhedraActualPos,"Print list of vertices sorted according to polyhedrons facets.");
+ 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("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");
+}
+
+#endif // YADE_CGAL
=== added file 'py/polyhedra_utils.py'
--- py/polyhedra_utils.py 1970-01-01 00:00:00 +0000
+++ py/polyhedra_utils.py 2013-10-15 16:14:46 +0000
@@ -0,0 +1,102 @@
+# 2013 Jan Elias, http://www.fce.vutbr.cz/STM/elias.j/, elias.j@xxxxxxxxxxxx
+# https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
+
+"""
+Auxiliary functions for polyhedra
+"""
+
+
+import math,random,doctest,geom,numpy
+from yade import Vector3
+from yade.wrapper import *
+#from miniEigen import *
+try: # use psyco if available
+ import psyco
+ psyco.full()
+except ImportError: pass
+
+
+# c++ implementations for performance reasons
+from yade._polyhedra_utils import *
+
+#**********************************************************************************
+def randomColor(seed=None):
+ random.seed(seed);
+ #Return random Vector3 with each component in interval 0...1 (uniform distribution)
+ return Vector3(random.random(),random.random(),random.random())
+
+#**********************************************************************************
+#create polyhedra, one can specify vertices directly, or leave it empty for random shape
+def polyhedra(material,size=Vector3(1,1,1),seed=None,v=[],mask=1,fixed=False, color=[-1,-1,-1]):
+ """create polyhedra, one can specify vertices directly, or leave it empty for random shape.
+
+ :param Material material: material of new body
+ :param Vector3 size: size of new body (see Polyhedra docs)
+ :param float seed: seed for random operations
+ :param [Vector3] v: list of body vertices (see Polyhedra docs)
+ """
+ b=Body()
+ random.seed(seed);
+ b.aspherical = True
+ if len(v)>0:
+ b.shape = Polyhedra(v=v)
+ else:
+ b.shape = Polyhedra(size = size, seed=random.randint(0,1E6))
+ if color[0] == -1:
+ b.shape.color = randomColor(seed=random.randint(0,1E6))
+ else:
+ b.shape.color = color
+ b.mat = material
+ b.state.mass = b.mat.density*b.shape.GetVolume()
+ b.state.inertia = b.shape.GetInertia()*b.mat.density
+ b.state.ori = b.shape.GetOri()
+ b.state.pos = b.shape.GetCentroid()
+ b.mask=mask
+ if fixed:
+ b.state.blockedDOFs = 'xyzXYZ'
+ return b
+
+#**********************************************************************************
+#creates polyhedra having N vertices and resembling sphere
+def polyhedralBall(radius, N, material, center,mask=1):
+ """creates polyhedra having N vertices and resembling sphere
+
+ :param float radius: ball radius
+ :param int N: number of vertices
+ :param Material material: material of new body
+ :param Vector3 center: center of the new body
+ """
+ pts = []
+
+ inc = math.pi * (3. - math.sqrt(5.))
+ off = 2. / float(N)
+ for k in range(0, N):
+ y = k * off - 1. + (off / 2.)
+ r = math.sqrt(1. - y*y)
+ phi = k * inc
+ pts.append([math.cos(phi)*r*radius, y*radius, math.sin(phi)*r*radius])
+
+ ball = polyhedra(material,v=pts)
+ ball.state.pos = center
+ 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)
+ :param Vector3 mincoord: first corner
+ :param Vector3 maxcoord: second corner
+ :param Vector3 sizemin: minimal size of bodies
+ :param Vector3 sizemax: maximal size of bodies
+ :param Vector3 ratio: scaling ratio
+ :param float seed: random seed
+ """
+ random.seed(seed);
+ v = fillBox_cpp(mincoord, maxcoord, sizemin,sizemax, ratio, random.randint(0,1E6), material)
+ #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
+
+
Follow ups
