← Back to branch summary

~ubuntu-branches/ubuntu/trusty/yade/trusty

« back to all changes in this revision

Viewing changes to examples/concrete/triax.py

  • Committer: Package Import Robot
  • Author(s): Anton Gladky, cf3f8d9
  • Date: 2013-10-30 20:56:33 UTC
  • mfrom: (20.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20131030205633-1f01r7hjce17d723
Tags: 1.05.0-2
http://bugs.debian.org/726009
[cf3f8d9] Pass -ftrack-macro-expansion=0 only if gcc>=4.8. (Closes: #726009)

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples/concrete/triax.py
 
1
################################################################################
 
2
#
 
3
# Triaxial test. Axial strain rate is prescribed and transverse prestress.
 
4
# Test is possible on prism or cylinder
 
5
# An independent c++ engine may be created from this script in the future.
 
6
#
 
7
################################################################################
 
8
from yade import pack, plot
 
9
import os
 
10
 
 
11
# default parameters or from table
 
12
readParamsFromTable(noTableOk=True,
 
13
        # type of test ['cyl','cube']
 
14
        testType = 'cyl',
 
15
 
 
16
        # material parameters
 
17
        young = 20e9,
 
18
        poisson = .2,
 
19
        frictionAngle = 1.2,
 
20
        sigmaT = 1.5e6,
 
21
        epsCrackOnset = 1e-4,
 
22
        relDuctility = 30,
 
23
 
 
24
        # prestress
 
25
        preStress = -3e6,
 
26
        # axial strain rate
 
27
        strainRate = -20,
 
28
 
 
29
        # assamlby parameters
 
30
        rParticle = .075e-3, #
 
31
        width = 2e-3,
 
32
        height = 5e-3,
 
33
        bcCoeff = 5,
 
34
 
 
35
        # facets division
 
36
        nw = 24,
 
37
        nh = 15,
 
38
 
 
39
        # output specifications
 
40
        fileName = 'test',
 
41
        exportDir = '/tmp',
 
42
        runGnuplot = False,
 
43
        runInGui = True,
 
44
)
 
45
from yade.params.table import *
 
46
assert testType in ['cyl','cube']
 
47
 
 
48
# materials
 
49
concMat = O.materials.append(CpmMat(
 
50
        young=young,frictionAngle=frictionAngle,poisson=poisson,sigmaT=sigmaT,
 
51
        epsCrackOnset=epsCrackOnset,relDuctility=relDuctility
 
52
))
 
53
frictMat = O.materials.append(FrictMat(
 
54
        young=young,poisson=poisson,frictionAngle=frictionAngle
 
55
))
 
56
 
 
57
# spheres
 
58
pred = pack.inCylinder((0,0,0),(0,0,height),.5*width) if testType=='cyl' else pack.inAlignedBox((-.5*width,-.5*width,0),(.5*width,.5*width,height)) if testType=='cube' else None
 
59
sp=SpherePack()
 
60
sp = pack.randomDensePack(pred,spheresInCell=2000,radius=rParticle,memoizeDb='/tmp/triaxTestOnCylinder.sqlite',material=concMat,returnSpherePack=True)
 
61
spheres=sp.toSimulation(color=(0,1,1))
 
62
# bottom and top of specimen. Will have prescribed velocity
 
63
bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
 
64
top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
 
65
vel = strainRate*(height-rParticle*2*bcCoeff)
 
66
for s in bot:
 
67
        s.shape.color = (1,0,0)
 
68
        s.state.blockedDOFs = 'xyzXYZ'
 
69
        s.state.vel = (0,0,-vel)
 
70
for s in top:
 
71
        s.shape.color = Vector3(0,1,0)
 
72
        s.state.blockedDOFs = 'xyzXYZ'
 
73
        s.state.vel = (0,0,vel)
 
74
 
 
75
# facets
 
76
facets = []
 
77
if testType == 'cyl':
 
78
        rCyl2 = .5*width / cos(pi/float(nw))
 
79
        for r in xrange(nw):
 
80
                for h in xrange(nh):
 
81
                        v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
 
82
                        v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
 
83
                        v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
 
84
                        v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
 
85
                        f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
 
86
                        f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
 
87
                        facets.extend((f1,f2))
 
88
elif testType == 'cube':
 
89
        nw2 = nw/4
 
90
        for r in xrange(nw2):
 
91
                for h in xrange(nh):
 
92
                        v11 = Vector3( -.5*width + (r+0)*width/nw2, -.5*width, height*(h+0)/float(nh) )
 
93
                        v12 = Vector3( -.5*width + (r+1)*width/nw2, -.5*width, height*(h+0)/float(nh) )
 
94
                        v13 = Vector3( -.5*width + (r+1)*width/nw2, -.5*width, height*(h+1)/float(nh) )
 
95
                        v14 = Vector3( -.5*width + (r+0)*width/nw2, -.5*width, height*(h+1)/float(nh) )
 
96
                        f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
 
97
                        f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
 
98
                        v21 = Vector3( +.5*width, -.5*width + (r+0)*width/nw2, height*(h+0)/float(nh) )
 
99
                        v22 = Vector3( +.5*width, -.5*width + (r+1)*width/nw2, height*(h+0)/float(nh) )
 
100
                        v23 = Vector3( +.5*width, -.5*width + (r+1)*width/nw2, height*(h+1)/float(nh) )
 
101
                        v24 = Vector3( +.5*width, -.5*width + (r+0)*width/nw2, height*(h+1)/float(nh) )
 
102
                        f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
 
103
                        f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
 
104
                        v31 = Vector3( +.5*width - (r+0)*width/nw2, +.5*width, height*(h+0)/float(nh) )
 
105
                        v32 = Vector3( +.5*width - (r+1)*width/nw2, +.5*width, height*(h+0)/float(nh) )
 
106
                        v33 = Vector3( +.5*width - (r+1)*width/nw2, +.5*width, height*(h+1)/float(nh) )
 
107
                        v34 = Vector3( +.5*width - (r+0)*width/nw2, +.5*width, height*(h+1)/float(nh) )
 
108
                        f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
 
109
                        f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
 
110
                        v41 = Vector3( -.5*width, +.5*width - (r+0)*width/nw2, height*(h+0)/float(nh) )
 
111
                        v42 = Vector3( -.5*width, +.5*width - (r+1)*width/nw2, height*(h+0)/float(nh) )
 
112
                        v43 = Vector3( -.5*width, +.5*width - (r+1)*width/nw2, height*(h+1)/float(nh) )
 
113
                        v44 = Vector3( -.5*width, +.5*width - (r+0)*width/nw2, height*(h+1)/float(nh) )
 
114
                        f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
 
115
                        f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
 
116
                        facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
 
117
O.bodies.append(facets)
 
118
mass = O.bodies[0].state.mass
 
119
for f in facets:
 
120
        f.state.mass = mass
 
121
        f.state.blockedDOFs = 'XYZz'
 
122
 
 
123
# plots 
 
124
plot.plots = { 'e':('s',), }
 
125
def plotAddData():
 
126
        f1 = sum(O.forces.f(b.id)[2] for b in top)
 
127
        f2 = sum(O.forces.f(b.id)[2] for b in bot)
 
128
        f = .5*(f2-f1)
 
129
        s = f/(pi*.25*width*width) if testType=='cyl' else f/(width*width) if testType=='cube' else None
 
130
        e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
 
131
        plot.addData(
 
132
                i = O.iter,
 
133
                s = s,
 
134
                e = e,
 
135
        )
 
136
 
 
137
# apply prestress to facets
 
138
def addForces():
 
139
        for f in facets:
 
140
                n = f.shape.normal
 
141
                a = f.shape.area
 
142
                O.forces.addF(f.id,preStress*a*n)
 
143
 
 
144
# stop condition and exit of the simulation
 
145
def stopIfDamaged(maxEps=5e-3):
 
146
        extremum = max(abs(s) for s in plot.data['s'])
 
147
        s = abs(plot.data['s'][-1])
 
148
        e = abs(plot.data['e'][-1])
 
149
        if O.iter < 1000 or s > .5*extremum and e < maxEps:
 
150
                return
 
151
        f = os.path.join(exportDir,fileName)
 
152
        print 'gnuplot',plot.saveGnuplot(f,term='png')
 
153
        if runGnuplot:
 
154
                import subprocess
 
155
                os.chdir(exportDir)
 
156
                subprocess.Popen(['gnuplot',f+'.gnuplot']).wait()
 
157
        print 'Simulation finished'
 
158
        O.pause()
 
159
        #sys.exit(0) # results in some threading exception
 
160
 
 
161
O.dt=.5*utils.PWaveTimeStep()
 
162
enlargeFactor=1.5
 
163
O.engines=[
 
164
        ForceResetter(),
 
165
        InsertionSortCollider([
 
166
                Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
 
167
                Bo1_Facet_Aabb()
 
168
        ]),
 
169
        InteractionLoop(
 
170
                [
 
171
                        Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
 
172
                        Ig2_Facet_Sphere_ScGeom(),
 
173
                ],
 
174
                [
 
175
                        Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5),
 
176
                        Ip2_FrictMat_CpmMat_FrictPhys(),
 
177
                        Ip2_FrictMat_FrictMat_FrictPhys(),
 
178
                ],
 
179
                [
 
180
                        Law2_ScGeom_CpmPhys_Cpm(),
 
181
                        Law2_ScGeom_FrictPhys_CundallStrack(),
 
182
                ],
 
183
        ),
 
184
        PyRunner(iterPeriod=1,command="addForces()"),
 
185
        NewtonIntegrator(damping=.3),
 
186
        CpmStateUpdater(iterPeriod=50,label='cpmStateUpdater'),
 
187
        PyRunner(command='plotAddData()',iterPeriod=10),
 
188
        PyRunner(iterPeriod=50,command='stopIfDamaged()'),
 
189
]
 
190
 
 
191
# run one step
 
192
O.step()
 
193
 
 
194
# reset interaction detection enlargement
 
195
bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0
 
196
 
 
197
# initialize auto-updated plot
 
198
if runInGui:
 
199
        plot.plot()
 
200
        try:
 
201
                from yade import qt
 
202
                renderer=qt.Renderer()
 
203
                # uncomment following line to exagerate displacement
 
204
                #renderer.dispScale=(100,100,100)
 
205
        except:
 
206
                pass
 
207
 
 
208
# run
 
209
O.run()