The small particles generated in PotentialBlocks have more vertices than the initial polyhedra.
Hi,
When I converted the polyhedron to PotentialBlocks, the conversion was successful when the initial polyhedron size was larger, but I found that the polyhedron with a smaller particle size (three-dimensional length is less than about 0.005) converted more vertices than the initial polyhedron.This problem is shown in the code below.
Vasileios[1] (answer#15) thinks 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.
[1]https://answers.launchpad.net/yade/+question/689434
#####################
from yade import polyhedra_utils,pack,plot,utils,export,qt,ymport
import numpy as np
import math
import random
##########################################################################
# Material
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
#########################################################################
# Generate polyhedra
t= polyhedra_utils.polyhedra(m,size=(0.005,0.005,0.005),seed =2)
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
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.
aa.append(vects[0])
bb.append(vects[1])
cc.append(vects[2])
dd.append(dvalue2)
id1=id1+1
chosenR=min(dd)/2
print(dd)
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= t.state.ori
bbb.state.pos = t.state.pos+Vector3(0.01,0,0)
O.bodies.append(bbb)
#######################################################################################
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')
]
###################
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)
Edited by Vasileios Angelidakis