]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement VectorTools::compute_global_error, document norms 2602/head
authorTimo Heister <timo.heister@gmail.com>
Sun, 15 May 2016 15:16:18 +0000 (16:16 +0100)
committerTimo Heister <timo.heister@gmail.com>
Sun, 22 May 2016 18:45:42 +0000 (19:45 +0100)
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.

20 files changed:
doc/news/changes.h
examples/step-11/step-11.cc
examples/step-20/step-20.cc
examples/step-34/step-34.cc
examples/step-38/step-38.cc
examples/step-48/step-48.cc
examples/step-51/step-51.cc
examples/step-7/step-7.cc
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_integrate_difference.inst.in
tests/mpi/integrate_difference.cc
tests/vector_tools/integrate_difference_02.cc [new file with mode: 0644]
tests/vector_tools/integrate_difference_02.output [new file with mode: 0644]
tests/vector_tools/integrate_difference_03.cc [new file with mode: 0644]
tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/vector_tools/integrate_difference_03.with_trilinos=true.with_p4est=truempirun=3.output [new file with mode: 0644]
tests/vector_tools/integrate_difference_04.cc [new file with mode: 0644]
tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=1.output [new file with mode: 0644]
tests/vector_tools/integrate_difference_04.with_trilinos=true.with_metis=true.mpirun=3.output [new file with mode: 0644]

index 88ce9b37a2a32360f71b6d6ed56e22f4b0cd5dbb..08159b146423b88b626ca26ca487027fba000b0c 100644 (file)
@@ -120,6 +120,13 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
+ <li> New: Add VectorTools::compute_global_error that computes global
+ errors from cellwise errors obtained by VectorTools::integrate_difference()
+ and do MPI collectives if necessary.
+ <br>
+ (Timo Heister, 2016/05/15)
+ </li>
+
  <li> New: Add functions to transform Cartesian coordinates to spherical and back:
  GeometricUtilities::Coordinates::to_spherical and
  GeometricUtilities::Coordinates::from_spherical.
index 8b392417869944d06021a386d8ab923065d19270..ea9b17a48f8cf427d754b4c64685e0a875149238 100644 (file)
@@ -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());
index fb1be4ff9567a60f7bf58e19fc62107b91236dad..72bb761cb08918e2e72b3c55e446fc95029066b0 100644 (file)
@@ -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
-  // <code>VectorTools::integrate_difference</code> 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
+  // <code>VectorTools::integrate_difference</code> 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
index dd54cbb4d73c5e2bdbb638fbd60d6d7ad450dbfc..87e4caa76381fab5bf84699a4d47cbe146d6c631 100644 (file)
@@ -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
index 70a1958ee2f51d318786c2419a20f75ca34128a5..9035706e217d453650ec032024f4a87140bbe16f 100644 (file)
@@ -525,8 +525,11 @@ namespace Step38
                                        QGauss<dim>(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;
   }
 
