def setUp(self):
		woo.master.reset()
		woo.master.scene=S=Scene(fields=[DemField()])
		r1,r2,p0,p1=1,.5,Vector3.Zero,Vector3(0,0,3)
		S.dem.par.addClumped([
			utils.sphere(p0,r1),
			utils.sphere(p1,r2)
		])
		for n in (S.dem.par[0].shape.nodes[0],S.dem.par[1].shape.nodes[0]): S.dem.nodesAppend(n);
		self.bC,self.b1,self.b2=S.dem.nodes

def randomPeriPack(radius,initSize,rRelFuzz=0.0,memoizeDb=None):
	"""Generate periodic dense packing.

	A cell of initSize is stuffed with as many spheres as possible, then we run periodic compression with PeriIsoCompressor, just like with randomDensePack.

	:param radius: mean sphere radius
	:param rRelFuzz: relative fuzz of sphere radius (equal distribution); see the same param for randomDensePack.
	:param initSize: initial size of the periodic cell.

	:return: SpherePack object, which also contains periodicity information.
	"""
	from math import pi
	from woo import core, dem
	sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,initSize[0],initSize[1],initSize[2],fullDim=Vector3(0,0,0),wantPeri=True,fillPeriodic=False,spheresInCell=-1,memoDbg=True)
	if sp: return sp
	#oldScene=O.scene
	S=core.Scene(fields=[dem.DemField()])
	sp=SpherePack()
	S.periodic=True
	S.cell.setBox(initSize)
	sp.makeCloud(Vector3().Zero,S.cell.size0,radius,rRelFuzz,-1,True)
	from woo import log
	log.setLevel('PeriIsoCompressor',log.DEBUG)
	S.engines=[dem.ForceResetter(),dem.InsertionSortCollider([dem.Bo1_Sphere_Aabb()],verletDist=.05*radius),dem.ContactLoop([dem.Cg2_Sphere_Sphere_L6Geom()],[dem.Cp2_FrictMat_FrictPhys()],[dem.Law2_L6Geom_FrictPhys_IdealElPl()],applyForces=True),dem.PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='print "done"; S.stop();',globalUpdateInt=20,keepProportions=True),dem.Leapfrog(damping=.8)]
	mat=dem.FrictMat(young=30e9,tanPhi=.1,ktDivKn=.3,density=1e3)
	for s in sp: S.dem.par.add(utils.sphere(s[0],s[1],mat=mat))
	S.dt=utils.pWaveDt(S)
	#O.timingEnabled=True
	S.run(); S.wait()
	ret=SpherePack()
	ret.fromDem(S,S.dem)
	_memoizePacking(memoizeDb,ret,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0)) # fullDim unused
	return ret

	def setUp(self):
		woo.master.scene=S=Scene(fields=[DemField()])
		S.periodic=True;
		S.cell.setBox(2.5,2.5,3)
		self.cellDist=Vector3i(0,0,10) # how many cells away we go
		self.relDist=Vector3(0,.999999999999999999,0) # rel position of the 2nd ball within the cell
		self.initVel=Vector3(0,0,5)
		S.dem.par.add(utils.sphere((1,1,1),.5))
		self.initPos=Vector3([S.dem.par[0].pos[i]+self.relDist[i]+self.cellDist[i]*S.cell.hSize0.col(i).norm() for i in (0,1,2)])
		S.dem.par.add(utils.sphere(self.initPos,.5))
		S.dem.par[1].vel=self.initVel
		S.engines=[Leapfrog(reset=True)]
		S.cell.nextGradV=Matrix3(0,0,0, 0,0,0, 0,0,-1)
		S.cell.homoDeform=Cell.HomoVel2
		S.dem.collectNodes() # avoid msg from Leapfrog
		S.dt=0 # do not change positions with dt=0 in NewtonIntegrator, but still update velocities from velGrad

 def setUp(self):
     woo.master.reset()
     woo.master.scene.fields=[DemField()]
     S=woo.master.scene
     self.count=100
     S.dem.par.add([utils.sphere([random.random(),random.random(),random.random()],random.random()) for i in range(0,self.count)])
     random.seed()

	def testNoCollide(self):
		"Clump: particles inside one clump don't collide with each other"
		# use a new scene, with a different clump in this test
		S=Scene(fields=[DemField()])
		S.dem.par.addClumped([utils.sphere(c,r) for c,r in [((1,0,0),1),((0,1,0),1),((0,0,1),1)]])
		S.engines=[InsertionSortCollider([Bo1_Sphere_Aabb()])]
		S.one()
		for i,j in [(0,1),(1,2),(0,2)]:
			self.assert_(not woo.dem.Collider.mayCollide(S.dem,S.dem.par[i],S.dem.par[j]))
		self.assert_(len(S.dem.con)==0)

