From 2a0e91226953a4eaebb86a4e1344689175c0f786 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 3 Jul 2023 20:47:23 -0600 Subject: [PATCH] Use balanced braces for doxygen formulas. --- examples/step-11/doc/intro.dox | 2 +- examples/step-14/doc/intro.dox | 8 +-- examples/step-15/doc/intro.dox | 4 +- examples/step-15/doc/results.dox | 4 +- examples/step-18/doc/intro.dox | 14 ++--- examples/step-20/doc/intro.dox | 30 +++++----- examples/step-21/doc/intro.dox | 30 +++++----- examples/step-22/doc/intro.dox | 58 +++++++++---------- examples/step-22/doc/results.dox | 6 +- examples/step-23/doc/intro.dox | 8 +-- examples/step-24/doc/intro.dox | 12 ++-- examples/step-25/doc/intro.dox | 2 +- examples/step-26/doc/intro.dox | 36 ++++++------ examples/step-26/doc/results.dox | 10 ++-- examples/step-28/doc/intro.dox | 56 +++++++++--------- examples/step-29/doc/intro.dox | 8 +-- examples/step-3/doc/intro.dox | 20 +++---- examples/step-31/doc/intro.dox | 46 +++++++-------- examples/step-31/doc/results.dox | 4 +- examples/step-32/doc/intro.dox | 44 +++++++------- examples/step-32/doc/results.dox | 4 +- examples/step-32/step-32.cc | 2 +- examples/step-33/doc/intro.dox | 16 ++--- examples/step-33/step-33.cc | 2 +- examples/step-34/doc/intro.dox | 16 ++--- examples/step-35/doc/intro.dox | 4 +- examples/step-36/doc/intro.dox | 6 +- examples/step-37/doc/intro.dox | 12 ++-- examples/step-37/doc/results.dox | 4 +- examples/step-38/doc/intro.dox | 6 +- examples/step-39/doc/intro.dox | 4 +- examples/step-40/doc/intro.dox | 4 +- examples/step-41/doc/intro.dox | 50 ++++++++-------- examples/step-42/doc/intro.dox | 26 ++++----- examples/step-43/doc/intro.dox | 32 +++++----- examples/step-44/doc/intro.dox | 20 +++---- examples/step-44/step-44.cc | 10 ++-- examples/step-45/doc/intro.dox | 4 +- examples/step-45/step-45.cc | 4 +- examples/step-46/doc/intro.dox | 28 ++++----- examples/step-5/doc/intro.dox | 4 +- examples/step-50/doc/intro.dox | 14 ++--- examples/step-51/doc/intro.dox | 20 +++---- examples/step-52/doc/intro.dox | 32 +++++----- examples/step-55/doc/intro.dox | 6 +- examples/step-56/doc/intro.dox | 14 ++--- examples/step-57/doc/intro.dox | 34 +++++------ examples/step-59/doc/intro.dox | 2 +- examples/step-61/doc/intro.dox | 24 ++++---- examples/step-62/doc/intro.dox | 6 +- examples/step-62/step-62.cc | 2 +- examples/step-64/step-64.cc | 2 +- examples/step-66/doc/intro.dox | 2 +- examples/step-68/doc/intro.dox | 2 +- examples/step-69/doc/intro.dox | 30 +++++----- examples/step-69/step-69.cc | 8 +-- examples/step-7/doc/intro.dox | 4 +- examples/step-70/doc/intro.dox | 8 +-- examples/step-71/doc/intro.dox | 4 +- examples/step-71/step-71.cc | 10 ++-- examples/step-72/doc/intro.dox | 6 +- examples/step-74/doc/intro.dox | 2 +- examples/step-79/doc/intro.dox | 2 +- examples/step-8/doc/intro.dox | 8 +-- examples/step-81/doc/intro.dox | 14 ++--- examples/step-85/doc/intro.dox | 42 +++++++------- examples/step-85/step-85.cc | 4 +- examples/step-9/doc/intro.dox | 4 +- include/deal.II/lac/scalapack.templates.h | 2 +- .../deal.II/numerics/smoothness_estimator.h | 8 +-- 70 files changed, 488 insertions(+), 488 deletions(-) diff --git a/examples/step-11/doc/intro.dox b/examples/step-11/doc/intro.dox index ab72e1fcaa..d092c787d2 100644 --- a/examples/step-11/doc/intro.dox +++ b/examples/step-11/doc/intro.dox @@ -3,7 +3,7 @@ The problem we will be considering is the solution of Laplace's problem with Neumann boundary conditions only: -@f{eqnarray*} +@f{eqnarray*}{ -\Delta u &=& f \qquad \mathrm{in}\ \Omega, \\ \partial_n u &=& g \qquad \mathrm{on}\ \partial\Omega. diff --git a/examples/step-14/doc/intro.dox b/examples/step-14/doc/intro.dox index 308d8d4db2..93405873ae 100644 --- a/examples/step-14/doc/intro.dox +++ b/examples/step-14/doc/intro.dox @@ -79,7 +79,7 @@ sum over cells where each cell's contribution can then be used as an error indicator for this cell. Thus, we split the scalar products into terms for each cell, and integrate by parts on each of them: -@f{eqnarray*} +@f{eqnarray*}{ J(e) &=& \sum_K (\nabla (u-u_h), \nabla (z-\varphi_h))_K @@ -98,7 +98,7 @@ with which this term could cancel, the weight $z-\varphi_h$ can be chosen as zero, and the whole term disappears. Thus, we have -@f{eqnarray*} +@f{eqnarray*}{ J(e) &=& \sum_K (f+\Delta u_h, z-\varphi_h)_K @@ -109,7 +109,7 @@ the value of this quantity as taken from this side of the cell (for the usual Lagrange elements, derivatives are not continuous across edges). We then rewrite the above formula by exchanging half of the edge integral of cell $K$ with the neighbor cell $K'$, to obtain -@f{eqnarray*} +@f{eqnarray*}{ J(e) &=& \sum_K (f+\Delta u_h, z-\varphi_h)_K @@ -126,7 +126,7 @@ the jump of the normal derivative by and get the final form after setting the discrete function $\varphi_h$, which is by now still arbitrary, to the point interpolation of the dual solution, $\varphi_h=I_h z$: -@f{eqnarray*} +@f{eqnarray*}{ J(e) &=& \sum_K (f+\Delta u_h, z-I_h z)_K diff --git a/examples/step-15/doc/intro.dox b/examples/step-15/doc/intro.dox index 109bd83f83..8347c9cb94 100644 --- a/examples/step-15/doc/intro.dox +++ b/examples/step-15/doc/intro.dox @@ -40,7 +40,7 @@ have to use Newton's method to compute the solution iteratively. In a classical sense, the problem is given in the following form: - @f{align*} + @f{align*}{ -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right) &= 0 \qquad \qquad &&\textrm{in} ~ \Omega \\ @@ -53,7 +53,7 @@ space. In this example, we choose $\Omega$ as the unit disk. As described above, we solve this equation using Newton's method in which we compute the $n$th approximate solution from the $(n-1)$th one, and use a damping parameter $\alpha^n$ to get better global convergence behavior: - @f{align*} + @f{align*}{ F'(u^{n},\delta u^{n})&=- F(u^{n}) \\ u^{n+1}&=u^{n}+\alpha^n \delta u^{n} diff --git a/examples/step-15/doc/results.dox b/examples/step-15/doc/results.dox index 858f39b3fa..7ac1adb84f 100644 --- a/examples/step-15/doc/results.dox +++ b/examples/step-15/doc/results.dox @@ -221,11 +221,11 @@ the form @f] we use a Newton iteration that requires us to repeatedly solve the linear partial differential equation - @f{align*} + @f{align*}{ F'(u^{n},\delta u^{n}) &=- F(u^{n}) @f} so that we can compute the update - @f{align*} + @f{align*}{ u^{n+1}&=u^{n}+\alpha^n \delta u^{n} @f} with the solution $\delta u^{n}$ of the Newton step. For the problem diff --git a/examples/step-18/doc/intro.dox b/examples/step-18/doc/intro.dox index 91df0426b4..3a84a62c60 100644 --- a/examples/step-18/doc/intro.dox +++ b/examples/step-18/doc/intro.dox @@ -48,7 +48,7 @@ In addition, initial conditions @f] and Dirichlet (displacement) or Neumann (traction) boundary conditions need to be specified for a unique solution: -@f{eqnarray*} +@f{eqnarray*}{ \mathbf{u}(\mathbf{x},t) &=& \mathbf{d}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_D\subset\partial\Omega, @@ -82,7 +82,7 @@ that all changes in external forcing happen on times scales that are much larger than $\tau$. In that case, the dynamic nature of the change is unimportant: we can consider the body to always be in static equilibrium, i.e. we can assume that at all times the body satisfies -@f{eqnarray*} +@f{eqnarray*}{ - \textrm{div}\ ( C \varepsilon(\mathbf{u})) &=& \mathbf{f}(\mathbf{x},t) \qquad \textrm{in}\ \Omega, @@ -127,7 +127,7 @@ implementing a real model, as we will see in step-44. To come back to defining our "artificial" model, let us first introduce a tensorial stress variable $\sigma$, and write the differential equations in terms of the stress: -@f{eqnarray*} +@f{eqnarray*}{ - \textrm{div}\ \sigma &=& \mathbf{f}(\mathbf{x},t) \qquad \textrm{in}\ \Omega(t), @@ -172,7 +172,7 @@ and $\Delta \mathbf{u}^n$ the incremental displacement for time step $n$. In addition, we have to specify initial data $\mathbf{u}(\cdot,0)=\mathbf{u}_0$. This way, if we want to solve for the displacement increment, we have to solve the following system: -@f{align*} +@f{align*}{ - \textrm{div}\ C \varepsilon(\Delta\mathbf{u}^n) &= \mathbf{f} + \textrm{div}\ \sigma^{n-1} \qquad &&\textrm{in}\ \Omega(t_{n-1}), @@ -190,7 +190,7 @@ finite element formulation, reads as follows: find $\Delta \mathbf{u}^n \in \{v\in H^1(\Omega(t_{n-1}))^d: v|_{\Gamma_D}=\mathbf{d}(\cdot,t_n) - \mathbf{d}(\cdot,t_{n-1})\}$ such that -@f{align*} +@f{align*}{ (C \varepsilon(\Delta\mathbf{u}^n), \varepsilon(\varphi) )_{\Omega(t_{n-1})} &= (\mathbf{f}, \varphi)_{\Omega(t_{n-1})} @@ -207,7 +207,7 @@ Using that $\sigma^{n-1} \mathbf{n} = [C \varepsilon(\mathbf{u}^{n-1})] \mathbf{n} = \mathbf{b}(\mathbf x, t_{n-1})$, these equations can be simplified to -@f{align*} +@f{align*}{ (C \varepsilon(\Delta\mathbf{u}^n), \varepsilon(\varphi) )_{\Omega(t_{n-1})} &= (\mathbf{f}, \varphi)_{\Omega(t_{n-1})} @@ -529,7 +529,7 @@ for (unsigned int i=0; imixed formulation: -@f{eqnarray*} +@f{eqnarray*}{ K^{-1} {\mathbf u} + \nabla p &=& 0 \qquad {\textrm{in}\ } \Omega, \\ -{\textrm{div}}\ {\mathbf u} &=& -f \qquad {\textrm{in}\ }\Omega, \\ p &=& g \qquad {\textrm{on}\ } \partial\Omega. @@ -103,11 +103,11 @@ under the common assumption that $K$ is a symmetric tensor. The weak formulation of this problem is found by multiplying the two equations with test functions and integrating some terms by parts: -@f{eqnarray*} +@f{eqnarray*}{ A(\{{\mathbf u},p\},\{{\mathbf v},q\}) = F(\{{\mathbf v},q\}), @f} where -@f{eqnarray*} +@f{eqnarray*}{ A(\{{\mathbf u},p\},\{{\mathbf v},q\}) &=& ({\mathbf v}, K^{-1}{\mathbf u})_\Omega - ({\textrm{div}}\ {\mathbf v}, p)_\Omega @@ -141,7 +141,7 @@ elements, for example the classic book by Brenner and Scott, also state the relevant results. In any case, with appropriate choices of function spaces, the discrete formulation reads as follows: Find ${\mathbf u}_h,p_h$ so that -@f{eqnarray*} +@f{eqnarray*}{ A(\{{\mathbf u}_h,p_h\},\{{\mathbf v}_h,q_h\}) = F(\{{\mathbf v}_h,q_h\}) \qquad\qquad \forall {\mathbf v}_h,q_h. @f} @@ -156,20 +156,20 @@ zero everywhere else. If we also choose ${\mathbf v}=0$ everywhere (remember that the weak form above has to hold for all discrete test functions $q,v$), then putting these choices of test functions into the weak formulation above implies in particular that -@f{eqnarray*} +@f{eqnarray*}{ - (1,{\textrm{div}}\ {\mathbf u}_h)_K = -(1,f)_K, @f} which we can of course write in more explicit form as -@f{eqnarray*} +@f{eqnarray*}{ \int_K {\textrm{div}}\ {\mathbf u}_h = \int_K f. @f} Applying the divergence theorem results in the fact that ${\mathbf u}_h$ has to satisfy, for every choice of cell $K$, the relationship -@f{eqnarray*} +@f{eqnarray*}{ \int_{\partial K} {\mathbf u}_h\cdot{\mathbf n} = \int_K f. @@ -202,7 +202,7 @@ The deal.II library (of course) implements Raviart-Thomas elements $RT(k)$ of arbitrary order $k$, as well as discontinuous elements $DG(k)$. If we forget about their particular properties for a second, we then have to solve a discrete problem -@f{eqnarray*} +@f{eqnarray*}{ A(x_h,w_h) = F(w_h), @f} with the bilinear form and right hand side as stated above, and $x_h=\{{\mathbf u}_h,p_h\}$, $w_h=\{{\mathbf v}_h,q_h\}$. Both $x_h$ and $w_h$ are from the space @@ -263,7 +263,7 @@ return a non-zero value for more than just one component. We could now attempt to rewrite the bilinear form above in terms of vector components. For example, in 2d, the first term could be rewritten like this (note that $u_0=x_0, u_1=x_1, p=x_2$): -@f{eqnarray*} +@f{eqnarray*}{ ({\mathbf u}_h^i, K^{-1}{\mathbf u}_h^j) = &\left((x_h^i)_0, K^{-1}_{00} (x_h^j)_0\right) + @@ -463,7 +463,7 @@ In view of the difficulties using standard solvers and preconditioners mentioned above, let us take another look at the matrix. If we sort our degrees of freedom so that all velocity come before all pressure variables, then we can subdivide the linear system $Ax=b$ into the following blocks: -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} M & B \\ B^T & 0 \end{array}\right) @@ -483,7 +483,7 @@ to the gradient. By block elimination, we can then re-order this system in the following way (multiply the first row of the system by $B^TM^{-1}$ and then subtract the second row from it): -@f{eqnarray*} +@f{eqnarray*}{ B^TM^{-1}B P &=& B^TM^{-1} F - G, \\ MU &=& F - BP. @f} @@ -713,7 +713,7 @@ We will try something along the second approach, as much to improve the performance of the program as to demonstrate some techniques. To this end, let us recall that the ideal preconditioner is, of course, $S^{-1}$, but that is unattainable. However, how about -@f{eqnarray*} +@f{eqnarray*}{ \tilde S^{-1} = [B^T ({\textrm{diag}\ }M)^{-1}B]^{-1} @f} as a preconditioner? That would mean that every time we have to do one @@ -782,12 +782,12 @@ formulation as stated above. Since we want to monitor convergence of the solution inside the program, we choose right hand side, boundary conditions, and the coefficient so that we recover a solution function known to us. In particular, we choose the pressure solution -@f{eqnarray*} +@f{eqnarray*}{ p = -\left(\frac \alpha 2 xy^2 + \beta x - \frac \alpha 6 x^3\right), @f} and for the coefficient we choose the unit matrix $K_{ij}=\delta_{ij}$ for simplicity. Consequently, the exact velocity satisfies -@f{eqnarray*} +@f{eqnarray*}{ {\mathbf u} = \left(\begin{array}{cc} \frac \alpha 2 y^2 + \beta - \frac \alpha 2 x^2 \\ diff --git a/examples/step-21/doc/intro.dox b/examples/step-21/doc/intro.dox index 96877da161..3068d8c028 100644 --- a/examples/step-21/doc/intro.dox +++ b/examples/step-21/doc/intro.dox @@ -38,7 +38,7 @@ however. The velocity with which molecules of each of the two phases move is determined by Darcy's law that states that the velocity is proportional to the pressure gradient: -@f{eqnarray*} +@f{eqnarray*}{ \mathbf{u}_{j} = -\frac{k_{rj}(S)}{\mu_{j}} \mathbf{K} \cdot \nabla p @@ -61,7 +61,7 @@ each phase, with a source term for each phase. By summing over the two phases, we can express the governing equations in terms of the so-called pressure equation: -@f{eqnarray*} +@f{eqnarray*}{ - \nabla \cdot (\mathbf{K}\lambda(S) \nabla p)= q. @f} Here, $q$ is the sum source term, and @@ -80,12 +80,12 @@ The second part of the equations is the description of the dynamics of the saturation, i.e., how the relative concentration of the two fluids changes with time. The saturation equation for the displacing fluid (water) is given by the following conservation law: -@f{eqnarray*} +@f{eqnarray*}{ S_{t} + \nabla \cdot (F(S) \mathbf{u}) = q_{w}, @f} which can be rewritten by using the product rule of the divergence operator in the previous equation: -@f{eqnarray*} +@f{eqnarray*}{ S_{t} + F(S) \left[\nabla \cdot \mathbf{u}\right] + \mathbf{u} \cdot \left[ \nabla F(S)\right] = S_{t} + F(S) q + \mathbf{u} \cdot \nabla F(S) = q_{w}. @@ -104,7 +104,7 @@ where the fractional flow is often parameterized via the (heuristic) expression @f] Putting it all together yields the saturation equation in the following, advected form: -@f{eqnarray*} +@f{eqnarray*}{ S_{t} + \mathbf{u} \cdot \nabla F(S) = 0, @f} where $\mathbf u$ is the total velocity @@ -124,7 +124,7 @@ $\mathbf u (1-F'(S))$. $F(S)$ is consequently often referred to as the fractional flow. In summary, what we get are the following two equations: -@f{eqnarray*} +@f{eqnarray*}{ - \nabla \cdot (\mathbf{K}\lambda(S) \nabla p) &=& q \qquad \textrm{in}\ \Omega\times[0,T], \\ @@ -169,7 +169,7 @@ In the reservoir simulation community, it is common to solve the equations derived above by going back to the first order, mixed formulation. To this end, we re-introduce the total velocity $\mathbf u$ and write the equations in the following form: -@f{eqnarray*} +@f{eqnarray*}{ \mathbf{u}+\mathbf{K}\lambda(S) \nabla p&=&0 \\ \nabla \cdot\mathbf{u} &=& q \\ S_{t} + \mathbf{u} \cdot \nabla F(S) &=& 0. @@ -194,7 +194,7 @@ L. Stone and A. O. Gardner Jr: Analysis of gas-cap or dissolved-gas reservoirs, Trans. SPE AIME, 222 (1961), pp. 92-104). In a slightly modified form, this algorithm can be written as follows: for each time step, solve -@f{eqnarray*} +@f{eqnarray*}{ \mathbf{u}^{n+1}+\mathbf{K}\lambda(S^n) \nabla p^{n+1}&=&0 \\ \nabla \cdot\mathbf{u}^{n+1} &=& q^{n+1} \\ \frac {S^{n+1}-S^n}{\triangle t} + \mathbf{u}^{n+1} \cdot \nabla F(S^n) &=& 0, @@ -211,7 +211,7 @@ splitting" method. step-58 has a long description of the idea behind this.) We can then state the problem in weak form as follows, by multiplying each equation with test functions $\mathbf v$, $\phi$, and $\sigma$ and integrating terms by parts: -@f{eqnarray*} +@f{eqnarray*}{ \left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{u}^{n+1},\mathbf v\right)_\Omega - (p^{n+1}, \nabla\cdot\mathbf v)_\Omega &=& - (p^{n+1}, \mathbf v)_{\partial\Omega} @@ -223,7 +223,7 @@ the boundary $\partial\Omega$ as boundary values for our problem. $\mathbf n$ denotes the unit outward normal vector to $\partial K$, as usual. For the saturation equation, we obtain after integrating by parts -@f{eqnarray*} +@f{eqnarray*}{ (S^{n+1}, \sigma)_\Omega - \triangle t @@ -238,7 +238,7 @@ For the saturation equation, we obtain after integrating by parts @f} Using the fact that $\nabla \cdot \mathbf{u}^{n+1}=q^{n+1}$, we can rewrite the cell term to get an equation as follows: -@f{eqnarray*} +@f{eqnarray*}{ (S^{n+1}, \sigma)_\Omega - \triangle t @@ -272,7 +272,7 @@ terms on the interfaces between cells, since discontinuous functions are not really defined there. In particular, we have to give a meaning to the last term on the left hand side of the saturation equation. To this end, let us define that we want to evaluate it in the following sense: -@f{eqnarray*} +@f{eqnarray*}{ &&\left(F(S^n) (\mathbf n \cdot \mathbf{u}^{n+1}), \sigma\right)_{\partial K} \\ &&\qquad = @@ -321,7 +321,7 @@ B & 0 & 0\\ where the individual matrices and vectors are defined as follows using shape functions $\mathbf v_i$ (of type Raviart Thomas $RT_k$) for velocities and $\phi_i$ (of type $DGQ_k$) for both pressures and saturations: -@f{eqnarray*} +@f{eqnarray*}{ M^u(S^n)_{ij} &=& \left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{v}_i,\mathbf v_j\right)_\Omega, @@ -471,7 +471,7 @@ boundary conditions for the saturation on the inflow part of the boundary, \mathbf{n} \cdot \mathbf{u}(\mathbf{x},t) < 0\}. @f] On this inflow boundary, we impose the following saturation values: -@f{eqnarray} +@f{eqnarray}{ S(\mathbf{x},t) = 1 & \textrm{on}\ \Gamma_{in}\cap\{x_1=0\}, \\ S(\mathbf{x},t) = 0 & \textrm{on}\ \Gamma_{in}\backslash \{x_1=0\}. @@ -540,7 +540,7 @@ functions introduced at the end of the results section of @ref step_20 linear solvers will no longer converge properly.
  • A function that models a somewhat random medium. Here, we choose - @f{eqnarray*} + @f{eqnarray*}{ k(\mathbf x) &=& \min \left\{ \max \left\{ \sum_{i=1}^N \sigma_i(\mathbf{x}), 0.01 \right\}, 4\right\}, diff --git a/examples/step-22/doc/intro.dox b/examples/step-22/doc/intro.dox index ec07c2bd89..4e6cb841e6 100644 --- a/examples/step-22/doc/intro.dox +++ b/examples/step-22/doc/intro.dox @@ -18,7 +18,7 @@ California Institute of Technology. This program deals with the Stokes system of equations which reads as follows in non-dimensionalized form: -@f{eqnarray*} +@f{eqnarray*}{ -2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& \textbf{f}, \\ -\textrm{div}\; \textbf{u} &=& 0, @@ -42,14 +42,14 @@ we will focus on the simpler Stokes system. Note that when deriving the more general compressible Navier-Stokes equations, the diffusion is modeled as the divergence of the stress tensor -@f{eqnarray*} +@f{eqnarray*}{ \tau = - \mu \left(2\varepsilon(\textbf{u}) - \frac{2}{3}\nabla \cdot \textbf{u} I\right), @f} where $\mu$ is the viscosity of the fluid. With the assumption of $\mu=1$ (assume constant viscosity and non-dimensionalize the equation by dividing out $\mu$) and assuming incompressibility ($\textrm{div}\; \textbf{u}=0$), we arrive at the formulation from above: -@f{eqnarray*} +@f{eqnarray*}{ \textrm{div}\; \tau = -2\textrm{div}\;\varepsilon(\textbf{u}). @f} A different formulation uses the Laplace operator ($-\triangle \textbf{u}$) @@ -57,14 +57,14 @@ instead of the symmetrized gradient. A big difference here is that the different components of the velocity do not couple. If you assume additional regularity of the solution $\textbf{u}$ (second partial derivatives exist and are continuous), the formulations are equivalent: -@f{eqnarray*} +@f{eqnarray*}{ \textrm{div}\; \tau = -2\textrm{div}\;\varepsilon(\textbf{u}) = -\triangle \textbf{u} - \nabla \cdot (\nabla\textbf{u})^T = -\triangle \textbf{u}. @f} This is because the $i$th entry of $\nabla \cdot (\nabla\textbf{u})^T$ is given by: -@f{eqnarray*} +@f{eqnarray*}{ [\nabla \cdot (\nabla\textbf{u})^T]_i = \sum_j \frac{\partial}{\partial x_j} [(\nabla\textbf{u})^T]_{i,j} = \sum_j \frac{\partial}{\partial x_j} [(\nabla\textbf{u})]_{j,i} @@ -92,7 +92,7 @@ toplevel overview of this topic can be found in the @ref vector_valued module. The weak form of the equations is obtained by writing it in vector form as -@f{eqnarray*} +@f{eqnarray*}{ \begin{pmatrix} {-2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p} \\ @@ -108,7 +108,7 @@ form as forming the dot product from the left with a vector-valued test function $\phi = \begin{pmatrix}\textbf{v} \\ q\end{pmatrix}$ and integrating over the domain $\Omega$, yielding the following set of equations: -@f{eqnarray*} +@f{eqnarray*}{ (\mathrm v, -2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p)_{\Omega} - @@ -132,7 +132,7 @@ weak formulation that has only one derivative on both test and trial function. In the current context, we integrate by parts the second term: -@f{eqnarray*} +@f{eqnarray*}{ (\textbf{v}, -2\; \textrm{div}\; \varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega} @@ -142,7 +142,7 @@ In the current context, we integrate by parts the second term: (\textbf{v}, \textbf{f})_\Omega. @f} Likewise, we integrate by parts the first term to obtain -@f{eqnarray*} +@f{eqnarray*}{ (\nabla \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\Omega} - (\textbf{n} \otimes \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\partial\Omega} @@ -155,7 +155,7 @@ Likewise, we integrate by parts the first term to obtain @f} where the scalar product between two tensor-valued quantities is here defined as -@f{eqnarray*} +@f{eqnarray*}{ (\nabla \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\Omega} = 2 \int_\Omega \sum_{i,j=1}^d \frac{\partial v_j}{\partial x_i} @@ -170,7 +170,7 @@ $\nabla\textbf{v}$ and a symmetric tensor like $\varepsilon(\textbf{u})$ equals the scalar product between the symmetrized forms of the two, we can also write the bilinear form above as follows: -@f{eqnarray*} +@f{eqnarray*}{ (\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega} - (\textbf{n} \otimes \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\partial\Omega} @@ -183,7 +183,7 @@ above as follows: @f} We will deal with the boundary terms in the next section, but it is already clear from the domain terms -@f{eqnarray*} +@f{eqnarray*}{ (\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - @@ -205,13 +205,13 @@ possibilities for imposing boundary conditions: $\Gamma_D\subset\partial\Omega$ we may impose Dirichlet conditions on the velocity $\textbf u$: - @f{eqnarray*} + @f{eqnarray*}{ \textbf u = \textbf g_D \qquad\qquad \textrm{on}\ \Gamma_D. @f} Because test functions $\textbf{v}$ come from the tangent space of the solution variable, we have that $\textbf{v}=0$ on $\Gamma_D$ and consequently that - @f{eqnarray*} + @f{eqnarray*}{ -(\textbf{n} \otimes \mathrm v, 2\; \varepsilon(\textbf{u}))_{\Gamma_D} + @@ -229,7 +229,7 @@ possibilities for imposing boundary conditions:
  • Neumann-type or natural boundary conditions: On the rest of the boundary $\Gamma_N=\partial\Omega\backslash\Gamma_D$, let us re-write the boundary terms as follows: - @f{eqnarray*} + @f{eqnarray*}{ -(\textbf{n} \otimes \mathrm v, 2\; \varepsilon(\textbf{u}))_{\Gamma_N} + @@ -262,14 +262,14 @@ possibilities for imposing boundary conditions: @f} In other words, on the Neumann part of the boundary we can prescribe values for the total stress: - @f{eqnarray*} + @f{eqnarray*}{ \textbf{n}\cdot [p \textbf{I} - 2\; \varepsilon(\textbf{u})] = \textbf g_N \qquad\qquad \textrm{on}\ \Gamma_N. @f} If the boundary is subdivided into Dirichlet and Neumann parts $\Gamma_D,\Gamma_N$, this then leads to the following weak form: - @f{eqnarray*} + @f{eqnarray*}{ (\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - @@ -283,13 +283,13 @@ possibilities for imposing boundary conditions:
  • Robin-type boundary conditions: Robin boundary conditions are a mixture of Dirichlet and Neumann boundary conditions. They would read - @f{eqnarray*} + @f{eqnarray*}{ \textbf{n}\cdot [p \textbf{I} - 2\; \varepsilon(\textbf{u})] = \textbf S \textbf u \qquad\qquad \textrm{on}\ \Gamma_R, @f} with a rank-2 tensor (matrix) $\textbf S$. The associated weak form is - @f{eqnarray*} + @f{eqnarray*}{ (\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - @@ -309,7 +309,7 @@ possibilities for imposing boundary conditions: dim-1 components of the velocity. The remaining component can be constrained by requiring that the normal component of the normal stress be zero, yielding the following set of boundary conditions: - @f{eqnarray*} + @f{eqnarray*}{ \textbf u_{\textbf t} &=& 0, \\ \textbf n \cdot \left(\textbf{n}\cdot [p \textbf{I} - 2\; @@ -329,7 +329,7 @@ possibilities for imposing boundary conditions: stratified by different densities but that both have small enough viscosities to not introduce much tangential stress on each other). In formulas, this means that - @f{eqnarray*} + @f{eqnarray*}{ \textbf{n}\cdot\textbf u &=& 0, \\ (\textbf 1-\textbf n\otimes\textbf n) @@ -353,7 +353,7 @@ As developed above, the weak form of the equations with Dirichlet and Neumann boundary conditions on $\Gamma_D$ and $\Gamma_N$ reads like this: find $\textbf u\in \textbf V_g = \{\varphi \in H^1(\Omega)^d: \varphi_{\Gamma_D}=\textbf g_D\}, p\in Q=L^2(\Omega)$ so that -@f{eqnarray*} +@f{eqnarray*}{ (\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - @@ -387,7 +387,7 @@ pressures. This then leads to the following discrete problem: find $\textbf u_h,p_h$ so that -@f{eqnarray*} +@f{eqnarray*}{ (\varepsilon(\textbf{v}_h), 2\; \varepsilon(\textbf u_h))_{\Omega} - (\textrm{div}\; \textbf{v}_h, p_h)_{\Omega} - @@ -408,7 +408,7 @@ vector_valued module. The weak form of the discrete equations naturally leads to the following linear system for the nodal values of the velocity and pressure fields: -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) @@ -423,11 +423,11 @@ linear system for the nodal values of the velocity and pressure fields: Like in step-20 and step-21, we will solve this system of equations by forming the Schur complement, i.e. we will first find the solution $P$ of -@f{eqnarray*} +@f{eqnarray*}{ BA^{-1}B^T P &=& BA^{-1} F - G, \\ @f} and then -@f{eqnarray*} +@f{eqnarray*}{ AU &=& F - B^TP. @f} The way we do this is pretty much exactly like we did in these previous @@ -611,7 +611,7 @@ complementation preconditioner in the results section of this program.

    A note on the structure of the linear system

    Above, we have claimed that the linear system has the form -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) @@ -645,7 +645,7 @@ a few entries on the diagonal, one for each constrained pressure degree of freedom, and a correct description of the linear system we have to solve is that it has the form -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & D_c \end{array}\right) @@ -732,7 +732,7 @@ year at most), leaving a crack in the earth crust that is filled with magma from below. Without trying to be entirely realistic, we model this situation by solving the following set of equations and boundary conditions on the domain $\Omega=[-2,2]\times[0,1]\times[-1,0]$: -@f{eqnarray*} +@f{eqnarray*}{ -2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& 0, \\ -\textrm{div}\; \textbf{u} &=& 0, diff --git a/examples/step-22/doc/results.dox b/examples/step-22/doc/results.dox index d9f1a109e6..b11108f4f3 100644 --- a/examples/step-22/doc/results.dox +++ b/examples/step-22/doc/results.dox @@ -357,14 +357,14 @@ complement, there is no other possibility. The alternative is to attack the block system at once and use an approximate Schur complement as efficient preconditioner. The idea is as follows: If we find a block preconditioner $P$ such that the matrix -@f{eqnarray*} +@f{eqnarray*}{ P^{-1}\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) @f} is simple, then an iterative solver with that preconditioner will converge in a few iterations. Using the Schur complement $S = B A^{-1} B^T$, one finds that -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} = \left(\begin{array}{cc} @@ -372,7 +372,7 @@ few iterations. Using the Schur complement $S = B A^{-1} B^T$, one finds that \end{array}\right) @f} would appear to be a good choice since -@f{eqnarray*} +@f{eqnarray*}{ P^{-1}\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) diff --git a/examples/step-23/doc/intro.dox b/examples/step-23/doc/intro.dox index 1310890811..6c52ea052c 100644 --- a/examples/step-23/doc/intro.dox +++ b/examples/step-23/doc/intro.dox @@ -13,7 +13,7 @@ solutions called solitons. The wave equation in its prototypical form reads as follows: find $u(x,t), x\in\Omega, t\in[0,T]$ that satisfies -@f{eqnarray*} +@f{eqnarray*}{ \frac{\partial^2 u}{\partial t^2} - \Delta u &=& f @@ -127,7 +127,7 @@ this end, let us introduce @f] and call this variable the velocity for obvious reasons. We can then reformulate the original wave equation as follows: -@f{eqnarray*} +@f{eqnarray*}{ \frac{\partial u}{\partial t} - v @@ -289,7 +289,7 @@ system: \right], \f} where -@f{eqnarray*} +@f{eqnarray*}{ M^n_{ij} &=& (\phi_i^n, \phi_j^n), \\ A^n_{ij} &=& (\nabla\phi_i^n, \nabla\phi_j^n), @@ -398,7 +398,7 @@ step_24 "step-24" shows a better way how to keep these things in sync. Although the program has all the hooks to deal with nonzero initial and boundary conditions and body forces, we take a simple case where the domain is a square $[-1,1]^2$ and -@f{eqnarray*} +@f{eqnarray*}{ f &=& 0, \\ u_0 &=& 0, diff --git a/examples/step-24/doc/intro.dox b/examples/step-24/doc/intro.dox index 9d4f12848b..1f8e2de6aa 100644 --- a/examples/step-24/doc/intro.dox +++ b/examples/step-24/doc/intro.dox @@ -71,7 +71,7 @@ where $\lambda = - \frac{\beta}{C_p}$. This somewhat strange equation with the derivative of a Dirac delta function on the right hand side can be rewritten as an initial value problem as follows: -@f{eqnarray*} +@f{eqnarray*}{ \Delta \bar{p}- \frac{1}{c_0^2} \frac{\partial^2 \bar{p}}{\partial t^2} & = & 0 \\ \bar{p}(0,\mathbf r) &=& c_0^2 \lambda a(\mathbf r) = b(\mathbf r) \\ @@ -120,12 +120,12 @@ v = \frac{\partial\bar{p}}{\partial t} With the second variable, one then transforms the forward problem into two separate equations: -@f{eqnarray*} +@f{eqnarray*}{ \bar{p}_{t} - v & = & 0 \\ \Delta\bar{p} - \frac{1}{c_0^2}\,v_{t} & = & f @f} with initial conditions: -@f{eqnarray*} +@f{eqnarray*}{ \bar{p}(0,\mathbf r) & = & b(r) \\ v(0,\mathbf r)=\bar{p}_t(0,\mathbf r) & = & 0. @f} @@ -135,7 +135,7 @@ to the thermoacoustic problem $f=0$. The semi-discretized, weak version of this model, using the general $\theta$ scheme introduced in step-23 is then: -@f{eqnarray*} +@f{eqnarray*}{ \left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_\Omega- \left(\theta v^{n}+(1-\theta)v^{n-1},\phi\right)_\Omega & = & 0 \\ -\left(\nabla((\theta\bar{p}^n+(1-\theta)\bar{p}^{n-1})),\nabla\phi\right)_\Omega- @@ -154,7 +154,7 @@ absorbing boundary conditions are incorporated into the weak form by using From this we obtain the discrete model by introducing a finite number of shape functions, and get -@f{eqnarray*} +@f{eqnarray*}{ M\bar{p}^{n}-k \theta M v^n & = & M\bar{p}^{n-1}+k (1-\theta)Mv^{n-1},\\ (-c_0^2k \theta A-c_0 B)\bar{p}^n-Mv^{n} & = & @@ -198,7 +198,7 @@ G_2 \\ By simple transformations, one then obtains two equations for the pressure potential and its derivative, just as in the previous tutorial program: -@f{eqnarray*} +@f{eqnarray*}{ (M+(k\,\theta\,c_{0})^{2}A+c_0k\theta B)\bar{p}^{n} & = & G_{1}+(k\, \theta)G_{2}-(c_0k)^2\theta (\theta F^{n}+(1-\theta)F^{n-1}) \\ Mv^n & = & -(c_0^2\,k\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 - diff --git a/examples/step-25/doc/intro.dox b/examples/step-25/doc/intro.dox index db22da07a2..15d6b222d3 100644 --- a/examples/step-25/doc/intro.dox +++ b/examples/step-25/doc/intro.dox @@ -151,7 +151,7 @@ letter the vector of coefficients (in the nodal basis) of a function denoted by the same letter in lower case; e.g., $u^n = \sum_{i=1}^N U^n_i \varphi_i$ where $U^n \in {R}^N$ and $u^n \in H^1(\Omega)$. Thus, the finite-dimensional version of the variational formulation requires that we solve the following matrix equations at each time step: -@f{eqnarray*} +@f{eqnarray*}{ F_h'(U^{n,l})\delta U^{n,l} &=& -F_h(U^{n,l}), \qquad U^{n,l+1} = U^{n,l} + \delta U^{n,l}, \qquad U^{n,0} = U^{n-1}; \\ MV^n &=& MV^{n-1} - k \theta AU^n -k (1-\theta) AU^{n-1} - k S(u^n,u^{n-1}). diff --git a/examples/step-26/doc/intro.dox b/examples/step-26/doc/intro.dox index 0fbd2b74d4..d7c7fbefd8 100644 --- a/examples/step-26/doc/intro.dox +++ b/examples/step-26/doc/intro.dox @@ -6,7 +6,7 @@ This program implements the heat equation -@f{align*} +@f{align*}{ \frac{\partial u(\mathbf x, t)}{\partial t} - \Delta u(\mathbf x, t) @@ -32,7 +32,7 @@ mesh every few time steps, etc. Our goal here will be to solve the equations above using the theta-scheme that discretizes the equation in time using the following approach, where we would like $u^n(\mathbf x)$ to approximate $u(\mathbf x, t_n)$ at some time $t_n$: -@f{align*} +@f{align*}{ \frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n} - \left[ @@ -59,7 +59,7 @@ does, by multiplying with test functions, integrating by parts, and then restricting everything to a finite dimensional subspace. This yields the following set of fully discrete equations after multiplying through with $k_n$: -@f{align*} +@f{align*}{ M U^n-MU^{n-1} + k_n \left[ @@ -78,7 +78,7 @@ $k_n$: where $M$ is the @ref GlossMassMatrix "mass matrix" and $A$ is the @ref GlossStiffnessMatrix "stiffness matrix" that results from discretizing the Laplacian. Bringing all known quantities to the right hand side yields the linear system we have to solve in every step: -@f{align*} +@f{align*}{ (M + k_n \theta A) U^n @@ -155,7 +155,7 @@ compared to the stationary case. Let us go through them in turn:
  • Test functions from different meshes: Let us consider again the semi-discrete equations we have written down above: - @f{align*} + @f{align*}{ \frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n} - \left[ @@ -172,14 +172,14 @@ compared to the stationary case. Let us go through them in turn: @f} We can here consider $u^{n-1}$ as data since it has presumably been computed before. Now, let us replace - @f{align*} + @f{align*}{ u^n(\mathbf x)\approx u_h^n(\mathbf x) = \sum_j U^n \varphi_j(\mathbf x), @f} multiply with test functions $\varphi_i(\mathbf x)$ and integrate by parts where necessary. In a process as outlined above, this would yield - @f{align*} + @f{align*}{ \sum_j (M + @@ -202,7 +202,7 @@ compared to the stationary case. Let us go through them in turn: $u^{n-1}$ are different! This pertains to the terms on the right hand side, the first of which we could more clearly write as (the second follows the same pattern) - @f{align*} + @f{align*}{ (\varphi_i, u_h^{n-1}) = (\varphi_i^n, u_h^{n-1}) @@ -232,7 +232,7 @@ compared to the stationary case. Let us go through them in turn: whole situation by interpolating the solution from the old to the new mesh every time we adapt the mesh. In other words, rather than solving the equations above, we instead solve the problem - @f{align*} + @f{align*}{ \sum_j (M + @@ -289,7 +289,7 @@ equations without coefficients (or constant coefficients) like the one here, there is a fairly standard protocol that rests on the following observation: if you choose as your domain a square $[0,1]^2$ (or, with slight modifications, a rectangle), then the exact solution can be written as -@f{align*} +@f{align*}{ u(x,y,t) = a(t) \sin(n_x \pi x) \sin(n_y \pi y) @f} (with integer constants $n_x,n_y$) @@ -303,7 +303,7 @@ As an example, let us consider the situation where we have $u_0(x,y)=\sin(n_x \pi x) \sin(n_x \pi y)$ and $f(x,y,t)=0$. With the claim (ansatz) of the form for $u(x,y,t)$ above, we get that -@f{align*} +@f{align*}{ \left(\frac{\partial}{\partial t} -\Delta\right) u(x,y,t) &= @@ -314,12 +314,12 @@ $u(x,y,t)$ above, we get that \left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y). @f} For this to be equal to $f(x,y,t)=0$, we need that -@f{align*} +@f{align*}{ a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 0 @f} and due to the initial conditions, $a(0)=1$. This differential equation can be integrated to yield -@f{align*} +@f{align*}{ a(t) = - e^{-(n_x^2+n_y^2)\pi^2 t}. @f} In other words, if the initial condition is a product of sines, then the @@ -341,7 +341,7 @@ If you have so verified that the time integrator is correct, take the situation where the right hand side is nonzero but the initial conditions are zero: $u_0(x,y)=0$ and $f(x,y,t)=\sin(n_x \pi x) \sin(n_x \pi y)$. Again, -@f{align*} +@f{align*}{ \left(\frac{\partial}{\partial t} -\Delta\right) u(x,y,t) &= @@ -352,12 +352,12 @@ $f(x,y,t)=\sin(n_x \pi x) \sin(n_x \pi y)$. Again, \left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y), @f} and for this to be equal to $f(x,y,t)$, we need that -@f{align*} +@f{align*}{ a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 1 @f} and due to the initial conditions, $a(0)=0$. Integrating this equation in time yields -@f{align*} +@f{align*}{ a(t) = \frac{1}{(n_x^2+n_y^2)\pi^2} \left[ 1 - e^{-(n_x^2+n_y^2)\pi^2 t} \right]. @f} @@ -380,7 +380,7 @@ become very smooth very quickly and then do not move very much any more. Rather, we here solve the equation on the L-shaped domain with zero Dirichlet boundary values and zero initial conditions, but as right hand side we choose -@f{align*} +@f{align*}{ f(\mathbf x, t) = \left\{ @@ -399,7 +399,7 @@ we choose \right. @f} Here, -@f{align*} +@f{align*}{ \chi_1(\mathbf x) &= \left\{ \begin{array}{ll} diff --git a/examples/step-26/doc/results.dox b/examples/step-26/doc/results.dox index db8010651f..3b4a2c3961 100644 --- a/examples/step-26/doc/results.dox +++ b/examples/step-26/doc/results.dox @@ -150,11 +150,11 @@ This is not what we would expect to happen (in nature). To get an idea of this behavior mathematically, let us consider a general, fully discrete problem: -@f{align*} +@f{align*}{ A u^{n} = B u^{n-1}. @f} The general form of the $i$th equation then reads: -@f{align*} +@f{align*}{ a_{ii} u^{n}_i &= b_{ii} u^{n-1}_i + \sum\limits_{j \in S_i} \left( b_{ij} u^{n-1}_j - a_{ij} u^{n}_j \right), @f} @@ -162,7 +162,7 @@ where $S_i$ is the set of degrees of freedom that DoF $i$ couples with (i.e., for which either the matrix $A$ or matrix $B$ has a nonzero entry at position $(i,j)$). If all coefficients fulfill the following conditions: -@f{align*} +@f{align*}{ a_{ii} &> 0, & b_{ii} &\geq 0, & a_{ij} &\leq 0, & b_{ij} &\geq 0, & \forall j &\in S_i, @@ -177,7 +177,7 @@ able to deduce conditions for the time step $k_n$. For the heat equation with the Crank-Nicolson scheme, Schatz et. al. have translated it to the following ones: -@f{align*} +@f{align*}{ (1 - \theta) k a_{ii} &\leq m_{ii},\qquad \forall i, & \theta k \left| a_{ij} \right| &\geq m_{ij},\qquad j \neq i, @@ -186,7 +186,7 @@ where $M = m_{ij}$ denotes the @ref GlossMassMatrix "mass matrix" and $A = a_{ij matrix with $a_{ij} \leq 0$ for $j \neq i$, respectively. With $a_{ij} \leq 0$, we can formulate bounds for the global time step $k$ as follows: -@f{align*} +@f{align*}{ k_{\text{max}} &= \frac{ 1 }{ 1 - \theta } \min\left( \frac{ m_{ii} }{ a_{ii} } \right),~ \forall i, & diff --git a/examples/step-28/doc/intro.dox b/examples/step-28/doc/intro.dox index 7745ee3d6b..400c426992 100644 --- a/examples/step-28/doc/intro.dox +++ b/examples/step-28/doc/intro.dox @@ -39,7 +39,7 @@ fluxes. Assume we have energy groups $g=1,\ldots,G$, where by convention we assume that the neutrons with the highest energy are in group 1 and those with the lowest energy in group $G$. Then the neutron flux of each group satisfies the following equations: -@f{eqnarray*} +@f{eqnarray*}{ \frac 1{v_g}\frac{\partial \phi_g(x,t)}{\partial t} &=& \nabla \cdot(D_g(x) \nabla \phi_g(x,t)) @@ -108,7 +108,7 @@ below). If we consider all energy groups at once, we may write above equations in the following operator form: -@f{eqnarray*} +@f{eqnarray*}{ \frac 1v \frac{\partial \phi}{\partial t} = -L\phi @@ -126,7 +126,7 @@ that $L$ is symmetric, whereas $F$ and $X$ are not. It is well known that this equation admits a stable solution if all eigenvalues of the operator $-L+F+X$ are negative. This can be readily seen by multiplying the equation by $\phi$ and integrating over the domain, leading to -@f{eqnarray*} +@f{eqnarray*}{ \frac 1{2v} \frac{\partial}{\partial t} \|\phi\|^2 = ((-L+F+X)\phi,\phi). @f} Stability means that the solution does not grow, i.e. we want the left hand @@ -136,7 +136,7 @@ not very desirable if a nuclear reactor produces neutron fluxes that grow exponentially, so eigenvalue analyses are the bread-and-butter of nuclear engineers. The main point of the program is therefore to consider the eigenvalue problem -@f{eqnarray*} +@f{eqnarray*}{ (L-F-X) \phi = \lambda \phi, @f} where we want to make sure that all eigenvalues are positive. Note that $L$, @@ -149,20 +149,20 @@ In nuclear engineering, one typically looks at a slightly different formulation of the eigenvalue problem. To this end, we do not just multiply with $\phi$ and integrate, but rather multiply with $\phi(L-X)^{-1}$. We then get the following evolution equation: -@f{eqnarray*} +@f{eqnarray*}{ \frac 1{2v} \frac{\partial}{\partial t} \|\phi\|^2_{(L-X)^{-1}} = ((L-X)^{-1}(-L+F+X)\phi,\phi). @f} Stability is then guaranteed if the eigenvalues of the following problem are all negative: -@f{eqnarray*} +@f{eqnarray*}{ (L-X)^{-1}(-L+F+X)\phi = \lambda_F \phi, @f} which is equivalent to the eigenvalue problem -@f{eqnarray*} +@f{eqnarray*}{ (L-X)\phi = \frac 1{\lambda_F+1} F \phi. @f} The typical formulation in nuclear engineering is to write this as -@f{eqnarray*} +@f{eqnarray*}{ (L-X) \phi = \frac 1{k_{\mathrm{eff}}} F \phi, @f} where $k_{\mathrm{eff}}=\frac 1{\lambda^F+1}$. @@ -211,7 +211,7 @@ like this: and $k_{\mathrm{eff}}^{(0)}$ and let $n=1$.
  • Define the so-called fission source by - @f{eqnarray*} + @f{eqnarray*}{ s_f^{(n-1)}(x) = \frac{1}{k_{\mathrm{eff}}^{(n-1)}} @@ -219,7 +219,7 @@ like this: @f}
  • Solve for all group fluxes $\phi_g,g=1,\ldots,G$ using - @f{eqnarray*} + @f{eqnarray*}{ -\nabla \cdot D_g\nabla \phi_g^{(n)} + \Sigma_{r,g}\phi_g^{(n)} @@ -232,7 +232,7 @@ like this: @f}
  • Update - @f{eqnarray*} + @f{eqnarray*}{ k_{\mathrm{eff}}^{(n)} = \sum_{g'=1}^G @@ -292,7 +292,7 @@ again use the a posteriori error estimator by Kelly, Gago, Zienkiewicz and Babuska which approximates the error per cell by integrating the jump of the gradient of the solution along the faces of each cell. Using this, we obtain indicators -@f{eqnarray*} +@f{eqnarray*}{ \eta_{g,K}, \qquad g=1,2,\ldots,G,\qquad K\in{\cal T}_g, @f} where ${\cal T}_g$ is the triangulation used in the solution of @@ -315,7 +315,7 @@ there, however, and simply assume that all energy groups are equally important, and will therefore normalize the error indicators $\eta_{g,K}$ for group $g$ by the maximum of the solution $\phi_g$. We then refine the cells whose errors satisfy -@f{eqnarray*} +@f{eqnarray*}{ \frac{\eta_{g,K}}{\|\phi_g\|_\infty} > \alpha_1 @@ -323,7 +323,7 @@ whose errors satisfy \frac{\eta_{g,K}}{\|\phi_g\|_\infty}} @f} and coarsen the cells where -@f{eqnarray*} +@f{eqnarray*}{ \frac{\eta_{g,K}}{\|\phi_g\|_\infty} < \alpha_2 @@ -352,7 +352,7 @@ form the weak form for the equation to compute $\phi_g^{(n)}$ as usual by multiplication with test functions $\varphi_g^i$ defined on the mesh for energy group $g$; in the process, we have to compute the right hand side vector that contains terms of the following form: -@f{eqnarray*} +@f{eqnarray*}{ F_i = \int_\Omega f(x) \varphi_g^i(x) \phi_{g'}(x) \ dx, @f} where $f(x)$ is one of the coefficient functions $\Sigma_{s,g'\to g}$ or @@ -362,7 +362,7 @@ the mesh for energy group $g'$, i.e. it can be expanded as $\phi_{g'}(x)=\sum_j\phi_{g'}^j \varphi_{g'}^j(x)$, with basis functions $\varphi_{g'}^j(x)$ defined on mesh $g'$. The contribution to the right hand side can therefore be written as -@f{eqnarray*} +@f{eqnarray*}{ F_i = \sum_j \left\{\int_\Omega f(x) \varphi_g^i(x) \varphi_{g'}^j(x) \ dx \right\} \phi_{g'}^j , @f} @@ -389,7 +389,7 @@ GridTools::get_finest_common_cells that computes exactly this set of cells that are active on at least one of two meshes. With this, we can write above integral as follows: -@f{eqnarray*} +@f{eqnarray*}{ F_i = \sum_{K \in {\cal T}_g \cap {\cal T}_{g'}} @@ -436,7 +436,7 @@ these three cases, as follows: of $\phi_g^i$ to child cell $K_c$ into the basis functions defined on that child cell (i.e. on cells on which the basis functions $\varphi_{g'}^l$ are defined): - @f{eqnarray*} + @f{eqnarray*}{ \phi_g^i|_{K_c} = B_c^{il} \varphi_{g'}^l|_{K_c}. @f} Here, and in the following, summation over indices appearing twice is @@ -445,7 +445,7 @@ these three cases, as follows: Then we can write the contribution of cell $K$ to the right hand side component $F_i$ as - @f{eqnarray*} + @f{eqnarray*}{ F_i|_K &=& \left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x) @@ -458,7 +458,7 @@ these three cases, as follows: \ dx \right\} \phi_{g'}^j. @f} In matrix notation, this can be written as - @f{eqnarray*} + @f{eqnarray*}{ F_i|_K = \sum_{0\le c<2^{\texttt{dim}}} @@ -478,11 +478,11 @@ these three cases, as follows: until we find an active cell. We then have to sum up all the contributions from all the children, grandchildren, etc, of cell $K$, with contributions of the form - @f{eqnarray*} + @f{eqnarray*}{ F_i|_{K_{cc'}} = (B_cB_{c'} M_{K_{cc'}})^{ij} \phi_{g'}^j, @f} or - @f{eqnarray*} + @f{eqnarray*}{ F_i|_{K_{cc'c''}} = (B_c B_{c'} B_{c''}M_{K_{cc'c''}})^{ij} \phi_{g'}^j, @f} @@ -499,7 +499,7 @@ these three cases, as follows: than $\varphi_g^i$ as before. This of course works in exactly the same way. If the children of $K$ are active on mesh $g$, then leading to the expression - @f{eqnarray*} + @f{eqnarray*}{ F_i|_K &=& \left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x) @@ -512,7 +512,7 @@ these three cases, as follows: \ dx \right\} \phi_{g'}^j. @f} In matrix notation, this expression now reads as - @f{eqnarray*} + @f{eqnarray*}{ F_i|_K = \sum_{0\le c<2^{\texttt{dim}}} @@ -525,11 +525,11 @@ these three cases, as follows: @f} and correspondingly for cases where cell $K$ is refined more than once on mesh $g$: - @f{eqnarray*} + @f{eqnarray*}{ F_i|_{K_{cc'}} = (M_{K_{cc'}} B_{c'}^T B_c^T)^{ij} \phi_{g'}^j, @f} or - @f{eqnarray*} + @f{eqnarray*}{ F_i|_{K_{cc'c''}} = (M_{K_{cc'c''}} B_{c''}^T B_{c'}^T B_c^T)^{ij} \phi_{g'}^j, @f} @@ -545,13 +545,13 @@ $(f \varphi_g^i, \varphi_{g'}^j)_K$ onto child cells, and then finally forming the inner product (the mass matrix) on the final cell. To make the symmetry in these cases more obvious, we can write them like this: for case (ii), we have -@f{eqnarray*} +@f{eqnarray*}{ F_i|_{K_{cc'\cdots c^{(k)}}} = [B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}}]^{ij} \phi_{g'}^j, @f} whereas for case (iii) we get -@f{eqnarray*} +@f{eqnarray*}{ F_i|_{K_{cc'\cdots c^{(k)}}} = [(B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}})^T]^{ij} \phi_{g'}^j, diff --git a/examples/step-29/doc/intro.dox b/examples/step-29/doc/intro.dox index 163ab0510e..43969c28de 100644 --- a/examples/step-29/doc/intro.dox +++ b/examples/step-29/doc/intro.dox @@ -73,7 +73,7 @@ dependency of amplitude and phase (relative to the source) of the waves of frequency ${\omega}$, with the amplitude being the quantity that we are interested in. By plugging this form of the solution into the wave equation, we see that for $u$ we have -@f{eqnarray*} +@f{eqnarray*}{ -\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\ u(x) &=& 1, \qquad x\in\Gamma_1. @f} @@ -101,7 +101,7 @@ fraction of the wave will be reflected back into the domain. If we are willing to accept this as a sufficient approximation to an absorbing boundary we finally arrive at the following problem for $u$: -@f{eqnarray*} +@f{eqnarray*}{ -\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\ c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\ u &=& 1, \qquad x\in\Gamma_1. @@ -115,7 +115,7 @@ imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system: -@f{eqnarray*} +@f{eqnarray*}{ \left.\begin{array}{ccc} -\omega^2 v - c^2\Delta v &=& 0 \quad\\ -\omega^2 w - c^2\Delta w &=& 0 \quad @@ -135,7 +135,7 @@ of $v$ and $w$ we have the following system: For test functions $\phi,\psi$ with $\phi|_{\Gamma_1}=\psi|_{\Gamma_1}=0$, after the usual multiplication, integration over $\Omega$ and applying integration by parts, we get the weak formulation -@f{eqnarray*} +@f{eqnarray*}{ -\omega^2 \langle \phi, v \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \phi, \nabla v \rangle_{\mathrm{L}^2(\Omega)} - c \omega \langle \phi, w \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0, \\ diff --git a/examples/step-3/doc/intro.dox b/examples/step-3/doc/intro.dox index f61fc92daf..0bbafd7184 100644 --- a/examples/step-3/doc/intro.dox +++ b/examples/step-3/doc/intro.dox @@ -9,7 +9,7 @@ This is the first example where we actually use finite elements to compute something. We will solve a simple version of Poisson's equation with zero boundary values, but a nonzero right hand side: -@f{align*} +@f{align*}{ -\Delta u &= f \qquad\qquad & \text{in}\ \Omega, \\ u &= 0 \qquad\qquad & \text{on}\ \partial\Omega. @@ -27,11 +27,11 @@ of the equation above, which we obtain by multiplying the equation by a test function $\varphi$ from the left (we will come back to the reason for multiplying from the left and not from the right below) and integrating over the domain $\Omega$: -@f{align*} +@f{align*}{ -\int_\Omega \varphi \Delta u = \int_\Omega \varphi f. @f} This can be integrated by parts: -@f{align*} +@f{align*}{ \int_\Omega \nabla\varphi \cdot \nabla u - \int_{\partial\Omega} \varphi \mathbf{n}\cdot \nabla u @@ -41,7 +41,7 @@ The test function $\varphi$ has to satisfy the same kind of boundary conditions (in mathematical terms: it needs to come from the tangent space of the set in which we seek the solution), so on the boundary $\varphi=0$ and consequently the weak form we are looking for reads -@f{align*} +@f{align*}{ (\nabla\varphi, \nabla u) = (\varphi, f), @f} @@ -83,7 +83,7 @@ we need the following: Through these steps, we now have a set of functions $\varphi_i$, and we can define the weak form of the discrete problem: Find a function $u_h$, i.e., find the expansion coefficients $U_j$ mentioned above, so that -@f{align*} +@f{align*}{ (\nabla\varphi_i, \nabla u_h) = (\varphi_i, f), \qquad\qquad @@ -106,7 +106,7 @@ With this, the problem reads: Find a vector $U$ so that A U = F, @f} where the matrix $A$ and the right hand side $F$ are defined as -@f{align*} +@f{align*}{ A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j), \\ F_i &= (\varphi_i, f). @@ -119,12 +119,12 @@ Before we move on with describing how these quantities can be computed, note that if we had multiplied the original equation from the right by a test function rather than from the left, then we would have obtained a linear system of the form -@f{align*} +@f{align*}{ U^T A = F^T @f} with a row vector $F^T$. By transposing this system, this is of course equivalent to solving -@f{align*} +@f{align*}{ A^T U = F @f} which here is the same as above since $A=A^T$. But in general is not, @@ -151,7 +151,7 @@ takes to make that happen: most commonly done using quadrature, i.e. the integrals are replaced by a weighted sum over a set of *quadrature points* on each cell. That is, we first split the integral over $\Omega$ into integrals over all cells, - @f{align*} + @f{align*}{ A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j) = \sum_{K \in {\mathbb T}} \int_K \nabla\varphi_i \cdot \nabla \varphi_j, \\ @@ -159,7 +159,7 @@ takes to make that happen: = \sum_{K \in {\mathbb T}} \int_K \varphi_i f, @f} and then approximate each cell's contribution by quadrature: - @f{align*} + @f{align*}{ A^K_{ij} &= \int_K \nabla\varphi_i \cdot \nabla \varphi_j \approx diff --git a/examples/step-31/doc/intro.dox b/examples/step-31/doc/intro.dox index cd3b7d74c6..88bb3b561e 100644 --- a/examples/step-31/doc/intro.dox +++ b/examples/step-31/doc/intro.dox @@ -27,7 +27,7 @@ down with gravity. In cases where the fluid moves slowly enough such that inertial effects can be neglected, the equations that describe such behavior are the Boussinesq equations that read as follows: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& -\rho\; \beta \; T\; \mathbf{g}, \\ @@ -162,12 +162,12 @@ time-derivatives and is therefore of the sort of an algebraic constraint that has to hold at each time instant. The main difference to step-21 is that the algebraic constraint there was a mixed Laplace system of the form -@f{eqnarray*} +@f{eqnarray*}{ \mathbf u + {\mathbf K}\lambda \nabla p &=& 0, \\ \nabla\cdot \mathbf u &=& f, @f} where now we have a Stokes system -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& f, \\ \nabla\cdot \mathbf u &=& 0, @f} @@ -192,7 +192,7 @@ solve the Stokes equations for velocity and pressure using the temperature field from the previous time step, which means that we get the velocity for the previous time step. In other words, we first solve the Stokes system for time step $n - 1$ as -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2\eta \varepsilon ({\mathbf u}^{n-1})) + \nabla p^{n-1} &=& -\rho\; \beta \; T^{n-1} \mathbf{g}, \\ @@ -209,7 +209,7 @@ the time derivative $\frac{\partial T}{\partial t}$ by the (one-sided) difference quotient $\frac{\frac 32 T^{n}-2T^{n-1}+\frac 12 T^{n-2}}{k}$ with $k$ the time step size. This gives the discretized-in-time temperature equation -@f{eqnarray*} +@f{eqnarray*}{ \frac 32 T^n - k\nabla \cdot \kappa \nabla T^n @@ -268,7 +268,7 @@ T^n \approx \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}, @f} and above equation is generalized as follows: -@f{eqnarray*} +@f{eqnarray*}{ \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n - k_n\nabla \cdot \kappa \nabla T^n @@ -312,7 +312,7 @@ stable pair $Q_{p+1}^d \times Q_p, p\ge 1$. These are continuous elements, so we can form the weak form of the Stokes equation without problem by integrating by parts and substituting continuous functions by their discrete counterparts: -@f{eqnarray*} +@f{eqnarray*}{ (\nabla {\mathbf v}_h, 2\eta \varepsilon ({\mathbf u}^{n-1}_h)) - (\nabla \cdot {\mathbf v}_h, p^{n-1}_h) @@ -331,7 +331,7 @@ anti-symmetric component of $\nabla {\mathbf v}_h$ plays no role and it leads to the entirely same form if we use the symmetric gradient of $\mathbf v_h$ instead. Consequently, the formulation we consider and that we implement is -@f{eqnarray*} +@f{eqnarray*}{ (\varepsilon({\mathbf v}_h), 2\eta \varepsilon ({\mathbf u}^{n-1}_h)) - (\nabla \cdot {\mathbf v}_h, p^{n-1}_h) @@ -370,7 +370,7 @@ necessary to stabilize the scheme. A better alternative is therefore to add some nonlinear viscosity to the model. Essentially, what this does is to transform the temperature equation from the form -@f{eqnarray*} +@f{eqnarray*}{ \frac{\partial T}{\partial t} + {\mathbf u} \cdot \nabla T @@ -378,7 +378,7 @@ equation from the form \nabla \cdot \kappa \nabla T &=& \gamma @f} to something like -@f{eqnarray*} +@f{eqnarray*}{ \frac{\partial T}{\partial t} + {\mathbf u} \cdot \nabla T @@ -394,7 +394,7 @@ To achieve this, the literature contains a number of approaches. We will here follow one developed by Guermond and Popov that builds on a suitably defined residual and a limiting procedure for the additional viscosity. To this end, let us define a residual $R_\alpha(T)$ as follows: -@f{eqnarray*} +@f{eqnarray*}{ R_\alpha(T) = \left( @@ -411,7 +411,7 @@ within the range $[1,2]$. Note that $R_\alpha(T)$ will be zero if $T$ satisfies the temperature equation, since then the term in parentheses will be zero. Multiplying terms out, we get the following, entirely equivalent form: -@f{eqnarray*} +@f{eqnarray*}{ R_\alpha(T) = \frac 1\alpha @@ -434,7 +434,7 @@ With this residual, we can now define the artificial viscosity as a piecewise constant function defined on each cell $K$ with diameter $h_K$ separately as follows: -@f{eqnarray*} +@f{eqnarray*}{ \nu_\alpha(T)|_K = \beta @@ -521,7 +521,7 @@ discretization error in smooth regions. Using the BDF-2 scheme introduced above, this yields for the simpler case of uniform time steps of size $k$: -@f{eqnarray*} +@f{eqnarray*}{ \frac 32 T^n - k\nabla \cdot \kappa \nabla T^n @@ -563,7 +563,7 @@ the advection works in the section on time stepping. The form for nonuniform time steps that we will have to use in reality is a bit more complicated (which is why we showed the simpler form above first) and reads: -@f{eqnarray*} +@f{eqnarray*}{ \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n - k_n\nabla \cdot \kappa \nabla T^n @@ -602,7 +602,7 @@ form above first) and reads: After settling all these issues, the weak form follows naturally from the strong form shown in the last equation, and we immediately arrive at the weak form of the discretized equations: -@f{eqnarray*} +@f{eqnarray*}{ \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} (\tau_h,T_h^n) + k_n (\nabla \tau_h, \kappa \nabla T_h^n) @@ -641,7 +641,7 @@ $\mathbf{n}\cdot\kappa\nabla T|_{\partial\Omega}=0$. This then results in a matrix equation of form -@f{eqnarray*} +@f{eqnarray*}{ \left( \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} M+k_n A_T\right) T_h^n = F(U_h^{n-1}, U_h^{n-2},T_h^{n-1},T_h^{n-2}), @f} @@ -671,13 +671,13 @@ efficient than the original approach. The best alternative identified there we to use a GMRES solver preconditioned by a block matrix involving the Schur complement. Specifically, the Stokes operator leads to a block structured matrix -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) @f} and as discussed there a good preconditioner is -@f{eqnarray*} +@f{eqnarray*}{ P = \left(\begin{array}{cc} @@ -696,7 +696,7 @@ where $S$ is the Schur complement of the Stokes operator $S=B^TA^{-1}B$. Of course, this preconditioner is not useful because we can't form the various inverses of matrices, but we can use the following as a preconditioner: -@f{eqnarray*} +@f{eqnarray*}{ \tilde P^{-1} = \left(\begin{array}{cc} @@ -742,7 +742,7 @@ a way that first all $x$-components of the velocity are numbered, then the $y$-components, and then the $z$-components, then the matrix $\hat A$ that is associated with this slightly different bilinear form has the form -@f{eqnarray*} +@f{eqnarray*}{ \hat A = \left(\begin{array}{ccc} A_s & 0 & 0 \\ 0 & A_s & 0 \\ 0 & 0 & A_s @@ -752,7 +752,7 @@ where $A_s$ is a Laplace matrix of size equal to the number of shape functions associated with each component of the vector-valued velocity. With this matrix, one could be tempted to define our preconditioner for the velocity matrix $A$ as follows: -@f{eqnarray*} +@f{eqnarray*}{ \tilde A^{-1} = \left(\begin{array}{ccc} \tilde A_s^{-1} & 0 & 0 \\ @@ -856,7 +856,7 @@ program is based. We could of course do the same here. The linear system that we would get would look like this: -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{ccc} A & B^T & 0 \\ B & 0 &0 \\ C & 0 & K \end{array}\right) diff --git a/examples/step-31/doc/results.dox b/examples/step-31/doc/results.dox index e2d8a0307a..a10f5e2488 100644 --- a/examples/step-31/doc/results.dox +++ b/examples/step-31/doc/results.dox @@ -304,7 +304,7 @@ theoretical handle on how to choose in an optimal way. These are: $k\le \min_K \frac{c_kh_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. Here, $c_k$ is dimensionless, but what is the right value?
  • In the computation of the artificial viscosity, -@f{eqnarray*} +@f{eqnarray*}{ \nu_\alpha(T)|_K = \beta @@ -548,7 +548,7 @@ factor of $\frac 78$). As a consequence, a formula that reconciles 2d, 3d, and variable polynomial degree and takes all factors in account reads as follows: -@f{eqnarray*} +@f{eqnarray*}{ k = \frac 1{2 \cdot 1.7} \frac 1{\sqrt{d}} \frac 2d diff --git a/examples/step-32/doc/intro.dox b/examples/step-32/doc/intro.dox index 598378f451..635a49203c 100644 --- a/examples/step-32/doc/intro.dox +++ b/examples/step-32/doc/intro.dox @@ -68,7 +68,7 @@ parallelize things, and finally the actual testcase we will consider. In step-31, we used the following Stokes model for the velocity and pressure field: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& -\rho \; \beta \; T \mathbf{g}, \\ @@ -89,13 +89,13 @@ expression of the form $\rho(T) = \rho_{\text{ref}} $\rho_{\text{ref}}$ at reference temperature and decreases linearly as the temperature increases (as the material expands). The force balance equation then looks properly written like this: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho_{\text{ref}} [1-\beta(T-T_{\text{ref}})] \mathbf{g}. @f} Now note that the gravity force results from a gravity potential as $\mathbf g=-\nabla \varphi$, so that we can re-write this as follows: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& -\rho_{\text{ref}} \; \beta\; T\; \mathbf{g} -\rho_{\text{ref}} [1+\beta T_{\text{ref}}] \nabla\varphi. @@ -104,7 +104,7 @@ The second term on the right is time independent, and so we could introduce a new "dynamic" pressure $p_{\text{dyn}}=p+\rho_{\text{ref}} [1+\beta T_{\text{ref}}] \varphi=p_{\text{total}}-p_{\text{static}}$ with which the Stokes equations would read: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p_{\text{dyn}} &=& -\rho_{\text{ref}} \; \beta \; T \; \mathbf{g}, \\ @@ -119,7 +119,7 @@ on the velocity field.) On the other hand, we will here use the form of the Stokes equations that considers the total pressure instead: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho(T)\; \mathbf{g}, \\ @@ -176,7 +176,7 @@ this turns out not to be a limiting factor.

    The scaling of discretized equations

    Remember that we want to solve the following set of equations: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho(T) \mathbf{g}, \\ @@ -195,7 +195,7 @@ the temperature equation forward by one time interval. The problem under consideration in this current section is with the Stokes problem: if we discretize it as usual, we get a linear system -@f{eqnarray*} +@f{eqnarray*}{ M \; X = \left(\begin{array}{cc} @@ -272,7 +272,7 @@ in our domain (which experiments show is best chosen to be the diameter of plumes — around 10 km — rather than the diameter of the domain). Using these %numbers for $\eta$ and $L$, this factor is around $10^{17}$. So, we now get this for the Stokes system: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho(T) \; \mathbf{g}, \\ @@ -282,7 +282,7 @@ The trouble with this is that the result is not symmetric any more (we have $\frac{\eta}{L} \nabla \cdot$ at the bottom left, but not its transpose operator at the top right). This, however, can be cured by introducing a scaled pressure $\hat p = \frac{L}{\eta}p$, and we get the scaled equations -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla \left(\frac{\eta}{L} \hat p\right) &=& \rho(T) \; \mathbf{g}, @@ -310,7 +310,7 @@ also have to scale the pressure immediately before solving. In this tutorial program, we apply a variant of the preconditioner used in step-31. That preconditioner was built to operate on the system matrix $M$ in block form such that the product matrix -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} M = \left(\begin{array}{cc} @@ -430,7 +430,7 @@ Similarly to step-31, we will use an artificial viscosity for stabilization based on a residual of the equation. As a difference to step-31, we will provide two slightly different definitions of the stabilization parameter. For $\alpha=1$, we use the same definition as in step-31: -@f{eqnarray*} +@f{eqnarray*}{ \nu_\alpha(T)|_K = \nu_1(T)|_K @@ -449,16 +449,16 @@ regions where the residual is large (around steep gradients). This definition has been shown to work well for the given case, $\alpha = 1$ in step-31, but it is usually less effective as the diffusion for $\alpha=2$. For that case, we choose a slightly more readable definition of the viscosity, -@f{eqnarray*} +@f{eqnarray*}{ \nu_2(T)|_K = \min (\nu_h^\mathrm{max}|_K,\nu_h^\mathrm{E}|_K) @f} where the first term gives again the maximum dissipation (similarly to a first order upwind scheme), -@f{eqnarray*} +@f{eqnarray*}{ \nu^\mathrm{max}_h|_K = \beta h_K \|\mathbf {u}\|_{L^\infty(K)} @f} and the entropy viscosity is defined as -@f{eqnarray*} +@f{eqnarray*}{ \nu^\mathrm{E}_h|_K = c_R \frac{h_K^2 \|R_\mathrm{2,E}(T)\|_{L^\infty(K)}} {\|E(T) - \bar{E}(T)\|_{L^\infty(\Omega)} }. @f} @@ -469,7 +469,7 @@ Comput. Phys., 230, 4248--4267. Compared to the case $\alpha = 1$, the residual is computed from the temperature entropy, $E(T) = \frac12 (T-T_m)^2$ with $T_m$ an average temperature (we choose the mean between the maximum and minimum temperature in the computation), which gives the following formula -@f{eqnarray*} +@f{eqnarray*}{ R_\mathrm{E}(T) = \frac{\partial E(T)}{\partial t} + (T-T_\mathrm{m}) \left(\mathbf{u} \cdot \nabla T - \kappa \nabla^2 T - \gamma\right). @f} @@ -504,7 +504,7 @@ decreasing $c_R$. As we did in step-31, we first determine an optimal value of $\beta$ by using the "traditional" formula -@f{eqnarray*} +@f{eqnarray*}{ \nu_\alpha(T)|_K = \beta @@ -564,7 +564,7 @@ The standard Taylor-Hood discretization for Stokes, using the $Q_{k+1}^d the divergence equation reads $(q_h, \textrm{div}\; \mathbf u_h)=0, \forall q_h\in Q_h$. Because the pressure space does contain the function $q_h=1$, we get -@f{align*} +@f{align*}{ 0 = (1, \textrm{div}\; \mathbf u_h)_\Omega = \int_\Omega \textrm{div}\; \mathbf u_h = \int_{\partial\Omega} \mathbf n \cdot \mathbf u_h @@ -589,7 +589,7 @@ This space turns out to be stable for the Stokes equation. Because the space is discontinuous, we can now in particular choose the test function $q_h(\mathbf x)=\chi_K(\mathbf x)$, i.e. the characteristic function of cell $K$. We then get in a similar fashion as above -@f{align*} +@f{align*}{ 0 = (q_h, \textrm{div}\; \mathbf u_h)_\Omega = (1, \textrm{div}\; \mathbf u_h)_K @@ -851,7 +851,7 @@ realistic description is a goal of the ASPECT code in development at the time of writing this. As a reminder, let us again state the equations we want to solve are these: -@f{eqnarray*} +@f{eqnarray*}{ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla \left( \frac{\eta}{L} \hat p\right) &=& \rho(T) \mathbf{g}, @@ -900,14 +900,14 @@ the following quantities: by looking at the pictures shown in the results section below. The initial temperature field we use here is given in terms of the radius by - @f{align*} + @f{align*}{ s &= \frac{\|\mathbf x\|-R_0}{R_1-R_0}, \\ \varphi &= \arctan \frac{y}{x}, \\ \tau &= s + \frac 15 s(1-s) \sin(6\varphi) q(z), \\ T(\mathbf x) &= T_0(1-\tau) + T_1\tau, @f} where - @f{align*} + @f{align*}{ q(z) = \left\{ \begin{array}{ll} 1 & \text{in 2d} \\ @@ -986,7 +986,7 @@ the following quantities: the medium on the other side is in motion as well, so the shear stress would, in the simplest case, be proportional to the velocity difference, leading to a boundary condition of the form - @f{align*} + @f{align*}{ \mathbf{n}\cdot [2\eta \varepsilon(\mathbf v)] &= s \mathbf{n} \times [\mathbf v - \mathbf v_0], diff --git a/examples/step-32/doc/results.dox b/examples/step-32/doc/results.dox index 93b1b06cee..a989c9cd2e 100644 --- a/examples/step-32/doc/results.dox +++ b/examples/step-32/doc/results.dox @@ -306,14 +306,14 @@ address to make the program more useful: that rises up cools adiabatically, and cold material that sinks down heats adiabatically. The correct temperature equation would therefore look somewhat like this: - @f{eqnarray*} + @f{eqnarray*}{ \frac{D T}{Dt} - \nabla \cdot \kappa \nabla T &=& \gamma + \tau\frac{Dp}{Dt}, @f} or, expanding the advected derivative $\frac{D}{Dt} = \frac{\partial}{\partial t} + \mathbf u \cdot \nabla$: - @f{eqnarray*} + @f{eqnarray*}{ \frac{\partial T}{\partial t} + {\mathbf u} \cdot \nabla T diff --git a/examples/step-32/step-32.cc b/examples/step-32/step-32.cc index e6f3c5af77..04f45cd483 100644 --- a/examples/step-32/step-32.cc +++ b/examples/step-32/step-32.cc @@ -213,7 +213,7 @@ namespace Step32 // introduction, this preconditioner differs in a number of key portions // from the one used in step-31. Specifically, it is a right preconditioner, // implementing the matrix - // @f{align*} + // @f{align*}{ // \left(\begin{array}{cc}A^{-1} & B^T // \\0 & S^{-1} // \end{array}\right) diff --git a/examples/step-33/doc/intro.dox b/examples/step-33/doc/intro.dox index 4f5b09ccec..89ee0634c8 100644 --- a/examples/step-33/doc/intro.dox +++ b/examples/step-33/doc/intro.dox @@ -52,7 +52,7 @@ G_i(\mathbf w)$, $i=1,\ldots,dim+2$. For the Euler equations, the flux matrix $\mathbf F$ (or system of flux functions) is defined as (shown here for the case $d=3$) -@f{eqnarray*} +@f{eqnarray*}{ \mathbf F(\mathbf w) = \left( @@ -67,7 +67,7 @@ is defined as (shown here for the case $d=3$) @f} and we will choose as particular right hand side forcing only the effects of gravity, described by -@f{eqnarray*} +@f{eqnarray*}{ \mathbf G(\mathbf w) = \left( @@ -82,7 +82,7 @@ gravity, described by @f} where $\mathbf g=(g_1,g_2,g_3)^T$ denotes the gravity vector. With this, the entire system of equations reads: -@f{eqnarray*} +@f{eqnarray*}{ \partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s + \delta_{is} p)}{\partial x_s} &=& g_i \rho, \qquad i=1,\dots,d, \\ \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\ @@ -108,7 +108,7 @@ in step-12: We choose a finite element space $V_h$, and integrate our conservation law against our (vector-valued) test function $\mathbf{z} \in V_h$. We then integrate by parts and approximate the boundary flux with a numerical flux $\mathbf{H}$, -@f{eqnarray*} +@f{eqnarray*}{ &&\int_{\Omega} (\partial_t \mathbf{w}, \mathbf{z}) + (\nabla \cdot \mathbf{F}(\mathbf{w}), \mathbf{z}) \\ &\approx &\int_{\Omega} (\partial_t \mathbf{w}, \mathbf{z}) - (\mathbf{F}(\mathbf{w}), \nabla \mathbf{z}) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}) + \int_{\partial \Omega} (\mathbf{H}(\mathbf{w}^+, \mathbf{w}^-, \mathbf{n}), \mathbf{z}^+), @f} @@ -136,7 +136,7 @@ Compressible Euler Equations", PhD thesis, University of Heidelberg, 2002). We use a time stepping scheme to substitute the time derivative in the above equations. For simplicity, we define $ \mathbf{B}({\mathbf{w}_{n}})(\mathbf z) $ as the spatial residual at time step $n$ : -@f{eqnarray*} +@f{eqnarray*}{ \mathbf{B}(\mathbf{w}_{n})(\mathbf z) &=& - \int_{\Omega} \left(\mathbf{F}(\mathbf{w}_n), \nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w}_n , \nabla \mathbf{z}) \\ @@ -151,7 +151,7 @@ above equations. For simplicity, we define $ \mathbf{B}({\mathbf{w}_{n}})(\mathb At each time step, our full discretization is thus that the residual applied to any test function $\mathbf z$ equals zero: -@f{eqnarray*} +@f{eqnarray*}{ R(\mathbf{W}_{n+1})(\mathbf z) &=& \int_{\Omega} \left(\frac{{\mathbf w}_{n+1} - \mathbf{w}_n}{\delta t}, \mathbf{z}\right)+ @@ -176,7 +176,7 @@ applied, and $\delta T$ the current time step. With these choices, equating the residual to zero results in a nonlinear system of equations $R(\mathbf{W}_{n+1})=0$. We solve this nonlinear system by a Newton iteration (in the same way as explained in step-15), i.e. by iterating -@f{eqnarray*} +@f{eqnarray*}{ R'(\mathbf{W}^k_{n+1},\delta \mathbf{W}_{n+1}^k)(\mathbf z) & = & - R(\mathbf{W}^{k}_{n+1})(\mathbf z) \qquad \qquad \forall \mathbf z\in V_h \\ \mathbf{W}^{k+1}_{n+1} &=& \mathbf{W}^k_{n+1} + \delta \mathbf{W}^k_{n+1}, @@ -184,7 +184,7 @@ R(\mathbf{W}^{k}_{n+1})(\mathbf z) \qquad \qquad \forall \mathbf z\in V_h \\ until $|R(\mathbf{W}^k_{n+1})|$ (the residual) is sufficiently small. By testing with the nodal basis of a finite element space instead of all $\mathbf z$, we arrive at a linear system for $\delta \mathbf W$: -@f{eqnarray*} +@f{eqnarray*}{ \mathbf R'(\mathbf{W}^k_{n+1})\delta \mathbf{W}^k_{n+1} & = & - \mathbf R(\mathbf{W}^{k}_{n+1}). @f} diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 15ac9725dd..d81985ff72 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -1809,7 +1809,7 @@ namespace Step33 // vector_valued module). It will be represented by the variable // component_i below. With this, the residual term can be // re-written as - // @f{eqnarray*} + // @f{eqnarray*}{ // R_i &=& // \left(\frac{(\mathbf{w}_{n+1} - // \mathbf{w}_n)_{\text{component\_i}}}{\delta diff --git a/examples/step-34/doc/intro.dox b/examples/step-34/doc/intro.dox index 704b9b4bcf..304af2b7a2 100644 --- a/examples/step-34/doc/intro.dox +++ b/examples/step-34/doc/intro.dox @@ -277,7 +277,7 @@ arbitrary constant $c$, which means that, if we set the Neumann data to be zero, then any constant $\phi = \phi_\infty$ will be a solution. Inserting the constant solution and the Neumann boundary condition in the boundary integral equation, we have -@f{align*} +@f{align*}{ \alpha\left(\mathbf{x}\right)\phi\left(\mathbf{x}\right) &=\int_{\Omega}\phi\left(\mathbf{y}\right)\delta\left(\mathbf{y}-\mathbf{x}\right)\, dy\\ \Rightarrow @@ -323,7 +323,7 @@ v}+\mathbf{v}_\infty=\nabla\phi+\mathbf{v}_\infty$ for some (unknown) pressure $p$ and a given constant $\rho$. In other words, we would like to verify that Bernoulli's law as stated above indeed holds. To show this, we use that the left hand side of this equation equates to -@f{align*} +@f{align*}{ \mathbf{v}\cdot\nabla\mathbf{v} &= [(\nabla\phi+\mathbf{v}_\infty)\cdot\nabla] (\nabla\phi+\mathbf{v}_\infty) @@ -336,7 +336,7 @@ write this expression as the gradient of something (remember that $\rho$ is a constant). The next step is more convenient if we consider the components of the equation individually (summation over indices that appear twice is implied): -@f{align*} +@f{align*}{ [\mathbf{v}\cdot\nabla\mathbf{v}]_i &= (\partial_j\phi+v_{\infty,j}) \partial_j \partial_i\phi @@ -351,7 +351,7 @@ convenient if we consider the components of the equation individually @f} because $\partial_j \partial_j\phi = \Delta \phi = 0$ and $\textrm{div} \ \mathbf{v}_\infty=0$. Next, -@f{align*} +@f{align*}{ [\mathbf{v}\cdot\nabla\mathbf{v}]_i &= \partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi] @@ -384,7 +384,7 @@ because $\partial_j \partial_j\phi = \Delta \phi = 0$ and $\textrm{div} @f} Again, the last term disappears because $\mathbf{v}_\infty$ is constant and we can merge the first and third term into one: -@f{align*} +@f{align*}{ [\mathbf{v}\cdot\nabla\mathbf{v}]_i &= \partial_i (\partial_j [(\partial_j\phi) \phi + v_{\infty,j} \partial_j \phi]) @@ -399,7 +399,7 @@ can merge the first and third term into one: We now only need to massage that last term a bit more. Using the product rule, we get -@f{align*} +@f{align*}{ \partial_j [\partial_i (\partial_j\phi) \phi] &= \partial_i [\partial_j \partial_j\phi] \phi @@ -410,7 +410,7 @@ The first of these terms is zero (because, again, the summation over $j$ gives $\Delta\phi$, which is zero). The last term can be written as $\frac 12 \partial_i [(\partial_j\phi)(\partial_j\phi)]$ which is in the desired gradient form. As a consequence, we can now finally state that -@f{align*} +@f{align*}{ [\mathbf{v}\cdot\nabla\mathbf{v}]_i &= \partial_i (\partial_j [(\partial_j\phi) \phi + v_{\infty,j} \partial_j \phi]) @@ -676,7 +676,7 @@ that respects the continuous geometry behind the discrete initial mesh. For a sphere of radius $a$ translating at a velocity of $U$ in the $x$ direction, the potential reads -@f{align*} +@f{align*}{ \phi = -\frac{1}{2}U \left(\frac{a}{r}\right)3 r \cos\theta @f} see, e.g. J. N. Newman, Marine Hydrodynamics, 1977, diff --git a/examples/step-35/doc/intro.dox b/examples/step-35/doc/intro.dox index 208960e8f7..49f7e80490 100644 --- a/examples/step-35/doc/intro.dox +++ b/examples/step-35/doc/intro.dox @@ -12,7 +12,7 @@ University. Most of the work for this program is by him.

    Motivation

    The purpose of this program is to show how to effectively solve the incompressible time-dependent Navier-Stokes equations. These equations describe the flow of a viscous incompressible fluid and read -@f{align*} +@f{align*}{ u_t + u \cdot \nabla u - \nu \Delta u + \nabla p = f, \\ \nabla \cdot u = 0, @f} @@ -42,7 +42,7 @@ In previous tutorial programs (see for instance step-20 and step-22) we have seen how to solve the time-independent Stokes equations using a Schur complement approach. For the time-dependent case, after time discretization, we would arrive at a system like -@f{align*} +@f{align*}{ \frac1\tau u^k - \nu \Delta u^k + \nabla p^k = F^k, \\ \nabla \cdot u^k = 0, @f} diff --git a/examples/step-36/doc/intro.dox b/examples/step-36/doc/intro.dox index 92816cefbe..47d66f0750 100644 --- a/examples/step-36/doc/intro.dox +++ b/examples/step-36/doc/intro.dox @@ -17,7 +17,7 @@ here, though the general techniques outlined in this program are of course equally applicable to the other applications above. Eigenspectrum problems have the general form -@f{align*} +@f{align*}{ L \Psi &= \varepsilon \Psi \qquad &&\text{in}\ \Omega, \\ \Psi &= 0 &&\text{on}\ \partial\Omega, @f} @@ -124,7 +124,7 @@ If you assume for a moment that we had renumbered degrees of freedom in such a way that all of those on the Dirichlet boundary come last, then the linear system we would get when solving a regular PDE with a right hand side would look like this: -@f{align*} +@f{align*}{ \begin{pmatrix} A_i & 0 \\ 0 & D_b \end{pmatrix} @@ -162,7 +162,7 @@ There, eliminating boundary values affects both matrices $A$ and $M$ that we will solve with in the current tutorial program. After elimination of boundary values, we then receive an eigenvalue problem that can be partitioned like this: -@f{align*} +@f{align*}{ \begin{pmatrix} A_i & 0 \\ 0 & D_A \end{pmatrix} diff --git a/examples/step-37/doc/intro.dox b/examples/step-37/doc/intro.dox index e50ea2c32b..1d9859bac6 100644 --- a/examples/step-37/doc/intro.dox +++ b/examples/step-37/doc/intro.dox @@ -64,7 +64,7 @@ solution. In order to find out how we can write a code that performs a matrix-vector product, but does not need to store the matrix elements, let us start at looking how a finite element matrix A is assembled: -@f{eqnarray*} +@f{eqnarray*}{ A = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_{\mathrm{cell,{loc-glob}}}^T A_{\mathrm{cell}} P_{\mathrm{cell,{loc-glob}}}. @f} @@ -76,7 +76,7 @@ variable and is used in the assembly calls filling matrices in deal.II. Here, Acell denotes the cell matrix associated with A. If we are to perform a matrix-vector product, we can hence use that -@f{eqnarray*} +@f{eqnarray*}{ y &=& A\cdot u = \left(\sum_{\text{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T A_\mathrm{cell} P_\mathrm{cell,{loc-glob}}\right) \cdot u \\ @@ -162,7 +162,7 @@ is much cheaper. One way to improve this is to realize that conceptually the local matrix can be thought of as the product of three matrices, -@f{eqnarray*} +@f{eqnarray*}{ A_\mathrm{cell} = B_\mathrm{cell}^T D_\mathrm{cell} B_\mathrm{cell}, @f} where for the example of the Laplace operator the (q*dim+d,i)-th @@ -175,7 +175,7 @@ dim copies of each of these values). This kind of representation of finite element matrices can often be found in the engineering literature. When the cell matrix is applied to a vector, -@f{eqnarray*} +@f{eqnarray*}{ A_\mathrm{cell}\cdot u_\mathrm{cell} = B_\mathrm{cell}^T D_\mathrm{cell} B_\mathrm{cell} \cdot u_\mathrm{cell}, @f} @@ -240,14 +240,14 @@ problem. Each diagonal block is the Jacobian transformation that goes from the reference cell to the real cell. Putting things together, we find that -@f{eqnarray*} +@f{eqnarray*}{ A_\mathrm{cell} = B_\mathrm{cell}^T D B_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^{-1} D_\mathrm{cell} J_\mathrm{cell}^{-\mathrm T} B_\mathrm{ref\_cell}, @f} so we calculate the product (starting the local product from the right) -@f{eqnarray*} +@f{eqnarray*}{ v_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^{-1} D J_\mathrm{cell}^{-\mathrm T} B_\mathrm{ref\_cell} u_\mathrm{cell}, \quad v = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T diff --git a/examples/step-37/doc/results.dox b/examples/step-37/doc/results.dox index 776ca401cb..d09c6ffe97 100644 --- a/examples/step-37/doc/results.dox +++ b/examples/step-37/doc/results.dox @@ -457,7 +457,7 @@ arise in the mathematical formulation and how they are implemented in a matrix-based variant. In essence, an inhomogeneous Dirichlet condition sets some of the nodal values in the solution to given values rather than determining them through the variational principles, -@f{eqnarray*} +@f{eqnarray*}{ u_h(\mathbf{x}) = \sum_{i\in \mathcal N} \varphi_i(\mathbf{x}) u_i = \sum_{i\in \mathcal N \setminus \mathcal N_D} \varphi_i(\mathbf{x}) u_i + \sum_{i\in \mathcal N_D} \varphi_i(\mathbf{x}) g_i, @@ -470,7 +470,7 @@ of boundary values on the Dirichlet-constrained node points $i\in \mathcal N_D$. We then insert this solution representation into the weak form, e.g. the Laplacian shown above, and move the known quantities to the right hand side: -@f{eqnarray*} +@f{eqnarray*}{ (\nabla \varphi_i, \nabla u_h)_\Omega &=& (\varphi_i, f)_\Omega \quad \Rightarrow \\ \sum_{j\in \mathcal N \setminus \mathcal N_D}(\nabla \varphi_i,\nabla \varphi_j)_\Omega \, u_j &=& (\varphi_i, f)_\Omega diff --git a/examples/step-38/doc/intro.dox b/examples/step-38/doc/intro.dox index ecbf3ab1b9..096ad856b7 100644 --- a/examples/step-38/doc/intro.dox +++ b/examples/step-38/doc/intro.dox @@ -18,7 +18,7 @@ In this example, we show how to solve a partial differential equation (PDE) on a codimension one surface $\Gamma \subset \mathbb R^3$ made of quadrilaterals, i.e. on a surface in 3d or a line in 2d. We focus on the following elliptic second order PDE -@f{align*} +@f{align*}{ -\Delta_\Gamma u &= f \qquad \text{on } \qquad \Gamma,\\ u &= g \qquad \text{on} \qquad \partial \Gamma, @f} @@ -78,7 +78,7 @@ and take advantage of the partition ${\mathbb T}$ to further write Moreover, each integral in the above expression is computed in the reference element $\hat K \dealcoloneq [0,1]^2$ so that -@f{align*} +@f{align*}{ \int_{K} \nabla_{K} u \cdot \nabla_{K} v &= \int_{\hat K} \nabla (u \circ \mathbf x_K)^T G_K^{-1} (D \mathbf @@ -165,7 +165,7 @@ We produce one test case for a 2d problem and another one for 3d: Taking into account that the transformation is length preserving, i.e. a segment of length $dt$ is mapped onto a piece of curve of exactly the same length, the tangential Laplacian then satisfies - @f{align*} + @f{align*}{ \Delta_\Gamma u &= \frac{d^2}{dt^2}(-2\cos t \sin t) = -2 \frac{d}{dt}(-\sin^2 t + \cos^2 t) diff --git a/examples/step-39/doc/intro.dox b/examples/step-39/doc/intro.dox index 924519948c..5a8d582931 100644 --- a/examples/step-39/doc/intro.dox +++ b/examples/step-39/doc/intro.dox @@ -25,7 +25,7 @@ averaging operator turns into a jump. The interior penalty method for the proble -\Delta u = f \text{ in }\Omega \qquad u = u^D \text{ on } \partial\Omega @f] becomes -@f{multline*} +@f{multline*}{ \sum_{K\in \mathbb T_h} (\nabla u, \nabla v)_K \\ + \sum_{F \in F_h^i} \biggl\{4\sigma_F (\average{ u \mathbf n}, \average{ v \mathbf n })_F @@ -64,7 +64,7 @@ right. As we will see below, even the error estimate is of the same structure, since it can be written as -@f{align*} +@f{align*}{ \eta^2 &= \eta_K^2 + \eta_F^2 + \eta_B^2 \\ \eta_K^2 &= \sum_{K\in \mathbb T_h} h^2 \|f + \Delta u_h\|^2 diff --git a/examples/step-40/doc/intro.dox b/examples/step-40/doc/intro.dox index c26d797eb2..cbd60bafd4 100644 --- a/examples/step-40/doc/intro.dox +++ b/examples/step-40/doc/intro.dox @@ -97,7 +97,7 @@ things work internally in this program. This program essentially re-solves what we already do in step-6, i.e. it solves the Laplace equation -@f{align*} +@f{align*}{ -\Delta u &= f \qquad &&\text{in}\ \Omega=[0,1]^2, \\ u &= 0 \qquad &&\text{on}\ \partial\Omega. @f} @@ -109,7 +109,7 @@ not to do something useful but to show how useful programs can be implemented using deal.II. Be that as it may, to make things at least a tiny bit interesting, we choose the right hand side as a discontinuous function, -@f{align*} +@f{align*}{ f(x,y) = \left\{ diff --git a/examples/step-41/doc/intro.dox b/examples/step-41/doc/intro.dox index 04eb08e134..3470ad73e3 100644 --- a/examples/step-41/doc/intro.dox +++ b/examples/step-41/doc/intro.dox @@ -49,7 +49,7 @@ is active there.

    Classical formulation

    The classical formulation of the problem possesses the following form: -@f{align*} +@f{align*}{ -\textrm{div}\ \sigma &\geq f & &\quad\text{in } \Omega,\\ \sigma &= \nabla u & &\quad\text{in } \Omega,\\ u(\mathbf x) &= 0 & &\quad\text{on }\partial\Omega,\\ @@ -91,15 +91,15 @@ obstacle).

    Derivation of the variational inequality

    An obvious way to obtain the variational formulation of the obstacle problem is to consider the total potential energy: -@f{equation*} +@f{equation*}{ E(u) \dealcoloneq \dfrac{1}{2}\int\limits_{\Omega} \nabla u \cdot \nabla u - \int\limits_{\Omega} fu. @f} We have to find a solution $u\in G$ of the following minimization problem: -@f{equation*} +@f{equation*}{ E(u)\leq E(v)\quad \forall v\in G, @f} with the convex set of admissible displacements: -@f{equation*} +@f{equation*}{ G \dealcoloneq \lbrace v\in V: v\geq g \text{ a.e. in } \Omega\rbrace,\quad V\dealcoloneq H^1_0(\Omega). @f} This set takes care of the third and fifth conditions above (the boundary @@ -107,7 +107,7 @@ values and the complementarity condition). Consider now the minimizer $u\in G$ of $E$ and any other function $v\in G$. Then the function -@f{equation*} +@f{equation*}{ F(\varepsilon) \dealcoloneq E(u+\varepsilon(v-u)),\quad\varepsilon\in\left[0,1\right], @f} takes its minimum at $\varepsilon = 0$ (because $u$ is a minimizer of the @@ -118,7 +118,7 @@ convexity of $G$. If we compute $F'(\varepsilon)\vert_{\varepsilon=0}$ it yields the variational formulation we are searching for: Find a function $u\in G$ with -@f{equation*} +@f{equation*}{ \left(\nabla u, \nabla(v-u)\right) \geq \left(f,v-u\right) \quad \forall v\in G. @f} @@ -151,12 +151,12 @@ besides regular functions as well. This yields: Find $u\in V$ and $\lambda\in K$ such that -@f{align*} +@f{align*}{ a(u,v) + b(v,\lambda) &= f(v),\quad &&v\in V\\ b(u,\mu - \lambda) &\leq \langle g,\mu - \lambda\rangle,\quad&&\mu\in K, @f} with -@f{align*} +@f{align*}{ a(u,v) &\dealcoloneq \left(\nabla u, \nabla v\right),\quad &&u,v\in V\\ b(u,\mu) &\dealcoloneq \langle u,\mu\rangle,\quad &&u\in V,\quad\mu\in V'. @f} @@ -180,7 +180,7 @@ inequality constraints. To get there, let us assume that we discretize both $u$ and $\lambda$ with the same finite element space, for example the usual $Q_k$ spaces. We would then get the equations -@f{eqnarray*} +@f{eqnarray*}{ &A U + B\Lambda = F,&\\ &[BU-G]_i \geq 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[BU-G]_i = 0 \qquad \forall i.& @@ -194,27 +194,27 @@ assembling all terms that yield this mass matrix, namely a quadrature formula where quadrature points are only located at the interpolation points at which shape functions are defined; since all but one shape function are zero at these locations, we get a diagonal mass matrix with -@f{align*} +@f{align*}{ B_{ii} = \int_\Omega \varphi_i(\mathbf x)^2\ \textrm{d}x, \qquad B_{ij}=0 \ \text{for } i\neq j. @f} To define $G$ we use the same technique as for $B$. In other words, we define -@f{align*} +@f{align*}{ G_{i} = \int_\Omega g_h(x) \varphi_i(\mathbf x)\ \textrm{d}x, @f} where $g_h$ is a suitable approximation of $g$. The integral in the definition of $B_{ii}$ and $G_i$ are then approximated by the trapezoidal rule. With this, the equations above can be restated as -@f{eqnarray*} +@f{eqnarray*}{ &A U + B\Lambda = F,&\\ &U_i-B_{ii}^{-1}G_i \ge 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[U_i-B_{ii}^{-1}G_i] = 0 \qquad \forall i\in{\cal S}.& @f} Now we define for each degree of freedom $i$ the function -@f{equation*} +@f{equation*}{ C([BU]_i,\Lambda_i) \dealcoloneq -\Lambda_i + \min\lbrace 0, \Lambda_i + c([BU]_i - G_i) \rbrace, @f} with some $c>0$. (In this program we choose $c = 100$. It is a kind of a @@ -224,7 +224,7 @@ current program if we use 7 global refinements.) After some head-scratching one can then convince oneself that the inequalities above can equivalently be rewritten as -@f{equation*} +@f{equation*}{ C([BU]_i,\Lambda_i) = 0, \qquad \forall i\in{\cal S}. @f} The primal-dual active set strategy we will use here is an iterative scheme which is based on @@ -243,7 +243,7 @@ The algorithm for the primal-dual active set method works as follows (NOTE: $B = $\mathcal{S}=\mathcal{A}_k\cup\mathcal{F}_k$ and $\mathcal{A}_k\cap\mathcal{F}_k=\emptyset$ and set $k=1$. 2. Find the primal-dual pair $(U^k,\Lambda^k)$ that satisfies - @f{align*} + @f{align*}{ AU^k + B\Lambda^k &= F,\\ [BU^k]_i &= G_i\quad&&\forall i\in\mathcal{A}_k,\\ \Lambda_i^k &= 0\quad&&\forall i\in\mathcal{F}_k. @@ -252,7 +252,7 @@ The algorithm for the primal-dual active set method works as follows (NOTE: $B = are fixed, with the first condition yielding the remaining $|S|$ equations necessary to determine both $U$ and $\Lambda$. 3. Define the new active and inactive sets by - @f{equation*} + @f{equation*}{ \begin{split} \mathcal{A}_{k+1} \dealcoloneq \lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)< 0\rbrace,\\ \mathcal{F}_{k+1} \dealcoloneq \lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)\geq 0\rbrace. @@ -280,7 +280,7 @@ condition, we can distinguish the following cases: Second, the method above appears intuitively correct and useful but a bit ad hoc. However, it can be derived in a concisely in the following way. To this end, note that we'd like to solve the nonlinear system -@f{eqnarray*} +@f{eqnarray*}{ &A U + B\Lambda = F,&\\ &C([BU-G]_i, \Lambda_i) = 0, \qquad \forall i.& @@ -289,20 +289,20 @@ We can iteratively solve this by always linearizing around the previous iterate (i.e., applying a Newton method), but for this we need to linearize the function $C(\cdot,\cdot)$ that is not differentiable. That said, it is slantly differentiable, and in fact we have -@f{equation*} +@f{equation*}{ \dfrac{\partial}{\partial U^k_i}C([BU^k]_i,\Lambda^k_i) = \begin{cases} cB_{ii},& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)< 0\\ 0,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)\geq 0. \end{cases} @f} -@f{equation*} +@f{equation*}{ \dfrac{\partial}{\partial\Lambda^k_i}C([BU^k]_i,\Lambda^k_i) = \begin{cases} 0,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)< 0\\ -1,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)\geq 0. \end{cases} @f} This suggest a semismooth Newton step of the form -@f{equation*} +@f{equation*}{ \begin{pmatrix} A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & B_{\mathcal{F}_k} & 0\\ A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & 0 & B_{\mathcal{A}_k}\\ @@ -325,7 +325,7 @@ Rather than solving for updates $\delta U, \delta \Lambda$, we can also solve for the variables we are interested in right away by setting $\delta U^k \dealcoloneq U^{k+1} - U^k$ and $\delta \Lambda^k \dealcoloneq \Lambda^{k+1} - \Lambda^k$ and bringing all known terms to the right hand side. This yields -@f{equation*} +@f{equation*}{ \begin{pmatrix} A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & B_{\mathcal{F}_k} & 0\\ A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & 0 & B_{\mathcal{A}_k}\\ @@ -345,7 +345,7 @@ These are the equations outlined above in the description of the basic algorithm We could even drive this a bit further. It's easy to see that we can eliminate the third row and the third column because it implies $\Lambda_{\mathcal{F}_k} = 0$: -@f{equation*} +@f{equation*}{ \begin{pmatrix} A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & 0\\ A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & B_{\mathcal{A}_k}\\ @@ -362,13 +362,13 @@ because it implies $\Lambda_{\mathcal{F}_k} = 0$: This shows that one in fact only needs to solve for the Lagrange multipliers located on the active set. By considering the second row one would then recover the full Lagrange multiplier vector through -@f{equation*} +@f{equation*}{ \Lambda^k_S = B^{-1}\left(f_{\mathcal{S}} - A_{\mathcal{S}}U^k_{\mathcal{S}}\right). @f} Because of the third row and the fact that $B_{\mathcal{A}_k}$ is a diagonal matrix we are able to calculate $U^k_{\mathcal{A}_k}=B^{-1}_{\mathcal{A}_k}G_{\mathcal{A}_k}$ directly. We can therefore also write the linear system as follows: -@f{equation*} +@f{equation*}{ \begin{pmatrix} A_{\mathcal{F}_k\mathcal{F}_k} & 0\\ 0 & Id_{\mathcal{A}_k} \\ @@ -385,7 +385,7 @@ linear system as follows: @f} Fortunately, this form is easy to arrive at: we simply build the usual Laplace linear system -@f{equation*} +@f{equation*}{ \begin{pmatrix} A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} \\ A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} diff --git a/examples/step-42/doc/intro.dox b/examples/step-42/doc/intro.dox index 527516706d..d09e30d270 100644 --- a/examples/step-42/doc/intro.dox +++ b/examples/step-42/doc/intro.dox @@ -53,7 +53,7 @@ step-32.

    Classical formulation

    The classical formulation of the problem possesses the following form: -@f{align*} +@f{align*}{ \varepsilon(\mathbf u) &= A\sigma + \varepsilon^p & &\quad\text{in } \Omega,\\ -\textrm{div}\ \sigma &= \mathbf f & &\quad\text{in } \Omega,\\ \varepsilon^p:(\tau - \sigma) &\geq 0\quad\forall\tau\text{ with @@ -141,11 +141,11 @@ hint at how difficult it is to construct good solvers and preconditioners. With this said, let us simply state the problem we obtain after reformulation (again, details can be found in the paper): Find a displacement $\mathbf u \in V^+$ so that -@f{align*} +@f{align*}{ \left(P_{\Pi}(C\varepsilon(\mathbf u)),\varepsilon(\varphi) - \varepsilon(\mathbf u)\right) \geq 0,\quad \forall \varphi\in V^+. @f} where the projector $P_\Pi$ is defined as -@f{align*} +@f{align*}{ P_{\Pi}(\tau) \dealcoloneq \begin{cases} \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\ \left[ @@ -158,7 +158,7 @@ where the projector $P_\Pi$ is defined as @f} and the space $V^+$ is the space of all displacements that satisfy the contact condition: -@f{align*} +@f{align*}{ V &= \left\{ \mathbf u\in \left[H^1(\Omega)\right]^{d}: @@ -196,7 +196,7 @@ a semi-smooth function with an appropriately chosen "derivative". In the current case, we will run our iteration by solving in each iteration $i$ the following equation (still an inequality, but linearized): -@f{align*} +@f{align*}{ \label{eq:linearization} \left(I_{\Pi}\varepsilon(\tilde {\mathbf u}^{i}), \varepsilon(\varphi) - \varepsilon(\tilde {\mathbf u}^{i})\right) \geq @@ -207,7 +207,7 @@ the following equation (still an inequality, but linearized): \quad \forall \varphi\in V^+, @f} where the rank-4 tensor $I_\Pi=I_\Pi(\varepsilon^D(\mathbf u^{i-1}))$ given by -@f{align} +@f{align}{ I_\Pi = \begin{cases} C_{\mu} + C_{\kappa}, & \hspace{-8em} \text{if } \vert C\varepsilon^D(\mathbf u^{i-1}) \vert \leq \sigma_0, \\ @@ -219,7 +219,7 @@ where the rank-4 tensor $I_\Pi=I_\Pi(\varepsilon^D(\mathbf u^{i-1}))$ given by This tensor is the (formal) linearization of $P_\Pi(C\cdot)$ around $\varepsilon^D(\mathbf u^{i-1})$. For the linear isotropic material we consider here, the bulk and shear components of the projector are given by -@f{gather*} +@f{gather*}{ C_{\kappa} = \kappa I\otimes I, \qquad\qquad\qquad\qquad C_{\mu} = 2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes @@ -241,7 +241,7 @@ As a final note about the Newton method let us mention that as is common with Newton methods we need to globalize it by controlling the step length. In other words, while the system above solves for $\tilde {\mathbf u}^{i}$, the final iterate will rather be -@f{align*} +@f{align*}{ {\mathbf u}^{i} = {\mathbf u}^{i-1} + \alpha_i (\tilde {\mathbf u}^{i} - {\mathbf u}^{i-1}) @f} where the difference in parentheses on the right takes the role of the @@ -285,7 +285,7 @@ method for the contact. It works as follows: the contact inequality.
  • Find the primal-dual pair $(\tilde U^i,\Lambda^i)$ that satisfies - @f{align*} + @f{align*}{ A\tilde U^i + B\Lambda^i & = F, &\\ \left[B^T\tilde U^i\right]_p & = G_p & \forall p\in\mathcal{A}_i,\\ \Lambda^i_p & = 0 & \forall p\in\mathcal{F}_i. @@ -303,7 +303,7 @@ method for the contact. It works as follows: @f{gather*}U^i \dealcoloneq \alpha^i_l\bar U^i + (1-\alpha^i_l)U^{i-1}@f} satisfies - @f{gather*} + @f{gather*}{ \vert {\hat R}\left({\mathbf u}^{i}\right) \vert < \vert {\hat R}\left({\mathbf u}^{i-1}\right) \vert. \f} with ${\hat R}\left({\mathbf u}\right)=\left(P_{Pi}(C\varepsilon(u)),\varepsilon(\varphi^{i}_p\right)$ with @@ -317,11 +317,11 @@ method for the contact. It works as follows: c\left(\left[B^TU^i\right]_p - G_p\right) \leq 0\rbrace.@f}
  • Project $U^i$ so that it satisfies the contact inequality, - @f{gather*}\hat U^i \dealcoloneq P_{\mathcal{A}_{i+1}}(U^i).@f} + @f{gather*}\hat U^i \dealcoloneq P_{\mathcal{A}_{i+1}}(U^i).@f}{ Here, $P_{\mathcal{A}}(U)$ is the projection of the active components in $\mathcal{A}$ to the gap - @f{gather*}P_{\mathcal{A}}(U)_p \dealcoloneq \begin{cases} + @f{gather*}P_{\mathcal{A}}(U)_p \dealcoloneq \begin{cases}{ U_p, & \textrm{if}\quad p\notin\mathcal{A}\\ g_{h,p}, & \textrm{if}\quad p\in\mathcal{A}, @@ -346,7 +346,7 @@ we can choose $B$ to be a matrix that has only one entry per row, set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg. 194, 2005, pp. 3147-3166). The vector $G$ is defined by a suitable approximation $g_h$ of the gap $g$ -@f{gather*}G_p = \begin{cases} +@f{gather*}G_p = \begin{cases}{ g_{h,p}, & \text{if}\quad p\in\mathcal{S}\\ 0, & \text{if}\quad p\notin\mathcal{S}. \end{cases}@f} diff --git a/examples/step-43/doc/intro.dox b/examples/step-43/doc/intro.dox index 2ec5bde432..8d90311422 100644 --- a/examples/step-43/doc/intro.dox +++ b/examples/step-43/doc/intro.dox @@ -66,13 +66,13 @@ We consider the flow of a two-phase immiscible, incompressible fluid. Capillary and gravity effects are neglected, and viscous effects are assumed dominant. The governing equations for such a flow that are identical to those used in step-21 and are -@f{align*} +@f{align*}{ \mathbf{u}_t &= - \mathbf{K} \lambda_t \left(S\right) \nabla p, \\ \nabla \cdot \mathbf{u}_t &= q, \\ \epsilon \frac{\partial S}{\partial t} + \nabla \cdot \left( \mathbf{u}_t F\left( S \right) \right)&=0, @f} where $S$ is the saturation (volume fraction between zero and one) of the second (wetting) phase, $p$ is the pressure, $\mathbf{K}$ is the permeability tensor, $\lambda_t$ is the total mobility, $\epsilon$ is the porosity, $F$ is the fractional flow of the wetting phase, $q$ is the source term and $\mathbf{u}_t$ is the total velocity. The total mobility, fractional flow of the wetting phase and total velocity are respectively given by -@f{align*} +@f{align*}{ \lambda_t(S)&= \lambda_w + \lambda_{nw} = \frac{k_{rw}(S)}{\mu_w} + \frac{k_{rnw}(S)}{\mu_{nw}}, \\ F(S) &= \frac{\lambda_w}{\lambda_t} = \frac{\lambda_w}{\lambda_w + \lambda_{nw}} = \frac{k_{rw}(S)/\mu_w}{k_{rw}(S)/\mu_w + k_{rnw}(S)/\mu_{nw}}, \\ \mathbf{u}_t &= \mathbf{u}_w + \mathbf{u}_{nw} = -\lambda_t(S)\mathbf{K} \cdot \nabla p, @@ -85,7 +85,7 @@ porosity $\epsilon$ in the saturation equation, which can be considered a scaling factor for the time variable, is set to one. Following a commonly used prescription for the dependence of the relative permeabilities $k_{rw}$ and $k_{rnw}$ on saturation, we use -@f{align*} +@f{align*}{ k_{rw} &= S^2, \qquad&\qquad k_{rnw} &= \left( 1-S \right)^2. @f} @@ -118,7 +118,7 @@ operator splitting" scheme. Here, we use the following a posteriori criterion to decide when to re-compute pressure and velocity variables (detailed derivations and descriptions can be found in @cite Chueh2013): -@f{align*} +@f{align*}{ \theta(n,n_p) = \max_{\kappa\in{\mathbb T}} @@ -159,7 +159,7 @@ steps are of uniform length, and both are chosen adaptively. Using this time discretization, we obtain the following set of equations for each time step from the IMPES approach (see step-21): -@f{align*} +@f{align*}{ \mathbf{u}^{(n)}_t + \lambda_t\left(S^{(n-1)}\right) \mathbf{K} \nabla p^{(n)} =0, \\ \nabla \cdot \mathbf{u}^{(n)}_t = q, \\ \epsilon \left( \frac{S^{(n-1)}-S^{(n)}}{\Delta t^{(n)}_c} \right) + \mathbf{u}^{(n)}_t \cdot \nabla F\left(S^{(n-1)}\right) + F\left(S^{(n-1)}\right) \nabla \cdot \mathbf{u}^{(n)}_t =0. @@ -168,7 +168,7 @@ each time step from the IMPES approach (see step-21): Using the fact that $\nabla \cdot \mathbf{u}_t = q$, the time discrete saturation equation becomes -@f{align*} +@f{align*}{ &\epsilon \left( \frac{S^{(n)}-S^{(n-1)}}{\Delta t^{(n)}_c} \right) + \mathbf{u}^{(n)}_t \cdot \nabla F\left(S^{(n-1)}\right) + F\left(S^{(n-1)}\right)q=0. @f} @@ -181,7 +181,7 @@ functions $\mathbf{v}$ and $w$ respectively and then integrating terms by parts as necessary, the weak form of the problem reads: Find $\mathbf u, p$ so that for all test functions $\mathbf{v}, w$ there holds -@f{gather*} +@f{gather*}{ \left( \left( \mathbf{K} \lambda_t\left(S^{(n-1)}\right) \right)^{-1} \mathbf{u}^{(n)}_t, \mathbf{v}\right)_{\Omega} - \left(p^{(n)}, \nabla \cdot \mathbf{v}\right)_{\Omega} = -\left(p^{(n)}, \mathbf{n} \cdot \mathbf{v} \right)_{\partial \Omega}, \\ - \left( \nabla \cdot \mathbf{u}^{(n)}_t,w\right)_{\Omega} = - \big(q,w\big)_{\Omega}. @f} @@ -220,7 +220,7 @@ validated in @cite Chueh2013 and This method modifies the (discrete) weak form of the saturation equation to read -@f{align*} +@f{align*}{ \left(\epsilon \frac{\partial S_h}{\partial t},\sigma_h\right) - \left(\mathbf{u}_t F\left( S_h \right), @@ -271,7 +271,7 @@ therefore provides for a higher degree of accuracy. On the other hand, it is nonlinear since $\nu$ depends on the saturation $S$. We avoid this difficulty by treating all nonlinear terms explicitly, which leads to the following fully discrete problem at time step $n$: -@f{align*} +@f{align*}{ &\left( \epsilon S_h^{(n)},\sigma_h\right)_{\Omega} - \Delta t^{(n)}_c \Big(F\left(S_h^{(n-1)}\right)\mathbf{u}^{*}_t,\nabla\sigma_h\Big)_{\Omega} + \Delta t^{(n)}_c \Big(F\left(S_h^{(n-1)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{*}_t\right),\sigma_h\Big)_{\partial\Omega} \nonumber \\ & \quad = \left( \epsilon S_h^{(n-1)},\sigma_h\right)_{\Omega} - \Delta t^{(n)}_c \bigg(\nu\left(S_h^{(n-1)}\right)\nabla S_h^{(n-1)},\nabla\sigma_h\bigg)_{\Omega} \nonumber \\ & \qquad + \Delta t^{(n)}_c \bigg(\mathbf{n}\cdot\nu\left(S_h^{(n-1)}\right)\nabla S^{(n-1)},\sigma_h\bigg)_{\partial\Omega} @@ -284,7 +284,7 @@ is to solve with a @ref GlossMassMatrix "mass matrix" on the saturation space. Since the Dirichlet boundary conditions for saturation are only imposed on the inflow boundaries, the third term on the left hand side of the equation above needs to be split further into two parts: -@f{align*} +@f{align*}{ &\Delta t^{(n)}_c \Big(F\left(S_h^{(n-1)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{(n)}_t\right),\sigma_h\Big)_{\partial\Omega} \nonumber \\ &\qquad= \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}_{(+)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{(n)}_{t(+)}\right),\sigma_h\Big)_{\partial\Omega_{(+)}} + \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}_{(-)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{(n)}_{t(-)}\right),\sigma_h\Big)_{\partial\Omega_{(-)}} @f} @@ -349,7 +349,7 @@ obtain a linear system of equations in time step $(n)$ of the following form: \right) @f] where the individual matrices and vectors are defined as follows using shape functions $\mathbf{v}_i$ for velocity, and $\phi_i$ for both pressure and saturation: -@f{align*} +@f{align*}{ \mathbf{M}^{\mathbf{u}}_{ij} &= \left( \left( \mathbf{K} \lambda_t\left(S^{(n-1)}\right) \right)^{-1} \mathbf{v}_{i},\mathbf{v}_{j}\right)_{\Omega}, @@ -399,7 +399,7 @@ We apply the Generalized Minimal Residual (GMRES) method @cite Saad1986 to this linear system. The ideal preconditioner for the velocity-pressure system is -@f{align*} +@f{align*}{ \mathbf{P} = \left( \begin{array}{cc} @@ -420,7 +420,7 @@ where $\mathbf{S}=\mathbf{B}\left(\mathbf{M}^{\mathbf{u}}\right)^{-1}\mathbf{B}^T$ is the Schur complement @cite Zhang2005 of the system. This preconditioner is optimal since -@f{align*} +@f{align*}{ \mathbf{P}^{-1} \left( \begin{array}{cc} @@ -444,7 +444,7 @@ approach by @cite SW94 originally proposed for the Stokes system. (See also the note in the "Possibilities for extensions" section of step-22.) Adapting it to the current set of equations yield the preconditioner -@f{align*} +@f{align*}{ \mathbf{\tilde{P}}^{-1} = \left( \begin{array}{cc} @@ -475,7 +475,7 @@ to be built with Dirichlet boundary conditions to ensure its invertibility. Once the velocity $\mathbf{U}^{(n)} \equiv \mathbf{u}^*_t$ is available, we can assemble $\mathbf{H}$ and $\mathbf{F}_{3}$ and solve for the saturations using -@f{align*} +@f{align*}{ \mathbf{M}^{S} \mathbf{S}^{(n)} = \mathbf{F}_{3} - \mathbf{H} \mathbf{U}^{(n)}. @f} where the mass matrix $\mathbf{M}^{S}$ is solved by the conjugate gradient @@ -523,7 +523,7 @@ Pressure and saturation uniquely determine a velocity, and the velocity determines whether a boundary segment is an inflow or outflow boundary. On the inflow part of the boundary, $\mathbf{\Gamma}_{in}(t)$, we impose -@f{align*} +@f{align*}{ S(\mathbf{x},t) = 1 \qquad & \textrm{on} \quad \mathbf{\Gamma}_{in}(t) \cap \left\{x = 0\right\}, \\ S(\mathbf{x},t) = 0 \qquad & \textrm{on} \quad \mathbf{\Gamma}_{in}(t) \backslash \left\{x = 0\right\}. @f} diff --git a/examples/step-44/doc/intro.dox b/examples/step-44/doc/intro.dox index edd216beb7..46e1fe52a7 100644 --- a/examples/step-44/doc/intro.dox +++ b/examples/step-44/doc/intro.dox @@ -251,7 +251,7 @@ Following the multiplicative decomposition of the deformation gradient, the Helm \Psi(\mathbf{b}) = \Psi_{\text{vol}}(J) + \Psi_{\text{iso}}(\overline{\mathbf{b}}) \, . @f] Similarly, the Kirchhoff stress can be decomposed into volumetric and isochoric parts as $\boldsymbol{\tau} = \boldsymbol{\tau}_{\text{vol}} + \boldsymbol{\tau}_{\text{iso}}$ where: -@f{align*} +@f{align*}{ \boldsymbol{\tau}_{\text{vol}} &= 2 \mathbf{b} \dfrac{\partial \Psi_{\textrm{vol}}(J)}{\partial \mathbf{b}} \\ @@ -325,7 +325,7 @@ can be written in the following decoupled form: \mathfrak{c} = \mathfrak{c}_{\text{vol}} + \mathfrak{c}_{\text{iso}} \, , @f] where -@f{align*} +@f{align*}{ J \mathfrak{c}_{\text{vol}} &= 4 \mathbf{b} \dfrac{\partial^2 \Psi_{\text{vol}}(J)} {\partial \mathbf{b} \partial \mathbf{b}} \mathbf{b} \\ @@ -382,7 +382,7 @@ The body force per unit current volume is denoted $\mathbf{b}^\text{p}$. The stationarity of the potential follows as -@f{align*} +@f{align*}{ R(\mathbf\Xi;\delta \mathbf{\Xi}) &= D_{\delta \mathbf{\Xi}}\Pi(\mathbf{\Xi}) \\ @@ -420,7 +420,7 @@ It should be noted however that they are equivalent. The Euler-Lagrange equations corresponding to the residual are: -@f{align*} +@f{align*}{ &\textrm{div}\ \boldsymbol{\sigma} + \mathbf{b}^\text{p} = \mathbf{0} && \textrm{[equilibrium]} \\ &J(\mathbf{u}) = \widetilde{J} && \textrm{[dilatation]} @@ -432,7 +432,7 @@ The second is the constraint that $J(\mathbf{u}) = \widetilde{J}$. The third is the definition of the pressure $\widetilde{p}$. @note The simplified single-field derivation ($\mathbf{u}$ is the only primary variable) below makes it clear how we transform the limits of integration to the reference domain: -@f{align*} +@f{align*}{ \int_{\Omega}\delta \mathbf{u} \cdot [\textrm{div}\ \boldsymbol{\sigma} + \mathbf{b}^\text{p}]~\mathrm{d}v &= \int_{\Omega} [-\mathrm{grad}\delta \mathbf{u}:\boldsymbol{\sigma} + \delta \mathbf{u} \cdot\mathbf{b}^\text{p}]~\mathrm{d}v @@ -483,7 +483,7 @@ The tangent is given by =: K(\mathbf{\Xi}_{\mathsf{i}}; d \mathbf{\Xi}, \delta \mathbf{\Xi}) \, . @f] Thus, -@f{align*} +@f{align*}{ K(\mathbf{\Xi}_{\mathsf{i}}; d \mathbf{\Xi}, \delta \mathbf{\Xi}) &= D_{d \mathbf{u}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) \cdot d \mathbf{u} @@ -495,7 +495,7 @@ Thus, D_{d \widetilde{J}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) d \widetilde{J} \, , @f} where -@f{align*} +@f{align*}{ D_{d \mathbf{u}} R( \mathbf{\Xi}; \delta \mathbf{\Xi}) &= \int_{\Omega_0} \bigl[ \textrm{grad}\ \delta \mathbf{u} : @@ -518,7 +518,7 @@ where @f} Note that the following terms are termed the geometrical stress and the material contributions to the tangent matrix: -@f{align*} +@f{align*}{ & \int_{\Omega_0} \textrm{grad}\ \delta \mathbf{u} : \textrm{grad}\ d \mathbf{u} [\boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}]~\textrm{d}V && \quad {[\textrm{Geometrical stress}]} \, , @@ -558,7 +558,7 @@ The linearised problem can be written as \mathbf{ \mathsf{F}}(\mathbf{\Xi}_{\textrm{i}}) @f] where -@f{align*} +@f{align*}{ \underbrace{\begin{bmatrix} \mathbf{\mathsf{K}}_{uu} & \mathbf{\mathsf{K}}_{u\widetilde{p}} & \mathbf{0} \\ @@ -596,7 +596,7 @@ by inverting a local matrix and multiplying it by the local right hand side. We can then insert the result into the remaining equations and recover a classical displacement-based method. In order to condense out the pressure and dilatation contributions at the element level we need the following results: -@f{align*} +@f{align*}{ d \widetilde{\mathbf{\mathsf{p}}} & = \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}^{-1} \bigl[ \mathbf{\mathsf{F}}_{\widetilde{J}} diff --git a/examples/step-44/step-44.cc b/examples/step-44/step-44.cc index 0e9fcba5c2..842ce33891 100644 --- a/examples/step-44/step-44.cc +++ b/examples/step-44/step-44.cc @@ -1452,7 +1452,7 @@ namespace Step44 BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); // The global system matrix initially has the following structure - // @f{align*} + // @f{align*}{ // \underbrace{\begin{bmatrix} // \mathsf{\mathbf{K}}_{uu} & \mathsf{\mathbf{K}}_{u\widetilde{p}} & // \mathbf{0} @@ -2503,7 +2503,7 @@ namespace Step44 // the dof associated with the current element // (denoted somewhat loosely as $\mathsf{\mathbf{k}}$) // is of the form: - // @f{align*} + // @f{align*}{ // \begin{bmatrix} // \mathsf{\mathbf{k}}_{uu} & \mathsf{\mathbf{k}}_{u\widetilde{p}} // & \mathbf{0} @@ -2514,7 +2514,7 @@ namespace Step44 // @f} // // We now need to modify it such that it appear as - // @f{align*} + // @f{align*}{ // \begin{bmatrix} // \mathsf{\mathbf{k}}_{\textrm{con}} & // \mathsf{\mathbf{k}}_{u\widetilde{p}} & \mathbf{0} @@ -2684,7 +2684,7 @@ namespace Step44 { // Firstly, here is the approach using the (permanent) augmentation of // the tangent matrix. For the following, recall that - // @f{align*} + // @f{align*}{ // \mathsf{\mathbf{K}}_{\textrm{store}} //\dealcoloneq // \begin{bmatrix} @@ -2697,7 +2697,7 @@ namespace Step44 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}} \end{bmatrix} \, . // @f} // and - // @f{align*} + // @f{align*}{ // d \widetilde{\mathsf{\mathbf{p}}} // & = // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1} diff --git a/examples/step-45/doc/intro.dox b/examples/step-45/doc/intro.dox index 11c6572797..2f7e3c293f 100644 --- a/examples/step-45/doc/intro.dox +++ b/examples/step-45/doc/intro.dox @@ -152,14 +152,14 @@ velocity component of a Stokes flow. On a quarter-circle defined by $\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}$ we are going to solve the Stokes problem -@f{eqnarray*} +@f{eqnarray*}{ -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\ -\textrm{div}\; \textbf{u}&=&0,\\ \textbf{u}|_{\Gamma_1}&=&{\bf 0}, @f} where the boundary $\Gamma_1$ is defined as $\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}$. For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e. -@f{align*} +@f{align*}{ u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\ u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1]. @f} diff --git a/examples/step-45/step-45.cc b/examples/step-45/step-45.cc index f7a44ecd1f..2aabdaaf8d 100644 --- a/examples/step-45/step-45.cc +++ b/examples/step-45/step-45.cc @@ -318,7 +318,7 @@ namespace Step45 // on the lower boundary given by $\text{vertices}_2=R\cdot // \text{vertices}_1+b$ where the rotation matrix $R$ and the offset $b$ are // given by - // @f{align*} + // @f{align*}{ // R=\begin{pmatrix} // 0&1\\-1&0 // \end{pmatrix}, @@ -405,7 +405,7 @@ namespace Step45 // $\text{vertices}_1$ of a face on the lower boundary given by // $\text{vertices}_2=R\cdot \text{vertices}_1+b$ where the rotation // matrix $R$ and the offset $b$ are given by - // @f{align*} + // @f{align*}{ // R=\begin{pmatrix} // 0&1\\-1&0 // \end{pmatrix}, diff --git a/examples/step-46/doc/intro.dox b/examples/step-46/doc/intro.dox index b198c68b93..b99feef93b 100644 --- a/examples/step-46/doc/intro.dox +++ b/examples/step-46/doc/intro.dox @@ -22,20 +22,20 @@ where you may want to read up on the individual equations): - In a part $\Omega_f$ of $\Omega$, we have a fluid flowing that satisfies the time independent Stokes equations (in the form that involves the strain tensor): - @f{align*} + @f{align*}{ -2\eta\nabla \cdot \varepsilon(\mathbf v) + \nabla p &= 0, \qquad \qquad && \text{in}\ \Omega_f\\ -\nabla \cdot \mathbf v &= 0 && \text{in}\ \Omega_f. @f} Here, $\mathbf v, p$ are the fluid velocity and pressure, respectively. We prescribe the velocity on part of the external boundary, - @f{align*} + @f{align*}{ \mathbf v = \mathbf v_0 \qquad\qquad \text{on}\ \Gamma_{f,1} \subset \partial\Omega \cap \partial\Omega_f @f} while we assume free-flow conditions on the remainder of the external boundary, - @f{align*} + @f{align*}{ (2\eta \varepsilon(\mathbf v) - p \mathbf 1) \cdot \mathbf n = 0 \qquad\qquad \text{on}\ \Gamma_{f,2} = \partial\Omega \cap \partial\Omega_f \backslash @@ -44,7 +44,7 @@ where you may want to read up on the individual equations): - The remainder of the domain, $\Omega_s = \Omega \backslash \Omega_f$ is occupied by a solid whose deformation field $\mathbf u$ satisfies the elasticity equation, - @f{align*} + @f{align*}{ -\nabla \cdot C \varepsilon(\mathbf u) = 0 \qquad\qquad & \text{in}\ \Omega_s, @f} @@ -55,19 +55,19 @@ where you may want to read up on the individual equations): be so small that it has no feedback effect on the fluid, i.e. the coupling is only in one direction. For simplicity, we will assume that the solid's external boundary is clamped, i.e. - @f{align*} + @f{align*}{ \mathbf u = \mathbf 0 \qquad\qquad \text{on}\ \Gamma_{s,1} = \partial\Omega \cap \partial\Omega_s @f} - As a consequence of the small displacement assumption, we will pose the following boundary conditions on the interface between the fluid and solid: first, we have no slip boundary conditions for the fluid, - @f{align*} + @f{align*}{ \mathbf v = \mathbf 0 \qquad\qquad \text{on}\ \Gamma_{i} = \partial\Omega_s \cap \partial\Omega_f. @f} Secondly, the forces (traction) on the solid equal the normal stress from the fluid, - @f{align*} + @f{align*}{ (C \varepsilon(\mathbf u)) \mathbf n = (2 \eta \varepsilon(\mathbf v) - p \mathbf 1) \mathbf n \qquad\qquad \text{on}\ \Gamma_{i} = \partial\Omega_s \cap \partial\Omega_f, @@ -80,7 +80,7 @@ multiplying from the left by a test function and integrating over the domain. It then looks like this: Find $y = \{\mathbf v, p, \mathbf u\} \in Y \subset H^1(\Omega_f)^d \times L_2(\Omega_f) \times H^1(\Omega_s)^d$ such that -@f{align*} +@f{align*}{ 2 \eta (\varepsilon(\mathbf a), \varepsilon(\mathbf v))_{\Omega_f} - (\nabla \cdot \mathbf a, p)_{\Omega_f} - (q, \nabla \cdot \mathbf v)_{\Omega_f} & @@ -130,7 +130,7 @@ appropriate function space for these variables? We know that on $\Omega_f$ we should require $\mathbf v \in H^1(\Omega_f)^d, p \in L_2(\Omega_f)$, so for the extensions $\tilde{\mathbf v}, \tilde p$ to the whole domain the following appears a useful set of function spaces: -@f{align*} +@f{align*}{ \tilde {\mathbf v} &\in V = \{\tilde {\mathbf v}|_{\Omega_f} \in H^1(\Omega_f)^d, \quad \tilde {\mathbf v}|_{\Omega_s} = 0 \} @@ -156,7 +156,7 @@ $V, P$. For Stokes, we know from step-22 that an appropriate choice is $Q_{p+1}^d\times Q_P$ but this only holds for that part of the domain occupied by the fluid. For the extended field, let's use the following subspaces defined on the triangulation $\mathbb T$: -@f{align*} +@f{align*}{ V_h &= \{{\mathbf v}_h \quad | \quad \forall K \in {\mathbb T}: @@ -193,13 +193,13 @@ zero. Obviously, we choose $Z_h=Z$. This entire discussion above can be repeated for the variables we use to describe the elasticity equation. Here, for the extended variables, we have -@f{align*} +@f{align*}{ \tilde {\mathbf u} &\in U = \{\tilde {\mathbf u}|_{\Omega_s} \in H^1(\Omega_f)^d, \quad \tilde {\mathbf u}|_{\Omega_f} \in Z(\Omega_s)^d \}, @f} and we will typically use a finite element space of the kind -@f{align*} +@f{align*}{ U_h &= \{{\mathbf u}_h \quad | \quad \forall K \in {\mathbb T}: @@ -213,7 +213,7 @@ of polynomial degree $r$. So to sum up, we are going to look for a discrete vector-valued solution $y_h = \{\mathbf v_h, p_h, \mathbf u_h\}$ in the following space: -@f{align*} +@f{align*}{ Y_h = \{ & y_h = \{\mathbf v_h, p_h, \mathbf u_h\} : \\ & y_h|_{\Omega_f} \in Q_{p+1}^d \times Q_p \times Z^d, \\ @@ -302,7 +302,7 @@ points: Let us first discuss implementing the bilinear form, which at the discrete level we recall to be -@f{align*} +@f{align*}{ 2 \eta (\varepsilon(\mathbf a_h), \varepsilon(\mathbf v_h))_{\Omega_f} - (\nabla \cdot \mathbf a_h, p_h)_{\Omega_f} - (q_h, \nabla \cdot \mathbf v_h)_{\Omega_f} & diff --git a/examples/step-5/doc/intro.dox b/examples/step-5/doc/intro.dox index 9e02ee8249..c8982a34f3 100644 --- a/examples/step-5/doc/intro.dox +++ b/examples/step-5/doc/intro.dox @@ -42,7 +42,7 @@ them are: The equation to solve here is as follows: -@f{align*} +@f{align*}{ -\nabla \cdot a(\mathbf x) \nabla u(\mathbf x) &= 1 \qquad\qquad & \text{in}\ \Omega, \\ u &= 0 \qquad\qquad & \text{on}\ \partial\Omega. @@ -67,7 +67,7 @@ more interpretations than just the two listed above. When assembling the linear system for this equation, we need the weak form which here reads as follows: -@f{align*} +@f{align*}{ (a \nabla \varphi, \nabla u) &= (\varphi, 1) \qquad \qquad \forall \varphi. @f} The implementation in the assemble_system function follows diff --git a/examples/step-50/doc/intro.dox b/examples/step-50/doc/intro.dox index 10c6c0dad8..e95f5d2a85 100644 --- a/examples/step-50/doc/intro.dox +++ b/examples/step-50/doc/intro.dox @@ -55,7 +55,7 @@ strongly when the number of unknowns per processor becomes smaller than

    The testcase

    We consider the variable-coefficient Laplacian weak formulation -@f{align*} +@f{align*}{ (\epsilon \nabla u, \nabla v) = (f,v) \quad \forall v \in V_h @f} on the domain $\Omega = [-1,1]^\text{dim} \setminus [0,1]^\text{dim}$ (an L-shaped domain @@ -69,7 +69,7 @@ the right-hand side is $f=1$. We use continuous $Q_2$ elements for the discrete finite element space $V_h$, and use a residual-based, cell-wise a posteriori error estimator $e(K) = e_{\text{cell}}(K) + e_{\text{face}}(K)$ from @cite karakashian2003posteriori with -@f{align*} +@f{align*}{ e_{\text{cell}}(K) &= h^2 \| f + \epsilon \triangle u \|_K^2, \\ e_{\text{face}}(K) &= \sum_F h_F \| \jump{ \epsilon \nabla u \cdot n } \|_F^2, @f} @@ -144,26 +144,26 @@ level. Now, let $N_{\ell}$ be the number of cells on level $\ell$ $p$. We will also denote by $P$ the total number of processors. Assuming that the workload for any one processor is proportional to the number of cells owned by that processor, the optimal workload per processor is given by -@f{align*} +@f{align*}{ W_{\text{opt}} = \frac1{P}\sum_{\ell} N_{\ell} = \sum_{\ell}\left(\frac1{P}\sum_{p}N_{\ell,p}\right). @f} Next, assuming a synchronization of work on each level (i.e., on each level of a V-cycle, work must be completed by all processors before moving on to the next level), the limiting effort on each level is given by -@f{align*} +@f{align*}{ W_\ell = \max_{p} N_{\ell,p}, @f} and the total parallel complexity -@f{align*} +@f{align*}{ W = \sum_{\ell} W_\ell. @f} Then we define $\mathbb{E}$ as a ratio of the optimal partition to the parallel complexity of the current partition -@f{align*} +@f{align*}{ \mathbb{E} = \frac{W_{\text{opt}}}{W}. @f} For the example distribution above, we have -@f{align*} +@f{align*}{ W_{\text{opt}}&=\frac{1}{P}\sum_{\ell} N_{\ell} = \frac{1}{3} \left(1+4+4\right)= 3 \qquad \\ W &= \sum_\ell W_\ell = 1 + 2 + 3 = 6 diff --git a/examples/step-51/doc/intro.dox b/examples/step-51/doc/intro.dox index 85b4bbce94..9b9466b699 100644 --- a/examples/step-51/doc/intro.dox +++ b/examples/step-51/doc/intro.dox @@ -71,7 +71,7 @@ the hybridizable discontinuous Galerkin (HDG) method. Let us write the complete linear system associated to the HDG problem as a block system with the discrete DG (cell interior) variables $U$ as first block and the skeleton (face) variables $\Lambda$ as the second block: -@f{eqnarray*} +@f{eqnarray*}{ \begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} U \\ \Lambda \end{pmatrix} = @@ -79,7 +79,7 @@ and the skeleton (face) variables $\Lambda$ as the second block: @f} Our aim is now to eliminate the $U$ block with a Schur complement approach similar to step-20, which results in the following two steps: -@f{eqnarray*} +@f{eqnarray*}{ (D - C A^{-1} B) \Lambda &=& G - C A^{-1} F, \\ A U &=& F - B \Lambda. @f} @@ -132,7 +132,7 @@ The HDG formulation used for this example is taken from We consider the convection-diffusion equation over the domain $\Omega$ with Dirichlet boundary $\partial \Omega_D$ and Neumann boundary $\partial \Omega_N$: -@f{eqnarray*} +@f{eqnarray*}{ \nabla \cdot (\mathbf{c} u) - \nabla \cdot (\kappa \nabla u) &=& f, \quad \text{ in } \Omega, \\ u &=& g_D, \quad \text{ on } \partial \Omega_D, \\ @@ -142,7 +142,7 @@ $\partial \Omega_N$: Introduce the auxiliary variable $\mathbf{q}=-\kappa \nabla u$ and rewrite the above equation as the first order system: -@f{eqnarray*} +@f{eqnarray*}{ \mathbf{q} + \kappa \nabla u &=& 0, \quad \text{ in } \Omega, \\ \nabla \cdot (\mathbf{c} u + \mathbf{q}) &=& f, \quad \text{ in } \Omega, \\ u &=& g_D, \quad \text{ on } \partial \Omega_D, \\ @@ -152,7 +152,7 @@ the above equation as the first order system: We multiply these equations by the weight functions $\mathbf{v}, w$ and integrate by parts over every element $K$ to obtain: -@f{eqnarray*} +@f{eqnarray*}{ (\mathbf{v}, \kappa^{-1} \mathbf{q})_K - (\nabla\cdot\mathbf{v}, u)_K + \left<\mathbf{v}\cdot\mathbf{n}, {\hat{u}}\right>_{\partial K} &=& 0, \\ - (\nabla w, \mathbf{c} u + \mathbf{q})_K @@ -167,7 +167,7 @@ these terms must be single-valued on any given element edge $\partial K$ even though, with discontinuous shape functions, there may of course be multiple values coming from the cells adjacent to an interface. We eliminate the numerical trace $\hat{\mathbf{q}}$ by using traces of the form: -@f{eqnarray*} +@f{eqnarray*}{ \widehat{\mathbf{c} u}+\hat{\mathbf{q}} = \mathbf{c}\hat{u} + \mathbf{q} + \tau(u - \hat{u})\mathbf{n} \quad \text{ on } \partial K. @f} @@ -188,7 +188,7 @@ stabilization parameter $\tau$ that tends to infinity prohibits jumps in the solution over the element boundaries, making the HDG solution approach the approximation with continuous finite elements. In the program below, we choose the stabilization parameter as -@f{eqnarray*} +@f{eqnarray*}{ \tau = \frac{\kappa}{\ell} + |\mathbf{c} \cdot \mathbf{n}| @f} where we set the diffusion $\kappa=1$ and the diffusion length scale to @@ -197,7 +197,7 @@ $\ell = \frac{1}{5}$. The trace/skeleton variables in HDG methods are single-valued on element faces. As such, they must strongly represent the Dirichlet data on $\partial\Omega_D$. This means that -@f{equation*} +@f{equation*}{ \hat{u}|_{\partial \Omega_D} = g_D, @f} where the equal sign actually means an $L_2$ projection of the boundary @@ -212,7 +212,7 @@ component of the numerical flux, and integrating by parts on the equation weighted by $w$, we arrive at the final form of the problem: Find $(\mathbf{q}_h, u_h, \hat{u}_h) \in \mathcal{V}_h^p \times \mathcal{W}_h^p \times \mathcal{M}_h^p$ such that -@f{align*} +@f{align*}{ (\mathbf{v}, \kappa^{-1} \mathbf{q}_h)_{\mathcal{T}} - ( \nabla\cdot\mathbf{v}, u_h)_{\mathcal{T}} + \left<\mathbf{v}\cdot\mathbf{n}, \hat{u}_h\right>_{\partial\mathcal{T}} @@ -301,7 +301,7 @@ $K$ under the constraint $\left(1, u_h^*\right)_K = \left(1, u_h\right)_K$. The constraint is necessary because the minimization functional does not determine the constant part of $u_h^*$. This translates to the following system of equations: -@f{eqnarray*} +@f{eqnarray*}{ \left(1, u_h^*\right)_K &=& \left(1, u_h\right)_K\\ \left(\nabla w_h^*, \kappa \nabla u_h^*\right)_K &=& -\left(\nabla w_h^*, \mathbf{q}_h\right)_K diff --git a/examples/step-52/doc/intro.dox b/examples/step-52/doc/intro.dox index 13c47faeb3..6097ff2a4c 100644 --- a/examples/step-52/doc/intro.dox +++ b/examples/step-52/doc/intro.dox @@ -26,7 +26,7 @@ time-dependent diffusion equation -- which is just a different name for the heat equation discussed in step-26, plus some extra terms. We assume that the medium is not fissible and therefore, the neutron flux satisfies the following equation: -@f{eqnarray*} +@f{eqnarray*}{ \frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla \cdot D(x) \nabla \phi(x,t) - \Sigma_a(x) \phi(x,t) + S(x,t) @f} @@ -41,7 +41,7 @@ Since this program only intends to demonstrate how to use advanced time stepping algorithms, we will only look for the solutions of relatively simple problems. Specifically, we are looking for a solution on a square domain $[0,b]\times[0,b]$ of the form -@f{eqnarray*} +@f{eqnarray*}{ \phi(x,t) = A\sin(\omega t)(bx-x^2). @f} By using quadratic finite elements, we can represent this function exactly at @@ -54,7 +54,7 @@ We impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and $x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. We choose the source term so that the corresponding solution is in fact of the form stated above: -@f{eqnarray*} +@f{eqnarray*}{ S=A\left(\frac{1}{v}\omega \cos(\omega t)(bx -x^2) + \sin(\omega t) \left(\Sigma_a (bx-x^2)+2D\right) \right). @f} @@ -73,25 +73,25 @@ consequently observe the temporal convergence order with either. The Runge-Kutta methods implemented in deal.II assume that the equation to be solved can be written as: -@f{eqnarray*} +@f{eqnarray*}{ \frac{dy}{dt} = g(t,y). @f} On the other hand, when using finite elements, discretized time derivatives always result in the presence of a @ref GlossMassMatrix "mass matrix" on the left hand side. This can easily be seen by considering that if the solution vector $y(t)$ in the equation above is in fact the vector of nodal coefficients $U(t)$ for a variable of the form -@f{eqnarray*} +@f{eqnarray*}{ u_h(x,t) = \sum_j U_j(t) \varphi_j(x) @f} with spatial shape functions $\varphi_j(x)$, then multiplying an equation of the form -@f{eqnarray*} +@f{eqnarray*}{ \frac{\partial u(x,t)}{\partial t} = q(t,u(x,t)) @f} by test functions, integrating over $\Omega$, substituting $u\rightarrow u_h$ and restricting the test functions to the $\varphi_i(x)$ from above, then this spatially discretized equation has the form -@f{eqnarray*} +@f{eqnarray*}{ M\frac{dU}{dt} = f(t,U), @f} where $M$ is the mass matrix and $f(t,U)$ is the spatially discretized version @@ -99,17 +99,17 @@ of $q(t,u(x,t))$ (where $q$ is typically the place where spatial derivatives appear, but this is not of much concern for the moment given that we only consider time derivatives). In other words, this form fits the general scheme above if we write -@f{eqnarray*} +@f{eqnarray*}{ \frac{dy}{dt} = g(t,y) = M^{-1}f(t,y). @f} Runke-Kutta methods are time stepping schemes that approximate $y(t_n)\approx y_{n}$ through a particular one-step approach. They are typically written in the form -@f{eqnarray*} +@f{eqnarray*}{ y_{n+1} = y_n + \sum_{i=1}^s b_i k_i @f} where for the form of the right hand side above -@f{eqnarray*} +@f{eqnarray*}{ k_i = \Delta t \, M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right). @f} Here $a_{ij}$, $b_i$, and $c_i$ are known coefficients that identify which @@ -180,7 +180,7 @@ $\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M$. The particular form of this operator results from the fact that each Newton step requires the solution of an equation of the form -@f{align*} +@f{align*}{ \left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y = -M h(t,y) @f} @@ -201,16 +201,16 @@ be re-used between stages because $\tau$ is the same every time). By expanding the solution of our model problem as always using shape functions $\psi_j$ and writing -@f{eqnarray*} +@f{eqnarray*}{ \phi_h(x,t) = \sum_j U_j(t) \psi_j(x), @f} we immediately get the spatially discretized version of the diffusion equation as -@f{eqnarray*} +@f{eqnarray*}{ M \frac{dU(t)}{dt} = -{\cal D} U(t) - {\cal A} U(t) + {\cal S}(t) @f} where -@f{eqnarray*} +@f{eqnarray*}{ M_{ij} &=& (\psi_i,\psi_j), \\ {\cal D}_{ij} &=& (D\nabla\psi_i,\nabla\psi_j)_\Omega, \\ {\cal A}_{ij} &=& (\Sigma_a\psi_i,\psi_j)_\Omega, \\ @@ -220,12 +220,12 @@ See also step-24 and step-26 to understand how we arrive here. Boundary terms are not necessary due to the chosen boundary conditions for the current problem. To use the Runge-Kutta methods, we recast this as follows: -@f{eqnarray*} +@f{eqnarray*}{ f(y) = -{\cal D}y - {\cal A}y + {\cal S}. @f} In the code, we will need to be able to evaluate this function $f(U)$ along with its derivative, -@f{eqnarray*} +@f{eqnarray*}{ \frac{\partial f}{\partial y} = -{\cal D} - {\cal A}. @f} diff --git a/examples/step-55/doc/intro.dox b/examples/step-55/doc/intro.dox index 4877a22847..32bc359021 100644 --- a/examples/step-55/doc/intro.dox +++ b/examples/step-55/doc/intro.dox @@ -50,7 +50,7 @@ The learning outcomes for this tutorial are: We are solving for a velocity $\textbf{u}$ and pressure $p$ that satisfy the Stokes equation, which reads -@f{eqnarray*} +@f{eqnarray*}{ - \triangle \textbf{u} + \nabla p &=& \textbf{f}, \\ -\textrm{div}\; \textbf{u} &=& 0. @f} @@ -93,7 +93,7 @@ bigger, the former will eventually beat the latter.

    The solver and preconditioner

    We precondition the linear system -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) @@ -107,7 +107,7 @@ We precondition the linear system @f} with the block diagonal preconditioner -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} = \left(\begin{array}{cc} diff --git a/examples/step-56/doc/intro.dox b/examples/step-56/doc/intro.dox index ee002f95cc..d1a6f70f00 100644 --- a/examples/step-56/doc/intro.dox +++ b/examples/step-56/doc/intro.dox @@ -33,7 +33,7 @@ Let $u \in H_0^1 = \{ u \in H^1(\Omega), u|_{\partial \Omega} = 0 \}$ and $p \in L_*^2 = \{ p \in L^2(\Omega), \int_\Omega p = 0 \}$. The Stokes equations read as follows in non-dimensionalized form: -@f{eqnarray*} +@f{eqnarray*}{ - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p & =& f \\ - \nabla \cdot u &=& 0 @@ -49,7 +49,7 @@ well as more expensive). The weak form of the discrete equations naturally leads to the following linear system for the nodal values of the velocity and pressure fields: -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) = \left(\begin{array}{c} F \\ 0 \end{array}\right). @@ -63,7 +63,7 @@ preconditioner, in the spirit of the approach outlined in the "Results" section of step-22. The idea is as follows: if we find a block preconditioner $P$ such that the matrix -@f{eqnarray*} +@f{eqnarray*}{ \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1} @f} @@ -72,14 +72,14 @@ converge in a few iterations. Notice that we are doing right preconditioning here. Using the Schur complement $S=BA^{-1}B^T$, we find that -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} = \left(\begin{array}{cc} A & B^T \\ 0 & S \end{array}\right)^{-1} @f} is a good choice. Let $\widetilde{A^{-1}}$ be an approximation of $A^{-1}$ and $\widetilde{S^{-1}}$ of $S^{-1}$, we see -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} = \left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right) \left(\begin{array}{cc} I & B^T \\ 0 & -I \end{array}\right) @@ -136,13 +136,13 @@ addition, the pressure is chosen to have a mean value of zero. For the "Method of Manufactured Solutions" of step-7, we need to find $\bf f$ such that: -@f{align*} +@f{align*}{ {\bf f} = - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p. @f} Using the reference solution above, we obtain: -@f{eqnarray*} +@f{eqnarray*}{ {\bf f} &=& (2 \pi^2 \sin (\pi x),- \pi^3 y \cos(\pi x),- \pi^3 z \cos(\pi x))\\ & & + (\pi \cos(\pi x) \cos(\pi y) \sin(\pi z) ,- \pi \sin(\pi y) \sin(\pi x) \sin(\pi z), \pi \cos(\pi diff --git a/examples/step-57/doc/intro.dox b/examples/step-57/doc/intro.dox index ea54a8fe17..8d42d53153 100644 --- a/examples/step-57/doc/intro.dox +++ b/examples/step-57/doc/intro.dox @@ -22,7 +22,7 @@ is assumed to be steady. In a domain $\Omega \subset $\partial \Omega$, and a given force field $\textbf{f}$, we seek a velocity field $\textbf{u}$ and a pressure field $\textbf{p}$ satisfying -@f{eqnarray*} +@f{eqnarray*}{ - \nu \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f}\\ - \nabla \cdot \textbf{u} &=& 0. @f} @@ -38,7 +38,7 @@ solve is necessary.

    Linearization of Navier-Stokes Equations

    We define a nonlinear function whose root is a solution to the NSE by -@f{eqnarray*} +@f{eqnarray*}{ F(\mathbf{u}, p) = \begin{pmatrix} - \nu \Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p - \mathbf{f} \\ @@ -50,7 +50,7 @@ Assuming the initial guess is good enough to guarantee the convergence of Newton's iteration and denoting $\textbf{x} = (\textbf{u}, p)$, Newton's iteration on a vector function can be defined as -@f{eqnarray*} +@f{eqnarray*}{ \textbf{x}^{k+1} = \textbf{x}^{k} - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}), @f} @@ -64,21 +64,21 @@ The Newton iteration formula implies the new solution is obtained by adding an update term to the old solution. Instead of evaluating the Jacobian matrix and taking its inverse, we consider the update term as a whole, that is -@f{eqnarray*} +@f{eqnarray*}{ \delta \textbf{x}^{k} = - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}), @f} where $\textbf{x}^{k+1}=\textbf{x}^{k}+\delta \textbf{x}^{k}$. We can find the update term by solving the system -@f{eqnarray*} +@f{eqnarray*}{ \nabla F(\textbf{x}^{k}) \delta \textbf{x}^{k} = -F(\textbf{x}^{k}). @f} Here, the left of the previous equation represents the directional gradient of $F(\textbf{x})$ along $\delta \textbf{x}^{k}$ at $\textbf{x}^{k}$. By definition, the directional gradient is given by -@f{eqnarray*} +@f{eqnarray*}{ & &\nabla F(\mathbf{u}^{k}, p^{k}) (\delta \mathbf{u}^{k}, \delta p^{k}) \\ \\ &=& \lim_{\epsilon \to 0} \frac{1}{\epsilon} @@ -108,7 +108,7 @@ directional gradient of $F(\textbf{x})$ along $\delta @f} Therefore, we arrive at the linearized system: -@f{eqnarray*} +@f{eqnarray*}{ -\nu \Delta \delta \mathbf{u}^{k} + \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k} + \delta \mathbf{u}^{k} \cdot \nabla \mathbf{u}^{k} @@ -166,13 +166,13 @@ $\nu_{i}$ as the initial guess of the NSE with $\nu_{i+1}$. This can be thought as a staircase from the Stokes equations to the NSE we want to solve. That is, we first solve a Stokes problem -@f{eqnarray*} +@f{eqnarray*}{ -\nu_{1} \Delta \textbf{u} + \nabla p &=& \textbf{f}\\ -\nabla \cdot \textbf{u} &=& 0 @f} to get the initial guess for -@f{eqnarray*} +@f{eqnarray*}{ -\nu_{1} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ -\nabla \cdot \textbf{u} &=& 0, @f} @@ -182,13 +182,13 @@ Here $\nu_{1}$ is relatively large so that the solution to the Stokes problem wi can be used as an initial guess for the NSE in Newton's iteration. Then the solution to -@f{eqnarray*} +@f{eqnarray*}{ -\nu_{i} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ -\nabla \cdot \textbf{u} &=& 0. @f} acts as the initial guess for -@f{eqnarray*} +@f{eqnarray*}{ -\nu_{i+1} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ -\nabla \cdot \textbf{u} &=& 0. @f} @@ -201,7 +201,7 @@ guess for the Newton iteration. At each step of Newton's iteration, the problem results in solving a saddle point systems of the form -@f{eqnarray*} +@f{eqnarray*}{ \begin{pmatrix} A & B^{T} \\ B & 0 @@ -220,7 +220,7 @@ saddle point systems of the form This system matrix has the same block structure as the one in step-22. However, the matrix $A$ at the top left corner is not symmetric because of the nonlinear term. Instead of solving the above system, we can solve the equivalent system -@f{eqnarray*} +@f{eqnarray*}{ \begin{pmatrix} A + \gamma B^TW^{-1}B & B^{T} \\ B & 0 @@ -242,7 +242,7 @@ $\gamma B^TW^{-1}B$ is the Augmented Lagrangian term; see @cite Benzi2006 for de Denoting the system matrix of the new system by $G$ and the right-hand side by $b$, we solve it iteratively with right preconditioning $P^{-1}$ as $GP^{-1}y = b$, where -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} = \begin{pmatrix} \tilde{A} & B^T \\ @@ -254,14 +254,14 @@ with $\tilde{A} = A + \gamma B^TW^{-1}B$ and $\tilde{S}$ is the corresponding Schur complement $\tilde{S} = B^T \tilde{A}^{-1} B$. We let $W = M_p$ where $M_p$ is the pressure @ref GlossMassMatrix "mass matrix", then $\tilde{S}^{-1}$ can be approximated by -@f{eqnarray*} +@f{eqnarray*}{ \tilde{S}^{-1} \approx -(\nu+\gamma)M_p^{-1}. @f} See @cite Benzi2006 for details. We decompose $P^{-1}$ as -@f{eqnarray*} +@f{eqnarray*}{ P^{-1} = \begin{pmatrix} \tilde{A}^{-1} & 0 \\ @@ -293,7 +293,7 @@ We use the lid driven cavity flow as our test case; see [this page](http://www.cfd-online.com/Wiki/Lid-driven_cavity_problem) for details. The computational domain is the unit square and the right-hand side is $f=0$. The boundary condition is -@f{eqnarray*} +@f{eqnarray*}{ (u(x, y), v(x,y)) &=& (1,0) \qquad\qquad \textrm{if}\ y = 1 \\ (u(x, y), v(x,y)) &=& (0,0) \qquad\qquad \textrm{otherwise}. @f} diff --git a/examples/step-59/doc/intro.dox b/examples/step-59/doc/intro.dox index f78b7930f5..6bfb94959c 100644 --- a/examples/step-59/doc/intro.dox +++ b/examples/step-59/doc/intro.dox @@ -49,7 +49,7 @@ For this tutorial program, we exemplify the matrix-free DG framework for the interior penalty discretization of the Laplacian, i.e., the same scheme as the one used for the step-39 tutorial program. The discretization of the Laplacian is given by the following weak form -@f{align*} +@f{align*}{ &\sum_{K\in\text{cells}} \left(\nabla v_h, \nabla u_h\right)_{K}+\\ &\sum_{F\in\text{faces}}\Big(-\left<\jump{v_h}, \average{\nabla u_h}\right>_{F} - \left<\average{\nabla v_h}, \jump{u_h}\right>_{F} + \left<\jump{v_h}, \sigma \jump{u_h}\right>_{F}\Big) \\ &= \sum_{K\in\text{cells}}\left(v_h, f\right)_{K}, diff --git a/examples/step-61/doc/intro.dox b/examples/step-61/doc/intro.dox index 04ab0a81f9..dc2f339add 100644 --- a/examples/step-61/doc/intro.dox +++ b/examples/step-61/doc/intro.dox @@ -98,16 +98,16 @@ Then we will calculate $L_2$ errors of pressure, velocity, and flux. The Poisson equation above has a solution $p$ that needs to satisfy the weak formulation of the problem, -@f{equation*} +@f{equation*}{ \mathcal{A}\left(p,q \right) = \mathcal{F} \left(q \right), @f} for all test functions $q$, where -@f{equation*} +@f{equation*}{ \mathcal{A}\left(p,q\right) \dealcoloneq \int_\Omega \left(\mathbf{K} \nabla p\right) \cdot \nabla q \;\mathrm{d}x, @f} and -@f{equation*} +@f{equation*}{ \mathcal{F}\left(q\right) \dealcoloneq \int_\Omega f \, q \;\mathrm{d}x - \int_{\Gamma^N} u_N q \; \mathrm{d}x. @@ -135,17 +135,17 @@ test function $q_h$ on the boundary (where we would simply take its interface part $q_h^\partial$) but we have to be careful with the gradient because that is only defined in cell interiors. Consequently, the weak Galerkin scheme for the Poisson equation is defined by -@f{equation*} +@f{equation*}{ \mathcal{A}_h\left(p_h,q \right) = \mathcal{F} \left(q_h \right), @f} for all discrete test functions $q_h$, where -@f{equation*} +@f{equation*}{ \mathcal{A}_h\left(p_h,q_h\right) \dealcoloneq \sum_{K \in \mathbb{T}} \int_K \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q_h \;\mathrm{d}x, @f} and -@f{equation*} +@f{equation*}{ \mathcal{F}\left(q_h\right) \dealcoloneq \sum_{K \in \mathbb{T}} \int_K f \, q_h^\circ \;\mathrm{d}x - \sum_{\gamma \in \Gamma_h^N} \int_\gamma u_N q_h^\partial \;\mathrm{d}x, @@ -245,7 +245,7 @@ shape functions will be of either one or the other kind. That is not important, here, however. What is important is that we need to wonder how we can represent $\nabla_{w,d} \varphi_j$ because that is clearly what will appear in the problem when we want to implement the bilinear form -@f{equation*} +@f{equation*}{ \mathcal{A}_h\left(p_h,q_h\right) = \sum_{K \in \mathbb{T}} \int_K \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q_h \;\mathrm{d}x, @@ -254,7 +254,7 @@ problem when we want to implement the bilinear form The key point is that $\nabla_{w,d} \varphi_j$ is known to be a member of the "broken" Raviart-Thomas space $DGRT_s$. What this means is that we can represent (on each cell $K$ separately) -@f{equation*} +@f{equation*}{ \nabla_{w,d} \varphi_j|_K = \sum_k C_{jk}^K \mathbf v_k|_K @f} @@ -402,7 +402,7 @@ use the discrete weak gradient of $p_h$ to calculate the velocity on each elemen As discussed above, on each element the gradient of the numerical pressure $\nabla p$ can be approximated by discrete weak gradients $ \nabla_{w,d}\phi_i$: -@f{equation*} +@f{equation*}{ \nabla_{w,d} p_h = \nabla_{w,d} \left(\sum_{i} P_i \phi_i\right) = \sum_{i} P_i \nabla_{w,d}\phi_i. @@ -433,12 +433,12 @@ $\left( \mathbf{Q}_h \left( \mathbf{Kv}_j \right),\mathbf{v}_k \right)_K = \left( \mathbf{Kv}_j,\mathbf{v}_k \right)_K.$ So, rather than the formula shown above, the numerical velocity on cell $K$ instead becomes -@f{equation*} +@f{equation*}{ \mathbf{u}_h = \mathbf{Q}_h \left( -\mathbf{K}\nabla_{w,d}p_h \right) = -\sum_i \sum_j P_i B^K_{ij}\mathbf{Q}_h \left( \mathbf{K}\mathbf{v}_j \right), @f} and we have the following system to solve for the coefficients $d_{jk}$: -@f{equation*} +@f{equation*}{ \sum_j \left(\mathbf{v}_i,\mathbf{v}_j\right) d_{jk} @@ -457,7 +457,7 @@ $ is called cell_matrix_E. Then the elementwise velocity is -@f{equation*} +@f{equation*}{ \mathbf{u}_h = -\sum_{i} \sum_{j}P_ic_{ij}\sum_{k}d_{jk}\mathbf{v}_k = \sum_{k}- \left(\sum_{j} \sum_{i} P_ic_{ij}d_{jk} \right)\mathbf{v}_k, @f} diff --git a/examples/step-62/doc/intro.dox b/examples/step-62/doc/intro.dox index e0164ca8fd..b2fb90ec0b 100644 --- a/examples/step-62/doc/intro.dox +++ b/examples/step-62/doc/intro.dox @@ -106,7 +106,7 @@ Instead of a time domain approach, this tutorial program converts the equations above into the frequency domain by performing a Fourier transform with regard to the time variable. The elastic equations in the frequency domain then read as follows -@f{eqnarray*} +@f{eqnarray*}{ \nabla\cdot(\boldsymbol{\bar\sigma} \xi \boldsymbol{\Lambda})&=&-\omega^2\rho\xi\mathbf{\bar u}\\ \boldsymbol{\bar \sigma} &=&\mathbf{C}\boldsymbol{\bar\varepsilon}\\ \boldsymbol{\bar\varepsilon}&=&\frac{1}{2}[(\nabla\mathbf{\bar{u}}\boldsymbol{\Lambda}+\boldsymbol{\Lambda}^\mathrm{T}(\nabla\mathbf{\bar{u}})^\mathrm{T})]\\ @@ -152,14 +152,14 @@ It is useful to introduce the tensors $\alpha_{mnkl}$ and $\beta_{mnkl}$. @f] We can multiply by $\varphi_m$ and integrate over the domain $\Omega$ and integrate by parts. -@f{eqnarray*} +@f{eqnarray*}{ -\omega^2\int_\Omega\rho\xi\varphi_m u_m + \int_\Omega\partial_n\varphi_m \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k u_l + \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l u_k\right) = \int_\Omega\varphi_m f_m @f} It is this set of equations we want to solve for a set of frequencies $\omega$ in order to compute the transmission coefficient as function of frequency. The linear system becomes -@f{eqnarray*} +@f{eqnarray*}{ AU&=&F\\ A_{ij} &=& -\omega^2\int_\Omega\rho \xi\varphi_m^i \varphi_m^j + \int_\Omega\partial_n\varphi_m^i \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k \varphi_l^j + \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l \varphi_k^j\right)\\ diff --git a/examples/step-62/step-62.cc b/examples/step-62/step-62.cc index 551ac1b4e3..8d5e178a10 100644 --- a/examples/step-62/step-62.cc +++ b/examples/step-62/step-62.cc @@ -378,7 +378,7 @@ namespace step62 // This function defines the spatial shape of the force vector pulse which // takes the form of a Gaussian function - // @f{align*} + // @f{align*}{ // F_x &= // \left\{ // \begin{array}{ll} diff --git a/examples/step-64/step-64.cc b/examples/step-64/step-64.cc index dccfe9f9a6..ae596a979f 100644 --- a/examples/step-64/step-64.cc +++ b/examples/step-64/step-64.cc @@ -152,7 +152,7 @@ namespace Step64 // The Helmholtz problem we want to solve here reads in weak form as follows: - // @f{eqnarray*} + // @f{eqnarray*}{ // (\nabla v, \nabla u)+ (v, a(\mathbf x) u) &=&(v,1) \quad \forall v. // @f} // If you have seen step-37, then it will be obvious that diff --git a/examples/step-66/doc/intro.dox b/examples/step-66/doc/intro.dox index b8516abd8b..5536ac0101 100644 --- a/examples/step-66/doc/intro.dox +++ b/examples/step-66/doc/intro.dox @@ -33,7 +33,7 @@ On the unit circle $\Omega = \bigl\{ x \in \mathbb{R}^2 : \|x\| \leq 1 \bigr\}$ we consider the following nonlinear elliptic boundary value problem subject to a homogeneous Dirichlet boundary condition: Find a function $u\colon\Omega\to\mathbb{R}$ such that it holds: -@f{align*} +@f{align*}{ - \Delta u &= \exp(u) & \quad & \text{in } \Omega,\\ u &= 0 & \quad & \text{on } \partial\Omega. @f} diff --git a/examples/step-68/doc/intro.dox b/examples/step-68/doc/intro.dox index 777a3994c0..5cf44cf2ee 100644 --- a/examples/step-68/doc/intro.dox +++ b/examples/step-68/doc/intro.dox @@ -194,7 +194,7 @@ The stream function $\Psi$ of this Rayleigh-Kothe vortex is defined as: \Psi = \frac{1}{\pi} \sin^2 (\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T} \right) @f] where $T$ is half the period of the flow. The velocity profile in 2D ($\textbf{u}=[u,v]^T$) is : -@f{eqnarray*} +@f{eqnarray*}{ u &=& - \frac{\partial\Psi}{\partial y} = -2 \sin^2 (\pi x) \sin (\pi y) \cos (\pi y) \cos \left( \pi \frac{t}{T} \right)\\ v &=& \frac{\partial\Psi}{\partial x} = 2 \cos(\pi x) \sin(\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T} \right) @f} diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index b6fd8d6efa..8d2e6c687b 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -56,7 +56,7 @@ the second-order scheme discussed in @cite GuermondEtAl2018. The compressible Euler's equations of gas dynamics are written in conservative form as follows: -@f{align} +@f{align}{ \mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0} , @f} where $\mathbf{u}(\textbf{x},t):\mathbb{R}^{d} \times \mathbb{R} @@ -69,7 +69,7 @@ $\textbf{u} = [\rho, \textbf{m}^\top,E]^{\top}$: where $\rho \in \mathbb{R}^+$ denotes the density, $\textbf{m} \in \mathbb{R}^d$ is the momentum, and $E \in \mathbb{R}^+$ is the total energy of the system. The flux of the system $\mathbb{f}(\mathbf{u})$ is defined as -@f{align*} +@f{align*}{ \mathbb{f}(\textbf{u}) = \begin{bmatrix} @@ -83,7 +83,7 @@ $\otimes$ denotes the tensor product. Here, we have introduced the pressure $p$ that, in general, is defined by a closed-form equation of state. In this tutorial we limit the discussion to the class of polytropic ideal gases for which the pressure is given by -@f{align*} +@f{align*}{ p = p(\textbf{u}) := (\gamma -1) \Big(E - \tfrac{|\textbf{m}|^2}{2\,\rho} \Big), @@ -96,7 +96,7 @@ heats.

    Solution theory

    Hyperbolic conservation laws, such as -@f{align*} +@f{align*}{ \mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}, @f} pose a significant challenge with respect to solution theory. An evident @@ -112,7 +112,7 @@ understanding of hyperbolic conservation laws was to assume that the solution is formally defined as $\mathbf{u} := \lim_{\epsilon \rightarrow 0^+} \mathbf{u}^{\epsilon}$ where $\mathbf{u}^{\epsilon}$ is the solution of the parabolic regularization -@f{align} +@f{align}{ \mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon}) - {\epsilon} \Delta \mathbf{u}^{\epsilon} = 0. @f} @@ -127,7 +127,7 @@ Global existence and uniqueness of such solutions is an open issue. However, we know at least that if such viscosity solutions exists they have to satisfy the constraint $\textbf{u}(\mathbf{x},t) \in \mathcal{B}$ for all $\mathbf{x} \in \Omega$ and $t \geq 0$ where -@f{align} +@f{align}{ \mathcal{B} = \big\{ \textbf{u} = [\rho, \textbf{m}^\top,E]^{\top} \in \mathbb{R}^{d+2} \, \big | \ @@ -139,7 +139,7 @@ all $\mathbf{x} \in \Omega$ and $t \geq 0$ where \big\}. @f} Here, $s(\mathbf{u})$ denotes the specific entropy -@f{align} +@f{align}{ s(\mathbf{u}) = \ln \Big(\frac{p(\mathbf{u})}{\rho^{\gamma}}\Big). @f} We will refer to $\mathcal{B}$ as the invariant set of Euler's equations. @@ -166,7 +166,7 @@ $\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$. Following Step-9, Step-12, Step-33, and Step-67, at this point it might look tempting to base a discretization of Euler's equations on a (semi-discrete) variational formulation: -@f{align*} +@f{align*}{ (\partial_t\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} - ( \mathbb{f}(\mathbf{u}_{h}) ,\text{grad} \, \textbf{v}_{h})_{L^2(\Omega)} + s_h(\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} = \boldsymbol{0} @@ -181,7 +181,7 @@ are based on such a (semi-discrete) variational approach. Fundamentally, from an analysis perspective, variational discretizations are conceived to provide some notion of global (integral) stability, meaning an estimate of the form -@f{align*} +@f{align*}{ |\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!| @f} holds true, where $|\!|\!| \cdot |\!|\!| $ could represent the @@ -211,7 +211,7 @@ have an edge in many regards. In this tutorial step we therefore depart from variational schemes. We will present a completely algebraic formulation (with the flavor of a collocation-type scheme) that preserves constraints pointwise, i.e., -@f{align*} +@f{align*}{ \textbf{u}_h(\mathbf{x}_i,t) \in \mathcal{B} \;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}. @f} @@ -320,7 +320,7 @@ Consider the following pseudo-code, illustrating a possible straight forward strategy for computing the solution $\textbf{U}^{n+1}$ at a new time $t_{n+1} = t_n + \tau_n$ given a known state $\textbf{U}^{n}$ at time $t_n$: -@f{align*} +@f{align*}{ &\textbf{for } i \in \mathcal{V} \\ &\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow \mathtt{gather\_cij\_vectors} (\textbf{c}, \mathcal{I}(i)) \\ @@ -358,7 +358,7 @@ patch). The actual implementation will deviate from above code in one key aspect: the time-step size $\tau$ has to be chosen subject to a CFL condition -@f{align*} +@f{align*}{ \tau_n = c_{\text{cfl}}\,\min_{ i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right), @f} @@ -414,7 +414,7 @@ For the case of the reflecting boundary conditions we will proceed as follows: enforce reflecting boundary conditions. At the end of the time step we enforce reflecting boundary conditions strongly in a post-processing step where we execute the projection - @f{align*} + @f{align*}{ \mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\widehat{\boldsymbol{\nu}}_i \cdot \mathbf{m}_i) \widehat{\boldsymbol{\nu}}_i \ \ \text{where} \ \ @@ -458,7 +458,7 @@ If $\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}$ represents Euler's equation with reflecting boundary conditions on the entirety of the boundary (i.e. $\partial\Omega^r \equiv \partial\Omega$) and we integrate in space and time $\int_{\Omega}\int_{t_1}^{t_2}$ we would obtain -@f{align*} +@f{align*}{ \int_{\Omega} \rho(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} = \int_{\Omega} \rho(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \ \int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} @@ -477,7 +477,7 @@ the domain, we would like to know that our implementation of reflecting boundary conditions is consistent with the conservation properties mentioned above. In particular, if we use the projection $\boldsymbol{(1)}$ in the entirety of the domain the following discrete mass-balance can be guaranteed: -@f{align*} +@f{align*}{ \sum_{i \in \mathcal{V}} m_i \rho_i^{n+1} = \sum_{i \in \mathcal{V}} m_i \rho_i^{n} \ , \ \ \sum_{i \in \mathcal{V}} m_i \mathbf{m}_i^{n+1} diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index b89f9db87b..8d282059a5 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -1008,7 +1008,7 @@ namespace Step69 // Finally, assuming that $\mathbf{x}_i$ is a support point at the boundary, // the (nodal) normals are defined as: // - // @f{align*} + // @f{align*}{ // \widehat{\boldsymbol{\nu}}_i \dealcoloneq // \frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \, // \, \mathrm{d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i @@ -1474,7 +1474,7 @@ namespace Step69 // Next, we need two local wavenumbers that are defined in terms of a // primitive state $[\rho, u, p, a]$ and a given pressure $p^\ast$ // @cite GuermondPopov2016 Eqn. (3.7): - // @f{align*} + // @f{align*}{ // \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} // \left(\frac{p^\ast-p}{p}\right)_+} // @f} @@ -1498,7 +1498,7 @@ namespace Step69 } // Analougously @cite GuermondPopov2016 Eqn. (3.8): - // @f{align*} + // @f{align*}{ // \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} // \left(\frac{p^\ast-p}{p}\right)_+} // @f} @@ -1557,7 +1557,7 @@ namespace Step69 // general, not as sharp as the two-rarefaction estimate. But it will // save the day in the context of near vacuum conditions when the // two-rarefaction approximation might attain extreme values: - // @f{align*} + // @f{align*}{ // \lambda_{\text{exp}} = \max(u_i,u_j) + 5. \max(a_i, a_j). // @f} // @note The constant 5.0 multiplying the maximum of the sound speeds diff --git a/examples/step-7/doc/intro.dox b/examples/step-7/doc/intro.dox index f70cf04e0e..494069fa42 100644 --- a/examples/step-7/doc/intro.dox +++ b/examples/step-7/doc/intro.dox @@ -127,7 +127,7 @@ boundary term ${(v,g_2)}_{\Gamma_2}$ has appeared by integration by parts and using $\partial_n u=g_2$ on $\Gamma_2$ and $v=0$ on $\Gamma_1$. The cell matrices and vectors which we use to build the global matrices and right hand side vectors in the discrete formulation therefore look like this: -@f{eqnarray*} +@f{eqnarray*}{ A_{ij}^K &=& \left(\nabla \varphi_i, \nabla \varphi_j\right)_K +\left(\varphi_i, \varphi_j\right)_K, \\ @@ -167,7 +167,7 @@ where the centers $x_i$ of the exponentials are $\mathbf x_3=(\frac 12,-\frac 12)$, and the half width is set to $\sigma=\frac {1}{8}$. The method of manufactured solution then says: choose -@f{align*} +@f{align*}{ f &= -\Delta \bar u + \bar u, \\ g_1 &= \bar u|_{\Gamma_1}, \\ g_2 &= {\mathbf n}\cdot \nabla\bar u|_{\Gamma_2}. diff --git a/examples/step-70/doc/intro.dox b/examples/step-70/doc/intro.dox index a7cacbed04..ca5be5c677 100644 --- a/examples/step-70/doc/intro.dox +++ b/examples/step-70/doc/intro.dox @@ -108,7 +108,7 @@ one or equal with respect to the dimension of the embedding domain $\Omega$ We are going to solve the following differential problem: given a sufficiently regular function $g$ on $\Gamma$, find the solution $(\textbf{u},p)$ to -@f{eqnarray*} +@f{eqnarray*}{ -\Delta \mathbf{u} + \nabla p &=& 0,\\ -\nabla \cdot \textbf{u} &=& 0,\\ \textbf{u} &=& \textbf{g} \text{ in } \Gamma,\\ @@ -131,7 +131,7 @@ on $\partial\Omega$. The weak form of the Stokes equations is obtained by first writing it in vector form as -@f{eqnarray*} +@f{eqnarray*}{ \begin{pmatrix} {-\Delta \textbf{u} + \nabla p} \\ @@ -147,7 +147,7 @@ form as forming the dot product from the left with a vector-valued test function $\phi = \begin{pmatrix}\textbf{v} \\ q\end{pmatrix}$, and integrating over the domain $\Omega$, yielding the following set of equations: -@f{eqnarray*} +@f{eqnarray*}{ (\mathrm v, -\Delta \textbf{u} + \nabla p)_{\Omega} - @@ -314,7 +314,7 @@ $\Omega$. We can think of this as a thick object moving around in the fluid. In the case of $\mathcal{L}^2$ penalization, the additional penalization term can be interpreted as a Darcy term within $\Gamma$, resulting in: -@f{eqnarray*} +@f{eqnarray*}{ (\nabla \textbf{v}, \nabla \textbf{u})_{\Omega} - & (\textrm{div}\; \textbf{v}, p)_{\Omega} - (q, \textrm{div}\; \textbf{u})_{\Omega} + \beta (\textbf{v},\textbf{u})_{\Gamma} = \beta (\textbf{v},\textbf{g})_{\Gamma}. diff --git a/examples/step-71/doc/intro.dox b/examples/step-71/doc/intro.dox index d4a2aba9c3..49a8a7a64a 100644 --- a/examples/step-71/doc/intro.dox +++ b/examples/step-71/doc/intro.dox @@ -339,7 +339,7 @@ if one of the input variables is a function of another, it is also held constant and the chain rule does not propagate any further, while the computing total derivative would imply judicious use of the chain rule. This can be better understood by comparing the following two statements: -@f{align*} +@f{align*}{ \frac{\partial f\left(x, y\left(x\right)\right)}{\partial x} &= \frac{d f\left(x, y\left(x\right)\right)}{d x} \Big\vert_{y} \\ \frac{d f\left(x, y\left(x\right)\right)}{d x} @@ -636,7 +636,7 @@ Under the assumptions that an incompressible medium is being tested, and that the deformation profile through the sample thickness is linear, then the displacement at some measurement point $\mathbf{X}$ within the sample, expressed in radial coordinates, is -@f{align*} +@f{align*}{ r(\mathbf{X}) &= \frac{R(X_{1}, X_{2})}{\sqrt{\lambda_{3}}} , \\ \theta(\mathbf{X}) diff --git a/examples/step-71/step-71.cc b/examples/step-71/step-71.cc index 9c8310e330..1b3f7bb36d 100644 --- a/examples/step-71/step-71.cc +++ b/examples/step-71/step-71.cc @@ -2001,7 +2001,7 @@ namespace Step71 // \boldsymbol{\mathbb{H}} // \right] // @f] - // @f{align} + // @f{align}{ // \mathbf{S}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} // \right) // \dealcoloneq 2 \frac{d \psi_{0} \left( \mathbf{C}, @@ -2083,7 +2083,7 @@ namespace Step71 // \right)}{d \boldsymbol{\mathbb{H}} \otimes d \boldsymbol{\mathbb{H}}} // + \mu_{0} \mu_{r} \text{det}(\mathbf{F}) \mathbf{C}^{-1} // @f] - // @f{align} + // @f{align}{ // \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} // \right) = - \frac{d \mathbf{S}^{\text{tot}}}{d \boldsymbol{\mathbb{H}}} // &= - \mu_{e} @@ -2114,7 +2114,7 @@ namespace Step71 // \cdot \boldsymbol{\mathbb{H}} // \right]}{d \mathbf{C} \otimes \mathbf{C} \boldsymbol{\mathbb{H}}} // @f} - // @f{align} + // @f{align}{ // \mathcal{H}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} // \right) = 2 \frac{d \mathbf{S}^{\text{tot}}}{d \mathbf{C}} // &= 2 \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) @@ -2198,7 +2198,7 @@ namespace Step71 // = - C^{-1}_{ac} C^{-1}_{be} \mathbb{H}_{e} - C^{-1}_{ae} \mathbb{H}_{e} // C^{-1}_{bc} // @f] - // @f{align} + // @f{align}{ // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot // \boldsymbol{\mathbb{H}}\right]}{d \mathbf{C} \otimes d \mathbf{C}} // &= -\frac{d \left[\left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} @@ -2712,7 +2712,7 @@ namespace Step71 // \right] \otimes \frac{d f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} // \right)}{d \boldsymbol{\mathbb{H}}} // @f] - // @f{align} + // @f{align}{ // \mathcal{H}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, // \boldsymbol{\mathbb{H}} \right) // &= 2 \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right) diff --git a/examples/step-72/doc/intro.dox b/examples/step-72/doc/intro.dox index 49009b1c25..354946d227 100644 --- a/examples/step-72/doc/intro.dox +++ b/examples/step-72/doc/intro.dox @@ -17,7 +17,7 @@ Award No. EAR-1550901 and The University of California-Davis. This program solves the same problem as step-15, that is, it solves for the [minimal surface equation](https://en.wikipedia.org/wiki/Minimal_surface) - @f{align*} + @f{align*}{ F(u) \dealcoloneq -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right) &= 0 \qquad \qquad &&\textrm{in} ~ \Omega \\ @@ -70,11 +70,11 @@ might actually be enough. Our goal then is this: When using a Newton iteration, we need to repeatedly solve the linear partial differential equation - @f{align*} + @f{align*}{ F'(u^{n},\delta u^{n}) &=- F(u^{n}) @f} so that we can compute the update - @f{align*} + @f{align*}{ u^{n+1}&=u^{n}+\alpha^n \delta u^{n} @f} with the solution $\delta u^{n}$ of the Newton step. As discussed in step-15, diff --git a/examples/step-74/doc/intro.dox b/examples/step-74/doc/intro.dox index 7957c1e964..127e1b61a5 100644 --- a/examples/step-74/doc/intro.dox +++ b/examples/step-74/doc/intro.dox @@ -56,7 +56,7 @@ respectively. Note that when $f\subset \partial \Omega$, we define $\jump{v} = v $\average{v}=v$. The discretization using the SIPG is given by the following weak formula (more details can be found in @cite di2011mathematical and the references therein) -@f{align*} +@f{align*}{ &\sum_{K\in {\mathbb T}_h} (\nabla v_h, \nu \nabla u_h)_K\\ &-\sum_{F \in F_h^i} \left\{ \left< \jump{v_h}, \nu\average{ \nabla u_h} \cdot \mathbf n \right>_F diff --git a/examples/step-79/doc/intro.dox b/examples/step-79/doc/intro.dox index 2d737e5bf3..8f3eedf6ae 100644 --- a/examples/step-79/doc/intro.dox +++ b/examples/step-79/doc/intro.dox @@ -125,7 +125,7 @@ formulations). Furthermore, we will make the assumption that the material is linear isotropic, in which case the stress-strain tensor can be expressed in terms of the Lamé parameters $\lambda,\mu$ such that -@f{align} +@f{align}{ \boldsymbol{\sigma} &= \rho^p (\lambda \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2 \mu \boldsymbol{\varepsilon}) , \\ \sigma_{i,j} &= \rho^p (\lambda \varepsilon_{k,k} \delta_{i,j} + 2 \mu \varepsilon_{i,j}) . @f} diff --git a/examples/step-8/doc/intro.dox b/examples/step-8/doc/intro.dox index e454fc4ef5..855f087c00 100644 --- a/examples/step-8/doc/intro.dox +++ b/examples/step-8/doc/intro.dox @@ -190,7 +190,7 @@ respectively. For example (though this sequence of shape functions is not guaranteed, and you should not rely on it), the following layout could be used by the library: -@f{eqnarray*} +@f{eqnarray*}{ \Phi_0({\mathbf x}) &=& \left(\begin{array}{c} \varphi_0({\mathbf x}) \\ 0 @@ -250,7 +250,7 @@ $U_i$ such that If we insert the definition of the bilinear form and the representation of ${\mathbf u}_h$ and ${\mathbf v}_h$ into this formula: -@f{eqnarray*} +@f{eqnarray*}{ \sum_{i,j} U_i V_j \sum_{k,l} @@ -307,7 +307,7 @@ $\Phi_i$, which are of course given as follows (see their definition): @f] with the Kronecker symbol $\delta_{nm}$. Due to this, we can delete some of the sums over $k$ and $l$: -@f{eqnarray*} +@f{eqnarray*}{ A^K_{ij} &=& \sum_{k,l} @@ -366,7 +366,7 @@ the sums over $k$ and $l$: @f} Likewise, the contribution of cell $K$ to the right hand side vector is -@f{eqnarray*} +@f{eqnarray*}{ f^K_j &=& \sum_l diff --git a/examples/step-81/doc/intro.dox b/examples/step-81/doc/intro.dox index 4c46d4bd7e..fb314effa0 100644 --- a/examples/step-81/doc/intro.dox +++ b/examples/step-81/doc/intro.dox @@ -47,7 +47,7 @@ is described by Maxwell's equations @cite Schwartz1972, @cite Monk2003 : -@f{align*} +@f{align*}{ \frac{\partial}{\partial t} \mathbf{H} + \nabla \times \mathbf{E} &= -\mathbf{M}_a, \\ \nabla \cdot \mathbf{H} &= \rho_m, @@ -87,7 +87,7 @@ $\tilde{\mathbf{F}}(\mathbf{x})$ is a corresponding complex-valued vector field (or density). Inserting this ansatz into Maxwell's equations, substituting the charge conservation equations and some minor algebra then yields the so-called time-harmonic Maxwell's equations: -@f{align*} +@f{align*}{ -i\omega \tilde{\mathbf{H}} + \nabla \times \tilde{\mathbf{E}} &= -\tilde{\mathbf{M}}_a, \\ @@ -144,7 +144,7 @@ and then taking the limit of the upper and lower part of the line integral approaching the sheet. In contrast, the tangential part of the electric field is continuous. By fixing a unit normal $\mathbf{\nu}$ on the hypersurface $\Sigma$ both jump conditions are -@f{align*} +@f{align*}{ \mathbf{\nu} \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma} &= \sigma^{\Sigma}\left[(\mathbf{\nu}\times \mathbf{E}\times \mathbf{\nu})\right]|_{\Sigma}, \\ @@ -201,7 +201,7 @@ Finally, the interface conductivity is rescaled as @f] Accordingly, our rescaled equations are -@f{align*} +@f{align*}{ -i\mu_r \hat{\mathbf{H}} + \hat{\nabla} \times \hat{\mathbf{E}} &= -\hat{\mathbf{M}}_a, \\ @@ -422,7 +422,7 @@ in which $r = e_r \cdot x$ and $s(\tau)$ is an appropriately chosen, nonnegative scaling function. We introduce the following $2\times2$ matrices -@f{align*} +@f{align*}{ A &= T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2}, \frac{1}{d\bar{d}}\right)T_{e_xe_r} \\ @@ -434,14 +434,14 @@ We introduce the following $2\times2$ matrices in which -@f{align*} +@f{align*}{ d &= 1 + is(r) \\ \bar{d} &= 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau @f} and $T_{e_xe_r}$ is the rotation matrix which rotates $e_r$ onto $e_x$. Thus, after applying the rescaling, we get the following modified parameters -@f{align*} +@f{align*}{ \bar{\mu}_r^{-1} &= \frac{\mu_r^{-1}}{d}, \\ \bar{\varepsilon}_r &= A^{-1} \varepsilon_r B^{-1}, \text{ and } diff --git a/examples/step-85/doc/intro.dox b/examples/step-85/doc/intro.dox index cc8d17a5ae..bdc5acf9b8 100644 --- a/examples/step-85/doc/intro.dox +++ b/examples/step-85/doc/intro.dox @@ -14,7 +14,7 @@ under grants 2014-6088 (Kreiss) and 2017-05038 (Massing). In this example, we show how to use the cut finite element method (CutFEM) in deal.II. For illustration, we want to solve the simplest possible problem, so we again consider Poisson's equation: -@f{align*} +@f{align*}{ -\Delta u &= f \qquad && \text{in }\, \Omega, \\ u &= u_D \qquad && \text{on }\, \Gamma = \partial \Omega, @@ -29,12 +29,12 @@ Since we no longer use the mesh to describe the geometry of the domain, we need some other way to represent it. This can be done in several ways but here we assume that $\Omega$ is described by a level set function, $\psi : \mathbb{R}^{\text{dim}} \to \mathbb{R}$ such that -@f{align*} +@f{align*}{ \Omega &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) < 0 \}, \\ \Gamma &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) = 0 \}. @f} For simplicity, we choose $\Omega$ to be a unit disk, so that -@f{equation*} +@f{equation*}{ \psi(x) = \| x \| - 1. @f} As can be seen from the figure below, @@ -44,17 +44,17 @@ zero on the boundary, and positive everywhere else. To solve this problem, we want to distribute degrees of freedom over the smallest submesh, $\mathcal{T}_\Omega^h$, that completely covers the domain: -@f{equation*} +@f{equation*}{ \mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Omega \neq \emptyset \}. @f} This is usually referred to as the "active mesh". @image html step-85-active-mesh.svg The finite element space where we want to find our numerical solution, $u_h$, is now -@f{equation*} +@f{equation*}{ V_\Omega^h = \{ v \in C(\mathcal{N}_\Omega^h) : v \in Q_p(T), \, T \in \mathcal{T}_\Omega^h \}, @f} where -@f{equation*} +@f{equation*}{ \mathcal{N}_\Omega^h = \bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T}, @f} and $\overline{T}$ denotes the closure of $T$. @@ -66,23 +66,23 @@ In this type of immersed finite element method, the standard way to apply boundary conditions is using Nitsche's method. Multiplying the PDE with a test function, $v_h \in V_\Omega^h$, and integrating by parts over $\Omega$, as usual, gives us -@f{equation*} +@f{equation*}{ (\nabla u_h, \nabla v_h)_\Omega - (\partial_n u_h, v_h)_\Gamma = (f,v)_\Omega. @f} Let $\gamma_D > 0$ be a scalar penalty parameter and let $h$ be some measure of the local cell size. We now note that the following terms are consistent with the Dirichlet boundary condition: -@f{align*} +@f{align*}{ -(u_h, \partial_n v_h)_\Gamma &= -(u_D, \partial_n v_h)_\Gamma, \\ \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma &= \left (\frac{\gamma_D}{h}u_D, v_h \right )_\Gamma. @f} Thus, we can add these to the weak formulation to enforce the boundary condition. This leads to the following weak formulation: Find $u_h \in V_\Omega^h$ such that -@f{equation*} +@f{equation*}{ a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h, @f} where -@f{align*} +@f{align*}{ a_h(u_h, v_h) &= (\nabla u_h, \nabla v_h)_\Omega - (\partial_n u_h, v_h)_\Gamma - (u_h, \partial_n v_h)_\Gamma @@ -101,7 +101,7 @@ that is, over $T \cap \Omega$ and $T \cap \Gamma$. @image html immersed_quadratures.svg Since $\Omega \cap T$ is the part of the cell that lies inside the domain, we shall refer to the following regions -@f{align*} +@f{align*}{ \{x \in T : \psi(x) < 0 \}, \\ \{x \in T : \psi(x) > 0 \}, \\ \{x \in T : \psi(x) = 0 \}, @@ -126,11 +126,11 @@ One way to avoid this problem is to add a so-called ghost penalty term, $g_h$, to the weak formulation (see e.g. @cite burman_hansbo_2012 and @cite cutfem_2015). This leads to the stabilized cut finite element method, which reads: Find $u_h \in V_\Omega^h$ such that -@f{equation*} +@f{equation*}{ A_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h, @f} where -@f{equation*} +@f{equation*}{ A_h(u_h,v_h) = a_h(u_h,v_h) + g_h(u_h, v_h). @f} The point of this ghost penalty is that it makes the numerical method essentially independent @@ -138,11 +138,11 @@ of how $\Omega$ relates to the background mesh. In particular, $A_h$ can be shown to be continuous and coercive, with constants that do not depend on how $\Omega$ intersects $\mathcal{T}^h$. To define the ghost penalty, let $\mathcal{T}_\Gamma^h$ be the set of intersected cells: -@f{equation*} +@f{equation*}{ \mathcal{T}_{\Gamma}^h = \{ T \in \mathcal{T}_{\Omega}^{h} : T \cap \Gamma \neq \emptyset \}, @f} and let $\mathcal{F}_h$ denote the interior faces of the intersected cells in the active mesh: -@f{equation*} +@f{equation*}{ \mathcal{F}_h = \{ F = \overline{T}_+ \cap \overline{T}_- : \, T_+ \in \mathcal{T}_{\Gamma}^h, \, T_- \in \mathcal{T}_{\Omega}^h @@ -150,11 +150,11 @@ and let $\mathcal{F}_h$ denote the interior faces of the intersected cells in th @f} @image html step-85-ghost-faces.svg The ghost penalty acts on these faces and reads -@f{equation*} +@f{equation*}{ g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} g_F(u_h, v_h), @f} where $g_F$ is the face-wise ghost penalty: -@f{equation*} +@f{equation*}{ g_F(u_h, v_h) = \gamma_A \sum_{k=0}^p \left(\frac{h_F^{2k-1}}{k!^2}[\partial_n^k u_h], [\partial_n^k v_h] \right)_F. @f} Here, $\gamma_A$ is a penalty parameter and $h_F$ is some measure of the face size. @@ -166,7 +166,7 @@ Hand-wavingly speaking, this is the reason why we obtain a cut-independent method when we enforce $g_F(u_h, v_h) = 0$ over the faces in $\mathcal{F}_h$. Here, we shall use a continuous space of $Q_1$-elements, so the ghost penalty is reduced to -@f{equation*} +@f{equation*}{ g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} (h_F [\partial_n u_h], [\partial_n v_h])_F. @f} @@ -176,7 +176,7 @@ such that the domain deforms with time. For such a problem, one would typically solve for an approximation of the level set function, $\psi_h \in V^h$, in a separate finite element space over the whole background mesh: -@f{equation*} +@f{equation*}{ V^h = \{ v \in C(\mathcal{N}^h) : v \in Q_p(T), \, T \in \mathcal{T}^h \}, @f} where $\mathcal{N}^h = \bigcup_{T \in \mathcal{T}^h} \overline{T}$. @@ -184,12 +184,12 @@ Even if we solve a much simpler problem with a stationary domain in this tutoria we shall, just to illustrate, still use a discrete level set function for the Poisson problem. Technically, this is a so-called "variational crime" because we are actually not using the bilinear form $a_h$ but instead -@f{equation*} +@f{equation*}{ a_h^\star(u_h, v_h) = (\nabla u_h, \nabla v_h)_{\Omega_h} - (\partial_n u_h, v_h)_{\Gamma_h} + \ldots @f} This is an approximation of $a_h$ since we integrate over the approximations of the geometry that we get via the discrete level set function: -@f{align*} +@f{align*}{ \Omega_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) < 0 \}, \\ \Gamma_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) = 0 \}. @f} diff --git a/examples/step-85/step-85.cc b/examples/step-85/step-85.cc index abcee95b07..82188c139b 100644 --- a/examples/step-85/step-85.cc +++ b/examples/step-85/step-85.cc @@ -427,7 +427,7 @@ namespace Step85 // etc., one can use its normal_vector(..)-function to get an outward // normal to the immersed surface, $\Gamma$. In terms of the level set // function, this normal reads - // @f{equation*} + // @f{equation*}{ // n = \frac{\nabla \psi}{\| \nabla \psi \|}. // @f} // An additional benefit of std::optional is that we do not need any @@ -573,7 +573,7 @@ namespace Step85 // To test that the implementation works as expected, we want to compute the // error in the solution in the $L^2$-norm. The analytical solution to the // Poisson problem stated in the introduction reads - // @f{align*} + // @f{align*}{ // u(x) = 1 - \frac{2}{\text{dim}}(\| x \|^2 - 1) , \qquad x \in // \overline{\Omega}. // @f} diff --git a/examples/step-9/doc/intro.dox b/examples/step-9/doc/intro.dox index ded3a8db95..e0f63bd7ad 100644 --- a/examples/step-9/doc/intro.dox +++ b/examples/step-9/doc/intro.dox @@ -303,7 +303,7 @@ for problems of the kind we have here. For the problem which we will solve in this tutorial program, we use the following domain and functions (in $d=2$ space dimensions): -@f{eqnarray*} +@f{eqnarray*}{ \Omega &=& [-1,1]^d \\ \beta({\mathbf x}) &=& @@ -436,7 +436,7 @@ criterion: \eta_K = h^{1+d/2} |\nabla_h u_h(K)|, @f] which is inspired by the following (not rigorous) argument: -@f{eqnarray*} +@f{eqnarray*}{ \|u-u_h\|^2_{L_2} &\le& C h^2 \|\nabla u\|^2_{L_2} diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index c0ba5835b3..413577eca0 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -807,7 +807,7 @@ extern "C" /* * Perform matrix sum: - * @f{equation*} + * @f{equation*}{ * C \dealcoloneq \beta C + \alpha op(A), * @f * where $op(A)$ denotes either $op(A) = A$ or $op(A)=A^T$. diff --git a/include/deal.II/numerics/smoothness_estimator.h b/include/deal.II/numerics/smoothness_estimator.h index 87d85c648f..44aaa043e0 100644 --- a/include/deal.II/numerics/smoothness_estimator.h +++ b/include/deal.II/numerics/smoothness_estimator.h @@ -74,7 +74,7 @@ namespace SmoothnessEstimator * * In one dimension, the finite element solution on cell $K$ with polynomial * degree $p$ can be written as - * @f{eqnarray*} + * @f{eqnarray*}{ * u_h(x) &=& \sum_j u_j \varphi_j (x) \\ * u_{h, k}(x) &=& \sum_{k=0}^{p} a_k \widetilde P_k (x), * \quad a_k = \sum_j {\cal L}_{k,j} u_j @@ -113,7 +113,7 @@ namespace SmoothnessEstimator * coefficients at once. If there are multiple coefficients corresponding to * the same absolute value of modes $\|{\bf k}\|_1$, we take the maximum * among those. Thus, the least-squares fit is performed on - * @f{eqnarray*} + * @f{eqnarray*}{ * \widetilde P_{\bf k}({\bf x}) &=& * \widetilde P_{k_1} (x_1) \ldots \widetilde P_{k_d} (x_d) \\ * \ln \left( \max\limits_{\|{\bf k}\|_1} |a_{\bf k}| \right) &\sim& @@ -260,7 +260,7 @@ namespace SmoothnessEstimator * From the definition, we can write our Fourier series expansion * $a_{\bf k}$ of the finite element solution on cell $K$ with polynomial * degree $p$ as a matrix product - * @f{eqnarray*} + * @f{eqnarray*}{ * u_h({\bf x}) &=& \sum_j u_j \varphi_j ({\bf x}) \\ * u_{h, {\bf k}}({\bf x}) &=& * \sum_{{\bf k}, \|{\bf k}\|\le p} a_{\bf k} \phi_{\bf k}({\bf x}), @@ -280,7 +280,7 @@ namespace SmoothnessEstimator * If the finite element approximation on cell $K$ is part of the Hilbert * space $H^s(K)$, then the following integral must exist for both the finite * element and spectral representation of our solution - * @f{eqnarray*} + * @f{eqnarray*}{ * \| \nabla^s u_h({\bf x}) \|_{L^2(K)}^2 &=& * \int\limits_K \left| \nabla^s u_h({\bf x}) \right|^2 d{\bf x} < * \infty \\ -- 2.39.5