# How to generate many polyhedra randomly without specific data in Potential Blocks code？

Hello,all
How to generate many polyhedra randomly without specific data in Potential Blocks code？

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-05-26
2020-05-26

## This question was reopened

 Vasileios Angelidakis (vsangelidakis) said on 2020-03-22: #1

Hi weijie,

Currently, I have not written a dedicated function to do something like this for the Potential Blocks. What you need to define a particle using this code, are mainly the equations of the planes corresponding to the particle faces, i.e. the coefficients "a,b,c,d". And then, you can choose an appropriate radius "r", which will be used in the calculation of the contact normals when a contact occurs.

To calculate the plane coefficients a,b,c,d, you need to have a mesh of the particle surface (vertices and connectivity), which this code currently does not give you. You can though use a different approach to calculate such a mesh, and then transform it into a Potential Block. E.g, the polyhedra code has a very nice function to generate random polyhedra with controlled size  and it also gives you the connectivity of the vertices making the particle surface .

All the best, Vasileios

 weijie (amandajoe) said on 2020-03-26: #2

Hi Vasileios,
I tried to make a random polyhedron with reference to the information you provided. However, I encountered a difficulty. A part of the converted polyhedron is displayed in black and it does not feel closed.How can I solve it？
Following is my code：
#####################
import numpy as np
import math
import random
#########################################################################Define physical parameters
powderDensity = 2500
young = 1e9
FricDegree1= atan(.6)
m = PolyhedraMat()
m.density = powderDensity
m.young = young
m.poisson = 0.25
m.frictionAngle = FricDegree1
#########################################################################Forming polyhedra with polyhedra function
t = polyhedra_utils.polyhedra(m,size = (0.06,0.06,0.06),seed =5)
#t.state.pos = (0.,0.,0.5)
O.bodies.append(t)
########################################################################Polyhedron transformation
def dvalue(vecn1,pp1):
dd1=-1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1
chosenR=0.0001
for b in O.bodies:
aa=[]
bb=[]
cc=[]
dd=[]
if not isinstance(b.shape,Polyhedra): # skip non-polyhedra bodies
continue
vs = [b.state.pos + b.state.ori*v for v in b.shape.v] # vertices in global coords
face2=b.shape.GetSurfaces()

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=vs[face11]-vs[face11]; vec1.normalize()
vec2=vs[face11]-vs[face11]; vec2.normalize() #Normalize this object in-place.
vects=vec1.cross(vec2)

dvalue2=dvalue(vects,vs[face11]) # Some dvalue2 values return equal to 0. Check this part of your script once more.
dv=dvalue2-chosenR

if dv<0:
vects=[-1*vects,-1*vects,-1*vects]
dv=-1*dv

aa.append(vects)
bb.append(vects)
cc.append(vects)
dd.append(dv)

id1=id1+1
#######################################################################################
O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, useFaceProperties=False, viscousDamping=0.2)],
),
#GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
#PotentialBlockVTKRecorder(fileName='./vtk/cubePBscaled',iterPeriod=10000,twoDimension=False,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2,label='vtkRecorder')
]
#print(aa,bb,cc,dd)
bbb=Body()
wire=False
color=[125,2,1]
highlight=True
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=aa, b=bb, c=cc, d=dd, id=len(O.bodies)-1)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos = [0,0.1,0]
# bbb.state.ori=Quaternion((0,1,0),pi)
lidID = O.bodies.append(bbb)

 Vasileios Angelidakis (vsangelidakis) said on 2020-03-26: #3

Hi weijie,

The black part is only a visualisation issue. We cull (hide) the back faces of the particles and for some reason one of the faces is oriented inside out. I will have a look on how to resolve it soon, but it is nothing to worry about. The particle is defined correctly.

I would advise using a larger chosenR value for numerical stability reasons (but still chosenR<min(dd)). This also happens to be solving your visualisation issue here.

Also, in your script, you haven't defined a contact law.

Regards,
Vasileios

 weijie (amandajoe) said on 2020-03-27: #4