def regularOrtho(predicate,radius,gap,**kw):
	"""Return set of spheres in regular orthogonal grid, clipped inside solid given by predicate.
	Created spheres will have given radius and will be separated by gap space."""
	ret=[]
	mn,mx=predicate.aabb()
	if(max([mx[i]-mn[i] for i in 0,1,2])==float('inf')): raise ValueError("Aabb of the predicate must not be infinite (didn't you use union | instead of intersection & for unbounded predicate such as notInNotch?");
	xx,yy,zz=[arange(mn[i]+radius,mx[i]-radius,2*radius+gap) for i in 0,1,2]
	for xyz in itertools.product(xx,yy,zz):
		if predicate(xyz,radius): ret+=[utils.sphere(xyz,radius=radius,**kw)]
	return ret

def fill_cylinder_with_spheres(sphereRadius,cylinderRadius,cylinderHeight,cylinderOrigin,cylinderSlope):
    spheresCount=0
    for h in range(0,cylinderHeight/sphereRadius/2):
            for r in range(1,cylinderRadius/sphereRadius/2):
                dfi = asin(0.5/r)*2
                for a in range(0,6.28/dfi):
                    x = cylinderOrigin[0]+2*r*sphereRadius*cos(dfi*a)
                    y = cylinderOrigin[1]+2*r*sphereRadius*sin(dfi*a)
                    z = cylinderOrigin[2]+h*2*sphereRadius
                    o.bodies.append(utils.sphere([x,y*cos(cylinderSlope)+z*sin(cylinderSlope),z*cos(cylinderSlope)-y*sin(cylinderSlope)],sphereRadius))
                    spheresCount+=1
    return spheresCount

def regularHexa(predicate,radius,gap,**kw):
	"""Return set of spheres in regular hexagonal grid, clipped inside solid given by predicate.
	Created spheres will have given radius and will be separated by gap space."""
	ret=[]
	a=2*radius+gap
	# thanks to Nasibeh Moradi for finding bug here:
	# http://www.mail-archive.com/[email protected]/msg01424.html
	hy,hz=a*sqrt(3)/2.,a*sqrt(6)/3.
	mn,mx=predicate.aabb()
	dim=mx-mn
	if(max(dim)==float('inf')): raise ValueError("Aabb of the predicate must not be infinite (didn't you use union | instead of intersection & for unbounded predicate such as notInNotch?");
	ii,jj,kk=[range(0,int(dim[0]/a)+1),range(0,int(dim[1]/hy)+1),range(0,int(dim[2]/hz)+1)]
	for i,j,k in itertools.product(ii,jj,kk):
		x,y,z=mn[0]+radius+i*a,mn[1]+radius+j*hy,mn[2]+radius+k*hz
		if j%2==0: x+= a/2. if k%2==0 else -a/2.
		if k%2!=0: x+=a/2.; y+=hy/2.
		if predicate((x,y,z),radius): ret+=[utils.sphere((x,y,z),radius=radius,**kw)]
	return ret

import woo
from woo import utils,pack,plot
from woo.dem import *
from woo.core import *

woo.master.scene=S=Scene(fields=[DemField(gravity=(0,0,-10))])
mat=woo.dem.PelletMat(young=1e6,tanPhi=.5,ktDivKn=.2,density=1000)
if 1:
	sp=pack.SpherePack()
	sp.makeCloud((0,0,0),(10,10,10),.4,rRelFuzz=.5)
	sp.toSimulation(S,mat=mat)
else:
	S.dem.par.add(utils.sphere((0,0,1),.5,mat=mat))

S.dem.par.add(utils.wall(0,axis=2,sense=1,mat=mat))
S.engines=utils.defaultEngines(damping=0.,cp2=Cp2_PelletMat_PelletPhys(),law=Law2_L6Geom_PelletPhys_Pellet(plastSplit=True))+[
	PyRunner(1,'S.plot.addData(i=S.step,t=S.time,Eerr=(S.energy.relErr() if S.step>100 else 0),**S.energy)'),
]
	
S.dt=.3*utils.pWaveDt(S)

