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.
32     radius [.5]
33       Cut-out radius.
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, '''
114      eNpjaGDADhlwiOEU1z8HZusbgukkg5BzRJqKFRoa1oD1HzfceA5NH9FmgKC10SuwOdONpM7DxDYa77gM
115      MueoMQPDEePzV2Hic42XXmoynnQRxvc3dryQbnz3Aoj91Mj3vJnx6fOqxjnnAQzkV94=''')
116    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
117      eNoNzE8og3EcBvC3uUo5rNUOnBSK9/19n0Ic0Eo5oJBmRxcaB04kUnPgoETmT2w7LVrtMBy4auMw+35/
118      7/vaykFSFEopKTnIe/jU01PPU6FNWcQIn+Or5CBfSqCGD1uDYhi7/KbW+dma5aK65gX6Y8Po8HSzZQ7y
119      vBniHyvFV9aq17V7TK42O9kwFS9YUzxhjXIcZxLCnIzjTsfxah/BMFJotjUlZYz6xYeoPqEPKaigbKhb
120      9lOj9NGa9KgtVmqJH9UT36gcp71dEr6HaVS5GS8f46AcQ9itx739SQXdBL8dRqeTo1odox35poh2yJVh
121      apEueucsRWWPgpJFoLKPNzeHC/fU+yl48pDyMi6dCFbsBNJODNu2iawOoE4PoVdP4kH/UkZeaEDaUJQG
122      zMg/DouRUg==''')