# Converting polyhedron body to Potential Particle or Potential Block

Asked by Ali Rafiee on 2019-10-24

Dear all,

I would like to know if there is a method for converting polyhedron objects to geometric format of potential particle or not.

Thanks a lot
Ali

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
2019-10-29
Last query:
2019-10-29
2019-10-29

## This question was reopened

 Vasileios Angelidakis (vsangelidakis) said on 2019-10-24: #1

Hi Ali,

To answer your question briefly, there is no built in function to make the transformation, but you can develop a function to do this. In more detail:

- Converting a polyhedron to a Potential Particle (PP) is not an exact transformation, since the edges of a PP are rounded and its faces are curved, so there is not an exact analogy between the shape of a PP and the polyhedral shape it is based upon. Though, you can approximate a polyhedron closely enough with a Potential Particle, ending up with a particle shape which is based on a polyhedron, and is a rounded version of it. If you would like to model a PP with fairly sharp edges (still rounded, but with a small radius "r") and with practically flat faces (very small value for k; close to zero), you can get something close enough to the original polyhedron.

- Converting a polyhedron to a Potential Block (PB) on the other hand is less complicated, since both codes (Polyhedra & PotentialBlock) simulate polyhedral particle shapes with sharp edges and flat faces.

To make the transformation from polyhedra to PP or PB:
After you have defined a polyhedron, you can use its attribute b.shape.getSurfaces()  to get the indices of the vertices lying on each face. You can calculate the normal vector for each face (coefficients: "a,b,c" of the Potential Blocks) by defining two vectors on each face and considering their cross product. Then, using a vertex on each face, you can then calculate the coefficients "d". The only thing left to define a Potential Block is to assume a radius "r", which is rather a modelling choice (i.e. there is no right/wrong value, as long as the chosen value is fairly small in regards to the particle dimensions). The selected value of "r" should then be deducted from each parameter d (i.e. d[i]=d[i]-r), to achieve the exact size of the initial polyhedral shape.

Hope this helps.
All the best,
Vasileios

 Ali Rafiee (ali-rafiee) said on 2019-10-24: #2

Thanks Vasileios Angelidakis, that solved my question.

 Ali Rafiee (ali-rafiee) said on 2019-10-28: #3

Thank you very much Vasileios,

I tried to do this converting method, and that is my code:

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

import numpy as np
import math
import random

######################################################################
def dotproduct(v1, v2):
return sum((a*b) for a, b in zip(v1, v2))

def length(v):
return math.sqrt(dotproduct(v, v))

def veccc(pp1,pp2):
pp3=[pp2-pp1,pp2-pp1,pp2-pp1]
if length(pp3)>0:
pp3=[pp3/length(pp3),pp3/length(pp3),pp3/length(pp3)]

return pp3
####################################
### calculate d value of a surface with normal vector and one point
def dvalue(vecn1,pp1):
dd1=-1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1

################################################################
powderDensity = 2000
distanceToCentre = 0.5
meanSize = 1.0
wallThickness = 0.5*meanSize
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless')) #The normal and shear stifness values are determined in the IPhys functor, thus the young,

########################################
young = 1e9
density=2500
FricDegree1= atan(.6)
m = PolyhedraMat()
m.density = density
m.young = young
m.poisson = 0.25
m.frictionAngle = FricDegree1

c1=polyhedron(((0.62500000,0.00000000,0.00000000),(0.62500000,0.50000000,0.00000000),(0.65382210,0.50000000,0.29263550),(0.65382210,0.00000000,0.29263550),(0.16342940,0.50000000,0.39018060),(0.16342940,0.00000000,0.39018060),(0.00000000,0.50000000,0.39018060),(0.00000000,0.00000000,0.39018060),(0.00000000,0.50000000,0.00000000),(0.00000000,0.00000000,0.00000000)),material=m)
O.bodies.append(c1)

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
#print(len(vs))
face2=[]
face2=b.shape.GetSurfaces()
#print(face2)
#print(vs)

