# platewithhole.py¶

In this script we solve the linear plane strain elasticity problem for an infinite plate with a circular hole under tension. We do this by placing the circle in the origin of a unit square, imposing symmetry conditions on the left and bottom, and Dirichlet conditions constraining the displacements to the analytical solution to the right and top. The traction-free circle is removed by means of the Finite Cell Method (FCM).

```10from nutils import mesh, function, solver, export, cli, testing
11import numpy, treelog
12
13def main(nelems:int, etype:str, btype:str, degree:int, traction:float, maxrefine:int, radius:float, poisson:float):
14  '''
15  Horizontally loaded linear elastic plate with FCM hole.
16
17  .. arguments::
18
19     nelems [9]
20       Number of elements along edge.
21     etype [square]
22       Type of elements (square/triangle/mixed).
23     btype [std]
24       Type of basis function (std/spline), with availability depending on the
25       selected element type.
26     degree [2]
27       Polynomial degree.
28     traction [.1]
29       Far field traction (relative to Young's modulus).
30     maxrefine [2]
31       Number or refinement levels used for the finite cell method.
34     poisson [.3]
35       Poisson's ratio, nonnegative and strictly smaller than 1/2.
36  '''
37
38  domain0, geom = mesh.unitsquare(nelems, etype)
39  domain = domain0.trim(function.norm2(geom) - radius, maxrefine=maxrefine)
40
41  ns = function.Namespace()
42  ns.x = geom
43  ns.lmbda = 2 * poisson
44  ns.mu = 1 - poisson
45  ns.ubasis = domain.basis(btype, degree=degree).vector(2)
46  ns.u_i = 'ubasis_ni ?lhs_n'
47  ns.X_i = 'x_i + u_i'
48  ns.strain_ij = '(u_i,j + u_j,i) / 2'
49  ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
50  ns.r2 = 'x_k x_k'
51  ns.R2 = radius**2 / ns.r2
52  ns.k = (3-poisson) / (1+poisson) # plane stress parameter
53  ns.scale = traction * (1+poisson) / 2
54  ns.uexact_i = 'scale (x_i ((k + 1) (0.5 + R2) + (1 - R2) R2 (x_0^2 - 3 x_1^2) / r2) - 2 δ_i1 x_1 (1 + (k - 1 + R2) R2))'
55  ns.du_i = 'u_i - uexact_i'
56
57  sqr = domain.boundary['left,bottom'].integral('(u_i n_i)^2 d:x' @ ns, degree=degree*2)
58  cons = solver.optimize('lhs', sqr, droptol=1e-15)
59  sqr = domain.boundary['top,right'].integral('du_k du_k d:x' @ ns, degree=20)
60  cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
61
62  res = domain.integral('ubasis_ni,j stress_ij d:x' @ ns, degree=degree*2)
63  lhs = solver.solve_linear('lhs', res, constrain=cons)
64
65  bezier = domain.sample('bezier', 5)
66  X, stressxx = bezier.eval(['X_i', 'stress_00'] @ ns, lhs=lhs)
67  export.triplot('stressxx.png', X, stressxx, tri=bezier.tri, hull=bezier.hull)
68
69  err = domain.integral('<du_k du_k, du_i,j du_i,j>_n d:x' @ ns, degree=max(degree,3)*2).eval(lhs=lhs)**.5
70  treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
71
72  return err, 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 platewithhole.py``` (view log). To select mixed elements and quadratic basis functions add `python3 platewithhole.py etype=mixed degree=2` (view log).

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

``` 89class test(testing.TestCase):
90
91  @testing.requires('matplotlib')
92  def test_spline(self):
93    err, cons, lhs = main(nelems=4, etype='square', btype='spline', degree=2, traction=.1, maxrefine=2, radius=.5, poisson=.3)
94    with self.subTest('l2-error'):
95      self.assertAlmostEqual(err[0], .00033, places=5)
96    with self.subTest('h1-error'):
97      self.assertAlmostEqual(err[1], .00672, places=5)
98    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
99      eNpjaGBoYGBAxvrnGBow4X89g3NQFSjQwLAGq7i10Wus4k+NfM8fNWZgOGL89upc47WX0ozvXjAzPn1e
100      1TjnPACrACoJ''')
101    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
102      eNpbZHbajIHhxzkGBhMgtgdi/XPypyRPvjFxO/PccPq5Vn2vcxr6luf+6xmcm2LMwLDQePf5c0bTzx8x
103      5D7vaTjnnIFhzbmlQPH5xhV39Y3vXlxtJHoh2EjvvLXR63MbgOIbjRdfrTXeecnUeO+Fn0Yrzj818j1/
104      FCh+xPjt1bnGay+lGd+9YGZ8+ryqcc55AK+AP/0=''')
105
106  @testing.requires('matplotlib')
107  def test_mixed(self):
108    err, cons, lhs = main(nelems=4, etype='mixed', btype='std', degree=2, traction=.1, maxrefine=2, radius=.5, poisson=.3)
109    with self.subTest('l2-error'):
110      self.assertAlmostEqual(err[0], .00024, places=5)
111    with self.subTest('h1-error'):
112      self.assertAlmostEqual(err[1], .00739, places=5)
113    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''