yade-dev team mailing list archive

Thread
Date

[Branch ~yade-dev/yade/trunk] Rev 2742: 1. implement packing for the WireMatPM

To: Yade developers <yade-dev@xxxxxxxxxxxxxxxxxxx>
From: noreply@xxxxxxxxxxxxx
Date: Wed, 16 Feb 2011 05:46:25 -0000
Reply-to: noreply@xxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx

------------------------------------------------------------
revno: 2742
committer: Klaus Thoeni <klaus.thoeni@xxxxxxxxx>
branch nick: yade
timestamp: Wed 2011-02-16 16:40:01 +1100
message:
  1. implement packing for the WireMatPM
  2. add some examples for the WireMatPM
added:
  examples/WireMatPM/
  examples/WireMatPM/wirepackings.py
  examples/WireMatPM/wiretensiltest.py
modified:
  py/pack/pack.py


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== added directory 'examples/WireMatPM'
=== added file 'examples/WireMatPM/wirepackings.py'
--- examples/WireMatPM/wirepackings.py	1970-01-01 00:00:00 +0000
+++ examples/WireMatPM/wirepackings.py	2011-02-16 05:40:01 +0000
@@ -0,0 +1,92 @@
+# encoding: utf-8
+from yade import pack, qt
+
+#### short description of script
+print 'This script shows the use of the function pack.hexaNet (interactions are not initialised)'
+
+#### define parameters for the net
+# mesh geometry
+mos = 0.08
+a = 0.04
+b = 0.04
+# wire diameter
+d = 2.7/1000.
+# net dimension
+Lx = 0.2
+Ly = 0.15
+
+# properties of particles
+radius = d*4.
+
+#### startAtCorner=True
+kw = {'color':[1,1,0],'wire':True,'highlight':False,'fixed':False,'material':-1}
+
+#### packing 1
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[0,0,0], xLength=Lx, yLength=Ly, mos=mos, a=a, b=b, startAtCorner=True, isSymmetric=True, **kw )
+O.bodies.append(netpack)
+print 'Packing 1:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### packing 2
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[0.4,0,0], xLength=Lx, yLength=Ly, mos=mos, a=a, b=b, startAtCorner=True, isSymmetric=False, **kw )
+O.bodies.append(netpack)
+print 'Packing 2:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### packing 3
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[0,-0.4,0], xLength=Lx, yLength=Ly+0.05, mos=mos, a=a, b=b, startAtCorner=True, isSymmetric=True, **kw )
+O.bodies.append(netpack)
+print 'Packing 3:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### packing 4
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[0.4,-0.4,0], xLength=Lx, yLength=Ly+0.05, mos=mos, a=a, b=b, startAtCorner=True, isSymmetric=False, **kw )
+O.bodies.append(netpack)
+print 'Packing 4:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### startAtCorner=False
+kw = {'color':[1,0,0],'wire':True,'highlight':False,'fixed':False,'material':-1}
+
+#### packing 1
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[1,0,0], xLength=Lx, yLength=Ly, mos=mos, a=a, b=b, startAtCorner=False, isSymmetric=True, **kw )
+O.bodies.append(netpack)
+print 'Packing 1:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### packing 2
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[1.4,0,0], xLength=Lx, yLength=Ly, mos=mos, a=a, b=b, startAtCorner=False, isSymmetric=False, **kw )
+O.bodies.append(netpack)
+print 'Packing 2:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### packing 3
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[1,-0.4,0], xLength=Lx, yLength=Ly+0.05, mos=mos, a=a, b=b, startAtCorner=False, isSymmetric=True, **kw )
+O.bodies.append(netpack)
+print 'Packing 3:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### packing 4
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[1.4,-0.4,0], xLength=Lx, yLength=Ly+0.05, mos=mos, a=a, b=b, startAtCorner=False, isSymmetric=False, **kw )
+O.bodies.append(netpack)
+print 'Packing 4:'
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+## to see it
+v=qt.Controller()
+v=qt.View()

