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.
35 radius [.5]
36 Cut-out radius.
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, '''
121 eNpjaGDADhlwiOEU1z8HZusbgukkg5BzRJqKFRoa1oD1HzfceA5NH9FmgKC10SuwOdONpM7DxDYa77gM
122 MueoMQPDEePzV2Hic42XXmoynnQRxvc3dryQbnz3Aoj91Mj3vJnx6fOqxjnnAQzkV94=''')
123 with self.subTest('left-hand side'):
124 self.assertAlmostEqual64(u, '''
125 eNoNzE8og3EcBvC3uUo5rNUOnBSK9/19n0Ic0Eo5oJBmRxcaB04kUnPgoETmT2w7LVrtMBy4auMw+35/
126 7/vaykFSFEopKTnIe/jU01PPU6FNWcQIn+Or5CBfSqCGD1uDYhi7/KbW+dma5aK65gX6Y8Po8HSzZQ7y
127 vBniHyvFV9aq17V7TK42O9kwFS9YUzxhjXIcZxLCnIzjTsfxah/BMFJotjUlZYz6xYeoPqEPKaigbKhb
128 9lOj9NGa9KgtVmqJH9UT36gcp71dEr6HaVS5GS8f46AcQ9itx739SQXdBL8dRqeTo1odox35poh2yJVh
129 apEueucsRWWPgpJFoLKPNzeHC/fU+yl48pDyMi6dCFbsBNJODNu2iawOoE4PoVdP4kH/UkZeaEDaUJQG
130 zMg/DouRUg==''')