← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-pkg/yade/git-trunk] Rev 3499: PolyhedraMat and PolyhedraPhys are now inherited from FrictMat and FrictPhys, modified related ex...

 

------------------------------------------------------------
revno: 3499
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Wed 2014-10-22 12:13:18 +0200
message:
  PolyhedraMat and PolyhedraPhys are now inherited from FrictMat and FrictPhys, modified related examples
  Added Ig2_Sphere_Polyhedra_ScGeom
modified:
  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
  py/_polyhedra_utils.cpp
  py/utils.py


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'examples/polyhedra/ball.py'
--- examples/polyhedra/ball.py	2014-09-22 12:28:50 +0000
+++ examples/polyhedra/ball.py	2014-10-22 10:13:18 +0000
@@ -3,8 +3,8 @@
 
 m = PolyhedraMat()
 m.density = 2600 #kg/m^3 
-m.Ks = 20000
-m.Kn = 1E6 #Pa
+m.young = 1E6 #Pa
+m.poisson = 20000/1E6
 m.frictionAngle = 0.6 #rad
 
 maxLoad = 3E6

=== modified file 'examples/polyhedra/free-fall.py'
--- examples/polyhedra/free-fall.py	2014-09-22 12:28:50 +0000
+++ examples/polyhedra/free-fall.py	2014-10-22 10:13:18 +0000
@@ -3,20 +3,20 @@
 
 gravel = PolyhedraMat()
 gravel.density = 2600 #kg/m^3 
-gravel.Ks = 20000
-gravel.Kn = 1E7 #Pa
+gravel.young = 1E7 #Pa
+gravel.poisson = 20000/1E7
 gravel.frictionAngle = 0.5 #rad
 
 steel = PolyhedraMat()
 steel.density = 7850 #kg/m^3 
-steel.Ks = 10*gravel.Ks
-steel.Kn = 10*gravel.Kn
+steel.young = 10*gravel.young
+steel.poisson = gravel.poisson
 steel.frictionAngle = 0.4 #rad
 
 rubber = PolyhedraMat()
 rubber.density = 1000 #kg/m^3 
-rubber.Ks = gravel.Ks/10
-rubber.Kn = gravel.Kn/10
+rubber.young = gravel.young/10
+rubber.poisson = gravel.poisson
 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)))

=== modified file 'examples/polyhedra/irregular.py'
--- examples/polyhedra/irregular.py	2014-09-22 12:28:50 +0000
+++ examples/polyhedra/irregular.py	2014-10-22 10:13:18 +0000
@@ -16,8 +16,8 @@
 
 m = PolyhedraMat()
 m.density = 2600 #kg/m^3 
-m.Ks = 20000
-m.Kn = 1E6 #Pa
+m.young = 1E6 #Pa
+m.poisson = 20000/1E6
 m.frictionAngle = 0.6 #rad
 
 O.bodies.append(utils.wall(0,axis=2,sense=1, material = m))

=== modified file 'examples/polyhedra/wall.py'
--- examples/polyhedra/wall.py	2014-09-22 12:28:50 +0000
+++ examples/polyhedra/wall.py	2014-10-22 10:13:18 +0000
@@ -5,8 +5,8 @@
 
 m = PolyhedraMat()
 m.density = 1000  
-m.Ks = 5E6
-m.Kn = 5E8 
+m.young = 5E8 
+m.poisson = 5E6/5E8
 m.frictionAngle = 0.7 
 
 size = 0.1;

=== modified file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp	2014-10-20 08:28:15 +0000
+++ pkg/dem/Polyhedra.cpp	2014-10-22 10:13:18 +0000
@@ -13,7 +13,7 @@
 
 #define _USE_MATH_DEFINES
 
-YADE_PLUGIN(/* self-contained in hpp: */ (Polyhedra) (PolyhedraGeom) (Bo1_Polyhedra_Aabb) (PolyhedraPhys) (PolyhedraMat) (Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys) (Law2_PolyhedraGeom_PolyhedraPhys_Volumetric)
+YADE_PLUGIN(/* self-contained in hpp: */ (Polyhedra) (PolyhedraGeom) (Bo1_Polyhedra_Aabb) (PolyhedraPhys) (PolyhedraMat) (Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys) (Ip2_FrictMat_PolyhedraMat_FrictPhys) (Law2_PolyhedraGeom_PolyhedraPhys_Volumetric)
 	/* some code in cpp (this file): */ 
 	#ifdef YADE_OPENGL
 		(Gl1_Polyhedra) (Gl1_PolyhedraGeom) (Gl1_PolyhedraPhys)
@@ -398,16 +398,22 @@
 	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 Kna 	= mat1->young;
+	Real Knb 	= mat2->young;
+	Real Ksa 	= mat1->young*mat1->poisson;
+	Real Ksb 	= mat2->young*mat2->poisson;
 	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);
 };
 
