/
finitestrain.py
109 lines (92 loc) · 5.27 KB
/
finitestrain.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#! /usr/bin/env python3
#
# 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.
import nutils, numpy
# The main function defines the parameter space for the script. Configurable
# parameters are the mesh density (in number of elements along an edge),
# element type (square, triangle, or mixed), type of basis function (std or
# spline, with availability depending on element type), polynomial degree,
# Poisson's ratio, wedge angle, Newton tolerance, and a boolean flag for a
# circular cutout.
def main(nelems: 'number of elements along edge' = 10,
etype: 'type of elements (square/triangle/mixed)' = 'square',
btype: 'type of basis function (std/spline)' = 'std',
degree: 'polynomial degree' = 1,
poisson: 'poisson ratio < 0.5' = .25,
angle: 'bend angle (degrees)' = 20,
restol: 'residual tolerance' = 1e-10,
trim: 'create circular-shaped hole' = False):
domain, geom = nutils.mesh.unitsquare(nelems, etype)
if trim:
domain = domain.trim(nutils.function.norm2(geom-.5)-.2, maxrefine=2)
bezier = domain.sample('bezier', 5)
ns = nutils.function.Namespace()
ns.x = geom
ns.angle = angle * numpy.pi / 180
ns.lmbda = 2 * poisson
ns.mu = 1 - 2 * poisson
ns.ubasis = domain.basis(btype, degree=degree).vector(2)
ns.u_i = 'ubasis_ki ?lhs_k'
ns.X_i = 'x_i + u_i'
ns.strain_ij = '.5 (u_i,j + u_j,i)'
ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
sqr = domain.boundary['left'].integral('u_k u_k d:x' @ ns, degree=degree*2)
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)
cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
energy = domain.integral('energy d:x' @ ns, degree=degree*2)
lhs0 = nutils.solver.optimize('lhs', energy, constrain=cons)
X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs0)
nutils.export.triplot('linear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
ns.strain_ij = '.5 (u_i,j + u_j,i + u_k,i u_k,j)'
ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
energy = domain.integral('energy d:x' @ ns, degree=degree*2)
lhs1 = nutils.solver.minimize('lhs', energy, lhs0=lhs0, constrain=cons).solve(restol)
X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs1)
nutils.export.triplot('nonlinear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
return lhs0, lhs1
# If the script is executed (as opposed to imported), :func:`nutils.cli.run`
# calls the main function with arguments provided from the command line. For
# example, to keep with the default arguments simply run :sh:`python3
# finitestrain.py`. To select quadratic splines and a cutout add :sh:`python3
# finitestrain.py btype=spline degree=2 trim`.
if __name__ == '__main__':
nutils.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 :mod:`nutils.testing`
# module, which builds on the standard :mod:`unittest` framework, facilitates
# this by providing :func:`nutils.testing.TestCase.assertAlmostEqual64` for the
# embedding of desired results as compressed base64 data.
class test(nutils.testing.TestCase):
@nutils.testing.requires('matplotlib')
def test_default(self):
lhs0, lhs1 = main(nelems=4, angle=10)
with self.subTest('linear'): self.assertAlmostEqual64(lhs0, '''
eNpjYICB8ku8+icMthvOM+K42G1ga6Rv/Mh42YVcQwnj/8bzTW5fUDbaaNxtomwK18CQfCnxkuPFL+f7
zt06d/Rc1rnbZ73Pyp4VPvvwzOwz7mckz3w4ffL0stMtpwGSOirA''')
with self.subTest('non-linear'): self.assertAlmostEqual64(lhs1, '''
eNpjYICBMu1b+jKGFw2bjdy1LICkk/Fx4+bLjwxdjAVM2k1uX1A22mjcbaJsCtfAoHz53sXiC27nGc6p
nD94Tutc5dlLZyLOSpw9fab4DOsZyTMfTp88vex0y2kA6e4nVQ==''')
@nutils.testing.requires('matplotlib')
def test_mixed(self):
lhs0, lhs1 = main(nelems=4, angle=10, etype='mixed')
with self.subTest('linear'): self.assertAlmostEqual64(lhs0, '''
eNoBZACb/wAAAADV0WwvAAChMAAAtjEAAKgyXjBl0UUyMjPeMyXQjzESM/ozqDQjMtvQsTOLNCM1AAAA
AIfS7NEAAM/RAADQzwAAmc7czsvOU871zUrNMs0NzenMk8xQzPDLGczJy6bLhMsZ2Sx5''')
with self.subTest('non-linear'): self.assertAlmostEqual64(lhs1, '''
eNoBZACb/wAAAAAr3xowAAD9MAAA1DEAAI8yHzGKLIEySDPKM6fS9TFCMwM0mzQjMtvQsTOLNCM1AAAA
AD/TYNEAAN7QAAA3zwAACc7SzgnPEc6TzdjMZ80TzdHMa8wXzPDLGczJy6bLhMthnih2''')
@nutils.testing.requires('matplotlib')
def test_spline(self):
lhs0, lhs1 = main(nelems=4, angle=10, degree=2, btype='spline')
with self.subTest('linear'): self.assertAlmostEqual64(lhs0, '''
eNpjYECAa1e+aE3Qu6Nfa9BlmHoxU/eHgbIRs3Gs8bwLr/S4jayNfxn7mGy/sEz/qNFz4wUmL0xuX/Az
EDDWMrlromyKZAxDlg6bbppOw1WXi2nnqy8svSBxwf980Ln3Z9+ffXP2+Nm8s6xnT59pOdNzJveM3Rnm
M/dOS55hOXPn9PbTU0+3nAYAZeQ9sA==''')
with self.subTest('non-linear'): self.assertAlmostEqual64(lhs1, '''
eNpjYEAAZ21dXWF9WYNug3RDPu1i/RzDYKNfRi7Gn5V9DVKNkoy/G+uaiF/qM/hi9NN4pckZk9sX/AwE
jLVM7poomyIZw3BIp0/H/a7qpf4LD85tvTD1wtrz+87tPRt8Vvuc0Lm1Z43PLjmTfGbXmQVn0s/onHl7
euNpyTMsZ+6c3n566umW0wB4sDra''')