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.

10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from nutils import mesh, function, solver, export, cli, testing
import numpy, treelog

def main(nrefine:int, traction:float, radius:float, poisson:float):
  '''
  Horizontally loaded linear elastic plate with IGA hole.

  .. arguments::

     nrefine [2]
       Number of uniform refinements starting from 1x2 base mesh.
     traction [.1]
       Far field traction (relative to Young's modulus).
     radius [.5]
       Cut-out radius.
     poisson [.3]
       Poisson's ratio, nonnegative and strictly smaller than 1/2.
  '''

create the coarsest level parameter domain

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

create geometry function

38
39
40
41
42
43
44
45
46
  indices = [0,2], [1,2], [2,1], [2,0]
  controlpoints = numpy.concatenate([
    numpy.take([0, 2**.5-1, 1], indices) * radius,
    numpy.take([0, .3, 1], indices) * (radius+1) / 2,
    numpy.take([0, 1, 1], indices)])
  geom = (nurbsbasis[:,numpy.newaxis] * controlpoints).sum(0)

  radiuserr = domain.boundary['left'].integral((function.norm2(geom) - radius)**2 * function.J(geom0), degree=9).eval()**.5
  treelog.info('hole radius exact up to L2 error {:.2e}'.format(radiuserr))

refine domain

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
  if nrefine:
    domain = domain.refine(nrefine)
    bsplinebasis = domain.basis('spline', degree=2)
    controlweights = domain.project(weightfunc, onto=bsplinebasis, geometry=geom0, ischeme='gauss9')
    nurbsbasis = bsplinebasis * controlweights / weightfunc

  ns = function.Namespace()
  ns.x = geom
  ns.lmbda = 2 * poisson
  ns.mu = 1 - poisson
  ns.ubasis = nurbsbasis.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['top,bottom'].integral('(u_i n_i)^2 d:x' @ ns, degree=9)
  cons = solver.optimize('lhs', sqr, droptol=1e-15)
  sqr = domain.boundary['right'].integral('du_k du_k d:x' @ ns, degree=20)
  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
84
85
  bezier = domain.sample('bezier', 9)
  X, stressxx = bezier.eval(['X_i', 'stress_00'] @ ns, lhs=lhs)
  export.triplot('stressxx.png', X, stressxx, tri=bezier.tri, hull=bezier.hull, clim=(numpy.nanmin(stressxx), numpy.nanmax(stressxx)))

evaluate error

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

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

107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
class test(testing.TestCase):

  @testing.requires('matplotlib')
  def test0(self):
    err, cons, lhs = main(nrefine=0, traction=.1, radius=.5, poisson=.3)
    with self.subTest('l2-error'):
      self.assertAlmostEqual(err[0], .00199, places=5)
    with self.subTest('h1-error'):
      self.assertAlmostEqual(err[1], .02269, places=5)
    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
      eNpjYGBoQIIggMZXOKdmnHRe3vjh+cvGDAwA6w0LgQ==''')
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNpjYJh07qLhhnOTjb0vTDdmAAKVcy/1u85lGYforQDzFc6pGSedlzd+eP4ykA8AvkQRaA==''')

  @testing.requires('matplotlib')
  def test2(self):
    err, cons, lhs = main(nrefine=2, traction=.1, radius=.5, poisson=.3)
    with self.subTest('l2-error'):
      self.assertAlmostEqual(err[0], .00009, places=5)
    with self.subTest('h1-error'):
      self.assertAlmostEqual(err[1], .00286, places=5)
    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
      eNpjYGBoIAKCwCBXp3kuysDjnLXR+3NPjTzPqxrnAnHeeQvjk+dTjZ9d2GG85soJYwYGAPkhPtE=''')
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNpjYOg890mv85yM4axz0kYHz+00Yj6vZJxzPtWY+0KPMffFucaml+caMwBB5LlCvYhzCw0qzu0wPHyu
      0sjlPIsx14VoY/6LvcaxlxYZz7myCKzO+dwWPZdzBwzqz20z/Hguxmj2+TtGHRdsjHdfbDB2v7zUeMXV
      pWB1VucC9B3OORmuOCdhZHR+ktGu87eNbC6oGstfLDA+eWm1seG19WB1Buf+6ruce2p469wco9Dzb4wm
      n2c23nZe3djqQqpx88XNxrOv7gOr0zwXZeBxztro/bmnRp7nVY1zgTjvvIXxSaBfnl3YYbzmygmgOgDU
      Imlr''')