]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use balanced braces for doxygen formulas. 15640/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Tue, 4 Jul 2023 02:47:23 +0000 (20:47 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Tue, 4 Jul 2023 02:47:23 +0000 (20:47 -0600)
70 files changed:
examples/step-11/doc/intro.dox
examples/step-14/doc/intro.dox
examples/step-15/doc/intro.dox
examples/step-15/doc/results.dox
examples/step-18/doc/intro.dox
examples/step-20/doc/intro.dox
examples/step-21/doc/intro.dox
examples/step-22/doc/intro.dox
examples/step-22/doc/results.dox
examples/step-23/doc/intro.dox
examples/step-24/doc/intro.dox
examples/step-25/doc/intro.dox
examples/step-26/doc/intro.dox
examples/step-26/doc/results.dox
examples/step-28/doc/intro.dox
examples/step-29/doc/intro.dox
examples/step-3/doc/intro.dox
examples/step-31/doc/intro.dox
examples/step-31/doc/results.dox
examples/step-32/doc/intro.dox
examples/step-32/doc/results.dox
examples/step-32/step-32.cc
examples/step-33/doc/intro.dox
examples/step-33/step-33.cc
examples/step-34/doc/intro.dox
examples/step-35/doc/intro.dox
examples/step-36/doc/intro.dox
examples/step-37/doc/intro.dox
examples/step-37/doc/results.dox
examples/step-38/doc/intro.dox
examples/step-39/doc/intro.dox
examples/step-40/doc/intro.dox
examples/step-41/doc/intro.dox
examples/step-42/doc/intro.dox
examples/step-43/doc/intro.dox
examples/step-44/doc/intro.dox
examples/step-44/step-44.cc
examples/step-45/doc/intro.dox
examples/step-45/step-45.cc
examples/step-46/doc/intro.dox
examples/step-5/doc/intro.dox
examples/step-50/doc/intro.dox
examples/step-51/doc/intro.dox
examples/step-52/doc/intro.dox
examples/step-55/doc/intro.dox
examples/step-56/doc/intro.dox
examples/step-57/doc/intro.dox
examples/step-59/doc/intro.dox
examples/step-61/doc/intro.dox
examples/step-62/doc/intro.dox
examples/step-62/step-62.cc
examples/step-64/step-64.cc
examples/step-66/doc/intro.dox
examples/step-68/doc/intro.dox
examples/step-69/doc/intro.dox
examples/step-69/step-69.cc
examples/step-7/doc/intro.dox
examples/step-70/doc/intro.dox
examples/step-71/doc/intro.dox
examples/step-71/step-71.cc
examples/step-72/doc/intro.dox
examples/step-74/doc/intro.dox
examples/step-79/doc/intro.dox
examples/step-8/doc/intro.dox
examples/step-81/doc/intro.dox
examples/step-85/doc/intro.dox
examples/step-85/step-85.cc
examples/step-9/doc/intro.dox
include/deal.II/lac/scalapack.templates.h
include/deal.II/numerics/smoothness_estimator.h

index ab72e1fcaa330444c9a0af82024cca82b9f19e00..d092c787d2f25ab37f1c9bb7633b9ee96ffadacc 100644 (file)
@@ -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.
index 308d8d4db2fdc69750d514e8daf58a0ed2584f30..93405873ae9f87f5984ec990eaa1cc8cb127a672 100644 (file)
@@ -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
index 109bd83f83d899cc19a57d4ab874f6e393f832b6..8347c9cb943b771d52b5c1c388ad4276811698a4 100644 (file)
@@ -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}
index 858f39b3fa2c2f9b7f347417b992185c29c07f16..7ac1adb84f93af3bf743feb3753b95ffdca44673 100644 (file)
@@ -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
index 91df0426b44935329ced4ada3046ce376cf28283..3a84a62c604583afba2980600052ac5c83fc1def 100644 (file)
@@ -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 <i>real</i> 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
 <a name="step_18.linear-system"></a>
