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:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
2020-03-30
Last query:
2020-03-30
Last reply:
2020-03-27

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.