Question #693455 on Yade changed:
https://answers.launchpad.net/yade/+question/693455
Status: Answered => Open
onyourself is still having a problem:
hi,
1) - it just works, just you wrongly interpret it (typically it is not
"visible in GUI" but observable according to numerical results)
if so, maybe i just need more time to see the difference? i didnt print
anything so the only way i can be sure that my code works is to check 3D
view by my eyes…i thought it should be visible as my specimen is actual
incompact.
2)got it, thanks, below is my code.
####################################
from yade.gridpfacet import *
from yade import pack, plot
from random import random
import numpy as np
from numpy import *
import math
import os
#parameters
rParticle = 0.14e-3#diameter0.28
rRelFuzz=.0005
mi,ma = (-500e-3,0.,0.),(500e-3,0,250e-3)
#nCyls,nSphs = 8000,28216
frictionAngleSph=30
frictionAngleCyl=30
width = 1000.e-3,
length = 10e-3,
height = 250.e-3,
# default parameters or from table
readParamsFromTable(noTableOk=True,
# type of test ['cyl','cube']
testType = 'cube',
# material parameters
young = 20e9,
poisson = .2,
frictionAngle = 1.2,
sigmaT = 1.5e6,
epsCrackOnset = 1e-4,
relDuctility = 30,
# prestress
preStress = -3e9,
# axial strain rate
strainRate = -100,
# assamlby parameters
rParticle = .075e-3, #
width = 1000e-3,
height = 250e-3,
bcCoeff = 5,
# facets division
nw = 240,
nh = 150,
# output specifications
fileName = 'test',
exportDir = '/tmp',
runGnuplot = False,
runInGui = True,
)
from yade.params.table import *
#meterials########
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=20,density=910e+6,label='sphcylMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6,label='cylinderMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.15,frictionAngle=30,density=2630e+6,label='sphereMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.15,frictionAngle=0,density=910e+6,label='walls'))
#O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6,label='frictMat'))
frictMat = O.materials.append(FrictMat(
young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6
))
##sphere
sp = yade.pack.SpherePack()
sp.makeCloud(mi,ma,porosity=0.7/1.7,psdSizes=[.00014,0.00016,0.00022,.0003,.00035],psdCumm=[0.,0.1,0.3,0.6,1.],num=35000)
#psdSizes=[0.01,0.0125,0.015,0.018,0.02,0.023,0.025,0.027,0.029,0.03],psdCumm=[0,0.05,0.15,0.38,0.54,0.77,0.9,0.97,0.99,1])
# makeCloud for spheres
spheres=sp.toSimulation(color=(0,0.5,0.7),material='sphereMat')
print "cheers"
num=0
vol=0
mass=0
area=0
for b in O.bodies:
if isinstance(b.shape,Sphere):
r=b.shape.radius
x=b.state.pos[0]
y=b.state.pos[1]
z=b.state.pos[2]
ID=b.id
num+=1
#vol+=4./3.*np.pi*r**3
mass+=2.63*4./3.*np.pi*r**3
area+=np.pi*r**2
print 'num,mass,voidRatio =',num,mass,((1000e-3*250e-3)-area)/area
length=6e-3
rfiber=0.023e-3/2
fibre_mass=np.pi*(rfiber)**2*length*0.91
Xw=0.5*1e-2 # fiber content
numFiber=Xw*mass/fibre_mass
print "the numer of fibres =",int(numFiber)
Ne=int(length/(0.28e-3/2))
print 'Ne=',Ne
nodesIds=[]
cylIds=[]
numFibre=0
#target_num=int(numFiber)
target_num=1200
deleted_sphere=0
####fiberfilewrite
fiberwr=open('fibers.txt','w')
cylNodes=[]
for n in range(2000):
if numFibre<target_num:
random.seed()
l_x=random.uniform(-500e-3,500e-3)
l_z=random.uniform(0,250e-3)
x0=l_x
z0=l_z
u1 = random.random()
u2 = random.random()*pi
u3 = random.random()*pi
a = sqrt(1-u1)
b = sqrt(u1)
t=((a*cos(u2))**2+(b*sin(u3))**2)**0.5
cosphi=a*cos(u2)/t
sinphi=b*sin(u3)/t
hx=length*cosphi
hz=length*sinphi
x1=x0+hx
y1=0
z1=z0+hz
if -500e-3<x1<500e-3 and 0<z1<250e-3 and -500e-3<x0<500e-3 and
0<z0<250e-3:
if len(cylNodes)<=1:
fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(0)+'\t'+str(z1)+'\n')
cylNodes.append([x0,z0,x1,z1])
numFibre+=1
vertices=[]
for i in range(0, Ne+1):
px=float(i)*hx/float(Ne)+x0; py=0;
pz=float(i)*hz/float(Ne)+z0;
vertices.append([px,py,pz])
cylinderConnection(vertices,0.023e-3/2,nodesIds,cylIds,color=[0.5,0.5,0],highlight=False,intMaterial='cylinderMat',extMaterial='fMat')
print "+1"
else:
valid=True
#print len(cylNodes)
for c in cylNodes:
if(max(c[0],c[2])>=min(x0,x1)
and max(x0,x1)>=min(c[0],c[2])
and max(c[1],c[3])>=min(z0,z1)
and max(z0,z1)>=min(c[1],c[3])):
x12=(x1-x0)
z12=(z1-z0)
x13=(c[0]-x0)
z13=(c[1]-z0)
x14=(c[2]-x0)
z14=(c[3]-z0)
cross1=x12*z13-x13*z12
cross2=x12*z14-x14*z12
x34=c[2]-c[0]
z34=c[3]-c[1]
#x23=c[0]-x1
#z23=c[1]-z1
x32=x1-c[0]
z32=z1-c[1]
#x24=c[2]-x1
#z24=c[3]-z1
x31=x0-c[0]
z31=z0-c[1]
cross3=x34*z32-x32*z34
cross4=x34*z31-x31*z34
if(cross1*cross2<=0
and cross3*cross4<=0):
valid=False
print "next"
break
if valid:
fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(0)+'\t'+str(z1)+'\n')
cylNodes.append([x0,z0,x1,z1])
numFibre+=1
#print "++1"
vertices=[]
for i in range(0, Ne+1):
px=float(i)*hx/float(Ne)+x0;
py=0; pz=float(i)*hz/float(Ne)+z0;
vertices.append([px,py,pz])
cylinderConnection(vertices,0.023e-3/2,nodesIds,cylIds,color=[0.5,0.5,0],highlight=False,intMaterial='cylinderMat',extMaterial='fMat')
###spherefileWrite
##################################################################
sphwr=open('spheres.txt','w')
for b in O.bodies:
if isinstance(b.shape,Sphere):
r=b.shape.radius
if r>0.012e-3:
x=b.state.pos[0]
y=b.state.pos[1]
z=b.state.pos[2]
ID=b.id
k=(z1-z0)/(x1-x0)
#kx-z+z0-kx0=0
A=k
B=-1
C=z0-k*x0
#if
(A*x+B*z+C)/(A^2+B^2)^0.5<=(0.14e-3+0.023e-3/2):
if
abs(k*x-1*z+C)/(k**2+1)**0.5<=(r+0.0115e-3):
O.bodies.erase(b.id)
deleted_sphere+=1
#else:
sphwr.write(str(ID)+'\t'+str(x)+'\t'+str(y)+'\t'+str(z)+'\t'+str(r)+'\n')
sphwr.close()
##################################################################
fiberwr.close()
print 'numFibre , deleted_sphere =',numFibre, deleted_sphere
for i in O.bodies:
if isinstance(i.shape,Sphere):
i.state.blockedDOFs="yXZ"
#i.shape.color=(0.,0.5,0.7)
facets = []
nw2 = nw/40
nw2 = 240
for r in xrange(nw2):
for h in xrange(nh):
v11 = Vector3( -.5*width + (r+0)*width/nw2, 5e-3, 0 )
v12 = Vector3( -.5*width + (r+1)*width/nw2, -5e-3, 0 )
v13 = Vector3( -.5*width + (r+1)*width/nw2, 5e-3, 0 )
v14 = Vector3( -.5*width + (r+0)*width/nw2, -5e-3, 0 )
f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
v21 = Vector3( +.5*width, -5e-3, height*(h+0)/float(nh) )
v22 = Vector3( +.5*width, 5e-3, height*(h+0)/float(nh) )
v23 = Vector3( +.5*width, 5e-3, height*(h+1)/float(nh) )
v24 = Vector3( +.5*width, -5e-3, height*(h+1)/float(nh) )
f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
v31 = Vector3( +.5*width - (r+0)*width/nw2, +5e-3, 250e-3 )
v32 = Vector3( +.5*width - (r+1)*width/nw2, +5e-3, 250e-3 )
v33 = Vector3( +.5*width - (r+1)*width/nw2, -5e-3, 250e-3 )
v34 = Vector3( +.5*width - (r+0)*width/nw2, -5e-3, 250e-3 )
f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
v41 = Vector3( -.5*width, -5e-3, height*(h+0)/float(nh) )
v42 = Vector3( -.5*width, 5e-3, height*(h+0)/float(nh) )
v43 = Vector3( -.5*width, 5e-3, height*(h+1)/float(nh) )
v44 = Vector3( -.5*width, -5e-3, height*(h+1)/float(nh) )
f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
f.state.mass = mass
f.shape.color=(200,255,240)
f.state.blockedDOFs = 'XYZz'
'''
#preStress
# apply prestress to facets
def addForces():
for f in facets:
n = f.shape.normal
a = f.shape.area
O.forces.addF(f.id,preStress*a*n)
'''
for f in facets:
n = f.shape.normal
a = f.shape.area
#O.forces.addF(f.id,preStress*a*n)
O.forces.setPermF(f.id,preStress*a*n)
O.engines=[
ForceResetter(),
InsertionSortCollider([
#Bo1_Box_Aabb(),
Bo1_Sphere_Aabb(),
Bo1_GridConnection_Aabb(),
Bo1_Facet_Aabb(),
]),
InteractionLoop([
Ig2_Sphere_Sphere_ScGeom(),
#Ig2_Box_Sphere_ScGeom(),
Ig2_GridNode_GridNode_GridNodeGeom6D(),
Ig2_Sphere_GridConnection_ScGridCoGeom(),
Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
Ig2_Facet_Sphere_ScGeom(),
],
[
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
# internal cylinder physics
Ip2_FrictMat_FrictMat_FrictPhys() # physics for external
interactions, i.e., cylinder-cylinder, sphere-sphere, cylinder-sphere
#Ip2_FrictMat_CpmMat_FrictPhys(),
],
[
Law2_ScGeom_FrictPhys_CundallStrack(), # contact law for
sphere-sphere
Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law
for cylinder-sphere
Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law
for "internal" cylinder forces
Law2_GridCoGridCoGeom_FrictPhys_CundallStrack() # contact law
for cylinder-cylinder interaction
#Law2_ScGeom_CpmPhys_Cpm(),
]
),
GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
#triax,
#TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.1,label='newton'),
#PyRunner(iterPeriod=200,command='targetSpecimen()',label='target'),
#PyRunner(iterPeriod=20,command='history()',label='recorder'),
]
######################################################
--
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
