1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
import random
N = 30
def filterPos(X):
# the first condition create a concavity (for fun), the other crops to extract a cylinder from a cube
return (X - Vector3(N / 2, N / 2, N)).norm() < N / 2 or (X - Vector3(N / 2, X[1], N / 2)).norm() > N / 2
random.seed(2)
for i in range(N):
for j in range(N):
for k in range(N):
X = Vector3(i + 0.01 * random.random(), j + 0.01 * random.random(), k + 0.01 * random.random())
if not filterPos(X):
O.bodies.append(sphere(X, 0.49, color=(0.7, 0.7, 0.1)))
boxes = O.bodies.append([box((N / 2, -0.51, N / 2), (N, 0, N), fixed=True), box((N / 2, N - 0.49, N / 2), (N, 0, N), fixed=True)])
O.bodies[boxes[1]].state.vel = (0, -0.1, 0)
TW = TesselationWrapper()
TW.triangulate()
TW.addBoundingPlane(1, True)
TW.addBoundingPlane(1, False)
a = 6
shrinkedA = (sqrt(a) - 0.5)**2
fixedA = True
ag = TW.getAlphaGraph(alpha=a, shrinkedAlpha=shrinkedA, fixedAlpha=fixedA)
graph = GlExtra_AlphaGraph(tesselationWrapper=TW, wire=True)
from yade import qt
rr = qt.Renderer()
rr.extraDrawers = [graph]
TW.applyAlphaForces(stress=-Matrix3.Identity, alpha=a, shrinkedAlpha=shrinkedA, fixedAlpha=fixedA, reset=False)
def updateMembraneForces():
TW.triangulate(reset=True)
TW.addBoundingPlane(1, True)
TW.addBoundingPlane(1, False)
TW.applyAlphaForces(stress=-1000 * Matrix3.Identity, alpha=a, shrinkedAlpha=shrinkedA, fixedAlpha=fixedA, reset=False)
def displayMembrane():
#FIXME: suboptimal, we can't re-use the existing triangulation because it would try to insert fictious spheres again
# as a consequence we rebuild the whole thing just for display :-/
TW.triangulate(reset=True)
TW.addBoundingPlane(1, True)
TW.addBoundingPlane(1, False)
ag = TW.getAlphaGraph(alpha=a, shrinkedAlpha=shrinkedA, fixedAlpha=fixedA)
graph.refresh()
for b in O.bodies:
if isinstance(b.shape, Sphere):
if O.forces.permF(b.id) != Vector3.Zero:
b.shape.color = (0.8, 0.4, 0.4)
else:
b.shape.color = (0.4, 0.7, 0.4)
else:
b.shape.color = (1, 0, 0)
v = qt.View()
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()], verletDist=0.3),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8),
#triax,
NewtonIntegrator(damping=0.3),
PyRunner(iterPeriod=50, initRun=True, command="updateMembraneForces()", label="membrane"),
PyRunner(iterPeriod=100, initRun=True, command="displayMembrane()", label="GuiMembrane")
]
O.timingEnabled = True
O.step()
# make Video
#O.engines = O.engines + [qt.SnapshotEngine(fileBase=O.tmpFilename(), iterPeriod=100, label='snapshotter')]
#def done():
#makeVideo(snapshotter.snapshots,out='alphaWithWalls.mpeg', fps=20, kbps=5000)
|
