]> https://gitweb.dealii.org/ - dealii.git/commitdiff
VectorTools::project uses a matrix-free implemementation if supported and can deal...
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 4 Nov 2016 10:40:11 +0000 (11:40 +0100)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 17 Nov 2016 16:36:22 +0000 (17:36 +0100)
20 files changed:
doc/news/changes.h
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_project.inst.in
tests/numerics/project_parallel_common.h [new file with mode: 0644]
tests/numerics/project_parallel_q_01.cc [new file with mode: 0644]
tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=1.output [new file with mode: 0644]
tests/numerics/project_parallel_q_01.with_trilinos=true.mpirun=4.output [new file with mode: 0644]
tests/numerics/project_parallel_q_02.cc [new file with mode: 0644]
tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=1.output [new file with mode: 0644]
tests/numerics/project_parallel_q_02.with_trilinos=true.mpirun=4.output [new file with mode: 0644]
tests/numerics/project_parallel_q_03.cc [new file with mode: 0644]
tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=1.output [new file with mode: 0644]
tests/numerics/project_parallel_q_03.with_trilinos=true.mpirun=4.output [new file with mode: 0644]
tests/numerics/project_parallel_q_04.cc [new file with mode: 0644]
tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=1.output [new file with mode: 0644]
tests/numerics/project_parallel_q_04.with_trilinos=true.mpirun=4.output [new file with mode: 0644]
tests/numerics/project_parallel_q_05.cc [new file with mode: 0644]
tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=1.output [new file with mode: 0644]
tests/numerics/project_parallel_q_05.with_trilinos=true.mpirun=4.output [new file with mode: 0644]

index c2fc5e7416546d60b8567ade82dcabb2ba190a59..ab13226ffefe9b897a480689a2e921b657544187 100644 (file)
@@ -465,7 +465,7 @@ inconvenience this causes.
  <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>
@@ -486,6 +486,12 @@ inconvenience this causes.
  (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)
@@ -535,7 +541,6 @@ inconvenience this causes.
  (Mayank Sabharwal, Wolfgang Bangerth, 2016/10/25)
  </li>
 
-
  <li> New: Add ArpackSolver::set_shift() to set the shift value in spectral
  transformation.
  <br>
index 8509ed417b7c0fe86259f52b114a731c1b7d6a85..f7a099628976bf398ccac2c7ff19e183fe2cdb52 100644 (file)
@@ -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<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
index b6bfa29e043012d97b9961500efb06872742d924..23ed9d56a46b813accfc2e4880730ca8921fa5a1 100644 (file)
@@ -19,6 +19,7 @@
 
 #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,
@@ -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 <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,
@@ -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<LinearAlgebra::distributed::Vector<Number> > cg(control);
       typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
       PreconditionJacobi<MatrixType> 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<LinearAlgebra::distributed::Vector<Number> > cg(control);
       typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
       PreconditionJacobi<MatrixType> preconditioner;
@@ -1163,6 +1356,87 @@ namespace VectorTools
 
 
 
+  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,
@@ -1174,10 +1448,27 @@ namespace VectorTools
                 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);
+      }
   }
 
 
@@ -1209,6 +1500,9 @@ namespace VectorTools
                 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,
@@ -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<Vector<number> > memory;
       SolverCG<Vector<number> >            cg(control,memory);
index 27e1397b54843734fc4121212359ee58071381e2..7c045dd3f5d19ea3d5cb13701319ce86f62c87d7 100644 (file)
@@ -20,7 +20,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
     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                &,
@@ -32,7 +32,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
      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>   &,
diff --git a/tests/numerics/project_parallel_common.h b/tests/numerics/project_parallel_common.h
new file mode 100644 (file)
index 0000000..348daea
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/numerics/project_parallel_q_01.cc b/tests/numerics/project_parallel_q_01.cc
new file mode 100644 (file)
index 0000000..2c5b11c
--- /dev/null
@@ -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 <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));
+}
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 (file)
index 0000000..9037bc6
--- /dev/null
@@ -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 (file)
index 0000000..9037bc6
--- /dev/null
@@ -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 (file)
index 0000000..ad66046
--- /dev/null
@@ -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 <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));
+}
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 (file)
index 0000000..c31b059
--- /dev/null
@@ -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 (file)
index 0000000..c31b059
--- /dev/null
@@ -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 (file)
index 0000000..fb2f597
--- /dev/null
@@ -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 <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));
+}
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 (file)
index 0000000..3f7b8d7
--- /dev/null
@@ -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 (file)
index 0000000..3f7b8d7
--- /dev/null
@@ -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 (file)
index 0000000..e605985
--- /dev/null
@@ -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 <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));
+}
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 (file)
index 0000000..e0be52d
--- /dev/null
@@ -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 (file)
index 0000000..e0be52d
--- /dev/null
@@ -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 (file)
index 0000000..a3ed5dc
--- /dev/null
@@ -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 <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));
+}
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 (file)
index 0000000..3136b53
--- /dev/null
@@ -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 (file)
index 0000000..3136b53
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.