yade-users team mailing list archive

[onyourself <question693455@xxxxxxxxxxxxxxxxxxxxx>]

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.