\documentclass[12pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TCIDATA{OutputFilter=LATEX.DLL}
%TCIDATA{Version=4.10.0.2363}
%TCIDATA{Created=Sunday, September 19, 2004 15:07:04}
%TCIDATA{LastRevised=Wednesday, November 10, 2004 15:38:12}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Standard LaTeX\Blank - Standard LaTeX Article">}
%TCIDATA{CSTFile=40 LaTeX article.cst}

\newtheorem{theorem}{Theorem}
\newtheorem{acknowledgement}[theorem]{Acknowledgement}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{axiom}[theorem]{Axiom}
\newtheorem{case}[theorem]{Case}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conclusion}[theorem]{Conclusion}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{criterion}[theorem]{Criterion}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{solution}[theorem]{Solution}
\newtheorem{summary}[theorem]{Summary}
\newenvironment{proof}[1][Proof]{\noindent\textbf{#1.} }{\ \rule{0.5em}{0.5em}}
\input{tcilatex}

\begin{document}


\section{Integrating the stiffness matrix}

Index after comma mean differentiation. Repeated index implies summation.

\subsection{Laplace equation\label{sec:laplace}}

The basis functions $f$ are defined on a reference element $R$ with
coordinates $\xi _{j}$. The functions $x=x(\xi )$, $x=\left[ x_{i}\right] $
define the element $E$ as the image $x(R)$. Denote 
\[
x_{i,j}=\frac{dx_{i}}{d\xi _{j}} 
\]%
Basis functions $g$ on $E$ are defined via this mapping: 
\[
g_{a}(x)=f_{a}(\xi ) 
\]%
where $a$ is the index of the basis function. We wish to compute the entries
of the local stiffness matrix of element $E$, given by 
\begin{equation}
K_{ab}=\int_{E}cg_{a,j}g_{b,j}dx
\end{equation}%
where $a$ and $b$ are the numbers of the shape functions (or nodes), $j$ is
the index of the partial derivatives, and $\dot{c}>0$ is constant
coefficient.

Using the substitution $x=x(\xi )$ and the chain rule $g_{a,j}(x)=f_{a,p}(%
\xi )\xi _{p,j}(x)$ we get from (\ref{eq:stiff})

\[
K_{ab}=\int_{R}v_{ab}(\xi )d\xi 
\]%
where

\[
v_{ab}=cf_{a,p}\xi _{p,j}f_{b,q}\xi _{j}\det \left[ x_{m,n}\right] 
\]%
We do not know the partial derivatives 
\[
\xi _{p,j}=\frac{\partial \xi _{p}}{\partial x_{j}} 
\]%
but we are given the mapping $x=x(\xi )$ so $\xi =\xi (x)$ is the inverse
mapping and its partial derivatives are the entries of the inverse matrix to
the Jacobian matrix of $x$, 
\[
\left[ \xi _{p,j}\right] =\left[ x_{j,p}\right] ^{-1} 
\]%
To compute the stiffness matrix numerically we need the quadrature nodes $%
z_{n}$ and weights $w_{n}$ on the reference element $R,$%
\[
\int_{R}v(\xi )d\xi \approx w_{n}v(z_{n}) 
\]%
(again, with summation over the repeated index $n$). In conclusion, we get
the formula for evaluating the stiffness matrix%
\[
K_{ab}=w_{n}\left[ cf_{a,p}\xi _{p,j}f_{b,q}\xi _{q,j}\det \left[ x_{m,n}%
\right] \right] _{\xi =z_{n}}\text{ where }\left[ \xi _{p,j}\right] =\left[
x_{j,p}\right] ^{-1} 
\]%
So, to evaluate the stiffness matrix we need: quadrature nodes $z_{n}$ and
weights $w_{n}$ on the reference element $R,$partial derivatives $f_{a,p}$of
the reference basis functions at the quadrature nodes, partial derivatives $%
x_{m,n}$ of the mapping that defines the element $E=x(R)$ and the
coefficient $c$.

For quadrature of the right-hand side the situation is simpler:\ then we get%
\begin{eqnarray*}
F_{a} &=&\int_{E}F(x)g_{a}(x)dx \\
&=&\int_{R}\underbrace{f(x(\xi ))f_{a}(\xi )\det \left[ x_{m,n}\right] }%
_{u_{a}(\xi )}d\xi \\
&=&w_{n}u_{a}(z_{n})
\end{eqnarray*}

\subsection{General formula}

The basis functions $f$ are defined on a reference element $R$ with
coordinates $\xi _{j}$. The functions $x=x(\xi )$, $x=\left[ x_{i}\right] $
define the element $E$ as the image $x(R)$. Note that 
\[
x_{i,j}=\frac{dx_{i}}{d\xi _{j}}
\]%
Basis functions $g$ on $E$ are defined via this mapping: 
\[
g_{a}(x)=f_{a}(\xi )
\]%
where $a$ is the index of the basis function. We wish to compute the entries
of the local stiffness matrix of element $E$, given by 
\begin{equation}
K_{iakb}=\int_{E}c_{ijkl}g_{a,j}g_{b,l}dx  \label{eq:stiff}
\end{equation}%
where $a$ and $b$ are the numbers of the shape functions (or nodes), $j$ and 
$l$ are the indices of the partial derivatives, and $c_{ijkl}^{\prime }$ are
constant coefficients.

Note that the formula (\ref{eq:stiff}) is quite general: it works for
elasticity where $i$ and $k$ are the directions of the displacement ($%
c_{ijkl}$ are the same as the coefficients from the Hooke's law, see Sec. %
\ref{sec:elasticity}), and in applications where the range of the indices $i$
and $k$ is just $1$ unrelated to the dimension of the space (Sec. \ref%
{sec:laplace}). Using the substitution $x=x(\xi )$ and the chain rule $%
g_{a,j}(x)=f_{a,p}(\xi )\xi _{p,j}(x)$ we get from (\ref{eq:stiff})

\[
K_{iakb}=\int_{R}v_{iakb}(\xi )d\xi 
\]%
where

\[
v_{iakb}=c_{ijkl}^{\prime }f_{a,p}\xi _{p,j}f_{b,q}\xi _{q,l}\det \left[
x_{m,n}\right] 
\]%
We do not know the partial derivatives 
\[
\xi _{p,j}=\frac{\partial \xi _{p}}{\partial x_{j}} 
\]%
but we but we are given the mapping $x=x(\xi )$ so $\xi =\xi (x)$ is the
inverse mapping so the partials are the entries of the inverse matrix 
\[
\left[ \xi _{p,j}\right] =\left[ x_{m,n}\right] ^{-1} 
\]%
To compute the stiffness matrix numerically we need the quadrature nodes $%
z_{n}$ and weights $w_{n}$ on the reference element $R,$%
\[
\int_{R}v(\xi )d\xi \approx w_{n}v(z_{n}) 
\]%
(again, with summation over the repeated index $n$). In conclusion, we get
the formula for evaluating the stiffness matrix%
\begin{eqnarray*}
K_{iakb} &=&w_{n}v_{iakb}(z_{n}) \\
v_{iakb} &=&c_{ijkl}^{\prime }f_{a,p}\xi _{p,j}f_{b,q}\xi _{q,l}\det \left[
x_{m,n}\right] \\
\left[ \xi _{p,j}\right] &=&\left[ x_{m,n}\right] ^{-1}
\end{eqnarray*}%
So, to evaluate the stiffness matrix we need: quadrature nodes $z_{n}$ and
weights $w_{n}$ on the reference element $R,$partial derivatives $f_{a,p}$of
the reference basis functions at the quadrature nodes, partial derivatives $%
x_{m,n}$ of the mapping that defines the element $E=x(R)$ and the
coefficients $c_{ijkl}^{\prime }$.

Note that an efficient and stable way of evaluating determinant is as well
as the inverse matrix are from a decomposition such as the QR decomposition:%
\[
A=QR 
\]%
where $Q$ is an orthogonal matrix and $R$ is an upper triangular matrix. Then%
\begin{eqnarray*}
A^{-1} &=&R^{-1}Q^{-1}=R^{-1}Q^{T} \\
\det A &=&\det Q\det R=\pm r_{11}\cdots r_{nn}
\end{eqnarray*}%
See Numerical Analysis for details.

\subsection{Elasticity}

The bilinear form for elasticity is%
\[
a(u,v)=\int_{\Omega }u_{\left( i,j\right) }c_{ijkl}v_{\left( k,l\right) }dx
\]%
where the subscript in parenthesis means symmetric part:

\[
u_{\left( i,j\right) }=\frac{1}{2}\left( u_{i,j}+u_{j,i}\right) 
\]%
Because of the symmetry $c_{ijkl}=c_{jikl}$ we have 
\[
u_{\left( i,j\right) }c_{ijkl}=u_{i,j}c_{ijkl}
\]

and the stiffness matrix on element $E$ is obtained by integrating over just
one element $E$,%
\[
a_{E}(u,v)=\int_{\Omega }u_{i,j}c_{ijkl}v_{k,l}dx
\]%
For elasticity, each shape function is identified with a node and with one
of the displacement fields, so the basis functions are%
\[
g_{ai}=e_{i}g_{a}(x)
\]%
where $e_{i}$ is the unit vector of displacement in coordinate direction $i$%
. Then (\ref{eq:stiff})\ becomes%
\[
K_{iakb}=\int_{E}c_{ijkl}g_{a,j}g_{b,l}dx
\]

\end{document}