Hi Vasileios, and thank you again.
I found that the direction of the transformed polyhedron is different from the original polyhedron. What if I want to make the two polyhedrons have the same direction? At the same time, what should I do to remove the previous polyhedron?
Best regards
weijie

 Vasileios Angelidakis (vsangelidakis) said on 2020-03-27: #5

> I found that the direction of the transformed polyhedron is different from the original polyhedron.

Yes, this has been reported before and is probably a bug in the calculation of the principal orientations of the PotentialBlocks. In general, when a particle has b.state.ori=Quaternion((1,0,0),0) it should be oriented to its principal axes and the two codes should lead to the same result. Will have another look if I can spot this bug.

> At the same time, what should I do to remove the previous polyhedron?

This is easier :) Do not append the body in the body container (i.e. do not use O.bodies.append() for the polyhedron). Just define its shape to calculate its surface triangulation.

Best regards,
Vasileios

 weijie (amandajoe) said on 2020-03-30: #6

Thanks Vasileios Angelidakis, that solved my question.

 Vasileios Angelidakis (vsangelidakis) said on 2020-04-22: #7

Hi again,

> I found that the direction of the transformed polyhedron is different from the original polyhedron.

This is now fixed in the latest trunk, where you can achieve the same orientation between PBs and polyhedra.

In your script above, comment these lines:

# if dv<0:
# vects=[-1*vects,-1*vects,-1*vects]
# dv=-1*dv

Best regards,
Vasileios

 weijie (amandajoe) said on 2020-04-23: #8

Hi Vasileios, and thank you again.

I have updated the latest trunk and also commented out these three lines, but I found that the orientation of the resulting polyhedron is still different. Can you please help to see what went wrong?

Following is my code：

#####################
import numpy as np
import math
import random
#########################################################################Define physical parameters
powderDensity = 2500
young = 1e9
FricDegree1= atan(.6)
m = PolyhedraMat()
m.density = powderDensity
m.young = young
m.poisson = 0.25
m.frictionAngle = FricDegree1
#########################################################################Forming polyhedra with polyhedra function
t = polyhedra_utils.polyhedra(m,size = (0.06,0.06,0.06),seed =5)
t.state.ori=Quaternion((1,0,0),0)
#t.state.pos = (0.,0.,0.5)
O.bodies.append(t)
########################################################################Polyhedron transformation
def dvalue(vecn1,pp1):
dd1=-1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1
chosenR=0.001
for b in O.bodies:
aa=[]
bb=[]
cc=[]
dd=[]
if not isinstance(b.shape,Polyhedra): # skip non-polyhedra bodies
continue
vs = [b.state.pos + b.state.ori*v for v in b.shape.v] # vertices in global coords
face2=b.shape.GetSurfaces()

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=vs[face11]-vs[face11]; vec1.normalize()
vec2=vs[face11]-vs[face11]; vec2.normalize() #Normalize this object in-place.
vects=vec1.cross(vec2)

dvalue2=dvalue(vects,vs[face11]) # Some dvalue2 values return equal to 0. Check this part of your script once more.
dv=dvalue2-chosenR

#if dv<0:
#vects=[-1*vects,-1*vects,-1*vects]
#dv=-1*dv

aa.append(vects)
bb.append(vects)
cc.append(vects)
dd.append(dv)

id1=id1+1
#######################################################################################
O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, useFaceProperties=False, viscousDamping=0.2)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False, allowViscousAttraction=True, traceEnergy=False)]
),
#GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
#PotentialBlockVTKRecorder(fileName='./vtk/cubePBscaled',iterPeriod=10000,twoDimension=False,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2,label='vtkRecorder')
]
#print(aa,bb,cc,dd)
bbb=Body()
wire=False
color=[125,2,1]
highlight=True
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=aa, b=bb, c=cc, d=dd, id=len(O.bodies)-1)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos = [0,0.1,0]
bbb.state.ori=Quaternion((1,0,0),0)
O.bodies.append(bbb)
###################

Best regards,
Jie

 Vasileios Angelidakis (vsangelidakis) said on 2020-04-23: #9

Hi Jie,

You have actually revealed a different bug in the code; many thanks for that! I see that if size=(0.06,0.06,0.06), we do not get correct volume/centroid, due to undetected duplicate vertices in the PBs. I will go through the source code again to find a permanent fix for this. In the meanwhile, if you scale the particle e.g. with size=(0.6,0.6,0.6), like below, you get the correct volume/centroid/inertia/orientation/vertices.

