~ubuntu-branches/debian/jessie/yade/jessie : revision 36

1

# encoding: utf-8

2

#

3

# utility functions for yade

4

#

5

# 2008-2009 © Václav Šmilauer <eudoxos@arcig.cz>

6

7

"""Heap of functions that don't (yet) fit anywhere else.

8

9

Devs: please DO NOT ADD more functions here, it is getting too crowded!

10

"""

11

12

import math,random,doctest,geom,numpy

13

from yade import *

14

from yade.wrapper import *

15

try: # use psyco if available

16

import psyco

17

psyco.full()

18

except ImportError: pass

19

20

try:

21

from minieigen import *

22

except ImportError:

23

from miniEigen import *

24

25

26

# c++ implementations for performance reasons

27

from yade._utils import *

28

29

def saveVars(mark='',loadNow=True,**kw):

30

"""Save passed variables into the simulation so that it can be recovered when the simulation is loaded again.

31

32

For example, variables *a*, *b* and *c* are defined. To save them, use::

33

34

>>> saveVars('something',a=1,b=2,c=3)

35

>>> from yade.params.something import *

36

>>> a,b,c

37

(1, 2, 3)

38

39

those variables will be save in the .xml file, when the simulation itself is saved. To recover those variables once the .xml is loaded again, use

40

41

>>> loadVars('something')

42

43

and they will be defined in the yade.params.\ *mark* module. The *loadNow* parameter calls :yref:`yade.utils.loadVars` after saving automatically.

44

45

If 'something' already exists, given variables will be inserted.

46

"""

47

import cPickle

48

try:

49

d=cPickle.loads(Omega().tags['pickledPythonVariablesDictionary'+mark]) #load dictionary d

50

for key in kw.keys():

51

d[key]=kw[key] #insert new variables into d

52

except KeyError:

53

d = kw

54

Omega().tags['pickledPythonVariablesDictionary'+mark]=cPickle.dumps(d)

55

if loadNow: loadVars(mark)

56

57

58

def loadVars(mark=None):

59

"""Load variables from :yref:`yade.utils.saveVars`, which are saved inside the simulation.

60

If ``mark==None``, all save variables are loaded. Otherwise only those with

61

the mark passed."""

62

import cPickle, types, sys, warnings

63

def loadOne(d,mark=None):

64

"""Load given dictionary into a synthesized module yade.params.name (or yade.params if *name* is not given). Update yade.params.__all__ as well."""

65

import yade.params

66

if mark:

67

if mark in yade.params.__dict__: warnings.warn('Overwriting yade.params.%s which already exists.'%mark)

68

modName='yade.params.'+mark

69

mod=types.ModuleType(modName)

70

mod.__dict__.update(d)

71

mod.__all__=list(d.keys()) # otherwise params starting with underscore would not be imported

72

sys.modules[modName]=mod

73

yade.params.__all__.append(mark)

74

yade.params.__dict__[mark]=mod

75

else:

76

yade.params.__all__+=list(d.keys())

77

yade.params.__dict__.update(d)

78

if mark!=None:

79

d=cPickle.loads(Omega().tags['pickledPythonVariablesDictionary'+mark])

80

loadOne(d,mark)

81

else: # load everything one by one

82

for m in Omega().tags.keys():

83

if m.startswith('pickledPythonVariablesDictionary'):

84

loadVars(m[len('pickledPythonVariableDictionary')+1:])

85

86

87

def SpherePWaveTimeStep(radius,density,young):

88

r"""Compute P-wave critical timestep for a single (presumably representative) sphere, using formula for P-Wave propagation speed $\Delta t_{c}=\frac{r}{\sqrt{E/\rho}}$.

89

If you want to compute minimum critical timestep for all spheres in the simulation, use :yref:`yade.utils.PWaveTimeStep` instead.

90

91

>>> SpherePWaveTimeStep(1e-3,2400,30e9)

92

2.8284271247461903e-07

93

"""

94

from math import sqrt

95

return radius/sqrt(young/density)

96

97

def randomColor():

98

"""Return random Vector3 with each component in interval 0…1 (uniform distribution)"""

99

return Vector3(random.random(),random.random(),random.random())

100

101

def typedEngine(name):

102

"""Return first engine from current O.engines, identified by its type (as string). For example:

103

104

>>> from yade import utils

105

>>> O.engines=[InsertionSortCollider(),NewtonIntegrator(),GravityEngine()]

106

>>> utils.typedEngine("NewtonIntegrator") == O.engines[1]

107

True

108

"""

109

return [e for e in Omega().engines if e.__class__.__name__==name][0]

110

111

def defaultMaterial():

112

"""Return default material, when creating bodies with :yref:`yade.utils.sphere` and friends, material is unspecified and there is no shared material defined yet. By default, this function returns::

113

114

.. code-block:: python

115

116

FrictMat(density=1e3,young=1e7,poisson=.3,frictionAngle=.5,label='defaultMat')

117

"""

118

return FrictMat(density=1e3,young=1e7,poisson=.3,frictionAngle=.5,label='defaultMat')

119

120

def _commonBodySetup(b,volume,geomInertia,material,pos,noBound=False,resetState=True,dynamic=None,fixed=False):

121

"""Assign common body parameters."""

122

if isinstance(material,int):

123

if material<0 and len(O.materials)==0: O.materials.append(defaultMaterial());

124

b.mat=O.materials[material]

125

elif isinstance(material,str): b.mat=O.materials[material]

126

elif isinstance(material,Material): b.mat=material

127

elif callable(material): b.mat=material()

128

else: raise TypeError("The 'material' argument must be None (for defaultMaterial), string (for shared material label), int (for shared material id) or Material instance.");

129

## resets state (!!)

130

if resetState: b.state=b.mat.newAssocState()

131

mass=volume*b.mat.density

132

b.state.mass,b.state.inertia=mass,geomInertia*b.mat.density

133

b.state.pos=b.state.refPos=pos

134

b.bounded=(not noBound)

135

if dynamic!=None:

136

import warnings

137

warnings.warn('dynamic=%s is deprecated, use fixed=%s instead'%(str(dynamic),str(not dynamic)),category=DeprecationWarning,stacklevel=2)

138

fixed=not dynamic

139

b.state.blockedDOFs=('xyzXYZ' if fixed else '')

140

141

142

def sphere(center,radius,dynamic=None,fixed=False,wire=False,color=None,highlight=False,material=-1,mask=1):

143

