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

```59if __name__ == '__main__':
60  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.

``` 68class test(testing.TestCase):
69
70  @testing.requires('matplotlib')
71  def test_default(self):
72    cons, lhs = main(nelems=4, etype='square', btype='std', degree=1, poisson=.25)
73    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
74      eNpjYMACGsiHP0wxMQBKlBdi''')
75    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
76      eNpjYMAEKcaiRmLGQQZCxgwMYsbrzqcYvz672KTMaIKJimG7CQPDBJM75xabdJ3NMO0xSjG1MUw0Beox
77      PXIuw7Tk7A/TXqMfQLEfQLEfQLEfpsVnAUzzHtI=''')
78
79  @testing.requires('matplotlib')
80  def test_mixed(self):
81    cons, lhs = main(nelems=4, etype='mixed', btype='std', degree=1, poisson=.25)
82    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
83      eNpjYICCBiiEsdFpIuEPU0wMAG6UF2I=''')
84    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
85      eNpjYICAJGMOI3ljcQMwx3i/JohSMr51HkQnGP8422eiYrjcJM+o3aToWq/Jy3PLTKafzTDtM0oxtTRM
86      MF2okmJ67lyGacnZH6aOhj9Mu41+mMZq/DA9dO6HaflZAAMdIls=''')
87
88  @testing.requires('matplotlib')
90    cons, lhs = main(nelems=4, etype='square', btype='std', degree=2, poisson=.25)
91    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
92      eNpjYCACNIxc+MOUMAYA/+NOFg==''')
93    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
94      eNqFzL9KA0EQx/HlLI5wprBJCol/rtfN7MxobZEXOQIJQdBCwfgAItwVStQmZSAvcOmtVW6z5wP4D2yE