polyhedral code, segmentation fault & memory leak

[Sway]

YADE version:1.11.1
I simulated a packing with Polyhedra shape recently. When running about 50,000 (not constant ) steps, a bug occurred casting information as follow:
    Error in `/usr/bin/python': malloc(): memory corruption (fast): 0x00007fc6780762d8
So give me some hints to find the reason, please.

gdb:

Program received signal SIGABRT, Aborted.
[Switching to Thread 0x7f3257fff700 (LWP 3066)]
0x00007f32b7760bb9 in __GI_raise (sig=sig@entry=6)
    at ../nptl/sysdeps/unix/sysv/linux/raise.c:56
56 ../nptl/sysdeps/unix/sysv/linux/raise.c: No such file or directory.
(gdb) bt
#0 0x00007f32b7760bb9 in __GI_raise (sig=sig@entry=6)
    at ../nptl/sysdeps/unix/sysv/linux/raise.c:56
#1 0x00007f32b7763fc8 in __GI_abort () at abort.c:89
#2 0x00007f32b779de14 in __libc_message (do_abort=do_abort@entry=1,
    fmt=fmt@entry=0x7f32b78ac668 "*** Error in `%s': %s: 0x%s ***\n")
    at ../sysdeps/posix/libc_fatal.c:175
#3 0x00007f32b77a8b77 in malloc_printerr (action=<optimized out>,
    str=0x7f32b78aca08 "malloc(): memory corruption (fast)",
    ptr=<optimized out>) at malloc.c:4996
#4 0x00007f32b77ab884 in _int_malloc (av=0x7f3250000020, bytes=48)
    at malloc.c:3359
#5 0x00007f32b77ad230 in __GI___libc_malloc (bytes=48) at malloc.c:2891
#6 0x00007f32b1639f2d in operator new(unsigned long) ()
   from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#7 0x00007f32b441cb23 in void CGAL::In_place_list<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > >, false, std::allocator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > > > >::insert<CGAL::internal::In_place_list_const_iterator<CGAL::HalfedgeDS_in_place_list_vert---Type <return> to continue, or q <return> to quit---
ex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > >, std::allocator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > > > > >(CGAL::internal::In_place_list_iterator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > >, std::allocator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > > > >, CGAL::internal::In_place_list_const_iterator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > >, std::allocator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > > > >, CGAL::internal::In_place_list_const_iterator<CGAL::HalfedgeDS_in_---Type <return> to continue, or q <return> to quit---
place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > >, std::allocator<CGAL::HalfedgeDS_in_place_list_vertex<CGAL::I_Polyhedron_vertex<CGAL::HalfedgeDS_vertex_base<CGAL::HalfedgeDS_list_types<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Point_3<CGAL::Epick> > > > > >) ()
   from .../temp-install/lib/x86_64-linux-gnu/yade-1.11.1/libyade.so
#8 0x00007f32b442abc4 in CGAL::HalfedgeDS_list<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >::HalfedgeDS_list(CGAL::HalfedgeDS_list<CGAL::Epick, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> > const&) ()
   from .../temp-install/lib/x86_64-linux-gnu/yade-1.11.1/libyade.so
#9 0x00007f32b49118e2 in Simplify(CGAL::Polyhedron_3<CGAL::Epick, CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator<int> >, double) ()
   from .../temp-install/lib/x86_64-linux-gnu/yade-1.11.1/libyade.so
#10 0x00007f32b4915e21 in Polyhedron_Plane_intersection(CGAL::Polyhedron_3<CGAL::Epick, CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator<int> >, CGAL::Plane_3<CGAL::Epick>, CGAL::Point_3<CGAL::Epick>, CGAL::Point_3<CGAL::Epick>) ()
   from .../temp-install/lib/x86_64-linux-gnu/yade-1.11.1/libyade.so
#11 0x00007f32b48ea0f8 in Ig2_Wall_Polyhedra_PolyhedraGeom::go(boost::shared_ptr<Shape> const&, boost::shared_ptr<Shape> const&, State const&, State const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, bool const&, boost::shared_ptr<Interaction> const&) ()