S.dem.collectNodes()
S.trackEnergy=True
S.saveTmp()
S.plot.plots={'i':(S.energy,None,('Eerr','g--'))}
S.plot.plot()
S.run(500)
#from woo import gl
#gl.Gl1_Wall.div=10
#gl.Gl1_InfCylinder.wire=True

#!/usr/local/bin/woo-trunk -x
# -*- coding: utf-8 -*-

o=Omega() 
o.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),]),
    IGeomDispatcher([Ig2_Sphere_Sphere_Dem3DofGeom()]),
    IPhysDispatcher([Ip2_2xFrictMat_CSPhys()]),
    LawDispatcher([Law2_Dem3Dof_CSPhys_CundallStrack()]),
    GravityEngine(gravity=[0,0,-9.81]),
    NewtonIntegrator(damping = 0.01)
]

from woo import utils

o.bodies.append(utils.sphere([0,0,6],1,dynamic=True, color=[0,1,0]))
o.bodies.append(utils.sphere([0,0,0],1,dynamic = False, color=[0,0,1]))
o.dt=.2*utils.PWaveTimeStep()

from woo import qt
qt.Controller()
qt.View()

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]

from woo.core import*
from woo.dem import *
import woo.gl
from woo import utils
S=woo.master.scene=Scene(fields=[DemField(gravity=(0,0,-10))])
S.dem.par.add([
	utils.wall((0,0,0),axis=2,sense=1),
	utils.sphere((0,0,1),.2)
])
S.dem.par[1].vel=(0,1,0)
S.dt=.7*utils.pWaveDt()
S.engines=utils.defaultEngines()+[woo.dem.Tracer(num=512,compress=2,stepPeriod=10,compSkip=1)]
S.saveTmp()
import woo.qt
woo.qt.View()
S.run(2000)

import sys
sys.path.append('.')
import clDem

from minieigen import *

from woo import utils
m=utils.defaultMaterial()
O.dem.par.append([
	utils.wall((0,0,0),axis=2,material=m,fixed=True),
	utils.sphere((0,0,1.3),radius=1,material=m)
])
O.dem.par[-1].vel=(-1.,0,0)
O.dem.par[-1].angVel=(0,1,0)
# O.dem.par[-1].material.tanPhi=0. # no friction
O.scene.dt=.003*utils.pWaveDt()
O.scene.trackEnergy=True
from woo.dem import *
from woo.core import*
O.scene.engines=[
	#PyRunner('if len(O.dem.con)>0 and O.dem.con[0].real: O.pause()'),
	Gravity(gravity=(0,0,-10)),
	Leapfrog(damping=.05,reset=True),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
	ContactLoop([Cg2_Sphere_Sphere_L6Geom(),Cg2_Wall_Sphere_L6Geom()],[Cp2_FrictMat_FrictPhys()],[Law2_L6Geom_FrictPhys_IdealElPl()],applyForces=True),
]
import woo.qt
woo.qt.View()

O.scene.clDev=(1,0) # intel
sim=woo.cld.CLDemField.wooToClDem(O.scene,stepPeriod=1,relTol=-1)

