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.

 7from nutils import mesh, function, solver, export, cli, testing
 8from nutils.expression_v2 import Namespace
 9
10
11def main(nelems: int, etype: str, btype: str, degree: int, poisson: float):
12    '''
13    Horizontally loaded linear elastic plate.
14
15    .. arguments::
16
17       nelems [10]
18         Number of elements along edge.
19       etype [square]
20         Type of elements (square/triangle/mixed).
21       btype [std]
22         Type of basis function (std/spline), with availability depending on the
23         configured element type.
24       degree [1]
25         Polynomial degree.
26       poisson [.25]
27         Poisson's ratio, nonnegative and strictly smaller than 1/2.
28    '''
29
30    domain, geom = mesh.unitsquare(nelems, etype)
31
32    ns = Namespace()
33    ns.δ = function.eye(domain.ndims)
34    ns.x = geom
35    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
36    ns.basis = domain.basis(btype, degree=degree).vector(2)
37    ns.u = function.dotarg('lhs', ns.basis)
38    ns.X_i = 'x_i + u_i'
39    ns.lmbda = 2 * poisson
40    ns.mu = 1 - 2 * poisson
41    ns.strain_ij = '(∇_j(u_i) + ∇_i(u_j)) / 2'
42    ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
43
44    sqr = domain.boundary['left'].integral('u_k u_k dS' @ ns, degree=degree*2)
45    sqr += domain.boundary['right'].integral('(u_0 - .5)^2 dS' @ ns, degree=degree*2)
46    cons = solver.optimize('lhs', sqr, droptol=1e-15)
47
48    res = domain.integral('∇_j(basis_ni) stress_ij dV' @ ns, degree=degree*2)
49    lhs = solver.solve_linear('lhs', res, constrain=cons)
50
51    bezier = domain.sample('bezier', 5)
52    X, sxy = bezier.eval(['X_i', 'stress_01'] @ ns, lhs=lhs)
53    export.triplot('shear.png', X, sxy, tri=bezier.tri, hull=bezier.hull)
54
55    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).

64if __name__ == '__main__':
65    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.

 74class test(testing.TestCase):
 75
 76    def test_default(self):
 77        cons, lhs = main(nelems=4, etype='square', btype='std', degree=1, poisson=.25)
 78        with self.subTest('constraints'):
 79            self.assertAlmostEqual64(cons, '''
 80                eNpjYMACGsiHP0wxMQBKlBdi''')
 81        with self.subTest('left-hand side'):
 82            self.assertAlmostEqual64(lhs, '''
 83                eNpjYMAEKcaiRmLGQQZCxgwMYsbrzqcYvz672KTMaIKJimG7CQPDBJM75xabdJ3NMO0xSjG1MUw0Beox
 84                PXIuw7Tk7A/TXqMfQLEfQLEfQLEfpsVnAUzzHtI=''')
 85
 86    def test_mixed(self):
 87        cons, lhs = main(nelems=4, etype='mixed', btype='std', degree=1, poisson=.25)
 88        with self.subTest('constraints'):
 89            self.assertAlmostEqual64(cons, '''
 90                eNpjYICCBiiEsdFpIuEPU0wMAG6UF2I=''')
 91        with self.subTest('left-hand side'):
 92            self.assertAlmostEqual64(lhs, '''
 93                eNpjYICAJGMOI3ljcQMwx3i/JohSMr51HkQnGP8422eiYrjcJM+o3aToWq/Jy3PLTKafzTDtM0oxtTRM
 94                MF2okmJ67lyGacnZH6aOhj9Mu41+mMZq/DA9dO6HaflZAAMdIls=''')
 95
 96    def test_quadratic(self):
 97        cons, lhs = main(nelems=4, etype='square', btype='std', degree=2, poisson=.25)
 98        with self.subTest('constraints'):
 99            self.assertAlmostEqual64(cons, '''
100                eNpjYCACNIxc+MOUMAYA/+NOFg==''')
101        with self.subTest('left-hand side'):
102            self.assertAlmostEqual64(lhs, '''
103                eNqFzL9KA0EQx/HlLI5wprBJCol/rtfN7MxobZEXOQIJQdBCwfgAItwVStQmZSAvcOmtVW6z5wP4D2yE
104                aKOwEhTnDRz4VvPhp9T/1zeP0ILF5hhSnUK5cQlKpaDvx3DoWvA57Zt128PIMO5CjHvNOn5s1lCpOi6V
105                MZ5PGS/k/1U0qGcqVMIcQ5jhmX4XM8N9N8dvWyFtG3RVjOjADOkNBrQMGV3rlJTKaMcN6NUOqWZHlBVV
106                PjER/0DIDAE/6ICVCjh2Id/ZiBdslY+LrpiOmLaYhJ90IibhNdcW0xHTFTPhUzPhX8h5W3rRuZicV1zO
107                N3bCgXRUeDFedjxvSc/ai/G86jzfWi87Xswfg5Nx3Q==''')
108
109    def test_poisson(self):
110        cons, lhs = main(nelems=4, etype='square', btype='std', degree=1, poisson=.4)
111        with self.subTest('constraints'):
112            self.assertAlmostEqual64(cons, '''
113                eNpjYMACGsiHP0wxMQBKlBdi''')
114        with self.subTest('left-hand side'):
115            self.assertAlmostEqual64(lhs, '''
116                eNpjYMAEFsaTjdcYvTFcasTAsMZI5JyFce6ZKSavjbNMFhhFmDAwZJkknJ1iInom0ZTJJNx0q1GgKQND
117                uKn32UTTf6d/mLKY/DDdZvQDKPbD1OvsD9M/pwGZyh9l''')