+void Ip2_FrictMat_PolyhedraMat_FrictPhys::go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction){
+	const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(pp1);
+	const shared_ptr<PolyhedraMat>& mat2 = YADE_PTR_CAST<PolyhedraMat>(pp2);
+	Ip2_FrictMat_FrictMat_FrictPhys().go(mat1,mat2,interaction);
+}
+
 //**************************************************************************************
 #if 1
 Real Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::getPlasticDissipation() {return (Real) plasticDissipation;}

=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp	2014-10-20 08:28:15 +0000
+++ pkg/dem/Polyhedra.hpp	2014-10-22 10:13:18 +0000
@@ -25,6 +25,10 @@
 #include <CGAL/convex_hull_3.h>
 #include <CGAL/Tetrahedron_3.h>
 #include <CGAL/linear_least_squares_fitting_3.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/squared_distance_3.h>
 
 #include<time.h>
 
@@ -36,6 +40,7 @@
 typedef CGAL::Polyhedron_3<K>	Polyhedron;
 typedef CGAL::Delaunay_triangulation_3<K> Triangulation;
 typedef K::Point_3 CGALpoint;
+typedef K::Triangle_3 CGALtriangle;
 typedef K::Vector_3 CGALvector;
 typedef CGAL::Aff_transformation_3<K> Transformation;
 typedef K::Segment_3 Segment;
@@ -43,6 +48,8 @@
 typedef CGAL::Plane_3<K> Plane;
 typedef CGAL::Line_3<K> Line;
 typedef CGAL::Origin CGAL_ORIGIN;
+typedef CGAL::AABB_tree<CGAL::AABB_traits<K,CGAL::AABB_triangle_primitive<K,std::vector<Triangle>::iterator>>> CGAL_AABB_tree;
+
 
 //**********************************************************************************
 class Polyhedra: public Shape{
@@ -147,36 +154,29 @@
 
 //***************************************************************************
 /*! Elastic material */
-class PolyhedraMat: public Material{
+class PolyhedraMat: public FrictMat{
 	public:
-		 PolyhedraMat(Real N, Real S, Real F){Kn=N; Ks=S; frictionAngle=F;};
 		 Real GetStrength(){return strength;};
 	virtual ~PolyhedraMat(){};
-	YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraMat,Material,"Elastic material with Coulomb friction.",
-		((Real,Kn,1e8,,"Normal 'stiffness' (N/m3 for Law_.._Volumetric, N/m for Law2_.._Simple)."))
-		((Real,Ks,1e5,,"Shear stiffness (N/m)."))
-		((Real,frictionAngle,.5,,"Contact friction angle (in radians)."))
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraMat,FrictMat,"Elastic material with Coulomb friction.",
 		((bool,IsSplitable,0,,"To be splitted ... or not"))
 		((Real,strength,100,,"Stress at whis polyhedra of volume 4/3*pi [mm] breaks.")),
 		/*ctor*/ createIndex();
 	);
-	REGISTER_CLASS_INDEX(PolyhedraMat,Material);
+	REGISTER_CLASS_INDEX(PolyhedraMat,FrictMat);
 };
 REGISTER_SERIALIZABLE(PolyhedraMat);
 
 //***************************************************************************
-class PolyhedraPhys: public IPhys{
+class PolyhedraPhys: public FrictPhys{
 	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")),
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraPhys,FrictPhys,"Simple elastic material with friction for volumetric constitutive laws",
+		/*attrs*/
+		,
 		/*ctor*/ createIndex();	
 	);
-	REGISTER_CLASS_INDEX(PolyhedraPhys,IPhys);
+	REGISTER_CLASS_INDEX(PolyhedraPhys,FrictPhys);
 };
 REGISTER_SERIALIZABLE(PolyhedraPhys);
 
@@ -239,6 +239,17 @@
 };
 REGISTER_SERIALIZABLE(Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys);
 
