'''
Various computations affected by the periodic boundary conditions.
'''
import unittest
import random
from minieigen import *
from woo import utils
from woo.dem import *
from woo.core import *
from math import *
import woo
[docs]class TestSimpleClump(unittest.TestCase):
"Test things on a simple clump composed of 2 spheres."
[docs] 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
#print 100*'#'
#print self.b1,self.b1.dem.master,id(self.b1.dem.master),self.b1.dem.master._cxxAddr
#print self.b2,self.b2.dem.master,id(self.b2.dem.master),self.b1.dem.master._cxxAddr
#print self.bC,id(self.bC),self.bC._cxxAddr
[docs] def testConsistency(self):
"Clump: ids and flags consistency"
S=woo.master.scene
b1,b2,bC=self.b1,self.b2,self.bC
#self.assertEqual(b1.clumpId,bC.id)
#self.assertEqual(b2.clumpId,bC.id)
#self.assertEqual(bC.clumpId,bC.id)
self.assertTrue(b1 in bC.dem.nodes)
self.assertTrue(b2 in bC.dem.nodes)
self.assertTrue(bC.dem.clump)
self.assertTrue(b1.dem.clumped)
self.assertTrue(b2.dem.clumped)
self.assertTrue(b1.dem.master==bC)
self.assertTrue(b2.dem.master==bC)
[docs] def testStaticProperties(self):
"Clump: mass, centroid, intertia"
S=woo.master.scene
b1,b2,bC=self.b1,self.b2,self.bC
# mass
self.assertEqual(bC.dem.mass,b1.dem.mass+b2.dem.mass)
# centroid
SS=b1.dem.mass*b1.pos+b2.dem.mass*b2.pos
c=SS/bC.dem.mass
self.assertEqual(bC.pos,c);
# inertia
i1,i2=(8./15)*pi*S.dem.par[0].mat.density*S.dem.par[0].shape.radius**5, (8./15)*pi*S.dem.par[1].mat.density*S.dem.par[1].shape.radius**5 # inertia of spheres
iMax=i1+i2+b1.dem.mass*(b1.pos-c).norm()**2+b2.dem.mass*(b2.pos-c).norm()**2 # minimum principal inertia
iMin=i1+i2 # perpendicular to the
# the order of bC.state.inertia is arbitrary (though must match the orientation)
iC=list(bC.dem.inertia); iC.sort()
self.assertAlmostEqual(iC[0],iMin)
self.assertAlmostEqual(iC[1],iMax)
self.assertAlmostEqual(iC[2],iMax)
# check orientation...?
#self.assertAlmostEqual
[docs] def testVelocity(self):
"Clump: velocities of member assigned by Leapfrog"
S=woo.master.scene
b1,b2,bC=self.b1,self.b2,self.bC
S.dt=0
#print bC.dem.vel,bC.dem.angVel
bC.dem.vel=(1.,.2,.4)
bC.dem.angVel=(0,.4,.1)
self.assertTrue(self.b1.dem.master==self.bC)
S.engines=[Leapfrog(reset=True)]; S.one() # update velocities
# linear velocities
self.assertEqual(b1.dem.vel,bC.dem.vel+bC.dem.angVel.cross(b1.pos-bC.pos))
self.assertEqual(b2.dem.vel,bC.dem.vel+bC.dem.angVel.cross(b2.pos-bC.pos))
# angular velocities
self.assertEqual(b1.dem.angVel,bC.dem.angVel);
self.assertEqual(b2.dem.angVel,bC.dem.angVel);
[docs] 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.assertTrue(not woo.dem.Collider.mayCollide(S.dem,S.dem.par[i],S.dem.par[j]))
self.assertTrue(len(S.dem.con)==0)
[docs]def sphereClumpPrincipalAxes(cc,rr):
'Clump: return vol,pos,ori,inertia of sphere clump defined by centers and radii of spheres'
ii=range(len(rr))
Ii=[Matrix3(Vector3.Ones*((8/15.)*pi*rr[i]**5)) for i in ii]
Vi=[(4/3.)*pi*rr[i]**3 for i in ii]
V=sum(Vi)
Sg=sum([Vi[i]*cc[i] for i in ii],Vector3.Zero)
Ig=sum([Ii[i]+Vi[i]*(cc[i].dot(cc[i])*Matrix3.Identity-cc[i].outer(cc[i])) for i in ii],Matrix3.Zero)
pos=Sg/V
import numpy, numpy.linalg
Ic_orientG=Ig-V*(pos.dot(pos)*Matrix3.Identity-pos.outer(pos))
eigval,eigvec=numpy.linalg.eigh(numpy.array(Ic_orientG))
ori=Quaternion(Matrix3(*eigvec.tolist()))
inertia=Vector3(tuple(eigval))
return V,pos,ori,inertia
[docs]class TestSphereClumpGeom(unittest.TestCase):
"Clump: geometry of clumps composed of spheres only"
[docs] def setUp(self):
# c has no self-intersections
self.c=SphereClumpGeom(centers=[(0,0,0),(0,0,3)],radii=(1,.5),div=-1)
# c2 has several identical (overlapping) spheres
self.c2=SphereClumpGeom(centers=[(0,0,0),(0,0,0)],radii=(1,1),div=-1)
[docs] def testSteiner(self):
"SphereClumpGeom: properties via Steiner's theorem"
vol,pos,ori,inertia=sphereClumpPrincipalAxes(self.c.centers,self.c.radii)
self.c.recompute(div=0)
self.assertAlmostEqual(min(inertia),min(self.c.inertia))
self.assertAlmostEqual(max(inertia),max(self.c.inertia))
for i in (0,1,2): self.assertAlmostEqual(pos[i],self.c.pos[i])
self.assertAlmostEqual(vol,self.c.volume)
[docs] def testSampling(self):
"SphereClumpGeom: properties via grid sampling"
self.c2.recompute(div=10)
exactV=(4/3.)*pi*self.c.radii[0]**3
exactI=(8/15.)*pi*self.c.radii[0]**5
self.assertAlmostEqual(self.c2.volume,exactV,delta=exactV*1e-2)
self.assertAlmostEqual(self.c2.inertia[0],exactI,delta=exactI*2e-2)
#print 'volume',sum([(4/3.)*pi*self.c.radii[i]**3 for i in range(len(self.c.radii))]),self.c.volume
#self.assertAlmostEqual()