-@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; i<dofs_per_cell; ++i)
   between symmetric tensors of even rank here.
 
   Assembling the local contributions
-  @f{eqnarray*}
+  @f{eqnarray*}{
       f^K_i &=&
       (\mathbf{f}, \varphi_i)_K -(\sigma^{n-1},\varepsilon(\varphi_i))_K
       \\
index de94225f50b65c4042b90822cc4a81032852f2c9..7cf1e41a3b18e9b0977fa3807caf0ec400e39d16 100644 (file)
@@ -8,7 +8,7 @@ particular Raviart-Thomas elements -- and using block matrices to define
 solvers, preconditioners, and nested versions of those that use the
 substructure of the system matrix. The equation we are going to solve is again
 the Poisson equation, though with a matrix-valued coefficient:
-@f{eqnarray*}
+@f{eqnarray*}{
   -\nabla \cdot K({\mathbf x}) \nabla p &=& f \qquad {\textrm{in}\ } \Omega, \\
   p &=& g \qquad {\textrm{on}\ }\partial\Omega.
 @f}
@@ -89,7 +89,7 @@ areas of high pressure to areas of low pressure (thus the negative sign).
 
 With this second variable, one then finds an alternative version of the
 Laplace equation, called the <i>mixed formulation</i>:
-@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 <i>all</i> 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 \\
index 96877da161774f17bc4ae3c1bde4304abef4ff3a..3068d8c028805adc09940439ba6e83234b6678d0 100644 (file)
@@ -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
 <i>fractional flow</i>.
 
 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: <i>Analysis of gas-cap or dissolved-gas
 reservoirs</i>, 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.
 
   <li>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\},
index ec07c2bd89732e19b4a8c6d02ea04a5ba7b8ec8e..4e6cb841e6035917d2024410e1b9f83f4f26ad36 100644 (file)
@@ -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:
 <li>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:
 
 <li>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:
     <code>dim</code>-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</a> in the results section of this program.
 <h4> A note on the structure of the linear system </h4>
 
 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,
index d9f1a109e6083b06c34a3a71506dada85475f9cb..b11108f4f3aeaa1a657b4ed9a873a920e7c1e42b 100644 (file)
@@ -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)
index 1310890811ad77a8977aaa636feae9f458a2144f..6c52ea052c862f0347bff253ee7a8de1e4f56552 100644 (file)
@@ -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 <i>velocity</i> 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,
index 9d4f12848bb7820d5abbcbab1fba3ff4185d19a6..1f8e2de6aa60387f183fdbbe9d1df5333110ad7d 100644 (file)
@@ -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 -
index db22da07a2327c469777b46e291d78d8eda51a3b..15d6b222d3ea626d51194a30f3020a63c9ad74dc 100644 (file)
@@ -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}).
index 0fbd2b74d4c5ee1c93b31db6b35b1badf7094b71..d7c7fbefd86dc16a0a8ab26124bafaf6925d3bad 100644 (file)
@@ -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:
 
   <li><i>Test functions from different meshes</i>: 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}
index db8010651ff0fe20a03adb29e1c6d95c13949f10..3b4a2c3961a170a054469d02ec2a294bc118c30a 100644 (file)
@@ -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,
 <a href="https://doi.org/10.2478/cmam-2010-0025">Schatz et. al.</a> 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,
   &
index 7745ee3d6bed4d6bac80494e33c33638eb998a4f..400c42699202d3c4c2299bc9b8916fdbd62c5a52 100644 (file)
@@ -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$.
 
 <li> Define the so-called <i>fission source</i> by
-  @f{eqnarray*}
+  @f{eqnarray*}{
     s_f^{(n-1)}(x)
     =
     \frac{1}{k_{\mathrm{eff}}^{(n-1)}}
@@ -219,7 +219,7 @@ like this:
   @f}
 
 <li> 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}
 
 <li> 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,
