<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
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:
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+<!-- MATH
+ \begin{displaymath}
-
-\partial_j (c_{ijkl} \partial_k u_l)
+ \partial_j (c_{ijkl} \partial_k u_l)
=
- 0,
+ f_i,
\qquad
i=1\ldots d,
-\end{displaymath} -->
-
-
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="237" HEIGHT="29"
- SRC="step-8.data/intro/img1.gif"
- ALT="\begin{displaymath}-
+ WIDTH="248" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img1.png"
+ ALT="$\displaystyle -
\partial_j (c_{ijkl} \partial_k u_l)
=
-0,
+f_i,
\qquad
i=1\ldots d,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-where the values <I>c</I><SUB><I>ijkl</I></SUB> are the stiffness coefficients and
+$">
+</DIV><P></P>
+where the values <IMG
+ WIDTH="33" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="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/img2.gif"
- ALT="$\lambda$">
-and <IMG
- WIDTH="13" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img3.gif"
- ALT="$\mu$">
-the coefficient tensor reduces to
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+ SRC="intro/img3.png"
+ ALT="$ \lambda$"> and <IMG
+ WIDTH="14" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="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} -->
-
-
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="236" HEIGHT="29"
- SRC="step-8.data/intro/img4.gif"
- ALT="\begin{displaymath}c_{ijkl}
+ WIDTH="241" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img5.png"
+ ALT="$\displaystyle c_{ijkl}
=
\lambda \delta_{ij} \delta_{kl} +
\mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}).
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
+$">
+</DIV><P></P>
+
<P>
The elastic equations can then be rewritten in much simpler a form:
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+<!-- MATH
+ \begin{displaymath}
-
-\nabla (\lambda+\mu) (\div \vec u)
+ \nabla \lambda (\div\vec u)
-
(\nabla \cdot \mu \nabla) \vec u
+ -
+ \div\mu (\nabla \vec u)^T
=
- 0,
-\end{displaymath} -->
-
-
+ \vec f,
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="246" HEIGHT="28"
- SRC="step-8.data/intro/img5.gif"
- ALT="\begin{displaymath}-
-\nabla (\lambda+\mu) (\div \vec u)
+ WIDTH="309" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img6.png"
+ ALT="$\displaystyle -
+\nabla \lambda (\div\vec u)
-
(\nabla \cdot \mu \nabla) \vec u
+-
+\div\mu (\nabla \vec u)^T
=
-0,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
+\vec f,
+$">
+</DIV><P></P>
and the respective bilinear form is then
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+<!-- MATH
+ \begin{displaymath}
a(\vec u, \vec v) =
-\left(
- (\lambda+\mu) \div \vec u, \div \vec v
+ \left(
+ \lambda \div\vec u, \div\vec v
\right)_\Omega
+
- \sum_k
+ \sum_{i,j}
\left(
- \mu \nabla u_k, \nabla v_k
+ \mu \partial_i u_j, \partial_i v_j
\right)_\Omega,
-\end{displaymath} -->
-
-
+ +
+ \sum_{i,j}
+ \left(
+ \mu \partial_i u_j, \partial_j v_i
+ \right)_\Omega,
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="369" HEIGHT="45"
- SRC="step-8.data/intro/img6.gif"
- ALT="\begin{displaymath}a(\vec u, \vec v) =
+ WIDTH="477" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img7.png"
+ ALT="$\displaystyle a(\vec u, \vec v) =
\left(
-(\lambda+\mu) \div \vec u, \div ...
-...+
-\sum_k
+\lambda \div\vec u, \div\vec v
+\right)...
+...\Omega,
++
+\sum_{i,j}
\left(
-\mu \nabla u_k, \nabla v_k
+\mu \partial_i u_j, \partial_j v_i
\right)_\Omega,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-or as a sum over components:
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+$">
+</DIV><P></P>
+or also writing the first term a sum over components:
+<!-- MATH
+ \begin{displaymath}
a(\vec u, \vec v) =
-\sum_{k,l}
+ \sum_{i,j}
\left(
- (\lambda+\mu) \partial_l u_l, \partial_k v_k
+ \lambda \partial_l u_l, \partial_k v_k
\right)_\Omega
+
\sum_{k,l}
\left(
- \mu \partial_l u_k, \partial_l v_k
+ \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} -->
-
-
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="383" HEIGHT="47"
- SRC="step-8.data/intro/img7.gif"
- ALT="\begin{displaymath}a(\vec u, \vec v) =
-\sum_{k,l}
+ WIDTH="492" HEIGHT="53" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img8.png"
+ ALT="$\displaystyle a(\vec u, \vec v) =
+\sum_{i,j}
\left(
-(\lambda+\mu) \parti...
-...}
+\lambda \partial_l u_l, \pa...
+...\Omega,
++
+\sum_{i,j}
\left(
-\mu \partial_l u_k, \partial_l v_k
+\mu \partial_i u_j, \partial_j v_i
\right)_\Omega.
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
+$">
+</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 <I>n</I> be the number of shape functions for the scalar finite element of
+let <IMG
+ WIDTH="14" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ SRC="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 <I>n</I>=4 in two space dimensions). Further, let <I>N</I> be the number of
-shape functions for the vector element; in two space dimensions, we need <I>n</I>shape functions for each component of the vector, so <I>N</I>=2<I>n</I>. Then, the <I>i</I>th
+already, and <IMG
+ WIDTH="43" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img10.png"
+ ALT="$ n=4$"> in two space dimensions). Further, let <IMG
+ WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ SRC="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="intro/img9.png"
+ ALT="$ n$">
+shape functions for each component of the vector, so <IMG
+ WIDTH="57" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img12.png"
+ ALT="$ N=2n$">. Then, the <IMG
+ WIDTH="10" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img13.png"
+ ALT="$ i$">th
shape function of the vector element has the form
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+<!-- MATH
+ \begin{displaymath}
\Phi_i(\vec x) = \varphi_{base(i)}(\vec x)\ \vec e_{comp(i)},
-\end{displaymath} -->
-
-
-<IMG
- WIDTH="194" HEIGHT="30"
- SRC="step-8.data/intro/img8.gif"
- ALT="\begin{displaymath}\Phi_i(\vec x) = \varphi_{base(i)}(\vec x)\ \vec e_{comp(i)},
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
+\end{displaymath}
+ -->
<P></P>
-where <I>e</I><SUB><I>l</I></SUB> is the <I>l</I>th unit vector, <I>comp</I>(<I>i</I>) is the function that tells
+<DIV ALIGN="CENTER">
+<IMG
+ WIDTH="200" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="intro/img15.png"
+ ALT="$ e_l$"> is the <IMG
+ WIDTH="9" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img16.png"
+ ALT="$ l$">th unit vector, <IMG
+ WIDTH="58" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img17.png"
+ ALT="$ comp(i)$"> is the function that tells
us which component of <IMG
- WIDTH="20" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img9.gif"
- ALT="$\Phi_i$">
-is the one that is nonzero (for
+ WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ SRC="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)$ -->
+zero). <!-- MATH
+ $\varphi_{base(i)}(x)$
+ -->
+<IMG
+ WIDTH="76" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="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="intro/img13.png"
+ ALT="$ i$"> is in the range <!-- MATH
+ $0,\ldots,N-1$
+ -->
<IMG
- WIDTH="75" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img10.gif"
- ALT="$\varphi_{base(i)}(x)$">
-describes the space dependence of the shape
-function, which is taken to be the <I>base</I>(<I>i</I>)-th shape function of the scalar
-element. Of course, while <I>i</I> is in the range
-<!-- MATH: $0,\ldots,N-1$ -->
+ WIDTH="89" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img21.png"
+ ALT="$ 0,\ldots,N-1$">, the functions
<IMG
- WIDTH="91" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img11.gif"
- ALT="$0,\ldots,N-1$">,
-the functions
-<I>comp</I>(<I>i</I>) and <I>base</I>(<I>i</I>) have the ranges 0,1 (in 2D) and
-<!-- MATH: $0,\ldots,n-1$ -->
+ WIDTH="58" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img17.png"
+ ALT="$ comp(i)$"> and <IMG
+ WIDTH="52" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img20.png"
+ ALT="$ base(i)$"> have the ranges <IMG
+ WIDTH="27" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img22.png"
+ ALT="$ 0,1$"> (in 2D) and <!-- MATH
+ $0,\ldots,n-1$
+ -->
<IMG
- WIDTH="86" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img12.gif"
- ALT="$0,\ldots,n-1$">,
+ WIDTH="84" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img23.png"
+ ALT="$ 0,\ldots,n-1$">,
respectively.
<P>
For example (but this layout is not mandatory, and you should not rely on it),
the following layout could be used by the library:
<DIV ALIGN="CENTER">
-<BR>
+</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="519" HEIGHT="84" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img13.gif"
- ALT="\begin{multline*}\Phi_0(\vec x) =
+ WIDTH="522" HEIGHT="94" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img24.png"
+ ALT="\begin{multline*}
+\Phi_0(\vec x) =
\begin{pmatrix}
-\varphi_0(\vec x) \\ 0
-\en...
+\varphi_0(\vec x) 0
+\...
... \begin{pmatrix}
-0 \\ \varphi_1(\vec x)
+0 \varphi_1(\vec x)
\end{pmatrix},
\ldots
-\end{multline*}">
-<BR></DIV>
-where here
-<BR><P></P>
+\end{multline*}"></DIV>
+<BR CLEAR="ALL">
+<P><P></P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-comp(0)=0, \quad comp(1)=1, \quad comp(2)=0, \quad comp(3)=1, \quad \ldots
-\end{displaymath} -->
-
-
-<IMG
- WIDTH="455" HEIGHT="28"
- SRC="step-8.data/intro/img14.gif"
- ALT="\begin{displaymath}comp(0)=0, \quad comp(1)=1, \quad comp(2)=0, \quad comp(3)=1, \quad \ldots
-\end{displaymath}">
</DIV>
-<BR CLEAR="ALL">
+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>
-<BR><P></P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-base(0)=0, \quad base(1)=0, \quad base(2)=1, \quad base(3)=1, \quad \ldots
-\end{displaymath} -->
-
-
<IMG
- WIDTH="429" HEIGHT="28"
- SRC="step-8.data/intro/img15.gif"
- ALT="\begin{displaymath}base(0)=0, \quad base(1)=0, \quad base(2)=1, \quad base(3)=1, \quad \ldots
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
+ WIDTH="459" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="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)}$ -->
+<!-- MATH
+ $\varphi_{base(i)}$
+ -->
<IMG
- WIDTH="53" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img16.gif"
- ALT="$\varphi_{base(i)}$">
-of the scalar element belongs to a shape function <IMG
- WIDTH="20" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img9.gif"
- ALT="$\Phi_i$">of the vector element. Let us therefore define
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+ WIDTH="54" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="intro/img18.png"
+ ALT="$ \Phi_i$">
+of the vector element. Let us therefore define
+<!-- MATH
+ \begin{displaymath}
\phi_i = \varphi_{base(i)}
-\end{displaymath} -->
-
-
-<IMG
- WIDTH="84" HEIGHT="30"
- SRC="step-8.data/intro/img17.gif"
- ALT="\begin{displaymath}\phi_i = \varphi_{base(i)}
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
+\end{displaymath}
+ -->
<P></P>
-by which we can write the vector shape function as
-<BR><P></P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-\Phi_i(\vec x) = \phi_{i}(\vec x)\ \vec e_{comp(i)}.
-\end{displaymath} -->
-
-
<IMG
- WIDTH="158" HEIGHT="30"
- SRC="step-8.data/intro/img18.gif"
- ALT="\begin{displaymath}\Phi_i(\vec x) = \phi_{i}(\vec x)\ \vec e_{comp(i)}.
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
+ WIDTH="90" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ SRC="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>
-You can now safely forget about the function <I>base</I>(<I>i</I>), at least for the rest
+<DIV ALIGN="CENTER">
+<IMG
+ WIDTH="164" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="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
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+<!-- MATH
+ \begin{displaymath}
\vec u_h(\vec x) =
-\sum_i \Phi_i(\vec x)\ u_i
-
-\end{displaymath} -->
-
-
-<IMG
- WIDTH="144" HEIGHT="44"
- SRC="step-8.data/intro/img19.gif"
- ALT="\begin{displaymath}\vec u_h(\vec x) =
-\sum_i \Phi_i(\vec x)\ u_i
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
+ \sum_i \Phi_i(\vec x)\ u_i
+\end{displaymath}
+ -->
<P></P>
-with scalar coefficients <I>u</I><SUB><I>i</I></SUB>. If we define an analog function <IMG
- WIDTH="21" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img20.gif"
- ALT="$\vec v_h$">
-as
-test function, we can write the discrete problem as follows: Find coefficients
-<I>u</I><SUB><I>i</I></SUB> such that
-<BR><P></P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+<IMG
+ WIDTH="150" HEIGHT="48" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="intro/img31.png"
+ ALT="$ u_i$">. If we define an analog function <IMG
+ WIDTH="22" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="intro/img31.png"
+ ALT="$ u_i$"> such that
+<!-- MATH
+ \begin{displaymath}
a(\vec u_h, \vec v_h) = (\vec f, \vec v_h)
-\qquad
+ \qquad
\forall \vec v_h.
-\end{displaymath} -->
-
-
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="193" HEIGHT="28"
- SRC="step-8.data/intro/img21.gif"
- ALT="\begin{displaymath}a(\vec u_h, \vec v_h) = (\vec f, \vec v_h)
+ WIDTH="197" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img33.png"
+ ALT="$\displaystyle a(\vec u_h, \vec v_h) = (\vec f, \vec v_h)
\qquad
\forall \vec v_h.
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
+$">
+</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/img22.gif"
- ALT="$\vec u_h$">
-and <IMG
- WIDTH="21" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img20.gif"
- ALT="$\vec v_h$">
-into this formula:
+ SRC="intro/img34.png"
+ ALT="$ \vec u_h$"> and <IMG
+ WIDTH="22" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img32.png"
+ ALT="$ \vec v_h$"> into this formula:
<DIV ALIGN="CENTER">
-<BR>
+</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="520" HEIGHT="83" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img23.gif"
- ALT="\begin{multline*}\sum_{i,j}
+ WIDTH="560" HEIGHT="97" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img35.png"
+ ALT="\begin{multline*}
+\sum_{i,j}
u_i v_j
\sum_{k,l}
\left\{
\left(
-(\lambda+\mu)...
+\lambda \pa...
...=
\sum_j v_j
\sum_l
f_l,
(\Phi_j)_l
\right)_\Omega.
-\end{multline*}">
-<BR></DIV>
-We note that here and in the following, the indices <I>k</I>,<I>l</I> run over spatial
-directions, i.e.
-<!-- MATH: $0\le k,l \le d-1$ -->
+\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="intro/img36.png"
+ ALT="$ k,l$"> run over spatial
+directions, i.e. <!-- MATH
+ $0\le k,l < d$
+ -->
<IMG
- WIDTH="110" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img24.gif"
- ALT="$0\le k,l \le d-1$">,
-and that indices <I>i</I>,<I>j</I> run over degrees
+ WIDTH="84" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img37.png"
+ ALT="$ 0\le k,l < d$">, and that indices <IMG
+ WIDTH="24" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img38.png"
+ ALT="$ i,j$"> run over degrees
of freedoms.
<P>
-The local stiffness matrix on cell <I>K</I> therefore has the following entries:
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+The local stiffness matrix on cell <IMG
+ WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img39.png"
+ ALT="$ K$"> therefore has the following entries:
+<!-- MATH
+ \begin{displaymath}
A^K_{ij}
-=
+ =
\sum_{k,l}
\left\{
\left(
- (\lambda+\mu) \partial_l (\Phi_i)_l, \partial_k (\Phi_j)_k
+ \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} -->
-
-
+\end{displaymath}
+ -->
+<P></P>
+<DIV ALIGN="CENTER">
<IMG
- WIDTH="443" HEIGHT="47"
- SRC="step-8.data/intro/img25.gif"
- ALT="\begin{displaymath}A^K_{ij}
+ WIDTH="567" HEIGHT="54" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img40.png"
+ ALT="$\displaystyle A^K_{ij}
=
\sum_{k,l}
\left\{
\left(
-(\lambda+\mu) \part...
-...tial_l (\Phi_i)_k, \partial_l (\Phi_j)_k
+\lambda \partial_l (\Phi_i)_...
+...
++
+\left(
+\mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l
\right)_K
\right\},
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-where <I>i</I>,<I>j</I> now are local degrees of freedom and therefore
-<!-- MATH: $0\le i,j \le N-1$ -->
+$">
+</DIV><P></P>
+where <IMG
+ WIDTH="24" HEIGHT="28" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img38.png"
+ ALT="$ i,j$"> now are local degrees of freedom and therefore <!-- MATH
+ $0\le i,j < N$
+ -->
<IMG
- WIDTH="116" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img26.gif"
- ALT="$0\le i,j \le N-1$">.
+ WIDTH="89" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ SRC="intro/img41.png"
+ ALT="$ 0\le i,j < N$">.
In these formulas, we always take some component of the vector shape functions
<IMG
- WIDTH="20" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img9.gif"
- ALT="$\Phi_i$">,
-which are of course given as follows (see their definition):
-<BR><P></P>
-<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
+ WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ SRC="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} -->
-
-
-<IMG
- WIDTH="134" HEIGHT="30"
- SRC="step-8.data/intro/img27.gif"
- ALT="\begin{displaymath}(\Phi_i)_l = \phi_i \delta_{l,comp(i)},
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
+\end{displaymath}
+ -->
<P></P>
-with the Kronecker symbol
-<!-- MATH: $\delta_{nm}$ -->
-<IMG
- WIDTH="30" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
- SRC="step-8.data/intro/img28.gif"
- ALT="$\delta_{nm}$">.
-Due to this, we can delete some of
-the sums over <I>k</I> and <I>l</I>:
<DIV ALIGN="CENTER">
-<BR>
<IMG
- WIDTH="521" HEIGHT="147" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img29.gif"
- ALT="\begin{align*}A^K_{ij}
-&=
-\sum_{k,l}
-\Bigl\{
-\left(
-(\lambda+\mu) \partial_...
-...abla \phi_i,
-\nabla \phi_j
-\right)_K
-\ \delta_{comp(i),comp(j)}.
-\end{align*}">
-<BR></DIV>
+ WIDTH="139" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="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="intro/img44.png"
+ ALT="$ k$"> and <IMG
+ WIDTH="9" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="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="intro/img45.png"
+ ALT="$\displaystyle A^K_{ij}$"></TD>
+<TD NOWRAP ALIGN="LEFT"><IMG
+ WIDTH="293" HEIGHT="53" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="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="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="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 <I>K</I> to the right hand side vector is
+Likewise, the contribution of cell <IMG
+ WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ SRC="intro/img39.png"
+ ALT="$ K$"> to the right hand side vector is
<DIV ALIGN="CENTER">
-<BR>
-<IMG
- WIDTH="191" HEIGHT="106" ALIGN="BOTTOM" BORDER="0"
- SRC="step-8.data/intro/img30.gif"
- ALT="\begin{align*}f^K_j
-&=
-\sum_l
-\left(
-f_l,
-(\Phi_j)_l
-\right)_K
-\\
-&=
-\s...
-...mp(j)}
-\right)_K
-\\
-&=
-\left(
-f_{comp(j)},
-\phi_j
-\right)_K.
-\end{align*}">
-<BR></DIV>
+</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="intro/img50.png"
+ ALT="$\displaystyle f^K_j$"></TD>
+<TD NOWRAP ALIGN="LEFT"><IMG
+ WIDTH="126" HEIGHT="49" ALIGN="MIDDLE" BORDER="0"
+ SRC="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="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="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.