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

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''')
```