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

10
import nutils, numpy

The main function defines the parameter space for the script. Configurable parameters are the mesh density (in number of elements along an edge), element type (square, triangle, or mixed), type of basis function (std or spline, with availability depending on element type), polynomial degree, far field traction, number of refinement levels for FCM, the cutout radius and Poisson’s ratio.

19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def main(nelems: 'number of elementsa long edge' = 9,
         etype: 'type of elements (square/triangle/mixed)' = 'square',
         btype: 'type of basis function (std/spline)' = 'std',
         degree: 'polynomial degree' = 2,
         traction: "far field traction (relative to Young's modulus)" = .1,
         maxrefine: 'maxrefine level for trimming' = 2,
         radius: 'cut-out radius' = .5,
         poisson: 'poisson ratio' = .3):

  domain0, geom = nutils.mesh.unitsquare(nelems, etype)
  domain = domain0.trim(nutils.function.norm2(geom) - radius, maxrefine=maxrefine)

  ns = nutils.function.Namespace()
  ns.x = geom
  ns.lmbda = 2 * poisson
  ns.mu = 1 - poisson
  ns.ubasis = domain.basis(btype, degree=degree).vector(2)
  ns.u_i = 'ubasis_ni ?lhs_n'
  ns.X_i = 'x_i + u_i'
  ns.strain_ij = '(u_i,j + u_j,i) / 2'
  ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
  ns.r2 = 'x_k x_k'
  ns.R2 = radius**2 / ns.r2
  ns.k = (3-poisson) / (1+poisson) # plane stress parameter
  ns.scale = traction * (1+poisson) / 2
  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))'
  ns.du_i = 'u_i - uexact_i'

  sqr = domain.boundary['left,bottom'].integral('(u_i n_i)^2 d:x' @ ns, degree=degree*2)
  sqr += domain.boundary['top,right'].integral('du_k du_k d:x' @ ns, degree=20)
  cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)

  res = domain.integral('ubasis_ni,j stress_ij d:x' @ ns, degree=degree*2)
  lhs = nutils.solver.solve_linear('lhs', res, constrain=cons)

  bezier = domain.sample('bezier', 5)
  X, stressxx = bezier.eval(['X_i', 'stress_00'] @ ns, lhs=lhs)
  nutils.export.triplot('stressxx.png', X, stressxx, tri=bezier.tri, hull=bezier.hull)

  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
  nutils.log.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))

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

69
70
if __name__ == '__main__':
  nutils.cli.run(main)

Once a simulation is developed and tested, it is good practice to save a few strategicly chosen return values for routine regression testing. Here we use the standard unittest framework, with nutils.numeric.assert_allclose64() facilitating the embedding of desired results as compressed base64 data.

 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
class test(nutils.testing.TestCase):

  @nutils.testing.requires('matplotlib')
  def test_spline(self):
    err, cons, lhs = main(nelems=4, etype='square', degree=2, btype='spline')
    nutils.numeric.assert_allclose64(err, 'eNpbpd6jCwAEZgGL')
    nutils.numeric.assert_allclose64(cons, 'eNpjaPC5XybfdX+dIkMDDP7TQ7ANDBFsayME+'
      '6nRUeMjxnON04zNjFWNYaL655B0nrNUgrFrzrHeh7Ffn/sNt8v3/Nk7X66uuXT3wunzOecBJ0s'
      'yCg==')
    nutils.numeric.assert_allclose64(lhs, 'eNoBjABz/6I2TN92H4rfriEeyuw05zGFLykv/i'
      '6UM6EzzjLEMUkxMDGlM58zLzOrMlMyOzKwM7EzfTM1M/ky5TLFM8QznTNmMzYzJTPLNvjONM4/'
      'zi/OGclHzJfOSs45zjDOOSK7z5fPC8+cznzOBd/D3d3RFdAuz+vO+yGg1bnSvdCoz03Pzdz01a'
      'zS3dDLz2zPaQdIRw==')

  @nutils.testing.requires('matplotlib')
  def test_mixed(self):
    err, cons, lhs = main(nelems=4, etype='mixed', degree=2, btype='std')
    nutils.numeric.assert_allclose64(err, 'eNpjU9+jCwACNgEX')
    nutils.numeric.assert_allclose64(cons, 'eNpjaGCAwx4pGMv/8UYZGFvrgagCkNZnaEgyY'
      'GjABw0NGRqOG+JXY23E0DDdCMTaaMzQcNT4iDGIPde4CUz7G6cD6adGZsaqxvjNgUD9c0BbgTj'
      'kHEwk+jE2dTVA+Y3nTsmB2GYPsZv1CqhG6jyItePye8XLd69dBbGXXZp0EUQ7Xrh7gaHB9/zp8'
      'znnAW7uYcc=')
    nutils.numeric.assert_allclose64(lhs, 'eNoNzcsrRGEYB2CxlbKY1CSXhUJxzvf+Co0Fml'
      'IWTCExdjaEBSuTSI0FiymRaxgrl9QsBgu2mqFc3vc75zCliGmQUaKkZCH+gKcnQaM4gI11rFaG'
      '3Gn1aJ6rAPlS0XzTGDG+zWOz/MFVlG1kGAGzx1yAF11YwBo2oKmDMrFDcRVSLmqkeqXUvTpVmw'
      'hjALvYRhCF+KAydCJKQfoim1qpliK0RBEsI4o9xBHDOPz/exAG8uBDL37oiapQghlp48/L2GUO'
      'u2WRp3mIT/iXa7iOW9jLGzzJ1WywhxX3cTvvy7Bc6RerO1VuhaVJ+vWbuOWCS2VKZnmMkxzls4'
      'Ln2yynKrly3encWHHtsjx2rp4Xv3akQl65/1+4E2nn0Hkvdu4S10f2hLVlz1kRXaAb9J3elWY5'
      'l0H5AxDbnCE=')