"""Create sphere with given parameters; mass and inertia computed automatically.

144

145

Last assigned material is used by default (*material* = -1), and utils.defaultMaterial() will be used if no material is defined at all.

146

147

:param Vector3 center: center

148

:param float radius: radius

149

:param float dynamic: deprecated, see "fixed"

150

:param float fixed: generate the body with all DOFs blocked?

151

:param material:

152

specify :yref:`Body.material`; different types are accepted:

153

* int: O.materials[material] will be used; as a special case, if material==-1 and there is no shared materials defined, utils.defaultMaterial() will be assigned to O.materials[0]

154

* string: label of an existing material that will be used

155

* :yref:`Material` instance: this instance will be used

156

* callable: will be called without arguments; returned Material value will be used (Material factory object, if you like)

157

:param int mask: :yref:`Body.mask` for the body

158

:param wire: display as wire sphere?

159

:param highlight: highlight this body in the viewer?

160

:param Vector3-or-None: body's color, as normalized RGB; random color will be assigned if ``None``.

161

162

:return:

163

A Body instance with desired characteristics.

164

165

166

Creating default shared material if none exists neither is given::

167

168

>>> O.reset()

169

>>> from yade import utils

170

>>> len(O.materials)

171

0

172

>>> s0=utils.sphere([2,0,0],1)

173

>>> len(O.materials)

174

1

175

176

Instance of material can be given::

177

178

>>> s1=utils.sphere([0,0,0],1,wire=False,color=(0,1,0),material=ElastMat(young=30e9,density=2e3))

179

>>> s1.shape.wire

180

False

181

>>> s1.shape.color

182

Vector3(0,1,0)

183

>>> s1.mat.density

184

2000.0

185

186

Material can be given by label::

187

188

>>> O.materials.append(FrictMat(young=10e9,poisson=.11,label='myMaterial'))

189

1

190

>>> s2=utils.sphere([0,0,2],1,material='myMaterial')

191

>>> s2.mat.label

192

'myMaterial'

193

>>> s2.mat.poisson

194

0.11

195

196

Finally, material can be a callable object (taking no arguments), which returns a Material instance.

197

Use this if you don't call this function directly (for instance, through yade.pack.randomDensePack), passing

198

only 1 *material* parameter, but you don't want material to be shared.

199

200

For instance, randomized material properties can be created like this:

201

202

>>> import random

203

>>> def matFactory(): return ElastMat(young=1e10*random.random(),density=1e3+1e3*random.random())

204

...

205

>>> s3=utils.sphere([0,2,0],1,material=matFactory)

206

>>> s4=utils.sphere([1,2,0],1,material=matFactory)

207

208

"""

209

b=Body()

210

b.shape=Sphere(radius=radius,color=color if color else randomColor(),wire=wire,highlight=highlight)

211

V=(4./3)*math.pi*radius**3

212

geomInert=(2./5.)*V*radius**2

213

_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,pos=center,dynamic=dynamic,fixed=fixed)

214

b.aspherical=False

215

b.mask=mask

216

return b

217

218

def box(center,extents,orientation=Quaternion(1,0,0,0),dynamic=None,fixed=False,wire=False,color=None,highlight=False,material=-1,mask=1):

219

"""Create box (cuboid) with given parameters.

220

221

:param Vector3 extents: half-sizes along x,y,z axes

222

223

See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""

224

b=Body()

225

b.shape=Box(extents=extents,color=color if color else randomColor(),wire=wire,highlight=highlight)

226

V=8*extents[0]*extents[1]*extents[2]

227

geomInert=Vector3(4*(extents[1]**2+extents[2]**2),4*(extents[0]**2+extents[2]**2),4*(extents[0]**2+extents[1]**2))

228

_commonBodySetup(b,V,geomInert,material,pos=center,dynamic=dynamic,fixed=fixed)

229

b.state.ori=orientation

230

b.mask=mask

231

b.aspherical=True

232

return b

233

234

235

def chainedCylinder(begin=Vector3(0,0,0),end=Vector3(1.,0.,0.),radius=0.2,dynamic=None,fixed=False,wire=False,color=None,highlight=False,material=-1,mask=1):

236

"""Create and connect a chainedCylinder with given parameters. The shape generated by repeted calls of this function is the Minkowski sum of polyline and sphere.

237

238

:param Real radius: radius of sphere in the Minkowski sum.

239

:param Vector3 begin: first point positioning the line in the Minkowski sum

240

:param Vector3 last: last point positioning the line in the Minkowski sum

241

242

