yade-users team mailing list archive

[Vasileios Angelidakis <question685375@xxxxxxxxxxxxxxxxxxxxx>]

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.