core dumped - Eigen - Triangulation

Asked by Son Pham Thai

Hi YADE community,

I want to generate a pore network by using regular triangulation method implemented in YADE.

The following properties are given:

1. Radius and center coordinate of solid particle
2. Cylindrical wall (facet)
3. The material properties of the solid particle and cylindrical wall

When I run initialization(), I got the error:

python: /usr/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:408: Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::Index) [with Derived = Eigen::Matrix<double, 3, 1>; Eigen::DenseCoeffsBase<Derived, 1>::Scalar = double; Eigen::Index = long int]: Assertion `index >= 0 && index < size()' failed.
SIGSEGV/SIGABRT handler called; gdb batch file is `/tmp/yade-sQh1QT/tmp-0'
GNU gdb (Ubuntu 7.7.1-0ubuntu5~14.04.3) 7.7.1
Copyright (C) 2014 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law. Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word".
Could not attach to process. If your uid matches the uid of the target
process, check the setting of /proc/sys/kernel/yama/ptrace_scope, or try
again as the root user. For more details, see /etc/sysctl.d/10-ptrace.conf
/tmp/yade-sQh1QT/tmp-0:1: Error in sourced command file:
ptrace: Operation not permitted.
(gdb) quit
Aborted (core dumped)

# ------------------------------------------------------------ #
A mini script to reproduce the above error:
# ------------------------------------------------------------ #

# encoding: utf-8
finalFricDegree = 30
damping0=0.4
young0= 30e9
poisson0= 0.2
density0 = 2530
tensileStrength = 40e5
shearStrength = 40e6
compFricDegree = 1
youngDrying= 30e9
poissonDrying= 0.2
densityDrying = 2530*1e3
surfaceTension= 0.07274

#####################################
### setup the material properties ###
#####################################

O.materials.append(CohFrictMat(young=young0,poisson=poisson0,density=density0,frictionAngle=radians(compFricDegree), normalCohesion=tensileStrength, shearCohesion=shearStrength, momentRotationLaw=True,label='CohSpheres'))

O.materials.append(FrictMat(young=young0,poisson=poisson0,frictionAngle=0,density=0,label='frictionlessWalls'))

O.materials.append(FrictMat(young=youngDrying,poisson=poissonDrying,frictionAngle=radians(finalFricDegree),density=densityDrying,label='frictionSphere'))

##################################################
### --- reloading the sphere and wall data --- ###
##################################################
import numpy as np

#-----------------------------------#
### the spheres ### of main network #
#-----------------------------------#
particleInMainNetworkRad= [ 0.000859884 ,
 0.000861149 ,
 0.000842029 ,
 0.00085647 ]

particleInMainNetworkPos=[( 0.00317429 , 0.00550978 , 0.004819 ),
( 0.00531241 , 0.00159537 , 0.000953932 ),
( 0.00144734 , 0.00667752 , 0.00547688 ),
( 0.00521373 , 0.00203944 , 0.00291114 )]

for i in range(len(particleInMainNetworkRad)):
 O.bodies.append([utils.sphere(center=particleInMainNetworkPos[i],radius=particleInMainNetworkRad[i],material='CohSpheres')])

#----------------------------------#
### the walls ### cylinder facet ###
#----------------------------------#

### --- generate the cylindrical wall --- ###
centerCylinder=(0.004095328, 0.0041199985, 0.0023)
radiusCylinder=0.00418
heightCylinder=0.0085
segmentsNumberCylinder=5 # number of edges on the cylinder surface (>=5)
oriCylinder=Quaternion((1,0,0),0)
wallMaskCylinder=5
kwMeshes={'color':[1,0,0],'wire':True,'dynamic':False,'material':"frictionlessWalls"} # COLOR = red
O.bodies.append(geom.facetCylinder(centerCylinder, radius=radiusCylinder, height=heightCylinder, orientation=oriCylinder, segmentsNumber=segmentsNumberCylinder, wallMask=wallMaskCylinder, **kwMeshes))

#---------------------#
### Engine setup ### #
#---------------------#

gravityVector = (0.0,0,-9.81)
damping0=0.4
newton=NewtonIntegrator(damping=damping0, gravity=gravityVector)

unsat=UnsaturatedEngine(dead=1,label="unsat")
globals()['unsat'] = unsat

O.engines=[
 ForceResetter(),

 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),

 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],

  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True, label="cohesiveIp")],

  [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
   useIncrementalForm=True,
   always_use_moment_law=False,
   label='cohesiveLaw')]
 ),

 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.5),

 newton,
]
O.dt=.5*PWaveTimeStep()

################################
## Drying simulation ###
################################

from numpy import *

#set boundary conditions
unsat.bndCondIsPressure=[0,0,0,1,0,0]
unsat.bndCondValue=[0,0,0,0,0,0]

unsat.isPhaseTrapped=True # the W-phase can be disconnected from its reservoir
unsat.isInvadeBoundary=True # if False, the NW-phase can not invade the side pores
unsat.drainageFirst=True ### "If true, activate drainage first (initial saturated), then imbibition; if false, activate imbibition first (initial unsaturated), then drainage.
unsat.isDrainageActivated=True
unsat.isImbibitionActivated=False
unsat.computeForceActivated=True

unsat.initialization() # Regular triangulation, pore network generation, and initializeReservoirs() (set pore saturation for boundingCells[2] as W-reservoir and boundingCells[3] as NW-reservoir )

# ----------------------------------------------------------------- #

Looking forward to your help!
Thanks and best regards,
Son Pham

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Son Pham Thai
Solved:
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#1

Hi Son,
I'm unsure what the triangulation will do with a hollow cylinder, but most likely it's not doing anything good.
To go further please describe more basic steps, since I suppose you cases where it runs without the error. For instance, do you get the same error with a sphere packing without a cylinder?
Bruno

Revision history for this message
Son Pham Thai (pham-thai-son-987) said :
#2

Hi Bruno,

Thank you for the response!

I have tried running triangulation without appending the cylinder facet, but the error is still the same.

I have also tried replacing the cylinder facet with the boxes and still got the same error.

The error log:

negative volume for an ordinary pore (temp warning, should still be safe)

python: /usr/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:408: Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::Index) [with Derived = Eigen::Matrix<double, 3, 1>; Eigen::DenseCoeffsBase<Derived, 1>::Scalar = double; Eigen::Index = long int]: Assertion `index >= 0 && index < size()' failed.

Best regards,
Son Pham

Revision history for this message
Robert Caulk (rcaulk) said :
#3
Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#4

If you can't point out precisely what makes the script fail by comparison with a working script (as suggested by Robert in #3) it will be difficult to understand the problem.
Bruno

Revision history for this message
Son Pham Thai (pham-thai-son-987) said :
#5

Hi Bruno,

Sorry for the late reply, it took me a while to compare the scripts.

1. Regarding the triangulation step, my scripts use initialization() while the suggested script uses updateTriangulation=True

2. The appearance order of the bodies:
    - My scripts: the spheres are generated before the walls
    - The suggested script has the walls generated before the spheres

I have generated the spheres, then the walls, then I save these data. Then I reload the wall before the sphere, and initialization() works fine for me now.

Cheers,
Thai-Son

Revision history for this message
Son Pham Thai (pham-thai-son-987) said :
#6

Thank you Robert for the suggested script. My old working scripts also have the appearance order of the bodies like this script.

Cheers,
Thai Son