+class Ip2_FrictMat_PolyhedraMat_FrictPhys: public IPhysFunctor{
+	public:
+		virtual void go(const shared_ptr<Material>& b1,
+			const shared_ptr<Material>& b2,
+			const shared_ptr<Interaction>& interaction);
+	FUNCTOR2D(FrictMat,PolyhedraMat);
+	YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_PolyhedraMat_FrictPhys,IPhysFunctor,"",		
+	);
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_PolyhedraMat_FrictPhys);
+
 //***************************************************************************
 /*! Calculate physical response based on penetration configuration given by PolyhedraGeom. */
 

=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp	2014-10-20 08:28:15 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp	2014-10-22 10:13:18 +0000
@@ -5,21 +5,22 @@
 
 #include"Polyhedra.hpp"
 #include"Polyhedra_Ig2.hpp"
+#include<pkg/dem/ScGeom.hpp>
 
 #define _USE_MATH_DEFINES
 
-YADE_PLUGIN(/* self-contained in hpp: */ (Ig2_Polyhedra_Polyhedra_PolyhedraGeom) (Ig2_Wall_Polyhedra_PolyhedraGeom) (Ig2_Facet_Polyhedra_PolyhedraGeom)	);
+YADE_PLUGIN(/* self-contained in hpp: */ (Ig2_Polyhedra_Polyhedra_PolyhedraGeom) (Ig2_Wall_Polyhedra_PolyhedraGeom) (Ig2_Facet_Polyhedra_PolyhedraGeom) (Ig2_Sphere_Polyhedra_ScGeom) );
 
 //**********************************************************************************
 /*! 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){
+bool Ig2_Polyhedra_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& shape1,const shared_ptr<Shape>& shape2,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());
+	Polyhedra* A = static_cast<Polyhedra*>(shape1.get());		
+	Polyhedra* B = static_cast<Polyhedra*>(shape2.get());
 
 	bool isNew = !interaction->geom;
 
@@ -104,12 +105,12 @@
 //**********************************************************************************
 /*! 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){
+bool Ig2_Wall_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& shape1,const shared_ptr<Shape>& shape2,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 int& PA(shape1->cast<Wall>().axis); 
+	const int& sense(shape1->cast<Wall>().sense);
 	const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;	
-	Polyhedra* B = static_cast<Polyhedra*>(cm2.get());
+	Polyhedra* B = static_cast<Polyhedra*>(shape2.get());
 
 	bool isNew = !interaction->geom;
 
@@ -170,13 +171,13 @@
 //**********************************************************************************
 /*! 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){
+bool Ig2_Facet_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& shape1,const shared_ptr<Shape>& shape2,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());
+	Facet*   A = static_cast<Facet*>(shape1.get());	
+	Polyhedra* B = static_cast<Polyhedra*>(shape2.get());
 
 	bool isNew = !interaction->geom;
 
@@ -267,15 +268,137 @@
 }
 
 //**********************************************************************************
-/*! Create Polyhedra (collision geometry) from colliding Polyhedron and Wall. */
+/*! Create Polyhedra (collision geometry) from colliding Polyhedron and Sphere. */
+
+bool Ig2_Sphere_Polyhedra_ScGeom::go(const shared_ptr<Shape>& shape1,const shared_ptr<Shape>& shape2,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;
+
+	Sphere* A = static_cast<Sphere*>(shape1.get());	
+	Real radius = A->radius;
+	Real r2 = radius*radius;
+	Polyhedra* B = static_cast<Polyhedra*>(shape2.get());
+
+	bool isNew = !interaction->geom;
+	shared_ptr<ScGeom> geom;
+
+	// use polyhedron surface triangulation
+	const std::vector<int> faceTri = B->GetSurfaceTriangulation();
+	// get vertices in global coordinate system
+	std::vector<Vector3r> pts(B->v);
+	for (int i=0; i<pts.size(); i++) pts[i] = se32.position + se32.orientation*pts[i];
+
+	//******************************************************************************
+	// find the closest point of polyhedron to the sphere center, see following paper for notation
+	// http://automatica.dei.unipd.it/public/Schenato/PSC/2010_2011/gruppo4-Building_termo_identification/IdentificazioneTermodinamica20062007/eggpsc/Simulatore%20termico/Generazione%20Mesh/1995%20-%203D%20Distance%20from%20a%20Point%20to%20a%20Triangle.pdf
+	bool isInside = true; // if center is inside polyhedron
+	const Vector3r& p0 = se31.position;
+	const Vector3r& center = p0;
+	Vector3r closest; // closest point to be found
+	Real dst2min = DBL_MAX; // minimal squared distance (large number initially)
+	// auxiliary value
+	Vector3r p1,p2,p3,e1,e2,e3,e,n,p0a,p10,p20,p30,v1,v2,v3,pa,pb,r,p0aa,pp,ppp;
+	Real dst,f1,f2,f3,t,e1n,e2n,e3n,en,p10n,p20n,p30n,o1,o2,o3,dst2,edst2,edst2min;
+	PointTriangleRelation rel=none,relTemp,relTemp2;
+	for (int i=0; i<faceTri.size()/3; i++) { // iterate over all triangles...
+		// triangle vertices
+		p1 = pts[faceTri[3*i+0]];
+		p2 = pts[faceTri[3*i+1]];
+		p3 = pts[faceTri[3*i+2]];
+		// triangle edges vectors
+		e1 = p2-p1; e1n = e1.norm();
+		e2 = p3-p2; e2n = e2.norm();
+		e3 = p1-p3; e3n = e3.norm();
+		n = (e1.cross(-e3)).normalized(); // tirangle outer normal
+		dst = (p0-p1).dot(n); // oriented distance of p0 to triangle plane
+		if (dst>0) isInside = false; // p0 lies in positive halfspace of triangle, cannot be inside
+		p0a = p0 - dst*n; // p0 projected to triangle plane
+		p10 = p0a-p1;
+		p20 = p0a-p2;
+		p30 = p0a-p3;
+		// inside/outside with respect to individual edges
+		o1 = (p10.cross(p20)).dot(n);
+		o2 = (p20.cross(p30)).dot(n);
+		o3 = (p30.cross(p10)).dot(n);
+		if (o1>0 && o2>0 && o3>0) { // p0a is inside triangle
+			pp = p0a;
+			relTemp = inside;
+		} else {
+			edst2min = DBL_MAX;
+			for (int j=0; j<3; j++) { // iterate over 3 edges end find the closest point of them
+				pa = j==0? p1 : (j==1? p2 : p3);
+				pb = j==0? p2 : (j==1? p3 : p1);
+				r = (((pb-p0a).cross(pa-p0a)).cross(pb-pa)).normalized();
+				p0aa = p0a + (pa-p0a).dot(r)*r; // projection of p0a on the edge
+				e = pb-pa;
+				en = e.norm();
+				e /= en;
+				t = (p0aa-pa).dot(e)/en; // parameter of edge line
+				if      (t<0.) { ppp = pa; relTemp2=vertex; }
+				else if (t>1.) { ppp = pb; relTemp2=vertex; }
+				else { ppp = p0aa; relTemp2=edge; }
+				edst2 = (ppp-p0).squaredNorm();
+				if (edst2 < edst2min) {
+					edst2min = edst2;
+					pp = ppp;
+					relTemp = relTemp2;
+				}
+			}
+		}
+		dst2 = (pp-p0).squaredNorm();
+		if (dst2 < dst2min) { // dst2 is now the minimal squared distance and pp is now the closest found point so far
+			dst2min = dst2;
+			closest = pp;
+			rel = relTemp;
+		}
+	}
+	//******************************************************************************
+
+	if (isNew && !(isInside || dst2min<=r2) ) return false;
+
+	if (isNew) {
+		geom = shared_ptr<ScGeom>(new ScGeom());
+		geom->radius1 = geom->radius2 = radius;
+		interaction->geom = geom;
+	} else {
+		geom = YADE_PTR_CAST<ScGeom>(interaction->geom);
+	}
+
+	Real penetrationDepth;
+	Vector3r normal,contactPoint;
+	if (isInside) { // some artificial tricks for unlikely case
+		normal = se32.position - center;
+		dst = normal.norm();
+		penetrationDepth = radius;
+		normal /= dst;
+		contactPoint = center + normal*(.5*radius);
+	} else { 
+		if (rel==none) { LOG_FATAL("TODO"); }
+		Real coeff = rel==edge? edgeCoeff : (rel==vertex? vertexCoeff : 1.0);
+		normal = closest - center;
+		dst = normal.norm();
+		normal /= dst;
+		penetrationDepth = radius - dst;
+		contactPoint = center + normal*(radius-.5*penetrationDepth);
+		penetrationDepth *= coeff;
+	}
+	geom->normal = normal;
+	geom->penetrationDepth = penetrationDepth;
+	geom->contactPoint = contactPoint;
+
+	geom->precompute(state1,state2,scene,interaction,normal,isNew,shift2,false/*avoidGranularRatcheting only for sphere-sphere*/);
+
+	return true;
+}
 
 /*
-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){
+bool Ig2_Sphere_Polyhedra_PolyhedraGeom::go(const shared_ptr<Shape>& shape1,const shared_ptr<Shape>& shape2,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;
+	Polyhedra* B = static_cast<Polyhedra*>(shape2.get());
+	const Real& r=shape1->cast<Sphere>().radius;
 
 	//cout << "Sphere x Polyhedra" << endl;
 

=== modified file 'pkg/dem/Polyhedra_Ig2.hpp'
--- pkg/dem/Polyhedra_Ig2.hpp	2014-10-15 06:44:01 +0000
+++ pkg/dem/Polyhedra_Ig2.hpp	2014-10-22 10:13:18 +0000
@@ -12,8 +12,8 @@
 {
 	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);
-		virtual bool goReverse(	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) { return go(cm1,cm2,state2,state1,-shift2,force,c); }
+		virtual bool go(const shared_ptr<Shape>& shape1, const shared_ptr<Shape>& shape2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+		virtual bool goReverse(	const shared_ptr<Shape>& shape1, const shared_ptr<Shape>& shape2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c) { return go(shape1,shape2,state2,state1,-shift2,force,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");	
@@ -28,7 +28,7 @@
 {
 	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);
+		virtual bool go(const shared_ptr<Shape>& shape1, const shared_ptr<Shape>& shape2, 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");	
@@ -43,7 +43,7 @@
 {
 	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);
+		virtual bool go(const shared_ptr<Shape>& shape1, const shared_ptr<Shape>& shape2, 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");	
@@ -54,12 +54,31 @@
 
 //***************************************************************************
 /*! Create Polyhedra (collision geometry) from colliding Sphere & Polyhedra. */