from woo import utils
from woo.core import *
from woo.dem import *
from woo import plot
from woo import *
import woo.log
woo.log.setLevel('LawTester',woo.log.INFO)
m=utils.defaultMaterial()
S=woo.master.scene=Scene(fields=[DemField()])
S.dem.par.add([
	utils.sphere((0,0,0),.5,fixed=False,wire=True,mat=m),
	utils.sphere((0,1.01,0),.5,fixed=False,wire=True,mat=m)
])
S.dem.collectNodes()
S.engines=utils.defaultEngines(damping=.01)+[
	LawTester(ids=(0,1),abWeight=.3,smooth=1e-4,stages=[
			LawTesterStage(values=(-1,0,0,0,0,0),whats='v.....',until='bool(C)',done='print "Stage finished, at step",stage.step,", contact is",C'),
			LawTesterStage(values=(-.01,0,0,0,0,0),whats='v.....',until='C and C.geom.uN<-1e-2',done='print "Compressed to",C.geom.uN'),
			LawTesterStage(values=(.01,0,0,0,0,0),whats='v.....',until=('not C'),done='print "Contact broken",S.step;'),
			LawTesterStage(values=(-1e5,0,0,0,0,0),whats='fvv...',until=('C and tester.fErrRel[0]<1e-1'),done='print "Force-loaded contact stabilized";'),
			LawTesterStage(values=(-1e5,-.01,0,0,0,0),whats='fvvvvv',until=('"plast" in E and E["plast"]>0'),done='print "Plastic sliding reached";'),
			LawTesterStage(values=(0,-.01,0,0,0,0),whats=('vvvvvv'),until='stage.step>5000',done='print "5000 steps sliding done";'),
			LawTesterStage(values=(0,0,0,0,0,0),whats='vvvvvv',until='stage.step>100',done='E["plast"]=0.'),
			LawTesterStage(values=(0,0,-.01,0,0,0),whats='vvvvvv',until=('E["plast"]>1000'),done='print "sliding in the z-direction reached"'),
			LawTesterStage(values=(0,0,0,.1,0,0),whats='vvvvvv',until='stage.step>10000',done='print "Twist done"'),
			LawTesterStage(values=(0,0,0,0,.1,0),whats='vvvvvv',until='stage.step>10000',done='print "Bending done"'),
			LawTesterStage(values=(0,0,0,0,-.1,0),whats='vvvvvv',until='stage.step>10000',done='print "Bending back"'),
		],
		done='tester.dead=True; S.stop(); print "Everything done, making myself dead and pausing."',
		label='tester'
	),

# 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=False,diameter=d,strainStressValues=strainStressValues,lambdaEps=0.4,lambdak=0.21))


#### create boddies, default: dynamic=True
O.bodies.append( utils.sphere([0,0,0], radius, wire=False, color=[1,0,0], highlight=False, material=netMat) )
O.bodies.append( utils.sphere([0,a,0], radius, wire=False, color=[0,1,0], highlight=False, material=netMat) )

FixedSphere=O.bodies[0]
MovingSphere=O.bodies[1]

FixedSphere.dynamic=False
MovingSphere.dynamic=False

def addPlotData():
    if O.iter < 1:
        plot.addData( Fn=0., un=0. )
        #plot.saveGnuplot('net-2part-displ-unloading')
    else:
        try:
            i=O.interactions[FixedSphere.id,MovingSphere.id]

				sp.rotate(*orientation.toAxisAngle())
			return filterSpherePack(predicate,sp,mat=mat)
		else: print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
		sys.stdout.flush()
	S=core.Scene(fields=[dem.DemField()])
	if wantPeri:
		# x1,y1,z1 already computed above
		sp=SpherePack()
		S.periodic=True
		S.cell.setBox(x1,y1,z1)
		#print cloudPorosity,beta,gamma,N100,x1,y1,z1,S.cell.refSize
		#print x1,y1,z1,radius,rRelFuzz
		S.engines=[dem.ForceResetter(),dem.InsertionSortCollider([dem.Bo1_Sphere_Aabb()],verletDist=.05*radius),dem.ContactLoop([dem.Cg2_Sphere_Sphere_L6Geom()],[dem.Cp2_FrictMat_FrictPhys()],[dem.Law2_L6Geom_FrictPhys_IdealElPl()],applyForces=True),dem.Leapfrog(damping=.7,reset=False),dem.PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='print "DONE"; S.stop();',globalUpdateInt=5,keepProportions=True,label='compressor')]
		num=sp.makeCloud(Vector3().Zero,S.cell.size0,radius,rRelFuzz,spheresInCell,True)
		mat=dem.FrictMat(young=30e9,tanPhi=.5,density=1e3,ktDivKn=.2)
		for s in sp: S.dem.par.add(utils.sphere(s[0],s[1],mat=mat))
		S.dem.collectNodes()
		S.dt=.5*utils.pWaveDt(S)
		S.run(); S.wait()
		sp=SpherePack(); sp.fromDem(S,S.dem)
		#print 'Resulting cellSize',sp.cellSize,'proportions',sp.cellSize[1]/sp.cellSize[0],sp.cellSize[2]/sp.cellSize[0]
		# repetition to the required cell size will be done below, after memoizing the result
	else:
		raise RuntimError("Aperiodic compression not implemented.")
		assumedFinalDensity=0.6
		V=(4/3)*pi*radius**3; N=assumedFinalDensity*fullDim[0]*fullDim[1]*fullDim[2]/V;
		TriaxialTest(
			numberOfGrains=int(N),radiusMean=radius,radiusStdDev=rRelFuzz,
			# upperCorner is just size ratio, if radiusMean is specified
			upperCorner=fullDim,
			## no need to touch any the following

# The sphere couple is oriented randomly, but the result should always be the same.
#
from woo import utils,plot
import random, woo.log
random.seed()
#woo.log.setLevel('LawTester',woo.log.TRACE)

# sphere's radii
r1,r2=.1,.2
# place sphere 1 at the origin
pt1=Vector3(0,0,0)
# random orientation of the interaction
normal=Vector3(random.random()-.5,random.random()-.5,random.random()-.5)
normal=Vector3.UnitX
O.bodies.append([
    utils.sphere(pt1,r1,wire=True,color=(.7,.7,.7)),
    utils.sphere(pt1+.999999*(r1+r2)*normal.normalized(),r2,wire=True,color=(0,0,0))
])

O.engines=[
    ForceResetter(),
    PyRunner(iterPeriod=1,command='import time; time.sleep(.05)'),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        #[Ig2_Sphere_Sphere_ScGeom()],    [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()] # ScGeom
        #[Ig2_Sphere_Sphere_L3Geom(approxMask=63)],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True,noSlip=False)] # L3Geom
        [Ig2_Sphere_Sphere_L6Geom(approxMask=63)],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L6Geom_FrictPhys_Linear(charLen=1)] # L6Geom
    ),
    LawTester(ids=[0,1],disPath=[(0,0,0)]*7+[(-1e-5,0,0),(-1e-5,.1,.1)],rotPath=[(0,.2,0),(0,0,0),(0,0,.2),(0,0,0),(.2,0,0),(-.2,0,0),(0,0,0)],pathSteps=[10],doneHook='tester.dead=True; O.pause();',label='tester',rotWeight=0),
    #LawTester(ids=[0,1],path=[
    #    (-1e-5,0,0),(-.1,0,0),(-.1,.1,0),(-1e-5,.1,0), # towards, shear, back to intial normal distance