id1=0
while id1<len(face2):
face11=face2[id1]
if len(face11)>2:
vec1=veccc(vs[face11],vs[face11])
vec2=veccc(vs[face11],vs[face11])

vects=[]
vects = list(np.cross(vec1, vec2))

#print(vects)
dvalue2=dvalue(vects,vs[face11])

#print("*******************************")
#print(dvalue2)
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
#print(aa)
#print(bb)
#print(cc)
#print(dd)

#### after that how can put correct values for PotentialBlock() ?????
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.001, r=chosenR, R=chosenR, a=aa, b=bb, c=cc, d=dd, id=len(O.bodies))
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos = [-2,0,0]
lidID = O.bodies.append(bbb)

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

the values for a, b, c and d are obtained, but I can't complete all parameters required by PotentialBlock() command just with these 4 parameters.

in any case, I do not know how to complete the last phase of this geometric transformation.

Thanks again
Ali

 Vasileios Angelidakis (vsangelidakis) said on 2019-10-28: #4

Hi Ali,

In general your code works. Visualisation in qt fails due to a bug, which for now you can avoid. The temporary fix:
- Use len(O.bodies)-1 for the potential block. There is a bug in the visualisation code of the Potential Blocks, which assumes the ids of the Potential Blocks start from 0. Will fix this soon

- For Vector3 types we have built-in functions to do standard calculations among vectors:
- dot product: v1.dot(v2)
- cross product: v1.cross(v2)
- norm (magnitude): v1.norm()
- subtraction: pp3=pp2-pp1 (you needn't do it separately for each item of the Vector3). Same for addition.
- normalize: pp3.normalize()
You can find these and many more if you type: "Vector3." + [tab] in the iPython interface! :)

Regarding the syntax of the PotentialBlocks:

> r=chosenR, R=chosenR
- Better use: r=chosenR, R=0.0. If R=0.0, an automatic value will be calculated, as the distance of the farthest vertices of the particle. You want "r" to be sufficiently small, while "R" to be around the dimensions of the particle.

> k=0.001
- The parameter "k" is not actively used in the PotentialBlocks, so you can set k=0.0 safely (or not even set a value for it at all); just don't spend any time thinking of what value to assign. ;)

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

- Check the sign of the "d[i]" values before subtracting "r". If d[i] is negative and you subtract "r", you are actually adding it. For some reason you seem to get the right particle shape here, but I think you should check twice if the particle shape you get is exactly correct (e.g. compare the vertices for the two particles, after assigning for them the same orientation).
- Do not use large "r" values: if dv<0: don't change the sign of dv; instead reduce "r". A change of sign here can lead to a different particle shape

All the best,
Vasileios

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

import numpy as np
import math
import random

######################################################################
#def dotproduct(v1, v2):
# return sum((a*b) for a, b in zip(v1, v2))

#def length(v):
# return math.sqrt(dotproduct(v, v))

#def veccc(pp1,pp2):
# pp3=pp2-pp1 #pp3=[pp2-pp1,pp2-pp1,pp2-pp1]
# if pp3.norm()>0: # if length(pp3)>0:
# pp3.normalize() # pp3=[pp3/length(pp3),pp3/length(pp3),pp3/length(pp3)]
# return pp3

####################################
### calculate d value of a surface with normal vector and one point
def dvalue(vecn1,pp1):
dd1=-1*(vecn1*pp1+vecn1*pp1+vecn1*pp1)
return dd1

################################################################
powderDensity = 2000
distanceToCentre = 0.5
meanSize = 1.0
wallThickness = 0.5*meanSize
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless')) #The normal and shear stifness values are determined in the IPhys functor, thus the young,

########################################
young = 1e9
density=2500
FricDegree1= atan(.6)
m = PolyhedraMat()
m.density = density
m.young = young
m.poisson = 0.25
m.frictionAngle = FricDegree1

