elasticity.py¶
In this script we solve the linear elasticity problem on a unit square domain, clamped at the left boundary, and stretched at the right boundary while keeping vertical displacements free.
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 | from nutils import mesh, function, solver, export, cli, testing
def main(nelems:int, etype:str, btype:str, degree:int, poisson:float):
'''
Horizontally loaded linear elastic plate.
.. arguments::
nelems [10]
Number of elements along edge.
etype [square]
Type of elements (square/triangle/mixed).
btype [std]
Type of basis function (std/spline), with availability depending on the
configured element type.
degree [1]
Polynomial degree.
poisson [.25]
Poisson's ratio, nonnegative and strictly smaller than 1/2.
'''
domain, geom = mesh.unitsquare(nelems, etype)
ns = function.Namespace()
ns.x = geom
ns.basis = domain.basis(btype, degree=degree).vector(2)
ns.u_i = 'basis_ni ?lhs_n'
ns.X_i = 'x_i + u_i'
ns.lmbda = 2 * poisson
ns.mu = 1 - 2 * poisson
ns.strain_ij = '(u_i,j + u_j,i) / 2'
ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
sqr = domain.boundary['left'].integral('u_k u_k d:x' @ ns, degree=degree*2)
sqr += domain.boundary['right'].integral('(u_0 - .5)^2 d:x' @ ns, degree=degree*2)
cons = solver.optimize('lhs', sqr, droptol=1e-15)
res = domain.integral('basis_ni,j stress_ij d:x' @ ns, degree=degree*2)
lhs = solver.solve_linear('lhs', res, constrain=cons)
bezier = domain.sample('bezier', 5)
X, sxy = bezier.eval(['X_i', 'stress_01'] @ ns, lhs=lhs)
export.triplot('shear.png', X, sxy, tri=bezier.tri, hull=bezier.hull)
return cons, lhs
|
If the script is executed (as opposed to imported), nutils.cli.run()
calls the main function with arguments provided from the command line. For
example, to keep with the default arguments simply run python3
elasticity.py
(view log). To select mixed elements and quadratic basis functions add
python3 elasticity.py etype=mixed degree=2
(view log).
59 60 | if __name__ == '__main__':
cli.run(main)
|
Once a simulation is developed and tested, it is good practice to save a few
strategic return values for regression testing. The nutils.testing
module, which builds on the standard unittest
framework, facilitates
this by providing nutils.testing.TestCase.assertAlmostEqual64()
for the
embedding of desired results as compressed base64 data.
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 | class test(testing.TestCase):
@testing.requires('matplotlib')
def test_default(self):
cons, lhs = main(nelems=4, etype='square', btype='std', degree=1, poisson=.25)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNpjYMACGsiHP0wxMQBKlBdi''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNpjYMAEKcaiRmLGQQZCxgwMYsbrzqcYvz672KTMaIKJimG7CQPDBJM75xabdJ3NMO0xSjG1MUw0Beox
PXIuw7Tk7A/TXqMfQLEfQLEfQLEfpsVnAUzzHtI=''')
@testing.requires('matplotlib')
def test_mixed(self):
cons, lhs = main(nelems=4, etype='mixed', btype='std', degree=1, poisson=.25)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNpjYICCBiiEsdFpIuEPU0wMAG6UF2I=''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNpjYICAJGMOI3ljcQMwx3i/JohSMr51HkQnGP8422eiYrjcJM+o3aToWq/Jy3PLTKafzTDtM0oxtTRM
MF2okmJ67lyGacnZH6aOhj9Mu41+mMZq/DA9dO6HaflZAAMdIls=''')
@testing.requires('matplotlib')
def test_quadratic(self):
cons, lhs = main(nelems=4, etype='square', btype='std', degree=2, poisson=.25)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNpjYCACNIxc+MOUMAYA/+NOFg==''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNqFzL9KA0EQx/HlLI5wprBJCol/rtfN7MxobZEXOQIJQdBCwfgAItwVStQmZSAvcOmtVW6z5wP4D2yE
aKOwEhTnDRz4VvPhp9T/1zeP0ILF5hhSnUK5cQlKpaDvx3DoWvA57Zt128PIMO5CjHvNOn5s1lCpOi6V
MZ5PGS/k/1U0qGcqVMIcQ5jhmX4XM8N9N8dvWyFtG3RVjOjADOkNBrQMGV3rlJTKaMcN6NUOqWZHlBVV
PjER/0DIDAE/6ICVCjh2Id/ZiBdslY+LrpiOmLaYhJ90IibhNdcW0xHTFTPhUzPhX8h5W3rRuZicV1zO
N3bCgXRUeDFedjxvSc/ai/G86jzfWi87Xswfg5Nx3Q==''')
@testing.requires('matplotlib')
def test_poisson(self):
cons, lhs = main(nelems=4, etype='square', btype='std', degree=1, poisson=.4)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNpjYMACGsiHP0wxMQBKlBdi''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNpjYMAEFsaTjdcYvTFcasTAsMZI5JyFce6ZKSavjbNMFhhFmPz/n2WScHaKieiZRFMmk3DTrUaBpv//
h5t6n000/Xf6hymLyQ/TbUY/gGI/TL3O/jD9cxoASiglXw==''')
|