index 163ab0510e541409a1fa39e0f743388890b1aec3..43969c28de23ed5749a93e2739604dce9f32ebf5 100644 (file)
@@ -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, \\
index f61fc92dafbc8a5d15a2db653d02a88cb1ebbf2f..0bbafd7184f50c7026dd908f400db9a4510e4641 100644 (file)
@@ -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$ <i>from the left</i> (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 <i>right</i> 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
index cd3b7d74c60a4b0e877fdd604668987e0b503267..88bb3b561e26f2fd3282f8cf4e6dbfba874f8833 100644 (file)
@@ -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)
index e2d8a0307a52543864be94b4d6ee738dc6f8871a..a10f5e24882092c923a627e80652afece29c7e27 100644 (file)
@@ -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?
   <li>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
index 598378f45183db68ecc58d1940a60ba0d82a1d2e..635a49203c9bb7aea2fee713514dd6c2d8988c85 100644 (file)
@@ -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.
 <h3> The scaling of discretized equations </h3>
 
 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 &mdash; around 10 km &mdash; 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 <i>before</i> 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.</i> 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 <i>ASPECT</i> 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 <a href="#Results">results section
   below</a>. 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 <i>velocity
   difference</i>, 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],
index 93b1b06ceea6298244b2cd7948bbac6476d899bc..a989c9cd2e60de8e0d511315f501e9b9cc677c5b 100644 (file)
@@ -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
index e6f3c5af7796e9fb43d0a7d0f0cdd7414eb32a8a..04f45cd4836003e3d5a89ece88cd43d742e559dd 100644 (file)
@@ -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)
index 4f5b09cceceb26de668529596824d915589702e0..89ee0634c894e8e621fbca8bd3441ab975804c3f 100644 (file)
@@ -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 <i> numerical </i> 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}
index 15ac9725ddabb3168376cca0339dc66443c75a1f..d81985ff724132bf9148268f9e8f7812761a0a3b 100644 (file)
@@ -1809,7 +1809,7 @@ namespace Step33
     // vector_valued module). It will be represented by the variable
     // <code>component_i</code> 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
index 704b9b4bcf086605114143826a0b4e769642d540..304af2b7a2013310de94a0f9cbe7f1d7dd512b03 100644 (file)
@@ -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, <i>Marine Hydrodynamics</i>, 1977,
index 208960e8f71afe1f31fb7840d59967adad6331d5..49f7e804900ff7d69d6cedb59c5fe03730b8da3f 100644 (file)
@@ -12,7 +12,7 @@ University. Most of the work for this program is by him.
 <h3> Motivation </h3>
 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}
index 92816cefbe8b8d3b52e3075ce1e0c970a646f05e..47d66f07507c12af52c4b99ea2fee8e4e24ae414 100644 (file)
@@ -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}
index e50ea2c32b980588de28aa8c4517f59701443996..1d9859bac67ea74f17875c77f0a108e8a84755e7 100644 (file)
@@ -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 <i>A</i> 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,
 <i>A</i><sub>cell</sub> denotes the cell matrix associated with <i>A</i>.
 
 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 (<i>q</i>*dim+<i>d,i</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
index 776ca401cb290a88040decf571e667e69fa46f10..d09c6ffe97909380447304e875d7ddd757b4db6d 100644 (file)
@@ -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
index ecbf3ab1b93f012a4e1b9de8472af1f30ce52288..096ad856b7519bfd029f9aa87faa5efcf6c37783 100644 (file)
@@ -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)
index 924519948c0a922c952d120ca602b0339cb45157..5a8d582931f9740ff4fbdd348fcb26eea2b511ca 100644 (file)
@@ -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
index c26d797eb24ae31ac942666c28f4a7c4f1d6055d..cbd60bafd4594c1389e485aed4da367e97a10f4a 100644 (file)
@@ -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\{
index 04eb08e134b127ddecd9294d3f558822af6051a1..3470ad73e3f723a7b68bf6d44ab055a554f83371 100644 (file)
@@ -49,7 +49,7 @@ is active there.
 <h3>Classical formulation</h3>
 
 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).
 <h3>Derivation of the variational inequality</h3>
 
 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:
 
 <i>Find a function $u\in G$ with</i>