index 24e7c5694b81b21086010d63675ab35f09db26a0..df18653e277b2ad0360992cc620845e7ec28da7b 100644 (file)
@@ -467,7 +467,9 @@ namespace Step48
                                        QGauss<dim>(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
index 80434c891750fcb6c4f69655004b0cd44acd2c17..aa2d32a578a6210d239649702da381dace5dd6fb 100644 (file)
@@ -1052,7 +1052,9 @@ namespace Step51
                                        QGauss<dim>(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<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
                                                   dim+1);
@@ -1063,7 +1065,9 @@ namespace Step51
                                        QGauss<dim>(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<dim>(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());
index 69a19b9a97d0e3d9acfc2b74b52412766316f12b..f1c4a197ac7210747396c0d070934585b8430d5d 100644 (file)
@@ -867,22 +867,26 @@ namespace Step7
                                        difference_per_cell,
                                        QGauss<dim>(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
     // <code>difference_per_cell</code> vector since it is no longer used
     // after computing the <code>L2_error</code> 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<dim>(),
                                        difference_per_cell,
                                        QGauss<dim>(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<dim> 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
index c617a7ad7ba838321c7a37563e95da0e89f4de95..a6e8638a25f65249f3c0234c1a1602f03636345a 100644 (file)
@@ -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 <i>p</i>th power is integrated and the pth
-     * root is computed on each cell. The exponent <i>p</i> 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 <i>L<sup>p</sup></i>.
+     * The same as the #H1_norm but using <i>L<sup>p</sup></i>:
+     * @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 <i>L<sup>infty</sup></i>.
+     * 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 <i>global</i> $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<double> local_errors (tria.n_active_cells());
-   *    VectorTools::integrate_difference (mapping, dof,
-   *                                       solution, exact_solution,
-   *                                       local_errors,
-   *                                       QGauss<dim>(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<spacedim,double>           *weight = 0,
                              const double exponent = 2.);
 
+  /**
+   * Take a Vector @p cellwise_error of errors on each cell with
+   * <tt>tria.n_active_cells()</tt> 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<double> and Vector<float>.
+   */
+  template <int dim, int spacedim, class InVector>
+  double compute_global_error(const Triangulation<dim,spacedim> &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
index e1f47314f41ed85a24641bfec1cca76dc3f3d2db..3bc45d83253cfb01a1e1108ddc8a46e11af47467 100644 (file)
@@ -6662,7 +6662,97 @@ namespace VectorTools
                               norm, weight, exponent);
   }
 
+  template <int dim, int spacedim, class InVector>
+  double compute_global_error(const Triangulation<dim,spacedim> &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<dim,spacedim>::cell_iterator it = tria.begin_active();
+      for (; i<cellwise_error.size(); ++i, ++it)
+        if (!it->is_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<dim,spacedim> *ptria =
+          dynamic_cast<const parallel::Triangulation<dim,spacedim>*>(&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<double>::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<dim,spacedim>::cell_iterator it = tria.begin_active();
+        for (i = 0; i<cellwise_error.size(); ++i, ++it)
+          if (it->is_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 <int dim, typename VectorType, int spacedim>
   void
index 690974ffb6bbdb605e6803a92e2146802d044c7c..77830070631812c8cf0fbab031aac27a9e157e79 100644 (file)
@@ -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<deal_II_dimension,deal_II_space_dimension,Vector<float> >(
+    const Triangulation<deal_II_dimension,deal_II_space_dimension> &,
+        const Vector<float> &,
+        const NormType &,
+        const double);
+
+  template
+    double compute_global_error<deal_II_dimension,deal_II_space_dimension,Vector<double> >(
+    const Triangulation<deal_II_dimension,deal_II_space_dimension> &,
+        const Vector<double> &,
+        const NormType &,
+        const double);
+  \}
+#endif
+  }
index a181c25df9ba2ebcf4c73812a0a0f71a8c201cb8..cc9bf9d7d44fe2f7f6f0a32287a8173fc7197575 100644 (file)
@@ -87,15 +87,10 @@ void test()
                                      results,
                                      QGauss<dim>(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 (file)
index 0000000..9ca4376
--- /dev/null
@@ -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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+using namespace dealii;
+
+
+// x+y+z, x^2+y^2, z+xy
+// div = 1+2y+1
+template <int dim>
+class Ref : public Function<dim>
+{
+public:
+  Ref()
+    :Function<dim>(dim)
+  {}
+
+  double value (const Point<dim> &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 <int dim>
+void test(VectorTools::NormType norm, double value, double exp = 2.0)
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  FESystem<dim> fe(FE_Q<dim>(4),dim);
+  DoFHandler<dim> dofh(tria);
+  dofh.distribute_dofs(fe);
+
+  Vector<double> solution (dofh.n_dofs ());
+  VectorTools::interpolate(dofh, Ref<dim>(), solution);
+
+  Vector<double> cellwise_errors (tria.n_active_cells());
+  QIterated<dim> quadrature (QTrapez<1>(), 5);
+
+  const dealii::Function<dim,double> *w = 0;
+  VectorTools::integrate_difference (dofh,
+                                     solution,
+                                     ZeroFunction<dim>(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 <int dim>
+void test()
+{
+  deallog << "Hdiv_seminorm:" << std::endl;
+  // sqrt(\int (div f)^2 = sqrt(\int (1+2y+1)^2)
+  test<dim>(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<dim>(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<dim>(VectorTools::H1_seminorm, std::sqrt(22.0/3.0));
+
+  deallog << "H1_norm:" << std::endl;
+  test<dim>(VectorTools::H1_norm, std::sqrt(229.0/60.0+22.0/3.0));
+
+  deallog << "L1_norm:" << std::endl;
+  test<dim>(VectorTools::L1_norm, 35.0/12.0);
+
+  deallog << "Linfty_norm:" << std::endl;
+  test<dim>(VectorTools::Linfty_norm, 3.0);
+
+  deallog << "mean:" << std::endl;
+  // int -(x+y+z + x^2+y^2 + z+xy)
+  test<dim>(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<dim>(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<dim>(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 (file)
index 0000000..f666b9b
--- /dev/null
@@ -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 (file)
index 0000000..40a7673
--- /dev/null
@@ -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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+using namespace dealii;
+
+
+// x+y+z, x^2+y^2, z+xy
+// div = 1+2y+1
+template <int dim>
+class Ref : public Function<dim>
+{
+public:
+  Ref()
+    :Function<dim>(dim)
+  {}
+
+  double value (const Point<dim> &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 <int dim>
+void test(VectorTools::NormType norm, double value, double exp = 2.0)
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  FESystem<dim> fe(FE_Q<dim>(4),dim);
+  DoFHandler<dim> dofh(tria);
+  dofh.distribute_dofs(fe);
+
+  TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(),
+                                             MPI_COMM_WORLD);
+
+  VectorTools::interpolate(dofh, Ref<dim>(), interpolated);
+
+  IndexSet relevant_set;
+  DoFTools::extract_locally_relevant_dofs (dofh, relevant_set);
+  TrilinosWrappers::MPI::Vector solution(relevant_set, MPI_COMM_WORLD);
+  solution = interpolated;
+
+  Vector<double> cellwise_errors (tria.n_active_cells());
+  QIterated<dim> quadrature (QTrapez<1>(), 5);
+
+  const dealii::Function<dim,double> *w = 0;
+  VectorTools::integrate_difference (dofh,
+                                     solution,
+                                     ZeroFunction<dim>(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 <int dim>
+void test()
+{
+  deallog << "Hdiv_seminorm:" << std::endl;
+  // sqrt(\int (div f)^2 = sqrt(\int (1+2y+1)^2)
+  test<dim>(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<dim>(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<dim>(VectorTools::H1_seminorm, std::sqrt(22.0/3.0));
+
+  deallog << "H1_norm:" << std::endl;
+  test<dim>(VectorTools::H1_norm, std::sqrt(229.0/60.0+22.0/3.0));
+
+  deallog << "L1_norm:" << std::endl;
+  test<dim>(VectorTools::L1_norm, 35.0/12.0);
+
+  deallog << "Linfty_norm:" << std::endl;
+  test<dim>(VectorTools::Linfty_norm, 3.0);
+
+  deallog << "mean:" << std::endl;
+  // int -(x+y+z + x^2+y^2 + z+xy)
+  test<dim>(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<dim>(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<dim>(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 (file)
index 0000000..d8f0ac1
--- /dev/null
@@ -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 (file)
index 0000000..630918c
--- /dev/null
@@ -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 (file)
index 0000000..4ae2b2e
--- /dev/null
@@ -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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/distributed/shared_tria.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+using namespace dealii;
+
+
+// x+y+z, x^2+y^2, z+xy
+// div = 1+2y+1
+template <int dim>
+class Ref : public Function<dim>
+{
+public:
+  Ref()
+    :Function<dim>(dim)
+  {}
+
+  double value (const Point<dim> &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 <int dim>
+void test(VectorTools::NormType norm, double value, double exp = 2.0)
+{
+  parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  FESystem<dim> fe(FE_Q<dim>(4),dim);
+  DoFHandler<dim> dofh(tria);
+  dofh.distribute_dofs(fe);
+
+  TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(),
+                                             MPI_COMM_WORLD);
+
+  VectorTools::interpolate(dofh, Ref<dim>(), interpolated);
+
+  IndexSet relevant_set;
+  DoFTools::extract_locally_relevant_dofs (dofh, relevant_set);
+  TrilinosWrappers::MPI::Vector solution(relevant_set, MPI_COMM_WORLD);
+  solution = interpolated;
+
+  Vector<double> cellwise_errors (tria.n_active_cells());
+  QIterated<dim> quadrature (QTrapez<1>(), 5);
+
+  const dealii::Function<dim,double> *w = 0;
+  VectorTools::integrate_difference (dofh,
+                                     solution,
+                                     ZeroFunction<dim>(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 <int dim>
+void test()
+{
+  deallog << "Hdiv_seminorm:" << std::endl;
+  // sqrt(\int (div f)^2 = sqrt(\int (1+2y+1)^2)
+  test<dim>(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<dim>(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<dim>(VectorTools::H1_seminorm, std::sqrt(22.0/3.0));
+
+  deallog << "H1_norm:" << std::endl;
+  test<dim>(VectorTools::H1_norm, std::sqrt(229.0/60.0+22.0/3.0));
+
+  deallog << "L1_norm:" << std::endl;
+  test<dim>(VectorTools::L1_norm, 35.0/12.0);
+
+  deallog << "Linfty_norm:" << std::endl;
+  test<dim>(VectorTools::Linfty_norm, 3.0);
+
+  deallog << "mean:" << std::endl;
+  // int -(x+y+z + x^2+y^2 + z+xy)
+  test<dim>(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<dim>(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<dim>(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 (file)
index 0000000..d8f0ac1
--- /dev/null
@@ -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 (file)
index 0000000..630918c
--- /dev/null
@@ -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
+

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.