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.λ = 2 * poisson
49    ns.μ = 1 - poisson
50    ns.add_field(('u', 'v'), domain.basis(btype, degree=degree), shape=[2])
51    ns.x_i = 'X_i + u_i'
52    ns.ε_ij = '(∇_j(u_i) + ∇_i(u_j)) / 2'
53    ns.σ_ij = 'λ ε_kk δ_ij + 2 μ ε_ij'
54    ns.r2 = 'X_k X_k'
55    ns.R2 = radius**2 / ns.r2
56    ns.k = (3-poisson) / (1+poisson)  # plane stress parameter
57    ns.scale = traction * (1+poisson) / 2
58    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))'
59    ns.du_i = 'u_i - uexact_i'
60
61    sqr = domain.boundary['left,bottom'].integral('(u_i n_i)^2 dS' @ ns, degree=degree*2)
62    cons = solver.optimize('u,', sqr, droptol=1e-15)
63    sqr = domain.boundary['top,right'].integral('du_k du_k dS' @ ns, degree=20)
64    cons = solver.optimize('u,', sqr, droptol=1e-15, constrain=cons)
65
66    res = domain.integral('∇_j(v_i) σ_ij dV' @ ns, degree=degree*2)
67    args = solver.solve_linear('u:v', res, constrain=cons)
68
69    bezier = domain.sample('bezier', 5)
70    x, σxx = bezier.eval(['x_i', 'σ_00'] @ ns, **args)
71    export.triplot('stressxx.png', x, σxx, tri=bezier.tri, hull=bezier.hull)
72
73    err = domain.integral(function.stack(['du_k du_k dV', '∇_j(du_i) ∇_j(du_i) dV'] @ ns), degree=max(degree, 3)*2).eval(**args)**.5
74    treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
75
76    return err, cons['u'], args['u']
```

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).

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

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