All the best,
Vasileios

#####################
import numpy as np
import math
import random
#########################################################################Define physical parameters
powderDensity = 2500
young = 1e9
FricDegree1= atan(.6)
m = PolyhedraMat()
m.density = powderDensity
m.young = young
m.poisson = 0.25
m.frictionAngle = FricDegree1
#########################################################################Forming polyhedra with polyhedra function
t = polyhedra_utils.polyhedra(m,size=(0.6,0.6,0.6),seed =3) #size=(0.06,0.06,0.06)
#t.state.ori=Quaternion((1,0,0),0)
#t.state.pos = (0.,0.,0.5)
O.bodies.append(t)
########################################################################Polyhedron transformation
def dvalue(vecn1,pp1):
dd1=1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1

chosenR=0.001
for b in O.bodies:
aa=[]
bb=[]
cc=[]
dd=[]
if not isinstance(b.shape,Polyhedra): # skip non-polyhedra bodies
continue
vs = [b.state.ori*v for v in b.shape.v] # vertices in global coords
face2=b.shape.GetSurfaces()

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=vs[face11]-vs[face11]; vec1.normalize()
vec2=vs[face11]-vs[face11]; vec2.normalize() #Normalize this object in-place.
vects=vec1.cross(vec2); vects.normalize()

dvalue2=dvalue(vects,vs[face11]) # Some dvalue2 values return equal to 0. Check this part of your script once more.
dv=dvalue2-chosenR

# if dv<0:
# vects=[-1*vects,-1*vects,-1*vects]
# dv=-1*dv

aa.append(vects)
bb.append(vects)
cc.append(vects)
dd.append(dv)

id1=id1+1

#######################################################################################
O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, useFaceProperties=False, viscousDamping=0.2)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False, allowViscousAttraction=True, traceEnergy=False)]
),
#GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
#PotentialBlockVTKRecorder(fileName='./vtk/cubePBscaled',iterPeriod=10000,twoDimension=False,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2,label='vtkRecorder')
]
#print(aa,bb,cc,dd)
bbb=Body()
bbb.aspherical=True
wire=False
color=[125,2,1]
highlight=True
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=aa, b=bb, c=cc, d=dd)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=bbb.shape.position, fixed=False)
bbb.state.ori=bbb.shape.orientation
#bbb.state.pos = [0,0.1,0]
#bbb.shape.position+t.shape.GetCentroid()
#bbb.state.ori=Quaternion((1,0,0),0)
O.bodies.append(bbb)
###################

v=qt.View()

print("Volume")
print(t.shape.GetVolume())
print(bbb.shape.volume)

print("Centroid")
print(t.shape.GetCentroid())
print(bbb.shape.position)

print("Principal inertia tensor")
print(t.shape.GetInertia())
print(bbb.shape.inertia)

print("Principal orientations")
print(t.shape.GetOri())
print(bbb.shape.orientation)

 weijie (amandajoe) said on 2020-04-24: #10

Hi Vasileios,

Thanks a lot for your help, scale the particle makes the direction of the polyhedron consistent.

Best regards,
Jie

 weijie (amandajoe) said on 2020-05-21: #11

Hi Vasileios, and thank you again.

I have a new question now. What if I want to efficiently convert many different polyhedra generated in makeCloud?

Best regards,
Jie

Below is my MWS：
#################################
import random
import numpy as np

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e7,poisson=.2,density=2.5e3)
powderDensity = 2500
#-------------------------------------------
#Dimensions

meanSize = 0.05
wallThickness = 0.5*meanSize
distanceToCentre = 0.01
lengthOfBase = 0.250
heightOfBase = 0.600

#-------------------------------------------
#Make Cloud

sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))
R=sqrt(3.0)*distanceToCentre
sp.makeCloud(mn,mx,R,0,-1,False)

for s in sp:
b=Body()
color=Vector3(random.random(),random.random(),random.random())
b=polyhedra_utils.polyhedra(material=n,size=(2*distanceToCentre,distanceToCentre,distanceToCentre),color=color)
b.state.pos = s #s stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s
O.bodies.append(b)

