From: Timo Heister Date: Sat, 6 Feb 2016 23:59:42 +0000 (-0500) Subject: fix block formulas X-Git-Tag: v8.4.0-rc2~27^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=581715734f012af2332116f9fcd504a9daf1fb49;p=dealii.git fix block formulas --- diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h index 66942a3d71..3eaef2cf13 100644 --- a/include/deal.II/base/derivative_form.h +++ b/include/deal.II/base/derivative_form.h @@ -32,7 +32,7 @@ DEAL_II_NAMESPACE_OPEN * with $\mathbf x\in {\mathbb R}^{\text{dim}}$, in such a way that the * directional derivative in direction $\mathbf d\in {\mathbb R}^{\text{dim}}$ * so that - * @f{align*} + * @f{align*}{ * \nabla f(\mathbf x) \mathbf d * = \lim_{\varepsilon\rightarrow 0} * \frac{f(\mathbf x + \varepsilon \mathbf d) - f(\mathbf x)}{\varepsilon}, diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index fa03f4b3d6..36048ac800 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -560,9 +560,15 @@ namespace Polynomials * degree four with vanishing values and derivatives at x=0 and * x=1, then by the product of this fourth order polynomial with * Legendre polynomials of increasing order. The implementation is - * @f{align*}{ p_0(x) &= 2x^3-3x^2+1 \\ p_1(x) &= -2x^2+3x^2 \\ p_2(x) &= - * x^3-2x^2+x \\ p_3(x) &= x^3-x^2 \\ p_4(x) &= 16x^2(x-1)^2 \\ \ldots & - * \ldots \\ p_k(x) &= x^2(x-1)^2 L_{k-4}(x) @f} + * @f{align*}{ + * p_0(x) &= 2x^3-3x^2+1 \\ + * p_1(x) &= -2x^2+3x^2 \\ + * p_2(x) &= x^3-2x^2+x \\ + * p_3(x) &= x^3-x^2 \\ + * p_4(x) &= 16x^2(x-1)^2 \\ + * \ldots & \ldots \\ + * p_k(x) &= x^2(x-1)^2 L_{k-4}(x) + * @f} * * @author Guido Kanschat * @date 2012 diff --git a/include/deal.II/base/polynomials_bdm.h b/include/deal.II/base/polynomials_bdm.h index 99582d8971..9af1b22a02 100644 --- a/include/deal.II/base/polynomials_bdm.h +++ b/include/deal.II/base/polynomials_bdm.h @@ -54,7 +54,7 @@ DEAL_II_NAMESPACE_OPEN * \begin{pmatrix} f_{y}(x,y) \\ -f_{x}(x,y) \end{pmatrix}$. * * The basis used to construct the $BDM_{1}$ shape functions is - * @f{align*} + * @f{align*}{ * \phi_0 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, * \phi_1 = \begin{pmatrix} -\sqrt{3}+2\sqrt{3}x \\ 0 \end{pmatrix}, * \phi_2 = \begin{pmatrix} -\sqrt{3}+2\sqrt{3}y \\ 0 \end{pmatrix}, diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 293b8f3751..d6217ca021 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -457,16 +457,26 @@ public: * International Journal for Numerical Methods in Engineering, vol 24, * pages 959–973. year 1987, the author applies the transformation on the * reference cell $[-1, 1]$ getting - @f{align*} n(1) &= 1, \\ n(-1) &= -1, \\ \frac{dn}{dx} &= 0 \text{ at } - x = x_0, \\ \frac{d^2n}{dx^2} &= 0 \text{ at } x = x_0 @f} + * @f{align*}{ + * n(1) &= 1, \\ n(-1) &= -1, \\ \frac{dn}{dx} &= 0 \text{ at } + * x = x_0, \\ \frac{d^2n}{dx^2} &= 0 \text{ at } x = x_0 + * @f} * We get - @f{align*} a &= \frac{1}{q}, \\ b &= -3 \frac{\bar{\Gamma}}{q}, \\ - c &= 3\frac{\bar{\Gamma}}{q}, \\ d &= -b, @f} + * @f{align*}{ + * a &= \frac{1}{q}, \\ + * b &= -3 \frac{\bar{\Gamma}}{q}, \\ + * c &= 3 \frac{\bar{\Gamma}}{q}, \\ + * d &= -b, + * @f} * with - @f{align*} \eta^{*} &= \bar{\eta}^2 - 1, \\ \bar{\Gamma} &= \sqrt[3]{ - \bar{\eta} \eta^{*} + |\eta^{*} | } + \sqrt[3]{ \bar{\eta} \eta^{*} - - |\eta^{*} | } + \bar{\eta}, \\ q &= (\Gamma-\bar{\Gamma})^3 + \bar{\Gamma} - \frac{\bar{\Gamma}^2+3}{1+3\bar{\Gamma}^2} @f} + * @f{align*}{ + * \eta^{*} &= \bar{\eta}^2 - 1, \\ + * \bar{\Gamma} &= \sqrt[3]{\bar{\eta} \eta^{*} + |\eta^{*} | } + * + \sqrt[3]{ \bar{\eta} \eta^{*} - |\eta^{*} | } + * + \bar{\eta}, \\ + * q &= (\Gamma-\bar{\Gamma})^3 + \bar{\Gamma} + * \frac{\bar{\Gamma}^2+3}{1+3\bar{\Gamma}^2} + * @f} * Since the library assumes $[0,1]$ as reference interval, we will map * these values on the proper reference interval in the implementation. * diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index 7ef51ae428..b0b41524e3 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -490,7 +490,9 @@ namespace DoFTools * * For example, if you wanted to solve the Stokes equations, * - * @f{align*} -\Delta \mathbf u + \nabla p &= 0,\\ \text{div}\ u &= 0 @f} + * @f{align*}{ + * -\Delta \mathbf u + \nabla p &= 0,\\ \text{div}\ u &= 0 + * @f} * * in two space dimensions, using stable Q2/Q1 mixed elements (using * the FESystem class), then you don't want all degrees of freedom diff --git a/include/deal.II/dofs/function_map.h b/include/deal.II/dofs/function_map.h index 94dc275396..9d54d96421 100644 --- a/include/deal.II/dofs/function_map.h +++ b/include/deal.II/dofs/function_map.h @@ -34,15 +34,22 @@ template class Function; * function that provides the boundary values on this part of the boundary. * This type is required in many functions in the library where, for example, * we need to know about the functions $h_i(\mathbf x)$ used in boundary - * conditions @f{align*} \mathbf n \cdot \nabla u = h_i \qquad \qquad - * \text{on}\ \Gamma_i\subset\partial\Omega. @f} An example is the function + * conditions + * @f{align*}{ + * \mathbf n \cdot \nabla u = h_i \qquad \qquad + * \text{on}\ \Gamma_i\subset\partial\Omega. + * @f} + * An example is the function * KellyErrorEstimator::estimate() that allows us to provide a set of * functions $h_i$ for all those boundary indicators $i$ for which the * boundary condition is supposed to be of Neumann type. Of course, the same * kind of principle can be applied to cases where we care about Dirichlet * values, where one needs to provide a map from boundary indicator $i$ to - * Dirichlet function $h_i$ if the boundary conditions are given as @f{align*} - * u = h_i \qquad \qquad \text{on}\ \Gamma_i\subset\partial\Omega. @f} This + * Dirichlet function $h_i$ if the boundary conditions are given as + * @f{align*}{ + * u = h_i \qquad \qquad \text{on}\ \Gamma_i\subset\partial\Omega. + * @f} + * This * is, for example, the case for the VectorTools::interpolate() functions. * * Tutorial programs step-6, step-7 and step-8 show examples of how to use diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 728e0390fc..b88ae6e6eb 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -630,11 +630,15 @@ namespace FEValuesViews * for the shape function and quadrature point selected by the arguments. * For 1d this function does not make any sense. Thus it is not * implemented for spacedim=1. In 2d the curl is defined as - * @f{equation*} \operatorname{curl}(u):=\frac{du_2}{dx} -\frac{du_1}{dy}, - * @f} whereas in 3d it is given by @f{equation*} + * @f{equation*}{ + * \operatorname{curl}(u):=\frac{du_2}{dx} -\frac{du_1}{dy}, + * @f} + * whereas in 3d it is given by + * @f{equation*}{ * \operatorname{curl}(u):=\left( \begin{array}{c} * \frac{du_3}{dy}-\frac{du_2}{dz}\\ \frac{du_1}{dz}-\frac{du_3}{dx}\\ - * \frac{du_2}{dx}-\frac{du_1}{dy} \end{array} \right). @f} + * \frac{du_2}{dx}-\frac{du_1}{dy} \end{array} \right). + * @f} * * @note The meaning of the arguments is as documented for the value() * function. diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h index 9d8a1ac7b4..69579ee458 100644 --- a/include/deal.II/fe/mapping_cartesian.h +++ b/include/deal.II/fe/mapping_cartesian.h @@ -36,11 +36,17 @@ DEAL_II_NAMESPACE_OPEN * In other words, the mapping is meant for cells for which the mapping from * the reference to the real cell is a scaling along the coordinate * directions: The transformation from reference coordinates $\hat {\mathbf - * x}$ to real coordinates $\mathbf x$ on each cell is of the form @f{align*} + * x}$ to real coordinates $\mathbf x$ on each cell is of the form + * @f{align*}{ * {\mathbf x}(\hat {\mathbf x}) = \begin{pmatrix} h_x & 0 \\ 0 & h_y - * \end{pmatrix} \hat{\mathbf x} + {\mathbf v}_0 @f} in 2d, and @f{align*} + * \end{pmatrix} \hat{\mathbf x} + {\mathbf v}_0 + * @f} + * in 2d, and + * @f{align*}{ * {\mathbf x}(\hat {\mathbf x}) = \begin{pmatrix} h_x & 0 & 0 \\ 0 & h_y & 0 - * \\ 0 & 0 & h_z \end{pmatrix} \hat{\mathbf x} + {\mathbf v}_0 @f} in 3d, + * \\ 0 & 0 & h_z \end{pmatrix} \hat{\mathbf x} + {\mathbf v}_0 + * @f} + * in 3d, * where ${\mathbf v}_0$ is the bottom left vertex and $h_x,h_y,h_z$ are the * extents of the cell along the axes. * diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index 57d07443fc..3324fa1242 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -636,16 +636,21 @@ block_diagonal_operator(const LinearOperatorx0 = A00^-1 y0. Then, we can use x0 to recover x1: + * @code * x1 = A11^-1 ( y1 - A01 x0 ) + * @endcode * and therefore: + * @code * xn = Ann^-1 ( yn - A0n x0 - ... - A(n-1)n x(n-1) ) + * @endcode * * @note We are not using all blocks of the BlockLinearOperator arguments: * Just the lower triangular block matrix of @p block_operator is used as @@ -742,16 +747,21 @@ block_forward_substitution(const BlockLinearOperator &block_opera * * Let us assume we have a linear system with the following block structure: * + * @code * A00 x0 + A01 x1 + ... + A0n xn = yn * A11 x1 + ... = y1 * ... .. * Ann xn = yn + * @endcode * - * First of all, xn = Ann^-1 yn - * Then, we can use xn to recover x(n-1): + * First of all, xn = Ann^-1 yn. Then, we can use xn to recover x(n-1): + * @code * x(n-1) = A(n-1)(n-1)^-1 ( y(n-1) - A(n-1)n x(n-1) ) + * @endcode * and therefore: + * @code * x0 = A00^-1 ( y0 - A0n xn - ... - A01 x1 ) + * @endcode * * @note We are not using all blocks of the BlockLinearOperator arguments: * Just the upper triangular block matrix of @p block_operator is used as diff --git a/include/deal.II/lac/schur_complement.h b/include/deal.II/lac/schur_complement.h index 156501fb3c..9ee77775bb 100644 --- a/include/deal.II/lac/schur_complement.h +++ b/include/deal.II/lac/schur_complement.h @@ -44,50 +44,50 @@ DEAL_II_NAMESPACE_OPEN * * Consider a general system of linear equations that can be * decomposed into two major sets of equations: - * @f{eqnarray*} - \mathbf{K}\mathbf{d} = \mathbf{f} - \quad \Rightarrow\quad - \left(\begin{array}{cc} - A & B \\ C & D - \end{array}\right) - \left(\begin{array}{cc} - x \\ y - \end{array}\right) - = - \left(\begin{array}{cc} - f \\ g - \end{array}\right), + * @f{eqnarray*}{ + * \mathbf{K}\mathbf{d} = \mathbf{f} + * \quad \Rightarrow\quad + * \left(\begin{array}{cc} + * A & B \\ C & D + * \end{array}\right) + * \left(\begin{array}{cc} + * x \\ y + * \end{array}\right) + * = + * \left(\begin{array}{cc} + * f \\ g + * \end{array}\right), * @f} * where $ A,B,C,D $ represent general subblocks of the matrix * $ \mathbf{K} $ and, similarly, general subvectors of * $ \mathbf{d},\mathbf{f} $ are given by $ x,y,f,g $ . * * This is equivalent to the following two statements: - * @f{eqnarray*} - (1) \quad Ax + By &=& f \\ - (2) \quad Cx + Dy &=& g \quad . + * @f{eqnarray*}{ + * (1) \quad Ax + By &=& f \\ + * (2) \quad Cx + Dy &=& g \quad . * @f} * * Assuming that $ A,D $ are both square and invertible, we could * then perform one of two possible substitutions, - * @f{eqnarray*} - (3) \quad x &=& A^{-1}(f - By) \quad \text{from} \quad (1) \\ - (4) \quad y &=& D^{-1}(g - Cx) \quad \text{from} \quad (2) , + * @f{eqnarray*}{ + * (3) \quad x &=& A^{-1}(f - By) \quad \text{from} \quad (1) \\ + * (4) \quad y &=& D^{-1}(g - Cx) \quad \text{from} \quad (2) , * @f} * which amount to performing block Gaussian elimination on * this system of equations. * * For the purpose of the current implementation, we choose to * substitute (3) into (2) - * @f{eqnarray*} - C \: A^{-1}(f - By) + Dy &=& g \\ - -C \: A^{-1} \: By + Dy &=& g - C \: A^{-1} \: f \quad . - @f} + * @f{eqnarray*}{ + * C \: A^{-1}(f - By) + Dy &=& g \\ + * -C \: A^{-1} \: By + Dy &=& g - C \: A^{-1} \: f \quad . + * @f} * This leads to the result * @f[ - (5) \quad (D - C\: A^{-1} \:B)y = g - C \: A^{-1} f - \quad \Rightarrow \quad Sy = g' - @f] + * (5) \quad (D - C\: A^{-1} \:B)y = g - C \: A^{-1} f + * \quad \Rightarrow \quad Sy = g' + * @f] * with $ S = (D - C\: A^{-1} \:B) $ being the Schur complement * and the modified right-hand side vector $ g' = g - C \: A^{-1} f $ arising from * the condensation step. @@ -98,8 +98,8 @@ DEAL_II_NAMESPACE_OPEN * So for any arbitrary vector $ a $, the Schur complement * performs the following operation: * @f[ - (6) \quad Sa = (D - C \: A^{-1} \: B)a - @f] + * (6) \quad Sa = (D - C \: A^{-1} \: B)a + * @f] * * A typical set of steps needed the solve a linear system (1),(2) * would be: @@ -119,51 +119,51 @@ DEAL_II_NAMESPACE_OPEN * 4. Perform pre-processing step on the RHS of (5) using * condense_schur_rhs(): * @f[ - g' = g - C \: A^{-1} \: f - @f] + * g' = g - C \: A^{-1} \: f + * @f] * 5. Solve for $ y $ in (5): * @f[ - y = S^{-1} g' - @f] + * y = S^{-1} g' + * @f] * 6. Perform the post-processing step from (3) using * postprocess_schur_solution(): * @f[ - x = A^{-1} (f - By) - @f] + * x = A^{-1} (f - By) + * @f] * * An illustration of typical usage of this operator for a fully coupled * system is given below. * @code - #include - - // Given BlockMatrix K and BlockVectors d,F - - // Decomposition of tangent matrix - const auto A = linear_operator(K.block(0,0)); - const auto B = linear_operator(K.block(0,1)); - const auto C = linear_operator(K.block(1,0)); - const auto D = linear_operator(K.block(1,1)); - - // Decomposition of solution vector - auto x = d.block(0); - auto y = d.block(1); - - // Decomposition of RHS vector - auto f = F.block(0); - auto g = F.block(1); - - // Construction of inverse of Schur complement - const auto prec_A = PreconditionSelector<...>(A); - const auto A_inv = inverse_operator<...>(A,prec_A); - const auto S = schur_complement(A_inv,B,C,D); - const auto S_prec = PreconditionSelector<...>(D); // D and S operate on same space - const auto S_inv = inverse_operator<...>(S,...,prec_S); - - // Solve reduced block system - auto rhs = condense_schur_rhs (A_inv,C,f,g); // PackagedOperation that represents the condensed form of g - y = S_inv * rhs; // Solve for y - x = postprocess_schur_solution (A_inv,B,y,f); // Compute x using resolved solution y - @endcode + * #include + * + * // Given BlockMatrix K and BlockVectors d,F + * + * // Decomposition of tangent matrix + * const auto A = linear_operator(K.block(0,0)); + * const auto B = linear_operator(K.block(0,1)); + * const auto C = linear_operator(K.block(1,0)); + * const auto D = linear_operator(K.block(1,1)); + * + * // Decomposition of solution vector + * auto x = d.block(0); + * auto y = d.block(1); + * + * // Decomposition of RHS vector + * auto f = F.block(0); + * auto g = F.block(1); + * + * // Construction of inverse of Schur complement + * const auto prec_A = PreconditionSelector<...>(A); + * const auto A_inv = inverse_operator<...>(A,prec_A); + * const auto S = schur_complement(A_inv,B,C,D); + * const auto S_prec = PreconditionSelector<...>(D); // D and S operate on same space + * const auto S_inv = inverse_operator<...>(S,...,prec_S); + * + * // Solve reduced block system + * auto rhs = condense_schur_rhs (A_inv,C,f,g); // PackagedOperation that represents the condensed form of g + * y = S_inv * rhs; // Solve for y + * x = postprocess_schur_solution (A_inv,B,y,f); // Compute x using resolved solution y + * @endcode * * In the above example, the preconditioner for $ S $ was defined as the * preconditioner for $ D $, which is valid since they operate on the same @@ -184,21 +184,21 @@ DEAL_II_NAMESPACE_OPEN * Thereafter we construct the approximate inverse operator $ \tilde{S}^{-1} $ * which is then used as the preconditioner for computing $ S^{-1} $. * @code - // Construction of approximate inverse of Schur complement - const auto A_inv_approx = linear_operator(preconditioner_A); - const auto S_approx = schur_complement(A_inv_approx,B,C,D); - const auto S_approx_prec = PreconditionSelector<...>(D); // D and S_approx operate on same space - const auto S_inv_approx = inverse_operator(S_approx,...,S_approx_prec); // Inner solver: Typically limited to few iterations using IterationNumberControl - - // Construction of exact inverse of Schur complement - const auto S = schur_complement(A_inv,B,C,D); - const auto S_inv = inverse_operator(S,...,S_inv_approx); // Outer solver - - // Solve reduced block system - auto rhs = condense_schur_rhs (A_inv,C,f,g); - y = S_inv * rhs; // Solve for y - x = postprocess_schur_solution (A_inv,B,y,f); - @endcode + * // Construction of approximate inverse of Schur complement + * const auto A_inv_approx = linear_operator(preconditioner_A); + * const auto S_approx = schur_complement(A_inv_approx,B,C,D); + * const auto S_approx_prec = PreconditionSelector<...>(D); // D and S_approx operate on same space + * const auto S_inv_approx = inverse_operator(S_approx,...,S_approx_prec); // Inner solver: Typically limited to few iterations using IterationNumberControl + * + * // Construction of exact inverse of Schur complement + * const auto S = schur_complement(A_inv,B,C,D); + * const auto S_inv = inverse_operator(S,...,S_inv_approx); // Outer solver + * + * // Solve reduced block system + * auto rhs = condense_schur_rhs (A_inv,C,f,g); + * y = S_inv * rhs; // Solve for y + * x = postprocess_schur_solution (A_inv,B,y,f); + * @endcode * Note that due to the construction of @c S_inv_approx and subsequently * @c S_inv, there are a pair of nested iterative solvers which could * collectively consume a lot of resources. @@ -353,7 +353,7 @@ schur_complement(const LinearOperator &A_inv, * @relates PackagedOperation * * For the system of equations - * @f{eqnarray*} + * @f{eqnarray*}{ Ax + By &=& f \\ Cx + Dy &=& g \quad , * @f} @@ -436,7 +436,7 @@ condense_schur_rhs (const LinearOperator &A_inv, * @relates PackagedOperation * * For the system of equations - * @f{eqnarray*} + * @f{eqnarray*}{ Ax + By &=& f \\ Cx + Dy &=& g \quad , * @f} diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 5a14d4f555..ef07f6f6c6 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -906,9 +906,13 @@ namespace VectorTools * Project a function or a set of functions to the boundary of the domain. * In other words, compute the solution of the following problem: Find $u_h * \in V_h$ (where $V_h$ is the finite element space represented by the - * DoFHandler argument of this function) so that @f{align*} \int_{\Gamma} - * \varphi_i u_h = \sum_{k \in {\cal K}} \int_{\Gamma_k} \varphi_i f_k, - * \qquad \forall \varphi_i \in V_h @f} where $\Gamma = \bigcup_{k \in {\cal + * DoFHandler argument of this function) so that + * @f{align*}{ + * \int_{\Gamma} \varphi_i u_h + * = \sum_{k \in {\cal K}} \int_{\Gamma_k} \varphi_i f_k, + * \qquad \forall \varphi_i \in V_h + * @f} + * where $\Gamma = \bigcup_{k \in {\cal * K}} \Gamma_k$, $\Gamma_k \subset \partial\Omega$, $\cal K$ is the set of * indices and $f_k$ the corresponding boundary functions represented in the * function map argument @p boundary_values to this function, and the @@ -1812,7 +1816,11 @@ namespace VectorTools * function object, and a finite element function. The result of this * function is the vector @p difference that contains one value per active * cell $K$ of the triangulation. Each of the values of this vector $d$ - * equals @f{align*} d_K = \| u-u_h \|_X @f} where $X$ denotes the norm + * equals + * @f{align*}{ + * d_K = \| u-u_h \|_X + * @f} + * where $X$ denotes the norm * chosen and $u$ represents the exact solution. * * It is assumed that the number of components of the function @p @@ -1896,8 +1904,11 @@ namespace VectorTools * @endcode * The squaring and taking the square root is necessary in order to compute * the sum of squares of norms over all all cells in the definition of the - * $L_2$ norm: @f{align*} \textrm{error} = \sqrt{\sum_K - * \|u-u_h\|_{L_2(K)}^2} @f} Obviously, if you are interested in computing + * $L_2$ norm: + * @f{align*}{ + * \textrm{error} = \sqrt{\sum_K \|u-u_h\|_{L_2(K)}^2} + * @f} + * Obviously, if you are interested in computing * the $L_1$ norm of the error, the correct form of the last two lines would * have been * @code