Question #685375 on Yade changed:
https://answers.launchpad.net/yade/+question/685375
Status: Open => Answered
Vasileios Angelidakis proposed the following answer:
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
Some comments on your script:
- 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[0],-1*vects[1],-1*vects[2]]
> 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
Below some updates I would apply on your script.
All the best,
Vasileios
####################################################
from yade import pack,export,qt,ymport
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[0]-pp1[0],pp2[1]-pp1[1],pp2[2]-pp1[2]]
# if pp3.norm()>0: # if length(pp3)>0:
# pp3.normalize() #
pp3=[pp3[0]/length(pp3),pp3[1]/length(pp3),pp3[2]/length(pp3)]
# return pp3
####################################
### calculate d value of a surface with normal vector and one point
def dvalue(vecn1,pp1):
dd1=-1*(vecn1[0]*pp1[0]+vecn1[1]*pp1[1]+vecn1[2]*pp1[2])
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[2]]-vs[face11[1]]; vec1.normalize()
#vec1=veccc(vs[face11[2]],vs[face11[1]])
vec2=vs[face11[0]]-vs[face11[1]]; vec2.normalize()
#vec2=veccc(vs[face11[0]],vs[face11[1]])
# vects=[]
vects=vec1.cross(vec2) #vects = list(np.cross(vec1,
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
#### 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)
###########################################################################
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
