From 388f455b6a69c8ef077d3d75f24130860ff9dc2c Mon Sep 17 00:00:00 2001 From: hartmann Date: Thu, 23 Feb 2006 17:14:18 +0000 Subject: [PATCH] Update notation. git-svn-id: https://svn.dealii.org/trunk@12478 0785d39b-7218-0410-832d-ea1e28bc413d --- .../step-12.data/intro.tex | 196 ++++++++-------- deal.II/examples/step-12/step-12.cc | 212 +++++++++--------- 2 files changed, 210 insertions(+), 198 deletions(-) diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-12.data/intro.tex b/deal.II/doc/tutorial/chapter-2.step-by-step/step-12.data/intro.tex index 9a4aaf67e1..a465e22431 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/step-12.data/intro.tex +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-12.data/intro.tex @@ -27,128 +27,136 @@ this example's purpose is to illustrate some implementational issues of the DG discretization only, here we simply consider the linear transport equation \begin{equation}\label{transport-equation} - \nabla\cdot \left\{\beta u\right\}=f \qquad\mbox{in }\Omega, + \nabla\cdot \left\{\boldsymbol\beta u\right\}=f \qquad\mbox{in }\Omega, \end{equation} subject to the boundary conditions \[ -u=g\quad\mbox{on }\partial\Omega_-. +u=g\quad\mbox{on }\Gamma_-, \] -Here, $\beta=\beta(x)$ is a vector field, $f$ a source function, $u$ the -(scalar) solution function, $g$ a boundary value function and +on the inflow part $\Gamma_-$ of the boundary $\Gamma=\partial\Omega$ +of the domain. Here, $\boldsymbol\beta=\boldsymbol\beta(x)$ denotes a +vector field, $f$ a source function, $u$ the (scalar) solution +function, $g$ a boundary value function, \[ -\partial\Omega_-:=\{x\in\partial\Omega, \beta(x)\cdot n(x)<0\} +\Gamma_-:=\{x\in\Gamma, \boldsymbol\beta(x)\cdot{\bf n}(x)<0\} \] -the inflow part of the boundary of the domain, whereas $n$ denotes the -unit outer normal to the boundary $\partial\Omega$. This -equation, \eqref{transport-equation}, is the conservative version of -the transport equation already considered in step 9 of this tutorial. +the inflow part of the boundary of the domain and ${\bf n}$ denotes +the unit outward normal to the boundary $\Gamma$. Equation +\eqref{transport-equation} is the conservative version of the +transport equation already considered in step 9 of this tutorial. In particular, we consider problem \eqref{transport-equation} on -$\Omega=[0,1]^2$ with $\beta=\frac{1}{|x|}(-x_2, x_1)$ representing a -circular counterclockwise flow field, $f=0$ and $g=1$ on -$x\in\partial\Omega_-^1:=[0,0.5]\times\{0\}$ and $g=0$ on $x\in -\partial\Omega_-\setminus \partial\Omega_-^1$. +$\Omega=[0,1]^2$ with $\boldsymbol\beta=\frac{1}{|x|}(-x_2, x_1)$ +representing a circular counterclockwise flow field, $f=0$ and $g=1$ +on $x\in\Gamma_-^1:=[0,0.5]\times\{0\}$ and $g=0$ on $x\in +\Gamma_-\setminus \Gamma_-^1$. \subsection{Discretization} -To follow the general paradigm of deriving DG discretizations for purely -hyperbolic equations we first consider the general hyperbolic problem +Following the general paradigm of deriving DG discretizations for +purely hyperbolic equations, we first consider the general hyperbolic +problem \[ - \nabla\cdot F(u)=f \qquad\mbox{in }\Omega, + \nabla\cdot {\mathcal F}(u)=f \qquad\mbox{in }\Omega, \] -subject to appropriate boundary conditions. Here $F$ denotes the flux -function of the equation under consideration that in our case, see -equation \eqref{transport-equation}, is represented by $F(u)=\beta u$. -For deriving the DG discretization we start with a variational, -mesh--dependent formulation of the problem, +subject to appropriate boundary conditions. Here ${\mathcal F}$ +denotes the flux function of the equation under consideration that in +our case, see equation \eqref{transport-equation}, is represented by +${\mathcal F}(u)=\boldsymbol\beta u$. For deriving the DG +discretization we start with a variational, mesh--dependent +formulation of the problem, \[ - \sum_K\left\{-(F(u),\nabla v)_K+(F(u)\cdot n, v)_{\partial K}\right\}=(f,v)_\Omega, + \sum_\kappa\left\{-({\mathcal F}(u),\nabla v)_\kappa+({\mathcal + F}(u)\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega, \] that originates from \eqref{transport-equation} by multiplication with -a test function and integration by parts on each cell $K$ of the -triangulation. Here $(\cdot, \cdot)_K$ and $(\cdot, \cdot)_{\partial - K}$ simply denote the integrals over the cell and the boundary of -the cell, respectively. To discretize the problem, the functions $u$ -and $v$ are replaced by discrete functions $u_h$ and $v_h$ that in the -case of discontinuous Galerkin methods belong to the space $V_h$ of +a test function $v$ and integration by parts on each cell $\kappa$ of +the triangulation. Here $(\cdot, \cdot)_\kappa$ and $(\cdot, +\cdot)_{\partial\kappa}$ simply denote the integrals over the cell +$\kappa$ and the boundary $\partial\kappa$ of the cell, +respectively. To discretize the problem, the functions $u$ and $v$ are +replaced by discrete functions $u_h$ and $v_h$ that in the case of +discontinuous Galerkin methods belong to the space $V_h$ of discontinuous piecewise polynomial functions of some degree $p$. Due to the discontinuity of the discrete function $u_h$ on interelement -faces, the flux $F(u)\cdot n$ must be replaced by a \emph{numerical - flux} function $H(u_h, \hat u_h, n)$, where $u_h|_{\partial K}$ -denotes the inner trace (w.r.t. the cell $K$) of $u_h$ and $\hat -u_h|_{\partial K}$ the outer trace, i.e. the value of $u_h$ on the -neighboring cell. Furthermore the numerical flux function $H$, among -other things, must be consistent, i.e. +faces, the flux ${\mathcal F}(u)\cdot{\bf n}$ must be replaced by a +\emph{numerical flux} function ${\mathcal H}(u_h^+, u_h^-, {\bf n})$, +where $u_h^+|_{\partial\kappa}$ denotes the inner trace (w.r.t. the +cell $\kappa$) of $u_h$ and $u_h^-|_{\partial\kappa}$ the outer trace, +i.e. the value of $u_h$ on the neighboring cell. Furthermore the +numerical flux function ${\mathcal H}$, among other things, must be +consistent, i.e. \[ -H(u,u,n)=F(u)\cdot n, +{\mathcal H}(u,u,{\bf n})={\mathcal F}(u)\cdot{\bf n}, \] and conservative, i.e. \begin{equation}\label{conservative} -H(v,w,n)=-H(w,v,-n). +{\mathcal H}(v,w,{\bf n})=-{\mathcal H}(w,v,-{\bf n}). \end{equation} This yields the following \emph{discontinuous Galerkin discretization}: find $u_h\in V_h$ such that \begin{equation}\label{dg-scheme} - \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. + \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} -Boundary conditions are realized by replacing $\hat u_h$ on the inflow boundary $\partial\Omega_-$ by the boundary function $g$. +Boundary conditions are realized by replacing $u_h^-$ on the inflow boundary $\Gamma_-$ by the boundary function $g$. In the special case of the transport equation \eqref{transport-equation} the numerical flux in its simplest form is given by \begin{equation}\label{flux-transport-equation} - H(u_h,\hat u_h,n)(x)=\left\{\begin{array}{ll} - (\beta\cdot n\, \hat u_h)(x),&\mbox{for } \beta(x)\cdot n(x)<0,\\ - (\beta\cdot n\, u_h)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0, + {\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} where on the inflow part of the cell the value is taken from the -neighboring cell, $\hat u_h$, and on the outflow part the value is -taken from the current cell, $u_h$. Hence, the discontinuous Galerkin +neighboring cell, $u_h^-$, and on the outflow part the value is +taken from the current cell, $u_h^+$. Hence, the discontinuous Galerkin scheme for the transport equation \eqref{transport-equation} is given by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following equation holds: \begin{equation}\label{dg-transport} - \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_-}, + \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} -where $\partial K_-:=\{x\in\partial K, \beta(x)\cdot n(x)<0\}$ denotes -the inflow boundary and $\partial K_+=\partial K\setminus \partial -K_-$ the outflow part of cell $K$. Below, this equation will be -referred to as \emph{first version} of the DG method. We note that -after integrating by parts again, we end up with: find $u_h\in V_h$ -such that +where $\partial\kappa_-:=\{x\in\partial\kappa, +\boldsymbol\beta(x)\cdot{\bf n}(x)<0\}$ denotes the inflow boundary +and $\partial\kappa_+=\partial\kappa\setminus \partial \kappa_-$ the +outflow part of cell $\kappa$. Below, this equation will be referred +to as \emph{first version} of the DG method. We note that after a +second integration by parts, we obtain: find $u_h\in V_h$ such that \[ - \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, + \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, \] -where $[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 $[u_h]=u_h-g$. 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 \eqref{dg-scheme} and \eqref{flux-transport-equation}, or -\eqref{dg-transport}, respectively. +where $[u_h]=u_h^+-u_h^-$ denotes the jump of the discrete function +between two neighboring cells and is defined to be $[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 \eqref{dg-scheme} and \eqref{flux-transport-equation}, +or \eqref{dg-transport}, respectively. -Finally we rewrite \eqref{dg-scheme} in terms of a summation over all faces -where each face $\Gamma=\partial -K\cap\partial \hat K$ between two neighboring cells $K$ and $\hat K$ -occurs twice: Find $u_h\in V_h$ such that +Finally, we rewrite \eqref{dg-scheme} in terms of a summation over all +faces where each face $e=\partial \kappa\cap\partial \kappa'$ +between two neighboring cells $\kappa$ and $\kappa'$ occurs twice: +Find $u_h\in V_h$ such that \begin{equation}\label{dg-scheme-faces-long} - -\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, + -\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} By employing conservativity \eqref{conservative} of the numerical flux this equation simplifies to: find $u_h\in V_h$ such that \begin{equation}\label{dg-scheme-faces} - -\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. + -\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} -Whereas the outer unit normal $n|_{\partial K}$ is uniquely defined -this is not so for $n_\Gamma$ as the latter might be the normal from -either side of the face. Hence, we need to fix the normal $n$ on the -face to be one of the two normals and denote the other normal by $-n$. -This way we get $-n$ in the second face term in -\eqref{dg-scheme-faces-long} that finally produces the minus sign in -the jump $[v_h]$ in equation \eqref{dg-scheme-faces}. +Whereas the outer unit normal ${\bf n}|_{\partial\kappa}$ is uniquely +defined this is not so for ${\bf n}_e$ as the latter might be the +normal from either side of the face. Hence, we need to fix the normal +${\bf n}$ on the face to be one of the two normals and denote the +other normal by $-{\bf n}$. This way we get $-{\bf n}$ in the second +face term in \eqref{dg-scheme-faces-long} that finally produces the +minus sign in the jump $[v_h]$ in equation \eqref{dg-scheme-faces}. For the linear transport equation \eqref{transport-equation} equation \eqref{dg-scheme-faces} simplifies to \begin{equation}\label{dg-transport-gamma} - -\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_-}, + -\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} which will be refered to as \emph{second version} of the DG method. \subsection{Implementation} @@ -156,24 +164,24 @@ 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 \eqref{dg-transport} of the DG method that includes a sum of integrals over all cell -boundaries $\partial K$. This is realized by a loop over all cells and +boundaries $\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 -$\Gamma=\partial K\cap\partial \hat K$ is treated twice, the first -time when the outer loop treats cell $K$ and the second time when it -treats cell $\hat K$. This way some values like the shape function +$e=\partial\kappa\cap\partial \kappa'$ is treated twice, the first +time when the outer loop treats cell $\kappa$ and the second time when it +treats cell $\kappa'$. This way some values like the shape function values at quadrature points on faces need to be computed twice. 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 \eqref{dg-transport-gamma} that -includes a sum of integrals over all faces $\Gamma$. Here, several +includes a sum of integrals over all faces $e$. Here, several difficulties occurs. \begin{itemize} \item 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 \emph{deal.II} + like in the first implementation. This, because in {\tt deal.II} faces are accessible from cells but not visa versa. \item 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 @@ -190,10 +198,10 @@ difficulties occurs. \end{itemize} Before we start with the description of the code we first introduce its main ingredients. The main class is called -\emph{DGMethod}. It comprises all basic objects like the +{\tt DGMethod}. 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 \emph{assemble\_system1} and \emph{assemble\_system2} +which are the {\tt assemble\_system1} and {\tt assemble\_system2} 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 @@ -216,27 +224,27 @@ three cases. To be kind of educational within this tutorial we perform 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 -\emph{DGTransportEquation}. This class includes objects of all -equation--specific functions, the \emph{RHS} and the -\emph{BoundaryValues} class, both derived from the \emph{Function} -class, and the \emph{Beta} class representing the vector field. -Furthermore, the \emph{DGTransportEquation} class comprises member -functions \emph{assemble\_face\_terms1} and -\emph{assemble\_face\_terms2} that are invoked by the -\emph{assemble\_system1} and \emph{assemble\_system2} functions of the -\emph{DGMethod}, respectively, and the functions -\emph{assemble\_cell\_term} and \emph{assemble\_boundary\_term} that +{\tt DGTransportEquation}. This class includes objects of all +equation--specific functions, the {\tt RHS} and the +{\tt BoundaryValues} class, both derived from the {\tt Function} +class, and the {\tt Beta} class representing the vector field. +Furthermore, the {\tt DGTransportEquation} class comprises member +functions {\tt assemble\_face\_terms1} and +{\tt assemble\_face\_terms2} that are invoked by the +{\tt assemble\_system1} and {\tt assemble\_system2} functions of the +{\tt DGMethod}, respectively, and the functions +{\tt assemble\_cell\_term} and {\tt assemble\_boundary\_term} that are the same for both assembling routines. Due to the encapsulation of all equation- and scheme-specific functions, the -\emph{DGTransportEquation} class can easily be replaced by a similar +{\tt DGTransportEquation} class can easily be replaced by a similar class that implements a different equation and a different DG method. -Indeed, the implementation of the \emph{assemble\_system1} and -\emph{assemble\_system2} functions of the \emph{DGMethod} class will +Indeed, the implementation of the {\tt assemble\_system1} and +{\tt assemble\_system2} functions of the {\tt DGMethod} 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 -\emph{DGMethod<2>} by \emph{DGMethod<3>} the code runs in 3d. +{\tt DGMethod<2>} by {\tt DGMethod<3>} the code runs in 3d. \end{document} diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc index 8ba4184f18..c577f33f78 100644 --- a/deal.II/examples/step-12/step-12.cc +++ b/deal.II/examples/step-12/step-12.cc @@ -42,20 +42,21 @@ // due to the pre-knowledge of the // actual polynomial degree 1. #include - // Here the discontinuous finite elements are - // defined. They are used in the same way as - // all other finite elements, though -- as - // you have seen in previous tutorial + // Here the discontinuous finite + // elements are defined. They are + // used in the same way as all other + // finite elements, though -- as you + // have seen in previous tutorial // programs -- there isn't much user - // interaction with finite element classes at - // all: the are passed to ``DoFHandler'' and - // ``FEValues'' objects, and that is about - // it. + // interaction with finite element + // classes at all: the are passed to + // ``DoFHandler'' and ``FEValues'' + // objects, and that is about it. #include // We are going to use the simplest // possible solver, called Richardson - // iteration, that represents a simple - // defect correction. This, in + // iteration, that represents a + // simple defect correction. This, in // combination with a block SSOR // preconditioner (defined in // precondition_block.h), that uses @@ -231,24 +232,24 @@ class DGTransportEquation DGTransportEquation(); void assemble_cell_term(const FEValues& fe_v, - FullMatrix &u_v_matrix, + FullMatrix &ui_vi_matrix, Vector &cell_vector) const; void assemble_boundary_term(const FEFaceValues& fe_v, - FullMatrix &u_v_matrix, + FullMatrix &ui_vi_matrix, Vector &cell_vector) const; void assemble_face_term1(const FEFaceValuesBase& fe_v, const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix) const; + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const; void assemble_face_term2(const FEFaceValuesBase& fe_v, const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix, - FullMatrix &u_vn_matrix, - FullMatrix &un_vn_matrix) const; + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_matrix) const; private: const Beta beta_function; const RHS rhs_function; @@ -270,7 +271,7 @@ DGTransportEquation::DGTransportEquation () // The ``assemble_cell_term'' // function assembles the cell terms // of the discretization. - // ``u_v_matrix'' is a cell matrix, + // ``ui_vi_matrix'' is a cell matrix, // i.e. for a DG method of degree 1, // it is of size 4 times 4, and // ``cell_vector'' is of size 4. @@ -281,7 +282,7 @@ DGTransportEquation::DGTransportEquation () template void DGTransportEquation::assemble_cell_term( const FEValues &fe_v, - FullMatrix &u_v_matrix, + FullMatrix &ui_vi_matrix, Vector &cell_vector) const { // First we ask ``fe_v'' for the @@ -305,7 +306,7 @@ void DGTransportEquation::assemble_cell_term( for (unsigned int i=0; i::assemble_cell_term( template void DGTransportEquation::assemble_boundary_term( const FEFaceValues& fe_v, - FullMatrix &u_v_matrix, + FullMatrix &ui_vi_matrix, Vector &cell_vector) const { // Again, as in the previous @@ -356,18 +357,18 @@ void DGTransportEquation::assemble_boundary_term( const double beta_n=beta[point] * normals[point]; // We assemble the term // $(\beta\cdot n - // u,v)_{\partial K_+}$, + // u,v)_{\partial\kappa_+}$, if (beta_n>0) for (unsigned int i=0; i::assemble_boundary_term( // on the face. // // In addition to the cell matrix - // ``u_v_matrix'' this function gets - // a new argument ``un_v_matrix'', + // ``ui_vi_matrix'' this function gets + // a new argument ``ue_vi_matrix'', // that stores contributions to the // system matrix that are based on // outer values of u, see $\hat u_h$ @@ -411,8 +412,8 @@ template void DGTransportEquation::assemble_face_term1( const FEFaceValuesBase& fe_v, const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix) const + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const { // Again, as in the previous // function, we ask the FEValues @@ -435,11 +436,11 @@ void DGTransportEquation::assemble_face_term1( const double beta_n=beta[point] * normals[point]; // We assemble the term // $(\beta\cdot n - // u,v)_{\partial K_+}$, + // u,v)_{\partial\kappa_+}$, if (beta_n>0) for (unsigned int i=0; i::assemble_face_term1( // and the // term $(\beta\cdot n // \hat u,v)_{\partial - // K_-}$. + // \kappa_-}$. for (unsigned int i=0; i::assemble_face_term1( // terms are given as a sum of // integrals over all faces. Here we // need two additional cell matrices - // ``u_vn_matrix'' and - // ``un_vn_matrix'' that will store + // ``ui_ve_matrix'' and + // ``ue_ve_matrix'' that will store // contributions due to terms // involving u and vn as well as un // and vn. @@ -478,10 +479,10 @@ template void DGTransportEquation::assemble_face_term2( const FEFaceValuesBase& fe_v, const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix, - FullMatrix &u_vn_matrix, - FullMatrix &un_vn_matrix) const + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_matrix) const { // the first few lines are the same const std::vector &JxW = fe_v.get_JxW_values (); @@ -496,10 +497,11 @@ void DGTransportEquation::assemble_face_term2( const double beta_n=beta[point] * normals[point]; if (beta_n>0) { - // This terms we've already seen. + // This term we've already + // seen. for (unsigned int i=0; i::assemble_face_term2( // We additionally assemble // the term $(\beta\cdot n // u,\hat v)_{\partial - // K_+}, + // \kappa_+}, for (unsigned int k=0; k::assemble_face_term2( // seen, too. for (unsigned int i=0; i::assemble_face_term2( // And this is another new // one: $(\beta\cdot n \hat // u,\hat v)_{\partial - // K_-}$. + // \kappa_-}$. for (unsigned int k=0; k::assemble_system1 () // and vectors. Here we need two // cell matrices, both for face // terms that include test - // functions ``v'' (shape functions + // functions ``vi'' (internal shape + // functions, i.e. shape functions // of the current cell). To be more // precise, the first matrix will - // include the `u and v terms' and - // the second that will include the - // `un and v terms'. Here we recall - // the convention that `un' is - // the shortcut for `u_neighbor' - // and represents the $u_hat$, see + // include the `ui and vi terms' + // and the second will include the + // `ue and vi terms'. Here we + // recall the convention that `ui' + // is the shortcut for $u_h^+$ and + // `ue' represents $u_h_-$, see the // introduction. - FullMatrix u_v_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix un_v_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); Vector cell_vector (dofs_per_cell); @@ -835,10 +838,10 @@ void DGMethod::assemble_system1 () // cell vector are only // ADDED. Therefore on each // cell we need to reset the - // ``u_v_matrix'' and + // ``ui_vi_matrix'' and // ``cell_vector'' to zero, // before assembling the cell terms. - u_v_matrix = 0; + ui_vi_matrix = 0; cell_vector = 0; // Now we reinit the ``FEValues'' @@ -852,7 +855,7 @@ void DGMethod::assemble_system1 () // previously reinit'ed on the // current cell. dg.assemble_cell_term(fe_v, - u_v_matrix, + ui_vi_matrix, cell_vector); // As in previous examples the @@ -869,9 +872,9 @@ void DGMethod::assemble_system1 () typename DoFHandler::face_iterator face=cell->face(face_no); // and clear the - // ``un_v_matrix'' on each + // ``ue_vi_matrix'' on each // face. - un_v_matrix = 0; + ue_vi_matrix = 0; // Now we distinguish the // four different cases in @@ -891,7 +894,7 @@ void DGMethod::assemble_system1 () // corresponding face // terms. dg.assemble_boundary_term(fe_v_face, - u_v_matrix, + ui_vi_matrix, cell_vector); } else @@ -984,7 +987,7 @@ void DGMethod::assemble_system1 () // We need to // reset the - // ``un_v_matrix'' + // ``ue_vi_matrix'' // on each // subface // because on @@ -994,7 +997,7 @@ void DGMethod::assemble_system1 () // different // neighboring // cells. - un_v_matrix = 0; + ue_vi_matrix = 0; // As already // mentioned @@ -1022,8 +1025,8 @@ void DGMethod::assemble_system1 () dg.assemble_face_term1(fe_v_subface, fe_v_face_neighbor, - u_v_matrix, - un_v_matrix); + ui_vi_matrix, + ue_vi_matrix); // Then we get // the dof @@ -1035,13 +1038,13 @@ void DGMethod::assemble_system1 () // and // distribute - // ``un_v_matrix'' + // ``ue_vi_matrix'' // to the // system_matrix for (unsigned int i=0; ihas_children())'' @@ -1086,8 +1089,8 @@ void DGMethod::assemble_system1 () dg.assemble_face_term1(fe_v_face, fe_v_face_neighbor, - u_v_matrix, - un_v_matrix); + ui_vi_matrix, + ue_vi_matrix); // End of ``if // (neighbor->level() // == @@ -1142,8 +1145,8 @@ void DGMethod::assemble_system1 () dg.assemble_face_term1(fe_v_face, fe_v_subface_neighbor, - u_v_matrix, - un_v_matrix); + ui_vi_matrix, + ue_vi_matrix); } // Now we get the @@ -1154,11 +1157,11 @@ void DGMethod::assemble_system1 () neighbor->get_dof_indices (dofs_neighbor); // and distribute the - // ``un_v_matrix''. + // ``ue_vi_matrix''. for (unsigned int i=0; i::assemble_system1 () } // Finally we distribute the - // ``u_v_matrix'' + // ``ui_vi_matrix'' for (unsigned int i=0; i::assemble_system2 () mapping, fe, face_quadrature, neighbor_face_update_flags); - FullMatrix u_v_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix un_v_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); // Additionally we need the // following two cell matrices, // both for face term that include - // test function ``vn'' (shape + // test function ``ve'' (external + // shape functions, i.e. shape // functions of the neighboring // cell). To be more precise, the // first matrix will include the `u // and vn terms' and the second // that will include the `un and vn // terms'. - FullMatrix u_vn_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix un_vn_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ui_ve_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ue_ve_matrix (dofs_per_cell, dofs_per_cell); Vector cell_vector (dofs_per_cell); @@ -1272,13 +1276,13 @@ void DGMethod::assemble_system2 () endc = dof_handler.end(); for (;cell!=endc; ++cell) { - u_v_matrix = 0; + ui_vi_matrix = 0; cell_vector = 0; fe_v.reinit (cell); dg.assemble_cell_term(fe_v, - u_v_matrix, + ui_vi_matrix, cell_vector); cell->get_dof_indices (dofs); @@ -1294,7 +1298,7 @@ void DGMethod::assemble_system2 () fe_v_face.reinit (cell, face_no); dg.assemble_boundary_term(fe_v_face, - u_v_matrix, + ui_vi_matrix, cell_vector); } else @@ -1319,19 +1323,19 @@ void DGMethod::assemble_system2 () ExcInternalError()); Assert (!neighbor_child->has_children(), ExcInternalError()); - un_v_matrix = 0; - u_vn_matrix = 0; - un_vn_matrix = 0; + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; fe_v_subface.reinit (cell, face_no, subface_no); fe_v_face_neighbor.reinit (neighbor_child, neighbor2); dg.assemble_face_term2(fe_v_subface, fe_v_face_neighbor, - u_v_matrix, - un_v_matrix, - u_vn_matrix, - un_vn_matrix); + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); neighbor_child->get_dof_indices (dofs_neighbor); @@ -1339,11 +1343,11 @@ void DGMethod::assemble_system2 () for (unsigned int j=0; j::assemble_system2 () { const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); - un_v_matrix = 0; - u_vn_matrix = 0; - un_vn_matrix = 0; + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; fe_v_face.reinit (cell, face_no); fe_v_face_neighbor.reinit (neighbor, neighbor2); dg.assemble_face_term2(fe_v_face, fe_v_face_neighbor, - u_v_matrix, - un_v_matrix, - u_vn_matrix, - un_vn_matrix); + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); neighbor->get_dof_indices (dofs_neighbor); @@ -1377,11 +1381,11 @@ void DGMethod::assemble_system2 () for (unsigned int j=0; j::assemble_system2 () for (unsigned int i=0; i