=== added file 'examples/WireMatPM/wiretensiltest.py'
--- examples/WireMatPM/wiretensiltest.py	1970-01-01 00:00:00 +0000
+++ examples/WireMatPM/wiretensiltest.py	2011-02-16 05:40:01 +0000
@@ -0,0 +1,143 @@
+# -*- coding: utf-8 -*-
+from yade import utils, pack, plot
+
+#### short description of script
+print 'This script shows a tensile test of a net by using the UniaxialStrainer'
+
+#### define parameters for the net
+# wire diameter
+d = 2.7/1000.
+# particle radius
+radius = d*4.
+# define piecewise lineare stress-strain curve
+strainStressValues=[(0.0019230769,2.5e8),(0.0192,3.2195e8),(0.05,3.8292e8),(0.15,5.1219e8),(0.25,5.5854e8),(0.3,5.6585e8),(0.35,5.6585e8)]
+# elastic material properties
+particleVolume = 4./3.*pow(radius,3)*pi
+particleMass = 3.9/1000.
+density = particleMass/particleVolume
+young = strainStressValues[0][1] / strainStressValues[0][0]
+poisson = 0.3
+
+
+#### material definition
+netMat = O.materials.append( WireMat( young=young,poisson=poisson,density=density,isDoubleTwist=True,diameter=d,strainStressValues=strainStressValues,lambdaEps=0.4,lambdak=0.66) )
+
+wireMat = O.materials.append( WireMat( young=young,poisson=poisson,density=density,isDoubleTwist=False,diameter=3.4/1000,strainStressValues=strainStressValues ) )
+
+
+#### get net packing
+kw = {'color':[1,1,0],'wire':True,'highlight':False,'fixed':False,'material':netMat}
+[netpack,lx,ly] = pack.hexaNet( radius=radius, cornerCoord=[0,0,0], xLength=1.0, yLength=0.55, mos=0.08, a=0.04, b=0.04, startAtCorner=False, isSymmetric=True, **kw )
+O.bodies.append(netpack)
+print 'Real net length in x-direction [m]: ', lx
+print 'Real net length in y-direction [m]: ', ly
+
+
+#### get bodies for single wire at the boundary in y-direction and change properties
+bb = utils.uniaxialTestFeatures(axis=0)
+negIds,posIds=bb['negIds'],bb['posIds']
+
+for id in negIds:
+	O.bodies[id].material = O.materials[wireMat]
+	O.bodies[id].shape.color = [0,0,1]
+for id in posIds:
+	O.bodies[id].material = O.materials[wireMat]
+	O.bodies[id].shape.color = [0,0,1]
+
+
+#### define engines to create link
+interactionRadius=2.8
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='aabb')]), 
+	InteractionLoop(
+	[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='Ig2ssGeom')],
+	[Ip2_WireMat_WireMat_WirePhys(linkThresholdIteration=1,label='interactionPhys')],
+	[Law2_ScGeom_WirePhys_WirePM(linkThresholdIteration=1,label='interactionLaw')]
+	),
+	NewtonIntegrator(damping=0.),
+]
+
+#### define additional vertical interactions at the boundary
+utils.createInteraction(negIds[0],negIds[2])
+utils.createInteraction(negIds[3],negIds[4])
+utils.createInteraction(negIds[5],negIds[6])
+utils.createInteraction(negIds[7],negIds[1])
+utils.createInteraction(posIds[0],posIds[2])
+utils.createInteraction(posIds[3],posIds[4])
+utils.createInteraction(posIds[5],posIds[6])
+utils.createInteraction(posIds[7],posIds[1])
+
+#### time step definition for first time step to create links
+O.step()
+
+
+#### initialize values for UniaxialStrainer
+bb = utils.uniaxialTestFeatures(axis=1)
+negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
+strainRateTension = 0.1
+setSpeeds = True
+
+
+##### delete horizontal interactions for corner particles
+bb = utils.uniaxialTestFeatures(axis=1)
+negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
+
+
+##### delete some interactions
+O.interactions.erase(0,4)
+O.interactions.erase(0,5)
+O.interactions.erase(1,154)
+O.interactions.erase(1,155)
+O.interactions.erase(2,26)
+O.interactions.erase(2,27)
+O.interactions.erase(3,176)
+O.interactions.erase(3,177)
+
+#### time step definition for deleting some links which have been created by the Ig2 functor
+O.step()
+
+
+#### initializes now the interaction detection factor
+aabb.aabbEnlargeFactor=-1.
+Ig2ssGeom.interactionDetectionFactor=-1.
+
+
+#### define engines for simulation with UniaxialStrainer
+O.engines = O.engines[:3] + [ UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=1,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=True,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
+NewtonIntegrator(damping=0.5),
+PyRunner(initRun=True,iterPeriod=1000,command='addPlotData()'),
+]
+
+
+#### plot some results
+plot.plots={'un':('Fn',)}
+plot.plot(noShow=False, subPlots=False)
+
+def addPlotData():
+	Fn = 0.
+	for i in posIds:
+		try:
+			inter=O.interactions.withBody(i)[0]
+			F = abs(inter.phys.normalForce[1])
+		except:
+			F = 0
+		Fn += F
+	un = O.bodies[O.bodies[posIds[0]].id].state.pos[1] - O.bodies[O.bodies[posIds[0]].id].state.refPos[1]
+	if un > 0.10:
+		O.pause()
+	plot.addData( un=un*1000, Fn=Fn/1000 )
+
+
+#### time step definition for simulation
+## critical time step proposed by Bertrand
+kn = 16115042 # stiffness of single wire from code
+O.dt = 0.2*sqrt(particleMass/(2.*kn))
+
+
+#### to see it
+from yade import qt
+v = qt.Controller()
+v = qt.View()
+rr = qt.Renderer()
+rr.intrAllWire = True

