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

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

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