# 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
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