+class Ig2_Sphere_Polyhedra_ScGeom: public IGeomFunctor
+{
+	public:
+		enum PointTriangleRelation { inside, edge, vertex, none };
+		virtual ~Ig2_Sphere_Polyhedra_ScGeom(){};
+		virtual bool go(const shared_ptr<Shape>& shape1, const shared_ptr<Shape>& shape2, 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_ATTRS(Ig2_Sphere_Polyhedra_ScGeom,IGeomFunctor,"Create/update geometry of collision between Sphere and Polyhedra",
+			((Real,edgeCoeff,1.0,,"TODO"))
+			((Real,vertexCoeff,1.0,,"TODO"))
+		);	
+		DECLARE_LOGGER;	
+	private:
+};
+REGISTER_SERIALIZABLE(Ig2_Sphere_Polyhedra_ScGeom);
+
+//***************************************************************************
+/*! 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);
+		virtual bool go(const shared_ptr<Shape>& shape1, const shared_ptr<Shape>& shape2, 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");	

=== modified file 'py/_polyhedra_utils.cpp'
--- py/_polyhedra_utils.cpp	2014-10-15 06:44:01 +0000
+++ py/_polyhedra_utils.cpp	2014-10-22 10:13:18 +0000
@@ -100,7 +100,7 @@
 			Real density=b->state->mass/p->GetVolume();
 			//get equivalent radius and use same equation as for sphere
 			Real equi_radius=pow(p->GetVolume()/((4./3.)*Mathr::PI),1./3.);				
-			dt=min(dt,equi_radius/sqrt(ebp->Kn*equi_radius/density));
+			dt=min(dt,equi_radius/sqrt(ebp->young*equi_radius/density));
 		}
 	}
 	if (dt==std::numeric_limits<Real>::infinity()) {

=== modified file 'py/utils.py'
--- py/utils.py	2014-10-09 20:32:40 +0000
+++ py/utils.py	2014-10-22 10:13:18 +0000
@@ -393,7 +393,7 @@
 def polyhedron(vertices,dynamic=True,fixed=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):
 	"""Create polyhedron with given parameters.
 
-	:param [Vector3,Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.
+	:param [[Vector3]] vertices: coordinates of vertices in the global coordinate system.
 
 	See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""
 	b=Body()