Question #689434 on Yade changed:
https://answers.launchpad.net/yade/+question/689434
Status: Open => Answered
Vasileios Angelidakis proposed the following answer:
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)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
