finitestrain.py¶
In this script we solve the nonlinear Saint Venant-Kichhoff problem on a unit square domain (optionally with a circular cutout), clamped at both the left and right boundary in such a way that an arc is formed over a specified angle. The configuration is constructed such that a symmetric solution is expected.
9from nutils import mesh, function, solver, export, cli, testing
10from nutils.expression_v2 import Namespace
11import numpy
12
13
14def main(nelems: int, etype: str, btype: str, degree: int, poisson: float, angle: float, restol: float, trim: bool):
15 '''
16 Deformed hyperelastic plate.
17
18 .. arguments::
19
20 nelems [10]
21 Number of elements along edge.
22 etype [square]
23 Type of elements (square/triangle/mixed).
24 btype [std]
25 Type of basis function (std/spline).
26 degree [1]
27 Polynomial degree.
28 poisson [.25]
29 Poisson's ratio, nonnegative and stricly smaller than 1/2.
30 angle [20]
31 Rotation angle for right clamp (degrees).
32 restol [1e-8]
33 Newton tolerance.
34 trim [no]
35 Create circular-shaped hole.
36 '''
37
38 domain, geom = mesh.unitsquare(nelems, etype)
39 if trim:
40 domain = domain.trim(function.norm2(geom-.5)-.2, maxrefine=2)
41 bezier = domain.sample('bezier', 5)
42
43 ns = Namespace()
44 ns.x = geom
45 ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
46 ns.angle = angle * numpy.pi / 180
47 ns.lmbda = 2 * poisson
48 ns.mu = 1 - 2 * poisson
49 ns.ubasis = domain.basis(btype, degree=degree)
50 ns.u = function.dotarg('lhs', ns.ubasis, shape=(domain.ndims,))
51 ns.X_i = 'x_i + u_i'
52 ns.strain_ij = '.5 (∇_j(u_i) + ∇_i(u_j))'
53 ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
54
55 sqr = domain.boundary['left'].integral('u_k u_k dS' @ ns, degree=degree*2)
56 sqr += domain.boundary['right'].integral('((u_0 - x_1 sin(2 angle) - cos(angle) + 1)^2 + (u_1 - x_1 (cos(2 angle) - 1) + sin(angle))^2) dS' @ ns, degree=degree*2)
57 cons = solver.optimize('lhs', sqr, droptol=1e-15)
58
59 energy = domain.integral('energy dV' @ ns, degree=degree*2)
60 lhs0 = solver.optimize('lhs', energy, constrain=cons)
61 X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs0)
62 export.triplot('linear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
63
64 ns.strain_ij = '.5 (∇_j(u_i) + ∇_i(u_j) + ∇_i(u_k) ∇_j(u_k))'
65 ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
66
67 energy = domain.integral('energy dV' @ ns, degree=degree*2)
68 lhs1 = solver.minimize('lhs', energy, arguments=dict(lhs=lhs0), constrain=cons).solve(restol)
69 X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs1)
70 export.triplot('nonlinear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
71
72 return lhs0, lhs1
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
finitestrain.py
(view log). To select quadratic splines and a cutout add python3
finitestrain.py btype=spline degree=2 trim=yes
(view log).
81if __name__ == '__main__':
82 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.
91class test(testing.TestCase):
92
93 def test_default(self):
94 lhs0, lhs1 = main(nelems=4, etype='square', btype='std', degree=1, poisson=.25, angle=10, restol=1e-8, trim=False)
95 with self.subTest('linear'):
96 self.assertAlmostEqual64(lhs0, '''
97 eNpjYMAE5ZeSL/HqJ146YeB4cbvhl/PzjPrOcVy8da7b4Og5W6Osc/rGt88+MvY+u+yC7NlcQ+GzEsYP
98 z/w3nn1mvon7mdsXJM8oG304vdH45Oluk2WnlU1bTgMAv04qwA==''')
99 with self.subTest('non-linear'):
100 self.assertAlmostEqual64(lhs1, '''
101 eNpjYMAEZdrKl2/p37soY1h84aKh2/lmI4Zz7loq5y0MD55rNtI652Rcefa48aUzzZcjzj4ylDjrYnz6
102 jIBJ8Zl2E9Yzty9InlE2+nB6o/HJ090my04rm7acBgAKcSdV''')
103
104 def test_mixed(self):
105 lhs0, lhs1 = main(nelems=4, etype='mixed', btype='std', degree=1, poisson=.25, angle=10, restol=1e-8, trim=False)
106 with self.subTest('linear'):
107 self.assertAlmostEqual64(lhs0, '''
108 eNpjYICAqxfbL+Xov7kIYi80OA+mtxleOA+iVxjNPBdncOdc6sXT51yNgs8ZGX89e8/Y66zqBaOz/Ya8
109 Z4WMX575ZTz5zAqTgDPKRh9O374geWaj8cnT3SbLTiubtpwGAJ6hLHk=''')
110 with self.subTest('non-linear'):
111 self.assertAlmostEqual64(lhs1, '''
112 eNpjYIAA7fv2l6UMEi6C2H8N7l0A0VcMzc+D6H4jznPyhpfOdelwnm80EjznYTz57CnjG2eWX0o/+9VQ
113 +KyT8cUzzCbZZ2abiJ9RNvpw+vYFyTMbjU+e7jZZdlrZtOU0AJN4KHY=''')
114
115 def test_spline(self):
116 lhs0, lhs1 = main(nelems=4, etype='square', btype='spline', degree=2, poisson=.25, angle=10, restol=1e-8, trim=False)
117 with self.subTest('linear'):
118 self.assertAlmostEqual64(lhs0, '''
119 eNpjYMAOrl3J0vmixaY7QS9N545+w9VaA5eLXYZp51MvVl/I1F164YeBxAVlI//zzMZB52KN35+dd+H9
120 2Vd6b85yGx0/a22cd/aXMetZH5PTZ7ZfaDmzTL/nzFGj3DPPje3OLDBhPvPC5N7p2xckz/gZsJwRML5z
121 Wstk++m7JlNPK5u2nAYATqg9sA==''')
122 with self.subTest('non-linear'):
123 self.assertAlmostEqual64(lhs1, '''
124 eNpjYMAOnLUP6ejq9ukI67vflTVQvdRt0H8h3fDBOT7trReK9adeyDFcez7YaN+5X0Z7z7oYB5/9rKx9
125 ztdA6Fyq0dqzScbGZ78bLzmja5J8RvzSrjN9BgvOfDFKP/PTWOfMSpO3p8+YbDx9+4LkGT8DljMCxndO
126 a5lsP33XZOppZdOW0wApLzra''')