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
import nutils, numpy

The main function defines the parameter space for the script. Configurable parameters are the number of uniform refinements starting from a 1x2 base mesh, the far field traction, cutout radius and Poisson’s ratio.

16
17
18
19
def main(nrefine = 2,
         traction: "far field traction (relative to Young's modulus)" = .1,
         radius: 'hole radius' = 0.5,
         poisson: "poisson's ratio" = 0.3):

create the coarsest level parameter domain

22
23
24
25
26
27
  domain, geom0 = nutils.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

30
31
32
33
34
35
36
37
38
  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((nutils.function.norm2(geom) - radius)**2, degree=9).eval()**.5
  nutils.log.info('hole radius exact up to L2 error {:.2e}'.format(radiuserr))

refine domain

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
  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 = nutils.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' @ ns, degree=9)
  sqr += domain.boundary['right'].integral('du_k du_k d:x' @ ns, degree=20)
  cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)

construct residual

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

solve system

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

vizualize result

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

evaluate error

79
80
81
82
  err = domain.integral('<du_k du_k, du_i,j du_i,j>_n d:x' @ ns, degree=9).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-nurbs.py (view log).

89
90
if __name__ == '__main__':
  nutils.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.

 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
class test(nutils.testing.TestCase):

  @nutils.testing.requires('matplotlib')
  def test0(self):
    err, cons, lhs = main(nrefine=0)
    with self.subTest('l2-error'):
      self.assertAlmostEqual(err[0], .00189, places=5)
    with self.subTest('h1-error'):
      self.assertAlmostEqual(err[1], .02174, places=5)
    with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
      eNqLuMnQAIJnNSC08XUNYznjy8YQXoYmhC6+Insu7/ze873aAH5yESc=''')
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNoBMADP/1jZ0DGVM5YzzSjfL2kzqDMz1ygzHjPTM5LOr85F0GgpJc6GzrIuc9Qdzm7Pvc+NKyFrF1c=''')

  @nutils.testing.requires('matplotlib')
  def test2(self):
    err, cons, lhs = main(nrefine=2)
    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, '''
      eNrzec7QgA7NpTDFuJ5iiuXJYIpteIgpdkM+ysDa6KmRqrGqsYVxqvEO4xPGmKoSsZgm8whTbLEcptjT
      +5hifxU1z3mce3/O83zu+bzzJ88/u7DxSvUdAExTStA=''')
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNoB8AAP/0zn8i4cMRsyuTIiM2UzjDOdM50zNxpxLqEwuDF5MgQzWzONM6IzojMK5bQuwDC2MVwy3DI8
      M4AzpTOlM24cUC9CMRgykjLbMiUzcDOrM68zsOH9L+UxnDLsMgMzJzNlM7MzvjPYH1owOzLlMiUzJTM4
      M2UzuDPIM4nOic6azsHOA89szwvQC9E102EcWc5YznjOw85EzwrQD9Fd0pzUHOJDzkTOf87xzpvPiNC7
      0UfTqNWjHjrOQM6ozjLPus880B/RyNIz1uXfMM5EztrOVc+Tz7bPOdCE0ZfV/SEpzkjO785Jz23Pbs/J
      z+bQsdR73GyTdWc=''')