-@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:
 
 <i>Find $u\in V$ and $\lambda\in K$ such that</i>
-@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}
 <i>with</i>
-@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}
index 527516706dd3aaf9337155dfd34c781aa472bad3..d09e30d2706787b7ba8cb8e4c49667a2a0f584c6 100644 (file)
@@ -53,7 +53,7 @@ step-32.
 <h3>Classical formulation</h3>
 
 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.
 
  <li> 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}
 
  <li>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}
index 2ec5bde432dda039f70a2ed12436f897d4df2837..8d90311422c372e8139c112c64eef7bb0e5b4a9e 100644 (file)
@@ -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}
index edd216beb7fa6b4a38f991bdce83606b6fb8752a..46e1fe52a7c8fb66c5104225d8e4439a465c0024 100644 (file)
@@ -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}}
index 0e9fcba5c2d3b9297b2313038fee7917124840c7..842ce3389130a4f9c06acd94b9fbadbae92344a6 100644 (file)
@@ -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}
index 11c65727976812cb5dfa987ab17b9657f81d702f..2f7e3c293fdc277036e95a8ff1202590b3ee579a 100644 (file)
@@ -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}
index f7a44ecd1ffe986627242f26f4c49d08657476d2..2aabdaaf8dc33459b2c70346d11559fc7d094268 100644 (file)
@@ -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},
index b198c68b933d01821d92ffd7af2b8ccb16b6978e..b99feef93bd4ac0aa85ca8f90e605516d7021085 100644 (file)
@@ -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} &
index 9e02ee824992d910fb1a8c2ad4fc969ea3f58ba2..c8982a34f3037ad6573b211f3827325fc0f423df 100644 (file)
@@ -42,7 +42,7 @@ them are:
 </ul>
 
 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 <code>assemble_system</code> function follows
index 10c6c0dad8fd967c362ce37d9e0c34c00f03f709..e95f5d2a85967a4ef0fa285764f23c4dd671ab47 100644 (file)
@@ -55,7 +55,7 @@ strongly when the number of unknowns per processor becomes smaller than
 <h3>The testcase</h3>
 
 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
index 85b4bbce94fb4d1c3d831d34093617fd543f96b9..9b9466b6995e2770889316e3d931abf50bcd71a5 100644 (file)
@@ -71,7 +71,7 @@ the <i>hybridizable discontinuous Galerkin</i> (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
index 13c47faeb371b649c0e1eab1ca53f7b9cc446061..6097ff2a4c0a0a584d24a265ff52e6acbedea850 100644 (file)
@@ -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}
 
index 4877a22847707ee6f87639402f2fd5caa6c88ef4..32bc359021c7631d36fa88e7b0456c01cb32253d 100644 (file)
@@ -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.
 <h3>The solver and preconditioner</h3>
 
 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}
index ee002f95cc6fdb51c30f89c6d17599908d3a870f..d1a6f70f00e6051e27dd01beb62c536b0172713d 100644 (file)
@@ -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
index ea54a8fe1790b0f9b7ba1b2e529ef287576be0c3..8d42d531534211e3a003fbfd892318c87992e7ee 100644 (file)
@@ -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.
 <h3> Linearization of Navier-Stokes Equations </h3>
 
 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}
index f78b7930f5020c7b5ffcbb76cfc7174e0cd2a8ed..6bfb94959c2a380909a9c889a9a728a9d3af084d 100644 (file)
@@ -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},
index 04ab0a81f9f8065d886ece74185b86ba8a88f1d7..dc2f339add91a3e3e7607f0e725a8558b051b4dc 100644 (file)
@@ -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 <code>cell_matrix_E</code>.
 
 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}
