From 8459429fc05b4dfa8d0a493ad821e73e5c394eb6 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 4 Nov 2016 11:40:11 +0100 Subject: [PATCH] VectorTools::project uses a matrix-free implemementation if supported and can deal with parallel Triangulations in this case --- doc/news/changes.h | 9 +- include/deal.II/numerics/vector_tools.h | 58 ++- .../deal.II/numerics/vector_tools.templates.h | 317 +++++++++++++++- source/numerics/vector_tools_project.inst.in | 4 +- tests/numerics/project_parallel_common.h | 359 ++++++++++++++++++ tests/numerics/project_parallel_q_01.cc | 28 ++ ...el_q_01.with_trilinos=true.mpirun=1.output | 23 ++ ...el_q_01.with_trilinos=true.mpirun=4.output | 23 ++ tests/numerics/project_parallel_q_02.cc | 28 ++ ...el_q_02.with_trilinos=true.mpirun=1.output | 23 ++ ...el_q_02.with_trilinos=true.mpirun=4.output | 23 ++ tests/numerics/project_parallel_q_03.cc | 28 ++ ...el_q_03.with_trilinos=true.mpirun=1.output | 78 ++++ ...el_q_03.with_trilinos=true.mpirun=4.output | 78 ++++ tests/numerics/project_parallel_q_04.cc | 28 ++ ...el_q_04.with_trilinos=true.mpirun=1.output | 12 + ...el_q_04.with_trilinos=true.mpirun=4.output | 12 + tests/numerics/project_parallel_q_05.cc | 28 ++ ...el_q_05.with_trilinos=true.mpirun=1.output | 34 ++ ...el_q_05.with_trilinos=true.mpirun=4.output | 34 ++ 20 files changed, 1202 insertions(+), 25 deletions(-) create mode 100644 tests/numerics/project_parallel_common.h create mode 100644 tests/numerics/project_parallel_q_01.cc create mode 100644 tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=1.output create mode 100644 tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=4.output create mode 100644 tests/numerics/project_parallel_q_02.cc create mode 100644 tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=1.output create mode 100644 tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=4.output create mode 100644 tests/numerics/project_parallel_q_03.cc create mode 100644 tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=1.output create mode 100644 tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=4.output create mode 100644 tests/numerics/project_parallel_q_04.cc create mode 100644 tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=1.output create mode 100644 tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=4.output create mode 100644 tests/numerics/project_parallel_q_05.cc create mode 100644 tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=1.output create mode 100644 tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=4.output diff --git a/doc/news/changes.h b/doc/news/changes.h index c2fc5e7416..ab13226ffe 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -465,7 +465,7 @@ inconvenience this causes.
(Sebastian Gonzalez-Pintor, Wolfgang Bangerth, 2016/11/15) - +
  • Fixed: Objects of type TrilinosWrappers::SparsityPattern::const_iterator were entirely unusable due to a bug. This is now fixed.
    @@ -486,6 +486,12 @@ inconvenience this causes. (Denis Davydov, 2016/10/31)
    • +
  • New: VectorTools::project can be used for parallel Triangulations in + specific cases. +
    + (Daniel Arndt, 2016/10/30) +
    • +
  • New: Add MatrixFreeOperators::LaplaceOperator representing a Laplace matrix.
    (Denis Davydov, 2016/10/30) @@ -535,7 +541,6 @@ inconvenience this causes. (Mayank Sabharwal, Wolfgang Bangerth, 2016/10/25)
    • -
  • New: Add ArpackSolver::set_shift() to set the shift value in spectral transformation.
    diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 8509ed417b..f7a0996289 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -763,21 +763,31 @@ namespace VectorTools * 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. @@ -787,7 +797,8 @@ namespace VectorTools * 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 @@ -922,6 +933,33 @@ namespace VectorTools const std_cxx11::function< VectorizedArray (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 + void project_generic (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary = false, + const Quadrature &q_boundary = (dim > 1 ? + QGauss(2) : + Quadrature(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 diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index b6bfa29e04..23ed9d56a4 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -50,6 +51,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -77,10 +82,8 @@ DEAL_II_NAMESPACE_OPEN - namespace VectorTools { - template class DoFHandlerType> void interpolate (const Mapping &mapping, @@ -756,6 +759,7 @@ namespace VectorTools } + /** * Return whether the boundary values try to constrain a degree of freedom * that is already constrained to something else @@ -891,6 +895,187 @@ namespace VectorTools + /** + * 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 + void project_parallel_helper (const Mapping &mapping, + const dealii::DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) + { + switch (dof.get_fe().degree) + { + case 1: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 2: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 3: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 4: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 5: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 6: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 7: + dealii::VectorTools::project_generic + (mapping, dof, constraints, quadrature, function, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + break; + + case 8: + dealii::VectorTools::project_generic + (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 + void project_parallel (const Mapping &mapping, + const dealii::DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary, + const Quadrature &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 + void project_matrix_free (const Mapping &mapping, + const dealii::DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary, + const Quadrature &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 *fe_ptr = &dof.get_fe().base_element(0); + if (fe_ptr->n_components() != 1) + use_matrix_free = false; + else if ((dynamic_cast,dim,dim>*>(fe_ptr)==0) + && (dynamic_cast >,dim,dim>*>(fe_ptr)==0) + &&(dynamic_cast*>(fe_ptr)==0) + &&(dynamic_cast*>(fe_ptr)==0) + &&(dynamic_cast*>(fe_ptr)==0) + &&(dynamic_cast*>(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* > (&(dof.get_triangulation()))==0), + ExcNotImplemented()); + do_project (mapping, dof, constraints, quadrature, + function, vec_result, + enforce_zero_boundary, q_boundary, + project_to_boundary_first); + } + } + + + template void project_parallel (const Mapping &mapping, const DoFHandler &dof, @@ -970,8 +1155,12 @@ namespace VectorTools 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 > cg(control); typename PreconditionJacobi::AdditionalData data(0.8); PreconditionJacobi preconditioner; @@ -1040,8 +1229,12 @@ namespace VectorTools 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 > cg(control); typename PreconditionJacobi::AdditionalData data(0.8); PreconditionJacobi preconditioner; @@ -1163,6 +1356,87 @@ namespace VectorTools + template + void project_generic (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) + { + const parallel::Triangulation *parallel_tria = + dynamic_cast*>(&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::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::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 matrix_free; + matrix_free.reinit (mapping, dof, constraints, + QGauss<1>(fe_degree+2), additional_data); + typedef MatrixFreeOperators::MassOperator MatrixType; + MatrixType mass_matrix; + mass_matrix.initialize(matrix_free); + mass_matrix.compute_diagonal(); + + typedef LinearAlgebra::distributed::Vector 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 > cg(control); + PreconditionJacobi 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 void project (const Mapping &mapping, const DoFHandler &dof, @@ -1174,10 +1448,27 @@ namespace VectorTools const Quadrature &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 *const mapping_ptr = dynamic_cast*> (&mapping); + const DoFHandler *const dof_ptr = dynamic_cast*> (&dof); + const Function *const function_ptr + = dynamic_cast*> (&function); + Assert (mapping_ptr!=0, ExcInternalError()); + Assert (dof_ptr!=0, ExcInternalError()); + project_matrix_free + (*mapping_ptr, *dof_ptr, constraints, quadrature, *function_ptr, vec_result, + enforce_zero_boundary, q_boundary, project_to_boundary_first); + } + else + { + Assert((dynamic_cast* > (&(dof.get_triangulation()))==0), + ExcNotImplemented()); + do_project (mapping, dof, constraints, quadrature, + function, vec_result, + enforce_zero_boundary, q_boundary, + project_to_boundary_first); + } } @@ -1209,6 +1500,9 @@ namespace VectorTools const hp::QCollection &q_boundary, const bool project_to_boundary_first) { + Assert((dynamic_cast* > (&(dof.get_triangulation()))==0), + ExcNotImplemented()); + do_project (mapping, dof, constraints, quadrature, function, vec_result, enforce_zero_boundary, q_boundary, @@ -2504,7 +2798,8 @@ namespace VectorTools { // 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 > memory; SolverCG > cg(control,memory); diff --git a/source/numerics/vector_tools_project.inst.in b/source/numerics/vector_tools_project.inst.in index 27e1397b54..7c045dd3f5 100644 --- a/source/numerics/vector_tools_project.inst.in +++ b/source/numerics/vector_tools_project.inst.in @@ -20,7 +20,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens namespace VectorTools \{ template - void project + void project (const Mapping &, const DoFHandler &, const ConstraintMatrix &, @@ -32,7 +32,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const bool); template - void project + void project (const DoFHandler &, const ConstraintMatrix &, const Quadrature &, diff --git a/tests/numerics/project_parallel_common.h b/tests/numerics/project_parallel_common.h new file mode 100644 index 0000000000..348daeabfc --- /dev/null +++ b/tests/numerics/project_parallel_common.h @@ -0,0 +1,359 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + + +template +void test (); + + +template +class F : public Function +{ +public: + F (const unsigned int q, + const unsigned int n_components) + : + Function(n_components), + q(q) + {} + + virtual double value (const Point &p, + const unsigned int component) const + { + Assert ((component == 0) && (this->n_components == 1), + ExcInternalError()); + double val = 0; + for (unsigned int d=0; d &p, + Vector &v) const + { + for (unsigned int c=0; c +void do_project (const parallel::distributed::Triangulation &triangulation, + const FiniteElement &fe, + const unsigned int order_difference) +{ + const unsigned int p = fe_degree; + DoFHandler 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 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(p+2), + F (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 (q, fe.n_components()), + error, + QGauss(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 +void test_no_hanging_nodes (const FiniteElement &fe, + const unsigned int order_difference = 0) +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); + + do_project (triangulation, fe, order_difference); +} + + + +// same test as above, but this time with a mesh that has hanging nodes +template +void test_with_hanging_nodes (const FiniteElement &fe, + const unsigned int order_difference = 0) +{ + parallel::distributed::Triangulation 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 (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 +void test_with_wrong_face_orientation (const FiniteElement &fe, + const unsigned int order_difference = 0) +{ + if (dim != 3) + return; + + for (unsigned int i=0; i<7; ++i) + { + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_ball (triangulation); + typename parallel::distributed::Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + std::advance (cell, i); + cell->set_refine_flag (); + triangulation.execute_coarsening_and_refinement (); + + do_project (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 +void test_with_2d_deformed_mesh (const FiniteElement &fe, + const unsigned int order_difference = 0) +{ + if (dim != 2) + return; + + std::vector > points_glob; + std::vector > points; + + points_glob.push_back (Point (0.0, 0.0)); + points_glob.push_back (Point (1.0, 0.0)); + points_glob.push_back (Point (1.0, 0.5)); + points_glob.push_back (Point (1.0, 1.0)); + points_glob.push_back (Point (0.6, 0.5)); + points_glob.push_back (Point (0.5, 1.0)); + points_glob.push_back (Point (0.0, 1.0)); + + // Prepare cell data + std::vector > 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 triangulation(MPI_COMM_WORLD); + triangulation.create_triangulation (points_glob, cells, SubCellData()); + + do_project (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 +void test_with_2d_deformed_refined_mesh (const FiniteElement &fe, + const unsigned int order_difference = 0) +{ + if (dim != 2) + return; + + for (unsigned int i=0; i<3; ++i) + { + std::vector > points_glob; + std::vector > points; + + points_glob.push_back (Point (0.0, 0.0)); + points_glob.push_back (Point (1.0, 0.0)); + points_glob.push_back (Point (1.0, 0.5)); + points_glob.push_back (Point (1.0, 1.0)); + points_glob.push_back (Point (0.6, 0.5)); + points_glob.push_back (Point (0.5, 1.0)); + points_glob.push_back (Point (0.0, 1.0)); + + // Prepare cell data + std::vector > 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 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 (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>(); +} diff --git a/tests/numerics/project_parallel_q_01.cc b/tests/numerics/project_parallel_q_01.cc new file mode 100644 index 0000000000..2c5b11c418 --- /dev/null +++ b/tests/numerics/project_parallel_q_01.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// 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 +void test () +{ + test_no_hanging_nodes (FESystem(FE_Q(1), 2)); + test_no_hanging_nodes (FESystem(FE_Q(2), 1)); +} diff --git a/tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=1.output b/tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..9037bc6dd3 --- /dev/null +++ b/tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=1.output @@ -0,0 +1,23 @@ + +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 diff --git a/tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=4.output b/tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..9037bc6dd3 --- /dev/null +++ b/tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=4.output @@ -0,0 +1,23 @@ + +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 diff --git a/tests/numerics/project_parallel_q_02.cc b/tests/numerics/project_parallel_q_02.cc new file mode 100644 index 0000000000..ad6604690b --- /dev/null +++ b/tests/numerics/project_parallel_q_02.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// 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 +void test () +{ + test_with_hanging_nodes (FESystem(FE_Q(1), 2)); + test_with_hanging_nodes (FESystem(FE_Q(2), 1)); +} diff --git a/tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=1.output b/tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..c31b059144 --- /dev/null +++ b/tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=1.output @@ -0,0 +1,23 @@ + +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 diff --git a/tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=4.output b/tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..c31b059144 --- /dev/null +++ b/tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=4.output @@ -0,0 +1,23 @@ + +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 diff --git a/tests/numerics/project_parallel_q_03.cc b/tests/numerics/project_parallel_q_03.cc new file mode 100644 index 0000000000..fb2f597fea --- /dev/null +++ b/tests/numerics/project_parallel_q_03.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// 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 +void test () +{ + test_with_wrong_face_orientation (FESystem(FE_Q(1), 2)); + test_with_wrong_face_orientation (FESystem(FE_Q(2), 1)); +} diff --git a/tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=1.output b/tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..3f7b8d79df --- /dev/null +++ b/tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=1.output @@ -0,0 +1,78 @@ + +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 diff --git a/tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=4.output b/tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..3f7b8d79df --- /dev/null +++ b/tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=4.output @@ -0,0 +1,78 @@ + +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 diff --git a/tests/numerics/project_parallel_q_04.cc b/tests/numerics/project_parallel_q_04.cc new file mode 100644 index 0000000000..e605985ac8 --- /dev/null +++ b/tests/numerics/project_parallel_q_04.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// 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 +void test () +{ + test_with_2d_deformed_mesh (FESystem(FE_Q(1), 2)); + test_with_2d_deformed_mesh (FESystem(FE_Q(2), 1)); +} diff --git a/tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=1.output b/tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..e0be52d746 --- /dev/null +++ b/tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=1.output @@ -0,0 +1,12 @@ + +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 diff --git a/tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=4.output b/tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..e0be52d746 --- /dev/null +++ b/tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=4.output @@ -0,0 +1,12 @@ + +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 diff --git a/tests/numerics/project_parallel_q_05.cc b/tests/numerics/project_parallel_q_05.cc new file mode 100644 index 0000000000..a3ed5dc17e --- /dev/null +++ b/tests/numerics/project_parallel_q_05.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// 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 +void test () +{ + test_with_2d_deformed_refined_mesh (FESystem(FE_Q(1), 2)); + test_with_2d_deformed_refined_mesh (FESystem(FE_Q(2), 1)); +} diff --git a/tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=1.output b/tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..3136b5379b --- /dev/null +++ b/tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=1.output @@ -0,0 +1,34 @@ + +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 diff --git a/tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=4.output b/tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..3136b5379b --- /dev/null +++ b/tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=4.output @@ -0,0 +1,34 @@ + +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 -- 2.39.5