c1=polyhedron((
[0.62500000,0.00000000,0.00000000],[0.62500000,0.50000000,0.00000000],[0.65382210,0.50000000,0.29263550],[0.65382210,0.00000000,0.29263550],[0.16342940,0.50000000,0.39018060],[0.16342940,0.00000000,0.39018060],[0.00000000,0.50000000,0.39018060],[0.00000000,0.00000000,0.39018060],[0.00000000,0.50000000,0.00000000],[0.00000000,0.00000000,0.00000000]),material=m)
O.bodies.append(c1)

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() #vec1=veccc(vs[face11],vs[face11])
vec2=vs[face11]-vs[face11]; vec2.normalize() #vec2=veccc(vs[face11],vs[face11])

# vects=[]
vects=vec1.cross(vec2) #vects = list(np.cross(vec1, 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

#### after that how can put correct values for PotentialBlock() ?????
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
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 = [-1,0,0]
# bbb.state.ori=Quaternion((0,1,0),pi)
lidID = O.bodies.append(bbb)

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

 Ali Rafiee (ali-rafiee) said on 2019-10-28: #5

Hi Vasileios, and thank you again.

just without last rotation around Y axis, the obtained PB is not same as the input polyhedron regarding to its orientation.

is it related to the sign of d value?

I want to convert large number of polyhedra, it is not possible to check them one by one, how can I solve this issue automatically?

Merci

Best regards
Ali

 Vasileios Angelidakis (vsangelidakis) said on 2019-10-28: #6

> just without last rotation around Y axis, the obtained PB is not same as the input polyhedron regarding to its orientation.

Regarding the orientation:
If you assign bbb.state.ori=bbb.shape.orientation, you should get the orientation of the particle in the random (non-principal) system you assigned it at, initially. I see that this doesn't work in this case, which makes me infer that yes, the sign of "d" plays a role here. If you plot the "dv" values before subtracting "chosenR", you will see that two of them are zero. Please have a look at the calculation of "dv" and find out why that happens (since all "dv" values should be non-zero and positive, for them to make sense).

Regarding the positioning of the particle:
Bear in mind that the PotentialBlock is centred to its centroid upon its definition, while I think you assign the particle to a local coordinate system, which is not centred around the centroid initially, so you can expect a discrepancy there.

> how can I solve this issue automatically?
I feel your struggle for automation :D I think after you sort out the issue with the zero "dv" values, you shouldn't need to do any manual checks for the bulk of particle shapes you use. You can do one or two checks for selected particles, to ensure that everything works fine. To automate the procedure, consider creating a function to do the transformation, instead of putting the whole script inside a loop for each particle shape.

Regards,
Vasileios

 Ali Rafiee (ali-rafiee) said on 2019-10-29: #7

Hi Vasileios,

I solve the problems related to the geometric transformation, but I think I have another problem concerning parameter setting in Potential block part, because I do not know exactly how it works.

1 ## below, you can find my code for an arch structure, it should slide with 0 degree for friction angle without cohesion, but it remains without any movement.

2 ## in order to generate vtk files with PotentialBlockVTKRecorder, the input parameters such as sampleX, how are they used?

Thanks a lot
Best regards
Ali

################################################# code part

import numpy as np
import math
import random
######################################
O.materials.append(FrictMat(young=1e10,poisson=0.25,frictionAngle=radians(0.0),density=2500,label='frictionless')) #The normal and shear stifness values are determined in the IPhys functor, thus the young,

chosenR=0.0001

O.engines=[
ForceResetter(),
InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.01),
InteractionLoop(
[Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, useFaceProperties=False, calJointLength=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,0,-9.81]),
#PotentialBlockVTKRecorder(fileName='./vtk/cubePBscaled',iterPeriod=10000,twoDimension=False,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2,label='vtkRecorder')
]

