yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3816
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Thu 2016-03-24 23:08:07 +0100
message:
  Fix division by zero crash in Polyhedra
  
  One need to check, whether maxFS is not equal
  zero. In this case ratio=0 and is a divider
  for shearforce, which leads to NaNs.
  
  Clean some headers and indentions.
  
  One need also to implement a "fool proofing" from adding
  several equal IGs in InteractionLoop.
modified:
  pkg/dem/Polyhedra.cpp
  pkg/dem/Polyhedra.hpp
  pkg/dem/Polyhedra_Ig2.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	2016-03-08 09:18:20 +0000
+++ pkg/dem/Polyhedra.cpp	2016-03-24 22:08:07 +0000
@@ -3,13 +3,8 @@
 
 #ifdef YADE_CGAL
 
-#include<core/Interaction.hpp>
-#include<core/Omega.hpp>
-#include<core/Scene.hpp>
-#include<core/State.hpp>
-#include<pkg/common/ElastMat.hpp>
-#include<pkg/common/Aabb.hpp>
-#include"Polyhedra.hpp"
+#include <pkg/common/ElastMat.hpp>
+#include "Polyhedra.hpp"
 
 #define _USE_MATH_DEFINES
 
@@ -355,7 +350,7 @@
 
 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);
@@ -443,45 +438,55 @@
 
 		//erase the interaction when aAbB shows separation, otherwise keep it to be able to store previous separating plane for fast detection of separation 
 		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])  {
+		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;
 		}
-			
+		
 		//zero penetration depth means no interaction force
-		if(!(contactGeom->equivalentPenetrationDepth > 1E-18) || !(contactGeom->penetrationVolume > 0)) {phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.); return true;} 
+		if(!(contactGeom->equivalentPenetrationDepth > 1E-18) || !(contactGeom->penetrationVolume > 0)) {
+			phys->normalForce = Vector3r(0.,0.,0.);
+			phys->shearForce = Vector3r(0.,0.,0.);
+			return true;
+		}
+    
 		Real prop = std::pow(contactGeom->penetrationVolume,volumePower);
 		Vector3r normalForce=contactGeom->normal*prop*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);
+		//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);
+		const Real maxFs = normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
+
+		if(shearForce.squaredNorm() > maxFs && maxFs){
+			//PFC3d SlipModel, is using friction angle. CoulombCriterion
+			const Real ratio = sqrt(maxFs) / shearForce.norm();
+
+			//Store prev force for definition of plastic slip
+			//Define the plastic work input and increment the total plastic energy dissipated
+			const Vector3r trialForce=shearForce;
+			shearForce *= ratio;
+			
+			if (scene->trackEnergy && traceEnergy) {
+				const 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);
 			}
-			// compute elastic energy as well
-			scene->energy->add(0.5*(normalForce.squaredNorm()/phys->kn+shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,true);
+		} else {
+			shearForce = Vector3r::Zero();
 		}
 
 		Vector3r F = -normalForce-shearForce;	

=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp	2016-03-08 09:18:20 +0000
+++ pkg/dem/Polyhedra.hpp	2016-03-24 22:08:07 +0000
@@ -6,16 +6,15 @@
 
 #ifdef YADE_CGAL
 
-#include<core/Shape.hpp>
-#include<core/IGeom.hpp>
-#include<core/GlobalEngine.hpp>
-#include<core/Material.hpp>
-#include<pkg/common/Aabb.hpp>
-#include<pkg/common/Dispatching.hpp>
-#include<pkg/dem/FrictPhys.hpp>
-#include<pkg/common/Wall.hpp>
-#include<pkg/common/Facet.hpp>
-#include<lib/base/openmp-accu.hpp>
+#include <core/Omega.hpp>
+#include <core/Shape.hpp>
+#include <core/Interaction.hpp>
+#include <core/Material.hpp>
+#include <pkg/dem/ScGeom.hpp>
+#include <pkg/dem/FrictPhys.hpp>
+#include <pkg/common/Wall.hpp>
+#include <pkg/common/Facet.hpp>
+#include <pkg/common/Dispatching.hpp>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Delaunay_triangulation_3.h>
@@ -30,8 +29,6 @@
 #include <CGAL/AABB_triangle_primitive.h>
 #include <CGAL/squared_distance_3.h>
 
-#include<time.h>
-
 #define likely(x)       __builtin_expect((x),1)
 #define unlikely(x)     __builtin_expect((x),0)
 

=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp	2016-03-08 09:18:20 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp	2016-03-24 22:08:07 +0000
@@ -3,9 +3,8 @@
 
 #ifdef YADE_CGAL
 
-#include"Polyhedra.hpp"
-#include"Polyhedra_Ig2.hpp"
-#include<pkg/dem/ScGeom.hpp>
+#include "Polyhedra.hpp"
+#include "Polyhedra_Ig2.hpp"
 
 #define _USE_MATH_DEFINES
 
@@ -14,12 +13,18 @@
 //**********************************************************************************
 /*! Create Polyhedra (collision geometry) from colliding Polyhedras. */
 
-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){
-
+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*>(shape1.get());		
+	Polyhedra* A = static_cast<Polyhedra*>(shape1.get());
 	Polyhedra* B = static_cast<Polyhedra*>(shape2.get());
 
 	bool isNew = !interaction->geom;
@@ -181,7 +186,7 @@
 																					 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*>(shape1.get());	
@@ -265,7 +270,7 @@
 	bang->contactPoint=centroid;
 	bang->penetrationVolume=volume;
 	bang->equivalentPenetrationDepth=volume/area;
-	if (interaction && force) bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
+	bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
 	bang->normal=normal;
 	return true;
 }