# 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
10import numpy
11
12def main(nelems:int, etype:str, btype:str, degree:int, poisson:float, angle:float, restol:float, trim:bool):
13  '''
14  Deformed hyperelastic plate.
15
16  .. arguments::
17
18     nelems [10]
19       Number of elements along edge.
20     etype [square]
21       Type of elements (square/triangle/mixed).
22     btype [std]
23       Type of basis function (std/spline).
24     degree [1]
25       Polynomial degree.
26     poisson [.25]
27       Poisson's ratio, nonnegative and stricly smaller than 1/2.
28     angle [20]
29       Rotation angle for right clamp (degrees).
30     restol [1e-10]
31       Newton tolerance.
32     trim [no]
33       Create circular-shaped hole.
34  '''
35
36  domain, geom = mesh.unitsquare(nelems, etype)
37  if trim:
38    domain = domain.trim(function.norm2(geom-.5)-.2, maxrefine=2)
39  bezier = domain.sample('bezier', 5)
40
41  ns = function.Namespace()
42  ns.x = geom
43  ns.angle = angle * numpy.pi / 180
44  ns.lmbda = 2 * poisson
45  ns.mu = 1 - 2 * poisson
46  ns.ubasis = domain.basis(btype, degree=degree).vector(2)
47  ns.u_i = 'ubasis_ki ?lhs_k'
48  ns.X_i = 'x_i + u_i'
49  ns.strain_ij = '.5 (u_i,j + u_j,i)'
50  ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
51
52  sqr = domain.boundary['left'].integral('u_k u_k d:x' @ ns, degree=degree*2)
53  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) d:x' @ ns, degree=degree*2)
54  cons = solver.optimize('lhs', sqr, droptol=1e-15)
55
56  energy = domain.integral('energy d:x' @ ns, degree=degree*2)
57  lhs0 = solver.optimize('lhs', energy, constrain=cons)
58  X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs0)
59  export.triplot('linear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
60
61  ns.strain_ij = '.5 (u_i,j + u_j,i + u_k,i u_k,j)'
62  ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
63
64  energy = domain.integral('energy d:x' @ ns, degree=degree*2)
65  lhs1 = solver.minimize('lhs', energy, lhs0=lhs0, constrain=cons).solve(restol)
66  X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs1)
67  export.triplot('nonlinear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
68
69  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).

```77if __name__ == '__main__':
78  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.

``` 86class test(testing.TestCase):
87
88  @testing.requires('matplotlib')
89  def test_default(self):
90    lhs0, lhs1 = main(nelems=4, etype='square', btype='std', degree=1, poisson=.25, angle=10, restol=1e-10, trim=False)
91    with self.subTest('linear'): self.assertAlmostEqual64(lhs0, '''
92      eNpjYMAE5ZeSL/HqJ146YeB4cbvhl/PzjPrOcVy8da7b4Og5W6Osc/rGt88+MvY+u+yC7NlcQ+GzEsYP
93      z/w3nn1mvon7mdsXJM8oG304vdH45Oluk2WnlU1bTgMAv04qwA==''')
94    with self.subTest('non-linear'): self.assertAlmostEqual64(lhs1, '''
95      eNpjYMAEZdrKl2/p37soY1h84aKh2/lmI4Zz7loq5y0MD55rNtI652Rcefa48aUzzZcjzj4ylDjrYnz6
96      jIBJ8Zl2E9Yzty9InlE2+nB6o/HJ090my04rm7acBgAKcSdV''')
97
98  @testing.requires('matplotlib')
99  def test_mixed(self):
100    lhs0, lhs1 = main(nelems=4, etype='mixed', btype='std', degree=1, poisson=.25, angle=10, restol=1e-10, trim=False)
101    with self.subTest('linear'): self.assertAlmostEqual64(lhs0, '''
102      eNpjYICAqxfbL+Xov7kIYi80OA+mtxleOA+iVxjNPBdncOdc6sXT51yNgs8ZGX89e8/Y66zqBaOz/Ya8
103      Z4WMX575ZTz5zAqTgDPKRh9O374geWaj8cnT3SbLTiubtpwGAJ6hLHk=''')
104    with self.subTest('non-linear'): self.assertAlmostEqual64(lhs1, '''
105      eNpjYIAA7fv2l6UMEi6C2H8N7l0A0VcMzc+D6H4jznPyhpfOdelwnm80EjznYTz57CnjG2eWX0o/+9VQ
106      +KyT8cUzzCbZZ2abiJ9RNvpw+vYFyTMbjU+e7jZZdlrZtOU0AJN4KHY=''')
107
108  @testing.requires('matplotlib')
109  def test_spline(self):
110    lhs0, lhs1 = main(nelems=4, etype='square', btype='spline', degree=2, poisson=.25, angle=10, restol=1e-10, trim=False)
111    with self.subTest('linear'): self.assertAlmostEqual64(lhs0, '''
112      eNpjYMAOrl3J0vmixaY7QS9N545+w9VaA5eLXYZp51MvVl/I1F164YeBxAVlI//zzMZB52KN35+dd+H9