index e0164ca8fdfd64e8fb85586fb4564dfb1f240ae7..b2fb90ec0b4f8a9d9db30613cbdd3ad09d5a8447 100644 (file)
@@ -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)\\
index 551ac1b4e39b310c766a450f40226f5501ae2601..8d5e178a10df9a27ff805ac51c0a93e4d001c9a0 100644 (file)
@@ -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}
index dccfe9f9a6736b5a572dd300eefa43901387a2e2..ae596a979f985f07c6e62b38aaa1007012a775ef 100644 (file)
@@ -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
index b8516abd8b3ff722f6574fc8ac415bc90639c516..5536ac010196ee1731fe683664c75943bfeb2245 100644 (file)
@@ -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}
index 777a3994c03c7f0f2f0ca1035703fd41d15248c0..5cf44cf2ee2b84d76eee49917fa925dd0912f02c 100644 (file)
@@ -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}
index b6fd8d6efad5dba367b5b88a4e5ec7fdbef44520..8d2e6c687b1be366fc1253f1d85d544fb5ed1381 100644 (file)
@@ -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</a>.
 <h4>Solution theory</h4>
 
 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}
index b89f9db87b53639bfc6551822b2b255115326dd5..8d282059a51646a540e356d6525b2984894b8f85 100644 (file)
@@ -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
index f70cf04e0ea2f2525722b173aa9db08772582813..494069fa426cd3b2e4d2e24a398e938f8b64478a 100644 (file)
@@ -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}.
index a7cacbed045f16047d685661c74b45371912b6de..ca5be5c677e447be6b0cd5e0361784d60c55797f 100644 (file)
@@ -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}.
index d4a2aba9c3a3f3a6be6788f7a6f0d39678c3a333..49a8a7a64ac6f8354938f70f303463677948fa20 100644 (file)
@@ -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})
index 9c8310e3305e3fa0bae3e70eb3c1a726f4b7f45d..1b3f7bb36dd54e9a8f472f9e1d31539dcc9201e8 100644 (file)
@@ -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)
index 49009b1c25d014b72fb4402cc4da3e9df8753b0d..354946d227a1841217c12f64a97fbf1a018c88b4 100644 (file)
@@ -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,
index 7957c1e964824139defaa6bbb9ddb3c66d8ed9b8..127e1b61a568a1945cd2e1f2d09457bd32544844 100644 (file)
@@ -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
index 2d737e5bf3eea8265ee3eb4b344a493cb331e3dc..8f3eedf6aed3557630f02971505f7319a91ccf68 100644 (file)
@@ -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&eacute;
 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}
index e454fc4ef57bf23d0d71908cb8e2c6a353ccd8b0..855f087c0064d625ec6f6b5af8562aaa2b0946b3 100644 (file)
@@ -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
index 4c46d4bd7e717b1837a6de3f71cda0e77c9d8d70..fb314effa0d0d7f9f814ebc5d9bb9f20f495a690 100644 (file)
@@ -47,7 +47,7 @@ is described by
 <a href="https://en.wikipedia.org/wiki/Maxwell%27s_equations">Maxwell's
 equations</a>
 @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 <em>time-harmonic</em> 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 }
index cc8d17a5ae98a44e048f23b950ec7274fbd02148..bdc5acf9b835f0d8a1f6fa6f39be3d33178a4fb7 100644 (file)
@@ -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}
index abcee95b07d7f54a3312541eebedf7208ff1ae04..82188c139bd9771b5229352b489c48d85bb392aa 100644 (file)
@@ -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}
index ded3a8db95f22e6ca5315bd2cc5f517e23b4b0df..e0f63bd7adba9b94072332e7971de7a2b0b3ed6e 100644 (file)
@@ -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}
index c0ba5835b3b34df4a9cce4260d9b8a5d1722413a..413577eca019b0aae3b4f4506f41464b7418630f 100644 (file)
@@ -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$.
index 87d85c648f8bdd54fe5bbb7e02c5888dee9d569c..44aaa043e00f22ca1f7c577243f7d5952ff7c5cb 100644 (file)
@@ -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 \\

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.