yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 4171
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Mon 2014-09-22 14:28:50 +0200
message:
  enable periodic simulations with Polhedron and Tetra
modified:
  pkg/dem/Polyhedra.cpp
  pkg/dem/Polyhedra_Ig2.cpp
  pkg/dem/Polyhedra_Ig2.hpp
  pkg/dem/Tetra.cpp


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

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp	2014-07-18 18:18:50 +0000
+++ pkg/dem/Polyhedra.cpp	2014-09-22 12:28:50 +0000
@@ -220,6 +220,7 @@
 	}
 }
 
+
 //****************************************************************************************
 /* Destructor */
 
@@ -424,7 +425,8 @@
 		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])  {
+		Vector3r shift2=scene->cell->hSize*I->cellDist.cast<Real>();
+		if (A->bound->min[0] >= B->bound->max[0]+shift2[0] || B->bound->min[0]+shift2[0] >= A->bound->max[0] || A->bound->min[1] >= B->bound->max[1]+shift2[1] || B->bound->min[1]+shift2[1] >= A->bound->max[1] || A->bound->min[2] >= B->bound->max[2]+shift2[2] || B->bound->min[2]+shift2[2] >= A->bound->max[2])  {
 			return false;
 		}
 			

=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp	2014-05-16 11:58:59 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp	2014-09-22 12:28:50 +0000
@@ -34,7 +34,7 @@
 
 	//move and rotate 2nd the CGAL structure Polyhedron
 	rot_mat = (se32.orientation).toRotationMatrix();
-	trans_vec = se32.position;
+	trans_vec = se32.position + shift2;
 	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);
@@ -57,18 +57,23 @@
 
 	//find intersection Polyhedra
 	Polyhedron Int;
-	Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position), bang->sep_plane);	
+	Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position+shift2), bang->sep_plane);	
 
 	//volume and centroid of intersection
 	double volume;
 	Vector3r centroid;	
 	P_volume_centroid(Int, &volume, &centroid);
- 	if(isnan(volume) || volume<=1E-25 || volume > min(A->GetVolume(),B->GetVolume())) {bang->equivalentPenetrationDepth=0; 	return true;}
+ 	if(isnan(volume) || volume<=1E-25 || volume > min(A->GetVolume(),B->GetVolume())) {
+		bang->equivalentPenetrationDepth=0;
+		bang->penetrationVolume=min(A->GetVolume(),B->GetVolume());
+		bang->normal = (A->GetVolume()>B->GetVolume() ? 1 : -1)*(se32.position+shift2-se31.position);
+		return true;
+	}
 	if ( (!Is_inside_Polyhedron(PA, ToCGALPoint(centroid))) or (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid))))  {bang->equivalentPenetrationDepth=0; return true;}
 
 	//find normal direction
         Vector3r normal = FindNormal(Int, PA, PB);
-	if((se32.position-centroid).dot(normal)<0) normal*=-1;	
+	if((se32.position+shift2-centroid).dot(normal)<0) normal*=-1;	
 
 	//calculate area of projection of Intersection into the normal plane
 	//double area = CalculateProjectionArea(Int, ToCGALVector(normal));
@@ -93,7 +98,6 @@
 	PrintPolyhedron2File(Int,fin);
 	fclose(fin);
 	*/
-
 	return true;	
 }
 

=== modified file 'pkg/dem/Polyhedra_Ig2.hpp'
--- pkg/dem/Polyhedra_Ig2.hpp	2013-10-16 15:24:49 +0000
+++ pkg/dem/Polyhedra_Ig2.hpp	2014-09-22 12:28:50 +0000
@@ -13,6 +13,7 @@
 	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); }
 		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");	

=== modified file 'pkg/dem/Tetra.cpp'
--- pkg/dem/Tetra.cpp	2014-07-18 18:18:50 +0000
+++ pkg/dem/Tetra.cpp	2014-09-22 12:28:50 +0000
@@ -513,7 +513,7 @@
 	for (int i=0; i<4; i++) {
 		vTemp = se31.position + se31.orientation*shape1->v[i];
 		p1[i] = Point(vTemp[0],vTemp[1],vTemp[2]);
-		vTemp = se32.position + se32.orientation*shape2->v[i];
+		vTemp = se32.position + se32.orientation*shape2->v[i] + shift2;
 		p2[i] = Point(vTemp[0],vTemp[1],vTemp[2]);
 	}