<br>
(Sebastian Gonzalez-Pintor, Wolfgang Bangerth, 2016/11/15)
</li>
-
+
<li> Fixed: Objects of type TrilinosWrappers::SparsityPattern::const_iterator
were entirely unusable due to a bug. This is now fixed.
<br>
(Denis Davydov, 2016/10/31)
</li>
+ <li> New: VectorTools::project can be used for parallel Triangulations in
+ specific cases.
+ <br>
+ (Daniel Arndt, 2016/10/30)
+ </li>
+
<li> New: Add MatrixFreeOperators::LaplaceOperator representing a Laplace matrix.
<br>
(Denis Davydov, 2016/10/30)
(Mayank Sabharwal, Wolfgang Bangerth, 2016/10/25)
</li>
-
<li> New: Add ArpackSolver::set_shift() to set the shift value in spectral
transformation.
<br>
* that you need not give a second quadrature formula if you don't want to
* project to the boundary first, but that you must if you want to do so.
*
- * This function needs the mass matrix of the finite element space on the
- * present grid. To this end, the mass matrix is assembled exactly using
- * MatrixTools::create_mass_matrix. This function performs numerical
- * quadrature using the given quadrature rule; you should therefore make
- * sure that the given quadrature formula is also sufficient for the
- * integration of the mass matrix.
+ * A MatrixFree implementation is used if the following conditions are met:
+ * - @p enforce_zero_boundary is false,
+ * - @p project_to_boundary_first is false,
+ * - the FiniteElement is supported by the MatrixFree class,
+ * - the FiniteElement has less than five components
+ * - the degree of the FiniteElement is less than nine.
+ * - dim==spacedim
+ *
+ * In this case, this function performs numerical quadrature using the given
+ * quadrature formula for integration of the provided function while a
+ * QGauss(fe_degree+2) object is used for the mass operator. You should
+ * therefore make sure that the given quadrature formula is sufficient for
+ * creating the right-hand side.
+ *
+ * Otherwise, only serial Triangulations are supported and the mass matrix
+ * is assembled exactly using MatrixTools::create_mass_matrix and the same
+ * quadrature rule as for the right-hand side.
+ * You should therefore make sure that the given quadrature formula is also
+ * sufficient for creating the mass matrix.
*
* See the general documentation of this namespace for further information.
*
* In 1d, the default value of the boundary quadrature formula is an invalid
* object since integration on the boundary doesn't happen in 1d.
*
- * @note This function is not implemented for MPI parallel computations,
- * see step-32 for a way to do a projection in parallel.
- *
* @param[in] mapping The mapping object to use.
* @param[in] dof The DoFHandler the describes the finite element space to
* project into and that corresponds to @p vec.
* mass matrix.
* @param[in] function The function to project into the finite element space.
* @param[out] vec The output vector where the projected function will be
- * stored in. This vector is required to be already initialized.
+ * stored in. This vector is required to be already initialized and must not
+ * have ghost elements.
* @param[in] enforce_zero_boundary If true, @p vec will have zero boundary
* conditions.
* @param[in] q_boundary Quadrature rule to be used if @p project_to_boundary_first
const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
VectorType &vec_result);
+ /**
+ * Implementation for the project() function with finite elements
+ * are supported by the MatrixFree class for arbitrary number of
+ * components and degree of the FiniteElement.
+ *
+ * This function should be used if you have more than four components
+ * or the degree of your FiniteElement is higher than eight. For all the
+ * other cases project() is already instantiated.
+ *
+ * The first two template arguments have to be specified explicitly.
+ * @p vec_result is expected to not have any ghost entries.
+ * @p project_to_boundary_first and @p enforce_zero_boundary are not yet
+ * implemented.
+ */
+ template <int components, int fe_degree, int dim, typename VectorType, int spacedim>
+ void project_generic (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim, spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const Function<spacedim, typename VectorType::value_type> &function,
+ VectorType &vec_result,
+ const bool enforce_zero_boundary = false,
+ const Quadrature<dim-1> &q_boundary = (dim > 1 ?
+ QGauss<dim-1>(2) :
+ Quadrature<dim-1>(0)),
+ const bool project_to_boundary_first = false);
+
/**
* Compute Dirichlet boundary conditions. This function makes up a map of
* degrees of freedom subject to Dirichlet boundary conditions and the
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/function.h>
+#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/distributed/tria_base.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_system.h>
DEAL_II_NAMESPACE_OPEN
-
namespace VectorTools
{
-
template <int dim, int spacedim, typename VectorType,
template <int, int> class DoFHandlerType>
void interpolate (const Mapping<dim,spacedim> &mapping,
}
+
/**
* Return whether the boundary values try to constrain a degree of freedom
* that is already constrained to something else
+ /**
+ * Helper interface. After figuring out the number of components
+ * in project_parallel, we determine the degree of the FiniteElement
+ * and call project_generic with the appropriate template arguments.
+ */
+ template <int components, int dim, typename VectorType, int spacedim>
+ void project_parallel_helper (const Mapping<dim, spacedim> &mapping,
+ const dealii::DoFHandler<dim, spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const Function<spacedim, typename VectorType::value_type> &function,
+ VectorType &vec_result,
+ const bool enforce_zero_boundary,
+ const Quadrature<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ switch (dof.get_fe().degree)
+ {
+ case 1:
+ dealii::VectorTools::project_generic<components, 1>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 2:
+ dealii::VectorTools::project_generic<components, 2>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 3:
+ dealii::VectorTools::project_generic<components, 3>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 4:
+ dealii::VectorTools::project_generic<components, 4>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 5:
+ dealii::VectorTools::project_generic<components, 5>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 6:
+ dealii::VectorTools::project_generic<components, 6>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 7:
+ dealii::VectorTools::project_generic<components, 7>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 8:
+ dealii::VectorTools::project_generic<components, 8>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ default:
+ Assert(false,
+ ExcMessage("Call project_generic yourself specifying the "
+ "appropriate template arguments!"));
+ }
+ }
+
+
+
+ // Helper interface for the matrix-free implementation of project().
+ // Used to determine the number of components.
+ template <int dim, typename VectorType, int spacedim>
+ void project_parallel (const Mapping<dim, spacedim> &mapping,
+ const dealii::DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const Function<spacedim, typename VectorType::value_type> &function,
+ VectorType &vec_result,
+ const bool enforce_zero_boundary,
+ const Quadrature<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ switch (dof.get_fe().n_components())
+ {
+ case 1:
+ project_parallel_helper<1>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 2:
+ project_parallel_helper<2>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 3:
+ project_parallel_helper<3>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ case 4:
+ project_parallel_helper<4>
+ (mapping, dof, constraints, quadrature, function, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ break;
+
+ default:
+ Assert(false,
+ ExcMessage("Call project_generic yourself specifying the "
+ "appropriate template arguments!"));
+ }
+ }
+
+
+
+ template <int dim, typename VectorType>
+ void project_matrix_free (const Mapping<dim> &mapping,
+ const dealii::DoFHandler<dim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const Function<dim, typename VectorType::value_type> &function,
+ VectorType &vec_result,
+ const bool enforce_zero_boundary,
+ const Quadrature<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ // If we can, use the matrix-free implementation
+ bool use_matrix_free = true;
+ // enforce_zero_boundary and project_to_boundary_first
+ // are not yet supported
+ if (enforce_zero_boundary || project_to_boundary_first)
+ use_matrix_free = false;
+ else
+ {
+ // Find out if the FiniteElement is supported
+ // This is copied from matrix_free/shape_info.templates.h
+ // and matrix_free/matrix_free.templates.h
+ if (dof.get_fe().degree==0 || dof.get_fe().n_base_elements()!=1)
+ use_matrix_free = false;
+ else
+ {
+ const FiniteElement<dim> *fe_ptr = &dof.get_fe().base_element(0);
+ if (fe_ptr->n_components() != 1)
+ use_matrix_free = false;
+ else if ((dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(fe_ptr)==0)
+ && (dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
+ Polynomials::PiecewisePolynomial<double> >,dim,dim>*>(fe_ptr)==0)
+ &&(dynamic_cast<const FE_DGP<dim>*>(fe_ptr)==0)
+ &&(dynamic_cast<const FE_Q_DG0<dim>*>(fe_ptr)==0)
+ &&(dynamic_cast<const FE_Q<dim>*>(fe_ptr)==0)
+ &&(dynamic_cast<const FE_DGQ<dim>*>(fe_ptr)==0))
+ use_matrix_free = false;
+ }
+ }
+
+ if (use_matrix_free)
+ project_parallel (mapping, dof, constraints, quadrature,
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
+ else
+ {
+ Assert((dynamic_cast<const parallel::Triangulation<dim>* > (&(dof.get_triangulation()))==0),
+ ExcNotImplemented());
+ do_project (mapping, dof, constraints, quadrature,
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
+ }
+ }
+
+
+
template <int dim, typename VectorType, int spacedim, int fe_degree>
void project_parallel (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
mass_matrix.vmult_add(rhs, inhomogeneities);
- //now invert the matrix
- ReductionControl control(rhs.size(), 0., 1e-12, false, false);
+ // now invert the matrix
+ // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
+ // steps may not be sufficient, since roundoff errors may accumulate for
+ // badly conditioned matrices. This behavior can be observed, e.g. for
+ // FE_Q_Hierarchical for degree higher than three.
+ ReductionControl control(5.*rhs.size(), 0., 1e-12, false, false);
SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
PreconditionJacobi<MatrixType> preconditioner;
mass_matrix.vmult_add(rhs, inhomogeneities);
- //now invert the matrix
- ReductionControl control(rhs.size(), 0., 1e-12, false, false);
+ // now invert the matrix
+ // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
+ // steps may not be sufficient, since roundoff errors may accumulate for
+ // badly conditioned matrices. This behavior can be observed, e.g. for
+ // FE_Q_Hierarchical for degree higher than three.
+ ReductionControl control(5.*rhs.size(), 0., 1e-12, false, false);
SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
PreconditionJacobi<MatrixType> preconditioner;
+ template <int components, int fe_degree, int dim, typename VectorType, int spacedim>
+ void project_generic (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim, spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const Function<spacedim, typename VectorType::value_type> &function,
+ VectorType &vec_result,
+ const bool enforce_zero_boundary,
+ const Quadrature<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ const parallel::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::Triangulation<dim,spacedim>*>(&dof.get_tria());
+ Assert (project_to_boundary_first == false, ExcNotImplemented());
+ Assert (enforce_zero_boundary == false, ExcNotImplemented());
+ (void) enforce_zero_boundary;
+ (void) project_to_boundary_first;
+ (void) q_boundary;
+
+ const IndexSet locally_owned_dofs = dof.locally_owned_dofs();
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs);
+
+ typedef typename VectorType::value_type number;
+ Assert (dof.get_fe().n_components() == function.n_components,
+ ExcDimensionMismatch(dof.get_fe().n_components(),
+ function.n_components));
+ Assert (vec_result.size() == dof.n_dofs(),
+ ExcDimensionMismatch (vec_result.size(), dof.n_dofs()));
+ Assert (dof.get_fe().degree == fe_degree,
+ ExcDimensionMismatch(fe_degree, dof.get_fe().degree));
+ Assert (dof.get_fe().n_components() == components,
+ ExcDimensionMismatch(components, dof.get_fe().n_components()));
+
+ // set up mass matrix and right hand side
+ typename MatrixFree<dim,number>::AdditionalData additional_data;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim,number>::AdditionalData::partition_color;
+ if (parallel_tria)
+ additional_data.mpi_communicator = parallel_tria->get_communicator();
+ additional_data.mapping_update_flags = (update_values | update_JxW_values);
+ MatrixFree<dim, number> matrix_free;
+ matrix_free.reinit (mapping, dof, constraints,
+ QGauss<1>(fe_degree+2), additional_data);
+ typedef MatrixFreeOperators::MassOperator<dim, fe_degree, fe_degree+2, components, number> MatrixType;
+ MatrixType mass_matrix;
+ mass_matrix.initialize(matrix_free);
+ mass_matrix.compute_diagonal();
+
+ typedef LinearAlgebra::distributed::Vector<number> LocalVectorType;
+ LocalVectorType vec, rhs, inhomogeneities;
+ matrix_free.initialize_dof_vector(vec);
+ matrix_free.initialize_dof_vector(rhs);
+ matrix_free.initialize_dof_vector(inhomogeneities);
+ constraints.distribute(inhomogeneities);
+ inhomogeneities*=-1.;
+
+ create_right_hand_side (mapping, dof, quadrature, function, rhs, constraints);
+ mass_matrix.vmult_add(rhs, inhomogeneities);
+
+ // now invert the matrix
+ // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
+ // steps may not be sufficient, since roundoff errors may accumulate for
+ // badly conditioned matrices. This behavior can be observed, e.g. for
+ // FE_Q_Hierarchical for degree higher than three.
+ ReductionControl control(5.*rhs.size(), 0., 1e-12, false, false);
+ SolverCG<LinearAlgebra::distributed::Vector<number> > cg(control);
+ PreconditionJacobi<MatrixType> preconditioner;
+ preconditioner.initialize(mass_matrix, 1.);
+ cg.solve (mass_matrix, vec, rhs, preconditioner);
+ vec+=inhomogeneities;
+
+ constraints.distribute (vec);
+ IndexSet::ElementIterator it = locally_owned_dofs.begin();
+ for (; it!=locally_owned_dofs.end(); ++it)
+ vec_result(*it) = vec(*it);
+ vec_result.compress(VectorOperation::insert);
+ }
+
+
+
template <int dim, typename VectorType, int spacedim>
void project (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim-1> &q_boundary,
const bool project_to_boundary_first)
{
- do_project (mapping, dof, constraints, quadrature,
- function, vec_result,
- enforce_zero_boundary, q_boundary,
- project_to_boundary_first);
+ if (dim==spacedim)
+ {
+ const Mapping<dim> *const mapping_ptr = dynamic_cast<const Mapping<dim>*> (&mapping);
+ const DoFHandler<dim> *const dof_ptr = dynamic_cast<const DoFHandler<dim>*> (&dof);
+ const Function<dim, typename VectorType::value_type> *const function_ptr
+ = dynamic_cast<const Function<dim, typename VectorType::value_type>*> (&function);
+ Assert (mapping_ptr!=0, ExcInternalError());
+ Assert (dof_ptr!=0, ExcInternalError());
+ project_matrix_free<dim, VectorType>
+ (*mapping_ptr, *dof_ptr, constraints, quadrature, *function_ptr, vec_result,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ }
+ else
+ {
+ Assert((dynamic_cast<const parallel::Triangulation<dim,spacedim>* > (&(dof.get_triangulation()))==0),
+ ExcNotImplemented());
+ do_project (mapping, dof, constraints, quadrature,
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
+ }
}
const hp::QCollection<dim-1> &q_boundary,
const bool project_to_boundary_first)
{
+ Assert((dynamic_cast<const parallel::Triangulation<dim,spacedim>* > (&(dof.get_triangulation()))==0),
+ ExcNotImplemented());
+
do_project (mapping, dof, constraints, quadrature,
function, vec_result,
enforce_zero_boundary, q_boundary,
{
// Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
// steps may not be sufficient, since roundoff errors may accumulate for
- // badly conditioned matrices
+ // badly conditioned matrices. This behavior can be observed, e.g. for
+ // FE_Q_Hierarchical for degree higher than three.
ReductionControl control(5.*rhs.size(), 0., 1.e-12, false, false);
GrowingVectorMemory<Vector<number> > memory;
SolverCG<Vector<number> > cg(control,memory);
namespace VectorTools \{
template
- void project
+ void project<deal_II_dimension, VEC, deal_II_space_dimension>
(const Mapping<deal_II_dimension,deal_II_space_dimension> &,
const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
const ConstraintMatrix &,
const bool);
template
- void project
+ void project<deal_II_dimension, VEC, deal_II_space_dimension>
(const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
const ConstraintMatrix &,
const Quadrature<deal_II_dimension> &,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Common framework to check the projection error for functions of order q
+// using an element of polynomial order p. Thereby, check if a function
+// is represented exactyly if this is expected.
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/fe/fe_abf.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgp_monomial.h>
+#include <deal.II/fe/fe_dgp_nonparametric.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_hierarchical.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <fstream>
+#include <vector>
+
+
+template <int dim>
+void test ();
+
+
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ F (const unsigned int q,
+ const unsigned int n_components)
+ :
+ Function<dim>(n_components),
+ q(q)
+ {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert ((component == 0) && (this->n_components == 1),
+ ExcInternalError());
+ double val = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<=q; ++i)
+ val += (d+1)*(i+1)*std::pow (p[d], 1.*i);
+ return val;
+ }
+
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ for (unsigned int c=0; c<v.size(); ++c)
+ {
+ v(c) = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<=q; ++i)
+ v(c) += (d+1)*(i+1)*std::pow (p[d], 1.*i)+c;
+ }
+ }
+
+private:
+ const unsigned int q;
+};
+
+
+template <int dim, int components, int fe_degree>
+void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const unsigned int order_difference)
+{
+ const unsigned int p = fe_degree;
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs (fe);
+
+ deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl;
+
+ const MPI_Comm &mpi_communicator = triangulation.get_communicator();
+ const IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+ ConstraintMatrix constraints;
+ constraints.reinit(locally_relevant_dofs);
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ constraints.close ();
+
+ TrilinosWrappers::MPI::Vector projection_tmp (locally_owned_dofs,
+ mpi_communicator);
+ TrilinosWrappers::MPI::Vector projection (locally_owned_dofs,
+ locally_relevant_dofs,
+ mpi_communicator);
+ Vector<float> error (triangulation.n_active_cells());
+ for (unsigned int q=0; q<=p+2-order_difference; ++q)
+ {
+ // project the function
+ VectorTools::project
+ (dof_handler, constraints, QGauss<dim>(p+2),
+ F<dim> (q, fe.n_components()), projection_tmp);
+ // just to make sure it doesn't get
+ // forgotten: handle hanging node
+ // constraints
+ constraints.distribute (projection_tmp);
+ projection = projection_tmp;
+
+ // then compute the interpolation error
+ VectorTools::integrate_difference (dof_handler,
+ projection,
+ F<dim> (q, fe.n_components()),
+ error,
+ QGauss<dim>(std::max(p,q)+1),
+ VectorTools::L2_norm);
+ const double L2_error_local = error.l2_norm();
+ const double L2_error
+ = std::sqrt(Utilities::MPI::sum(L2_error_local * L2_error_local,
+ mpi_communicator));
+ const double projection_l2_norm = projection_tmp.l2_norm();
+ deallog << fe.get_name() << ", P_" << q
+ << ", rel. error=" << L2_error / projection_l2_norm
+ << std::endl;
+
+ if (q<=p-order_difference)
+ if (L2_error > 1e-10*projection_l2_norm)
+ deallog << "Projection failed with relative error "
+ << L2_error / projection_l2_norm
+ << std::endl;
+ }
+}
+
+
+
+// check the given element of polynomial order p. the last parameter, if
+// given, denotes a gap in convergence order; for example, the Nedelec element
+// of polynomial degree p has normal components of degree p-1 and therefore
+// can only represent polynomials of degree p-1 exactly. the gap is then 1.
+template <int dim, int components, int fe_degree>
+void test_no_hanging_nodes (const FiniteElement<dim> &fe,
+ const unsigned int order_difference = 0)
+{
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (3);
+
+ do_project<dim, components, fe_degree> (triangulation, fe, order_difference);
+}
+
+
+
+// same test as above, but this time with a mesh that has hanging nodes
+template <int dim, int components, int fe_degree>
+void test_with_hanging_nodes (const FiniteElement<dim> &fe,
+ const unsigned int order_difference = 0)
+{
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (1);
+ triangulation.begin_active()->set_refine_flag ();
+ triangulation.execute_coarsening_and_refinement ();
+ triangulation.refine_global (1);
+
+ do_project<dim, components, fe_degree> (triangulation, fe, order_difference);
+}
+
+
+
+// test with a 3d grid that has cells with face_orientation==false and hanging
+// nodes. this trips up all sorts of pieces of code, for example there was a
+// crash when computing hanging node constraints on such faces (see
+// bits/face_orientation_crash), and it triggers all sorts of other
+// assumptions that may be hidden in places
+//
+// the mesh we use is the 7 cells of the hyperball mesh in 3d, with each of
+// the cells refined in turn. that then makes 7 meshes with 14 active cells
+// each. this also cycles through all possibilities of coarser or finer cell
+// having face_orientation==false
+template <int dim, int components, int fe_degree>
+void test_with_wrong_face_orientation (const FiniteElement<dim> &fe,
+ const unsigned int order_difference = 0)
+{
+ if (dim != 3)
+ return;
+
+ for (unsigned int i=0; i<7; ++i)
+ {
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ GridGenerator::hyper_ball (triangulation);
+ typename parallel::distributed::Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ std::advance (cell, i);
+ cell->set_refine_flag ();
+ triangulation.execute_coarsening_and_refinement ();
+
+ do_project<dim, components, fe_degree> (triangulation, fe, order_difference);
+ }
+}
+
+
+
+
+// test with a 2d mesh that forms a square but subdivides it into 3
+// elements. this tests the case of the sign_change thingy in
+// fe_poly_tensor.cc
+template <int dim, int components, int fe_degree>
+void test_with_2d_deformed_mesh (const FiniteElement<dim> &fe,
+ const unsigned int order_difference = 0)
+{
+ if (dim != 2)
+ return;
+
+ std::vector<Point<dim> > points_glob;
+ std::vector<Point<dim> > points;
+
+ points_glob.push_back (Point<dim> (0.0, 0.0));
+ points_glob.push_back (Point<dim> (1.0, 0.0));
+ points_glob.push_back (Point<dim> (1.0, 0.5));
+ points_glob.push_back (Point<dim> (1.0, 1.0));
+ points_glob.push_back (Point<dim> (0.6, 0.5));
+ points_glob.push_back (Point<dim> (0.5, 1.0));
+ points_glob.push_back (Point<dim> (0.0, 1.0));
+
+ // Prepare cell data
+ std::vector<CellData<dim> > cells (3);
+
+ cells[0].vertices[0] = 0;
+ cells[0].vertices[1] = 1;
+ cells[0].vertices[2] = 4;
+ cells[0].vertices[3] = 2;
+ cells[0].material_id = 0;
+
+ cells[1].vertices[0] = 4;
+ cells[1].vertices[1] = 2;
+ cells[1].vertices[2] = 5;
+ cells[1].vertices[3] = 3;
+ cells[1].material_id = 0;
+
+ cells[2].vertices[0] = 0;
+ cells[2].vertices[1] = 4;
+ cells[2].vertices[2] = 6;
+ cells[2].vertices[3] = 5;
+ cells[2].material_id = 0;
+
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ triangulation.create_triangulation (points_glob, cells, SubCellData());
+
+ do_project<dim, components, fe_degree> (triangulation, fe, order_difference);
+}
+
+
+
+// same as test_with_2d_deformed_mesh, but refine each element in turn. this
+// makes sure we also check the sign_change thingy for refined cells
+template <int dim, int components, int fe_degree>
+void test_with_2d_deformed_refined_mesh (const FiniteElement<dim> &fe,
+ const unsigned int order_difference = 0)
+{
+ if (dim != 2)
+ return;
+
+ for (unsigned int i=0; i<3; ++i)
+ {
+ std::vector<Point<dim> > points_glob;
+ std::vector<Point<dim> > points;
+
+ points_glob.push_back (Point<dim> (0.0, 0.0));
+ points_glob.push_back (Point<dim> (1.0, 0.0));
+ points_glob.push_back (Point<dim> (1.0, 0.5));
+ points_glob.push_back (Point<dim> (1.0, 1.0));
+ points_glob.push_back (Point<dim> (0.6, 0.5));
+ points_glob.push_back (Point<dim> (0.5, 1.0));
+ points_glob.push_back (Point<dim> (0.0, 1.0));
+
+ // Prepare cell data
+ std::vector<CellData<dim> > cells (3);
+
+ cells[0].vertices[0] = 0;
+ cells[0].vertices[1] = 1;
+ cells[0].vertices[2] = 4;
+ cells[0].vertices[3] = 2;
+ cells[0].material_id = 0;
+
+ cells[1].vertices[0] = 4;
+ cells[1].vertices[1] = 2;
+ cells[1].vertices[2] = 5;
+ cells[1].vertices[3] = 3;
+ cells[1].material_id = 0;
+
+ cells[2].vertices[0] = 0;
+ cells[2].vertices[1] = 4;
+ cells[2].vertices[2] = 6;
+ cells[2].vertices[3] = 5;
+ cells[2].material_id = 0;
+
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ triangulation.create_triangulation (points_glob, cells, SubCellData());
+
+ switch (i)
+ {
+ case 0:
+ triangulation.begin_active()->set_refine_flag();
+ break;
+ case 1:
+ (++(triangulation.begin_active()))->set_refine_flag();
+ break;
+ case 2:
+ (++(++(triangulation.begin_active())))->set_refine_flag();
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ triangulation.execute_coarsening_and_refinement ();
+
+ do_project<dim, components, fe_degree> (triangulation, fe, order_difference);
+ }
+}
+
+
+
+int main (int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_init_finalize(argc, argv, 1);
+ mpi_initlog();
+ test<2>();
+ test<3>();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that VectorTools::project works for Q elements correctly
+
+#include "project_parallel_common.h"
+
+
+template <int dim>
+void test ()
+{
+ test_no_hanging_nodes<dim, 2, 1> (FESystem<dim>(FE_Q<dim>(1), 2));
+ test_no_hanging_nodes<dim, 1, 2> (FESystem<dim>(FE_Q<dim>(2), 1));
+}
--- /dev/null
+
+DEAL::n_dofs=162
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=1.81989e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=2.48720e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=6.65062e-05
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.000153713
+DEAL::n_dofs=289
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=2.31928e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=2.70614e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=3.56528e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=2.41461e-06
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=6.79248e-06
+DEAL::n_dofs=1458
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=7.35909e-18
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.08733e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=2.01238e-05
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=4.66822e-05
+DEAL::n_dofs=4913
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=5.92284e-18
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=5.33425e-18
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.45412e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=5.05504e-07
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=1.43736e-06
--- /dev/null
+
+DEAL::n_dofs=162
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=1.81989e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=2.48720e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=6.65062e-05
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.000153713
+DEAL::n_dofs=289
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=2.31928e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=2.70614e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=3.56528e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=2.41461e-06
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=6.79248e-06
+DEAL::n_dofs=1458
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=7.35909e-18
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.08733e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=2.01238e-05
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=4.66822e-05
+DEAL::n_dofs=4913
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=5.92284e-18
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=5.33425e-18
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.45412e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=5.05504e-07
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=1.43736e-06
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that VectorTools::project works for Q elements correctly
+
+#include "project_parallel_common.h"
+
+
+template <int dim>
+void test ()
+{
+ test_with_hanging_nodes<dim, 2, 1> (FESystem<dim>(FE_Q<dim>(1), 2));
+ test_with_hanging_nodes<dim, 1, 2> (FESystem<dim>(FE_Q<dim>(2), 1));
+}
--- /dev/null
+
+DEAL::n_dofs=82
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.06109e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=5.39290e-14
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.000373967
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.000951721
+DEAL::n_dofs=141
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=1.74177e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=3.86131e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.22920e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=2.61058e-05
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=7.91713e-05
+DEAL::n_dofs=446
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=2.14261e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=2.16569e-13
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.000158305
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.000390601
+DEAL::n_dofs=1375
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=1.65385e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=2.30619e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=4.86253e-14
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=8.13208e-06
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=2.39265e-05
--- /dev/null
+
+DEAL::n_dofs=82
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.06109e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=5.39290e-14
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.000373967
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.000951721
+DEAL::n_dofs=141
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=1.74177e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=3.86131e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.22920e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=2.61058e-05
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=7.91713e-05
+DEAL::n_dofs=446
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=2.14261e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=2.16569e-13
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.000158305
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.000390601
+DEAL::n_dofs=1375
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=1.65385e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=2.30619e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=4.86253e-14
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=8.13208e-06
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=2.39265e-05
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that VectorTools::project works for Q elements correctly
+
+#include "project_parallel_common.h"
+
+
+template <int dim>
+void test ()
+{
+ test_with_wrong_face_orientation<dim, 2, 1> (FESystem<dim>(FE_Q<dim>(1), 2));
+ test_with_wrong_face_orientation<dim, 1, 2> (FESystem<dim>(FE_Q<dim>(2), 1));
+}
--- /dev/null
+
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=8.46380e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.18133e-16
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0117673
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.0111866
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=5.48256e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=2.07334e-13
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0121021
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.0116930
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.39446e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.36014e-16
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.00926848
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00926122
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=7.42043e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=7.71262e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0100923
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00943746
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.73171e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=3.49454e-14
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.00984841
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00955841
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.03815e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.12255e-16
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0114470
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.0108987
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.84970e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=3.30955e-14
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0101389
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00977176
+DEAL::n_dofs=195
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=3.63837e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=2.80245e-14
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=2.68796e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00312284
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00305022
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=3.82484e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=2.22128e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.71631e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00316296
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00294735
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=2.98728e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.07050e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.14817e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00234293
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00227309
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=4.01977e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.46250e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.35756e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00210960
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00203856
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=2.84169e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.27722e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.18704e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00267340
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00253077
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=4.46008e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.80346e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.40395e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00288112
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00269959
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=3.12696e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.40254e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.10869e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00220695
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00214412
--- /dev/null
+
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=8.46380e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.18133e-16
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0117673
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.0111866
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=5.48256e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=2.07334e-13
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0121021
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.0116930
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.39446e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.36014e-16
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.00926848
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00926122
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=7.42043e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=7.71262e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0100923
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00943746
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.73171e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=3.49454e-14
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.00984841
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00955841
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.03815e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=1.12255e-16
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0114470
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.0108987
+DEAL::n_dofs=70
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_0, rel. error=6.84970e-17
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_1, rel. error=3.30955e-14
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_2, rel. error=0.0101389
+DEAL::FESystem<3>[FE_Q<3>(1)^2], P_3, rel. error=0.00977176
+DEAL::n_dofs=195
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=3.63837e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=2.80245e-14
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=2.68796e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00312284
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00305022
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=3.82484e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=2.22128e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.71631e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00316296
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00294735
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=2.98728e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.07050e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.14817e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00234293
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00227309
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=4.01977e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.46250e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.35756e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00210960
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00203856
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=2.84169e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.27722e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.18704e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00267340
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00253077
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=4.46008e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.80346e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.40395e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00288112
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00269959
+DEAL::n_dofs=194
+DEAL::FESystem<3>[FE_Q<3>(2)], P_0, rel. error=3.12696e-17
+DEAL::FESystem<3>[FE_Q<3>(2)], P_1, rel. error=1.40254e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_2, rel. error=1.10869e-13
+DEAL::FESystem<3>[FE_Q<3>(2)], P_3, rel. error=0.00220695
+DEAL::FESystem<3>[FE_Q<3>(2)], P_4, rel. error=0.00214412
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that VectorTools::project works for Q elements correctly
+
+#include "project_parallel_common.h"
+
+
+template <int dim>
+void test ()
+{
+ test_with_2d_deformed_mesh<dim, 2, 1> (FESystem<dim>(FE_Q<dim>(1), 2));
+ test_with_2d_deformed_mesh<dim, 1, 2> (FESystem<dim>(FE_Q<dim>(2), 1));
+}
--- /dev/null
+
+DEAL::n_dofs=14
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.44571e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=3.20140e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00667964
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.0144138
+DEAL::n_dofs=19
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=8.95863e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=1.73057e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.35999e-16
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000954360
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00254228
--- /dev/null
+
+DEAL::n_dofs=14
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.44571e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=3.20140e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00667964
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.0144138
+DEAL::n_dofs=19
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=8.95863e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=1.73057e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.35999e-16
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000954360
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00254228
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that VectorTools::project works for Q elements correctly
+
+#include "project_parallel_common.h"
+
+
+template <int dim>
+void test ()
+{
+ test_with_2d_deformed_refined_mesh <dim, 2, 1> (FESystem<dim>(FE_Q<dim>(1), 2));
+ test_with_2d_deformed_refined_mesh <dim, 1, 2> (FESystem<dim>(FE_Q<dim>(2), 1));
+}
--- /dev/null
+
+DEAL::n_dofs=24
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.89749e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=1.33525e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00544145
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.0123174
+DEAL::n_dofs=24
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=8.00269e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=1.25939e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00479825
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.0101146
+DEAL::n_dofs=24
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.80857e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=1.56761e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00315064
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.00683382
+DEAL::n_dofs=37
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=3.61277e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=1.76493e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.60801e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000777566
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00214199
+DEAL::n_dofs=37
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=3.64108e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=8.55604e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=5.15666e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000621354
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00166958
+DEAL::n_dofs=37
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=8.21833e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=4.54279e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.11092e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000395444
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00105516
--- /dev/null
+
+DEAL::n_dofs=24
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.89749e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=1.33525e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00544145
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.0123174
+DEAL::n_dofs=24
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=8.00269e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=1.25939e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00479825
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.0101146
+DEAL::n_dofs=24
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_0, rel. error=5.80857e-17
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_1, rel. error=1.56761e-15
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_2, rel. error=0.00315064
+DEAL::FESystem<2>[FE_Q<2>(1)^2], P_3, rel. error=0.00683382
+DEAL::n_dofs=37
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=3.61277e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=1.76493e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.60801e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000777566
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00214199
+DEAL::n_dofs=37
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=3.64108e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=8.55604e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=5.15666e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000621354
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00166958
+DEAL::n_dofs=37
+DEAL::FESystem<2>[FE_Q<2>(2)], P_0, rel. error=8.21833e-17
+DEAL::FESystem<2>[FE_Q<2>(2)], P_1, rel. error=4.54279e-14
+DEAL::FESystem<2>[FE_Q<2>(2)], P_2, rel. error=1.11092e-13
+DEAL::FESystem<2>[FE_Q<2>(2)], P_3, rel. error=0.000395444
+DEAL::FESystem<2>[FE_Q<2>(2)], P_4, rel. error=0.00105516
projects / dealii.git / commitdiff