################################## geometry part, by A.Rafiee
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[-0.0, -1.0, -0.0, -0.0, 1.0, -0.0], b=[-1.0, -0.0, -0.0, -0.0, -0.0, 1.0], c=[0.0, -0.0, -1.0, 1.0, -0.0, -0.0], d=[1.8749, 2.8124, 0.0624, 0.0624, 2.8124, 1.8749], id=0)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos =[2.125, 0.25, -0.0625]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[-0.995184720391626, 0.0, -0.0, -0.0, -0.0, -0.19509018176012233, 1.0], b=[2.7205239199536258e-18, -1.0, -5.551115123125783e-17, 1.0, -0.0, 1.0829680583419728e-17, -1.1102230246251564e-16], c=[0.09801720409724536, -0.0, 1.0, -0.0, -1.0, -0.9807853083018742, -2.8454080613571164e-16], d=[0.337181872405393, 0.2499, 0.17967108333414306, 0.2499, 0.21030951666585695, 0.17888318179899965, 0.303692080616296], id=1)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[0.30379208061629603, 0.25, 0.17977108333414304]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.1950901817601223, 0.9569403497890105, -1.4555161670882268e-16, -0.38268355158955414, 2.262943718606007e-16, -0.9569403375852945], b=[-4.537043464545097e-17, 1.3428361705711876e-17, -0.9951847174017722, -1.466534115884631e-17, 0.9951847174017722, 2.0392625979676485e-16], c=[0.9807853083018742, -0.2902846309153932, 1.1777448297131143e-16, -0.9238794831268877, 4.976255189356365e-18, 0.290284671145658], d=[0.1717652880563909, 0.23684874980903367, 0.14917770761026583, 0.17176528561465593, 0.14917770761026589, 0.26054366159508413], id=2)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[0.44708096679737336, 0.25000000000000006, 0.5089911606262629]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.38268355158955414, 0.8819212004273606, 8.127355190878246e-17, -0.5555703541431359, 1.958258197514367e-16, -0.8819212357727806], b=[-1.2821403426828246e-16, 6.54194554658077e-17, -0.9951847260300383, -1.7972755132828346e-16, 0.9951847260300383, -1.1957524611753766e-16], c=[0.9238794831268877, -0.47139685641374757, 1.9621163966280541e-16, -0.8314695313703745, -1.0467111550412906e-16, 0.4713967902871334], d=[0.17176529181696487, 0.2368487673900063, 0.1491777089045058, 0.17176529360968104, 0.14917770890450577, 0.2605436475477011], id=3)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[0.5786209907389809, 0.25, 0.8265567947287535]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.5555703541431359, 0.7730105349646802, -3.5803953569760556e-18, -0.7071067811865477, 8.582140233584773e-17, -0.7730106376137247], b=[3.8463192416736605e-17, -2.934658001475229e-17, -0.9951847535586482, -1.3084103822364795e-16, 0.995184753558648, -1.4303570135982168e-16], c=[0.8314695313703745, -0.634393184731377, 9.407783687275172e-17, -0.7071067811865475, -7.043178622195482e-17, 0.6343930596531008], d=[0.17176528321389758, 0.23684876708531535, 0.1491777130337973, 0.1717652803907167, 0.1491777130337972, 0.2605436269756058], id=4)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[0.7695874613517345, 0.25, 1.1123583282186593]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.7071067811865478, 0.6343931847313773, 5.936831281291228e-17, -0.8314695313703748, -1.298001029743927e-16, -0.6343930596531011], b=[-5.233641528945917e-16, 8.325647318886739e-16, -0.9951847140642529, 3.08054519409162e-16, 0.995184714064253, -9.716725512485364e-17], c=[0.7071067811865474, -0.7730105349646799, 3.718747074082228e-16, -0.5555703541431357, -2.8605330027204803e-16, 0.7730106376137245], d=[0.1717652803907161, 0.23684876708531505, 0.14917770710963785, 0.17176528321389797, 0.1491777071096379, 0.26054362697560574], id=5)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[1.0126416717813402, 0.25, 1.3554125386482654]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.8314695313703745, 0.47139685641374734, -1.6347077095121204e-16, -0.9238794831268876, -6.588657534994284e-18, -0.47139679028713377], b=[2.6049271698971687e-17, -3.7280529721843463e-16, -0.9951847256846824, -9.609762509852285e-16, 0.9951847256846824, -1.9134732646077375e-16], c=[0.5555703541431363, -0.8819212004273607, -4.4143378261607424e-16, -0.38268355158955436, 4.606860979300525e-16, 0.8819212357727804], d=[0.17176529360968132, 0.2368487673900062, 0.1491777088527024, 0.17176529181696482, 0.1491777088527025, 0.2605436475477015], id=6)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[1.2984432052712465, 0.25000000000000006, 1.5463790092610195]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.9238794831268876, 0.2902846309153935, 1.6303580705062413e-16, -0.9807853083018743, -4.9289839508535835e-17, -0.29028467114565765], b=[-3.1257340906439395e-16, 3.0042562302002276e-16, -0.9951847092768998, 2.356301414989026e-16, 0.9951847092768997, -2.118908447279992e-16], c=[0.382683551589554, -0.9569403497890105, -5.2059346508324875e-17, -0.19509018176012224, -8.021203425885159e-17, 0.9569403375852945], d=[0.1717652856146558, 0.2368487498090334, 0.149177706391535, 0.17176528805639094, 0.14917770639153488, 0.26054366159508413], id=7)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[1.6160088393737368, 0.24999999999999997, 1.6779190332026268]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.9807853083018742, 0.09801705227898266, 6.399136854325753e-17, -1.0, -8.031448525403647e-17, -0.09801720409724561], b=[-1.8049467639032894e-17, 9.20730839112388e-17, -0.9951847320167924, -1.8503717077085946e-16, 0.9951847320167923, -1.8136826133024222e-17], c=[0.1950901817601225, -0.995184735344418, 6.90547993979119e-17, 3.3306690738754706e-16, 9.667673614731509e-17, 0.9951847203916262], d=[0.17176527403820857, 0.23684876474180003, 0.14917770980251885, 0.17176527838164274, 0.14917770980251882, 0.2605436253595993], id=8)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[1.9531347216183572, 0.25, 1.7449775064827124]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[1.0, -0.09801705227898319, -6.153106884084172e-16, -0.9807853083018743, 4.791663943962904e-16, 0.09801720409724606], b=[5.551115123125779e-16, -6.068488975904394e-16, -0.995184735344418, -6.537290837710836e-16, 0.995184735344418, 4.408396346434073e-16], c=[5.551115123125779e-16, -0.995184735344418, -4.971946531206894e-16, 0.19509018176012194, 5.524385034674323e-16, 0.995184720391626], d=[0.17176527838164363, 0.23684876474179992, 0.14917771030166288, 0.17176527403820796, 0.14917771030166294, 0.2605436253595999], id=9)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[2.2968652783816434, 0.25000000000000006, 1.7449775064827127]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.9807853083018743, -0.2902846309153928, -5.3318352810194274e-17, -0.9238794831268876, -4.976255189356427e-18, 0.2902846711456571], b=[1.0928386802567227e-16, -4.078525166657175e-16, -0.9951847174017724, -0.0, 0.9951847174017725, 0.0], c=[-0.19509018176012255, -0.9569403497890108, 6.693180533084786e-17, 0.3826835515895541, 2.262943718606008e-16, 0.9569403375852947], d=[0.171765288056391, 0.23684874980903325, 0.14917770761026583, 0.1717652856146561, 0.14917770761026597, 0.26054366159508446], id=10)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[2.633991160626263, 0.25000000000000006, 1.677919033202627]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.9238794831268876, -0.47139685641374807, -2.6167778876032296e-17, -0.8314695313703748, -3.9236101011569525e-17, 0.47139679028713394], b=[-0.0, 0.0, -0.9951847260300383, 2.5874781661336978e-17, 0.9951847260300379, 2.06801174502733e-16], c=[-0.3826835515895542, -0.8819212004273601, -4.895645493785917e-17, 0.5555703541431356, 7.604763890745331e-17, 0.8819212357727804], d=[0.17176529181696437, 0.2368487673900065, 0.14917770890450574, 0.1717652936096813, 0.14917770890450574, 0.26054364754770065], id=11)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[2.951556794728753, 0.25000000000000006, 1.546379009261019]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.8314695313703749, -0.6343931847313768, -9.407783687275171e-17, -0.7071067811865477, 4.6785735571157934e-17, 0.6343930596531009], b=[-2.5874781661336954e-17, -2.017288423954374e-16, -0.9951847535586482, -3.271025955591196e-17, 0.9951847535586482, 2.8976857250342004e-16], c=[-0.5555703541431356, -0.7730105349646803, -3.580395356976111e-18, 0.7071067811865476, 1.6806240931471946e-16, 0.7730106376137246], d=[0.17176528321389806, 0.23684876708531497, 0.14917771303379723, 0.17176528039071642, 0.14917771303379732, 0.26054362697560607], id=12)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[3.2373583282186598, 0.25000000000000006, 1.3554125386482656]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.7071067811865476, -0.7730105349646805, 1.4294539431761725e-17, -0.5555703541431365, -1.0011594656793655e-16, 0.7730106376137243], b=[-3.271025955591198e-17, 3.575892059148325e-17, -0.9951847140642529, -1.281524780432438e-16, 0.9951847140642531, -1.5314522225483436e-16], c=[-0.7071067811865477, -0.6343931847313765, -1.7054773672941695e-16, 0.8314695313703743, 1.0011594656793657e-16, 0.6343930596531013], d=[0.17176528039071673, 0.2368487670853149, 0.1491777071096379, 0.1717652832138976, 0.14917770710963804, 0.2605436269756061], id=13)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[3.480412538648266, 0.25000000000000006, 1.1123583282186595]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.5555703541431364, -0.8819212004273601, -5.940828355224751e-17, -0.3826835515895541, -9.791291418020404e-17, 0.8819212357727808], b=[-0.0, 0.0, -0.9951847256846824, -7.081068168667528e-17, 0.9951847256846825, 1.6318821031013536e-16], c=[-0.8314695313703744, -0.4713968564137484, 2.877832969195519e-16, 0.9238794831268876, -5.233556005286094e-17, 0.4713967902871328], d=[0.1717652936096816, 0.23684876739000565, 0.14917770885270232, 0.1717652918169646, 0.14917770885270248, 0.26054364754770154], id=14)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[3.67137900926102, 0.25000000000000006, 0.8265567947287534]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.3826835515895541, -0.9569403497890109, 1.3120987351899408e-16, -0.1950901817601222, -5.205934650832478e-17, 0.9569403375852948], b=[2.063573865343809e-16, -3.482132063787019e-17, -0.9951847092768996, -1.814817385818039e-16, 0.9951847092768998, -5.371345426694166e-17], c=[-0.9238794831268877, -0.2902846309153925, 2.6066774258432397e-16, 0.9807853083018743, -1.6303580705062395e-16, 0.2902846711456568], d=[0.17176528561465604, 0.23684874980903342, 0.1491777063915349, 0.17176528805639077, 0.14917770639153502, 0.2605436615950846], id=15)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[3.802919033202627, 0.25, 0.5089911606262629]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#####################
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=chosenR, R=0.0, a=[0.9951847203916262, 0.19509018176012233, 0.0, -1.0, -0.0, -0.0, -8.881784197001253e-17], b=[-0.0, -0.0, -0.19509018176012233, -0.0, 0.19509018176012233, -0.0, -5.551115123125783e-17], c=[0.09801720409724574, -0.9807853083018744, -0.0, -0.0, 0.0, -1.0, 1.0], d=[0.33718187240539327, 0.17888318179899973, 0.04867254544003058, 0.30369208061629593, 0.04867254544003057, 0.21030951666585695, 0.17967108333414306], id=16)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
bbb.state.pos =[3.946207919383704, 0.25, 0.17977108333414304]
bbb.state.ori=Quaternion((0,1,0),pi)
O.bodies.append(bbb)
#########################################################################################################################

 Ali Rafiee (ali-rafiee) said on 2019-10-29: #8

