yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12591
[Branch ~yade-pkg/yade/git-trunk] Rev 3824: Check isnan in some coordinates before calling CGAL.
------------------------------------------------------------
revno: 3824
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Thu 2016-03-24 23:08:26 +0100
message:
Check isnan in some coordinates before calling CGAL.
It crashes, if isnan.
Minor indention fixes.
modified:
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
--
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-24 22:08:26 +0000
+++ pkg/dem/Polyhedra.cpp 2016-03-24 22:08:26 +0000
@@ -481,9 +481,9 @@
} else {
shearForce = Vector3r::Zero();
}
-
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));
@@ -501,9 +501,8 @@
fprintf(fin,"A\t%e\t%e\t%e\n",A->state->pos[0],A->state->pos[1],A->state->pos[2]);
fprintf(fin,"B\t%e\t%e\t%e\n",B->state->pos[0],B->state->pos[1],B->state->pos[2]);
fprintf(fin,"centroid\t%e\t%e\t%e\n",contactGeom->contactPoint[0],contactGeom->contactPoint[1],contactGeom->contactPoint[2]);
- fclose(fin);
- */
-
+ fclose(fin);
+ */
//needed to be able to acces interaction forces in other parts of yade
phys->normalForce = normalForce;
phys->shearForce = shearForce;
=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp 2016-03-24 22:08:26 +0000
+++ pkg/dem/Polyhedra.hpp 2016-03-24 22:08:26 +0000
@@ -5,6 +5,11 @@
#ifdef YADE_CGAL
+// NDEBUG causes crashes in CGAL sometimes. Anton
+#ifdef NDEBUG
+ #undef NDEBUG
+#endif
+
#include <core/Omega.hpp>
#include <core/Shape.hpp>
#include <core/Interaction.hpp>
@@ -60,7 +65,7 @@
Vector3r GetInertia(){Initialize(); return inertia;}
vector<int> GetSurfaceTriangulation(){Initialize(); return faceTri;}
vector<vector<int>> GetSurfaces() const;
- void Initialize();
+ void Initialize();
bool IsInitialized(){return init;}
std::vector<Vector3r> GetOriginalVertices();
Real GetVolume(){Initialize(); return volume;}
=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp 2016-03-24 22:08:26 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp 2016-03-24 22:08:26 +0000
@@ -29,21 +29,30 @@
//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.);
+
+ Transformation t_rot_trans(
+ rot_mat(0,0), rot_mat(0,1), rot_mat(0,2), trans_vec[0],
+ rot_mat(1,0), rot_mat(1,1), rot_mat(1,2), trans_vec[1],
+ rot_mat(2,0), rot_mat(2,1), rot_mat(2,2), trans_vec[2], 1.
+ );
+
Polyhedron PA = A->GetPolyhedron();
std::transform( PA.points_begin(), PA.points_end(), PA.points_begin(), t_rot_trans);
+
std::transform( PA.facets_begin(), PA.facets_end(), PA.planes_begin(),Plane_equation());
-
-
+
//move and rotate 2nd the CGAL structure Polyhedron
rot_mat = (se32.orientation).toRotationMatrix();
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);
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());
@@ -53,31 +62,37 @@
interaction->geom = bang;
}else{
// use data from old interaction
- bang=YADE_PTR_CAST<PolyhedraGeom>(interaction->geom);
+ 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+shift2), bang->sep_plane);
-
+
//volume and centroid of intersection
Real volume;
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
+
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;}
-
+
+ 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);
+ Vector3r normal = FindNormal(Int, PA, PB);
+
if((se32.position+shift2-centroid).dot(normal)<0) normal*=-1;
-
+
//calculate area of projection of Intersection into the normal plane
//Real area = CalculateProjectionArea(Int, ToCGALVector(normal));
//if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
@@ -132,7 +147,11 @@
//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.);
+
+ 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());
@@ -201,22 +220,25 @@
const Se3r& se32=state2.se3;
Facet* A = static_cast<Facet*>(shape1.get());
Polyhedra* B = static_cast<Polyhedra*>(shape2.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.);
+ Transformation t_rot_trans(
+ rot_mat(0,0),rot_mat(0,1),rot_mat(0,2),trans_vec[0],
+ rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),trans_vec[1],
+ rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),trans_vec[2],1.);
Polyhedron PB = B->GetPolyhedron();
std::transform( PB.points_begin(), PB.points_end(), PB.points_begin(), t_rot_trans);
std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());
-
+
//move and rotate facet
vector<CGALpoint> v;
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]));
@@ -229,7 +251,7 @@
v[0]=v[1];
v[1]=help;
}
-
+
Real 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
@@ -242,7 +264,7 @@
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
@@ -256,22 +278,22 @@
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
Real volume;
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
if (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid))) {bang->equivalentPenetrationDepth=0; return true;}
-
+
//find normal direction
Vector3r normal = FindNormal(Int, PA, PB);
if((se32.position-centroid).dot(normal)<0) normal*=-1;
-
+
//calculate area of projection of Intersection into the normal plane
Real area = volume/1E-8;
@@ -282,6 +304,7 @@
bang->equivalentPenetrationDepth=volume/area;
bang->precompute(state1,state2,scene,interaction,normal,bang->isShearNew,shift2);
bang->normal=normal;
+
return true;
}
=== modified file 'pkg/dem/Polyhedra_Ig2.hpp'
--- pkg/dem/Polyhedra_Ig2.hpp 2016-03-24 22:08:26 +0000
+++ pkg/dem/Polyhedra_Ig2.hpp 2016-03-24 22:08:26 +0000
@@ -2,7 +2,6 @@
// https://www.vutbr.cz/www_base/gigadisk.php?i=95194aa9a
#ifdef YADE_CGAL
-
#include "Polyhedra.hpp"
//***************************************************************************
=== modified file 'pkg/dem/Polyhedra_splitter.cpp'
--- pkg/dem/Polyhedra_splitter.cpp 2016-03-24 22:08:26 +0000
+++ pkg/dem/Polyhedra_splitter.cpp 2016-03-24 22:08:26 +0000
@@ -80,7 +80,7 @@
if(!b || !b->material || !b->shape) continue;
shared_ptr<Polyhedra> p=YADE_PTR_DYN_CAST<Polyhedra>(b->shape);
shared_ptr<PolyhedraMat> m=YADE_PTR_DYN_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[b->id];
=== modified file 'pkg/dem/Polyhedra_splitter.hpp'
--- pkg/dem/Polyhedra_splitter.hpp 2016-03-24 22:08:26 +0000
+++ pkg/dem/Polyhedra_splitter.hpp 2016-03-24 22:08:26 +0000
@@ -5,9 +5,9 @@
#pragma once
#ifdef YADE_CGAL
-
+#include "Polyhedra.hpp"
#include <pkg/common/PeriodicEngines.hpp>
-#include "Polyhedra.hpp"
+
//*********************************************************************************
/* Polyhedra Splitter */
=== modified file 'pkg/dem/Polyhedra_support.cpp'
--- pkg/dem/Polyhedra_support.cpp 2016-03-24 22:07:41 +0000
+++ pkg/dem/Polyhedra_support.cpp 2016-03-24 22:08:26 +0000
@@ -3,7 +3,7 @@
#ifdef YADE_CGAL
-#include"Polyhedra.hpp"
+#include "Polyhedra.hpp"
#define _USE_MATH_DEFINES
@@ -368,7 +368,7 @@
Vector3r SolveLinSys3x3(Matrix3r A, Vector3r y){
//only system 3x3 by Cramers rule
Real 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);}
+ if (det == 0) {LOG_WARN("error in linear solver"); 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);
}
@@ -378,9 +378,10 @@
*/
Polyhedron ConvexHull(vector<CGALpoint> &planes){
Polyhedron Int;
- //cout << "C"; cout.flush();
+ for (const auto p : planes) {
+ if (isnan(p.x()) || isnan(p.y()) || isnan(p.z())) return Int;
+ }
if (planes.size()>3) CGAL::convex_hull_3(planes.begin(), planes.end(), Int);
- //cout << "-C"<< endl; cout.flush();
return Int;
}
@@ -520,7 +521,7 @@
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;
+ LOG_WARN("Error in line-plane intersection");
}
}
if(Is_inside_Polyhedron(A,inside,DISTANCE_LIMIT) && Oriented_squared_distance(B,inside)<=-lim) intersection_found = true;
@@ -564,6 +565,7 @@
//compute convex hull of it
Intersection = ConvexHull(dual_planes);
+ if (Intersection.empty()) return Intersection;
//return to original position
std::transform( Intersection.points_begin(), Intersection.points_end(), Intersection.points_begin(), transl_back);
@@ -673,10 +675,20 @@
//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
+ for(Polyhedron::Plane_iterator pi =Intersection.planes_begin(); pi!=Intersection.planes_end(); ++pi) {
+ const auto pX = -pi->a()/pi->d();
+ const auto pY = -pi->b()/pi->d();
+ const auto pZ = -pi->c()/pi->d();
+ if (isnan(pX) || isnan(pY) || isnan(pZ)) {
+ Polyhedron IntersectionEmpty;
+ return IntersectionEmpty;
+ } else {
+ 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()) return Intersection;
//return to original position
std::transform( Intersection.points_begin(), Intersection.points_end(), Intersection.points_begin(), transl_back);
@@ -792,37 +804,37 @@
//**********************************************************************************
//split polyhedra
-shared_ptr<Body> SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction, Vector3r point){
-
- direction.normalize();
-
+shared_ptr<Body> SplitPolyhedra(const shared_ptr<Body>& body, const Vector3r direction, const Vector3r point){
+ Scene* scene=Omega::instance().getScene().get();
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;
+ const Vector3r OrigPos = X->pos;
+ const Vector3r OrigVel = X->vel;
+ const 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.);
+ 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 first splitted polyhedrons
Plane B(ToCGALPoint(point-direction*SPLITTER_GAP), ToCGALVector(direction));
Polyhedron S1 = Polyhedron_Plane_intersection(PA, B, ToCGALPoint(se3.position), B.projection(ToCGALPoint(OrigPos)) - 1E-6*ToCGALVector(direction));
+
B = Plane(ToCGALPoint(point+direction*SPLITTER_GAP), ToCGALVector((-1.)*direction));
Polyhedron S2 = Polyhedron_Plane_intersection(PA, B, ToCGALPoint(se3.position), B.projection(ToCGALPoint(OrigPos)) + 1E-6*ToCGALVector(direction));
- Scene* scene=Omega::instance().getScene().get();
//scene->bodies->erase(body->id);
-
//replace original polyhedron
- A->Clear();
- for(Polyhedron::Vertex_iterator vi = S1.vertices_begin(); vi != S1.vertices_end(); vi++){ A->v.push_back(FromCGALPoint(vi->point()));};
+ 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();
@@ -830,13 +842,12 @@
X->refPos[0] = X->refPos[0]+1.;
X->inertia =A->GetInertia()*body->material->density;
X->vel = OrigVel + OrigAngVel.cross(X->pos-OrigPos);
- X->angVel = OrigAngVel;
+ X->angVel = OrigAngVel;
//second polyhedron
vector<Vector3r> v2;
- for(Polyhedron::Vertex_iterator vi = S2.vertices_begin(); vi != S2.vertices_end(); vi++){ v2.push_back(FromCGALPoint(vi->point()));};
+ 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(Real(rand())/RAND_MAX,Real(rand())/RAND_MAX,Real(rand())/RAND_MAX);
scene->bodies->insert(BP);
//set proper state variables
