<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.
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());
// 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
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
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
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;
}
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
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);
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,
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());
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
//
// 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,
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
*
* 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
{
/**
* 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
//@{
/**
- * 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
* 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
* 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.
* 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
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
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
\}
#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
+ }
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
// 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());
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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
--- /dev/null
+
+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
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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
--- /dev/null
+
+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
+
projects / dealii.git / commitdiff