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

Asked by weijie on 2020-03-22

Hello,all
     How to generate many polyhedra randomly without specific data in Potential Blocks code?
     Thanks in advance.

Question information

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

This question was reopened

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 [1] and it also gives you the connectivity of the vertices making the particle surface [2].
You might find [3] helpful.

All the best, Vasileios

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/polyhedra/irregular.py#L26
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Polyhedra.GetSurfaces
[3] https://answers.launchpad.net/yade/+question/685375

weijie (amandajoe) said : #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:
#####################
from yade import polyhedra_utils,pack,plot,utils,export,qt,ymport
import numpy as np
import math
import random
#########################################################################Define physical parameters
powderDensity = 2500
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless'))
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[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
 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[2]]-vs[face11[1]]; vec1.normalize()
   vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize() #Normalize this object in-place.
   vects=vec1.cross(vec2)

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

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

   aa.append(vects[0])
   bb.append(vects[1])
   cc.append(vects[2])
   dd.append(dv)

  id1=id1+1
#######################################################################################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.01, avoidSelfInteractionMask=2),
 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)

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 : #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

> 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 : #6

Thanks Vasileios Angelidakis, that solved my question.

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[0],-1*vects[1],-1*vects[2]]
# dv=-1*dv

Best regards,
Vasileios

weijie (amandajoe) said : #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:

#####################
from yade import polyhedra_utils,pack,plot,utils,export,qt,ymport
import numpy as np
import math
import random
#########################################################################Define physical parameters
powderDensity = 2500
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless'))
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[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
 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[2]]-vs[face11[1]]; vec1.normalize()
   vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize() #Normalize this object in-place.
   vects=vec1.cross(vec2)

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

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

   aa.append(vects[0])
   bb.append(vects[1])
   cc.append(vects[2])
   dd.append(dv)

  id1=id1+1
#######################################################################################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.01, avoidSelfInteractionMask=2),
 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

Hi Jie,

Please find below some corrections to your script.
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

#####################
from yade import polyhedra_utils,pack,plot,utils,export,qt,ymport
import numpy as np
import math
import random
#########################################################################Define physical parameters
powderDensity = 2500
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless'))
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[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
 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[2]]-vs[face11[1]]; vec1.normalize()
   vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize() #Normalize this object in-place.
   vects=vec1.cross(vec2); vects.normalize()

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

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

   aa.append(vects[0])
   bb.append(vects[1])
   cc.append(vects[2])
   dd.append(dv)

  id1=id1+1

#######################################################################################
O.engines=[
        ForceResetter(),
        InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.01, avoidSelfInteractionMask=2),
        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)
###################

from yade import qt
v=qt.View()
v.sceneRadius=2.0

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 : #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 : #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:
#################################
from yade import polyhedra_utils,pack,plot,utils,qt
import random
import numpy as np

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e7,poisson=.2,density=2.5e3)
powderDensity = 2500
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless'))
#-------------------------------------------
#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()
    b.mask=1
    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[0] #s[0] stores center
    b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
    O.bodies.append(b)

def dvalue(vecn1,pp1):
 dd1=1*(vecn1[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
 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[2]]-vs[face11[1]]; vec1.normalize()
   vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize() #Normalize this object in-place.
   vects=vec1.cross(vec2); vects.normalize()
   dvalue2=dvalue(vects,vs[face11[0]]) # Some dvalue2 values return equal to 0. Check this part of your script once more.
   dv=dvalue2-chosenR
   aa.append(vects[0])
   bb.append(vects[1])
   cc.append(vects[2])
   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)

from yade import qt
v=qt.View()
O.saveTmp()

weijie (amandajoe) said : #12

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

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

#################################
from yade import polyhedra_utils,pack,plot,utils,qt
import random
import numpy as np

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e7,poisson=.2,density=2.5e3)
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=2500,label='frictionless'))
#-------------------------------------------
#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[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
 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[0] #s[0] stores center
 b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
 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[2]]-vs[face11[1]]; vec1.normalize()
   vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize() #Normalize this object in-place.
   vects=vec1.cross(vec2); vects.normalize()
   dvalue2=dvalue(vects,vs[face11[0]])
   aa.append(vects[0])
   bb.append(vects[1])
   cc.append(vects[2])
   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)

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

weijie (amandajoe) said : #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

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

from yade import polyhedra_utils,pack,plot,utils,qt
import random
import numpy as np

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e7,poisson=.2,density=2.5e3)
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=2500,label='frictionless'))
#-------------------------------------------
#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[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
 return dd1

for center,radius in sp:
 # Generate polyhedra
 color=Vector3(random.random(),random.random(),random.random())
 b=polyhedra_utils.polyhedra(material=n,size=(2*radius,radius,radius),color=color)
 b.state.pos = center #s[0] stores center
 b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
 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[2]]-vs[face11[1]]; vec1.normalize()
   vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize() #Normalize this object in-place.
   vects=vec1.cross(vec2); vects.normalize()
   dvalue2=dvalue(vects,vs[face11[0]])
   aa.append(vects[0])
   bb.append(vects[1])
   cc.append(vects[2])
   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)

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

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 [1], 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 [2], 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

[1] https://gitlab.com/groups/yade-dev/-/issues
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/PotentialBlock.cpp#L340

weijie (amandajoe) said : #16

Hi Vasileios, and thank you again.

I will try to do it.

Best regards,
Jie

Can you help with this problem?

Provide an answer of your own, or ask weijie for more information if necessary.

To post a message you must log in.