-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% intro.html was generated from this file
-%% with latex2html and some handwork
-%% (copying out the relevant parts from the
-%% generated html file, replacing IMG=
-%% by the proper path)
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\documentclass{article}
-\usepackage{amsmath}
-\usepackage{amsfonts}
+<a name="Intro"></a>
+<h1>Introduction</h1>
-\renewcommand{\vec}[1]{{\mathbf #1}}
-\renewcommand{\div}{\nabla \cdot}
-\begin{document}
In real life, most partial differential equations are really systems
of equations. Accordingly, the solutions are usually
which is subject to a force. Of course, the force is also
vector-valued, meaning that in each point it has a direction and an
absolute value. The elastic equations are the following:
-$$
+@f[
-
\partial_j (c_{ijkl} \partial_k u_l)
=
f_i,
\qquad
i=1\ldots d,
-$$
+@f]
where the values $c_{ijkl}$ are the stiffness coefficients and
will usually depend on the space coordinates. In
many cases, one knows that the material under consideration is
isotropic, in which case by introduction of the two coefficients
$\lambda$ and $\mu$ the coefficient tensor reduces to
-$$
+@f[
c_{ijkl}
=
\lambda \delta_{ij} \delta_{kl} +
\mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}).
-$$
+@f]
The elastic equations can then be rewritten in much simpler a form:
-$$
+@f[
-
- \nabla \lambda (\div \vec u)
+ \nabla \lambda (\nabla\cdot {\mathbf u})
-
- (\nabla \cdot \mu \nabla) \vec u
+ (\nabla \cdot \mu \nabla) {\mathbf u}
-
- \div \mu (\nabla \vec u)^T
+ \nabla\cdot \mu (\nabla {\mathbf u})^T
=
- \vec f,
-$$
+ {\mathbf f},
+@f]
and the respective bilinear form is then
-$$
- a(\vec u, \vec v) =
+@f[
+ a({\mathbf u}, {\mathbf v}) =
\left(
- \lambda \div \vec u, \div \vec v
+ \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v}
\right)_\Omega
+
\sum_{i,j}
\left(
\mu \partial_i u_j, \partial_j v_i
\right)_\Omega,
-$$
+@f]
or also writing the first term a sum over components:
-$$
- a(\vec u, \vec v) =
+@f[
+ a({\mathbf u}, {\mathbf v}) =
\sum_{i,j}
\left(
\lambda \partial_l u_l, \partial_k v_k
\left(
\mu \partial_i u_j, \partial_j v_i
\right)_\Omega.
-$$
+@f]
How do we now assemble the matrix for such an equation? The first thing we
let $n$ be the number of shape functions for the scalar finite element of
which we build the vector element (for example, we will use bilinear functions
for each component of the vector-valued finite element, so the scalar finite
-element is the \texttt{FEQ1} element which we have used in previous examples
+element is the <code>FE_Q(1)</code> element which we have used in previous examples
already, and $n=4$ in two space dimensions). Further, let $N$ be the number of
shape functions for the vector element; in two space dimensions, we need $n$
shape functions for each component of the vector, so $N=2n$. Then, the $i$th
shape function of the vector element has the form
-$$
- \Phi_i(\vec x) = \varphi_{base(i)}(\vec x)\ \vec e_{comp(i)},
-$$
+@f[
+ \Phi_i({\mathbf x}) = \varphi_{base(i)}({\mathbf x})\ {\mathbf e}_{comp(i)},
+@f]
where $e_l$ is the $l$th unit vector, $comp(i)$ is the function that tells
us which component of $\Phi_i$ is the one that is nonzero (for
each vector shape function, only one component is nonzero, and all others are
For example (though this sequence of shape functions is not
guaranteed, and you should not rely on it),
the following layout could be used by the library:
-\begin{center}
-\begin{multline*}
- \Phi_0(\vec x) =
- \begin{pmatrix}
- \varphi_0(\vec x) \\ 0
- \end{pmatrix},
- \qquad
- \Phi_1(\vec x) =
- \begin{pmatrix}
- 0 \\ \varphi_0(\vec x)
- \end{pmatrix},
+@f{eqnarray*}
+ \Phi_0({\mathbf x}) &=&
+ \left(\begin{array}{c}
+ \varphi_0({\mathbf x}) \\ 0
+ \end{array}\right),
\\
- \Phi_2(\vec x) =
- \begin{pmatrix}
- \varphi_1(\vec x) \\ 0
- \end{pmatrix},
- \qquad
- \Phi_3(\vec x) =
- \begin{pmatrix}
- 0 \\ \varphi_1(\vec x)
- \end{pmatrix},
+ \Phi_1({\mathbf x}) &=&
+ \left(\begin{array}{c}
+ 0 \\ \varphi_0({\mathbf x})
+ \end{array}\right),
+ \\
+ \Phi_2({\mathbf x}) &=&
+ \left(\begin{array}{c}
+ \varphi_1({\mathbf x}) \\ 0
+ \end{array}\right),
+ \\
+ \Phi_3({\mathbf x}) &=&
+ \left(\begin{array}{c}
+ 0 \\ \varphi_1({\mathbf x})
+ \end{array}\right),
\ldots
-\end{multline*}
-\end{center}
+@f}
where here
-$$
+@f[
comp(0)=0, \quad comp(1)=1, \quad comp(2)=0, \quad comp(3)=1, \quad \ldots
-$$
-$$
+@f]
+@f[
base(0)=0, \quad base(1)=0, \quad base(2)=1, \quad base(3)=1, \quad \ldots
-$$
+@f]
In all but very rare cases, you will not need to know which shape function
$\varphi_{base(i)}$ of the scalar element belongs to a shape function $\Phi_i$
of the vector element. Let us therefore define
-$$
+@f[
\phi_i = \varphi_{base(i)}
-$$
+@f]
by which we can write the vector shape function as
-$$
- \Phi_i(\vec x) = \phi_{i}(\vec x)\ \vec e_{comp(i)}.
-$$
+@f[
+ \Phi_i({\mathbf x}) = \phi_{i}({\mathbf x})\ {\mathbf e}_{comp(i)}.
+@f]
You can now safely forget about the function $base(i)$, at least for the rest
of this example program.
Now using this vector shape functions, we can write the discrete finite
element solution as
-$$
- \vec u_h(\vec x) =
- \sum_i \Phi_i(\vec x)\ u_i
-$$
-with scalar coefficients $u_i$. If we define an analog function $\vec v_h$ as
+@f[
+ {\mathbf u}_h({\mathbf x}) =
+ \sum_i \Phi_i({\mathbf x})\ u_i
+@f]
+with scalar coefficients $u_i$. If we define an analog function ${\mathbf v}_h$ as
test function, we can write the discrete problem as follows: Find coefficients
$u_i$ such that
-$$
- a(\vec u_h, \vec v_h) = (\vec f, \vec v_h)
+@f[
+ a({\mathbf u}_h, {\mathbf v}_h) = ({\mathbf f}, {\mathbf v}_h)
\qquad
- \forall \vec v_h.
-$$
+ \forall {\mathbf v}_h.
+@f]
If we insert the definition of the bilinear form and the representation of
-$\vec u_h$ and $\vec v_h$ into this formula:
-\begin{center}
-\begin{multline*}
+${\mathbf u}_h$ and ${\mathbf v}_h$ into this formula:
+@f{eqnarray*}
\sum_{i,j}
u_i v_j
\sum_{k,l}
f_l,
(\Phi_j)_l
\right)_\Omega.
-\end{multline*}
-\end{center}
+@f}
We note that here and in the following, the indices $k,l$ run over spatial
directions, i.e. $0\le k,l < d$, and that indices $i,j$ run over degrees
of freedoms.
The local stiffness matrix on cell $K$ therefore has the following entries:
-$$
+@f[
A^K_{ij}
=
\sum_{k,l}
\mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l
\right)_K
\right\},
-$$
+@f]
where $i,j$ now are local degrees of freedom and therefore $0\le i,j < N$.
In these formulas, we always take some component of the vector shape functions
$\Phi_i$, which are of course given as follows (see their definition):
-$$
+@f[
(\Phi_i)_l = \phi_i \delta_{l,comp(i)},
-$$
+@f]
with the Kronecker symbol $\delta_{nm}$. Due to this, we can delete some of
the sums over $k$ and $l$:
-\begin{center}
-\begin{align*}
+@f{eqnarray*}
A^K_{ij}
- &=
+ &=&
\sum_{k,l}
\Bigl\{
\left(
\partial_k \phi_j\ \delta_{k,comp(j)}
\right)_K
\\
- &\qquad\qquad +
+ &\qquad\qquad& +
\left(
\mu \partial_l \phi_i\ \delta_{k,comp(i)},
\partial_l \phi_j\ \delta_{k,comp(j)}
\right)_K
\Bigr\}
\\
- &=
+ &=&
\left(
\lambda \partial_{comp(i)} \phi_i,
\partial_{comp(j)} \phi_j
\partial_{comp(i)} \phi_j
\right)_K
\\
- &=
+ &=&
\left(
\lambda \partial_{comp(i)} \phi_i,
\partial_{comp(j)} \phi_j
\mu \partial_{comp(j)} \phi_i,
\partial_{comp(i)} \phi_j
\right)_K.
-\end{align*}
-\end{center}
+@f}
Likewise, the contribution of cell $K$ to the right hand side vector is
-\begin{center}
-\begin{align*}
+@f{eqnarray*}
f^K_j
- &=
+ &=&
\sum_l
\left(
f_l,
(\Phi_j)_l
\right)_K
\\
- &=
+ &=&
\sum_l
\left(
f_l,
\phi_j \delta_{l,comp(j)}
\right)_K
\\
- &=
+ &=&
\left(
f_{comp(j)},
\phi_j
\right)_K.
-\end{align*}
-\end{center}
+@f}
This is the form in which we will implement the local stiffness matrix and
right hand side vectors.
-As a final note: in the step-17 example program, we will revisit the elastic
+As a final note: in the @ref step_17 "step-17" example program, we will revisit the elastic
problem laid out here, and will show how to solve it in parallel on a cluster
of computers. The resulting program will thus be able to solve this problem to
-significantly higher accuracy, and more efficiently if this is required.
-
+significantly higher accuracy, and more efficiently if this is
+required. In addition, in @ref step_20 "step-20", we will revisit some
+vector-valued problems and show a few techniques that may make it
+simpler to actually go through all the stuff shown above, with
+<code>FiniteElement::system_to_component_index</code> etc.
-\end{document}
+++ /dev/null
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-<P>
-In real life, most partial differential equations are really systems
-of equations. Accordingly, the solutions are usually
-vector-valued. The deal.II library supports such problems, and we will show
-that that is mostly rather simple. The only more complicated problems
-are in assembling matrix and right hand side, but these are easily
-understood as well.
-
-<P>
-In the example, we will want to solve the elastic equations. They are
-an extension to Laplace's equation with a vector-valued solution that
-describes the displacement in each space direction of a rigid body
-which is subject to a force. Of course, the force is also
-vector-valued, meaning that in each point it has a direction and an
-absolute value. The elastic equations are the following:
-<!-- MATH
- \begin{displaymath}
--
- \partial_j (c_{ijkl} \partial_k u_l)
- =
- f_i,
- \qquad
- i=1\ldots d,
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="248" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img1.png"
- ALT="$\displaystyle -
-\partial_j (c_{ijkl} \partial_k u_l)
-=
-f_i,
-\qquad
-i=1\ldots d,
-$">
-</DIV><P></P>
-where the values <IMG
- WIDTH="33" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img2.png"
- ALT="$ c_{ijkl}$"> are the stiffness coefficients and
-will usually depend on the space coordinates. In
-many cases, one knows that the material under consideration is
-isotropic, in which case by introduction of the two coefficients
-<IMG
- WIDTH="13" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img3.png"
- ALT="$ \lambda$"> and <IMG
- WIDTH="14" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img4.png"
- ALT="$ \mu$"> the coefficient tensor reduces to
-<!-- MATH
- \begin{displaymath}
-c_{ijkl}
- =
- \lambda \delta_{ij} \delta_{kl} +
- \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}).
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="241" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img5.png"
- ALT="$\displaystyle c_{ijkl}
-=
-\lambda \delta_{ij} \delta_{kl} +
-\mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}).
-$">
-</DIV><P></P>
-
-<P>
-The elastic equations can then be rewritten in much simpler a form:
-<!-- MATH
- \begin{displaymath}
--
- \nabla \lambda (\div\vec u)
- -
- (\nabla \cdot \mu \nabla) \vec u
- -
- \div\mu (\nabla \vec u)^T
- =
- \vec f,
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="309" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img6.png"
- ALT="$\displaystyle -
-\nabla \lambda (\div\vec u)
--
-(\nabla \cdot \mu \nabla) \vec u
--
-\div\mu (\nabla \vec u)^T
-=
-\vec f,
-$">
-</DIV><P></P>
-and the respective bilinear form is then
-<!-- MATH
- \begin{displaymath}
-a(\vec u, \vec v) =
- \left(
- \lambda \div\vec u, \div\vec v
- \right)_\Omega
- +
- \sum_{i,j}
- \left(
- \mu \partial_i u_j, \partial_i v_j
- \right)_\Omega,
- +
- \sum_{i,j}
- \left(
- \mu \partial_i u_j, \partial_j v_i
- \right)_\Omega,
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="477" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img7.png"
- ALT="$\displaystyle a(\vec u, \vec v) =
-\left(
-\lambda \div\vec u, \div\vec v
-\right)...
-...\Omega,
-+
-\sum_{i,j}
-\left(
-\mu \partial_i u_j, \partial_j v_i
-\right)_\Omega,
-$">
-</DIV><P></P>
-or also writing the first term a sum over components:
-<!-- MATH
- \begin{displaymath}
-a(\vec u, \vec v) =
- \sum_{i,j}
- \left(
- \lambda \partial_l u_l, \partial_k v_k
- \right)_\Omega
- +
- \sum_{k,l}
- \left(
- \mu \partial_i u_j, \partial_i v_j
- \right)_\Omega,
- +
- \sum_{i,j}
- \left(
- \mu \partial_i u_j, \partial_j v_i
- \right)_\Omega.
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="492" HEIGHT="53" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img8.png"
- ALT="$\displaystyle a(\vec u, \vec v) =
-\sum_{i,j}
-\left(
-\lambda \partial_l u_l, \pa...
-...\Omega,
-+
-\sum_{i,j}
-\left(
-\mu \partial_i u_j, \partial_j v_i
-\right)_\Omega.
-$">
-</DIV><P></P>
-
-<P>
-How do we now assemble the matrix for such an equation? The first thing we
-need is some knowledge about how the shape functions work in the case of
-vector-valued finite elements. Basically, this comes down to the following:
-let <IMG
- WIDTH="14" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img9.png"
- ALT="$ n$"> be the number of shape functions for the scalar finite element of
-which we build the vector element (for example, we will use bilinear functions
-for each component of the vector-valued finite element, so the scalar finite
-element is the <TT>FEQ1</TT> element which we have used in previous examples
-already, and <IMG
- WIDTH="43" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img10.png"
- ALT="$ n=4$"> in two space dimensions). Further, let <IMG
- WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img11.png"
- ALT="$ N$"> be the number of
-shape functions for the vector element; in two space dimensions, we need <IMG
- WIDTH="14" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img9.png"
- ALT="$ n$">
-shape functions for each component of the vector, so <IMG
- WIDTH="57" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img12.png"
- ALT="$ N=2n$">. Then, the <IMG
- WIDTH="10" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img13.png"
- ALT="$ i$">th
-shape function of the vector element has the form
-<!-- MATH
- \begin{displaymath}
-\Phi_i(\vec x) = \varphi_{base(i)}(\vec x)\ \vec e_{comp(i)},
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="200" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img14.png"
- ALT="$\displaystyle \Phi_i(\vec x) = \varphi_{base(i)}(\vec x) \vec e_{comp(i)},
-$">
-</DIV><P></P>
-where <IMG
- WIDTH="16" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img15.png"
- ALT="$ e_l$"> is the <IMG
- WIDTH="9" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img16.png"
- ALT="$ l$">th unit vector, <IMG
- WIDTH="58" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img17.png"
- ALT="$ comp(i)$"> is the function that tells
-us which component of <IMG
- WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img18.png"
- ALT="$ \Phi_i$"> is the one that is nonzero (for
-each vector shape function, only one component is nonzero, and all others are
-zero). <!-- MATH
- $\varphi_{base(i)}(x)$
- -->
-<IMG
- WIDTH="76" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img19.png"
- ALT="$ \varphi_{base(i)}(x)$"> describes the space dependence of the shape
-function, which is taken to be the <IMG
- WIDTH="52" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img20.png"
- ALT="$ base(i)$">-th shape function of the scalar
-element. Of course, while <IMG
- WIDTH="10" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img13.png"
- ALT="$ i$"> is in the range <!-- MATH
- $0,\ldots,N-1$
- -->
-<IMG
- WIDTH="89" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img21.png"
- ALT="$ 0,\ldots,N-1$">, the functions
-<IMG
- WIDTH="58" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img17.png"
- ALT="$ comp(i)$"> and <IMG
- WIDTH="52" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img20.png"
- ALT="$ base(i)$"> have the ranges <IMG
- WIDTH="27" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img22.png"
- ALT="$ 0,1$"> (in 2D) and <!-- MATH
- $0,\ldots,n-1$
- -->
-<IMG
- WIDTH="84" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img23.png"
- ALT="$ 0,\ldots,n-1$">,
-respectively.
-
-<P>
-For example (though this sequence of shape functions is not
-guaranteed, and you should not rely on it),
-the following layout could be used by the library:
-<DIV ALIGN="CENTER">
-</DIV><P></P>
-<DIV ALIGN="CENTER"><!-- MATH
- \begin{multline*}
-\Phi_0(\vec x) =
- \begin{pmatrix}
- \varphi_0(\vec x) \\0
- \end{pmatrix},
- \qquad
- \Phi_1(\vec x) =
- \begin{pmatrix}
- 0 \\\varphi_0(\vec x)
- \end{pmatrix},
- \\
- \Phi_2(\vec x) =
- \begin{pmatrix}
- \varphi_1(\vec x) \\0
- \end{pmatrix},
- \qquad
- \Phi_3(\vec x) =
- \begin{pmatrix}
- 0 \\\varphi_1(\vec x)
- \end{pmatrix},
- \ldots
-\end{multline*}
- -->
-<IMG
- WIDTH="522" HEIGHT="94" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img24.png"
- ALT="\begin{multline*}
-\Phi_0(\vec x) =
-\begin{pmatrix}
-\varphi_0(\vec x) 0
-\...
-... \begin{pmatrix}
-0 \varphi_1(\vec x)
-\end{pmatrix},
-\ldots
-\end{multline*}"></DIV>
-<BR CLEAR="ALL">
-<P><P></P>
-<DIV ALIGN="CENTER">
-</DIV>
-where here
-<!-- MATH
- \begin{displaymath}
-comp(0)=0, \quad comp(1)=1, \quad comp(2)=0, \quad comp(3)=1, \quad \ldots
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="459" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img25.png"
- ALT="$\displaystyle comp(0)=0, \quad comp(1)=1, \quad comp(2)=0, \quad comp(3)=1, \quad \ldots
-$">
-</DIV><P></P>
-<!-- MATH
- \begin{displaymath}
-base(0)=0, \quad base(1)=0, \quad base(2)=1, \quad base(3)=1, \quad \ldots
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="433" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img26.png"
- ALT="$\displaystyle base(0)=0, \quad base(1)=0, \quad base(2)=1, \quad base(3)=1, \quad \ldots
-$">
-</DIV><P></P>
-
-<P>
-In all but very rare cases, you will not need to know which shape function
-<!-- MATH
- $\varphi_{base(i)}$
- -->
-<IMG
- WIDTH="54" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img27.png"
- ALT="$ \varphi_{base(i)}$"> of the scalar element belongs to a shape function <IMG
- WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img18.png"
- ALT="$ \Phi_i$">
-of the vector element. Let us therefore define
-<!-- MATH
- \begin{displaymath}
-\phi_i = \varphi_{base(i)}
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="90" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img28.png"
- ALT="$\displaystyle \phi_i = \varphi_{base(i)}
-$">
-</DIV><P></P>
-by which we can write the vector shape function as
-<!-- MATH
- \begin{displaymath}
-\Phi_i(\vec x) = \phi_{i}(\vec x)\ \vec e_{comp(i)}.
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="164" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img29.png"
- ALT="$\displaystyle \Phi_i(\vec x) = \phi_{i}(\vec x) \vec e_{comp(i)}.
-$">
-</DIV><P></P>
-You can now safely forget about the function <IMG
- WIDTH="52" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img20.png"
- ALT="$ base(i)$">, at least for the rest
-of this example program.
-
-<P>
-Now using this vector shape functions, we can write the discrete finite
-element solution as
-<!-- MATH
- \begin{displaymath}
-\vec u_h(\vec x) =
- \sum_i \Phi_i(\vec x)\ u_i
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="150" HEIGHT="48" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img30.png"
- ALT="$\displaystyle \vec u_h(\vec x) =
-\sum_i \Phi_i(\vec x) u_i
-$">
-</DIV><P></P>
-with scalar coefficients <IMG
- WIDTH="19" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img31.png"
- ALT="$ u_i$">. If we define an analog function <IMG
- WIDTH="22" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img32.png"
- ALT="$ \vec v_h$"> as
-test function, we can write the discrete problem as follows: Find coefficients
-<IMG
- WIDTH="19" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img31.png"
- ALT="$ u_i$"> such that
-<!-- MATH
- \begin{displaymath}
-a(\vec u_h, \vec v_h) = (\vec f, \vec v_h)
- \qquad
- \forall \vec v_h.
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="197" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img33.png"
- ALT="$\displaystyle a(\vec u_h, \vec v_h) = (\vec f, \vec v_h)
-\qquad
-\forall \vec v_h.
-$">
-</DIV><P></P>
-
-<P>
-If we insert the definition of the bilinear form and the representation of
-<IMG
- WIDTH="22" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img34.png"
- ALT="$ \vec u_h$"> and <IMG
- WIDTH="22" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img32.png"
- ALT="$ \vec v_h$"> into this formula:
-<DIV ALIGN="CENTER">
-</DIV><P></P>
-<DIV ALIGN="CENTER"><!-- MATH
- \begin{multline*}
-\sum_{i,j}
- u_i v_j
- \sum_{k,l}
- \left\{
- \left(
- \lambda \partial_l (\Phi_i)_l, \partial_k (\Phi_j)_k
- \right)_\Omega
- +
- \left(
- \mu \partial_l (\Phi_i)_k, \partial_l (\Phi_j)_k
- \right)_\Omega
- +
- \left(
- \mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l
- \right)_\Omega
- \right\}
-\\
-=
- \sum_j v_j
- \sum_l
- \left(
- f_l,
- (\Phi_j)_l
- \right)_\Omega.
-\end{multline*}
- -->
-<IMG
- WIDTH="560" HEIGHT="97" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img35.png"
- ALT="\begin{multline*}
-\sum_{i,j}
-u_i v_j
-\sum_{k,l}
-\left\{
-\left(
-\lambda \pa...
-...=
-\sum_j v_j
-\sum_l
-\left(
-f_l,
-(\Phi_j)_l
-\right)_\Omega.
-\end{multline*}"></DIV>
-<BR CLEAR="ALL">
-<P><P></P>
-<DIV ALIGN="CENTER">
-</DIV>
-We note that here and in the following, the indices <IMG
- WIDTH="25" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img36.png"
- ALT="$ k,l$"> run over spatial
-directions, i.e. <!-- MATH
- $0\le k,l < d$
- -->
-<IMG
- WIDTH="84" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img37.png"
- ALT="$ 0\le k,l < d$">, and that indices <IMG
- WIDTH="24" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img38.png"
- ALT="$ i,j$"> run over degrees
-of freedoms.
-
-<P>
-The local stiffness matrix on cell <IMG
- WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img39.png"
- ALT="$ K$"> therefore has the following entries:
-<!-- MATH
- \begin{displaymath}
-A^K_{ij}
- =
- \sum_{k,l}
- \left\{
- \left(
- \lambda \partial_l (\Phi_i)_l, \partial_k (\Phi_j)_k
- \right)_K
- +
- \left(
- \mu \partial_l (\Phi_i)_k, \partial_l (\Phi_j)_k
- \right)_K
- +
- \left(
- \mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l
- \right)_K
- \right\},
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="567" HEIGHT="54" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img40.png"
- ALT="$\displaystyle A^K_{ij}
-=
-\sum_{k,l}
-\left\{
-\left(
-\lambda \partial_l (\Phi_i)_...
-...
-+
-\left(
-\mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l
-\right)_K
-\right\},
-$">
-</DIV><P></P>
-where <IMG
- WIDTH="24" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img38.png"
- ALT="$ i,j$"> now are local degrees of freedom and therefore <!-- MATH
- $0\le i,j < N$
- -->
-<IMG
- WIDTH="89" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img41.png"
- ALT="$ 0\le i,j < N$">.
-In these formulas, we always take some component of the vector shape functions
-<IMG
- WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img18.png"
- ALT="$ \Phi_i$">, which are of course given as follows (see their definition):
-<!-- MATH
- \begin{displaymath}
-(\Phi_i)_l = \phi_i \delta_{l,comp(i)},
-\end{displaymath}
- -->
-<P></P>
-<DIV ALIGN="CENTER">
-<IMG
- WIDTH="139" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img42.png"
- ALT="$\displaystyle (\Phi_i)_l = \phi_i \delta_{l,comp(i)},
-$">
-</DIV><P></P>
-with the Kronecker symbol <!-- MATH
- $\delta_{nm}$
- -->
-<IMG
- WIDTH="31" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img43.png"
- ALT="$ \delta_{nm}$">. Due to this, we can delete some of
-the sums over <IMG
- WIDTH="13" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img44.png"
- ALT="$ k$"> and <IMG
- WIDTH="9" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img16.png"
- ALT="$ l$">:
-<DIV ALIGN="CENTER">
-</DIV><P></P>
-<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE">
-<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="28" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img45.png"
- ALT="$\displaystyle A^K_{ij}$"></TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="293" HEIGHT="53" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img46.png"
- ALT="$\displaystyle = \sum_{k,l} \Bigl\{ \left( \lambda \partial_l \phi_i \delta_{l,comp(i)}, \partial_k \phi_j \delta_{k,comp(j)} \right)_K$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-<TR VALIGN="MIDDLE">
-<TD> </TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="587" HEIGHT="45" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img47.png"
- ALT="$\displaystyle \qquad\qquad + \left( \mu \partial_l \phi_i \delta_{k,comp(i)}, ...
-..._i \delta_{k,comp(i)}, \partial_k \phi_j \delta_{l,comp(j)} \right)_K \Bigr\}$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-<TR VALIGN="MIDDLE">
-<TD> </TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="656" HEIGHT="50" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img48.png"
- ALT="$\displaystyle = \left( \lambda \partial_{comp(i)} \phi_i, \partial_{comp(j)} \p...
-...j)} + \left( \mu \partial_{comp(j)} \phi_i, \partial_{comp(i)} \phi_j \right)_K$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-<TR VALIGN="MIDDLE">
-<TD> </TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="638" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img49.png"
- ALT="$\displaystyle = \left( \lambda \partial_{comp(i)} \phi_i, \partial_{comp(j)} \p...
-...)} + \left( \mu \partial_{comp(j)} \phi_i, \partial_{comp(i)} \phi_j \right)_K.$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-</TABLE></DIV>
-<BR CLEAR="ALL"><P></P>
-<DIV ALIGN="CENTER">
-</DIV>
-
-<P>
-Likewise, the contribution of cell <IMG
- WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img39.png"
- ALT="$ K$"> to the right hand side vector is
-<DIV ALIGN="CENTER">
-</DIV><P></P>
-<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE">
-<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="26" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img50.png"
- ALT="$\displaystyle f^K_j$"></TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="126" HEIGHT="49" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img51.png"
- ALT="$\displaystyle = \sum_l \left( f_l, (\Phi_j)_l \right)_K$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-<TR VALIGN="MIDDLE">
-<TD> </TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="171" HEIGHT="49" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img52.png"
- ALT="$\displaystyle = \sum_l \left( f_l, \phi_j \delta_{l,comp(j)} \right)_K$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-<TR VALIGN="MIDDLE">
-<TD> </TD>
-<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="132" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img53.png"
- ALT="$\displaystyle = \left( f_{comp(j)}, \phi_j \right)_K.$"></TD>
-<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
- </TD></TR>
-</TABLE></DIV>
-<BR CLEAR="ALL"><P></P>
-<DIV ALIGN="CENTER">
-</DIV>
-
-<P>
-This is the form in which we will implement the local stiffness matrix and
-right hand side vectors.
-
-<P>
-As a final note: in the step-17 example program, we will revisit the elastic
-problem laid out here, and will show how to solve it in parallel on a cluster
-of computers. The resulting program will thus be able to solve this problem to
-significantly higher accuracy, and more efficiently if this is required.