yade-users team mailing list archive

[Ali Rafiee <question685375@xxxxxxxxxxxxxxxxxxxxx>]

Question #685375 on Yade changed:
https://answers.launchpad.net/yade/+question/685375

    Status: Solved => Open

Ali Rafiee is still having a problem:
Thank you very much Vasileios,

I tried to do this converting method, and that is my code:

####################################################
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[0]-pp1[0],pp2[1]-pp1[1],pp2[2]-pp1[2]]
    if length(pp3)>0:
        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.001

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
    #print(len(vs))
    face2=[]
    face2=b.shape.GetSurfaces()
    #print(face2)
    #print(vs)
    
    id1=0
    while id1<len(face2):
        face11=face2[id1]
        if len(face11)>2:
            vec1=veccc(vs[face11[2]],vs[face11[1]])
            vec2=veccc(vs[face11[0]],vs[face11[1]])          
            
            vects=[]
            vects = list(np.cross(vec1, vec2))
            
            #print(vects)
            dvalue2=dvalue(vects,vs[face11[0]])
            
                
            #print("*******************************")            
            #print(dvalue2)
            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
    #print(aa)
    #print(bb)
    #print(cc)
    #print(dd)
    
    
    #### 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.001, r=chosenR, R=chosenR, a=aa, b=bb, c=cc, d=dd, id=len(O.bodies))
    utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia,material='frictionless', pos=[0,0,0], fixed=False)
    bbb.state.pos = [-2,0,0]
    lidID = O.bodies.append(bbb)
    
    
 ###########################################################################

the values for a, b, c and d are obtained, but I can't complete all
parameters required by PotentialBlock() command  just with these 4
parameters.

in any case, I do not know how to complete the last phase of this
geometric transformation.

Thanks again
Ali

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.