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

create the coarsest level parameter domain

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

create geometry function

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

refine domain

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

construct residual

82    res = domain.integral('∇_j(ubasis_ni) stress_ij dV' @ ns, degree=9)

solve system

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

vizualize result

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

evaluate error

93    err = function.sqrt(domain.integral(['du_k du_k dV', '∇_j(du_i) ∇_j(du_i) dV'] @ ns, degree=9)).eval(lhs=lhs)
94    treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
95
96    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).

104if __name__ == '__main__':
105    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.

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