In above formulation, $\varepsilon(\mathbf{u})= \frac 12 (\nabla \mathbf{u} + \nabla
\mathbf{u}^T)$ is the symmetric gradient of the displacement, also called the
<em>strain</em>. $C$ is a tensor of rank 4, called the <em>stress-strain
- tensor</em> that contains knowledge of the elastic strength of the material; its
+ tensor</em> (the inverse of the <a
+ href="https://en.wikipedia.org/wiki/Hooke%27s_law#Hooke's_law_for_continuous_media"><em>compliance
+ tensor</em></a>)
+that contains knowledge of the elastic strength of the material; its
symmetry properties make sure that it maps symmetric tensors of rank 2
(“matrices” of dimension $d$, where $d$ is the spatial dimensionality) onto
symmetric tensors of the same rank. We will comment on the roles of the strain
(\mathbf{f}, \varphi)_{\Omega(t_{n-1})}
-(\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}
\\
- &\qquad +(\mathbf{b}(\mathbf{x},t_n)-\mathbf{b}(\mathbf{x},t_{n-1}), \varphi)_{\Gamma_N}
+ &\qquad
+ +(\mathbf{b}(\mathbf{x},t_n)-\mathbf{b}(\mathbf{x},t_{n-1}), \varphi)_{\Gamma_N}
+ +(\sigma^{n-1} \mathbf{n}, \varphi)_{\Gamma_N}
+ \\
+ &\qquad\qquad
+ \forall \varphi \in \{\mathbf{v}\in H^1(\Omega(t_{n-1}))^d: \mathbf{v}|_{\Gamma_D}=0\}.
+@f}
+Using that $\sigma^{n-1} \mathbf{n}
+ = [C \varepsilon(\mathbf{u}^{n-1})] \mathbf{n}
+ = \mathbf{b}(\mathbf x, t_{n-1})$,
+these equations can be simplified to
+@f{align*}
+ (C \varepsilon(\Delta\mathbf{u}^n), \varepsilon(\varphi) )_{\Omega(t_{n-1})}
+ &=
+ (\mathbf{f}, \varphi)_{\Omega(t_{n-1})}
+ -(\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}
+ +(\mathbf{b}(\mathbf{x},t_n),t_{n-1}), \varphi)_{\Gamma_N}
\\
&\qquad\qquad
\forall \varphi \in \{\mathbf{v}\in H^1(\Omega(t_{n-1}))^d: \mathbf{v}|_{\Gamma_D}=0\}.
\qquad
\textrm{[linear-system]}
@f}
+
We note that, for simplicity, in the program we will always assume that there
are no boundary forces, i.e. $\mathbf{b} = 0$, and that the deformation of the
body is driven by body forces $\mathbf{f}$ and prescribed boundary displacements
$\mathbf{d}$ alone. It is also worth noting that when integrating by parts, we
would get terms of the form $(C \varepsilon(\Delta\mathbf{u}^n), \nabla \varphi
-)_{\Omega(t_{n-1})}$, but that we replace it with the term involving the
+)_{\Omega(t_{n-1})}$, but that we replace them with the term involving the
symmetric gradient $\varepsilon(\varphi)$ instead of $\nabla\varphi$. Due to
-the symmetry of $C$, the two terms are equivalent, but the symmetric version
-avoids a potential for round-off to render the resulting matrix slightly
-non-symmetric.
+the symmetry of $C$, the two terms are mathematically equivalent, but
+the symmetric version avoids the potential for round-off errors making
+the resulting matrix slightly non-symmetric.
The system at time step $n$, to be solved on the old domain
$\Omega(t_{n-1})$, has exactly the form of a stationary elastic