=== modified file 'py/pack/pack.py'
--- py/pack/pack.py	2011-01-31 19:09:24 +0000
+++ py/pack/pack.py	2011-02-16 05:40:01 +0000
@@ -13,7 +13,7 @@
 * :ysrc:`scripts/test/pack-cloud.py`
 * :ysrc:`scripts/test/pack-predicates.py`
 * :ysrc:`examples/regular-sphere-pack/regular-sphere-pack.py`
-
+* :ysrc:`examples/WireMatPM/wirepackings.py`
 """
 
 import itertools,warnings
@@ -499,3 +499,93 @@
 	O.switchScene()
 	return ret
 
+def hexaNet( radius, cornerCoord=[0,0,0], xLength=1., yLength=0.5, mos=0.08, a=0.04, b=0.04, startAtCorner=True, isSymmetric=False, **kw ):
+	"""Definition of the particles for a hexagonal wire net in the x-y-plane for the WireMatPM.
+
+	:param radius: radius of the particle
+	:param cornerCoord: coordinates of the lower left corner of the net
+	:param xLenght: net length in x-direction
+	:param yLenght: net length in y-direction
+	:param mos: mesh opening size
+	:param a: length of double-twist 
+	:param b: height of single wire section
+	:param startAtCorner: if true the generation starts with a double-twist at the lower left corner
+	:param isSymmetric: defines if the net is symmetric with respect to the y-axis
+
+	:return: set of spheres which defines the net (net) and exact dimensions of the net (lx,ly).
+	
+	note::
+	This packing works for the WireMatPM only. The particles at the corner are always generated first.
+
+	"""
+	# check input dimension
+	if(xLength<mos): raise ValueError("xLength must be greather than mos!");
+	if(yLength<2*a+b): raise ValueError("yLength must be greather than 2*a+b!");
+	xstart = cornerCoord[0]
+	ystart = cornerCoord[1]
+	z = cornerCoord[2]
+	ab = a+b
+	# number of double twisted sections in y-direction and real length ly
+	ny = int( (yLength-a)/ab ) + 1
+	ly = ny*a+(ny-1)*b
+	jump=0
+	# number of sections in x-direction and real length lx
+	if isSymmetric:
+		nx = int( xLength/mos ) + 1
+		lx = (nx-1)*mos
+		if not startAtCorner:
+			nx+=-1
+	else:
+		nx = int( (xLength-0.5*mos)/mos ) + 1
+		lx = (nx-1)*mos+0.5*mos
+	net = []
+	# generate corner particles
+	if startAtCorner:
+		if (ny%2==0): # if ny even no symmetry in y-direction
+			net+=[utils.sphere((xstart,ystart+ly,z),radius=radius,**kw)] # upper left corner
+			if isSymmetric:
+				net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
+			else:
+				net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
+		else: # if ny odd symmetry in y-direction
+			if not isSymmetric:
+				net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
+				net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
+		jump=1
+	else: # do not start at corner
+		if (ny%2==0): # if ny even no symmetry in y-direction
+			net+=[utils.sphere((xstart,ystart,z),radius=radius,**kw)] # lower left corner
+			if isSymmetric:
+				net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
+			else:
+				net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
+		else: # if ny odd symmetry in y-direction
+			net+=[utils.sphere((xstart,ystart,z),radius=radius,**kw)] # lower left corner
+			net+=[utils.sphere((xstart,ystart+ly,z),radius=radius,**kw)] # upper left corner
+			if isSymmetric:
+				net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
+				net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
+		xstart+=0.5*mos
+	# generate other particles
+	if isSymmetric:
+		for i in range(ny):
+			y = ystart + i*ab
+			for j in range(nx):
+				x = xstart + j*mos
+				# add two particles of one vertical section (double-twist)
+				net+=[utils.sphere((x,y,z),radius=radius,**kw)]
+				net+=[utils.sphere((x,y+a,z),radius=radius,**kw)]
+			# set values for next section
+			xstart = xstart - 0.5*mos*pow(-1,i+jump)
+			nx = int(nx + 1*pow(-1,i+jump))
+	else:
+		for i in range(ny):
+			y = ystart + i*ab
+			for j in range(nx):
+				x = xstart + j*mos
+				# add two particles of one vertical section (double-twist)
+				net+=[utils.sphere((x,y,z),radius=radius,**kw)]
+				net+=[utils.sphere((x,y+a,z),radius=radius,**kw)]
+			# set values for next section
+			xstart = xstart - 0.5*mos*pow(-1,i+jump)
+	return [net,lx,ly]