Segmentation fault with PotentialBlocks

Dear all,

I am trying to use PotentialBlocks to model particles with nonregular shapes. In the following script, I define 2 PotentialBlocks from a sphere. I know there is a yade function, but it is a first example. The idea is then to have a nonregular shape. One grain is fixed, and the other is falling to it. The grains are defined with an angular discretization ('n_theta') related to the number of planes considered.

It appears when this discretization is equal to 7, the code is working well. But when this discretization is equal to at least 8, a Segmentation fault occurs. Grains are correctly generated but when the contact occurs the simulation stops.

What do you think about it? Is it because of the computational cost (increasing the number of planes)?

Thanks in advance for your answer.
Regards
Alexandre Sac-Morane

Here, my mwe :

from yade import utils
from potential_utils import *
import math

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Material Parameters
m = O.materials.append(FrictMat(young=-1, poisson=-1, frictionAngle=atan(0.5), density=2000, label='frictMat'))

Kn = 1e7 #Pa/m
Ks = Kn * 2 / 3 #Pa/m
kn_i = 5 * Kn
ks_i = 5 * Ks

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([PotentialBlock2AABB()], verletDist=0.00),
        InteractionLoop(
                [Ig2_PB_PB_ScGeom(calContactArea=True)],
                [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=kn_i, ks_i=ks_i, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
                [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', initialOverlapDistance = 1e-6, allowViscousAttraction=False)]
        ),
        NewtonIntegrator(damping=0.0, exactAsphericalRot=True, gravity=[0, 0, -9.81], label='newton')
]

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Parameters
radius = 1 # the radius of grains
r = 0.1
n_theta = 7 #7 max with this definition on my laptop
n_phi = int(n_theta/2+1)

# creation of a discretized sphere
a_L = []
b_L = []
c_L = []
d_L = []
for i_phi in range(n_phi):
 phi = math.pi*i_phi/(n_phi-1)
 # top plane
 if i_phi == 0:
  a_L.append(0) # x coordinate of the norm of the plane
  b_L.append(0) # y coordinate of the norm of the plane
  c_L.append(1) # z coordinate of the norm of the plane
  d_L.append(radius-r) # distance to the center of the plane
 # bottom plane
 elif i_phi == n_phi-1:
  a_L.append( 0) # x coordinate of the norm of the plane
  b_L.append( 0) # y coordinate of the norm of the plane
  c_L.append(-1) # z coordinate of the norm of the plane
  d_L.append(radius-r) # distance to the center of the plane
 # current plane
 else :
  if int(n_theta*math.sin(phi)) != 0 :
   for i_theta in range(int(n_theta*math.sin(phi))):
    theta = 2*math.pi*i_theta/int(n_theta*math.sin(phi))
    a_L.append(math.cos(theta)*math.sin(phi)) # x coordinate of the norm of the plane
    b_L.append(math.sin(theta)*math.sin(phi)) # y coordinate of the norm of the plane
    c_L.append( math.cos(phi)) # z coordinate of the norm of the plane
    d_L.append(radius-r) # distance to the center of the plane

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# potential_utils.potentialblock: Generic Potential Block defined by a,b,c,d

# grain 1
g1 = Body()
g1.aspherical = True
g1.shape = PotentialBlock(
        a=a_L,
        b=b_L,
        c=c_L,
        d=d_L,
        r=r,
        R=0.0,
        AabbMinMax=True,
        isBoundary=True
)
utils._commonBodySetup(g1, g1.shape.volume, g1.shape.inertia, material='frictMat', pos=[0, 0, 0], fixed=True)
g1.state.pos = [0, 0 , -radius]
g1.state.ori = g1.shape.orientation
O.bodies.append(g1)

# grain 2
g2 = Body()
g2.aspherical = True
g2.shape = PotentialBlock(
        a=a_L,
        b=b_L,
        c=c_L,
        d=d_L,
        r=r,
        R=0.0,
        AabbMinMax=True,
        isBoundary=False
        )
utils._commonBodySetup(g2, g2.shape.volume, g2.shape.inertia, material='frictMat', pos=[0, 0, 0], fixed=False)
g2.state.pos = [0, 0 , 1.1*radius]
g2.state.ori = g2.shape.orientation
O.bodies.append(g2)

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Timestep
O.dt = 5e-5

O.run()