From b498b7cbe7ff3a00148f5431461bda50fcc75953 Mon Sep 17 00:00:00 2001 From: maier Date: Mon, 8 Apr 2013 18:27:48 +0000 Subject: [PATCH] Improve VectorTools::subtract_mean_value git-svn-id: https://svn.dealii.org/trunk@29220 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 +- .../include/deal.II/numerics/vector_tools.h | 182 +++++++----------- .../deal.II/numerics/vector_tools.templates.h | 48 ++++- deal.II/source/numerics/vector_tools.cc | 36 ---- deal.II/source/numerics/vector_tools.inst.in | 9 + 5 files changed, 134 insertions(+), 150 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0d7ec5d70a..b89c937287 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -87,13 +87,20 @@ this function.

Specific improvements

+
    +
  1. Fixed: VectorTools::subtract_mean_value can now be called without a +boolean mask. The vector type is templatified and instantiated for all +non distributed vectors. +
    +(Matthias Maier, 2013/04/08) +
    2. +
  2. Fixed: It is now possible to call ConvergenceTable::evaluate_convergence_rates multiple times.
    (Matthias Maier, 2013/04/08)
    3. -
    1. Fixed: GridTools::distort_random (previously called Triangulation::distort_random) had a bug where points were only ever moved in positive coordinate directions rather than with uniform probability in either direction. The 1d diff --git a/deal.II/include/deal.II/numerics/vector_tools.h b/deal.II/include/deal.II/numerics/vector_tools.h index a61a1264cb..4b72071289 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@ -2173,117 +2173,80 @@ namespace VectorTools //@{ /** - * Subtract the (algebraic) mean value - * from a vector. This function is most - * frequently used as a mean-value filter - * for Stokes: The pressure in Stokes' - * equations with only Dirichlet - * boundaries for the velocities is only - * determined up to a constant. This - * function allows to subtract the mean - * value of the pressure. It is usually - * called in a preconditioner and - * generates updates with mean value - * zero. The mean value is computed as - * the mean value of the degrees of - * freedom values as given by the input - * vector; they are not weighted by the - * area of cells, i.e. the mean is - * computed as $\sum_i v_i$, rather than - * as $\int_\Omega v(x) = \int_\Omega - * \sum_i v_i \phi_i(x)$. The latter can - * be obtained from the - * VectorTools::compute_mean_function, - * however. - * - * Apart from the vector @p v to operate - * on, this function takes a boolean mask - * that has a true entry for every element - * of the vector for which the mean value - * shall be computed and later - * subtracted. The argument is used to - * denote which components of the solution - * vector correspond to the pressure, and - * avoid touching all other components of - * the vector, such as the velocity - * components. (Note, however, that the - * mask is not a @ref GlossComponentMask - * operating on the vector components of - * the finite element the solution vector - * @p v may be associated with; rather, it - * is a mask on the entire vector, without - * reference to what the vector elements - * mean.) - * - * @note In the context of using this - * function to filter out the kernel of - * an operator (such as the null space of - * the Stokes operator that consists of - * the constant pressures), this function - * only makes sense for finite elements - * for which the null space indeed - * consists of the vector - * $(1,1,\ldots,1)^T$. This is the case - * for example for the usual Lagrange - * elements where the sum of all shape - * functions equals the function that is - * constant one. However, it is not true - * for some other functions: for example, - * for the FE_DGP element (another valid - * choice for the pressure in Stokes - * discretizations), the first shape - * function on each cell is constant - * while further elements are $L_2$ - * orthogonal to it (on the reference - * cell); consequently, the sum of all - * shape functions is not equal to one, - * and the vector that is associated with - * the constant mode is not equal to - * $(1,1,\ldots,1)^T$. For such elements, - * a different procedure has to be used - * when subtracting the mean value. + * Subtract the (algebraic) mean value from a vector. + * + * This function is most frequently used as a mean-value filter for + * Stokes: The pressure in Stokes' equations with only Dirichlet + * boundaries for the velocities is only determined up to a constant. + * This function allows to subtract the mean value of the pressure. It is + * usually called in a preconditioner and generates updates with mean + * value zero. The mean value is computed as the mean value of the + * degrees of freedom values as given by the input vector; they are not + * weighted by the area of cells, i.e. the mean is computed as $\sum_i + * v_i$, rather than as $\int_\Omega v(x) = \int_\Omega \sum_i v_i + * \phi_i(x)$. The latter can be obtained from the + * VectorTools::compute_mean_function, however. + * + * Apart from the vector @p v to operate on, this function takes a + * boolean mask @p p_select that has a true entry for every element of + * the vector for which the mean value shall be computed and later + * subtracted. The argument is used to denote which components of the + * solution vector correspond to the pressure, and avoid touching all + * other components of the vector, such as the velocity components. + * (Note, however, that the mask is not a @ref GlossComponentMask + * operating on the vector components of the finite element the solution + * vector @p v may be associated with; rather, it is a mask on the entire + * vector, without reference to what the vector elements mean.) + * + * The boolean mask @p p_select has an empty vector as default value, + * which will be interpreted as selecting all vector elements, hence, + * subtracting the algebraic mean value on the whole vector. This allows + * to call this function without a boolean mask if the whole vector + * should be processed. + * + * @note In the context of using this function to filter out the kernel + * of an operator (such as the null space of the Stokes operator that + * consists of the constant pressures), this function only makes sense + * for finite elements for which the null space indeed consists of the + * vector $(1,1,\ldots,1)^T$. This is the case for example for the usual + * Lagrange elements where the sum of all shape functions equals the + * function that is constant one. However, it is not true for some other + * functions: for example, for the FE_DGP element (another valid choice + * for the pressure in Stokes discretizations), the first shape function + * on each cell is constant while further elements are $L_2$ orthogonal + * to it (on the reference cell); consequently, the sum of all shape + * functions is not equal to one, and the vector that is associated with + * the constant mode is not equal to $(1,1,\ldots,1)^T$. For such + * elements, a different procedure has to be used when subtracting the + * mean value. */ - void subtract_mean_value(Vector &v, - const std::vector &p_select); - - /** - * Compute the mean value of one - * component of the solution. - * - * This function integrates the - * chosen component over the - * whole domain and returns the - * result, i.e. it computes - * $\int_\Omega [u_h(x)]_c \; dx$ - * where $c$ is the vector component - * and $u_h$ is the function - * representation of the nodal - * vector given as fourth - * argument. The integral is evaluated - * numerically using the quadrature - * formula given as third argument. + template + void subtract_mean_value(VECTOR &v, + const std::vector &p_select = std::vector()); + + + /** + * Compute the mean value of one component of the solution. * - * This function is used in the - * "Possibilities for extensions" part of - * the results section of @ref step_3 - * "step-3". + * This function integrates the chosen component over the whole domain + * and returns the result, i.e. it computes $\int_\Omega [u_h(x)]_c \; + * dx$ where $c$ is the vector component and $u_h$ is the function + * representation of the nodal vector given as fourth argument. The + * integral is evaluated numerically using the quadrature formula given + * as third argument. * - * @note The function is most often used - * when solving a problem whose solution - * is only defined up to a constant, for - * example a pure Neumann problem or the - * pressure in a Stokes or Navier-Stokes - * problem. In both cases, subtracting - * the mean value as computed by the - * current function, from the nodal - * vector does not generally yield the - * desired result of a finite element - * function with mean value zero. In - * fact, it only works for Lagrangian - * elements. For all other elements, you - * will need to compute the mean value - * and subtract it right inside the - * evaluation routine. + * This function is used in the "Possibilities for extensions" part of + * the results section of @ref step_3 "step-3". + * + * @note The function is most often used when solving a problem whose + * solution is only defined up to a constant, for example a pure Neumann + * problem or the pressure in a Stokes or Navier-Stokes problem. In both + * cases, subtracting the mean value as computed by the current function, + * from the nodal vector does not generally yield the desired result of a + * finite element function with mean value zero. In fact, it only works + * for Lagrangian elements. For all other elements, you will need to + * compute the mean value and subtract it right inside the evaluation + * routine. */ template double compute_mean_value (const Mapping &mapping, @@ -2293,8 +2256,7 @@ namespace VectorTools const unsigned int component); /** - * Calls the other compute_mean_value() - * function, see above, with + * Calls the other compute_mean_value() function, see above, with * mapping=MappingQ1@(). */ template diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index b62546dbe2..33bde2ad2f 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -5635,6 +5635,49 @@ namespace VectorTools } + + template + void + subtract_mean_value(VECTOR &v, + const std::vector &p_select) + { + const unsigned int n = v.size(); + Assert(p_select.size() == 0 || p_select.size() == n, + ExcDimensionMismatch(p_select.size(), n)); + + typename VECTOR::value_type s = 0.; + + if(p_select.size() == 0) + { + // In case of an empty boolean mask operate on the whole vector: + for (unsigned int i=0; i 0, ExcNoComponentSelected()); + + s /= counter; + + for (unsigned int i=0; i double compute_mean_value (const Mapping &mapping, @@ -5672,9 +5715,8 @@ namespace VectorTools } #ifdef DEAL_II_WITH_P4EST - // if this was a distributed - // DoFHandler, we need to do the - // reduction over the entire domain + // if this was a distributed DoFHandler, we need to do the reduction + // over the entire domain if (const parallel::distributed::Triangulation * p_d_triangulation = dynamic_cast *>(&dof.get_tria())) diff --git a/deal.II/source/numerics/vector_tools.cc b/deal.II/source/numerics/vector_tools.cc index c4a0be5f6f..af44a3a505 100644 --- a/deal.II/source/numerics/vector_tools.cc +++ b/deal.II/source/numerics/vector_tools.cc @@ -11,47 +11,11 @@ // //--------------------------------------------------------------------------- -#include -#include - #include DEAL_II_NAMESPACE_OPEN - -namespace VectorTools -{ - - void - subtract_mean_value(Vector &v, - const std::vector &p_select) - { - const unsigned int n = v.size(); - Assert(n == p_select.size(), - ExcDimensionMismatch(n, p_select.size())); - - double s = 0; - unsigned int counter = 0; - - for (unsigned int i=0; i 0, ExcNoComponentSelected()); - - s /= counter; - - for (unsigned int i=0; i &, const std::vector &); + \} +} + + for (deal_II_dimension : DIMENSIONS) { namespace VectorTools \{ -- 2.39.5