In order to build a correct chain, last point of element of rank N must correspond to first point of element of rank N+1 in the same chain (with some tolerance, since bounding boxes will be used to create connections.

243

244

:return: Body object with the :yref:`ChainedCylinder` :yref:`shape<Body.shape>`.

245

"""

246

segment=end-begin

247

b=Body()

248

b.shape=ChainedCylinder(radius=radius,length=segment.norm(),color=color if color else randomColor(),wire=wire,highlight=highlight)

249

b.shape.segment=segment

250

V=2*(4./3)*math.pi*radius**3

251

geomInert=(2./5.)*V*radius**2+b.shape.length*b.shape.length*2*(4./3)*math.pi*radius**3

252

b.state=ChainedState()

253

b.state.addToChain(O.bodies.append(b))

254

_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,pos=begin,resetState=False,dynamic=dynamic,fixed=fixed)

255

b.mask=mask

256

b.bound=Aabb(color=[0,1,0])

257

b.state.ori.setFromTwoVectors(Vector3(0.,0.,1.),segment)

258

if (end == begin): b.state.ori = Quaternion((1,0,0),0)

259

return b

260

261

def gridNode(center,radius,dynamic=None,fixed=False,wire=False,color=None,highlight=False,material=-1):

262

b=Body()

263

b.shape=GridNode(radius=radius,color=color if color else randomColor(),wire=wire,highlight=highlight)

264

#V=(4./3)*math.pi*radius**3 # will be overriden by the connection

265

V=0.

266

geomInert=(2./5.)*V*radius**2 # will be overriden by the connection

267

_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,pos=center,dynamic=dynamic,fixed=fixed)

268

b.aspherical=False

269

b.bounded=False

270

b.mask=0 #avoid contact detection with the nodes. Manual interaction will be set for them in "gridConnection" below.

271

return b

272

273

def gridConnection(id1,id2,radius,wire=False,color=None,highlight=False,material=-1,mask=1,cellDist=None):

274

b=Body()

275

b.shape=GridConnection(radius=radius,color=color if color else randomColor(),wire=wire,highlight=highlight)

276

sph1=O.bodies[id1] ; sph2=O.bodies[id2]

277

i=createInteraction(id1,id2)

278

nodeMat=sph1.material

279

b.shape.node1=sph1 ; b.shape.node2=sph2

280

sph1.shape.addConnection(b) ; sph2.shape.addConnection(b)

281

if(O.periodic):

282

if(cellDist!=None):

283

i.cellDist=cellDist

284

segt=sph2.state.pos + O.cell.hSize*i.cellDist - sph1.state.pos

285

else: segt=sph2.state.pos - sph1.state.pos

286

L=segt.norm()

287

V=0.5*L*math.pi*radius**2

288

geomInert=(2./5.)*V*radius**2

289

_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,pos=sph1.state.pos,dynamic=False,fixed=True)

290

sph1.state.mass = sph1.state.mass + V*nodeMat.density

291

sph2.state.mass = sph2.state.mass + V*nodeMat.density

292

for k in [0,1,2]:

293

sph1.state.inertia[k] = sph1.state.inertia[k] + geomInert*nodeMat.density

294

sph2.state.inertia[k] = sph2.state.inertia[k] + geomInert*nodeMat.density

295

b.aspherical=False

296

if O.periodic:

297

i.phys.unp= -(sph2.state.pos + O.cell.hSize*i.cellDist - sph1.state.pos).norm() + sph1.shape.radius + sph2.shape.radius

298

b.shape.periodic=True

299

b.shape.cellDist=i.cellDist

300

else:

301

i.phys.unp= -(sph2.state.pos - sph1.state.pos).norm() + sph1.shape.radius + sph2.shape.radius

302

i.geom.connectionBody=b

303

I=math.pi*(2.*radius)**4/64.

304

E=nodeMat.young

305

i.phys.kn=E*math.pi*(radius**2)/L

306

i.phys.kr=E*I/L

307

i.phys.ks=12.*E*I/(L**3)

308

G=E/(2.*(1+nodeMat.poisson))

309

i.phys.ktw=2.*I*G/L

310

b.mask=mask

311

return b

312

313

314

def wall(position,axis,sense=0,color=None,material=-1,mask=1):

315

"""Return ready-made wall body.

316

317

:param float-or-Vector3 position: center of the wall. If float, it is the position along given axis, the other 2 components being zero

318

:param ∈{0,1,2} axis: orientation of the wall normal (0,1,2) for x,y,z (sc. planes yz, xz, xy)

319

:param ∈{-1,0,1} sense: sense in which to interact (0: both, -1: negative, +1: positive; see :yref:`Wall`)

320

321

See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""

322

b=Body()

323

b.shape=Wall(sense=sense,axis=axis,color=color if color else randomColor())

324

if isinstance(position,(int,long,float)):

325

pos2=Vector3(0,0,0); pos2[axis]=position

326

else: pos2=position

327

_commonBodySetup(b,0,Vector3(0,0,0),material,pos=pos2,fixed=True)

328

b.aspherical=False # wall never moves dynamically

329

b.mask=mask

330

return b

331

332

def facet(vertices,dynamic=None,fixed=True,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):

333

"""Create facet with given parameters.

334

335

:param [Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.

336

:param bool wire: if ``True``, facets are shown as skeleton; otherwise facets are filled

337

:param bool noBound: set :yref:`Body.bounded`

338

:param Vector3-or-None color: color of the facet; random color will be assigned if ``None``.

339

340

See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""

341

b=Body()

342

center=inscribedCircleCenter(vertices[0],vertices[1],vertices[2])

343

vertices=Vector3(vertices[0])-center,Vector3(vertices[1])-center,Vector3(vertices[2])-center

344

b.shape=Facet(color=color if color else randomColor(),wire=wire,highlight=highlight,vertices=vertices)

345

_commonBodySetup(b,0,Vector3(0,0,0),material,noBound=noBound,pos=center,fixed=fixed)

346

b.aspherical=False # mass and inertia are 0 anyway; fell free to change to ``True`` if needed

347

b.mask=mask

348

b.chain=chain

349

return b

350

351

def tetraPoly(vertices,strictCheck=True,dynamic=True,fixed=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):

352

"""Create tetrahedron (actually simple Polyhedra) with given parameters.

353

354

:param [Vector3,Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.

355

356

See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""

357

b=Body()

358

b.shape = Polyhedra(v=vertices,color=color if color else randomColor(),wire=wire,highlight=highlight)

359

volume = b.shape.GetVolume()

360

inertia = b.shape.GetInertia()

361

center = b.shape.GetCentroid()

362

_commonBodySetup(b,volume,inertia,material,noBound=noBound,pos=center,fixed=fixed)

363

b.aspherical=False # mass and inertia are 0 anyway; fell free to change to ``True`` if needed

364

b.state.ori = b.shape.GetOri()

365

b.mask=mask

366

b.chain=chain

367

return b

368

369

def tetra(vertices,strictCheck=True,dynamic=True,fixed=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):

370

"""Create tetrahedron with given parameters.

371

372

:param [Vector3,Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.

373

:param bool strictCheck: checks vertices order, raise RuntimeError for negative volume

374

375

See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""

376

b=Body()

377

center = .25*sum(vertices,Vector3.Zero)

378

volume = TetrahedronSignedVolume(vertices)

379

if volume < 0:

380

if strictCheck:

381

raise RuntimeError, "tetra: wrong order of vertices"

382

temp = vertices[3]

383

vertices[3] = vertices[2]

384

vertices[2] = temp

385

volume = TetrahedronSignedVolume(vertices)

386

assert(volume>0)

387

b.shape = Tetra(v=vertices,color=color if color else randomColor(),wire=wire,highlight=highlight)

388

# modifies pos, ori and inertia

389

ori = TetrahedronWithLocalAxesPrincipal(b)

390

_commonBodySetup(b,volume,b.state.inertia,material,noBound=noBound,pos=center,fixed=fixed)

391

b.state.ori = b.state.refOri = ori

392

b.aspherical = True

393

b.mask = mask

394

b.chain = chain

395

return b

396

397

398

399

400

401

402

#def setNewVerticesOfFacet(b,vertices):

403

# center = inscribedCircleCenter(vertices[0],vertices[1],vertices[2])

404

# vertices = Vector3(vertices[0])-center,Vector3(vertices[1])-center,Vector3(vertices[2])-center

405

# b.shape.vertices = vertices

406

# b.state.pos = center

407

408

409

def facetBox(*args,**kw):

410

"|ydeprecated|"

411

_deprecatedUtilsFunction('facetBox','geom.facetBox')

412

return geom.facetBox(*args,**kw)

413

414

def facetCylinder(*args,**kw):

415

"|ydeprecated|"

416

_deprecatedUtilsFunction('facetCylinder','geom.facetCylinder')

417

return geom.facetCylinder(*args,**kw)

418

419

def aabbWalls(extrema=None,thickness=0,oversizeFactor=1.5,**kw):

420

"""Return 6 boxes that will wrap existing packing as walls from all sides;

421

extrema are extremal points of the Aabb of the packing (will be calculated if not specified)

422

thickness is wall thickness (will be 1/10 of the X-dimension if not specified)

423

Walls will be enlarged in their plane by oversizeFactor.

424

returns list of 6 wall Bodies enclosing the packing, in the order minX,maxX,minY,maxY,minZ,maxZ.

425

"""

426

walls=[]

427

if not extrema: extrema=aabbExtrema()

428

#if not thickness: thickness=(extrema[1][0]-extrema[0][0])/10.

429

for axis in [0,1,2]:

430

mi,ma=extrema

431

center=[(mi[i]+ma[i])/2. for i in range(3)]

432

extents=[.5*oversizeFactor*(ma[i]-mi[i]) for i in range(3)]

433

extents[axis]=thickness/2.

434

for j in [0,1]:

435

center[axis]=extrema[j][axis]+(j-.5)*thickness

436

walls.append(box(center=center,extents=extents,fixed=True,**kw))

437

walls[-1].shape.wire=True

438

return walls

439

440

441

def aabbDim(cutoff=0.,centers=False):

442

"""Return dimensions of the axis-aligned bounding box, optionally with relative part *cutoff* cut away."""

443

a=aabbExtrema(cutoff,centers)

444

return (a[1][0]-a[0][0],a[1][1]-a[0][1],a[1][2]-a[0][2])

445

446

def aabbExtrema2d(pts):

447

"""Return 2d bounding box for a sequence of 2-tuples."""

448

inf=float('inf')

449

min,max=[inf,inf],[-inf,-inf]

450

for pt in pts:

451

if pt[0]<min[0]: min[0]=pt[0]

452

elif pt[0]>max[0]: max[0]=pt[0]

453

if pt[1]<min[1]: min[1]=pt[1]

454

elif pt[1]>max[1]: max[1]=pt[1]

455

return tuple(min),tuple(max)

456

457

def perpendicularArea(axis):

458

"""Return area perpendicular to given axis (0=x,1=y,2=z) generated by bodies

459

for which the function consider returns True (defaults to returning True always)

460

and which is of the type :yref:`Sphere`.

461

"""

462

ext=aabbExtrema()

463

other=((axis+1)%3,(axis+2)%3)

464

return (ext[1][other[0]]-ext[0][other[0]])*(ext[1][other[1]]-ext[0][other[1]])

465

466

def fractionalBox(fraction=1.,minMax=None):

467

"""retrurn (min,max) that is the original minMax box (or aabb of the whole simulation if not specified)

468

linearly scaled around its center to the fraction factor"""

469

if not minMax: minMax=aabbExtrema()

470

half=[.5*(minMax[1][i]-minMax[0][i]) for i in [0,1,2]]

471

return (tuple([minMax[0][i]+(1-fraction)*half[i] for i in [0,1,2]]),tuple([minMax[1][i]-(1-fraction)*half[i] for i in [0,1,2]]))

472

473

474

def randomizeColors(onlyDynamic=False):

475

"""Assign random colors to :yref:`Shape::color`.

476

477

If onlyDynamic is true, only dynamic bodies will have the color changed.

478

"""

479

for b in O.bodies:

480

color=(random.random(),random.random(),random.random())

481

if b.dynamic or not onlyDynamic: b.shape.color=color

482

483

def avgNumInteractions(cutoff=0.,skipFree=False,considerClumps=False):

484

r"""Return average numer of interactions per particle, also known as *coordination number* $Z$. This number is defined as

485

486

.. math:: Z=2C/N

487

488

where $C$ is number of contacts and $N$ is number of particles. When clumps are present, number of particles is the sum of standalone spheres plus the sum of clumps.

489

Clumps are considered in the calculation if cutoff != 0 or skipFree = True. If cutoff=0 (default) and skipFree=False (default) one needs to set considerClumps=True to consider clumps in the calculation.

490

491

With *skipFree*, particles not contributing to stable state of the packing are skipped, following equation (8) given in [Thornton2000]_:

492

493

.. math:: Z_m=\frac{2C-N_1}{N-N_0-N_1}

494

495

:param cutoff: cut some relative part of the sample's bounding box away.

496

:param skipFree: see above.

497

:param considerClumps: also consider clumps if cutoff=0 and skipFree=False; for further explanation see above.

498

499

"""

500

if cutoff==0 and not skipFree and not considerClumps: return 2*O.interactions.countReal()*1./len(O.bodies)

501

else:

502

nums,counts=bodyNumInteractionsHistogram(aabbExtrema(cutoff))

503

## CC is 2*C

504

CC=sum([nums[i]*counts[i] for i in range(len(nums))]); N=sum(counts)

505

if not skipFree: return CC*1./N if N>0 else float('nan')

506

## find bins with 0 and 1 spheres

507

N0=0 if (0 not in nums) else counts[nums.index(0)]

508

N1=0 if (1 not in nums) else counts[nums.index(1)]

509

NN=N-N0-N1

510

return (CC-N1)*1./NN if NN>0 else float('nan')

511

512

def plotNumInteractionsHistogram(cutoff=0.):

513

"Plot histogram with number of interactions per body, optionally cutting away *cutoff* relative axis-aligned box from specimen margin."

514

nums,counts=bodyNumInteractionsHistogram(aabbExtrema(cutoff))

515

import pylab

516

pylab.bar(nums,counts)

517

pylab.title('Number of interactions histogram, average %g (cutoff=%g)'%(avgNumInteractions(cutoff),cutoff))

518

pylab.xlabel('Number of interactions')

519

pylab.ylabel('Body count')

520

pylab.ion()

521

pylab.show()

522

523

def plotDirections(aabb=(),mask=0,bins=20,numHist=True,noShow=False):

524

"""Plot 3 histograms for distribution of interaction directions, in yz,xz and xy planes and

525

(optional but default) histogram of number of interactions per body.

526

527

:returns: If *noShow* is ``False``, displays the figure and returns nothing. If *noShow*, the figure object is returned without being displayed (works the same way as :yref:`yade.plot.plot`).

528

"""

529

import pylab,math

530

from yade import utils

531

for axis in [0,1,2]:

532

d=utils.interactionAnglesHistogram(axis,mask=mask,bins=bins,aabb=aabb)

533

fc=[0,0,0]; fc[axis]=1.

534

subp=pylab.subplot(220+axis+1,polar=True);

535

# 1.1 makes small gaps between values (but the column is a bit decentered)

536

pylab.bar(d[0],d[1],width=math.pi/(1.1*bins),fc=fc,alpha=.7,label=['yz','xz','xy'][axis])

537

#pylab.title(['yz','xz','xy'][axis]+' plane')

538

pylab.text(.5,.25,['yz','xz','xy'][axis],horizontalalignment='center',verticalalignment='center',transform=subp.transAxes,fontsize='xx-large')

539

if numHist:

540

pylab.subplot(224,polar=False)

541

nums,counts=utils.bodyNumInteractionsHistogram(aabb if len(aabb)>0 else utils.aabbExtrema())

542

avg=sum([nums[i]*counts[i] for i in range(len(nums))])/(1.*sum(counts))

543

pylab.bar(nums,counts,fc=[1,1,0],alpha=.7,align='center')

544

pylab.xlabel('Interactions per body (avg. %g)'%avg)

545

pylab.axvline(x=avg,linewidth=3,color='r')

546

pylab.ylabel('Body count')

547

if noShow: return pylab.gcf()

548

else:

549

pylab.ion()

550

pylab.show()

551

552

553

def encodeVideoFromFrames(*args,**kw):

554

"|ydeprecated|"

555

_deprecatedUtilsFunction('utils.encodeVideoFromFrames','utils.makeVideo')

556

return makeVideo(*args,**kw)

557

558

def makeVideo(frameSpec,out,renameNotOverwrite=True,fps=24,kbps=6000,bps=None):

559

"""Create a video from external image files using `mencoder <http://www.mplayerhq.hu>`__. Two-pass encoding using the default mencoder codec (mpeg4) is performed, running multi-threaded with number of threads equal to number of OpenMP threads allocated for Yade.

560

561

:param frameSpec: wildcard | sequence of filenames. If list or tuple, filenames to be encoded in given order; otherwise wildcard understood by mencoder's mf:// URI option (shell wildcards such as ``/tmp/snap-*.png`` or and printf-style pattern like ``/tmp/snap-%05d.png``)

562

:param str out: file to save video into

563

:param bool renameNotOverwrite: if True, existing same-named video file will have -*number* appended; will be overwritten otherwise.

564

:param int fps: Frames per second (``-mf fps=…``)

565

:param int kbps: Bitrate (``-lavcopts vbitrate=…``) in kb/s

566

"""

567

import os,os.path,subprocess,warnings

568

if bps!=None:

569

warnings.warn('plot.makeVideo: bps is deprecated, use kbps instead (the significance is the same, but the name is more precise)',stacklevel=2,category=DeprecationWarning)

570

kbps=bps

571

if renameNotOverwrite and os.path.exists(out):

572

i=0

573

while(os.path.exists(out+"~%d"%i)): i+=1

574

os.rename(out,out+"~%d"%i); print "Output file `%s' already existed, old file renamed to `%s'"%(out,out+"~%d"%i)

575

if isinstance(frameSpec,list) or isinstance(frameSpec,tuple): frameSpec=','.join(frameSpec)

576

for passNo in (1,2):

577

cmd=['mencoder','mf://%s'%frameSpec,'-mf','fps=%d'%int(fps),'-ovc','lavc','-lavcopts','vbitrate=%d:vpass=%d:threads=%d:%s'%(int(kbps),passNo,O.numThreads,'turbo' if passNo==1 else ''),'-o',('/dev/null' if passNo==1 else out)]

578

print 'Pass %d:'%passNo,' '.join(cmd)

579

ret=subprocess.call(cmd)

580

if ret!=0: raise RuntimeError("Error when running mencoder.")

581

582

def replaceCollider(colliderEngine):

583

"""Replaces collider (Collider) engine with the engine supplied. Raises error if no collider is in engines."""

584

colliderIdx=-1

585

for i,e in enumerate(O.engines):

586

if O.isChildClassOf(e.__class__.__name__,"Collider"):

587

colliderIdx=i

588

break

589

if colliderIdx<0: raise RuntimeError("No Collider found within O.engines.")

590

O.engines=O.engines[:colliderIdx]+[colliderEngine]+O.engines[colliderIdx+1:]

591

592

def _procStatus(name):

593

import os

594

for l in open('/proc/%d/status'%os.getpid()):

595

if l.split(':')[0]==name: return l

596

raise "No such line in /proc/[pid]/status: "+name

597

def vmData():

598

"Return memory usage data from Linux's /proc/[pid]/status, line VmData."

599

l=_procStatus('VmData'); ll=l.split(); assert(ll[2]=='kB')

600

return int(ll[1])

601

602

def uniaxialTestFeatures(filename=None,areaSections=10,axis=-1,distFactor=2.2,**kw):

603

"""Get some data about the current packing useful for uniaxial test:

604

605

#. Find the dimensions that is the longest (uniaxial loading axis)

606

607

#. Find the minimum cross-section area of the specimen by examining several (areaSections) sections perpendicular to axis, computing area of the convex hull for each one. This will work also for non-prismatic specimen.

608

609

#. Find the bodies that are on the negative/positive boundary, to which the straining condition should be applied.

610

611

:param filename: if given, spheres will be loaded from this file (ASCII format); if not, current simulation will be used.

612

:param float areaSection: number of section that will be used to estimate cross-section

613

:param ∈{0,1,2} axis: if given, force strained axis, rather than computing it from predominant length

614

:return: dictionary with keys ``negIds``, ``posIds``, ``axis``, ``area``.

615

616

.. warning::

617

The function :yref:`yade.utils.approxSectionArea` uses convex hull algorithm to find the area, but the implementation is reported to be *buggy* (bot works in some cases). Always check this number, or fix the convex hull algorithm (it is documented in the source, see :ysrc:`py/_utils.cpp`).

618

619

"""

620

if filename: ids=spheresFromFile(filename,**kw)

621

else: ids=[b.id for b in O.bodies]

622

mm,mx=aabbExtrema()

623

dim=aabbDim();

624

if axis<0: axis=list(dim).index(max(dim)) # list(dim) for compat with python 2.5 which didn't have index defined for tuples yet (appeared in 2.6 first)

625

assert(axis in (0,1,2))

626

import numpy

627

areas=[approxSectionArea(coord,axis) for coord in numpy.linspace(mm[axis],mx[axis],num=10)[1:-1]]

628

negIds,posIds=negPosExtremeIds(axis=axis,distFactor=distFactor)

629

return {'negIds':negIds,'posIds':posIds,'axis':axis,'area':min(areas)}

630

631

def voxelPorosityTriaxial(triax,resolution=200,offset=0):

632

"""

633

Calculate the porosity of a sample, given the TriaxialCompressionEngine.

634

635

A function :yref:`yade.utils.voxelPorosity` is invoked, with the volume of a box enclosed by TriaxialCompressionEngine walls.

636

The additional parameter offset allows using a smaller volume inside the box, where each side of the volume is at offset distance

637

from the walls. By this way it is possible to find a more precise porosity of the sample, since at walls' contact the porosity is usually reduced.

638

639

A recommended value of offset is bigger or equal to the average radius of spheres inside.

640

641

The value of resolution depends on size of spheres used. It can be calibrated by invoking voxelPorosityTriaxial with offset=0 and

642

comparing the result with TriaxialCompressionEngine.porosity. After calibration, the offset can be set to radius, or a bigger value, to get

643

the result.

644

645

:param triax: the TriaxialCompressionEngine handle

646

:param resolution: voxel grid resolution

647

:param offset: offset distance

648

:return: the porosity of the sample inside given volume

649

650

Example invocation::

651

652

from yade import utils

653

rAvg=0.03

654

TriaxialTest(numberOfGrains=200,radiusMean=rAvg).load()

655

O.dt=-1

656

O.run(1000)

657

O.engines[4].porosity

658

0.44007807740143889

659

utils.voxelPorosityTriaxial(O.engines[4],200,0)

660

0.44055412500000002

661

utils.voxelPorosityTriaxial(O.engines[4],200,rAvg)

662

0.36798199999999998

663

"""

664

p_bottom = O.bodies[triax.wall_bottom_id].state.se3[0]

665

p_top = O.bodies[triax.wall_top_id ].state.se3[0]

666

p_left = O.bodies[triax.wall_left_id ].state.se3[0]

667

p_right = O.bodies[triax.wall_right_id ].state.se3[0]

668

p_front = O.bodies[triax.wall_front_id ].state.se3[0]

669

p_back = O.bodies[triax.wall_back_id ].state.se3[0]

670

th = (triax.thickness)*0.5+offset

671

x_0 = p_left [0] + th

672

x_1 = p_right [0] - th

673

y_0 = p_bottom[1] + th

674

y_1 = p_top [1] - th

675

z_0 = p_back [2] + th

676

z_1 = p_front [2] - th

677

a=Vector3(x_0,y_0,z_0)

678

b=Vector3(x_1,y_1,z_1)

679

return voxelPorosity(resolution,a,b)

680

681

682

def trackPerfomance(updateTime=5):

683

"""

684

Track perfomance of a simulation. (Experimental)

685

Will create new thread to produce some plots.

686

Useful for track perfomance of long run simulations (in bath mode for example).

687

"""

688

689

def __track_perfomance(updateTime):

690

pid=os.getpid()

691

threadsCpu={}

692

lastTime,lastIter=-1,-1

693

while 1:

694

time.sleep(updateTime)

695

if not O.running:

696

lastTime,lastIter=-1,-1

697

continue

698

if lastTime==-1:

699

lastTime=time.time();lastIter=O.iter

700

plot.plots.update({'Iteration':('Perfomance',None,'Bodies','Interactions')})

701

continue

702

curTime=time.time();curIter=O.iter

703

perf=(curIter-lastIter)/(curTime-lastTime)

704

out=subprocess.Popen(['top','-bH','-n1', ''.join(['-p',str(pid)])],stdout=subprocess.PIPE).communicate()[0].splitlines()

705

for s in out[7:-1]:

706

w=s.split()

707

threadsCpu[w[0]]=float(w[8])

708

plot.addData(Iteration=curIter,Iter=curIter,Perfomance=perf,Bodies=len(O.bodies),Interactions=len(O.interactions),**threadsCpu)

709

plot.plots.update({'Iter':threadsCpu.keys()})

710

lastTime=time.time();lastIter=O.iter

711

712

thread.start_new_thread(__track_perfomance,(updateTime))

713

714

715

def NormalRestitution2DampingRate(en):

716

r"""Compute the normal damping rate as a function of the normal coefficient of restitution $e_n$. For $e_n\in\langle0,1\rangle$ damping rate equals

717

718

.. math:: -\frac{\log e_n}{\sqrt{e_n^2+\pi^2}}

719

720

"""

721

if en == 0.0: return 0.999999999

722

if en == 1.0: return 0.0

723

from math import sqrt,log,pi

724

ln_en = math.log(en)

725

return (-ln_en/math.sqrt((math.pow(ln_en,2) + math.pi*math.pi)))

726

727

def xMirror(half):

728

"""Mirror a sequence of 2d points around the x axis (changing sign on the y coord).

729

The sequence should start up and then it will wrap from y downwards (or vice versa).

730

If the last point's x coord is zero, it will not be duplicated."""

731

return list(half)+[(x,-y) for x,y in reversed(half[:-1] if half[-1][1]==0 else half)]

732

733

#############################

734

##### deprecated functions

735

736

737

def _deprecatedUtilsFunction(old,new):

738

"Wrapper for deprecated functions, example below."

739

import warnings

740

warnings.warn('Function utils.%s is deprecated, use %s instead.'%(old,new),stacklevel=2,category=DeprecationWarning)

741

742

# example of _deprecatedUtilsFunction usage:

743

#

744

# def import_mesh_geometry(*args,**kw):

745

# "|ydeprecated|"

746

# _deprecatedUtilsFunction('import_mesh_geometry','yade.import.gmsh')

747

# import yade.ymport

748

# return yade.ymport.stl(*args,**kw)

749

750

751

class TableParamReader():

752

"""Class for reading simulation parameters from text file.

753

754

Each parameter is represented by one column, each parameter set by one line. Colums are separated by blanks (no quoting).

755

756

First non-empty line contains column titles (without quotes).

757

You may use special column named 'description' to describe this parameter set;

758

if such colum is absent, description will be built by concatenating column names and corresponding values (``param1=34,param2=12.22,param4=foo``)

759

760

* from columns ending in ``!`` (the ``!`` is not included in the column name)

761

* from all columns, if no columns end in ``!``.

762

763

Empty lines within the file are ignored (although counted); ``#`` starts comment till the end of line. Number of blank-separated columns must be the same for all non-empty lines.

764

765

A special value ``=`` can be used instead of parameter value; value from the previous non-empty line will be used instead (works recursively).

766

767

This class is used by :yref:`yade.utils.readParamsFromTable`.

768

"""

769

def __init__(self,file):

770

"Setup the reader class, read data into memory."

771

import re

772

# read file in memory, remove newlines and comments; the [''] makes lines 1-indexed

773

ll=[re.sub('\s*#.*','',l[:-1]) for l in ['']+open(file,'r').readlines()]

774

# usable lines are those that contain something else than just spaces

775

usableLines=[i for i in range(len(ll)) if not re.match(r'^\s*(#.*)?$',ll[i])]

776

headings=ll[usableLines[0]].split()

777

# use all values of which heading has ! after its name to build up the description string

778

# if there are none, use all columns

779

if not 'description' in headings:

780

bangHeads=[h[:-1] for h in headings if h[-1]=='!'] or headings

781

headings=[(h[:-1] if h[-1]=='!' else h) for h in headings]

782

usableLines=usableLines[1:] # and remove headinds from usableLines

783

values={}

784

for l in usableLines:

785

val={}

786

for i in range(len(headings)):

787

val[headings[i]]=ll[l].split()[i]

788

values[l]=val

789

lines=values.keys(); lines.sort()

790

# replace '=' by previous value of the parameter

791

for i,l in enumerate(lines):

792

for j in values[l].keys():

793

if values[l][j]=='=':

794

try:

795

values[l][j]=values[lines[i-1]][j]

796

except IndexError,KeyError:

797

raise RuntimeError("The = specifier on line %d refers to nonexistent value on previous line?"%l)

798

#import pprint; pprint.pprint(headings); pprint.pprint(values)

799

# add descriptions, but if they repeat, append line number as well

800

if not 'description' in headings:

801

descs=set()

802

for l in lines:

803

dd=','.join(head.replace('!','')+'='+('%g'%values[head] if isinstance(values[l][head],float) else str(values[l][head])) for head in bangHeads).replace("'",'').replace('"','')

804

if dd in descs: dd+='__line=%d__'%l

805

values[l]['description']=dd

806

descs.add(dd)

807

self.values=values

808

809

def paramDict(self):

810

"""Return dictionary containing data from file given to constructor. Keys are line numbers (which might be non-contiguous and refer to real line numbers that one can see in text editors), values are dictionaries mapping parameter names to their values given in the file. The special value '=' has already been interpreted, ``!`` (bangs) (if any) were already removed from column titles, ``description`` column has already been added (if absent)."""

811

return self.values

812

813

if __name__=="__main__":

814

tryTable="""head1 important2! !OMP_NUM_THREADS! abcd

815

1 1.1 1.2 1.3

816

'a' 'b' 'c' 'd' ### comment

817

818

# empty line

819

1 = = g

820

"""

821

file='/tmp/try-tbl.txt'

822

f=open(file,'w')

823

f.write(tryTable)

824

f.close()

825

from pprint import *

826

pprint(TableParamReader(file).paramDict())

827

828

def runningInBatch():

829

'Tell whether we are running inside the batch or separately.'

830

import os

831

return 'YADE_BATCH' in os.environ

832

833

def waitIfBatch():

834

'Block the simulation if running inside a batch. Typically used at the end of script so that it does not finish prematurely in batch mode (the execution would be ended in such a case).'

835

if runningInBatch(): O.wait()

836

837

def readParamsFromTable(tableFileLine=None,noTableOk=True,unknownOk=False,**kw):

838

"""

839

Read parameters from a file and assign them to __builtin__ variables.

840

841

The format of the file is as follows (commens starting with # and empty lines allowed)::

842

843

# commented lines allowed anywhere

844

name1 name2 … # first non-blank line are column headings

845

# empty line is OK, with or without comment

846

val1 val2 … # 1st parameter set

847

val2 val2 … # 2nd

848

…

849

850

Assigned tags (the ``description`` column is synthesized if absent,see :yref:`yade.utils.TableParamReader`);

851

852

O.tags['description']=… # assigns the description column; might be synthesized

853

O.tags['params']="name1=val1,name2=val2,…" # all explicitly assigned parameters

854

O.tags['defaultParams']="unassignedName1=defaultValue1,…" # parameters that were left at their defaults

855

O.tags['d.id']=O.tags['id']+'.'+O.tags['description']

856

O.tags['id.d']=O.tags['description']+'.'+O.tags['id']

857

858

All parameters (default as well as settable) are saved using :yref:`yade.utils.saveVars`\ ``('table')``.

859

860

:param tableFile: text file (with one value per blank-separated columns)

861

:param int tableLine: number of line where to get the values from

862

:param bool noTableOk: if False, raise exception if the file cannot be open; use default values otherwise

863

:param bool unknownOk: do not raise exception if unknown column name is found in the file, and assign it as well

864

:return: number of assigned parameters

865

"""

866

tagsParams=[]

867

# dictParams is what eventually ends up in yade.params.table (default+specified values)

868

dictDefaults,dictParams,dictAssign={},{},{}

869

import os, __builtin__,re,math

870

if not tableFileLine and ('YADE_BATCH' not in os.environ or os.environ['YADE_BATCH']==''):

871

if not noTableOk: raise EnvironmentError("YADE_BATCH is not defined in the environment")

872

O.tags['line']='l!'

873

else:

874

if not tableFileLine: tableFileLine=os.environ['YADE_BATCH']

875

env=tableFileLine.split(':')

876

tableFile,tableLine=env[0],int(env[1])

877

allTab=TableParamReader(tableFile).paramDict()

878

if not allTab.has_key(tableLine): raise RuntimeError("Table %s doesn't contain valid line number %d"%(tableFile,tableLine))

879

vv=allTab[tableLine]

880

O.tags['line']='l%d'%tableLine

881

O.tags['description']=vv['description']

882

O.tags['id.d']=O.tags['id']+'.'+O.tags['description']; O.tags['d.id']=O.tags['description']+'.'+O.tags['id']

883

# assign values specified in the table to python vars

884

# !something cols are skipped, those are env vars we don't treat at all (they are contained in description, though)

885

for col in vv.keys():

886

if col=='description' or col[0]=='!': continue

887

if col not in kw.keys() and (not unknownOk): raise NameError("Parameter `%s' has no default value assigned"%col)

888

if vv[col]=='*': vv[col]=kw[col] # use default value for * in the table

889

elif col in kw.keys(): kw.pop(col) # remove the var from kw, so that it contains only those that were default at the end of this loop

890

#print 'ASSIGN',col,vv[col]

891

tagsParams+=['%s=%s'%(col,vv[col])];

892

dictParams[col]=eval(vv[col],math.__dict__)

893

# assign remaining (default) keys to python vars

894

defaults=[]

895

for k in kw.keys():

896

dictDefaults[k]=kw[k]

897

defaults+=["%s=%s"%(k,kw[k])];

898

O.tags['defaultParams']=",".join(defaults)

899

O.tags['params']=",".join(tagsParams)

900

dictParams.update(dictDefaults)

901

saveVars('table',loadNow=True,**dictParams)

902

return len(tagsParams)

903

904

def psd(bins=5, mass=True, mask=-1):

905

"""Calculates particle size distribution.

906

907

:param int bins: number of bins

908

:param bool mass: if true, the mass-PSD will be calculated

909

:param int mask: :yref:`Body.mask` for the body

910

:return:

911

* binsSizes: list of bin's sizes

912

* binsProc: how much material (in percents) are in the bin, cumulative

913

* binsSumCum: how much material (in units) are in the bin, cumulative

914

915

binsSizes, binsProc, binsSumCum

916

"""

917

maxD = 0.0

918

minD = 0.0

919

920

for b in O.bodies:

921

if (isinstance(b.shape,Sphere) and ((mask<0) or ((b.mask&mask)<>0))):

922

if ((2*b.shape.radius) > maxD) : maxD = 2*b.shape.radius

923

if (((2*b.shape.radius) < minD) or (minD==0.0)): minD = 2*b.shape.radius

924

925

if (minD==maxD): return false #All particles are having the same size

926

927

binsSizes = numpy.linspace(minD, maxD, bins+1)

928

929

deltaBinD = (maxD-minD)/bins

930

binsMass = numpy.zeros(bins)

931

binsNumbers = numpy.zeros(bins)

932

933

for b in O.bodies:

934

if (isinstance(b.shape,Sphere) and ((mask<0) or ((b.mask&mask)<>0))):

935

d=2*b.shape.radius

936

937

basketId = int(math.floor( (d-minD) / deltaBinD ) )

938

if (d == maxD): basketId = bins-1 #If the diameter equals the maximal diameter, put the particle into the last bin

939

binsMass[basketId] = binsMass[basketId] + b.state.mass #Put masses into the bin

940

binsNumbers[basketId] = binsNumbers[basketId] + 1 #Put numbers into the bin

941

942

binsProc = numpy.zeros(bins+1)

943

binsSumCum = numpy.zeros(bins+1)

944

945

i=1

946

for size in binsSizes[:-1]:

947

if (mass): binsSumCum[i]+=binsSumCum[i-1]+binsMass[i-1] #Calculate mass

948

else: binsSumCum[i]+=binsSumCum[i-1]+binsNumbers[i-1] #Calculate number of particles

949

i+=1

950

951

if (binsSumCum[len(binsSumCum)-1] > 0):

952

i=0

953

for l in binsSumCum:

954

binsProc[i] = binsSumCum[i]/binsSumCum[len(binsSumCum)-1]

955

i+=1

956

return binsSizes, binsProc, binsSumCum

957

958

class clumpTemplate:

959

"""Create a clump template by a list of relative radii and a list of relative positions. Both lists must have the same length.

960

961

:param [float,float,...] relRadii: list of relative radii (minimum length = 2)

962

:param [Vector3,Vector3,...] relPositions: list of relative positions (minimum length = 2)

963

964

"""

965

def __init__(self,relRadii=[],relPositions=[[],[]]):

966

if (len(relRadii) != len(relPositions)):

967

raise ValueError("Given lists must have same length! Given lists does not match clump template structure.");

968

if (len(relRadii) < 2):

969

raise ValueError("One or more of given lists for relative radii have length < 2! Given lists does not match clump template structure.");

970

for ii in range(0,len(relPositions)):

971

if len(relPositions[ii]) != 3:

972

raise ValueError("One or more of given lists for relative positions do not have length of 3! Given lists does not match clump template structure.");

973

for jj in range(ii,len(relPositions)):

974

if ii != jj:

975

if (relPositions[ii] == relPositions[jj]):

976

raise ValueError("Two or more of given lists for relative positions are equal! Given lists does not match clump template structure.");

977

self.numCM = len(relRadii)

978

self.relRadii = relRadii

979

self.relPositions = relPositions

980

981

class UnstructuredGrid:

982

"""EXPERIMENTAL.

983

Class representing triangulated FEM-like unstructured grid. It is used for transfereing data from ad to YADE and external FEM program. The main purpose of this class is to store information about individual grid vertices/nodes coords (since facets stores only coordinates of vertices in local coords) and to avaluate and/or apply nodal forces from contact forces (from actual contact force and contact point the force is distributed to nodes using linear approximation).

984

985

TODO rewrite to C++

986

TODO better docs

987

988

:param dict vertices: dict of {internal vertex label:vertex}, e.g. {5:(0,0,0),22:(0,1,0),23:(1,0,0)}

989

:param dict connectivityTable: dict of {internal element label:[indices of vertices]}, e.g. {88:[5,22,23]}

990

"""

991

def __init__(self,vertices={},connectivityTable={},**kw):

992

self.vertices = vertices

993

self.connectivityTable = connectivityTable

994

self.forces = dict( (i,Vector3(0,0,0)) for i in vertices)

995

self.build(**kw)

996

def setup(self,vertices,connectivityTable,toSimulation=False,bodies=None,**kw):

997

"""Sets new information to receiver

998

999

:param dict vertices: see constructor for explanation

1000

:param dict connectivityTable: see constructor for explanation

1001

:param bool toSimulation: if new information should be inserted to Yade simulation (create new bodies or not)

1002

:param [[int]]|None bodies: list of list of bodies indices to be appended as clumps (thus no contact detection is done within one body)

1003

"""

1004

self.vertices = dict( (i,v) for i,v in vertices.iteritems())

1005

self.connectivityTable = dict( (i,e) for i,e in connectivityTable.iteritems())

1006

self.build(**kw)

1007

if toSimulation:

1008

self.toSimulation(bodies)

1009

def build(self,**kw):

1010

self.elements = {}

1011

for i,c in self.connectivityTable.iteritems():

1012

if len(c) == 3:

1013

b = facet([self.vertices[j] for j in c],**kw)

1014

elif len(c) == 4:

1015

b = tetra([self.vertices[j] for j in c],**kw)

1016

else:

1017

raise RuntimeError, "Unsupported cell shape (should be triangle or tetrahedron)"

1018

self.elements[i] = b

1019

def resetForces(self):

1020

for i in self.vertices:

1021

self.forces[i] = Vector3(0,0,0)

1022

def getForcesOfNodes(self):

1023

"""Computes forces for each vertex/node. The nodal force is computed from contact force and contact point using linear approximation

1024

"""

1025

self.resetForces()

1026

for i,e in self.elements.iteritems():

1027

ie = self.connectivityTable[i]

1028

for i in e.intrs():

1029

fn = i.phys.normalForce

1030

fs = i.phys.shearForce if hasattr(i.phys,"shearForce") else Vector3.Zero

1031

f = (fn+fs) * (-1 if e.id == i.id1 else +1 if e.id == i.id2 else 'ERROR')

1032

cp = i.geom.contactPoint

1033

if isinstance(e.shape,Facet):

1034

v0,v1,v2 = [Vector3(self.vertices[j]) for j in ie]

1035

w0 = ((cp-v1).cross(v2-v1)).norm()

1036

w1 = ((cp-v2).cross(v0-v2)).norm()

1037

w2 = ((cp-v0).cross(v1-v0)).norm()

1038

ww = w0+w1+w2

1039

self.forces[ie[0]] += f*w0/ww

1040

self.forces[ie[1]] += f*w1/ww

1041

self.forces[ie[2]] += f*w2/ww

1042

elif isinstance(e.shape,Tetra):

1043

v0,v1,v2,v3 = [Vector3(self.vertices[j]) for j in ie]

1044

w0 = TetrahedronVolume((cp,v1,v2,v3))

1045

w1 = TetrahedronVolume((cp,v2,v3,v0))

1046

w2 = TetrahedronVolume((cp,v3,v0,v1))

1047

w3 = TetrahedronVolume((cp,v0,v1,v2))

1048

ww = w0+w1+w2+w3

1049

self.forces[ie[0]] += f*w0/ww

1050

self.forces[ie[1]] += f*w1/ww

1051

self.forces[ie[2]] += f*w2/ww

1052

self.forces[ie[3]] += f*w3/ww

1053

return self.forces

1054

def setPositionsOfNodes(self,newPoss):

1055

"""Sets new position of nodes and also updates all elements in the simulation

1056

1057

:param [Vector3] newPoss: list of new positions

1058

"""

1059

for i in self.vertices:

1060

self.vertices[i] = newPoss[i]

1061

self.updateElements()

1062

def updateElements(self):

1063

"""Updates positions of all elements in the simulation

1064

"""

1065

for i,c in self.connectivityTable.iteritems():

1066

e = self.elements[i]

1067

e.state.pos = Vector3(0,0,0)

1068

e.state.ori = Quaternion((1,0,0),0)

1069

if isinstance(e.shape,Facet):

1070

e.shape.vertices = [self.vertices[j] for j in c]

1071

elif isinstance(e.shape,Tetra):

1072

e.shape.v = [self.vertices[j] for j in c]

1073

def toSimulation(self,bodies=None):

1074

"""Insert all elements to Yade simulation

1075

"""

1076

if bodies:

1077

e = self.elements.values()

1078

for b in bodies:

1079

O.bodies.appendClumped([e[i] for i in b])

1080

else:

1081

O.bodies.append(self.elements.values())