platewithhole-nurbs.py

In this script we solve the same infinite plane strain problem as in platewithhole.py, but instead of using FCM to create the hole we use a NURBS-based mapping. A detailed description of the testcase can be found in Hughes et al., Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering, Elsevier, 2005, 194, 4135-4195.

10from nutils import mesh, function, solver, export, cli, testing
11import numpy, treelog
12
13def main(nrefine:int, traction:float, radius:float, poisson:float):
14  '''
15  Horizontally loaded linear elastic plate with IGA hole.
16
17  .. arguments::
18
19     nrefine [2]
20       Number of uniform refinements starting from 1x2 base mesh.
21     traction [.1]
22       Far field traction (relative to Young's modulus).
23     radius [.5]
24       Cut-out radius.
25     poisson [.3]
26       Poisson's ratio, nonnegative and strictly smaller than 1/2.
27  '''

create the coarsest level parameter domain

30  domain, geom0 = mesh.rectilinear([1, 2])
31  bsplinebasis = domain.basis('spline', degree=2)
32  controlweights = numpy.ones(12)
33  controlweights[1:3] = .5 + .25 * numpy.sqrt(2)
34  weightfunc = bsplinebasis.dot(controlweights)
35  nurbsbasis = bsplinebasis * controlweights / weightfunc

create geometry function

38  indices = [0,2], [1,2], [2,1], [2,0]
39  controlpoints = numpy.concatenate([
40    numpy.take([0, 2**.5-1, 1], indices) * radius,
41    numpy.take([0, .3, 1], indices) * (radius+1) / 2,
42    numpy.take([0, 1, 1], indices)])
43  geom = (nurbsbasis[:,numpy.newaxis] * controlpoints).sum(0)
44
45  radiuserr = domain.boundary['left'].integral((function.norm2(geom) - radius)**2 * function.J(geom0), degree=9).eval()**.5
46  treelog.info('hole radius exact up to L2 error {:.2e}'.format(radiuserr))

refine domain

49  if nrefine:
50    domain = domain.refine(nrefine)
51    bsplinebasis = domain.basis('spline', degree=2)
52    controlweights = domain.project(weightfunc, onto=bsplinebasis, geometry=geom0, ischeme='gauss9')
53    nurbsbasis = bsplinebasis * controlweights / weightfunc
54
55  ns = function.Namespace()
56  ns.x = geom
57  ns.lmbda = 2 * poisson
58  ns.mu = 1 - poisson
59  ns.ubasis = nurbsbasis.vector(2)
60  ns.u_i = 'ubasis_ni ?lhs_n'
61  ns.X_i = 'x_i + u_i'
62  ns.strain_ij = '(u_i,j + u_j,i) / 2'
63  ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
64  ns.r2 = 'x_k x_k'
65  ns.R2 = radius**2 / ns.r2
66  ns.k = (3-poisson) / (1+poisson) # plane stress parameter
67  ns.scale = traction * (1+poisson) / 2
68  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))'
69  ns.du_i = 'u_i - uexact_i'
70
71  sqr = domain.boundary['top,bottom'].integral('(u_i n_i)^2 d:x' @ ns, degree=9)
72  cons = solver.optimize('lhs', sqr, droptol=1e-15)
73  sqr = domain.boundary['right'].integral('du_k du_k d:x' @ ns, degree=20)
74  cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)

construct residual

77  res = domain.integral('ubasis_ni,j stress_ij d:x' @ ns, degree=9)

solve system

80  lhs = solver.solve_linear('lhs', res, constrain=cons)

vizualize result

83  bezier = domain.sample('bezier', 9)
84  X, stressxx = bezier.eval(['X_i', 'stress_00'] @ ns, lhs=lhs)
85  export.triplot('stressxx.png', X, stressxx, tri=bezier.tri, hull=bezier.hull, clim=(numpy.nanmin(stressxx), numpy.nanmax(stressxx)))

evaluate error

88  err = domain.integral('<du_k du_k, du_i,j du_i,j>_n d:x' @ ns, degree=9).eval(lhs=lhs)**.5
89  treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
90
91  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-nurbs.py (view log).

98if __name__ == '__main__':
99  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.

107class test(testing.TestCase):
108
109  @testing.requires('matplotlib')
110  def test0(self):
111    err, cons, lhs = main(nrefine=0, traction=.1, radius=.5, poisson=.3)
112    with self.subTest('l2-error'):
113      self.assertAlmostEqual(err[0], .00199, places=5)
114    with self.subTest('h1-error'):
115      self.assertAlmostEqual(err[1], .02269, places=5)
116    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
117      eNpjYGBoQIIggMZXOKdmnHRe3vjh+cvGDAwA6w0LgQ==''')
118    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
119      eNpjYJh07qLhhnOTjb0vTDdmAAKVcy/1u85lGYforQDzFc6pGSedlzd+eP4ykA8AvkQRaA==''')
120
121  @testing.requires('matplotlib')
122  def test2(self):
123    err, cons, lhs = main(nrefine=2, traction=.1, radius=.5, poisson=.3)
124    with self.subTest('l2-error'):
125      self.assertAlmostEqual(err[0], .00009, places=5)
126    with self.subTest('h1-error'):
127      self.assertAlmostEqual(err[1], .00286, places=5)
128    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
129      eNpjYGBoIAKCwCBXp3kuysDjnLXR+3NPjTzPqxrnAnHeeQvjk+dTjZ9d2GG85soJYwYGAPkhPtE=''')
130    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
131      eNpjYOg890mv85yM4axz0kYHz+00Yj6vZJxzPtWY+0KPMffFucaml+caMwBB5LlCvYhzCw0qzu0wPHyu
132      0sjlPIsx14VoY/6LvcaxlxYZz7myCKzO+dwWPZdzBwzqz20z/Hguxmj2+TtGHRdsjHdfbDB2v7zUeMXV
133      pWB1VucC9B3OORmuOCdhZHR+ktGu87eNbC6oGstfLDA+eWm1seG19WB1Buf+6ruce2p469wco9Dzb4wm
134      n2c23nZe3djqQqpx88XNxrOv7gOr0zwXZeBxztro/bmnRp7nVY1zgTjvvIXxSaBfnl3YYbzmygmgOgDU
135      Imlr''')