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