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
|
# UFL input for hyperleasticity
# =============================
#
# The first step is to define the variational problem at hand. We define
# the variational problem in UFL terms in a separate form file
# :download:`HyperElasticity.ufl`.
#
# We are interested in solving for a discrete vector field in three
# dimensions, so first we need the appropriate finite element space and
# trial and test functions on this space::
# Function spaces
element = VectorElement("Lagrange", tetrahedron, 1)
# Trial and test functions
du = TrialFunction(element) # Incremental displacement
v = TestFunction(element) # Test function
# Note that ``VectorElement`` creates a finite element space of vector
# fields. The dimension of the vector field (the number of components)
# is assumed to be the same as the spatial dimension (in this case 3),
# unless otherwise specified.
#
# Next, we will be needing functions for the boundary source ``B``, the
# traction ``T`` and the displacement solution itself ``u``::
# Functions
u = Coefficient(element) # Displacement from previous iteration
# B = Coefficient(element) # Body force per unit volume
# T = Coefficient(element) # Traction force on the boundary
# Now, we can define the kinematic quantities involved in the model::
# Kinematics
d = len(u)
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Before defining the energy density and thus the total potential
# energy, it only remains to specify constants for the elasticity
# parameters::
# Elasticity parameters
# mu = Constant(tetrahedron)
# lmbda = Constant(tetrahedron)
E = 10.0
nu = 0.3
mu = E/(2*(1 + nu))
lmbda = E*nu/((1 + nu)*(1 - 2*nu))
# Both the first variation of the potential energy, and the Jacobian of
# the variation, can be automatically computed by a call to
# ``derivative``::
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx # - inner(B, u)*dx - inner(T, u)*ds
# First variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(F, u, du)
# Note that ``derivative`` is here used with three arguments: the form
# to be differentiated, the variable (function) we are supposed to
# differentiate with respect too, and the direction the derivative is
# taken in.
#
# Before the form file can be used in the C++ program, it must be
# compiled using FFCX by running (on the command-line):
#
# .. code-block:: sh
#
# ffcx HyperElasticity.ufl
|