sphereRadius = 0.01
nbSpheres = (32, 11, 32)
print("Creating %d spheres..." % (nbSpheres[0] * nbSpheres[1] * nbSpheres[2]), end=" ")
for i in range(nbSpheres[0]):
    for j in range(nbSpheres[1]):
        for k in range(nbSpheres[2]):
            x = (i * 2 - nbSpheres[0]) * sphereRadius * 1.1 + sphereRadius * random.uniform(-0.1, 0.1)
            y = -j * sphereRadius * 2.2 - 0.01
            z = (k * 2 - nbSpheres[2]) * sphereRadius * 1.1 + sphereRadius * random.uniform(-0.1, 0.1)
            r = random.uniform(sphereRadius, sphereRadius * 0.9)
            dynamic = True
            color = [0.51, 0.52, 0.4]
            if i == 0 or i == nbSpheres[0] - 1 or j == nbSpheres[1] - 1 or k == 0 or k == nbSpheres[2] - 1:
                dynamic = False
                color = [0.21, 0.22, 0.1]
            O.bodies.append(utils.sphere([x, y, z], r, color=color, dynamic=dynamic))
print("done\n")

## Estimate time step
# O.dt=utils.PWaveTimeStep()
O.dt = 0.0001


## Engines
O.engines = [
    ## Resets forces and momenta the act on bodies
    ForceResetter(),
    ## Using bounding boxes find possible body collisions.
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
        # [Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],

 def setUp(self):
     woo.master.scene=S=Scene(fields=[DemField()])
     S.engines=utils.defaultEngines()
     S.dem.par.add(utils.sphere((0,0,0),radius=1))

	NewtonIntegrator(damping=0.0),
	###
	### NOTE this extra engine:
	###
	### You want snapshot to be taken every 1 sec (realTimeLim) or every 50 iterations (iterLim),
	### whichever comes soones. virtTimeLim attribute is unset, hence virtual time period is not taken into account.
	PyRunner(iterPeriod=1,command='myAddPlotData()')
]

## define and append material
mat=FrictMat(young=600.0e6,poisson=0.6,density=2.60e3,frictionAngle=26,label='Friction')
O.materials.append(mat)

## create two spheres (one will be fixed) and append them
from woo import utils
s0=utils.sphere([0,0,0],1,color=[0,1,0],fixed=False,wire=True,material='Friction')
s1=utils.sphere([2,0,0],1,color=[0,2,0],fixed=True,wire=True,material='Friction')
O.bodies.append(s0)
O.bodies.append(s1)

## time step
O.dt=.2*utils.PWaveTimeStep()
O.saveTmp('Mindlin')

from woo import qt
qt.View()
qt.Controller()

############################################
##### now the part pertaining to plots #####
############################################