Something wrong in function "Simplify" within Polyhedra_support.cpp:

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;
}

[Jan Stránský (honzik)]

Hello,

please always attach a minimal working script reproducing the error.

Thanks
Jan
Dne 11.10.2014 16:40 "Sway" <email address hidden> napsal(a):

> Public bug reported:
>
> YADE version:1.11.1
> I simulated a packing with Polyhedra shape recently. When running about
> 50,000 steps, a bug occurred casting information as follow:
> Error in `/usr/bin/python': malloc(): memory corruption (fast):
> 0x00007fc6780762d8
> So give me some hints to find the reason, please.
>
> ** Affects: yade
> Importance: Undecided
> Status: New
>
> --
> You received this bug notification because you are a member of Yade
> developers, which is subscribed to Yade.
> https://bugs.launchpad.net/bugs/1380103
>
> Title:
> memory leak?
>
> Status in Yet Another Dynamic Engine:
> New
>
> Bug description:
> YADE version:1.11.1
> I simulated a packing with Polyhedra shape recently. When running about
> 50,000 steps, a bug occurred casting information as follow:
> Error in `/usr/bin/python': malloc(): memory corruption (fast):
> 0x00007fc6780762d8
> So give me some hints to find the reason, please.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/yade/+bug/1380103/+subscriptions
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>

[Bruno Chareyre (bruno-chareyre)]

Dear Sway, for your information, attaching files to emails does not make it through launchpad answers.
So in case you attached a script as suggested by Jan we can't see it at the moment. There must be a button to upload things, but you can also paste your script directly in a message.

[Sway (zhswee)]

Hi,
There is a script to reproduce the bug. And the attachment file 2.dat includes vertices of particles. To reproduce the bug quickly, I have set the time step to 1e-5s.

Here is the script:

from yade._polyhedra_utils import *
from yade.wrapper import *
from yade import plot, polyhedra_utils
#from packingutils import *
import numpy as np

import os

gravel = PolyhedraMat()
gravel.density = 2600 #kg/m^3
gravel.Ks = 1.0e9
gravel.Kn = 5.0e10 #Pa
gravel.frictionAngle = 0.0 #rad

wall = PolyhedraMat()
wall.density = 2600 #kg/m^3
wall.Ks = 1.0e7
wall.Kn = 1.0e16 #Pa
wall.frictionAngle = 0.0 #rad
####input data
mn,mx=Vector3(0,0,0),Vector3(0.1,0.075,0.05) # corners of the initial packing

## create materials for spheres and plates
##O.materials.append(PolyhedraMat(frictionAngle=0.0,density=2600,Kn=5.0e10,Ks=1.0e9,label='gravel'))
##O.materials.append(PolyhedraMat(frictionAngle=0.0,density=2600,Kn=1e16,Ks=1e7,label='walls'))

###some functions###
def gettop(b):##maximum of y direction
 vv=[b.state.se3[1].toRotationMatrix()*v+b.state.se3[0] for v in b.shape.v]
 vvy=[i[1] for i in vv] ##y
 return max(vvy)

def loading():
 #O.bodies.append(utils.wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2))
 global plate1,plate2
 # without this line, the plate variable would only exist inside this function
 plate1=O.bodies[3]# the last particles is the plate
 #plate1.state.pos=(0,max([gettop(b) for b in O.bodies if isinstance(b.shape,Polyhedra)]),0)
 plate1.state.vel=(0,-0.1,0)
 #if plate1.state.pos[1]<0.05:
  #packing.dead=True
  #O.stopAtIter=O.iter+1

def outputPos(filename):
 fobj=open(filename,'a+')
 for b in O.bodies:
  if isinstance(b.shape,Polyhedra):
   vv=[b.state.se3[1].toRotationMatrix()*v+b.state.se3[0] for v in b.shape.v]
   vvv=[[v[0],v[1],v[2]] for v in vv]
   print>>fobj,vvv
        fobj.close()

