Entries for November 2017
Vectorial Finite Elements
Many initial boundary value problems require solving for unknown vector fields, such as solving for displacements in a mechanical problem. Discretization of weak forms of such problems leads to higher-order linear systems which need to be reshaped to be solved by regular linear solvers. There are also more indices involved than a scalar problem, which can be confusing. In this post, I’ll try to elucidate the procedure by deriving for a basic higher-order system and giving an example.
The weak formulation of a linear vectorial problem reads
Find $\Bu\in V$ such that
\[\begin{equation} a(\Bu, \Bv) = b(\Bv) \end{equation}\]for all $\Bv\in V$.
Discretizations of vectorial problems requires the expansion of vectorial quantities as linear combinations of the basis vectors $\Be_i$:
\[\begin{equation} \Bu = \suml{i=1}{\ndim} u_i \,\Be_i \label{eq:discrete6} \end{equation}\]where $\cbr{u_i}_{i=1}^{\ndim}$ are the components corresponding to the basis vectors and $\ndim=\dim V$. Here, we chose Cartesian basis vectors for simplicity.
We can therefore express its discretization as
\[\begin{equation} \Bu_h = \suml{I=1}{\nnode} \Bu^I N^I = \suml{I=1}{\nnode}\suml{i=1}{\ndim} u^I_i \Be_i N^I. \label{eq:discrete7} \end{equation}\]Substituting discretized functions in the weak formulation, we obtain
\[\begin{equation} \begin{aligned} a(\Bu_h, \Bv_h) &= \suml{i=1}{\ndim}\suml{j=1}{\ndim} \, a(u_{h,j}\,\Be_j, v_{h,i}\,\Be_i)\\ &= \suml{I=1}{\nnode}\suml{J=1}{\nnode} \suml{i=1}{\ndim}\suml{j=1}{\ndim} u^J_j v^I_i \,a(\Be_j N^J, \Be_i N^I) \end{aligned} \end{equation}\]and
\[\begin{equation} b(\Bv_h) = \suml{i=1}{\ndim} b(v_{h,i}\,\Be_i) = \suml{I=1}{\nnode} \suml{i=1}{\ndim} v^I_i b(\Be_i N^I). \end{equation}\]We define the following arrays
\[\begin{equation} \boxed{ \begin{aligned} A^{I\!J}_{ij} &= a(\Be_j N^J, \Be_i N^I) \\ b^{I}_{i} &= b(\Be_i N^I). \end{aligned} } \end{equation}\]Hence we can express the linear system
\[\begin{equation} a(\Bu_h,\Bv_h) = b(\Bv_h) \end{equation}\]as
\[\begin{equation} \suml{I=1}{\nnode}\suml{J=1}{\nnode} \suml{i=1}{\ndim}\suml{j=1}{\ndim} u^J_j v^I_i \,A^{I\!J}_{ij} = \suml{I=1}{\nnode} \suml{i=1}{\ndim} v^I_i b^{I}_{i}. \end{equation}\]For arbitrary $\Bv_h$, this yields the following system of equations
\[\begin{equation} \boxed{ \suml{J=1}{\nnode} \suml{j=1}{\ndim} A^{I\!J}_{ij} \,u^J_j = b^{I}_{i} } \label{eq:discrete8} \end{equation}\]for $i=1,\dots,\ndim$ and $I=1,\dots,\nnode$.
We reshape this higher-order system as shown in the previous post Reshaping Higher Order Linear Systems:
\[\begin{equation} \BAhat \Buhat = \Bbhat \end{equation}\]by defining a map $i_d$ that maps original indices to the reshaped indices
\[\begin{equation} i_d := \left\{ \begin{array}{rl} [1,\nnode]\times[1,\ndim] & \to [1,\nnode\ndim]\\[1ex] (I,i) & \mapsto \ndim(I-1)+i\\ \end{array} \right. \end{equation}\]where we used 1-based indexing of the arrays. We set
\[\begin{equation} \boxed{ \begin{alignedat}{3} \alpha &:= i_d(I,i) &&= \ndim(I-1) + i \\ \beta &:= i_d(J,j) &&= \ndim(J-1) + j \\ \end{alignedat} } \end{equation}\]and write
\[\begin{equation} \hat{A}_{\alpha\beta} = A^{I\!J}_{ij} \quad,\quad \hat{u}_{\beta} = u^{J}_{j} \eqand \hat{b}_{\alpha} = b^{I}_{i} \end{equation}\]The inverse index mapping can be obtained as shown in the previous post.
Example: Linear Elasticity
Our initial-boundary value problem is
\[\begin{equation} \begin{alignedat}{4} -\div\Bsigma &= \rho\Bgamma \qquad&& \text{in} \qquad&& \Omega\\ \Bu &= \bar{\Bu} && \text{on} && \del\Omega_u \\ \Bt &= \bar{\Bt} && \text{on} && \del\Omega_t \\ \end{alignedat} \end{equation}\]The weak formulation reads
Find $\Bu\in V$ such that
\[\begin{equation} -\int_\Omega \div\Bsigma \dtp \Bv \dv = \int_\Omega \rho \Bgamma\dtp\Bv \dv \end{equation}\]for all $\Bv\in V$ where $V=H^1(\Omega)$.
We apply integration by parts on the left-hand side
\[\begin{equation} \int_\Omega \div\Bsigma\dtp\Bv \dv = \int_\Omega \div(\Bsigma\Bv) \dv - \int_\Omega \Bsigma : \nabla\Bv \dv \end{equation}\]and apply the divergence theorem to the first resulting term:
\[\begin{equation} \int_\Omega \div(\Bsigma\Bv) \dv = \int_{\del\Omega_t} \bar{\Bt}\dtp\Bv \da \end{equation}\]Substituting the linear stress $\Bsigma=\IC:\Bvareps=\IC:\nabla\Bu$, we obtain the following variational forms:
\[\begin{align} \label{eq:linelastdiscretebilinear} a(\Bu,\Bv) &= \int_\Omega \nabla\Bv:\IC:\nabla\Bu \dv \\ \label{eq:linelastdiscretelinear} b(\Bv) &= \int_\Omega \rho \Bgamma\dtp\Bv \dv + \int_{\del\Omega_t} \bar{\Bt}\dtp\Bv \da \end{align}\]We have the following discretizations of the unknown function and test function
\[\begin{equation} \Bu_h = \suml{J=1}{\nnode} \Bu^J N^J \eqand \Bv_h = \suml{I=1}{\nnode} \Bv^I N^I. \end{equation}\]With the given discretizations, the matrix corresponding to \eqref{eq:linelastdiscretebilinear} can be calculated as
\[\begin{equation} \begin{aligned} A^{I\!J}_{ij} = a(\Be_j N^J, \Be_i N^I) &= \int_\Omega \nabla(\Be_iN^I) : \IC : \nabla(\Be_jN^J) \dv \\ &= \int_\Omega (\Be_i\dyd \nabla N^I) : \IC : (\Be_j \dyd \nabla N^J) \dv \\ &= \int_\Omega \partd{N^I}{x_k} \, C_{ikjl} \, \partd{N^J}{x_l} \dv, \end{aligned} \end{equation}\]and finally obtain
\[\begin{equation} \boxed{ A^{I\!J}_{ij} = \int_\Omega B^I_k \, C_{ikjl} \, B^J_l \dv \,. } \end{equation}\]The vector corresponding to \eqref{eq:linelastdiscretelinear} is calculated as
\[\begin{equation} b^{I}_{i} = b(\Be_i N^I) = \int_{\del\Omega_t} \bar{\Bt}\dtp(\Be_iN^I) \da + \int_\Omega\rho\Bgamma\dtp(\Be_iN^I) \dv \end{equation}\]which yields
\[\begin{equation} \boxed{ b^{I}_{i} = \int_{\del\Omega_t} \bar{t}_i N^I \da + \int_\Omega \rho \gamma_i N^I \dv } \end{equation}\]-
Reshaping Higher Order Linear Systems
In vectorial problems, we end up with linear systems of higher order, such as
\[\begin{equation} \suml{k=1}{N} \suml{l=1}{M} A_{ijkl} \, u_{kl} = b_{ij} \label{eq:1} \end{equation}\]for $i=1,\dots,N$ and $j=1,\dots,M$.
These systems cannot be solved readily with existing software. In order to be able to solve them with existing software, we need to reshape them by defining a matrix of matrices $\BAhat$ and vector of vectors $\Buhat$ and $\Bbhat$:
\[\begin{equation} \underbrace{ \left[ \begin{array}{ccc|c|ccc} A_{1111} & \cdots & A_{111M} & & A_{11N1} & \cdots & A_{11NM} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ A_{1M11} & \cdots & A_{1M1M} & & A_{1MN1} & \cdots & A_{1MNM} \\[1ex] \hline & \vdots & &\ddots & & \vdots & \\ \hline &&&&&& \\[-1.5ex] A_{N111} & \cdots & A_{N11M} & & A_{N1N1} & \cdots & A_{N1NM} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ A_{NM11} & \cdots & A_{NM1M} & & A_{NMN1} & \cdots & A_{NMNM} \end{array} \right] }_{\BAhat} \underbrace{ \left[ \begin{array}{c} u_{11}\\ \vdots \\ u_{1M} \\[1ex] \hline \vdots\\ \hline \\[-1.5ex] u_{N1} \\ \vdots \\ u_{NM} \end{array} \right] }_{\Buhat} = \underbrace{ \left[ \begin{array}{c} b_{11}\\ \vdots \\ b_{1M} \\[1ex] \hline \vdots\\ \hline \\[-1.5ex] b_{N1} \\ \vdots \\ b_{NM} \end{array} \right]}_{\Bbhat} \end{equation}\]This allows us to express the linear system as
\[\begin{equation} \BAhat \Buhat = \Bbhat \label{eq:3} \end{equation}\]Here, we reshape the system by defining a map $i_d$ that maps original indices to the reshaped indices
\[\begin{equation} i_d := \left\{ \begin{array}{rl} [1,N]\times[1,M] & \to [1,NM]\\[1ex] (i,j) & \mapsto M(i-1)+j\\ \end{array} \right. \end{equation}\]where we used 1-based indexing of the arrays. We set
\[\begin{equation} \boxed{ \begin{alignedat}{3} \alpha &:= i_d(i,j) &&= M(i-1) + j \\ \beta &:= i_d(k,l) &&= M(k-1) + l \\ \end{alignedat} } \end{equation}\]and write
\[\begin{equation} \hat{A}_{\alpha\beta} = A_{ijkl} \quad,\quad \hat{u}_{\beta} = u_{kl} \eqand \hat{b}_{\alpha} = b_{ij} \end{equation}\]For reference, the inverse of the index mapping reads
\[\begin{equation} i_d^{-1} := \left\{ \begin{array}{rl} [1,NM] & \to [1,N]\times[1,M] \\[1ex] \alpha & \mapsto (1+(\alpha-\modop(\alpha,M))/M\;,\; \modop(\alpha,M)) \end{array} \right. \end{equation}\]Thus, we have for our reshaped indices,
\[\begin{equation} \begin{aligned} j &= \modop(\alpha,M) \\ i &= 1+(\alpha-j)/M \end{aligned} \quad\eqand\quad \begin{aligned} l &= \modop(\beta,M) \\ k &= 1+(\beta-l)/M \end{aligned} \end{equation}\]Expressed as a regular linear system \eqref{eq:3}, the higher-order system \eqref{eq:1} can be solved with a linear solver such as LAPACK.
-
Linear Finite Elements
Beginning with this post, I’ll be publishing about the basics of finite element formulations, from personal notes that accumulated over the years. This one is about linear and scalar problems which came to be the “Hello World” for FE. Details regarding spaces and discretization are omitted for the sake of brevity. For those who want to delve into theory, I recommend “The Finite Element Method: Theory, Implementation, and Applications” by Larson and Bengzon.
The weak formulation of a canonical linear problem reads
Find $u\in V$ such that
\[\begin{equation} a(u, v) = b(v) \label{eq:femlinear1} \end{equation}\]for all $v \in V$ where $a(\cdot, \cdot)$ is a bilinear form and $b(\cdot)$ is a linear form.
We define the discretization of $u$ as
\[\begin{equation} u_h := \suml{J=1}{\nnode} u^J N^J ,\quad u_h \in V_h \quad\text{where}\quad V_h\subset V \end{equation}\]The discretization $u_h$ is a linear combination of basis functions $N^J$ and corresponding scalars $u^J$, $J=1,\dots,\nnode$ so that $V_h$ is a subset of $V$. The discretization of \eqref{eq:femlinear1} then reads
\[\begin{equation} a(u_h, v_h) = b(v_h) \quad \forall v_h \in V_h . \end{equation}\]We then have
\[\begin{equation} a\rbr{\suml{J=1}{\nnode} u^J N^J, \suml{I=1}{\nnode} v^I N^I} = b\rbr{\suml{I=1}{\nnode} v^I N^I} \end{equation}\]Using the linearity properties,
\[\begin{equation} a(\alpha u, \beta v) = \alpha\beta\, a(u,v) \eqand b(\alpha v) = \alpha b(v) \end{equation}\]we obtain
\[\begin{equation} \suml{I=1}{\nnode} \suml{J=1}{\nnode} u^J v^I a(N^J, N^I) = \suml{I=1}{\nnode} v^I b(N^I) . \label{eq:femlinear2} \end{equation}\]For arbitrary test function values $v^I$, we can express \eqref{eq:femlinear2} as a system of $\nnode$ equations
\[\begin{equation} \suml{J=1}{\nnode} u^J a(N^J, N^I) = b(N^I) \label{eq:femlinear3} \end{equation}\]for $I = 1,2,\dots,\nnode$. If we expand the summations as
\[\begin{alignat*}{6} & a(N^1, N^1) u^1 &&+ a(N^2, N^1) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^1) u^{\nnode} &&\quad=\quad b(N^1) \\ & a(N^1, N^2) u^1 &&+ a(N^2, N^2) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^2) u^{\nnode} &&\quad=\quad b(N^2) \\ & \qquad\vdots && \qquad\quad\;\vdots && \quad\;\;\vdots && \qquad\qquad\vdots && \qquad\qquad\vdots \\ & a(N^1, N^{\nnode}) u^1 &&+ a(N^2, N^{\nnode}) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^{\nnode}) u^{\nnode} &&\quad=\quad b(N^{\nnode}) \end{alignat*}\]we can see that the terms with $a$ constitute a matrix $\BA$ and the terms with $b$ constitute a vector $\Bb$, allowing us to write
\[\begin{equation} \BA\Bu = \Bb \label{eq:discrete9} \end{equation}\]where we chose to express the unknown coefficients $u^I$ as a vector $\Bu = [u^1,u^2,\dots,u^{\nnode}]\tra$.
\It can be seen that the components of the $\BA$ and $\Bb$ are defined as
\[\begin{equation} \boxed{ \Aelid{I\!J}{} = a(N^J,N^I) \eqand b^I = b(N^I), } \end{equation}\]we can express the linear system as
\[\begin{alignat*}{6} & \Aelid{11}{} u^1 &&+ \Aelid{12}{} u^2 &&+ \cdots &&+ \Aelid{1\nnode}{} u^{\nnode} &&\quad=\quad b^1 \\ & \Aelid{21}{} u^1 &&+ \Aelid{22}{} u^2 &&+ \cdots &&+ \Aelid{2\nnode}{} u^{\nnode} &&\quad=\quad b^2 \\ & \quad\vdots && \qquad\;\vdots && \quad\;\;\vdots && \qquad\;\vdots && \qquad\quad\;\;\vdots \\ & \Aelid{\nnode 1}{} u^1 &&+ \Aelid{\nnode 2}{} u^2 &&+ \cdots &&+ \Aelid{\nnode\nnode}{} u^{\nnode} &&\quad=\quad b^{\nnode} \end{alignat*}\]Note that with the given definitions, \eqref{eq:femlinear3} becomes
\[\begin{equation} \boxed{ \suml{J=1}{\nnode} \Aelid{I\!J}{} \,u^J = b^I \quad\text{for}\quad I=1,2,\dots\nnode. } \label{eq:discrete10} \end{equation}\]Example: Poisson’s Equation
In the weak form of Poisson’s equation
\[\begin{equation} \begin{alignedat}{4} - \Var u &= f \quad && \text{in} \quad && \Omega \\ u &= 0 \quad && \text{on} \quad && \del\Omega \end{alignedat} \end{equation}\]The weak formulation reads
Find $u\in V$ such that
\[\begin{equation} - \int_\Omega \Delta(u) v \dv= \int_\Omega f v \dv \end{equation}\]for all $v\in V$ where $V=H^1_0(\Omega)$.
Applying integration by parts and divergence theorem on the left-hand side
\[\begin{equation} \begin{aligned} \int_\Omega \Delta(u) v \dv &= \int_\Omega \nabla \dtp (\nabla (u) v) \dv - \int_\Omega \nabla u\dtp\nabla v \dv \\ &= \underbrace{\int_{\del\Omega} v (\nabla u\dtp\Bn) \da}_{v = 0 \text{ on } \del\Omega} - \int_\Omega \nabla u\dtp\nabla v \dv \\ \end{aligned} \end{equation}\]We have the following variational forms:
\[\begin{equation} \begin{aligned} a(u,v) &= \int_{\Omega} \nabla u \dtp \nabla v \dv\\ b(v) &= \int_{\Omega} f \, v \dv\\ \end{aligned} \end{equation}\]Following \eqref{eq:femlinear3}, we can calculate the stiffness matrix $\BA$ as
\[\begin{equation} \begin{aligned} \Aelid{I\!J}{} = a(N^J, N^I) &= \int_{\Omega} \nabla N^J \dtp \nabla N^I \dv \\ &= \int_{\Omega} \BB^J \dtp \BB^I \dv \end{aligned} \end{equation}\]where we have defined the gradient of the basis functions as
\[\begin{equation} \BB^I := \nabla N^I\,. \end{equation}\]Similarly, we integrate the force term into a vector $\Bb$ as
\[\begin{equation} \begin{aligned} b^I &= \int_{\Omega} f N^I \dv \end{aligned} \end{equation}\]