# 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