From: Timo Heister Date: Sun, 15 May 2016 15:16:18 +0000 (+0100) Subject: implement VectorTools::compute_global_error, document norms X-Git-Tag: v8.5.0-rc1~1022^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9dbe5f5080d55b0b7b0ead7ef48cc3f84c8d76f8;p=dealii.git implement VectorTools::compute_global_error, document norms Add a function that computes global errors from cellwise errors obtained by VectorTools::integrate_difference(). Do MPI collective if necessary. Also add a large amount of documentation to NormType. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 88ce9b37a2..08159b1464 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -120,6 +120,13 @@ inconvenience this causes.

General

    +
  1. New: Add VectorTools::compute_global_error that computes global + errors from cellwise errors obtained by VectorTools::integrate_difference() + and do MPI collectives if necessary. +
    + (Timo Heister, 2016/05/15) +
    2. +
  2. New: Add functions to transform Cartesian coordinates to spherical and back: GeometricUtilities::Coordinates::to_spherical and GeometricUtilities::Coordinates::from_spherical. diff --git a/examples/step-11/step-11.cc b/examples/step-11/step-11.cc index 8b39241786..ea9b17a48f 100644 --- a/examples/step-11/step-11.cc +++ b/examples/step-11/step-11.cc @@ -367,9 +367,10 @@ namespace Step11 VectorTools::H1_seminorm); // Then, the function just called returns its results as a vector of // values each of which denotes the norm on one cell. To get the global - // norm, a simple computation shows that we have to take the l2 norm of - // the vector: - const double norm = norm_per_cell.l2_norm(); + // norm, we do the following: + const double norm = VectorTools::compute_global_error(triangulation, + norm_per_cell, + VectorTools::H1_seminorm); // Last task -- generate output: output_table.add_value ("cells", triangulation.n_active_cells()); diff --git a/examples/step-20/step-20.cc b/examples/step-20/step-20.cc index fb1be4ff95..72bb761cb0 100644 --- a/examples/step-20/step-20.cc +++ b/examples/step-20/step-20.cc @@ -777,9 +777,9 @@ namespace Step20 // frequently the case in mixed finite element applications). What we // therefore have to do is to `mask' the components that we are interested // in. This is easily done: the - // VectorTools::integrate_difference function takes as its last - // argument a pointer to a weight function (the parameter defaults to the - // null pointer, meaning unit weights). What we simply have to do is to pass + // VectorTools::integrate_difference function takes as one of its + // arguments a pointer to a weight function (the parameter defaults to the + // null pointer, meaning unit weights). What we have to do is to pass // a function object that equals one in the components we are interested in, // and zero in the other ones. For example, to compute the pressure error, // we should pass a function that represents the constant vector with a unit @@ -829,13 +829,17 @@ namespace Step20 cellwise_errors, quadrature, VectorTools::L2_norm, &pressure_mask); - const double p_l2_error = cellwise_errors.l2_norm(); + const double p_l2_error = VectorTools::compute_global_error(triangulation, + cellwise_errors, + VectorTools::L2_norm); VectorTools::integrate_difference (dof_handler, solution, exact_solution, cellwise_errors, quadrature, VectorTools::L2_norm, &velocity_mask); - const double u_l2_error = cellwise_errors.l2_norm(); + const double u_l2_error = VectorTools::compute_global_error(triangulation, + cellwise_errors, + VectorTools::L2_norm); std::cout << "Errors: ||e_p||_L2 = " << p_l2_error << ", ||e_u||_L2 = " << u_l2_error diff --git a/examples/step-34/step-34.cc b/examples/step-34/step-34.cc index dd54cbb4d7..87e4caa763 100644 --- a/examples/step-34/step-34.cc +++ b/examples/step-34/step-34.cc @@ -768,8 +768,9 @@ namespace Step34 difference_per_cell, QGauss<(dim-1)>(2*fe.degree+1), VectorTools::L2_norm); - const double L2_error = difference_per_cell.l2_norm(); - + const double L2_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); // The error in the alpha vector can be computed directly using the // Vector::linfty_norm() function, since on each node, the value should be diff --git a/examples/step-38/step-38.cc b/examples/step-38/step-38.cc index 70a1958ee2..9035706e21 100644 --- a/examples/step-38/step-38.cc +++ b/examples/step-38/step-38.cc @@ -525,8 +525,11 @@ namespace Step38 QGauss(2*fe.degree+1), VectorTools::H1_norm); + double h1_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_norm); std::cout << "H1 error = " - << difference_per_cell.l2_norm() + << h1_error << std::endl; } diff --git a/examples/step-48/step-48.cc b/examples/step-48/step-48.cc index 24e7c5694b..df18653e27 100644 --- a/examples/step-48/step-48.cc +++ b/examples/step-48/step-48.cc @@ -467,7 +467,9 @@ namespace Step48 QGauss(fe_degree+1), VectorTools::L2_norm); const double solution_norm = - std::sqrt(Utilities::MPI::sum (norm_per_cell.norm_sqr(), MPI_COMM_WORLD)); + VectorTools::compute_global_error(triangulation, + norm_per_cell, + VectorTools::L2_norm); pcout << " Time:" << std::setw(8) << std::setprecision(3) << time diff --git a/examples/step-51/step-51.cc b/examples/step-51/step-51.cc index 80434c8917..aa2d32a578 100644 --- a/examples/step-51/step-51.cc +++ b/examples/step-51/step-51.cc @@ -1052,7 +1052,9 @@ namespace Step51 QGauss(fe.degree+2), VectorTools::L2_norm, &value_select); - const double L2_error = difference_per_cell.l2_norm(); + const double L2_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); ComponentSelectFunction gradient_select (std::pair(0, dim), dim+1); @@ -1063,7 +1065,9 @@ namespace Step51 QGauss(fe.degree+2), VectorTools::L2_norm, &gradient_select); - const double grad_error = difference_per_cell.l2_norm(); + const double grad_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); VectorTools::integrate_difference (dof_handler_u_post, solution_u_post, @@ -1071,7 +1075,9 @@ namespace Step51 difference_per_cell, QGauss(fe.degree+3), VectorTools::L2_norm); - const double post_error = difference_per_cell.l2_norm(); + const double post_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); convergence_table.add_value("cells", triangulation.n_active_cells()); convergence_table.add_value("dofs", dof_handler.n_dofs()); diff --git a/examples/step-7/step-7.cc b/examples/step-7/step-7.cc index 69a19b9a97..f1c4a197ac 100644 --- a/examples/step-7/step-7.cc +++ b/examples/step-7/step-7.cc @@ -867,22 +867,26 @@ namespace Step7 difference_per_cell, QGauss(3), VectorTools::L2_norm); - const double L2_error = difference_per_cell.l2_norm(); + const double L2_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); // By same procedure we get the H1 semi-norm. We re-use the // difference_per_cell vector since it is no longer used // after computing the L2_error variable above. The global // $H^1$ semi-norm error is then computed by taking the sum of squares // of the errors on each individual cell, and then the square root of - // it -- an operation that conveniently again coincides with taking - // the $l_2$ norm of the vector of error indicators. + // it -- an operation that is conveniently performed by + // VectorTools::compute_global_error. VectorTools::integrate_difference (dof_handler, solution, Solution(), difference_per_cell, QGauss(3), VectorTools::H1_seminorm); - const double H1_error = difference_per_cell.l2_norm(); + const double H1_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_seminorm); // Finally, we compute the maximum norm. Of course, we can't actually // compute the true maximum, but only the maximum at the quadrature @@ -896,10 +900,8 @@ namespace Step7 // // Using this special quadrature rule, we can then try to find the maximal // error on each cell. Finally, we compute the global L infinity error - // from the L infinite errors on each cell. Instead of summing squares, we - // now have to take the maximum value over all cell-wise entries, an - // operation that is conveniently done using the Vector::linfty() - // function: + // from the L infinity errors on each cell with a call to + // VectorTools::compute_global_error. const QTrapez<1> q_trapez; const QIterated q_iterated (q_trapez, 5); VectorTools::integrate_difference (dof_handler, @@ -908,7 +910,9 @@ namespace Step7 difference_per_cell, q_iterated, VectorTools::Linfty_norm); - const double Linfty_error = difference_per_cell.linfty_norm(); + const double Linfty_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::Linfty_norm); // After all these errors have been computed, we finally write some // output. In addition, we add the important data to the TableHandler by diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index c617a7ad7b..a6e8638a25 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -238,8 +238,10 @@ class ConstraintMatrix; * * This data, one number per active cell, can be used to generate graphical * output by directly passing it to the DataOut class through the - * DataOut::add_data_vector function. Alternatively, it can be interpolated to - * the nodal points of a finite element field using the + * DataOut::add_data_vector function. Alternatively, the global error can be + * computed using VectorTools::compute_global_error(). Finally, the output per + * cell from VectorTools::integrate_difference() can be interpolated to the + * nodal points of a finite element field using the * DoFTools::distribute_cell_to_dof_vector function. * * Presently, there is the possibility to compute the following values from @@ -317,61 +319,243 @@ namespace VectorTools { /** * Denote which norm/integral is to be computed by the - * integrate_difference() function of this namespace. The following - * possibilities are implemented: + * integrate_difference() function on each cell and compute_global_error() + * for the whole domain. + * Let $f:\Omega \rightarrow \mathbb{R}^c$ be a finite element function + * with $c$ components where component $c$ is denoted by $f_c$ and $\hat{f}$ + * be the reference function (the @p fe_function and @p exact_solution + * arguments to integrate_difference()). Let $e_c = \hat{f}_c - f_c$ + * be the difference or error between the two. Further, + * let $w:\Omega \rightarrow \mathbb{R}^c$ be the @p weight function of integrate_difference(), which is + * assumed to be equal to one if not supplied. Finally, let $p$ be the + * @p exponent argument (for $L_p$-norms). + * + * In the following,we denote by $E_K$ the local error computed by + * integrate_difference() on cell $K$, whereas $E$ is the global error + * computed by compute_global_error(). Note that integrals are + * approximated by quadrature in the usual way: + * @f[ + * \int_A f(x) dx \approx \sum_q f(x_q) \omega_q. + * @f] + * Similarly for suprema over a cell $T$: + * @f[ + * \sup_{x\in T} |f(x)| dx \approx \max_q |f(x_q)|. + * @f] */ enum NormType { /** - * The function or difference of functions is integrated on each cell. + * The function or difference of functions is integrated on each cell $K$: + * @f[ + * E_K + * = \int_K \sum_c (\hat{f}_c - f_c) \, w_c + * = \int_K \sum_c e_c \, w_c + * @f] + * and summed up to get + * @f[ + * E = \sum_K E_K + * = \int_\Omega \sum_c (\hat{f}_c - f_c) \, w_c + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \int_\Omega (\hat{f} - f) + * = \int_\Omega e. + * @f] + * + * Note: This differs from what is typically known as + * the mean of a function by a factor of $\frac{1}{|\Omega|}$. To + * compute the mean you can also use compute_mean_value(). Finally, + * pay attention to the sign: if $\hat{f}=0$, this will compute the + * negative of the mean of $f$. */ mean, + /** - * The absolute value of the function is integrated. + * The absolute value of the function is integrated: + * @f[ + * E_K = \int_K \sum_c |e_c| \, w_c + * @f] + * and + * @f[ + * E = \sum_K E_K = \int_\Omega \sum_c |e_c| w_c, + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \| e \|_{L^1}. + * @f] */ L1_norm, + /** * The square of the function is integrated and the the square root of the - * result is computed on each cell. + * result is computed on each cell: + * @f[ + * E_K = \sqrt{ \int_K \sum_c e_c^2 \, w_c } + * @f] + * and + * @f[ + * E = \sqrt{\sum_K E_K^2} = \sqrt{ \int_\Omega \sum_c e_c^2 \, w_c } + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \sqrt{ \int_\Omega e^2 } + * = \| e \|_{L^2} + * @f] */ L2_norm, + /** - * The absolute value to the pth power is integrated and the pth - * root is computed on each cell. The exponent p is the last - * parameter of the function. + * The absolute value to the $p$-th power is integrated and the $p$-th + * root is computed on each cell. The exponent $p$ is the @p + * exponent argument of integrate_difference() and compute_mean_value(): + * @f[ + * E_K = \left( \int_K \sum_c |e_c|^p \, w_c \right)^{1/p} + * @f] + * and + * @f[ + * E = \left( \sum_K E_K^p \right)^{1/p} + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \| e \|_{L^p}. + * @f] */ Lp_norm, + /** - * The maximum absolute value of the function. + * The maximum absolute value of the function: + * @f[ + * E_K = \sup_K \max_c |e_c| \, w_c + * @f] + * and + * @f[ + * E = \max_K E_K + * = \sup_\Omega \max_c |e_c| \, w_c + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \sup_\Omega \|e\|_\infty = \| e \|_{L^\infty}. + * @f] */ Linfty_norm, + /** - * #L2_norm of the gradient. + * #L2_norm of the gradient: + * @f[ + * E_K = \sqrt{ \int_K \sum_c (\nabla e_c)^2 \, w_c } + * @f] + * and + * @f[ + * E = \sqrt{\sum_K E_K^2} = \sqrt{ \int_\Omega \sum_c (\nabla e_c)^2 \, w_c } + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \| \nabla e \|_{L^2}. + * @f] */ H1_seminorm, + /** - * #L2_norm of the divergence of a vector field + * #L2_norm of the divergence of a vector field. The function $f$ is + * expected to have $c \geq \text{dim}$ components and the first @p dim + * will be used to compute the divergence: + * @f[ + * E_K = \sqrt{ \int_K \left( \sum_c \frac{\partial e_c}{\partial x_c} \, \sqrt{w_c} \right)^2 } + * @f] + * and + * @f[ + * E = \sqrt{\sum_K E_K^2} + * = \sqrt{ \int_\Omega \left( \sum_c \frac{\partial e_c}{\partial x_c} \, \sqrt{w_c} \right)^2 } + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \| \nabla \cdot e \|_{L^2}. + * @f] */ Hdiv_seminorm, + /** * The square of this norm is the square of the #L2_norm plus the square - * of the #H1_seminorm. + * of the #H1_seminorm: + * @f[ + * E_K = \sqrt{ \int_K \sum_c (e_c^2 + (\nabla e_c)^2) \, w_c } + * @f] + * and + * @f[ + * E = \sqrt{\sum_K E_K^2} = \sqrt{ \int_\Omega \sum_c (e_c^2 + (\nabla e_c)^2) \, w_c } + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \left( \| e \|_{L^2}^2 + \| \nabla e \|_{L^2}^2 \right)^{1/2}. + * @f] */ H1_norm, + /** - * #Lp_norm of the gradient. + * #Lp_norm of the gradient: + * @f[ + * E_K = \left( \int_K \sum_c |\nabla e_c|^p \, w_c \right)^{1/p} + * @f] + * and + * @f[ + * E = \left( \sum_K E_K^p \right)^{1/p} + * = \left( \int_\Omega \sum_c |\nabla e_c|^p \, w_c \right)^{1/p} + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \| \nabla e \|_{L^p}. + * @f] */ W1p_seminorm, + /** - * same as #H1_norm for Lp. + * The same as the #H1_norm but using Lp: + * @f[ + * E_K = \left( \int_K \sum_c (|e_c|^p + |\nabla e_c|^p) \, w_c \right)^{1/p} + * @f] + * and + * @f[ + * E = \left( \sum_K E_K^p \right)^{1/p} + * = \left( \int_\Omega \sum_c (|e_c|^p + |\nabla e_c|^p) \, w_c \right)^{1/p} + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \left( \| e \|_{L^p}^p + \| \nabla e \|_{L^p}^p \right)^{1/p}. + * @f] */ W1p_norm, + /** - * #Linfty_norm of the gradient. + * #Linfty_norm of the gradient: + * @f[ + * E_K = \sup_K \max_c |\nabla e_c| \, w_c + * @f] + * and + * @f[ + * E = \max_K E_K + * = \sup_\Omega \max_c |\nabla e_c| \, w_c + * @f] + * or, for $w \equiv 1$: + * @f[ + * E = \| \nabla e \|_{L^\infty}. + * @f] + * */ W1infty_seminorm, + /** - * same as #H1_norm for Linfty. + * The sum of #Linfty_norm and #W1infty_seminorm: + * @f[ + * E_K = \sup_K \max_c |e_c| \, w_c + \sup_K \max_c |\nabla e_c| \, w_c. + * @f] + * The global norm is not implemented in compute_global_error(), + * because it is impossible to compute the sum of the global + * norms from the values $E_K$. As a work-around, you can compute the + * global #Linfty_norm and #W1infty_seminorm separately and then add them + * to get (with $w \equiv 1$): + * @f[ + * E = \| e \|_{L^\infty} + \| \nabla e \|_{L^\infty}. + * @f] */ W1infty_norm @@ -1831,7 +2015,7 @@ namespace VectorTools //@{ /** - * Compute the error of the finite element solution. Integrate the + * Compute the cellwise error of the finite element solution. Integrate the * difference between a reference function which is given as a continuous * function object, and a finite element function. The result of this * function is the vector @p difference that contains one value per active @@ -1845,6 +2029,10 @@ namespace VectorTools * It is assumed that the number of components of the function @p * exact_solution matches that of the finite element used by @p dof. * + * To compute a global error norm of a finite element solution, use + * VectorTools::compute_global_error() with the output vector computed with + * this function. + * * @param[in] mapping The mapping that is used when integrating the * difference $u-u_h$. * @param[in] dof The DoFHandler object that describes the finite element @@ -1886,8 +2074,9 @@ namespace VectorTools * function, a null pointer, is interpreted as "no weighting function", * i.e., weight=1 in the whole domain for all vector components uniformly. * @param[in] exponent This value denotes the $p$ used in computing - * $L^p$-norms and $W^{1,p}$-norms. The value is ignores if a @p norm other - * than NormType::Lp_norm or NormType::W1p_norm is chosen. + * $L^p$-norms and $W^{1,p}$-norms. The value is ignored if a @p norm other + * than NormType::Lp_norm, NormType::W1p_norm, or NormType::W1p_seminorm + * is chosen. * * * See the general documentation of this namespace for more information. @@ -1904,38 +2093,8 @@ namespace VectorTools * The vector computed will, in the case of a distributed triangulation, * contain zeros for cells that are not locally owned. As a consequence, in * order to compute the global $L_2$ error (for example), the errors - * from different processors need to be combined, but this is simple because - * every processor only computes contributions for those cells of the global - * triangulation it locally owns (and these sets are, by definition, - * mutually disjoint). Consequently, the following piece of code computes - * the global $L_2$ error across multiple processors sharing a - * parallel::distribute::Triangulation: - * @code - * Vector local_errors (tria.n_active_cells()); - * VectorTools::integrate_difference (mapping, dof, - * solution, exact_solution, - * local_errors, - * QGauss(fe.degree+2), - * VectorTools::L2_norm); - * const double total_local_error = local_errors.l2_norm(); - * const double total_global_error - * = std::sqrt (Utilities::MPI::sum (total_local_error * total_local_error, MPI_COMM_WORLD)); - * @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 the $L_1$ norm of the - * error, the correct form of the last two lines would have been - * @code - * const double total_local_error = local_errors.l1_norm(); - * const double total_global_error - * = Utilities::MPI::sum (total_local_error, MPI_COMM_WORLD); - * @endcode - * instead, and similar considerations hold when computing the $L_\infty$ - * norm of the error. + * from different processors need to be combined, see + * VectorTools::compute_global_error(). * * Instantiations for this template are provided for some vector types (see * the general documentation of the namespace), but only for InVectors as in @@ -1995,6 +2154,37 @@ namespace VectorTools const Function *weight = 0, const double exponent = 2.); + /** + * Take a Vector @p cellwise_error of errors on each cell with + * tria.n_active_cells() entries and return the global + * error as given by @p norm. + * + * The @p cellwise_error vector is typically an output produced by + * VectorTools::integrate_difference() and you normally want to supply the + * same value for @p norm as you used in VectorTools::integrate_difference(). + * + * If the given Triangulation is a parallel::Triangulation, entries + * in @p cellwise_error that do not correspond to locally owned cells are + * assumed to be 0.0 and a parallel reduction using MPI is done to compute + * the global error. + * + * @param tria The Triangulation with active cells corresponding with the + * entries in @p cellwise_error. + * @param cellwise_error Vector of errors on each active cell. + * @param norm The type of norm to compute. + * @param exponent The exponent $p$ to use for $L^p$-norms and + * $W^{1,p}$-norms. The value is ignored if a @p norm other + * than NormType::Lp_norm, NormType::W1p_norm, or NormType::W1p_seminorm + * is chosen. + * + * @note Instantiated for type Vector and Vector. + */ + template + double compute_global_error(const Triangulation &tria, + const InVector &cellwise_error, + const NormType &norm, + const double exponent = 2.); + /** * Point error evaluation. Find the first cell containing the given point * and compute the difference of a (possibly vector-valued) finite element diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index e1f47314f4..3bc45d8325 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -6662,7 +6662,97 @@ namespace VectorTools norm, weight, exponent); } + template + double compute_global_error(const Triangulation &tria, + const InVector &cellwise_error, + const NormType &norm, + const double exponent) + { + Assert( cellwise_error.size() == tria.n_active_cells(), + ExcMessage("input vector cell_error has invalid size!")); +#ifdef DEBUG + { + // check that off-processor entries are zero. Otherwise we will compute + // wrong results below! + typename InVector::size_type i = 0; + typename Triangulation::cell_iterator it = tria.begin_active(); + for (; iis_locally_owned()) + Assert(std::fabs(cellwise_error[i]) < 1e-20, + ExcMessage("cellwise_error of cells that are not locally owned need to be zero!")); + } +#endif + + MPI_Comm comm = MPI_COMM_SELF; +#ifdef DEAL_II_WITH_MPI + if (const parallel::Triangulation *ptria = + dynamic_cast*>(&tria)) + comm = ptria->get_communicator(); +#endif + switch (norm) + { + case L2_norm: + case H1_seminorm: + case H1_norm: + case Hdiv_seminorm: + { + const double local = cellwise_error.l2_norm(); + return std::sqrt (Utilities::MPI::sum (local * local, comm)); + } + + case L1_norm: + { + const double local = cellwise_error.l1_norm(); + return Utilities::MPI::sum (local, comm); + } + + case Linfty_norm: + case W1infty_seminorm: + { + const double local = cellwise_error.linfty_norm(); + return Utilities::MPI::max (local, comm); + } + + case W1infty_norm: + { + AssertThrow(false, ExcMessage( + "compute_global_error() is impossible for " + "the W1infty_norm. See the documentation for " + "NormType::W1infty_norm for more information.")); + return std::numeric_limits::infinity(); + } + + case mean: + { + // Note: mean is defined as int_\Omega f = sum_K \int_K f, so we need + // the sum of the cellwise errors not the Euclidean mean value that + // is returned by Vector<>::mean_value(). + const double local = cellwise_error.mean_value() + * cellwise_error.size(); + return Utilities::MPI::sum (local, comm); + } + + case Lp_norm: + case W1p_norm: + case W1p_seminorm: + { + double local = 0; + typename InVector::size_type i; + typename Triangulation::cell_iterator it = tria.begin_active(); + for (i = 0; iis_locally_owned()) + local += std::pow(cellwise_error[i], exponent); + + return std::pow (Utilities::MPI::sum (local, comm), 1./exponent); + } + + default: + AssertThrow(false, ExcNotImplemented()); + break; + } + return 0.0; + } template void diff --git a/source/numerics/vector_tools_integrate_difference.inst.in b/source/numerics/vector_tools_integrate_difference.inst.in index 690974ffb6..7783007063 100644 --- a/source/numerics/vector_tools_integrate_difference.inst.in +++ b/source/numerics/vector_tools_integrate_difference.inst.in @@ -114,3 +114,24 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens \} #endif } + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + namespace VectorTools \{ + template + double compute_global_error >( + const Triangulation &, + const Vector &, + const NormType &, + const double); + + template + double compute_global_error >( + const Triangulation &, + const Vector &, + const NormType &, + const double); + \} +#endif + } diff --git a/tests/mpi/integrate_difference.cc b/tests/mpi/integrate_difference.cc index a181c25df9..cc9bf9d7d4 100644 --- a/tests/mpi/integrate_difference.cc +++ b/tests/mpi/integrate_difference.cc @@ -87,15 +87,10 @@ void test() results, QGauss(3), VectorTools::L2_norm); - double local = results.l2_norm() * results.l2_norm(); - double global; - - MPI_Allreduce (&local, &global, 1, MPI_DOUBLE, - MPI_SUM, - tr.get_communicator()); + double global = VectorTools::compute_global_error(tr, results, VectorTools::L2_norm); if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) - deallog << "difference = " << std::sqrt(global) + deallog << "difference = " << global << std::endl; // we have f(\vec x)=x, so the difference @@ -105,7 +100,7 @@ void test() // note that we have used a quadrature // formula of sufficient order to get exact // results - Assert (std::fabs(std::sqrt(global) - 1./std::sqrt(3.)) < 1e-6, + Assert (std::fabs(global - 1./std::sqrt(3.)) < 1e-6, ExcInternalError()); } diff --git a/tests/vector_tools/integrate_difference_02.cc b/tests/vector_tools/integrate_difference_02.cc new file mode 100644 index 0000000000..9ca4376744 --- /dev/null +++ b/tests/vector_tools/integrate_difference_02.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test integrate_difference and compute_global_error + +// see http://www.wolframalpha.com/input/?i=integrate+(x%2By%2Bz)%5E3%2B(x%5E2%2By%5E2)%5E3%2B(z%2Bxy)%5E3+from+x%3D0..1,y%3D0..1,z%3D0..1 + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +// x+y+z, x^2+y^2, z+xy +// div = 1+2y+1 +template +class Ref : public Function +{ +public: + Ref() + :Function(dim) + {} + + double value (const Point &p, const unsigned int c) const + { + if (c==0) + return p[0]+p[1]+((dim==3)?p[2]:0.0); + if (c==1) + return p[0]*p[0]+p[1]*p[1]; + if (c==2) + return p[2]+p[0]*p[1]; + return 0.0; + } +}; + + + +template +void test(VectorTools::NormType norm, double value, double exp = 2.0) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FESystem fe(FE_Q(4),dim); + DoFHandler dofh(tria); + dofh.distribute_dofs(fe); + + Vector solution (dofh.n_dofs ()); + VectorTools::interpolate(dofh, Ref(), solution); + + Vector cellwise_errors (tria.n_active_cells()); + QIterated quadrature (QTrapez<1>(), 5); + + const dealii::Function *w = 0; + VectorTools::integrate_difference (dofh, + solution, + ZeroFunction(dim), + cellwise_errors, + quadrature, + norm, + w, + exp); + + const double error + = VectorTools::compute_global_error(tria, cellwise_errors, norm, exp); + + const double difference = std::abs(error-value); + deallog << "computed: " << error + << " expected: " << value + << " difference: " << difference + << std::endl; + Assert(difference<2e-3, ExcMessage("Error in integrate_difference")); +} + +template +void test() +{ + deallog << "Hdiv_seminorm:" << std::endl; + // sqrt(\int (div f)^2 = sqrt(\int (1+2y+1)^2) + test(VectorTools::Hdiv_seminorm, 2.0*std::sqrt(7.0/3.0)); + + deallog << "L2_norm:" << std::endl; + // sqrt(\int_\Omega f^2) = sqrt(\int (x+y+z)^2+(x^2+y^2)^2+(z+xy)^2) + test(VectorTools::L2_norm, std::sqrt(229.0/60.0)); + + deallog << "H1_seminorm:" << std::endl; + // sqrt( \int sum_k | d/dxi f_k |_0^2 ) + // = sqrt( \int 3+ (2x)^2+(2y)^2 + y^2+x^2+1 + // = sqrt( 22/3 ) + test(VectorTools::H1_seminorm, std::sqrt(22.0/3.0)); + + deallog << "H1_norm:" << std::endl; + test(VectorTools::H1_norm, std::sqrt(229.0/60.0+22.0/3.0)); + + deallog << "L1_norm:" << std::endl; + test(VectorTools::L1_norm, 35.0/12.0); + + deallog << "Linfty_norm:" << std::endl; + test(VectorTools::Linfty_norm, 3.0); + + deallog << "mean:" << std::endl; + // int -(x+y+z + x^2+y^2 + z+xy) + test(VectorTools::mean, -35.0/12.0); + + deallog << "Lp_norm:" << std::endl; + // (int (x+y+z)^p+(x^2+y^2)^p+(z+xy)^p) ) ^ 1./p + // = std::pow(9937.0/1680.0, 1.0/3.0) + test(VectorTools::Lp_norm, std::pow(9937.0/1680.0, 1.0/3.0), 3.0); + + deallog << "W1p_seminorm:" << std::endl; + // ( \int_K sum_k | d/dxi f_k |_0^2^(p/2) )^1/p + // ( integrate 3^(3/2) + ((2x)^2+(2y)^2)^(3/2) + (y^2+x^2+1)^(3/2) ) ^1/p + // = (12.4164) ^1/3 = 2.31560 + test(VectorTools::W1p_seminorm, 2.31560, 3.0); + + deallog << "OK" << std::endl; +} + + +int main (int argc, char **argv) +{ + initlog(); + test<3>(); +} diff --git a/tests/vector_tools/integrate_difference_02.output b/tests/vector_tools/integrate_difference_02.output new file mode 100644 index 0000000000..f666b9bfd9 --- /dev/null +++ b/tests/vector_tools/integrate_difference_02.output @@ -0,0 +1,20 @@ + +DEAL::Hdiv_seminorm: +DEAL::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL::L2_norm: +DEAL::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL::H1_seminorm: +DEAL::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL::H1_norm: +DEAL::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL::L1_norm: +DEAL::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL::Linfty_norm: +DEAL::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL::mean: +DEAL::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL::Lp_norm: +DEAL::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL::W1p_seminorm: +DEAL::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL::OK diff --git a/tests/vector_tools/integrate_difference_03.cc b/tests/vector_tools/integrate_difference_03.cc new file mode 100644 index 0000000000..40a7673539 --- /dev/null +++ b/tests/vector_tools/integrate_difference_03.cc @@ -0,0 +1,157 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test integrate_difference and compute_global_error in parallel +// see integrate_difference_02.cc for the serial version + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +// x+y+z, x^2+y^2, z+xy +// div = 1+2y+1 +template +class Ref : public Function +{ +public: + Ref() + :Function(dim) + {} + + double value (const Point &p, const unsigned int c) const + { + if (c==0) + return p[0]+p[1]+((dim==3)?p[2]:0.0); + if (c==1) + return p[0]*p[0]+p[1]*p[1]; + if (c==2) + return p[2]+p[0]*p[1]; + return 0.0; + } +}; + + + +template +void test(VectorTools::NormType norm, double value, double exp = 2.0) +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FESystem fe(FE_Q(4),dim); + DoFHandler dofh(tria); + dofh.distribute_dofs(fe); + + TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(), + MPI_COMM_WORLD); + + VectorTools::interpolate(dofh, Ref(), interpolated); + + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dofh, relevant_set); + TrilinosWrappers::MPI::Vector solution(relevant_set, MPI_COMM_WORLD); + solution = interpolated; + + Vector cellwise_errors (tria.n_active_cells()); + QIterated quadrature (QTrapez<1>(), 5); + + const dealii::Function *w = 0; + VectorTools::integrate_difference (dofh, + solution, + ZeroFunction(dim), + cellwise_errors, + quadrature, + norm, + w, + exp); + + const double error + = VectorTools::compute_global_error(tria, cellwise_errors, norm, exp); + + const double difference = std::abs(error-value); + deallog << "computed: " << error + << " expected: " << value + << " difference: " << difference + << std::endl; + Assert(difference<2e-3, ExcMessage("Error in integrate_difference")); +} + +template +void test() +{ + deallog << "Hdiv_seminorm:" << std::endl; + // sqrt(\int (div f)^2 = sqrt(\int (1+2y+1)^2) + test(VectorTools::Hdiv_seminorm, 2.0*std::sqrt(7.0/3.0)); + + deallog << "L2_norm:" << std::endl; + // sqrt(\int_\Omega f^2) = sqrt(\int (x+y+z)^2+(x^2+y^2)^2+(z+xy)^2) + test(VectorTools::L2_norm, std::sqrt(229.0/60.0)); + + deallog << "H1_seminorm:" << std::endl; + // sqrt( \int sum_k | d/dxi f_k |_0^2 ) + // = sqrt( \int 3+ (2x)^2+(2y)^2 + y^2+x^2+1 + // = sqrt( 22/3 ) + test(VectorTools::H1_seminorm, std::sqrt(22.0/3.0)); + + deallog << "H1_norm:" << std::endl; + test(VectorTools::H1_norm, std::sqrt(229.0/60.0+22.0/3.0)); + + deallog << "L1_norm:" << std::endl; + test(VectorTools::L1_norm, 35.0/12.0); + + deallog << "Linfty_norm:" << std::endl; + test(VectorTools::Linfty_norm, 3.0); + + deallog << "mean:" << std::endl; + // int -(x+y+z + x^2+y^2 + z+xy) + test(VectorTools::mean, -35.0/12.0); + + deallog << "Lp_norm:" << std::endl; + // (int (x+y+z)^p+(x^2+y^2)^p+(z+xy)^p) ) ^ 1./p + // = std::pow(9937.0/1680.0, 1.0/3.0) + test(VectorTools::Lp_norm, std::pow(9937.0/1680.0, 1.0/3.0), 3.0); + + deallog << "W1p_seminorm:" << std::endl; + // ( \int_K sum_k | d/dxi f_k |_0^2^(p/2) )^1/p + // ( integrate 3^(3/2) + ((2x)^2+(2y)^2)^(3/2) + (y^2+x^2+1)^(3/2) ) ^1/p + // = (12.4164) ^1/3 = 2.31560 + test(VectorTools::W1p_seminorm, 2.31560, 3.0); + + deallog << "OK" << std::endl; +} + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + test<3>(); +} diff --git a/tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=true.mpirun=1.output b/tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..d8f0ac1c2a --- /dev/null +++ b/tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,20 @@ + +DEAL:0::Hdiv_seminorm: +DEAL:0::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:0::L2_norm: +DEAL:0::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:0::H1_seminorm: +DEAL:0::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:0::H1_norm: +DEAL:0::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:0::L1_norm: +DEAL:0::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:0::Linfty_norm: +DEAL:0::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:0::mean: +DEAL:0::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:0::Lp_norm: +DEAL:0::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:0::W1p_seminorm: +DEAL:0::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:0::OK diff --git a/tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=truempirun=3.output b/tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=truempirun=3.output new file mode 100644 index 0000000000..630918c014 --- /dev/null +++ b/tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=truempirun=3.output @@ -0,0 +1,62 @@ + +DEAL:0::Hdiv_seminorm: +DEAL:0::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:0::L2_norm: +DEAL:0::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:0::H1_seminorm: +DEAL:0::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:0::H1_norm: +DEAL:0::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:0::L1_norm: +DEAL:0::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:0::Linfty_norm: +DEAL:0::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:0::mean: +DEAL:0::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:0::Lp_norm: +DEAL:0::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:0::W1p_seminorm: +DEAL:0::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:0::OK + +DEAL:1::Hdiv_seminorm: +DEAL:1::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:1::L2_norm: +DEAL:1::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:1::H1_seminorm: +DEAL:1::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:1::H1_norm: +DEAL:1::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:1::L1_norm: +DEAL:1::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:1::Linfty_norm: +DEAL:1::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:1::mean: +DEAL:1::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:1::Lp_norm: +DEAL:1::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:1::W1p_seminorm: +DEAL:1::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:1::OK + + +DEAL:2::Hdiv_seminorm: +DEAL:2::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:2::L2_norm: +DEAL:2::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:2::H1_seminorm: +DEAL:2::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:2::H1_norm: +DEAL:2::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:2::L1_norm: +DEAL:2::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:2::Linfty_norm: +DEAL:2::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:2::mean: +DEAL:2::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:2::Lp_norm: +DEAL:2::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:2::W1p_seminorm: +DEAL:2::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:2::OK + diff --git a/tests/vector_tools/integrate_difference_04.cc b/tests/vector_tools/integrate_difference_04.cc new file mode 100644 index 0000000000..4ae2b2e72a --- /dev/null +++ b/tests/vector_tools/integrate_difference_04.cc @@ -0,0 +1,157 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test integrate_difference and compute_global_error in parallel +// see integrate_difference_02.cc for the serial version + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +// x+y+z, x^2+y^2, z+xy +// div = 1+2y+1 +template +class Ref : public Function +{ +public: + Ref() + :Function(dim) + {} + + double value (const Point &p, const unsigned int c) const + { + if (c==0) + return p[0]+p[1]+((dim==3)?p[2]:0.0); + if (c==1) + return p[0]*p[0]+p[1]*p[1]; + if (c==2) + return p[2]+p[0]*p[1]; + return 0.0; + } +}; + + + +template +void test(VectorTools::NormType norm, double value, double exp = 2.0) +{ + parallel::shared::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FESystem fe(FE_Q(4),dim); + DoFHandler dofh(tria); + dofh.distribute_dofs(fe); + + TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(), + MPI_COMM_WORLD); + + VectorTools::interpolate(dofh, Ref(), interpolated); + + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dofh, relevant_set); + TrilinosWrappers::MPI::Vector solution(relevant_set, MPI_COMM_WORLD); + solution = interpolated; + + Vector cellwise_errors (tria.n_active_cells()); + QIterated quadrature (QTrapez<1>(), 5); + + const dealii::Function *w = 0; + VectorTools::integrate_difference (dofh, + solution, + ZeroFunction(dim), + cellwise_errors, + quadrature, + norm, + w, + exp); + + const double error + = VectorTools::compute_global_error(tria, cellwise_errors, norm, exp); + + const double difference = std::abs(error-value); + deallog << "computed: " << error + << " expected: " << value + << " difference: " << difference + << std::endl; + Assert(difference<2e-3, ExcMessage("Error in integrate_difference")); +} + +template +void test() +{ + deallog << "Hdiv_seminorm:" << std::endl; + // sqrt(\int (div f)^2 = sqrt(\int (1+2y+1)^2) + test(VectorTools::Hdiv_seminorm, 2.0*std::sqrt(7.0/3.0)); + + deallog << "L2_norm:" << std::endl; + // sqrt(\int_\Omega f^2) = sqrt(\int (x+y+z)^2+(x^2+y^2)^2+(z+xy)^2) + test(VectorTools::L2_norm, std::sqrt(229.0/60.0)); + + deallog << "H1_seminorm:" << std::endl; + // sqrt( \int sum_k | d/dxi f_k |_0^2 ) + // = sqrt( \int 3+ (2x)^2+(2y)^2 + y^2+x^2+1 + // = sqrt( 22/3 ) + test(VectorTools::H1_seminorm, std::sqrt(22.0/3.0)); + + deallog << "H1_norm:" << std::endl; + test(VectorTools::H1_norm, std::sqrt(229.0/60.0+22.0/3.0)); + + deallog << "L1_norm:" << std::endl; + test(VectorTools::L1_norm, 35.0/12.0); + + deallog << "Linfty_norm:" << std::endl; + test(VectorTools::Linfty_norm, 3.0); + + deallog << "mean:" << std::endl; + // int -(x+y+z + x^2+y^2 + z+xy) + test(VectorTools::mean, -35.0/12.0); + + deallog << "Lp_norm:" << std::endl; + // (int (x+y+z)^p+(x^2+y^2)^p+(z+xy)^p) ) ^ 1./p + // = std::pow(9937.0/1680.0, 1.0/3.0) + test(VectorTools::Lp_norm, std::pow(9937.0/1680.0, 1.0/3.0), 3.0); + + deallog << "W1p_seminorm:" << std::endl; + // ( \int_K sum_k | d/dxi f_k |_0^2^(p/2) )^1/p + // ( integrate 3^(3/2) + ((2x)^2+(2y)^2)^(3/2) + (y^2+x^2+1)^(3/2) ) ^1/p + // = (12.4164) ^1/3 = 2.31560 + test(VectorTools::W1p_seminorm, 2.31560, 3.0); + + deallog << "OK" << std::endl; +} + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + test<3>(); +} diff --git a/tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=1.output b/tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=1.output new file mode 100644 index 0000000000..d8f0ac1c2a --- /dev/null +++ b/tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=1.output @@ -0,0 +1,20 @@ + +DEAL:0::Hdiv_seminorm: +DEAL:0::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:0::L2_norm: +DEAL:0::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:0::H1_seminorm: +DEAL:0::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:0::H1_norm: +DEAL:0::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:0::L1_norm: +DEAL:0::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:0::Linfty_norm: +DEAL:0::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:0::mean: +DEAL:0::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:0::Lp_norm: +DEAL:0::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:0::W1p_seminorm: +DEAL:0::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:0::OK diff --git a/tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=3.output b/tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=3.output new file mode 100644 index 0000000000..630918c014 --- /dev/null +++ b/tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=3.output @@ -0,0 +1,62 @@ + +DEAL:0::Hdiv_seminorm: +DEAL:0::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:0::L2_norm: +DEAL:0::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:0::H1_seminorm: +DEAL:0::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:0::H1_norm: +DEAL:0::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:0::L1_norm: +DEAL:0::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:0::Linfty_norm: +DEAL:0::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:0::mean: +DEAL:0::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:0::Lp_norm: +DEAL:0::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:0::W1p_seminorm: +DEAL:0::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:0::OK + +DEAL:1::Hdiv_seminorm: +DEAL:1::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:1::L2_norm: +DEAL:1::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:1::H1_seminorm: +DEAL:1::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:1::H1_norm: +DEAL:1::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:1::L1_norm: +DEAL:1::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:1::Linfty_norm: +DEAL:1::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:1::mean: +DEAL:1::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:1::Lp_norm: +DEAL:1::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:1::W1p_seminorm: +DEAL:1::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:1::OK + + +DEAL:2::Hdiv_seminorm: +DEAL:2::computed: 3.05532 expected: 3.05505 difference: 0.000272760 +DEAL:2::L2_norm: +DEAL:2::computed: 1.95470 expected: 1.95363 difference: 0.00106613 +DEAL:2::H1_seminorm: +DEAL:2::computed: 2.70878 expected: 2.70801 difference: 0.000769213 +DEAL:2::H1_norm: +DEAL:2::computed: 3.34041 expected: 3.33916 difference: 0.00124760 +DEAL:2::L1_norm: +DEAL:2::computed: 2.91750 expected: 2.91667 difference: 0.000833333 +DEAL:2::Linfty_norm: +DEAL:2::computed: 3.00000 expected: 3.00000 difference: 0.00000 +DEAL:2::mean: +DEAL:2::computed: -2.91750 expected: -2.91667 difference: 0.000833333 +DEAL:2::Lp_norm: +DEAL:2::computed: 1.80970 expected: 1.80849 difference: 0.00121796 +DEAL:2::W1p_seminorm: +DEAL:2::computed: 2.31644 expected: 2.31560 difference: 0.000840093 +DEAL:2::OK +