def readData(filename):
    if os.path.isfile(filename):
 fobj=open(filename)

 vv=list()
 while 1:
            line=fobj.readline()
            if not line:break
     line=line[2:-3]
     temp_v=list()
     s=line.split('], [')
     for i in s:
  sv=i.split(", ")
  temp_v.append([float(j) for j in sv])
     vv.append(temp_v)
 fobj.close()
 return vv
    else:
 print "Error:The file is not found!"
 return False

def addBodies(vv):
 for v in vv:
  obj1=polyhedra_utils.polyhedra(gravel,v=v,color=np.random.rand(3))
  O.bodies.append(obj1)

## create walls around the packing
O.bodies.append(utils.wall(0,axis=0,sense=1, material = wall))#left x
O.bodies.append(utils.wall(mx[0],axis=0,sense=-1, material = wall))#right x
O.bodies.append(utils.wall(0,axis=1,sense=1, material = wall))#bottom y
O.bodies.append(utils.wall(mx[1],axis=1,sense=-1, material = wall))#top y
O.bodies.append(utils.wall(0,axis=2,sense=1, material = wall))#back z
O.bodies.append(utils.wall(mx[2],axis=2,sense=-1, material = wall))#front z
##fill particles
vv=readData("2.dat")
addBodies(vv)
#polyhedra_utils.fillBox(mn, mx,gravel,sizemin=[0.003,0.003,0.003],sizemax=[0.008,0.008,0.008],seed=10002)