with setting dt=0.0001, my code works but I see this error message:
"analytic centre did not converge"

thanks
Ali Vasileios Angelidakis (vsangelidakis) said on 2019-10-29: #9

Hi Ali,

Nice arch :)
This thread is getting too lengthy, so I would advise to open a new question if you need to ask something else.

> the input parameters such as sampleX, how are they used?

See [1,2].
The PotentialBlockVTKRecorder does not record the actual polyhedral particle shape. What it records is a rounded version of the particle, close to the shape of the inner PP, which is used to calculate the contact normal. To do that, the PotentialBlockVTKRecorder considers a cubic drawing space, deriving its dimensions from the attributes: PotentialBlock.minAabb, PotentialBlock.maxAabb (it considers the maximum absolute x,y,z values). This cubic space is discretized into a grid, using the sampleX, sampleY, sampleZ parameters (i.e. the size of the grid.). Large sampleX,Y,Z values will give you better detail, but will make the output larger and slower to write.

So, what you can do is assign the values below, right after defining each particle shape:
bbb.shape.minAabb=Vector3(bbb.shape.R, bbb.shape.R, bbb.shape.R)
bbb.shape.maxAabb=Vector3(bbb.shape.R, bbb.shape.R, bbb.shape.R)

Having in mind how these values are used now, you can increase/decrease them per direction until you get a decent visualisation.
Also make sure to create a "vtk" folder, as I see you don't create one at the moment in your script.

