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