def checkUnbalancedI():
    print "iter %d, time elapsed %f, time step %.5e, unbalanced forces = %.5f"%(O.iter, O.realtime, O.dt, utils.u...

[Sway (zhswee)]

Hi,
Another problem is that the routine can't run when I use the data in the attachment 1.dat.
Sway

[Jan Elias (elias-j)]

Dear Sway,

I am the "developer" of the polyhedral code. The problem you are experiencing is related to CGAL. When two particles come into the contact, there is an algorithm that computes the volume of the intersection. The volume is computed based on two runs of CGAL method Convex_Hull (you can find it in paper Simulation of railway ballast using crushable polyhedral particles. Powder Technology, 246, pp. 458-465). However, that method may result in error as occured to you due to inconvenient input point (like points almost in plane). It would be nice if CGAL just return warning and one can handle it then, but instead it results in termination of the process.

In the YADE polyhedral code, there are few check points that try to identify the wrong input and handle it as no contact. But they are based on some chosen level, at which the points are rejected. You can set those values differently to be more strict in rejection of the inputs, but as a result you will get problems with stability (your contacts will imediatelly jump from no penetration to high penetration depth). As you now realize, these problems are typical for materials with high stiffness. I used 2x10^13 for gravel, I had these problems only few times. But for higher stiffness (you use 10^16, what material do you simulate), it will become more and more frequent.

Well, I do not have much help for you. You can try to decrease stiffness of your material (I tried on your code and it works just fine). Or you can play with controlling levels (in PolyhedraSupport.cpp), they are currently set to:

//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 2E-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-19 //19
//FIND_NORMAL_LIMIT - to determine which facet of intersection belongs to which polyhedron
#define FIND_NORMAL_LIMIT 1E-40

Best regards, Jan Elias

[Sway (zhswee)]

Dear Jan Elias,

Thank you for your detailed answer! Yeah, the problem is related to CGAL. especially function convex_hull_3. Do you mean that if the intersection is very closed to a plane then a crash occurs in convex_hull_3? I changed the ConvexHull as this:

Polyhedron ConvexHull(vector<CGALpoint> &planes){
 Polyhedron Int;
 CGAL::Object obj;//not sure the return of convex_hull_3 is a polyhedron exactly.
 if (planes.size() > 3){//are those points collinear?
   CGAL::convex_hull_3(planes.begin(),planes.end(),obj);
    if ( CGAL::assign (Int, obj)) {
      std::transform( Int.facets_begin(), Int.facets_end(), Int.planes_begin(),Plane_equation());
      if (CGAL::is_strongly_convex_3(Int)) return Int;//ConvexHull is a polyhedron
    }
 }
 //if (planes.size()>3) CGAL::convex_hull_3(planes.begin(), planes.end(), Int);
 //cout << "-C"<< endl; cout.flush();
 return Int;
}

If the convex hull is not a polyhedron, then return empty. Unfortunately, this didn't give a help. I am not sure whether the bug results from convex_hull_3 or not.
As you say, maybe I can decrease the stiffness to get a larger intersection which is helpful to avoid the crash.

Best,
Sway

[Bruno Chareyre (bruno-chareyre)]

Hello Jan,
Did you try with a different kernel (e.g.
Exact_predicates_exact_constructions_kernel)?
Did you report the issue to CGAL devs?

B

[Jan Elias (elias-j)]

Sway,
the crashes occur at row CGAL::convex_hull_3(planes.begin(),planes.end(),obj); - inside CGAL. Any check that you perform after it is useless. It must be done before. It might be also other issue, but all the problems I had with the polyhedral code were only there.

Bruno,
both of the suggestions are good. Yes, I tried different kernel. It has been already two years since I tried it so I am not sure why I didn't use it. I think it caused trouble in other parts of the code, but I do not really remember. I also think that it was much slower. Initially, I didn't know where the problem is. The bug is very hard to reproduce. Exact_predicates_exact_constructions_kernel would probably work.
I did not reported that to CGAL developers. I should do it, you are right. For very long time (one year), I wasn't sure about the source of the problem. When I finally realized where is the bug, I just tried to fix it somehow. It worked for me, for my stiffness. I will try to force myslef to report it to CGAL. Thanks.

Jan

[Sway (zhswee)]

Hello Jan,

The crashes occur at row "CGAL::convex_hull_3(planes.begin(),planes.end(),obj); " in the most case. See the original description I provided and the crash occurs at " for (Polyhedron::Edge_iterator hei = P.edges_begin(); hei!=P.edges_end(); ++hei)" in function SImpify().

I'm sure the kernel Exact_predicates_exact_constructions_kernel would make the routine very slow.

Decreasing the stiffness is not a good workaround for me, say the shaping particle in the attachment 1.dat, the crash is always there.

Best,
Sway

[Sway (zhswee)]

Hello Jan and Bruno,

A good news is that the bug seems to be fixed. I think the bug derives from the dual polyhedron. In the present polyhedral code, once we find a point in the intersection, we move the origin to the point. The problem is that there may be a plane which is very closed to the new origin, so the dual point of this plane will be far away from the new origin. This dual point is not a good input of convex_hull_3. A available workaround, I think, is scaling the potential intersection manually. We can scale the real intersection back to its real size after getting the intersection. So I added several lines as follows:

//after translating the polyhedra to the new coordinate system
Transformation scale(CGAL::SCALING, 10);//scale 10 times
Transformation scale_back(CGAL::SCALING,0.1);

//scaling A,B
std::transforms(..., scale);
...
std::transforms(...,scale_back);//before translate the new origin to the old origin.

-----------------------------------------------------------------------------------
This method can avoid the crash for me even if I adopt a high stiffness like 1e16.

Thanks!
Best,
Sway

[Jan Elias (elias-j)]

Sway,
your solution looks reasonable, two years ago I remember some scaling attempts as well. But then I moved somewhere else where finally everything worked for me. I am suprised by the position of the error, I have never seen it there. But I still do not understand what exactly the problem is.
Anyway, this scaling might help you also with the high stiffness. As I wrote before, there are several check points that returns no contact if the overlap is extremely small to prevent CGAL error. But for high stiffnesses, it would result in unstable behavior. When you scale the whole problem, you actually decrease also the limits for return and help to restore stable situation.
Thanks for solving it, I do not have much time now to go again deeply into it. I am very happy that someone uses my code, understand it and even improve it. Good luck with your work.
Jan