Entries for April 2018
Lumped L2 Projection
When utilizing Galerkin-type solutions for IBVPs, we often have to compute integrals using numerical methods such as Gauss quadrature. In such a solution, we solve for the values of a function at mesh nodes, whereas the integration takes place at the quadrature points. Depending on the case, we may need to compute the values of a function at mesh nodes, given their values at quadrature points, e.g. stress recovery for mechanical problems.
There are many ways of achieving this, such as superconvergent patch recovery. In this post, I wanted to document a widely used solution which is easy to implement, and which is used in research oriented codebases such as FEAP.
L2 Projection
Given a function $u \in L^2(\Omega)$, its projection into a finite element space $V_h\subset L^2(\Omega)$ is defined through the following optimization problem:
Find $u_h\in V_h$ such that
\[\begin{equation} \Pi(u_h) := \frac{1}{2}\lVert u_h-u \rVert^2_{L^2(\Omega)} \quad\rightarrow\quad \text{min} \end{equation}\]There is a unique solution to the problem since $\Pi(\cdot)$ is convex. Taking its variation, we have \(\begin{equation} D \Pi(u_h) \cdot v_h = \langle u_h-u, v_h \rangle = 0 \end{equation}\)
for all $v_h\in V_h$. Thus we have the following variational formulation
Find $u_h\in V_h$ such that
\[\begin{equation} \langle u_h,v_h\rangle = \langle u, v_h\rangle \end{equation}\]for all $v_h\in V_h$.
Here,
\[\begin{equation} \begin{alignedat}{3} m(u_h,v_h) &= \langle u_h,v_h\rangle && = \int_\Omega u_hv_h \,dx \quad\text{and} \\ b(v_h) &= \langle u, v_h\rangle && = \int_\Omega u v_h \,dx \end{alignedat} \end{equation}\]are our bilinear and linear forms respectively. Substituting FE discretizations $u_h = \sum_{J=1}^{\nnode} u^JN^J$ and $v_h = \sum_{I=1}^{\nnode} v^IN^I$, we have
\[\begin{equation} \suml{J=1}{\nnode} M^{I\!J} u^J = b^I \label{eq:projectionsystem1} \end{equation}\]for $I=1,\dots,\nnode$, where the FE matrix and vector are defined as
\[\begin{equation} \begin{alignedat}{3} M^{I\!J} &= m(N^J,N^I) &&= \int_\Omega N^JN^I \,dx \quad\text{and} \\ b^{I} &= b(N^I) &&= \int_\Omega u N^I \,dx \end{alignedat} \end{equation}\]Thus L2 projection requires the solution of a linear system
\[\boldsymbol{M}\boldsymbol{u}=\boldsymbol{b}\]which depending on the algorithm used, can have a complexity of at least $O(n^2)$ and at most $O(n^3)$.
Lumped L2 Projection
The L2 projection requires the solution of a system which can be computationally expensive. It is possible to convert the matrix—called the mass matrix in literature—to a diagonal form through a procedure called lumping.
The operator for row summation is defined as
\[\begin{equation} \rowsum{(\cdot)}_i := \suml{j=1}{\nnode} (\cdot)_{ij} \end{equation}\]For the mass matrix, we have
\[\begin{equation} \rowsum M^{I} = \suml{J=1}{\nnode} \int_\Omega N^JN^I \,dx = \int_\Omega N^I \,dx =: m^I \end{equation}\]since $\sum_{J=1}^{\nnode} N^J = 1$. Substituting the lumped mass matrix allows us to decouple the linear system of equations in \eqref{eq:projectionsystem1} and instead write
\[\begin{equation} m^I u^I = b^I \end{equation}\]for $I=1,\dots,\nnode$. The lumped L2 projection is then as simple as
\[\begin{equation} u^I = \frac{b^I}{m^I} = \frac{\displaystyle\int_\Omega u N^I\,dx}{\displaystyle\int_\Omega N^I \,dx} \end{equation}\]This results in a very efficient algorithm with $O(n)$ complexity.
Conclusion
Lumped L2 projection is a faster working approximation to L2 projection that is easy to implement for quick results. You can use it when developing a solution for an IBVP, and don’t want to wait too long when debugging, while not forgetting that it introduces some error.
-
Disadvantages of Engineering Notation in Finite Elements
Suppose we have the following stiffness matrix of linear elasticity:
\[\begin{equation} A^{I\!J}_{ij} = \int_\Omega B^I_k \, C_{ikjl} \, B^J_l \,dv \label{eq:engnot1} \end{equation}\]where $\boldsymbol{B}^I = \nabla N^I$ are the gradients of the shape functions $N^I$ and $\mathbb{C}$ is the linear elasticity tensor (you see the contraction of their components in the equation).
Despite being of the most explicit form, these types of indicial expressions are avoided in most texts on finite elements. There are two reasons for this:
- Engineers are not taught the Einstein summation convention.
- The presence of indices result in a seemingly cluttered expression.
They avoid the indicial expression by reshaping it into matrix multiplications. In engineering notation, the left- and right-hand sides are reshaped as
\[\begin{equation} A_{\alpha\beta} = \int_\Omega B_{\gamma\alpha}C_{\gamma\delta}B_{\delta\beta} \,dv \label{eq:engnot2} \end{equation}\]which allows us to write
\[\begin{equation} \boldsymbol{A} = \int_\Omega \tilde{\boldsymbol{B}}^T\tilde{\boldsymbol{C}}\tilde{\boldsymbol{B}} \,dv \label{eq:engnot3} \end{equation}\]The matrices $\tilde{\boldsymbol{B}}$ and $\tilde{\boldsymbol{C}}$ are set on with tildes in order to differentiate them from the boldface symbols used in the previous sections. Here,
- $\tilde{\boldsymbol{C}}$ is a matrix containing the unique components of the elasticity tensor $\mathbb{C}$, according to the Voigt notation. In this reshaping, only the minor symmetries are taken into account. If the dimension of the vectorial problem is $d$, then $\tilde{\boldsymbol{C}}$ is of the size $d(d+1)/2 \times d(d+1)/2$. For example, if the problem is 3 dimensional, $\tilde{\boldsymbol{C}}$ is of the size $6\times 6$:
- $\tilde{\boldsymbol{B}}$ is a $nd\times d(d+1)/2$ matrix whose components are adjusted so that \eqref{eq:engnot2} is equivalent to \eqref{eq:engnot1}. It has the components of $\boldsymbol{B}^I$ for $I=1,\dots,n$ where $n$ is the number of basis functions. Since $\tilde{\boldsymbol{B}}$ is adjusted to account for the reshaping of $\mathbb{C}$, it has many zero components. A 3d example:
Although \eqref{eq:engnot3} looks nice on paper, it is much less optimal for implementation. Implementing it requires the implementation of \eqref{eq:engnotB}, which adds another layer of complexity to the algorithm. The same cannot be said for \eqref{eq:engnotC}, because using Voigt notation might be more efficient in inelastic problems. In the most complex problems, the most efficient method is to implement \eqref{eq:engnot1} in conjunction with Voigt notation.
To prove the inefficiency of \eqref{eq:engnot3} we can readily compare it with \eqref{eq:engnot1} in terms of required number of iterations. Indices in \eqref{eq:engnot1} have the following ranges:
\[\begin{equation} I,J = 1,\dots,n \quad\text{and}\quad i,j,k,l = 1,\dots,d \end{equation}\]so $n^2d^4$ iterations are required. Indices in \eqref{eq:engnot2} have the following ranges:
\[\begin{equation} \alpha,\beta=1,\dots,nd \quad\text{and}\quad \gamma,\delta=1,\dots,d(d+1)/2 \end{equation}\]so
\[\begin{equation} (nd)^2\left(\frac{d(d+1)}{2}\right)^2 = n^2d^4\frac{(d+1)^2}{4} \end{equation}\]iterations are required. So engineering notation requires $(d+1)^2/4$ times more equations than index notation. For $d=2$, engineering notation is $2.25$ times slower and for $d=3$ it is $4$ times slower. For example, calculation of a stiffness matrix for $n=8$ and $d=3$ requires $20736$ iterations for engineering notation, whereas it only requires $5184$ iterations for index notation.
Although \eqref{eq:engnot3} seems less cluttered, what actually happens is that one trades off complexity in one expression for a much increased complexity in another one, in this case \eqref{eq:engnotB}. And to make it worse, it results in a slower algorithm.
The only obstacle to the widespread adoption of index notation seems to be its lack in undergraduate engineering curricula. If engineers were taught the index notation and summation convention as well as the formal notation, such expressions would not be as confusing at first sight. A good place would be in elementary calculus and physics courses, where one heavily uses vector calculus.
-
Variational Formulation of Elasticity
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.
Linearization of the Variational Formulation
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}$.
Discretization of the Lagrangian Form
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 Form
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.