def dvalue(vecn1,pp1):
dd1=1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1

chosenR=0.005
for b in O.bodies:
aa=[]
bb=[]
cc=[]
dd=[]
if not isinstance(b.shape,Polyhedra): # skip non-polyhedra bodies
continue
vs = [b.state.ori*v for v in b.shape.v] # vertices in global coords
face2=b.shape.GetSurfaces()

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=vs[face11]-vs[face11]; vec1.normalize()
vec2=vs[face11]-vs[face11]; vec2.normalize() #Normalize this object in-place.
vects=vec1.cross(vec2); vects.normalize()
dvalue2=dvalue(vects,vs[face11]) # Some dvalue2 values return equal to 0. Check this part of your script once more.
dv=dvalue2-chosenR
aa.append(vects)
bb.append(vects)
cc.append(vects)
dd.append(dv)
id1=id1+1

bbb=Body()
bbb.aspherical=True
wire=False
color=[125,2,1]
highlight=True
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=aa, b=bb, c=cc, d=dd)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=bbb.shape.position, fixed=False)
bbb.state.ori=bbb.shape.orientation
bbb.state.pos = [0,0.5,0]
O.bodies.append(bbb)

v=qt.View()
O.saveTmp()

 weijie (amandajoe) said on 2020-05-22: #12

In the above MWS, I found that only one polyhedron can be converted. How can I convert multiple polyhedrons?

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-25: #13

Hi,

To convert all the polyhedra into PBs, you have to define the generation of PBs inside the loop. You are currently creating only the last polyhedron. Then, for the newly created PBs, you have to assign the same orientation and position as the one of the polyhedral particles.

Below I generate the PBs shifted by Vector3(lengthOfBase,0,0), so that you can compare the polyhedra (on the left) and the PBs (on the right). Also, I choose different "r=chosenR" for each body, depending on their individual dd values.

I see some faces are still not rendered properly in the PBs, but this is only a visualisation issue.

#Best Regards,
#Vasileios

#################################
import random
import numpy as np

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e7,poisson=.2,density=2.5e3)
#-------------------------------------------
#Dimensions

meanSize = 0.05
wallThickness = 0.5*meanSize
distanceToCentre = 0.01
lengthOfBase = 5*meanSize
heightOfBase = 12*meanSize

#-------------------------------------------
#Make Cloud

sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))
R=sqrt(3.0)*distanceToCentre
sp.makeCloud(mn,mx,R,0,-1,False)

def dvalue(vecn1,pp1):
dd1=1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1

for s in sp:
# Generate polyhedra
color=Vector3(random.random(),random.random(),random.random())
b=polyhedra_utils.polyhedra(material=n,size=(2*distanceToCentre,distanceToCentre,distanceToCentre),color=color)
b.state.pos = s #s stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s
O.bodies.append(b)

# Generate PBs
aa=[]
bb=[]
cc=[]
dd=[]
vs=b.shape.v
face2=b.shape.GetSurfaces()

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=vs[face11]-vs[face11]; vec1.normalize()
vec2=vs[face11]-vs[face11]; vec2.normalize() #Normalize this object in-place.
vects=vec1.cross(vec2); vects.normalize()
dvalue2=dvalue(vects,vs[face11])
aa.append(vects)
bb.append(vects)
cc.append(vects)
dd.append(dvalue2)
id1=id1+1

chosenR=min(dd)/2

bbb=Body()
bbb.aspherical=True
wire=False
highlight=True
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=aa, b=bb, c=cc, d=np.array(dd)-chosenR, color=b.shape.color)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=bbb.shape.position, fixed=False)
bbb.state.ori= b.state.ori
bbb.state.pos = b.state.pos+Vector3(lengthOfBase,0,0)
O.bodies.append(bbb)

# Count number of bodies with b.shape=Polyhedra
countPol=0
for b in O.bodies:
if isinstance(b.shape,Polyhedra):
countPol=countPol+1
print("number of Polyhedra = ", countPol)

