-<H3><A NAME="SECTION00001000000000000000">
+<H2><A NAME="SECTION00001000000000000000">
Overview</A>
-</H3>
-This example is devoted to the <EM>discontinuous Galerkin method</EM>, or
+</H2>
+This example is devoted to the <I>discontinuous Galerkin method</I>, or
in short: DG method. It includes the following topics.
+
<UL>
<LI>Discretization of the linear transport equation with the DG method
+</LI>
<LI>Two different assembling routines for the system matrix based on
face terms given as a sum of integrals that
-<DL COMPACT>
-<DT>1.
-<DD>loops over all cell and all their faces, or that
-<DT>2.
-<DD>loops over all faces, whereas each face is treated only once.
-</DL>
+
+<OL>
+<LI>loops over all cell and all their faces, or that
+</LI>
+<LI>loops over all faces, whereas each face is treated only once.
+</LI>
+</OL>
+</LI>
<LI>Time comparison of the two assembling routines.
+</LI>
</UL>
+
<P>
-<H3><A NAME="SECTION00002000000000000000">
+<H2><A NAME="SECTION00002000000000000000">
Problem</A>
-</H3>
+</H2>
The DG method was first introduced to discretize simple transport
equations. Over the past years DG methods have been applied to a
variety of problems and many different schemes were introduced
this example's purpose is to illustrate some implementational issues
of the DG discretization only, here we simply consider the linear
transport equation
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
-\nabla\cdot \left\{\beta u\right\}=f \qquad\mbox{in }\Omega,
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="transport-equation"> </A><IMG
- WIDTH="170" HEIGHT="31"
- SRC="step-12.data/img2.gif"
- ALT="\begin{displaymath}
-\nabla\cdot \left\{\beta u\right\}=f \qquad\mbox{in }\Omega,
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="transport-equation"></A><!-- MATH
+ \begin{equation}
+\nabla\cdot \left\{\boldsymbol\beta u\right\}=f \qquad\mbox{in }\Omega,
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="105" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img1.png"
+ ALT="$\displaystyle \nabla\cdot \left\{\boldsymbol\beta u\right\}=f$"> in <IMG
+ WIDTH="22" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img2.png"
+ ALT="$\displaystyle \Omega,$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(1)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
subject to the boundary conditions
-<BR><P></P>
+<P><!-- MATH
+ \begin{displaymath}
+u=g\quad\mbox{on }\Gamma_-,
+\end{displaymath}
+ -->
+</P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-u=g\quad\mbox{on }\partial\Omega_-.
-\end{displaymath} -->
-
-
<IMG
- WIDTH="122" HEIGHT="31"
- SRC="step-12.data/img3.gif"
- ALT="\begin{displaymath}u=g\quad\mbox{on }\partial\Omega_-.
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-Here,
-<!-- MATH: $\beta=\beta(x)$ -->
-<IMG
- WIDTH="73" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img4.gif"
- ALT="$\beta=\beta(x)$">
-is a vector field, <I>f</I> a source function, <I>u</I> the
-(scalar) solution function, <I>g</I> a boundary value function and
-<BR><P></P>
+ WIDTH="46" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img3.png"
+ ALT="$\displaystyle u=g$"> on <IMG
+ WIDTH="31" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img4.png"
+ ALT="$\displaystyle \Gamma_-,
+$">
+</DIV><P>
+</P>
+on the inflow part <IMG
+ WIDTH="27" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img5.png"
+ ALT="$ \Gamma_-$">
+ of the boundary <!-- MATH
+ $\Gamma=\partial\Omega$
+ -->
+<IMG
+ WIDTH="61" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img6.png"
+ ALT="$ \Gamma=\partial\Omega$">
+
+of the domain. Here, <!-- MATH
+ $\boldsymbol\beta=\boldsymbol\beta(x)$
+ -->
+<IMG
+ WIDTH="75" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img7.png"
+ ALT="$ \boldsymbol\beta=\boldsymbol\beta(x)$">
+ denotes a
+vector field, <IMG
+ WIDTH="15" HEIGHT="33" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img8.png"
+ ALT="$ f$">
+ a source function, <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img9.png"
+ ALT="$ u$">
+ the (scalar) solution
+function, <IMG
+ WIDTH="13" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img10.png"
+ ALT="$ g$">
+ a boundary value function,
+<P><!-- MATH
+ \begin{displaymath}
+\Gamma_-:=\{x\in\Gamma, \boldsymbol\beta(x)\cdot{\bf n}(x)<0\}
+\end{displaymath}
+ -->
+</P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-\partial\Omega_-:=\{x\in\partial\Omega, \beta(x)\cdot n(x)<0\}
-\end{displaymath} -->
-
-
<IMG
- WIDTH="254" HEIGHT="31"
- SRC="step-12.data/img5.gif"
- ALT="\begin{displaymath}\partial\Omega_-:=\{x\in\partial\Omega, \beta(x)\cdot n(x)<0\}
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-the inflow part of the boundary of the domain, whereas <I>n</I> denotes the
-unit outer normal to the boundary
-<!-- MATH: $\partial\Omega$ -->
-<IMG
- WIDTH="27" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img6.gif"
- ALT="$\partial\Omega$">.
-This
-equation, (<A HREF="step-12.html#transport-equation">1</A>), is the conservative version of
-the transport equation already considered in <A HREF="step-9.html#SECTION00000010000000000000">step 9</A> of this tutorial.
+ WIDTH="237" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img11.png"
+ ALT="$\displaystyle \Gamma_-:=\{x\in\Gamma, \boldsymbol\beta(x)\cdot{\bf n}(x)<0\}
+$">
+</DIV><P>
+</P>
+the inflow part of the boundary of the domain and <IMG
+ WIDTH="15" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img12.png"
+ ALT="$ {\bf n}$">
+ denotes
+the unit outward normal to the boundary <IMG
+ WIDTH="15" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img13.png"
+ ALT="$ \Gamma$">
+. Equation
+(<A HREF="#transport-equation">1</A>) is the conservative version of the
+transport equation already considered in step 9 of this tutorial.
<P>
-In particular, we consider problem (<A HREF="step-12.html#transport-equation">1</A>) on
-
-<!-- MATH: $\Omega=[0,1]^2$ -->
-<IMG
- WIDTH="83" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img7.gif"
- ALT="$\Omega=[0,1]^2$">
-with
-<!-- MATH: $\beta=\frac{1}{|x|}(-x_2, x_1)$ -->
-<IMG
- WIDTH="127" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img8.gif"
- ALT="$\beta=\frac{1}{\vert x\vert}(-x_2, x_1)$">
-representing a
-circular counterclockwise flow field, <I>f</I>=0 and <I>g</I>=1 on
-
-<!-- MATH: $x\in\partial\Omega_-^1:=[0,0.5]\times\{0\}$ -->
-<IMG
- WIDTH="194" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img9.gif"
- ALT="$x\in\partial\Omega_-^1:=[0,0.5]\times\{0\}$">
-and <I>g</I>=0 on
-<!-- MATH: $x\in
-\partial\Omega_-\setminus \partial\Omega_-^1$ -->
-<IMG
- WIDTH="121" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img10.gif"
- ALT="$x\in
-\partial\Omega_-\setminus \partial\Omega_-^1$">.
-
-<H3><A NAME="SECTION00003000000000000000">
+In particular, we consider problem (<A HREF="#transport-equation">1</A>) on
+<!-- MATH
+ $\Omega=[0,1]^2$
+ -->
+<IMG
+ WIDTH="83" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img14.png"
+ ALT="$ \Omega=[0,1]^2$">
+ with <!-- MATH
+ $\boldsymbol\beta=\frac{1}{|x|}(-x_2, x_1)$
+ -->
+<IMG
+ WIDTH="129" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img15.png"
+ ALT="$ \boldsymbol\beta=\frac{1}{\vert x\vert}(-x_2, x_1)$">
+
+representing a circular counterclockwise flow field, <IMG
+ WIDTH="47" HEIGHT="33" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img16.png"
+ ALT="$ f=0$">
+ and <IMG
+ WIDTH="45" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img17.png"
+ ALT="$ g=1$">
+
+on <!-- MATH
+ $x\in\Gamma_-^1:=[0,0.5]\times\{0\}$
+ -->
+<IMG
+ WIDTH="182" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img18.png"
+ ALT="$ x\in\Gamma_-^1:=[0,0.5]\times\{0\}$">
+ and <IMG
+ WIDTH="45" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img19.png"
+ ALT="$ g=0$">
+ on <!-- MATH
+ $x\in
+\Gamma_-\setminus \Gamma_-^1$
+ -->
+<IMG
+ WIDTH="96" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img20.png"
+ ALT="$ x\in
+\Gamma_-\setminus \Gamma_-^1$">
+.
+
+<H2><A NAME="SECTION00003000000000000000">
Discretization</A>
-</H3>
-To follow the general paradigm of deriving DG discretizations for purely
-hyperbolic equations we first consider the general hyperbolic problem
-<BR><P></P>
+</H2>
+Following the general paradigm of deriving DG discretizations for
+purely hyperbolic equations, we first consider the general hyperbolic
+problem
+<P><!-- MATH
+ \begin{displaymath}
+\nabla\cdot {\mathcal F}(u)=f \qquad\mbox{in }\Omega,
+\end{displaymath}
+ -->
+</P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-\nabla\cdot F(u)=f \qquad\mbox{in }\Omega,
-\end{displaymath} -->
-
-
<IMG
- WIDTH="170" HEIGHT="31"
- SRC="step-12.data/img11.gif"
- ALT="\begin{displaymath}\nabla\cdot F(u)=f \qquad\mbox{in }\Omega,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-subject to appropriate boundary conditions. Here <I>F</I> denotes the flux
-function of the equation under consideration that in our case, see
-equation (<A HREF="step-12.html#transport-equation">1</A>), is represented by
-<!-- MATH: $F(u)=\beta u$ -->
-<IMG
- WIDTH="86" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img12.gif"
- ALT="$F(u)=\beta u$">.
-For deriving the DG discretization we start with a variational,
-mesh-dependent formulation of the problem,
-<BR><P></P>
+ WIDTH="103" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img21.png"
+ ALT="$\displaystyle \nabla\cdot {\mathcal F}(u)=f$"> in <IMG
+ WIDTH="22" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img2.png"
+ ALT="$\displaystyle \Omega,$">
+</DIV><P>
+</P>
+subject to appropriate boundary conditions. Here <!-- MATH
+ ${\mathcal F}$
+ -->
+<IMG
+ WIDTH="19" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img22.png"
+ ALT="$ {\mathcal F}$">
+
+denotes the flux function of the equation under consideration that in
+our case, see equation (<A HREF="#transport-equation">1</A>), is represented by
+<!-- MATH
+ ${\mathcal F}(u)=\boldsymbol\beta u$
+ -->
+<IMG
+ WIDTH="87" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img23.png"
+ ALT="$ {\mathcal F}(u)=\boldsymbol\beta u$">
+. For deriving the DG
+discretization we start with a variational, mesh-dependent
+formulation of the problem,
+<P><!-- MATH
+ \begin{displaymath}
+\sum_\kappa\left\{-({\mathcal F}(u),\nabla v)_\kappa+({\mathcal
+ F}(u)\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega,
+\end{displaymath}
+ -->
+</P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-\sum_K\left\{-(F(u),\nabla v)_K+(F(u)\cdot n, v)_{\partial K}\right\}=(f,v)_\Omega,
-\end{displaymath} -->
-
-
<IMG
- WIDTH="363" HEIGHT="50"
- SRC="step-12.data/img13.gif"
- ALT="\begin{displaymath}\sum_K\left\{-(F(u),\nabla v)_K+(F(u)\cdot n, v)_{\partial K}\right\}=(f,v)_\Omega,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-that originates from (<A HREF="step-12.html#transport-equation">1</A>) by multiplication with
-a test function and integration by parts on each cell <I>K</I> of the
-triangulation. Here
-<!-- MATH: $(\cdot, \cdot)_K$ -->
-<IMG
- WIDTH="49" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img14.gif"
- ALT="$(\cdot, \cdot)_K$">
-and
-<!-- MATH: $(\cdot, \cdot)_{\partial
-K}$ -->
-<IMG
- WIDTH="57" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img15.gif"
- ALT="$(\cdot, \cdot)_{\partial
-K}$">
-simply denote the integrals over the cell and the boundary of
-the cell, respectively. To discretize the problem, the functions <I>u</I>
-and <I>v</I> are replaced by discrete functions <I>u</I><SUB><I>h</I></SUB> and <I>v</I><SUB><I>h</I></SUB> that in the
-case of discontinuous Galerkin methods belong to the space <I>V</I><SUB><I>h</I></SUB> of
-discontinuous piecewise polynomial functions of some degree <I>p</I>. Due
-to the discontinuity of the discrete function <I>u</I><SUB><I>h</I></SUB> on interelement
-faces, the flux
-<!-- MATH: $F(u)\cdot n$ -->
-<IMG
- WIDTH="65" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img16.gif"
- ALT="$F(u)\cdot n$">
-must be replaced by a <EM>numerical
- flux</EM> function
-<!-- MATH: $H(u_h, \hat u_h, n)$ -->
-<IMG
- WIDTH="97" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img17.gif"
- ALT="$H(u_h, \hat u_h, n)$">,
-where
-<!-- MATH: $u_h|_{\partial K}$ -->
-<IMG
- WIDTH="49" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img18.gif"
- ALT="$u_h\vert _{\partial K}$">denotes the inner trace (w.r.t. the cell <I>K</I>) of <I>u</I><SUB><I>h</I></SUB> and
-<!-- MATH: $\hat
-u_h|_{\partial K}$ -->
-<IMG
- WIDTH="49" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img19.gif"
- ALT="$\hat
-u_h\vert _{\partial K}$">
-the outer trace, i.e. the value of <I>u</I><SUB><I>h</I></SUB> on the
-neighboring cell. Furthermore the numerical flux function <I>H</I>, among
-other things, must be consistent, i.e.
-<BR><P></P>
+ WIDTH="363" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img24.png"
+ ALT="$\displaystyle \sum_\kappa\left\{-({\mathcal F}(u),\nabla v)_\kappa+({\mathcal
+F}(u)\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega,
+$">
+</DIV><P>
+</P>
+that originates from (<A HREF="#transport-equation">1</A>) by multiplication with
+a test function <IMG
+ WIDTH="13" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img25.png"
+ ALT="$ v$">
+ and integration by parts on each cell <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img26.png"
+ ALT="$ \kappa$">
+ of
+the triangulation. Here <!-- MATH
+ $(\cdot, \cdot)_\kappa$
+ -->
+<IMG
+ WIDTH="44" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img27.png"
+ ALT="$ (\cdot, \cdot)_\kappa$">
+ and <!-- MATH
+ $(\cdot,
+\cdot)_{\partial\kappa}$
+ -->
+<IMG
+ WIDTH="52" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img28.png"
+ ALT="$ (\cdot,
+\cdot)_{\partial\kappa}$">
+ simply denote the integrals over the cell
+<IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img26.png"
+ ALT="$ \kappa$">
+ and the boundary <!-- MATH
+ $\partial\kappa$
+ -->
+<IMG
+ WIDTH="25" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img29.png"
+ ALT="$ \partial\kappa$">
+ of the cell,
+respectively. To discretize the problem, the functions <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img9.png"
+ ALT="$ u$">
+ and <IMG
+ WIDTH="13" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img25.png"
+ ALT="$ v$">
+ are
+replaced by discrete functions <IMG
+ WIDTH="23" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img30.png"
+ ALT="$ u_h$">
+ and <IMG
+ WIDTH="21" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img31.png"
+ ALT="$ v_h$">
+ that in the case of
+discontinuous Galerkin methods belong to the space <IMG
+ WIDTH="23" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img32.png"
+ ALT="$ V_h$">
+ of
+discontinuous piecewise polynomial functions of some degree <IMG
+ WIDTH="13" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img33.png"
+ ALT="$ p$">
+. Due
+to the discontinuity of the discrete function <IMG
+ WIDTH="23" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img30.png"
+ ALT="$ u_h$">
+ on interelement
+faces, the flux <!-- MATH
+ ${\mathcal F}(u)\cdot{\bf n}$
+ -->
+<IMG
+ WIDTH="66" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img34.png"
+ ALT="$ {\mathcal F}(u)\cdot{\bf n}$">
+ must be replaced by a
+<I>numerical flux</I> function <!-- MATH
+ ${\mathcal H}(u_h^+, u_h^-, {\bf n})$
+ -->
+<IMG
+ WIDTH="102" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img35.png"
+ ALT="$ {\mathcal H}(u_h^+, u_h^-, {\bf n})$">
+,
+where <!-- MATH
+ $u_h^+|_{\partial\kappa}$
+ -->
+<IMG
+ WIDTH="47" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img36.png"
+ ALT="$ u_h^+\vert _{\partial\kappa}$">
+ denotes the inner trace (w.r.t. the
+cell <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img26.png"
+ ALT="$ \kappa$">
+) of <IMG
+ WIDTH="23" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img30.png"
+ ALT="$ u_h$">
+ and <!-- MATH
+ $u_h^-|_{\partial\kappa}$
+ -->
+<IMG
+ WIDTH="47" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img37.png"
+ ALT="$ u_h^-\vert _{\partial\kappa}$">
+ the outer trace,
+i.e. the value of <IMG
+ WIDTH="23" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img30.png"
+ ALT="$ u_h$">
+ on the neighboring cell. Furthermore the
+numerical flux function <!-- MATH
+ ${\mathcal H}$
+ -->
+<IMG
+ WIDTH="19" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img38.png"
+ ALT="$ {\mathcal H}$">
+, among other things, must be
+consistent, i.e.
+<P><!-- MATH
+ \begin{displaymath}
+{\mathcal H}(u,u,{\bf n})={\mathcal F}(u)\cdot{\bf n},
+\end{displaymath}
+ -->
+</P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-H(u,u,n)=F(u)\cdot n,
-\end{displaymath} -->
-
-
<IMG
- WIDTH="162" HEIGHT="31"
- SRC="step-12.data/img20.gif"
- ALT="\begin{displaymath}H(u,u,n)=F(u)\cdot n,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
+ WIDTH="169" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img39.png"
+ ALT="$\displaystyle {\mathcal H}(u,u,{\bf n})={\mathcal F}(u)\cdot{\bf n},
+$">
+</DIV><P>
+</P>
and conservative, i.e.
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
-H(v,w,n)=-H(w,v,-n).
-\end{equation} -->
-<A NAME="conservative"> </A>
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP>
-<I>H</I>(<I>v</I>,<I>w</I>,<I>n</I>)=-<I>H</I>(<I>w</I>,<I>v</I>,-<I>n</I>).
-</TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="conservative"></A><!-- MATH
+ \begin{equation}
+{\mathcal H}(v,w,{\bf n})=-{\mathcal H}(w,v,-{\bf n}).
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="214" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img40.png"
+ ALT="$\displaystyle {\mathcal H}(v,w,{\bf n})=-{\mathcal H}(w,v,-{\bf n}).$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(2)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
-This yields the following <EM>discontinuous Galerkin
- discretization</EM>: find
-<!-- MATH: $u_h\in V_h$ -->
-<IMG
- WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img21.gif"
- ALT="$u_h\in V_h$">
-such that
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
-\sum_K\left\{-(F(u_h),\nabla v_h)_K+(H(u_h,\hat u_h,n), v_h)_{\partial K}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h.
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-scheme"> </A><IMG
- WIDTH="513" HEIGHT="50"
- SRC="step-12.data/img22.gif"
- ALT="\begin{displaymath}
-\sum_K\left\{-(F(u_h),\nabla v_h)_K+(H(u_h,\hat u_h,n), v_h)_{\partial K}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h.
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+This yields the following <I>discontinuous Galerkin
+ discretization</I>: find <!-- MATH
+ $u_h\in V_h$
+ -->
+<IMG
+ WIDTH="63" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img41.png"
+ ALT="$ u_h\in V_h$">
+ such that
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="dg-scheme"></A><!-- MATH
+ \begin{equation}
+\sum_\kappa\left\{-({\mathcal F}(u_h),\nabla v_h)_\kappa+({\mathcal H}(u_h^+,u_h^-,{\bf n}), v_h)_{\partial\kappa}\right\}=(f,v_h)_\Omega, \quad\forall v_h\in V_h.
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="526" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img42.png"
+ ALT="$\displaystyle \sum_\kappa\left\{-({\mathcal F}(u_h),\nabla v_h)_\kappa+({\mathc...
+...\bf n}), v_h)_{\partial\kappa}\right\}=(f,v_h)_\Omega, \quad\forall v_h\in V_h.$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(3)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
Boundary conditions are realized by replacing <IMG
- WIDTH="23" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img23.gif"
- ALT="$\hat u_h$">
-on the inflow boundary
-<!-- MATH: $\partial\Omega_-$ -->
-<IMG
- WIDTH="39" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img24.gif"
- ALT="$\partial\Omega_-$">
-by the boundary function <I>g</I>.
+ WIDTH="26" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img43.png"
+ ALT="$ u_h^-$">
+ on the inflow boundary <IMG
+ WIDTH="27" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img5.png"
+ ALT="$ \Gamma_-$">
+ by the boundary function <IMG
+ WIDTH="13" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img10.png"
+ ALT="$ g$">
+.
In the special case of the transport equation
-(<A HREF="step-12.html#transport-equation">1</A>) the numerical flux in its simplest form
+(<A HREF="#transport-equation">1</A>) the numerical flux in its simplest form
is given by
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
-H(u_h,\hat u_h,n)(x)=\left\{\begin{array}{ll}
- (\beta\cdot n\, u_h)(x),&\mbox{for } \beta(x)\cdot n(x)<0,\\
- (\beta\cdot n\, \hat u_h)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0,
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="flux-transport-equation"></A><!-- MATH
+ \begin{equation}
+{\mathcal H}(u_h^+,u_h^-,{\bf n})(x)=\left\{\begin{array}{ll}
+ (\boldsymbol\beta\cdot{\bf n}\, u_h^-)(x),&\mbox{for } \boldsymbol\beta(x)\cdot{\bf n}(x)<0,\\
+ (\boldsymbol\beta\cdot{\bf n}\, u_h^+)(x),&\mbox{for } \boldsymbol\beta(x)\cdot{\bf n}(x)\geq 0,
\end{array}
\right.
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="flux-transport-equation"> </A><IMG
- WIDTH="417" HEIGHT="50"
- SRC="step-12.data/img25.gif"
- ALT="\begin{displaymath}
-H(u_h,\hat u_h,n)(x)=\left\{\begin{array}{ll}
-(\beta\cdot ...
-...)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0,
-\end{array}\right.
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="444" HEIGHT="61" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img44.png"
+ ALT="$\displaystyle {\mathcal H}(u_h^+,u_h^-,{\bf n})(x)=\left\{\begin{array}{ll} (\b...
+...)(x),&\mbox{for } \boldsymbol\beta(x)\cdot{\bf n}(x)\geq 0, \end{array} \right.$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(4)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
where on the inflow part of the cell the value is taken from the
neighboring cell, <IMG
- WIDTH="23" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img23.gif"
- ALT="$\hat u_h$">,
-and on the outflow part the value is
-taken from the current cell, <I>u</I><SUB><I>h</I></SUB>. Hence, the discontinuous Galerkin
-scheme for the transport equation (<A HREF="step-12.html#transport-equation">1</A>) is given
-by: find
-<!-- MATH: $u_h\in V_h$ -->
-<IMG
- WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img21.gif"
- ALT="$u_h\in V_h$">
-such that for all
-<!-- MATH: $v_h\in V_h$ -->
-<IMG
- WIDTH="62" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img26.gif"
- ALT="$v_h\in V_h$">
-following
+ WIDTH="26" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img43.png"
+ ALT="$ u_h^-$">
+, and on the outflow part the value is
+taken from the current cell, <IMG
+ WIDTH="26" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img45.png"
+ ALT="$ u_h^+$">
+. Hence, the discontinuous Galerkin
+scheme for the transport equation (<A HREF="#transport-equation">1</A>) is given
+by: find <!-- MATH
+ $u_h\in V_h$
+ -->
+<IMG
+ WIDTH="63" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img41.png"
+ ALT="$ u_h\in V_h$">
+ such that for all <!-- MATH
+ $v_h\in V_h$
+ -->
+<IMG
+ WIDTH="61" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img46.png"
+ ALT="$ v_h\in V_h$">
+ following
equation holds:
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
-\sum_K\left\{-(u_h,\beta\cdot\nabla v_h)_K+(\beta\cdot n\, u_h, v_h)_{\partial K_+\setminus\partial\Omega}+(\beta\cdot n\, \hat u_h, v_h)_{\partial K_-\setminus\partial\Omega}\right\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-transport"> </A><IMG
- WIDTH="718" HEIGHT="50"
- SRC="step-12.data/img27.gif"
- ALT="\begin{displaymath}
-\sum_K\left\{-(u_h,\beta\cdot\nabla v_h)_K+(\beta\cdot n\, ...
-...ght\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="dg-transport"></A><!-- MATH
+ \begin{equation}
+\sum_\kappa\left\{-(u_h,\boldsymbol\beta\cdot\nabla v_h)_\kappa+(\boldsymbol\beta\cdot{\bf n}\, u_h, v_h)_{\partial\kappa_+\setminus\Gamma}+(\boldsymbol\beta\cdot{\bf n}\, u_h^-, v_h)_{\partial\kappa_-\setminus\Gamma}\right\}=(f,v_h)_\Omega-(\boldsymbol\beta\cdot{\bf n}\, g, v_h)_{\Gamma_-},
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="701" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img47.png"
+ ALT="$\displaystyle \sum_\kappa\left\{-(u_h,\boldsymbol\beta\cdot\nabla v_h)_\kappa+(...
+...mma}\right\}=(f,v_h)_\Omega-(\boldsymbol\beta\cdot{\bf n}\, g, v_h)_{\Gamma_-},$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(5)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
-where
-<!-- MATH: $\partial K_-:=\{x\in\partial K, \beta(x)\cdot n(x)<0\}$ -->
-<IMG
- WIDTH="263" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img28.gif"
- ALT="$\partial K_-:=\{x\in\partial K, \beta(x)\cdot n(x)<0\}$">
-denotes
-the inflow boundary and
-<!-- MATH: $\partial K_+=\partial K\setminus \partial
-K_-$ -->
-<IMG
- WIDTH="143" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img29.gif"
- ALT="$\partial K_+=\partial K\setminus \partial
-K_-$">
-the outflow part of cell <I>K</I>. Below, this equation will be
-referred to as <EM>first version</EM> of the DG method. We note that
-after integrating by parts again, we end up with: find
-<!-- MATH: $u_h\in V_h$ -->
-<IMG
- WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img21.gif"
- ALT="$u_h\in V_h$">such that
-<BR><P></P>
+where <!-- MATH
+ $\partial\kappa_-:=\{x\in\partial\kappa,
+\boldsymbol\beta(x)\cdot{\bf n}(x)<0\}$
+ -->
+<IMG
+ WIDTH="255" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img48.png"
+ ALT="$ \partial\kappa_-:=\{x\in\partial\kappa,
+\boldsymbol\beta(x)\cdot{\bf n}(x)<0\}$">
+ denotes the inflow boundary
+and <!-- MATH
+ $\partial\kappa_+=\partial\kappa\setminus \partial \kappa_-$
+ -->
+<IMG
+ WIDTH="127" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img49.png"
+ ALT="$ \partial\kappa_+=\partial\kappa\setminus \partial \kappa_-$">
+ the
+outflow part of cell <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img26.png"
+ ALT="$ \kappa$">
+. Below, this equation will be referred
+to as <I>first version</I> of the DG method. We note that after a
+second integration by parts, we obtain: find <!-- MATH
+ $u_h\in V_h$
+ -->
+<IMG
+ WIDTH="63" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img41.png"
+ ALT="$ u_h\in V_h$">
+ such that
+<P><!-- MATH
+ \begin{displaymath}
+\sum_\kappa\left\{(\nabla\cdot\{\boldsymbol\beta u_h\},\nabla v_h)_\kappa-(\boldsymbol\beta\cdot{\bf n} [u_h], v_h)_{\partial\kappa_-}\right\}=(f,v_h)_\Omega, \quad\forall v_h\in V_h,
+\end{displaymath}
+ -->
+</P>
<DIV ALIGN="CENTER">
-<!-- MATH: \begin{displaymath}
-\sum_K\left\{(\nabla\cdot\{\beta u_h\},\nabla v_h)_K-(\beta\cdot n [u_h], v_h)_{\partial K_-}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h,
-\end{displaymath} -->
-
-
<IMG
- WIDTH="509" HEIGHT="50"
- SRC="step-12.data/img30.gif"
- ALT="\begin{displaymath}\sum_K\left\{(\nabla\cdot\{\beta u_h\},\nabla v_h)_K-(\beta\c...
-..._{\partial K_-}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h,
-\end{displaymath}">
-</DIV>
-<BR CLEAR="ALL">
-<P></P>
-where
-<!-- MATH: $[u_h]=u_h-\hat u_h$ -->
-<IMG
- WIDTH="114" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img31.gif"
- ALT="$[u_h]=u_h-\hat u_h$">
-denotes the jump of the discrete function
-between two neighboring cells, that on the boundary of the domain is defined to be
-<!-- MATH: $[u_h]=u_h-g$ -->
-[<I>u</I><SUB><I>h</I></SUB>]=<I>u</I><SUB><I>h</I></SUB>-<I>g</I>. This is the discontinuous Galerkin
-scheme for the transport equation given in its original notation.
-Nevertheless, we will base the implementation of the scheme on the form
-given by (<A HREF="step-12.html#dg-scheme">3</A>) and (<A HREF="step-12.html#flux-transport-equation">4</A>), or
-(<A HREF="step-12.html#dg-transport">5</A>), respectively.
+ WIDTH="518" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img50.png"
+ ALT="$\displaystyle \sum_\kappa\left\{(\nabla\cdot\{\boldsymbol\beta u_h\},\nabla v_h...
+...u_h], v_h)_{\partial\kappa_-}\right\}=(f,v_h)_\Omega, \quad\forall v_h\in V_h,
+$">
+</DIV><P>
+</P>
+where <!-- MATH
+ $[u_h]=u_h^+-u_h^-$
+ -->
+<IMG
+ WIDTH="120" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img51.png"
+ ALT="$ [u_h]=u_h^+-u_h^-$">
+ denotes the jump of the discrete function
+between two neighboring cells and is defined to be <!-- MATH
+ $[u_h]=u_h^+-g$
+ -->
+<IMG
+ WIDTH="107" HEIGHT="39" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img52.png"
+ ALT="$ [u_h]=u_h^+-g$">
+ on
+the boundary of the domain. This is the discontinuous Galerkin scheme
+for the transport equation given in its original notation.
+Nevertheless, we will base the implementation of the scheme on the
+form given by (<A HREF="#dg-scheme">3</A>) and (<A HREF="#flux-transport-equation">4</A>),
+or (<A HREF="#dg-transport">5</A>), respectively.
<P>
-Finally we rewrite (<A HREF="step-12.html#dg-scheme">3</A>) in terms of a summation over all faces
-where each face
-<!-- MATH: $\Gamma=\partial
-K\cap\partial \hat K$ -->
-<IMG
- WIDTH="111" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img32.gif"
- ALT="$\Gamma=\partial
-K\cap\partial \hat K$">
-between two neighboring cells <I>K</I> and <IMG
- WIDTH="21" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img33.gif"
- ALT="$\hat K$">occurs twice: Find
-<!-- MATH: $u_h\in V_h$ -->
-<IMG
- WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img21.gif"
- ALT="$u_h\in V_h$">
-such that
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
--\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma\left\{(H(u_h,\hat u_h,n), v_h)_\Gamma+(H(\hat u_h, u_h,-n), \hat v_h)_{\Gamma\setminus\partial\Omega}\right\}=(f,v)_\Omega \quad\forall v_h\in V_h,
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-scheme-faces-long"> </A><IMG
- WIDTH="723" HEIGHT="50"
- SRC="step-12.data/img34.gif"
- ALT="\begin{displaymath}
--\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma\left\{(H(u_h,\hat ...
-...s\partial\Omega}\right\}=(f,v)_\Omega \quad\forall v_h\in V_h,
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+Finally, we rewrite (<A HREF="#dg-scheme">3</A>) in terms of a summation over all
+faces where each face <!-- MATH
+ $e=\partial \kappa\cap\partial \kappa'$
+ -->
+<IMG
+ WIDTH="100" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img53.png"
+ ALT="$ e=\partial \kappa\cap\partial \kappa'$">
+
+between two neighboring cells <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img26.png"
+ ALT="$ \kappa$">
+ and <IMG
+ WIDTH="19" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img54.png"
+ ALT="$ \kappa'$">
+ occurs twice:
+Find <!-- MATH
+ $u_h\in V_h$
+ -->
+<IMG
+ WIDTH="63" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img41.png"
+ ALT="$ u_h\in V_h$">
+ such that
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="dg-scheme-faces-long"></A><!-- MATH
+ \begin{equation}
+-\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e\left\{({\mathcal H}(u_h^+,u_h^-,{\bf n}), v_h)_e+({\mathcal H}(u_h^-, u_h^+,-{\bf n}), v_h^-)_{e\setminus\Gamma}\right\}=(f,v_h)_\Omega \quad\forall v_h\in V_h,
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="731" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img55.png"
+ ALT="$\displaystyle -\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e\left\{({...
+... n}), v_h^-)_{e\setminus\Gamma}\right\}=(f,v_h)_\Omega \quad\forall v_h\in V_h,$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(6)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
-By employing conservativity (<A HREF="step-12.html#conservative">2</A>) of the numerical flux
-this equation simplifies to: find
-<!-- MATH: $u_h\in V_h$ -->
-<IMG
- WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img21.gif"
- ALT="$u_h\in V_h$">
-such that
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
--\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n), [v_h])_{\Gamma\setminus\partial\Omega}+(H(u_h,g,n), v_h)_{\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-scheme-faces"> </A><IMG
- WIDTH="695" HEIGHT="50"
- SRC="step-12.data/img35.gif"
- ALT="\begin{displaymath}
--\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n),...
-...), v_h)_{\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+By employing conservativity (<A HREF="#conservative">2</A>) of the numerical flux
+this equation simplifies to: find <!-- MATH
+ $u_h\in V_h$
+ -->
+<IMG
+ WIDTH="63" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img41.png"
+ ALT="$ u_h\in V_h$">
+ such that
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="dg-scheme-faces"></A><!-- MATH
+ \begin{equation}
+-\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e({\mathcal H}(u_h^+,u_h^-,{\bf n}), [v_h])_{e\setminus\Gamma}+({\mathcal H}(u_h,g,{\bf n}), v_h)_{\Gamma}=(f,v_h)_\Omega \quad\forall v_h\in V_h.
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="687" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img56.png"
+ ALT="$\displaystyle -\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e({\mathca...
+...athcal H}(u_h,g,{\bf n}), v_h)_{\Gamma}=(f,v_h)_\Omega \quad\forall v_h\in V_h.$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(7)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
-Whereas the outer unit normal
-<!-- MATH: $n|_{\partial K}$ -->
-<IMG
- WIDTH="41" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img36.gif"
- ALT="$n\vert _{\partial K}$">
-is uniquely defined
-this is not so for <IMG
- WIDTH="24" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img37.gif"
- ALT="$n_\Gamma$">
-as the latter might be the normal from
-either side of the face. Hence, we need to fix the normal <I>n</I> on the
-face to be one of the two normals and denote the other normal by -<I>n</I>.
-This way we get -<I>n</I> in the second face term in
-(<A HREF="step-12.html#dg-scheme-faces-long">6</A>) that finally produces the minus sign in
-the jump [<I>v</I><SUB><I>h</I></SUB>] in equation (<A HREF="step-12.html#dg-scheme-faces">7</A>).
+Whereas the outer unit normal <!-- MATH
+ ${\bf n}|_{\partial\kappa}$
+ -->
+<IMG
+ WIDTH="37" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img57.png"
+ ALT="$ {\bf n}\vert _{\partial\kappa}$">
+ is uniquely
+defined this is not so for <IMG
+ WIDTH="23" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img58.png"
+ ALT="$ {\bf n}_e$">
+ as the latter might be the
+normal from either side of the face. Hence, we need to fix the normal
+<IMG
+ WIDTH="15" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img12.png"
+ ALT="$ {\bf n}$">
+ on the face to be one of the two normals and denote the
+other normal by <IMG
+ WIDTH="29" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img59.png"
+ ALT="$ -{\bf n}$">
+. This way we get <IMG
+ WIDTH="29" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img59.png"
+ ALT="$ -{\bf n}$">
+ in the second
+face term in (<A HREF="#dg-scheme-faces-long">6</A>) that finally produces the
+minus sign in the jump <IMG
+ WIDTH="31" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img60.png"
+ ALT="$ [v_h]$">
+ in equation (<A HREF="#dg-scheme-faces">7</A>).
<P>
-For the linear transport equation (<A HREF="step-12.html#transport-equation">1</A>)
-equation (<A HREF="step-12.html#dg-scheme-faces">7</A>) simplifies to
-<BR><P></P>
-<DIV ALIGN="CENTER">
-
-<!-- MATH: \begin{equation}
--\sum_K(u_h,\beta\cdot\nabla v_h)_K+\sum_\Gamma\left\{(\beta\cdot n\, u_h, [v_h])_{\Gamma_+\setminus\partial\Omega}+(\beta\cdot n\, \hat u_h, [v_h])_{\Gamma_-\setminus\partial\Omega}\right\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
-\end{equation} -->
-
-<TABLE WIDTH="100%" ALIGN="CENTER">
-<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-transport-gamma"> </A><IMG
- WIDTH="744" HEIGHT="50"
- SRC="step-12.data/img38.gif"
- ALT="\begin{displaymath}
--\sum_K(u_h,\beta\cdot\nabla v_h)_K+\sum_\Gamma\left\{(\bet...
-...ght\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
-\end{displaymath}"></TD>
-<TD WIDTH=10 ALIGN="RIGHT">
+For the linear transport equation (<A HREF="#transport-equation">1</A>)
+equation (<A HREF="#dg-scheme-faces">7</A>) simplifies to
+<P></P>
+<DIV ALIGN="CENTER"><A NAME="dg-transport-gamma"></A><!-- MATH
+ \begin{equation}
+-\sum_\kappa(u_h,\boldsymbol\beta\cdot\nabla v_h)_\kappa+\sum_e\left\{(\boldsymbol\beta\cdot{\bf n}\, u_h, [v_h])_{e_+\setminus\Gamma}+(\boldsymbol\beta\cdot{\bf n}\, u_h^-, [v_h])_{e_-\setminus\Gamma}\right\}=(f,v_h)_\Omega-(\boldsymbol\beta\cdot{\bf n}\, g, v_h)_{\Gamma_-},
+\end{equation}
+ -->
+<TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
+<TR VALIGN="MIDDLE">
+<TD NOWRAP ALIGN="CENTER"><IMG
+ WIDTH="730" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/intro/img61.png"
+ ALT="$\displaystyle -\sum_\kappa(u_h,\boldsymbol\beta\cdot\nabla v_h)_\kappa+\sum_e\l...
+...mma}\right\}=(f,v_h)_\Omega-(\boldsymbol\beta\cdot{\bf n}\, g, v_h)_{\Gamma_-},$"></TD>
+<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(8)</TD></TR>
-</TABLE>
-</DIV>
+</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
-which will be refered to as <EM>second version</EM> of the DG method.
+which will be refered to as <I>second version</I> of the DG method.
-<H3><A NAME="SECTION00004000000000000000">
+<H2><A NAME="SECTION00004000000000000000">
Implementation</A>
-</H3>
+</H2>
As already mentioned at the beginning of this example we will
implement assembling the system matrix in two different ways.
-The first one will be based on the first version (<A HREF="step-12.html#dg-transport">5</A>)
+The first one will be based on the first version (<A HREF="#dg-transport">5</A>)
of the DG method that includes a sum of integrals over all cell
-boundaries
-<!-- MATH: $\partial K$ -->
-<IMG
- WIDTH="31" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img39.gif"
- ALT="$\partial K$">.
-This is realized by a loop over all cells and
+boundaries <!-- MATH
+ $\partial\kappa$
+ -->
+<IMG
+ WIDTH="25" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img29.png"
+ ALT="$ \partial\kappa$">
+. This is realized by a loop over all cells and
a nested loop over all faces of each cell. Thereby each inner face
-
-<!-- MATH: $\Gamma=\partial K\cap\partial \hat K$ -->
-<IMG
- WIDTH="111" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img32.gif"
- ALT="$\Gamma=\partial
-K\cap\partial \hat K$">
-is treated twice, the first
-time when the outer loop treats cell <I>K</I> and the second time when it
+<!-- MATH
+ $e=\partial\kappa\cap\partial \kappa'$
+ -->
+<IMG
+ WIDTH="100" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img53.png"
+ ALT="$ e=\partial \kappa\cap\partial \kappa'$">
+ is treated twice, the first
+time when the outer loop treats cell <IMG
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img26.png"
+ ALT="$ \kappa$">
+ and the second time when it
treats cell <IMG
- WIDTH="21" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img33.gif"
- ALT="$\hat K$">.
-This way some values like the shape function
+ WIDTH="19" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img54.png"
+ ALT="$ \kappa'$">
+. This way some values like the shape function
values at quadrature points on faces need to be computed twice.
<P>
To overcome this overhead and for comparison, we implement
assembling of matrix also in a second and different way. This will
-be based on the second version (<A HREF="step-12.html#dg-transport-gamma">8</A>) that
+be based on the second version (<A HREF="#dg-transport-gamma">8</A>) that
includes a sum of integrals over all faces <IMG
- WIDTH="15" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img40.gif"
- ALT="$\Gamma$">.
-Here, several
+ WIDTH="12" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-12.data/intro/img62.png"
+ ALT="$ e$">
+. Here, several
difficulties occurs.
+
<UL>
<LI>As degrees of freedom are associated with cells (and not to faces)
and as a normal is only defined w.r.t. a cell adjacent to the face we
cannot simply run over all faces of the triangulation but need to
perform the nested loop over all cells and all faces of each cell
- like in the first implementation. This, because in <EM>deal.II</EM>
+ like in the first implementation. This, because in <TT>deal.II</TT>
faces are accessible from cells but not visa versa.
+</LI>
<LI>Due to the nested loop we arrive twice at each face. In order to
assemble face terms only once we either need to track which
faces we have treated before, or we introduce a simple rule that decides
treated from. Here, we employ the second approach and define the
following rule:
<DL COMPACT>
-<DT>a)
+<DT>a).</DT>
<DD>If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index on this level.
-<DT>b)
+
+</DD>
+<DT>b).</DT>
<DD>If the two cells are of different refinement levels we access
and treat the face from the coarser cell.
- </DL></UL>Before we start with the description of the code we first introduce
+
+</DD>
+</DL>
+</LI>
+</UL>
+Before we start with the description of the code we first introduce
its main ingredients. The main class is called
-<EM>DGMethod</EM>. It comprises all basic objects like the
+<TT>DGMethod</TT>. It comprises all basic objects like the
triangulation, the dofhandler, the system matrix and solution vectors.
Furthermore it has got some member functions, the most prominent of
-which are the <EM>assemble_system1</EM> and <EM>assemble_system2</EM>
+which are the <TT>assemble_system1</TT> and <TT>assemble_system2</TT>
functions that implement the two different ways mentioned above for
assembling the system matrix. Within these assembling routines several
different cases must be distinguished while performing the nested
loops over all cells and all faces of each cell and assembling the
respective face terms. While sitting on the current cell and looking
at a specific face there are the cases
-<DL COMPACT>
-<DT>1.
-<DD>face is at boundary,
-<DT>2.
-<DD>neighboring cell is finer,
-<DT>3.
-<DD>neighboring cell is of the same refinement level, and
-<DT>4.
-<DD>neighboring cell is coarser
-</DL>where the `neighboring cell' and the current cell have the mentioned
+
+<OL>
+<LI>face is at boundary,
+</LI>
+<LI>neighboring cell is finer,
+</LI>
+<LI>neighboring cell is of the same refinement level, and
+</LI>
+<LI>neighboring cell is coarser
+</LI>
+</OL>
+where the `neighboring cell' and the current cell have the mentioned
faces in common. In last three cases the assembling of the face terms
are almost the same. Hence, we can implement the assembling of the
face terms either by `copy and paste' (the lazy way, whose
the latter approach, of course. We go even further and encapsulate
this function and everything that is needed for assembling the
specific equation under consideration within a class called
-<EM>DGTransportEquation</EM>. This class includes objects of all
-equation-specific functions, the <EM>RHS</EM> and the
-<EM>BoundaryValues</EM> class, both derived from the <EM>Function</EM>
-class, and the <EM>Beta</EM> class representing the vector field.
-Furthermore, the <EM>DGTransportEquation</EM> class comprises member
-functions <EM>assemble_face_terms1</EM> and
-<EM>assemble_face_terms2</EM> that are invoked by the
-<EM>assemble_system1</EM> and <EM>assemble_system2</EM> functions of the
-<EM>DGMethod</EM>, respectively, and the functions
-<EM>assemble_cell_term</EM> and <EM>assemble_boundary_term</EM> that
+<TT>DGTransportEquation</TT>. This class includes objects of all
+equation-specific functions, the <TT>RHS</TT> and the
+<TT>BoundaryValues</TT> class, both derived from the <TT>Function</TT>
+class, and the <TT>Beta</TT> class representing the vector field.
+Furthermore, the <TT>DGTransportEquation</TT> class comprises member
+functions <TT>assemble_face_terms1</TT> and
+<TT>assemble_face_terms2</TT> that are invoked by the
+<TT>assemble_system1</TT> and <TT>assemble_system2</TT> functions of the
+<TT>DGMethod</TT>, respectively, and the functions
+<TT>assemble_cell_term</TT> and <TT>assemble_boundary_term</TT> that
are the same for both assembling routines. Due to the encapsulation of
all equation- and scheme-specific functions, the
-<EM>DGTransportEquation</EM> class can easily be replaced by a similar
+<TT>DGTransportEquation</TT> class can easily be replaced by a similar
class that implements a different equation and a different DG method.
-Indeed, the implementation of the <EM>assemble_system1</EM> and
-<EM>assemble_system2</EM> functions of the <EM>DGMethod</EM> class will
+Indeed, the implementation of the <TT>assemble_system1</TT> and
+<TT>assemble_system2</TT> functions of the <TT>DGMethod</TT> class will
be general enough to serve for different DG methods, different
equations, even for systems of equations (!) and, under small
-modifications, for nonlinear problems. Finally we note that the
+modifications, for nonlinear problems. Finally, we note that the
program is dimension independent, i.e. after replacing
-<EM>DGMethod<2></EM> by <EM>DGMethod<3></EM> the code runs in 3d.
+<TT>DGMethod<2></TT> by <TT>DGMethod<3></TT> the code runs in 3d.