$
\newcommand{\Ua}{\mathrm{a}}
\newcommand{\Ub}{\mathrm{b}}
\newcommand{\Uc}{\mathrm{c}}
\newcommand{\Ud}{\mathrm{d}}
\newcommand{\Ue}{\mathrm{e}}
\newcommand{\Uf}{\mathrm{f}}
\newcommand{\Ug}{\mathrm{g}}
\newcommand{\Uh}{\mathrm{h}}
\newcommand{\Ui}{\mathrm{i}}
\newcommand{\Uj}{\mathrm{j}}
\newcommand{\Uk}{\mathrm{k}}
\newcommand{\Ul}{\mathrm{l}}
\newcommand{\Um}{\mathrm{m}}
\newcommand{\Un}{\mathrm{n}}
\newcommand{\Uo}{\mathrm{o}}
\newcommand{\Up}{\mathrm{p}}
\newcommand{\Uq}{\mathrm{q}}
\newcommand{\Ur}{\mathrm{r}}
\newcommand{\Us}{\mathrm{s}}
\newcommand{\Ut}{\mathrm{t}}
\newcommand{\Uu}{\mathrm{u}}
\newcommand{\Uv}{\mathrm{v}}
\newcommand{\Uw}{\mathrm{w}}
\newcommand{\Ux}{\mathrm{x}}
\newcommand{\Uy}{\mathrm{y}}
\newcommand{\Uz}{\mathrm{z}}
\newcommand{\UA}{\mathrm{A}}
\newcommand{\UB}{\mathrm{B}}
\newcommand{\UC}{\mathrm{C}}
\newcommand{\UD}{\mathrm{D}}
\newcommand{\UE}{\mathrm{E}}
\newcommand{\UF}{\mathrm{F}}
\newcommand{\UG}{\mathrm{G}}
\newcommand{\UH}{\mathrm{H}}
\newcommand{\UI}{\mathrm{I}}
\newcommand{\UJ}{\mathrm{J}}
\newcommand{\UK}{\mathrm{K}}
\newcommand{\UL}{\mathrm{L}}
\newcommand{\UM}{\mathrm{M}}
\newcommand{\UN}{\mathrm{N}}
\newcommand{\UO}{\mathrm{O}}
\newcommand{\UP}{\mathrm{P}}
\newcommand{\UQ}{\mathrm{Q}}
\newcommand{\UR}{\mathrm{R}}
\newcommand{\US}{\mathrm{S}}
\newcommand{\UT}{\mathrm{T}}
\newcommand{\UU}{\mathrm{U}}
\newcommand{\UV}{\mathrm{V}}
\newcommand{\UW}{\mathrm{W}}
\newcommand{\UX}{\mathrm{X}}
\newcommand{\UY}{\mathrm{Y}}
\newcommand{\UZ}{\mathrm{Z}}
%
\newcommand{\Uzero }{\mathrm{0}}
\newcommand{\Uone }{\mathrm{1}}
\newcommand{\Utwo }{\mathrm{2}}
\newcommand{\Uthree}{\mathrm{3}}
\newcommand{\Ufour }{\mathrm{4}}
\newcommand{\Ufive }{\mathrm{5}}
\newcommand{\Usix }{\mathrm{6}}
\newcommand{\Useven}{\mathrm{7}}
\newcommand{\Ueight}{\mathrm{8}}
\newcommand{\Unine }{\mathrm{9}}
%
\newcommand{\Ja}{\mathit{a}}
\newcommand{\Jb}{\mathit{b}}
\newcommand{\Jc}{\mathit{c}}
\newcommand{\Jd}{\mathit{d}}
\newcommand{\Je}{\mathit{e}}
\newcommand{\Jf}{\mathit{f}}
\newcommand{\Jg}{\mathit{g}}
\newcommand{\Jh}{\mathit{h}}
\newcommand{\Ji}{\mathit{i}}
\newcommand{\Jj}{\mathit{j}}
\newcommand{\Jk}{\mathit{k}}
\newcommand{\Jl}{\mathit{l}}
\newcommand{\Jm}{\mathit{m}}
\newcommand{\Jn}{\mathit{n}}
\newcommand{\Jo}{\mathit{o}}
\newcommand{\Jp}{\mathit{p}}
\newcommand{\Jq}{\mathit{q}}
\newcommand{\Jr}{\mathit{r}}
\newcommand{\Js}{\mathit{s}}
\newcommand{\Jt}{\mathit{t}}
\newcommand{\Ju}{\mathit{u}}
\newcommand{\Jv}{\mathit{v}}
\newcommand{\Jw}{\mathit{w}}
\newcommand{\Jx}{\mathit{x}}
\newcommand{\Jy}{\mathit{y}}
\newcommand{\Jz}{\mathit{z}}
\newcommand{\JA}{\mathit{A}}
\newcommand{\JB}{\mathit{B}}
\newcommand{\JC}{\mathit{C}}
\newcommand{\JD}{\mathit{D}}
\newcommand{\JE}{\mathit{E}}
\newcommand{\JF}{\mathit{F}}
\newcommand{\JG}{\mathit{G}}
\newcommand{\JH}{\mathit{H}}
\newcommand{\JI}{\mathit{I}}
\newcommand{\JJ}{\mathit{J}}
\newcommand{\JK}{\mathit{K}}
\newcommand{\JL}{\mathit{L}}
\newcommand{\JM}{\mathit{M}}
\newcommand{\JN}{\mathit{N}}
\newcommand{\JO}{\mathit{O}}
\newcommand{\JP}{\mathit{P}}
\newcommand{\JQ}{\mathit{Q}}
\newcommand{\JR}{\mathit{R}}
\newcommand{\JS}{\mathit{S}}
\newcommand{\JT}{\mathit{T}}
\newcommand{\JU}{\mathit{U}}
\newcommand{\JV}{\mathit{V}}
\newcommand{\JW}{\mathit{W}}
\newcommand{\JX}{\mathit{X}}
\newcommand{\JY}{\mathit{Y}}
\newcommand{\JZ}{\mathit{Z}}
%
\newcommand{\Jzero }{\mathit{0}}
\newcommand{\Jone }{\mathit{1}}
\newcommand{\Jtwo }{\mathit{2}}
\newcommand{\Jthree}{\mathit{3}}
\newcommand{\Jfour }{\mathit{4}}
\newcommand{\Jfive }{\mathit{5}}
\newcommand{\Jsix }{\mathit{6}}
\newcommand{\Jseven}{\mathit{7}}
\newcommand{\Jeight}{\mathit{8}}
\newcommand{\Jnine }{\mathit{9}}
%
\newcommand{\BA}{\boldsymbol{A}}
\newcommand{\BB}{\boldsymbol{B}}
\newcommand{\BC}{\boldsymbol{C}}
\newcommand{\BD}{\boldsymbol{D}}
\newcommand{\BE}{\boldsymbol{E}}
\newcommand{\BF}{\boldsymbol{F}}
\newcommand{\BG}{\boldsymbol{G}}
\newcommand{\BH}{\boldsymbol{H}}
\newcommand{\BI}{\boldsymbol{I}}
\newcommand{\BJ}{\boldsymbol{J}}
\newcommand{\BK}{\boldsymbol{K}}
\newcommand{\BL}{\boldsymbol{L}}
\newcommand{\BM}{\boldsymbol{M}}
\newcommand{\BN}{\boldsymbol{N}}
\newcommand{\BO}{\boldsymbol{O}}
\newcommand{\BP}{\boldsymbol{P}}
\newcommand{\BQ}{\boldsymbol{Q}}
\newcommand{\BR}{\boldsymbol{R}}
\newcommand{\BS}{\boldsymbol{S}}
\newcommand{\BT}{\boldsymbol{T}}
\newcommand{\BU}{\boldsymbol{U}}
\newcommand{\BV}{\boldsymbol{V}}
\newcommand{\BW}{\boldsymbol{W}}
\newcommand{\BX}{\boldsymbol{X}}
\newcommand{\BY}{\boldsymbol{Y}}
\newcommand{\BZ}{\boldsymbol{Z}}
\newcommand{\Ba}{\boldsymbol{a}}
\newcommand{\Bb}{\boldsymbol{b}}
\newcommand{\Bc}{\boldsymbol{c}}
\newcommand{\Bd}{\boldsymbol{d}}
\newcommand{\Be}{\boldsymbol{e}}
\newcommand{\Bf}{\boldsymbol{f}}
\newcommand{\Bg}{\boldsymbol{g}}
\newcommand{\Bh}{\boldsymbol{h}}
\newcommand{\Bi}{\boldsymbol{i}}
\newcommand{\Bj}{\boldsymbol{j}}
\newcommand{\Bk}{\boldsymbol{k}}
\newcommand{\Bl}{\boldsymbol{l}}
\newcommand{\Bm}{\boldsymbol{m}}
\newcommand{\Bn}{\boldsymbol{n}}
\newcommand{\Bo}{\boldsymbol{o}}
\newcommand{\Bp}{\boldsymbol{p}}
\newcommand{\Bq}{\boldsymbol{q}}
\newcommand{\Br}{\boldsymbol{r}}
\newcommand{\Bs}{\boldsymbol{s}}
\newcommand{\Bt}{\boldsymbol{t}}
\newcommand{\Bu}{\boldsymbol{u}}
\newcommand{\Bv}{\boldsymbol{v}}
\newcommand{\Bw}{\boldsymbol{w}}
\newcommand{\Bx}{\boldsymbol{x}}
\newcommand{\By}{\boldsymbol{y}}
\newcommand{\Bz}{\boldsymbol{z}}
%
\newcommand{\Bzero }{\boldsymbol{0}}
\newcommand{\Bone }{\boldsymbol{1}}
\newcommand{\Btwo }{\boldsymbol{2}}
\newcommand{\Bthree}{\boldsymbol{3}}
\newcommand{\Bfour }{\boldsymbol{4}}
\newcommand{\Bfive }{\boldsymbol{5}}
\newcommand{\Bsix }{\boldsymbol{6}}
\newcommand{\Bseven}{\boldsymbol{7}}
\newcommand{\Beight}{\boldsymbol{8}}
\newcommand{\Bnine }{\boldsymbol{9}}
%
\newcommand{\Balpha }{\boldsymbol{\alpha} }
\newcommand{\Bbeta }{\boldsymbol{\beta} }
\newcommand{\Bgamma }{\boldsymbol{\gamma} }
\newcommand{\Bdelta }{\boldsymbol{\delta} }
\newcommand{\Bepsilon}{\boldsymbol{\epsilon} }
\newcommand{\Bvareps }{\boldsymbol{\varepsilon} }
\newcommand{\Bvarepsilon}{\boldsymbol{\varepsilon}}
\newcommand{\Bzeta }{\boldsymbol{\zeta} }
\newcommand{\Beta }{\boldsymbol{\eta} }
\newcommand{\Btheta }{\boldsymbol{\theta} }
\newcommand{\Bvarthe }{\boldsymbol{\vartheta} }
\newcommand{\Biota }{\boldsymbol{\iota} }
\newcommand{\Bkappa }{\boldsymbol{\kappa} }
\newcommand{\Blambda }{\boldsymbol{\lambda} }
\newcommand{\Bmu }{\boldsymbol{\mu} }
\newcommand{\Bnu }{\boldsymbol{\nu} }
\newcommand{\Bxi }{\boldsymbol{\xi} }
\newcommand{\Bpi }{\boldsymbol{\pi} }
\newcommand{\Brho }{\boldsymbol{\rho} }
\newcommand{\Bvrho }{\boldsymbol{\varrho} }
\newcommand{\Bsigma }{\boldsymbol{\sigma} }
\newcommand{\Bvsigma }{\boldsymbol{\varsigma} }
\newcommand{\Btau }{\boldsymbol{\tau} }
\newcommand{\Bupsilon}{\boldsymbol{\upsilon} }
\newcommand{\Bphi }{\boldsymbol{\phi} }
\newcommand{\Bvarphi }{\boldsymbol{\varphi} }
\newcommand{\Bchi }{\boldsymbol{\chi} }
\newcommand{\Bpsi }{\boldsymbol{\psi} }
\newcommand{\Bomega }{\boldsymbol{\omega} }
\newcommand{\BGamma }{\boldsymbol{\Gamma} }
\newcommand{\BDelta }{\boldsymbol{\Delta} }
\newcommand{\BTheta }{\boldsymbol{\Theta} }
\newcommand{\BLambda }{\boldsymbol{\Lambda} }
\newcommand{\BXi }{\boldsymbol{\Xi} }
\newcommand{\BPi }{\boldsymbol{\Pi} }
\newcommand{\BSigma }{\boldsymbol{\Sigma} }
\newcommand{\BUpsilon}{\boldsymbol{\Upsilon} }
\newcommand{\BPhi }{\boldsymbol{\Phi} }
\newcommand{\BPsi }{\boldsymbol{\Psi} }
\newcommand{\BOmega }{\boldsymbol{\Omega} }
%
\newcommand{\IA}{\mathbb{A}}
\newcommand{\IB}{\mathbb{B}}
\newcommand{\IC}{\mathbb{C}}
\newcommand{\ID}{\mathbb{D}}
\newcommand{\IE}{\mathbb{E}}
\newcommand{\IF}{\mathbb{F}}
\newcommand{\IG}{\mathbb{G}}
\newcommand{\IH}{\mathbb{H}}
\newcommand{\II}{\mathbb{I}}
\renewcommand{\IJ}{\mathbb{J}}
\newcommand{\IK}{\mathbb{K}}
\newcommand{\IL}{\mathbb{L}}
\newcommand{\IM}{\mathbb{M}}
\newcommand{\IN}{\mathbb{N}}
\newcommand{\IO}{\mathbb{O}}
\newcommand{\IP}{\mathbb{P}}
\newcommand{\IQ}{\mathbb{Q}}
\newcommand{\IR}{\mathbb{R}}
\newcommand{\IS}{\mathbb{S}}
\newcommand{\IT}{\mathbb{T}}
\newcommand{\IU}{\mathbb{U}}
\newcommand{\IV}{\mathbb{V}}
\newcommand{\IW}{\mathbb{W}}
\newcommand{\IX}{\mathbb{X}}
\newcommand{\IY}{\mathbb{Y}}
\newcommand{\IZ}{\mathbb{Z}}
%
\newcommand{\FA}{\mathsf{A}}
\newcommand{\FB}{\mathsf{B}}
\newcommand{\FC}{\mathsf{C}}
\newcommand{\FD}{\mathsf{D}}
\newcommand{\FE}{\mathsf{E}}
\newcommand{\FF}{\mathsf{F}}
\newcommand{\FG}{\mathsf{G}}
\newcommand{\FH}{\mathsf{H}}
\newcommand{\FI}{\mathsf{I}}
\newcommand{\FJ}{\mathsf{J}}
\newcommand{\FK}{\mathsf{K}}
\newcommand{\FL}{\mathsf{L}}
\newcommand{\FM}{\mathsf{M}}
\newcommand{\FN}{\mathsf{N}}
\newcommand{\FO}{\mathsf{O}}
\newcommand{\FP}{\mathsf{P}}
\newcommand{\FQ}{\mathsf{Q}}
\newcommand{\FR}{\mathsf{R}}
\newcommand{\FS}{\mathsf{S}}
\newcommand{\FT}{\mathsf{T}}
\newcommand{\FU}{\mathsf{U}}
\newcommand{\FV}{\mathsf{V}}
\newcommand{\FW}{\mathsf{W}}
\newcommand{\FX}{\mathsf{X}}
\newcommand{\FY}{\mathsf{Y}}
\newcommand{\FZ}{\mathsf{Z}}
\newcommand{\Fa}{\mathsf{a}}
\newcommand{\Fb}{\mathsf{b}}
\newcommand{\Fc}{\mathsf{c}}
\newcommand{\Fd}{\mathsf{d}}
\newcommand{\Fe}{\mathsf{e}}
\newcommand{\Ff}{\mathsf{f}}
\newcommand{\Fg}{\mathsf{g}}
\newcommand{\Fh}{\mathsf{h}}
\newcommand{\Fi}{\mathsf{i}}
\newcommand{\Fj}{\mathsf{j}}
\newcommand{\Fk}{\mathsf{k}}
\newcommand{\Fl}{\mathsf{l}}
\newcommand{\Fm}{\mathsf{m}}
\newcommand{\Fn}{\mathsf{n}}
\newcommand{\Fo}{\mathsf{o}}
\newcommand{\Fp}{\mathsf{p}}
\newcommand{\Fq}{\mathsf{q}}
\newcommand{\Fr}{\mathsf{r}}
\newcommand{\Fs}{\mathsf{s}}
\newcommand{\Ft}{\mathsf{t}}
\newcommand{\Fu}{\mathsf{u}}
\newcommand{\Fv}{\mathsf{v}}
\newcommand{\Fw}{\mathsf{w}}
\newcommand{\Fx}{\mathsf{x}}
\newcommand{\Fy}{\mathsf{y}}
\newcommand{\Fz}{\mathsf{z}}
%
\newcommand{\Fzero }{\mathsf{0}}
\newcommand{\Fone }{\mathsf{1}}
\newcommand{\Ftwo }{\mathsf{2}}
\newcommand{\Fthree}{\mathsf{3}}
\newcommand{\Ffour }{\mathsf{4}}
\newcommand{\Ffive }{\mathsf{5}}
\newcommand{\Fsix }{\mathsf{6}}
\newcommand{\Fseven}{\mathsf{7}}
\newcommand{\Feight}{\mathsf{8}}
\newcommand{\Fnine }{\mathsf{9}}
%
\newcommand{\CA}{\mathcal{A}}
\newcommand{\CB}{\mathcal{B}}
\newcommand{\CC}{\mathcal{C}}
\newcommand{\CD}{\mathcal{D}}
\newcommand{\CE}{\mathcal{E}}
\newcommand{\CF}{\mathcal{F}}
\newcommand{\CG}{\mathcal{G}}
\newcommand{\CH}{\mathcal{H}}
\newcommand{\CI}{\mathcal{I}}
\newcommand{\CJ}{\mathcal{J}}
\newcommand{\CK}{\mathcal{K}}
\newcommand{\CL}{\mathcal{L}}
\newcommand{\CM}{\mathcal{M}}
\newcommand{\CN}{\mathcal{N}}
\newcommand{\CO}{\mathcal{O}}
\newcommand{\CP}{\mathcal{P}}
\newcommand{\CQ}{\mathcal{Q}}
\newcommand{\CR}{\mathcal{R}}
\newcommand{\CS}{\mathcal{S}}
\newcommand{\CT}{\mathcal{T}}
\newcommand{\CU}{\mathcal{U}}
\newcommand{\CV}{\mathcal{V}}
\newcommand{\CW}{\mathcal{W}}
\newcommand{\CX}{\mathcal{X}}
\newcommand{\CY}{\mathcal{Y}}
\newcommand{\CZ}{\mathcal{Z}}
%
\newcommand{\KA}{\mathfrak{A}}
\newcommand{\KB}{\mathfrak{B}}
\newcommand{\KC}{\mathfrak{C}}
\newcommand{\KD}{\mathfrak{D}}
\newcommand{\KE}{\mathfrak{E}}
\newcommand{\KF}{\mathfrak{F}}
\newcommand{\KG}{\mathfrak{G}}
\newcommand{\KH}{\mathfrak{H}}
\newcommand{\KI}{\mathfrak{I}}
\newcommand{\KJ}{\mathfrak{J}}
\newcommand{\KK}{\mathfrak{K}}
\newcommand{\KL}{\mathfrak{L}}
\newcommand{\KM}{\mathfrak{M}}
\newcommand{\KN}{\mathfrak{N}}
\newcommand{\KO}{\mathfrak{O}}
\newcommand{\KP}{\mathfrak{P}}
\newcommand{\KQ}{\mathfrak{Q}}
\newcommand{\KR}{\mathfrak{R}}
\newcommand{\KS}{\mathfrak{S}}
\newcommand{\KT}{\mathfrak{T}}
\newcommand{\KU}{\mathfrak{U}}
\newcommand{\KV}{\mathfrak{V}}
\newcommand{\KW}{\mathfrak{W}}
\newcommand{\KX}{\mathfrak{X}}
\newcommand{\KY}{\mathfrak{Y}}
\newcommand{\KZ}{\mathfrak{Z}}
\newcommand{\Ka}{\mathfrak{a}}
\newcommand{\Kb}{\mathfrak{b}}
\newcommand{\Kc}{\mathfrak{c}}
\newcommand{\Kd}{\mathfrak{d}}
\newcommand{\Ke}{\mathfrak{e}}
\newcommand{\Kf}{\mathfrak{f}}
\newcommand{\Kg}{\mathfrak{g}}
\newcommand{\Kh}{\mathfrak{h}}
\newcommand{\Ki}{\mathfrak{i}}
\newcommand{\Kj}{\mathfrak{j}}
\newcommand{\Kk}{\mathfrak{k}}
\newcommand{\Kl}{\mathfrak{l}}
\newcommand{\Km}{\mathfrak{m}}
\newcommand{\Kn}{\mathfrak{n}}
\newcommand{\Ko}{\mathfrak{o}}
\newcommand{\Kp}{\mathfrak{p}}
\newcommand{\Kq}{\mathfrak{q}}
\newcommand{\Kr}{\mathfrak{r}}
\newcommand{\Ks}{\mathfrak{s}}
\newcommand{\Kt}{\mathfrak{t}}
\newcommand{\Ku}{\mathfrak{u}}
\newcommand{\Kv}{\mathfrak{v}}
\newcommand{\Kw}{\mathfrak{w}}
\newcommand{\Kx}{\mathfrak{x}}
\newcommand{\Ky}{\mathfrak{y}}
\newcommand{\Kz}{\mathfrak{z}}
%
\newcommand{\Kzero }{\mathfrak{0}}
\newcommand{\Kone }{\mathfrak{1}}
\newcommand{\Ktwo }{\mathfrak{2}}
\newcommand{\Kthree}{\mathfrak{3}}
\newcommand{\Kfour }{\mathfrak{4}}
\newcommand{\Kfive }{\mathfrak{5}}
\newcommand{\Ksix }{\mathfrak{6}}
\newcommand{\Kseven}{\mathfrak{7}}
\newcommand{\Keight}{\mathfrak{8}}
\newcommand{\Knine }{\mathfrak{9}}
%
$
$
\newcommand{\Lin}{\mathop{\rm Lin}\nolimits}
\newcommand{\modop}{\mathop{\rm mod}\nolimits}
\renewcommand{\div}{\mathop{\rm div}\nolimits}
\newcommand{\Var}{\Delta}
\newcommand{\evat}{\bigg|}
\newcommand\varn[3]{D_{#2}#1\cdot #3}
\newcommand{\dtp}{\cdot}
\newcommand{\dyd}{\otimes}
\newcommand{\tra}{^T}
\newcommand{\del}{\partial}
\newcommand{\dif}{d}
\newcommand{\rbr}[1]{\left(#1\right)}
\newcommand{\sbr}[1]{\left[#1\right]}
\newcommand{\cbr}[1]{\left\{#1\right\}}
\newcommand{\cbrn}[1]{\{#1\}}
\newcommand{\abr}[1]{\left\langle #1 \right\rangle}
\newcommand{\abrn}[1]{\langle #1 \rangle}
\newcommand{\deriv}[2]{\frac{d #1}{d #2}}
\newcommand{\dderiv}[2]{\frac{d^2 #1}{d {#2}^2}}
\newcommand{\partd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\nnode}{n_n}
\newcommand{\ndim}{n_d}
\newcommand{\suml}[2]{\sum\limits_{#1}^{#2}}
\newcommand{\Aelid}[2]{A^{#1}_{#2}}
\newcommand{\dv}{\, dv}
\newcommand{\dx}{\, dx}
\newcommand{\ds}{\, ds}
\newcommand{\da}{\, da}
\newcommand{\dV}{\, dV}
\newcommand{\dA}{\, dA}
\newcommand{\eqand}{\quad\text{and}\quad}
\newcommand{\eqor}{\quad\text{or}\quad}
\newcommand{\eqwith}{\quad\text{and}\quad}
\newcommand{\inv}{^{-1}}
\newcommand{\veci}[1]{#1_1,\ldots,#1_n}
\newcommand{\var}{\delta}
\newcommand{\Var}{\Delta}
\newcommand{\eps}{\epsilon}
\newcommand{\ddt}{\frac{d}{dt}}
\newcommand{\Norm}[1]{\left\lVert#1\right\rVert}
\newcommand{\Abs}[1]{\left|#1\right|}
\newcommand{\dabr}[1]{\left\langle\!\left\langle #1 \right\rangle\!\right\rangle}
\newcommand{\dabrn}[1]{\langle\!\langle #1 \rangle\!\rangle}
\newcommand{\idxsep}{\,}
$
$
\newcommand{\argmin}{\mathop{\rm argmin}\nolimits}
\newcommand{\cof}{\mathop{\rm cof}\nolimits}
\newcommand{\sym}{\mathop{\rm sym}\nolimits}
\newcommand{\invtra}{^{-T}}
\newcommand{\eps}{\epsilon}
\newcommand{\var}{\Delta}
\newcommand{\Vvphi}{\Delta\Bvarphi}
\newcommand{\vvphi}{\delta\Bvarphi}
\newcommand{\BFC}{\boldsymbol{\mathsf{C}}}
\newcommand{\BFc}{\boldsymbol{\mathsf{c}}}
\newcommand{\push}{\Bvarphi_\ast}
\newcommand{\pull}{\Bvarphi^\ast}
$
There are many books that give an outline of hyperelasticity, but there are few
that try to help the reader implement solutions, and even fewer that
manage to do it in a concise manner. Peter Wriggers’
Nonlinear Finite Element Methods
is a great reference for those who like to
roll up their sleeves and get lost in theory. It helped me understand a lot
about how solutions to hyperelastic and inelastic problems are implemented.
One thing did not quite fit my taste though—it was very formal in the way that
it didn’t give out indicial expressions. And if it wasn’t clear up until this
point, I love indicial expressions, because they pack enough information to
implement a solution in a single line. Almost all books skip these
because they seem cluttered and the professors who wrote them think they’re
trivial to derive. In fact, they are not.
So below, I’ll try to derive indicial expressions for the update equations of
hyperelasticity.
In the case of a hyperelastic material, there exists a strain energy function
\[\begin{equation}
\Psi: \BF \mapsto \Psi(\BF)
\end{equation}\]
which describes the elastic energy stored in the solid, i.e. energy
density per unit mass of the reference configuration.
The total energy stored in $\CB$ is described by the the stored energy
functional
\[\begin{equation}
E(\Bvarphi) := \int_\CB \Psi(\BF)\, \dif m = \int_\CB \rho_0 \Psi(\BF) \dV
\end{equation}\]
The loads acting on the body also form a potential:
\[\begin{equation}
L(\Bvarphi) := \int_\CB \rho_0\bar{\BGamma}\dtp\Bvarphi \dV + \int_{\del\CB_t} \bar{\BT}\dtp\Bvarphi \dA
\end{equation}\]
where $\bar{\BGamma}$ and $\bar{\BT}$ are the prescribed body forces per unit mass and surface tractions
respectively, where $\BT=\BP\BN$ with Cauchy’s stress theorem.
The potential energy of $\CB$ for deformation $\Bvarphi$ is defined as
\[\begin{equation}
\Pi(\Bvarphi) := E(\Bvarphi) - L(\Bvarphi)
\end{equation}\]
Thus the variational formulation reads
Find $\Bvarphi\in V$ such that the functional
\[\begin{equation}
\Pi(\Bvarphi)
= \int_\CB \rho_0\Psi(\BF) \dV
- \int_\CB \rho_0\bar{\BGamma}\dtp\Bvarphi \dV - \int_{\del\CB_t} \bar{\BT}\dtp\Bvarphi \dA
\end{equation}\]
is minimized for $\Bvarphi=\bar{\Bvarphi}$ on $\del\CB_u$.
The solution is one that minimizes the potential energy:
\[\begin{equation}
\Bvarphi^\ast = \argmin_{\Bvarphi\in V} \Pi(\Bvarphi)
\end{equation}\]
A stationary point for $\Pi$ means that its first variation vanishes: $\var\Pi=0$.
\[\begin{equation}
\begin{aligned}
\var\Pi &= \varn{\Pi}{\Bvarphi}{\vvphi}
=: G(\Bvarphi,\vvphi) \\
&= \int_\CB \rho_0\partd{\Psi}{\BF}: \nabla(\vvphi) \dV
- \int_\CB \rho_0\bar{\BGamma}\dtp\vvphi \dV
- \int_{\del\CB} \bar{\BT}\dtp\vvphi \dA
\end{aligned}
\end{equation}\]
Using $\BP=\BF\BS$ and $\BP = \rho_0\del\Psi/\del\BF$,
\[\begin{equation}
\rho_0\partd{\Psi}{\BF}: \nabla(\vvphi)
= \BF\BS:\nabla(\vvphi)
= \BS:\BF\tra\nabla(\vvphi)
\end{equation}\]
The symmetric part of the term on the right hand side of the contraction is
equal to the variation of the Green-Lagrange strain tensor:
\[\begin{equation}
\begin{aligned}
\var\BE = \varn{\BE}{\Bvarphi}{\vvphi}
&= \deriv{}{\eps} \frac{1}{2}
[\nabla(\Bvarphi+\eps\vvphi)\tra\nabla(\Bvarphi+\eps\vvphi) - \BI]\evat_{\eps=0}
\\
&= \frac{1}{2} [\nabla(\vvphi)\tra\BF + \BF\tra\nabla(\vvphi)]
\end{aligned}
\end{equation}\]
Substituting, we obtain the semilinear form $G$ in terms of the
second Piola-Kirchhoff stress tensor:
\[\begin{equation}
\boxed{
G(\Bvarphi,\vvphi)
= \int_\CB \BS: \var\BE \dV
- \int_\CB \rho_0\bar{\BGamma}\dtp\vvphi \dV
- \int_{\del\CB} \bar{\BT}\dtp\vvphi \dA = 0
}
\label{eq:lagrangianform1}
\end{equation}\]
We can write a Eulerian version of this form by pushing-forward the stresses and
strains.
The Almansi strain $\Be$ is the pull-back of the Green-Lagrange strain $\BE$ and
vice versa:
\[\begin{equation}
\Be = \push(\BE) = \BF\invtra \BE \BF\inv
\eqand
\BE = \pull(\Be) = \BF\tra \BE \BF
\end{equation}\]
Commutative diagram for the pull-back and push-forward relationships of the Green-Lagrange and
Almansi strain tensors.
Thus we can deduce the variation of the Almansi strain
\[\begin{equation}
\begin{aligned}
\var \Be
= \BF\invtra \var\BE\BF\inv
&= \frac{1}{2} [\nabla(\vvphi)\BF\inv+\BF\invtra \nabla(\vvphi)\tra] \\
&= \frac{1}{2} [\nabla_x(\vvphi)+ \nabla_x(\vvphi)\tra]
\end{aligned}
\end{equation}\]
where we have used the identity
\[\begin{equation}
\nabla_X(\cdot)\BF\inv = \nabla_x(\cdot).
\label{eq:defgradidentity1}
\end{equation}\]
The second Piola-Kirchhoff stress is the
pull-back of the Kirchhoff stress $\Btau$:
\[\begin{equation}
\BS = \pull(\Btau) = \BF\inv\Btau\BF\invtra
\end{equation}\]
Then it is evident that
\[\begin{equation}
\BS:\var\BE
= (\BF\inv\Btau\BF\invtra):(\BF\tra\var\Be\BF)
= \Btau:\var\Be
\end{equation}\]
We can thus write the Eulerian version of \eqref{eq:lagrangianform1}:
\[\begin{equation}
\boxed{
G(\Bvarphi,\vvphi)
= \int_\CB \Btau: \var\Be \dV
- \int_\CB \rho_0\bar{\BGamma}\dtp\vvphi \dV
- \int_{\del\CB} \bar{\BT}\dtp\vvphi \dA = 0
}
\end{equation}\]
Introducing the Cauchy stresses $\Bsigma=\Btau/J$, we can also transport the
integrals to the current configuration
\[\begin{equation}
\boxed{
G(\Bvarphi,\vvphi)
= \int_\CS \Bsigma:\var\Be \dv
- \int_\CS \rho\bar{\Bgamma}\dtp\vvphi \dv
- \int_{\del\CS_t} \bar{\Bt}\dtp\vvphi \da = 0
}
\end{equation}\]
Here, we substituted the following differential identities:
\[\begin{equation}
\rho_0\BGamma\dV = \rho\Bgamma\dv
\end{equation}\]
for the body forces, and
\[\begin{equation}
\BT\dA
= \BP\BN \dA
= \Bsigma J\BF\invtra\BN \dA
= \Bsigma\Bn \da
= \Bt \da
\end{equation}\]
for the surface tractions, where we used the Piola identity.
We linearize $G$:
\[\begin{equation}
\Lin G \evat_{\bar{\Bvarphi}}
= G(\bar{\Bvarphi}, \vvphi)
+ \Var G \evat_{\bar{\Bvarphi}}
= 0
\end{equation}\]
Then we have the variational setting
\[\begin{equation}
a(\Vvphi,\vvphi)=b(\vvphi)
\end{equation}\]
where
\[\begin{equation}
a(\Vvphi,\vvphi) = \Var G \evat_{\bar{\Bvarphi}}
\eqand
b(\vvphi) = -G(\bar{\Bvarphi}, \vvphi)
\end{equation}\]
Commutative diagram of the linearized solution procedure. Each
iteration brings the current iterate $\bar{\Bvarphi}$ closer to the optimum
value $\Bvarphi^\ast$.
Mappings between line elements belonging to the tangent spaces of
the linearization.
The variation $\Var G$ is calculated as
\[\begin{equation}
\Var G
= \varn{G}{\Bvarphi}{\Vvphi}
= \int_\CB [\Var\BS:\var\BE + \BS:\Var(\var\BE)] \dV
\end{equation}\]
Consecutive variations of the Green-Lagrange strain tensor is calculated as
\[\begin{equation}
\Var(\var\BE) = \varn{\var\BE}{\Bvarphi}{\Vvphi}
= \frac{1}{2}[\nabla(\vvphi)\tra\nabla(\Vvphi) + \nabla(\Vvphi)\tra\nabla(\vvphi)]
\end{equation}\]
The term on the left is calculated as
\[\begin{equation}
\Var\BS = \varn{\BS}{\Bvarphi}{\Vvphi}
= \partd{\BS}{\BC}:\Var\BC
= 2 \partd{\BS}{\BC}:\Var\BE
\end{equation}\]
where we substitute the Lagrangian elasticity tensor
\[\begin{equation}
\BFC := 2 \partd{\BS}{\BC} = 4\rho_0 \frac{\del^2\Psi}{\del\BC\del\BC}
\end{equation}\]
and $\Var\BE$ is calculated in the same manner as $\var\BE$:
\[\begin{equation}
\Var\BE = \frac{1}{2} [\nabla(\Vvphi)\tra\BF + \BF\tra\nabla(\Vvphi)]
\end{equation}\]
Then the variational forms of the linearized setting are
\[\begin{equation}
\boxed{
\begin{aligned}
a(\Vvphi,\vvphi)
&= \int_\CB \var\bar{\BE}:\bar{\BFC}:\Var\bar{\BE}
+ \bar{\BS} : [\nabla(\vvphi)\tra\nabla(\Vvphi)] \dV
\\
b(\vvphi)
&= - \int_\CB \bar{\BS}: \var\bar{\BE} \dV
+ \int_\CB \rho_0\bar{\BGamma}\dtp\vvphi \dV
+ \int_{\del\CB} \bar{\BT}\dtp\vvphi \dA
\end{aligned}
}
\end{equation}\]
where the bars denote evaluation $\Bvarphi=\bar{\Bvarphi}$ of dependent
variables.
Eulerian Version of the Linearization
We also have the following relationship between the Lagrangian and Eulerian
elasticity tensors
\[\begin{equation}
c_{abcd} = F_{aA}F_{bB}F_{cC}F_{dD} C_{ABCD}
\end{equation}\]
Substituting Eulerian expansions, we obtain the
following identity:
\[\begin{equation}
\begin{aligned}
\var\BE:\BFC:\Var\BE
&= (\BF\tra\var\Be\BF):\BFC:(\BF\tra\Var\Be\BF) \\
&=F_{aA}\var e_{ab} F_{bB} C_{ABCD} F_{cC}\Var e_{cd}F_{dD} \\
&=\var e_{ab} c_{abcd} \Var e_{cd} \\
&= \var\Be:\BFc:\Var\Be
\end{aligned}
\end{equation}\]
Thus we have
\[\begin{equation}
\begin{aligned}
\BS:[\nabla(\vvphi)\tra\nabla(\Vvphi)]
&= [\BF\inv\Btau\BF\invtra] :[\nabla(\vvphi)\tra\nabla(\Vvphi)] \\
&= \Btau : [(\nabla(\vvphi)\BF\inv)\tra\nabla(\Vvphi)\BF\inv] \\
&= \Btau : [\nabla_x(\vvphi)\tra\nabla_x(\Vvphi)] \\
\end{aligned}
\end{equation}\]
With these results at hand, we can write the Eulerian version of our variational
formulation:
\[\begin{equation}
\boxed{
\begin{aligned}
a(\Vvphi,\vvphi)
&= \int_\CB \var\bar{\Be}:\bar{\BFc}:\Var\bar{\Be}
+ \bar{\Btau} : [\nabla_{\bar{x}}(\vvphi)\tra\nabla_{\bar{x}}(\Vvphi)] \dV
\\
b(\vvphi)
&= - \int_\CB
\bar{\Btau}:\var\bar{\Be} \dV
+ \int_\CB \rho_0\bar{\BGamma}\dtp\vvphi \dV
+ \int_{\del\CB} \bar{\BT}\dtp\vvphi \dA
\end{aligned}
}
\end{equation}\]
If we introduce the Cauchy stress tensor $\Bsigma$ and corresponding elasticity tensor
$\BFc^\sigma = \BFc/J$,
our variational formulation can be expressed completely in terms of Eulerian quantities:
\[\begin{equation}
\boxed{
\begin{aligned}
a(\Vvphi,\vvphi)
&= \int_{\bar{\CS}} \var\bar{\Be}:\bar{\BFc}^\sigma:\Var\bar{\Be}
+ \bar{\Bsigma} : [\nabla_{\bar{x}}(\vvphi)\tra\nabla_{\bar{x}}(\Vvphi)] \,\dif\bar{v}
\\
b(\vvphi)
&= - \int_{\bar{\CS}}
\bar{\Bsigma}:\var\bar{\Be} \,\dif\bar{v}
+ \int_{\bar{\CS}} \rho\bar{\Bgamma}\dtp\vvphi \,\dif\bar{v}
+ \int_{\del\bar{\CS}_t} \bar{\Bt}\dtp\vvphi \,\dif\bar{a}
\end{aligned}
}
\end{equation}\]
We have the following relationships of the differential forms:
\[\begin{equation}
\dif \bar{v} = \bar{J}\dv
\eqand
\bar{\Bn} \,\dif \bar{a} = \cof \bar{\BF}\BN \dA
\end{equation}\]
where $\bar{\BF} = \nabla_X\bar{\Bvarphi}$ and $\bar{J} = \det\bar{\BF}$.
We use the following FE discretization:
\[\begin{equation}
\Bvarphi_h
= \suml{\gamma=1}{\nnode} \Bvarphi^\gamma N^\gamma
= \suml{\gamma=1}{\nnode}\suml{a=1}{\ndim} \varphi_a^\gamma \Be_a N^\gamma
\end{equation}\]
where $\nnode$ is the number of element nodes and $\ndim$ is the number of
spatial dimensions.
We use the same discretization for $\vvphi$ and $\Vvphi$. Then the linear system
at hand becomes
\[\begin{equation}
\suml{\delta=1}{\nnode}\suml{b=1}{\ndim}A_{ab}^{\gamma\delta} \Var\varphi_b^\delta = b_a^\gamma
\end{equation}\]
for $a=1,\dots,\ndim$ and $\gamma=1,\dots,\nnode$ where the $\BA$ and $\Bb$
are calculated from the variational forms as
\[\begin{equation}
\begin{aligned}
A_{ab}^{\gamma\delta} &= a(\Be_bN^\delta, \Be_aN^\gamma) \\
b_a^\gamma &= b(\Be_aN^\gamma)
\end{aligned}
\end{equation}\]
For detailed derivation, see the post
Vectorial Finite Elements .
For discretized gradients, we have the following relationship
\[\begin{equation}
\nabla_X(\Be_aN^\gamma) = (\Be_a\dyd\BB^\gamma)
\end{equation}\]
where $\BB^\gamma:= \nabla_X N^\gamma$. For the first term in $a$, we can get rid of the symmetries:
\[\begin{equation}
\begin{aligned}
\sym&(\bar{\BF}\tra\nabla(\Be_aN^\gamma)):\bar{\BFC}:
\sym(\bar{\BF}\tra\nabla(\Be_bN^\delta)) \\
&= (\bar{\BF}\tra(\Be_a\dyd\BB^\gamma)):\bar{\BFC}:
(\bar{\BF}\tra(\Be_b\dyd\BB^\delta)) \\
&= \bar{F}_{aA}B^\gamma_B\bar{C}_{ABCD}\bar{F}_{bC}B^\delta_D
\end{aligned}
\end{equation}\]
and for the second term, we have
\[\begin{equation}
\begin{aligned}
\bar{\BS}:[\nabla(\Be_aN^\gamma)\tra \nabla(\Be_bN^\delta)]
&= \bar{\BS}:[(\Be_a\dyd \BB^\gamma)\tra (\Be_b \dyd \BB^\delta)] \\
&= \bar{\BS}:[\BB^\gamma\dyd\BB^\delta] g_{ab} \\
&= B^\gamma_A \bar{S}_{AB}B^\delta_B g_{ab}
\end{aligned}
\end{equation}\]
where $g_{ab}$ are the components of the Eulerian metric tensor.
For the first term in $b$, we have
\[\begin{equation}
\bar{\BS} : \sym(\bar{\BF}\tra\nabla(\Be_aN^\gamma))
= \bar{\BS} : (\bar{\BF}\tra(\Be_a \dyd \BB^\gamma))
= \bar{S}_{AB} \bar{F}_{aA} B^\gamma_B
\end{equation}\]
Remaining terms can be calculated in a straightforward manner. We then have for
$\BA$ and $\Bb$:
\[\begin{equation}
\boxed{
\begin{aligned}
A_{ab}^{\gamma\delta}
&= \int_\CB
\bar{F}_{aA}B^\gamma_B\bar{C}_{ABCD}\bar{F}_{bC}B^\delta_D
+ B^\gamma_A \bar{S}_{AB}B^\delta_B g_{ab} \dV
\\
b_a^\gamma
&= -\int_\CB
\bar{S}_{AB} \bar{F}_{aA} B^\gamma_B \dV
+ \int_\CB\rho_0\bar{\Gamma}_aN^\gamma \dV
+ \int_{\del\CB_t}\bar{T}_aN^\gamma \dA
\end{aligned}
}
\end{equation}\]
The lowercase indices in $\bar{\Bgamma}$ and $\bar{\BT}$ might be confusing, but
in fact
\[\begin{equation}
\begin{aligned}
\Gamma_a(\BX,t) &= \gamma_a(\Bx, t) \circ \Bvarphi(\BX,t) \\
T_a(\BX,t) &= t_a(\Bx, t) \circ \Bvarphi(\BX,t) \\
\end{aligned}
\end{equation}\]
The system is solved for $\Vvphi$ at each Newton iteration with the following
update equation:
\[\begin{equation}
\Bvarphi \leftarrow \bar{\Bvarphi} + \Vvphi
\label{eq:lagrangianupdate1}
\end{equation}\]
Discretization of the Eulerian formulation parallels that of Lagrangian.
\[\begin{gather}
\boxed{
\begin{aligned}
A_{ab}^{\gamma\delta}
&= \int_\CB
\bar{B}^\gamma_c \bar{c}_{acbd}\bar{B}^\delta_d
+ \bar{B}^\gamma_e \bar{\tau}_{ef}\bar{B}^\delta_f g_{ab} \dV
\\
b_a^\gamma
&= -\int_\CB
\bar{\tau}_{ab} \bar{B}^\gamma_b \dV
+ \int_\CB\rho_0 \bar{\Gamma}_aN^\gamma \dV
+ \int_{\del\CB_t}\bar{T}_aN^\gamma \dA
\end{aligned}
}
\\
\text{or} \nonumber
\\
\boxed{
\begin{aligned}
A_{ab}^{\gamma\delta}
&= \int_{\bar{\CS}}
\bar{B}^\gamma_c \bar{c}^\sigma_{acbd}\bar{B}^\delta_d
+ \bar{B}^\gamma_e \bar{\sigma}_{ef}\bar{B}^\delta_f g_{ab} \,\dif\bar{v}
\\
b_a^\gamma
&= -\int_{\bar{\CS}}
\bar{\sigma}_{ab} \bar{B}^\gamma_b \,\dif\bar{v}
+ \int_{\bar{\CS}}\rho \bar{\gamma}_aN^\gamma \,\dif\bar{v}
+ \int_{\del\bar{\CS}_t}\bar{t}_aN^\gamma \,\dif\bar{a}
\end{aligned}
}
\end{gather}\]
Here, $\bar{\BB}^\gamma = \nabla_{\bar{x}} N^\gamma$ denote the spatial
gradients of the shape functions. One way of calculating is
$\bar{\BB}^\gamma = \bar{\BF}\invtra\BB^\gamma$, similar to
\eqref{eq:defgradidentity1}.
The update equation
\eqref{eq:lagrangianupdate1} holds for the Eulerian version.
Conclusion
The equations above in boxes contain all the information needed to implement the
nonlinear solution scheme of hyperelasticity.