New question #707162 on Yade:
https://answers.launchpad.net/yade/+question/707162
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()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