# Count number of bodies with b.shape=PotentialBlock
countPBs=0
for b in O.bodies:
if isinstance(b.shape,PotentialBlock):
countPBs=countPBs+1
print("number of PotentialBlocks = ", countPBs)

v=qt.View()
v.ortho=True # I activate orthotropic projection, to make visual comparisons easier
O.saveTmp()

 weijie (amandajoe) said on 2020-05-26: #14

Hi Vasileios, and thank you again.

I zoomed in and found that some of the converted particles have some overlapping parts, which seems to be a problem with the conversion. I changed the size of the particles generated in makeCloud and found that the smaller the particles, the more serious the problem. Could you please help to see what went wrong?

Best regards,
Jie

##################

import random
import numpy as np

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e7,poisson=.2,density=2.5e3)
#-------------------------------------------
#Dimensions

meanSize = 0.05
wallThickness = 0.5*meanSize
distanceToCentre = 0.01
lengthOfBase = 5*meanSize
heightOfBase = 12*meanSize

#-------------------------------------------
#Make Cloud

sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))
#R=sqrt(3.0)*distanceToCentre
#sp.makeCloud(mn,mx,R,0,-1,False)
sp.makeCloud(mn,mx,psdSizes=[distanceToCentre,2*distanceToCentre],psdCumm=(0.1,1),num=200)

def dvalue(vecn1,pp1):
dd1=1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1

# Generate polyhedra
color=Vector3(random.random(),random.random(),random.random())
b.state.pos = center #s stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s
O.bodies.append(b)

# Generate PBs
aa=[]
bb=[]
cc=[]
dd=[]
vs=b.shape.v
face2=b.shape.GetSurfaces()

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=vs[face11]-vs[face11]; vec1.normalize()
vec2=vs[face11]-vs[face11]; vec2.normalize() #Normalize this object in-place.
vects=vec1.cross(vec2); vects.normalize()
dvalue2=dvalue(vects,vs[face11])
aa.append(vects)
bb.append(vects)
cc.append(vects)
dd.append(dvalue2)
id1=id1+1

chosenR=min(dd)/2

bbb=Body()
bbb.aspherical=True
wire=False
highlight=True
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=aa, b=bb, c=cc, d=np.array(dd)-chosenR, color=b.shape.color)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=bbb.shape.position, fixed=False)
bbb.state.ori= b.state.ori
bbb.state.pos = b.state.pos+Vector3(lengthOfBase,0,0)
O.bodies.append(bbb)

# Count number of bodies with b.shape=Polyhedra
countPol=0
for b in O.bodies:
if isinstance(b.shape,Polyhedra):
countPol=countPol+1
print("number of Polyhedra = ", countPol)

# Count number of bodies with b.shape=PotentialBlock
countPBs=0
for b in O.bodies:
if isinstance(b.shape,PotentialBlock):
countPBs=countPBs+1
print("number of PotentialBlocks = ", countPBs)

v=qt.View()
v.ortho=True # I activate orthotropic projection, to make visual comparisons easier
O.saveTmp()

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-26: #15

Hi Jie,

> found that the smaller the particles, the more serious the problem
This is exactly the problem I have identified so far.

What we do here, is that from polyhedra we calculate the a,b,c,d of the PBs and using these, we recalculate the vertices, which are used to calculate the centroid/volume/inertia and also visualise the particles. I think the root of the problem is isolated in the calculation of vertices of the PBs, where some duplicate vertices are not detected properly. This is why the PBs we generate for these small particle sizes have more vertices than the initial polyhedra.

I will try to have a look in the coming days on how to solve this in a more robust way. Until then, you can:
- File a bug report , describing the issue.
- Consider using a different unit system (e.g. instead of meters/Newtons/Pascals use milimiters/Newtons/MegaPascals), to virtually upscale the size of the particles.
- Play with the tolerance of the duplicate vertices in the PB shape class , recompile, and see if you can achieve the same number of vertices between the generated PBs and the actual polyhedra.

These are messy solutions, but might help you progress with your work, until I find a permanent way to solve this.

Kind regards,
Vasileios

 weijie (amandajoe) said on 2020-05-27: #16

Hi Vasileios, and thank you again.

I will try to do it.

Best regards,
Jie