In general, the PotentialBlockVTKRecorder can be slow if very detailed visualisation is requested, but it is very helpful to assess visually how rounded the inner Potential Particle is (the one used to calculate the contact normals) and if the chosen "r" radius is appropriate in relation to the dimensions of each particle.

* A different question: Have you reduced the r values from the d arrays, i.e. d[i]=d[i]-r? If not, your particles will be oversized by r, in every dimension. In case you have done this manually, I would advise to do this in the working script, so in case you decide to change "r", the "d" values are updated too.

> "analytic centre did not converge"

This can happen for a number of reasons. The most probable:
1. The chosen radius "r" is too small
2. The particles are exactly touching, without practically any initial overlap (which is the case here).
This warning is not the end of the world, since you will see it happens for a few timesteps, when the particles do not have enough overlap for the analytic centre (i.e. the contact point) to be calculated, but then the code starts to converge.

To solve 2, there is a trick/workaround you can do: You can increase the size of the particles by a very very small value, which will then be deducted from the penetrationDepth in the contactlaw , to achieve the actual contact conditions.

To do this, set these variables:
initialOverlap=1e-8 # a very small distance in comparison with the dimensions of the problem
chosenR=0.0001+initialOverlap # particles are enlarged by initialOVerlap

and then inside the contact law, use:
initialOverlapDistance=initialOverlap # the contact law assumes a negative (tensile) gap equal to initialOverlap

These two practically cancel each other out, while contact is established from a geometrical perspective.

Please be more specific if you don't want me to start explaining the code line by line :) A good read of the documentation and the source code will help on this.

Best of luck,
Vasileios

 Ali Rafiee (ali-rafiee) said on 2019-10-29: #10

Thank you very much Vasileios for your time.

I think you're right, for the moment, I have a lot to do with your detailed explanations.
I'm finishing here for this question.

Best Regards
Ali

 Ali Rafiee (ali-rafiee) said on 2019-10-29: #11

Thanks Vasileios Angelidakis, that solved my question.