]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add VectorTools::project() for scalar-valued quadrature data 3290/head
authorDenis Davydov <davydden@gmail.com>
Sat, 29 Oct 2016 14:40:31 +0000 (16:40 +0200)
committerDenis Davydov <davydden@gmail.com>
Sun, 30 Oct 2016 09:29:54 +0000 (10:29 +0100)
15 files changed:
cmake/config/template-arguments.in
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_qp_common.h [new file with mode: 0644]
tests/numerics/project_parallel_qp_q_01.cc [new file with mode: 0644]
tests/numerics/project_parallel_qp_q_01.mpirun=3.with_cxx11=on.with_trilinos=true.output [new file with mode: 0644]
tests/numerics/project_parallel_qp_q_02.cc [new file with mode: 0644]
tests/numerics/project_parallel_qp_q_02.mpirun=3.with_cxx11=on.with_trilinos=true.output [new file with mode: 0644]
tests/numerics/project_parallel_qpmf_common.h [new file with mode: 0644]
tests/numerics/project_parallel_qpmf_q_01.cc [new file with mode: 0644]
tests/numerics/project_parallel_qpmf_q_01.mpirun=3.with_cxx11=on.output [new file with mode: 0644]
tests/numerics/project_parallel_qpmf_q_02.cc [new file with mode: 0644]
tests/numerics/project_parallel_qpmf_q_02.mpirun=3.with_cxx11=on.output [new file with mode: 0644]

index f2f37a3ecdb125b329e50636e1424c97393fa648..7830833e121e8444d6b7ed545bcc49c2d49bd225 100644 (file)
@@ -69,6 +69,21 @@ REAL_SERIAL_VECTORS := { Vector<double>;
                     @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR_REAL@;
                   }
 
+REAL_NONBLOCK_VECTORS := { Vector<double>;
+                    Vector<float> ;
+
+                    LinearAlgebra::Vector<double>;
+                    LinearAlgebra::Vector<float> ;
+
+                    LinearAlgebra::distributed::Vector<double>;
+                    LinearAlgebra::distributed::Vector<float> ;
+
+                    @DEAL_II_EXPAND_TRILINOS_VECTOR@;
+                    @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
+                    @DEAL_II_EXPAND_PETSC_VECTOR_REAL@;
+                    @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@;
+                  }
+
 EXTERNAL_SEQUENTIAL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_VECTOR@;
                                  @DEAL_II_EXPAND_TRILINOS_BLOCKVECTOR@;
                                  @DEAL_II_EXPAND_PETSC_VECTOR@;
index 7c697cdc63ff95b92e277632d3f025a6fdb84cbe..bfe83b245a03be78fb5490ad98373fa5659855d9 100644 (file)
@@ -405,6 +405,11 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> New: Add VectorTools::project() to do L2 projection
+ of scalar-valued quadrature point data in parallel.
+ <br>
+ (Denis Davydov, 2016/10/28)
+ <li>
 
  <li> Fixed: Increased precision of timesteps in DataOutInterface::write_pvd_record().
  <br>
index 4b24acebcafa173c2101cd55163a88a051612205..7071137f61ee188c84a105243dd14da25709aaa9 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/point.h>
+#include <deal.II/base/std_cxx11/function.h>
 #include <deal.II/dofs/function_map.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/hp/dof_handler.h>
@@ -35,6 +36,7 @@ DEAL_II_NAMESPACE_OPEN
 template <int dim, typename Number> class Function;
 template <int dim> class Quadrature;
 template <int dim> class QGauss;
+template <int dim, typename number> class MatrixFree;
 
 template <typename number> class Vector;
 template <typename number> class FullMatrix;
@@ -856,6 +858,71 @@ namespace VectorTools
                     Quadrature<dim-1>(0)),
                 const bool                          project_to_boundary_first = false);
 
+  /**
+   * The same as above for projection of scalar-valued quadrature data.
+   * The user provided function should return a value at the quadrature point
+   * based on the cell iterator and quadrature number and of course should be
+   * consistent with the provided @p quadrature object, which will be used
+   * to assemble the right-hand-side.
+   *
+   * This function can be used with lambdas:
+   * @code
+   * VectorTools::project
+   * (mapping,
+   *  dof_handler,
+   *  constraints,
+   *  quadrature_formula,
+   *  [=] (const typename DoFHandler<dim>::active_cell_iterator & cell, const unsigned int q) -> double
+   *  { return qp_data.get_data(cell)[q]->density; },
+   *  field);
+   * @endcode
+   * where <code>qp_data</code> is a CellDataStorage object, which stores
+   * quadrature point data.
+   */
+  template <int dim, typename VectorType, int spacedim>
+  void project (const Mapping<dim, spacedim>   &mapping,
+                const DoFHandler<dim,spacedim> &dof,
+                const ConstraintMatrix         &constraints,
+                const Quadrature<dim>          &quadrature,
+                const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::active_cell_iterator &, const unsigned int)> func,
+                VectorType                     &vec_result);
+
+  /**
+   * The same as above for projection of scalar-valued MatrixFree quadrature
+   * data.
+   * The user provided function @p func should return a VectorizedArray value
+   * at the quadrature point based on the cell number and quadrature number and
+   * should be consistent with the @p n_q_points_1d.
+   *
+   * This function can be used with lambdas:
+   * @code
+   * VectorTools::project
+   * (matrix_free_data,
+   *  constraints,
+   *  3,
+   *  [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray<double>
+   *  { return qp_data(cell,q); },
+   *  field);
+   * @endcode
+   * where <code>qp_data</code> is a an object of type Table<2, VectorizedArray<double> >,
+   * which stores quadrature point data.
+   */
+  template <int dim, typename VectorType>
+  void project (const MatrixFree<dim,typename VectorType::value_type> &data,
+                const ConstraintMatrix &constraints,
+                const unsigned int n_q_points_1d,
+                const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+                VectorType &vec_result);
+
+  /**
+   * Same as above but for <code>n_q_points_1d = matrix_free.get_dof_handler().get_fe().degree+1</code>.
+   */
+  template <int dim, typename VectorType>
+  void project (const MatrixFree<dim,typename VectorType::value_type> &data,
+                const ConstraintMatrix &constraints,
+                const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+                VectorType &vec_result);
+
   /**
    * Compute Dirichlet boundary conditions.  This function makes up a map of
    * degrees of freedom subject to Dirichlet boundary conditions and the
index eae1a44941277021a0866b62cbce11deb6911d6c..0040f2cacdd0f5a6aeabd374946b397e732aea8d 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/base/function.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/qprojector.h>
+#include <deal.II/distributed/tria_base.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/lac/la_vector.h>
@@ -37,6 +38,9 @@
 #include <deal.II/lac/vector_memory.h>
 #include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_boundary.h>
@@ -881,9 +885,282 @@ namespace VectorTools
       for (unsigned int i=0; i<vec.size(); ++i)
         vec_result(i) = vec(i);
     }
+
+
+
+    template <int dim, typename VectorType, int spacedim, int fe_degree>
+    void project_parallel (const Mapping<dim, spacedim>   &mapping,
+                           const DoFHandler<dim,spacedim> &dof,
+                           const ConstraintMatrix         &constraints,
+                           const Quadrature<dim>          &quadrature,
+                           const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::active_cell_iterator &, const unsigned int)> func,
+                           VectorType                     &vec_result)
+    {
+      const parallel::Triangulation<dim,spacedim> *parallel_tria =
+        dynamic_cast<const parallel::Triangulation<dim,spacedim>*>(&dof.get_tria());
+
+      typedef typename VectorType::value_type Number;
+      Assert (dof.get_fe().n_components() == 1,
+              ExcDimensionMismatch(dof.get_fe().n_components(),
+                                   1));
+      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));
+
+      // 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 !=0)
+        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+1), additional_data);
+      typedef MatrixFreeOperators::MassOperator<dim, fe_degree, 1, 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.;
+
+      // assemble right hand side:
+      {
+        FEValues<dim> fe_values (mapping, dof.get_fe(), quadrature,
+                                 update_values | update_JxW_values);
+
+        const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
+        const unsigned int   n_q_points    = quadrature.size();
+        Vector<Number>       cell_rhs (dofs_per_cell);
+        std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+        typename DoFHandler<dim, spacedim>::active_cell_iterator
+        cell = dof.begin_active(),
+        endc = dof.end();
+        for (; cell!=endc; ++cell)
+          if (cell->is_locally_owned())
+            {
+              cell_rhs = 0;
+              fe_values.reinit (cell);
+              for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                {
+                  const double val_q = func(cell, q_point);
+                  for (unsigned int i=0; i<dofs_per_cell; ++i)
+                    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                    val_q *
+                                    fe_values.JxW(q_point));
+                }
+
+              cell->get_dof_indices (local_dof_indices);
+              constraints.distribute_local_to_global (cell_rhs,
+                                                      local_dof_indices,
+                                                      rhs);
+            }
+        rhs.compress (VectorOperation::add);
+      }
+
+      mass_matrix.vmult_add(rhs, inhomogeneities);
+
+      //now invert the matrix
+      ReductionControl     control(rhs.size(), 0., 1e-12, false, false);
+      SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
+      typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
+      data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
+      PreconditionChebyshev<MatrixType, LocalVectorType> preconditioner;
+      preconditioner.initialize(mass_matrix, data);
+      cg.solve (mass_matrix, vec, rhs, preconditioner);
+      vec+=inhomogeneities;
+
+      constraints.distribute (vec);
+
+      const IndexSet locally_owned_dofs = dof.locally_owned_dofs();
+      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 fe_degree, int n_q_points_1d>
+    void project_parallel (const MatrixFree<dim,typename VectorType::value_type> &matrix_free,
+                           const ConstraintMatrix &constraints,
+                           const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+                           VectorType &vec_result)
+    {
+      const DoFHandler<dim> &dof = matrix_free.get_dof_handler();
+
+      typedef typename VectorType::value_type Number;
+      Assert (dof.get_fe().n_components() == 1,
+              ExcDimensionMismatch(dof.get_fe().n_components(),
+                                   1));
+      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));
+
+      typedef MatrixFreeOperators::MassOperator<dim, fe_degree, 1, 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.;
+
+      // assemble right hand side:
+      {
+        FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval(matrix_free);
+        const unsigned int n_cells = matrix_free.n_macro_cells();
+        const unsigned int n_q_points = fe_eval.n_q_points;
+
+        for (unsigned int cell=0; cell<n_cells; ++cell)
+          {
+            fe_eval.reinit(cell);
+            for (unsigned int q= 0; q < n_q_points; ++q)
+              fe_eval.submit_value(func(cell,q), q);
+
+            fe_eval.integrate (true,false);
+            fe_eval.distribute_local_to_global (rhs);
+          }
+        rhs.compress(VectorOperation::add);
+      }
+
+      mass_matrix.vmult_add(rhs, inhomogeneities);
+
+      //now invert the matrix
+      ReductionControl     control(rhs.size(), 0., 1e-12, false, false);
+      SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
+      typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
+      data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
+      PreconditionChebyshev<MatrixType, LocalVectorType> preconditioner;
+      preconditioner.initialize(mass_matrix, data);
+      cg.solve (mass_matrix, vec, rhs, preconditioner);
+      vec+=inhomogeneities;
+
+      constraints.distribute (vec);
+
+      const IndexSet locally_owned_dofs = dof.locally_owned_dofs();
+      IndexSet::ElementIterator it = locally_owned_dofs.begin();
+      for (; it!=locally_owned_dofs.end(); ++it)
+        vec_result(*it) = vec(*it);
+      vec_result.compress(VectorOperation::insert);
+    }
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  void project (const Mapping<dim, spacedim>   &mapping,
+                const DoFHandler<dim,spacedim> &dof,
+                const ConstraintMatrix         &constraints,
+                const Quadrature<dim>          &quadrature,
+                const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::active_cell_iterator &, const unsigned int)> func,
+                VectorType                     &vec_result)
+  {
+    switch (dof.get_fe().degree)
+      {
+      case 1:
+        project_parallel<dim, VectorType,spacedim,1> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 2:
+        project_parallel<dim, VectorType,spacedim,2> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 3:
+        project_parallel<dim, VectorType,spacedim,3> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 4:
+        project_parallel<dim, VectorType,spacedim,4> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 5:
+        project_parallel<dim, VectorType,spacedim,5> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 6:
+        project_parallel<dim, VectorType,spacedim,6> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 7:
+        project_parallel<dim, VectorType,spacedim,7> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      case 8:
+        project_parallel<dim, VectorType,spacedim,8> (mapping,dof,constraints,quadrature,func,vec_result);
+        break;
+      default:
+        Assert (false, ExcNotImplemented());
+      }
   }
 
 
+
+  template <int dim, typename VectorType>
+  void project (const MatrixFree<dim,typename VectorType::value_type> &matrix_free,
+                const ConstraintMatrix &constraints,
+                const unsigned int n_q_points_1d,
+                const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+                VectorType &vec_result)
+  {
+    const unsigned int fe_degree = matrix_free.get_dof_handler().get_fe().degree;
+
+    Assert (fe_degree+1 == n_q_points_1d,
+            ExcNotImplemented());
+
+    switch (fe_degree)
+      {
+      case 1:
+        project_parallel<dim, VectorType,1,2> (matrix_free,constraints,func,vec_result);
+        break;
+      case 2:
+        project_parallel<dim, VectorType,2,3> (matrix_free,constraints,func,vec_result);
+        break;
+      case 3:
+        project_parallel<dim, VectorType,3,4> (matrix_free,constraints,func,vec_result);
+        break;
+      case 4:
+        project_parallel<dim, VectorType,4,5> (matrix_free,constraints,func,vec_result);
+        break;
+      case 5:
+        project_parallel<dim, VectorType,5,6> (matrix_free,constraints,func,vec_result);
+        break;
+      case 6:
+        project_parallel<dim, VectorType,6,7> (matrix_free,constraints,func,vec_result);
+        break;
+      case 7:
+        project_parallel<dim, VectorType,7,8> (matrix_free,constraints,func,vec_result);
+        break;
+      case 8:
+        project_parallel<dim, VectorType,8,9> (matrix_free,constraints,func,vec_result);
+        break;
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+  }
+
+
+
+  template <int dim, typename VectorType>
+  void project (const MatrixFree<dim,typename VectorType::value_type> &matrix_free,
+                const ConstraintMatrix &constraints,
+                const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+                VectorType &vec_result)
+  {
+    project (matrix_free,
+             constraints,
+             matrix_free.get_dof_handler().get_fe().degree+1,
+             func,
+             vec_result);
+  }
+
+
+
   template <int dim, typename VectorType, int spacedim>
   void project (const Mapping<dim, spacedim>   &mapping,
                 const DoFHandler<dim,spacedim> &dof,
index 27e1397b54843734fc4121212359ee58071381e2..1c2d2f414129d67362a826ee3e2783bd9f2f5f8f 100644 (file)
@@ -67,3 +67,36 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
     \}
 #endif
 }
+
+for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS)
+{
+    namespace VectorTools \{
+
+    template
+    void project<deal_II_dimension,VEC,deal_II_dimension>(
+        const Mapping<deal_II_dimension, deal_II_dimension>   &,
+        const DoFHandler<deal_II_dimension,deal_II_dimension> &,
+        const ConstraintMatrix &,
+        const Quadrature<deal_II_dimension> &,
+        const std_cxx11::function<VEC::value_type (const DoFHandler<deal_II_dimension, deal_II_dimension>::active_cell_iterator &, const unsigned int)>,
+        VEC &);
+
+    template
+    void project<deal_II_dimension, VEC>(
+        const MatrixFree<deal_II_dimension,VEC::value_type> &matrix_free,
+        const ConstraintMatrix &constraints,
+        const std_cxx11::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)>,
+        VEC &);
+
+    template
+    void project<deal_II_dimension, VEC>(
+        const MatrixFree<deal_II_dimension,VEC::value_type> &matrix_free,
+        const ConstraintMatrix &constraints,
+        const unsigned int,
+        const std_cxx11::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)>,
+        VEC &);
+
+    \}
+
+}
+
diff --git a/tests/numerics/project_parallel_qp_common.h b/tests/numerics/project_parallel_qp_common.h
new file mode 100644 (file)
index 0000000..4d9cb56
--- /dev/null
@@ -0,0 +1,245 @@
+// ---------------------------------------------------------------------
+//
+// 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 whether an element of polynomial order p can
+// represent functions of order q for projection from quadrature points.
+
+#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/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/lac/generic_linear_algebra.h>
+#include <deal.II/base/quadrature_point_data.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+#include <fstream>
+#include <vector>
+
+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 = 0) 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;
+  }
+
+  VectorizedArray<double> value(const Point<dim, VectorizedArray<double>> &p_vec) const
+  {
+    VectorizedArray<double> res = make_vectorized_array (0.);
+    Point<dim> p;
+    for (unsigned int v = 0; v < VectorizedArray<double>::n_array_elements; ++v)
+      {
+        for (unsigned int d = 0; d < dim; d++)
+          p[d] = p_vec[d][v];
+        res[v] = value(p);
+      }
+    return res;
+  }
+
+
+private:
+  const unsigned int q;
+};
+
+struct QData
+{
+  double density;
+};
+
+
+template <typename VectorType, int dim>
+void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
+                 const FiniteElement<dim> &fe,
+                 const unsigned int        p)
+{
+  AssertThrow(fe.n_components() ==1,
+              ExcNotImplemented());
+  DoFHandler<dim>        dof_handler(triangulation);
+  dof_handler.distribute_dofs (fe);
+
+  deallog << "n_cells=" << triangulation.n_global_active_cells() << std::endl;
+  deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl;
+
+  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 ();
+
+  for (unsigned int q=0; q<=p; ++q)
+    {
+      // setup quadrature data:
+      F<dim> function(q, fe.n_components());
+      VectorType field_relevant;
+      double field_l2_norm = 0.;
+      {
+        QGauss<dim> quadrature_formula(p+2);
+        CellDataStorage<typename parallel::distributed::Triangulation<dim,dim>::cell_iterator,QData> qp_manager;
+        {
+          FEValues<dim> fe_values (fe, quadrature_formula,
+                                   update_quadrature_points);
+
+          const unsigned int n_q_points  = quadrature_formula.size();
+
+          typename DoFHandler<dim>::active_cell_iterator
+          cell = dof_handler.begin_active(),
+          endc = dof_handler.end();
+          for (; cell!=endc; ++cell)
+            if (cell->is_locally_owned())
+              {
+                fe_values.reinit (cell);
+                qp_manager.initialize(cell, quadrature_formula.size());
+                std::vector<std::shared_ptr<QData> > qpd = qp_manager.get_data(cell);
+                for (unsigned int q=0; q<n_q_points; ++q)
+                  qpd[q]->density = function.value(fe_values.quadrature_point(q));
+              }
+        }
+
+        VectorType field(dof_handler.locally_owned_dofs(),
+                         MPI_COMM_WORLD);
+        VectorTools::project
+        (MappingQGeneric<dim>(1),
+         dof_handler,
+         constraints,
+         quadrature_formula,
+         [=] (const typename DoFHandler<dim>::active_cell_iterator & cell, const unsigned int q) -> double { return qp_manager.get_data(cell)[q]->density; },
+         field);
+
+        field_relevant.reinit(dof_handler.locally_owned_dofs(),
+                              locally_relevant_dofs,
+                              MPI_COMM_WORLD);
+        field_relevant = field;
+
+        field_l2_norm = field.l2_norm();
+      }
+
+      // L2 norm of the difference between FE field and the function
+      double L2_norm = 0.;
+      {
+        QGauss<dim> quadrature_formula_error(std::max(p,q)+1);
+        FEValues<dim> fe_values (fe, quadrature_formula_error,
+                                 update_values  | update_quadrature_points | update_JxW_values);
+
+        const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+        const unsigned int   n_q_points    = quadrature_formula_error.size();
+        std::vector<double>  values (n_q_points);
+
+        typename DoFHandler<dim>::active_cell_iterator
+        cell = dof_handler.begin_active(),
+        endc = dof_handler.end();
+        for (; cell!=endc; ++cell)
+          if (cell->is_locally_owned())
+            {
+              fe_values.reinit (cell);
+              fe_values.get_function_values(field_relevant,
+                                            values);
+
+              for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                L2_norm += Utilities::fixed_power<2>(values[q_point] -  function.value(fe_values.quadrature_point(q_point))) *
+                           fe_values.JxW(q_point);
+
+            }
+      }
+
+      L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD);
+      L2_norm = std::sqrt(L2_norm);
+
+      deallog << fe.get_name() << ", P_" << q
+              << ", rel. error=" << L2_norm / field_l2_norm
+              << std::endl;
+
+      if (q<=p)
+        if (L2_norm > 1e-10*field_l2_norm)
+          deallog << "Projection failed with relative error "
+                  << L2_norm / field_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 <typename VectorType, int dim>
+void test_no_hanging_nodes (const FiniteElement<dim> &fe,
+                            const unsigned int        p)
+{
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (3);
+
+  do_project<VectorType> (triangulation, fe, p);
+}
+
+
+
+// same test as above, but this time with a mesh that has hanging nodes
+template <typename VectorType, int dim>
+void test_with_hanging_nodes (const FiniteElement<dim> &fe,
+                              const unsigned int        p)
+{
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (2);
+  if (triangulation.begin_active()->is_locally_owned())
+    triangulation.begin_active()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement ();
+  triangulation.refine_global (1);
+
+  do_project<VectorType> (triangulation, fe, p);
+}
+
diff --git a/tests/numerics/project_parallel_qp_q_01.cc b/tests/numerics/project_parallel_qp_q_01.cc
new file mode 100644 (file)
index 0000000..06027d6
--- /dev/null
@@ -0,0 +1,45 @@
+// ---------------------------------------------------------------------
+//
+// 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 VectorTools::project_parallel() of quadrature data for
+// FE_Q on a mesh without hanging nodes.
+
+#include "project_parallel_qp_common.h"
+
+namespace LA
+{
+  using namespace ::LinearAlgebraTrilinos;
+}
+
+template <int dim>
+void test ()
+{
+  for (unsigned int p=1; p<7-dim; ++p)
+    test_no_hanging_nodes<LA::MPI::Vector> (FE_Q<dim>(p), p);
+}
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/numerics/project_parallel_qp_q_01.mpirun=3.with_cxx11=on.with_trilinos=true.output b/tests/numerics/project_parallel_qp_q_01.mpirun=3.with_cxx11=on.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..d338fe1
--- /dev/null
@@ -0,0 +1,38 @@
+
+DEAL::n_cells=64
+DEAL::n_dofs=81
+DEAL::FE_Q<2>(1), P_0, rel. error=0
+DEAL::FE_Q<2>(1), P_1, rel. error=0
+DEAL::n_cells=64
+DEAL::n_dofs=289
+DEAL::FE_Q<2>(2), P_0, rel. error=0
+DEAL::FE_Q<2>(2), P_1, rel. error=0
+DEAL::FE_Q<2>(2), P_2, rel. error=0
+DEAL::n_cells=64
+DEAL::n_dofs=625
+DEAL::FE_Q<2>(3), P_0, rel. error=0
+DEAL::FE_Q<2>(3), P_1, rel. error=0
+DEAL::FE_Q<2>(3), P_2, rel. error=0
+DEAL::FE_Q<2>(3), P_3, rel. error=0
+DEAL::n_cells=64
+DEAL::n_dofs=1089
+DEAL::FE_Q<2>(4), P_0, rel. error=0
+DEAL::FE_Q<2>(4), P_1, rel. error=0
+DEAL::FE_Q<2>(4), P_2, rel. error=0
+DEAL::FE_Q<2>(4), P_3, rel. error=0
+DEAL::FE_Q<2>(4), P_4, rel. error=0
+DEAL::n_cells=512
+DEAL::n_dofs=729
+DEAL::FE_Q<3>(1), P_0, rel. error=0
+DEAL::FE_Q<3>(1), P_1, rel. error=0
+DEAL::n_cells=512
+DEAL::n_dofs=4913
+DEAL::FE_Q<3>(2), P_0, rel. error=0
+DEAL::FE_Q<3>(2), P_1, rel. error=0
+DEAL::FE_Q<3>(2), P_2, rel. error=0
+DEAL::n_cells=512
+DEAL::n_dofs=15625
+DEAL::FE_Q<3>(3), P_0, rel. error=0
+DEAL::FE_Q<3>(3), P_1, rel. error=0
+DEAL::FE_Q<3>(3), P_2, rel. error=0
+DEAL::FE_Q<3>(3), P_3, rel. error=0
diff --git a/tests/numerics/project_parallel_qp_q_02.cc b/tests/numerics/project_parallel_qp_q_02.cc
new file mode 100644 (file)
index 0000000..230793a
--- /dev/null
@@ -0,0 +1,45 @@
+// ---------------------------------------------------------------------
+//
+// 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 VectorTools::project_parallel() of quadrature data for
+// FE_Q on a mesh with hanging nodes.
+
+#include "project_parallel_qp_common.h"
+
+namespace LA
+{
+  using namespace ::LinearAlgebraTrilinos;
+}
+
+template <int dim>
+void test ()
+{
+  for (unsigned int p=1; p<7-dim; ++p)
+    test_with_hanging_nodes<LA::MPI::Vector> (FE_Q<dim>(p), p);
+}
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/numerics/project_parallel_qp_q_02.mpirun=3.with_cxx11=on.with_trilinos=true.output b/tests/numerics/project_parallel_qp_q_02.mpirun=3.with_cxx11=on.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..5a79d55
--- /dev/null
@@ -0,0 +1,38 @@
+
+DEAL::n_cells=76
+DEAL::n_dofs=97
+DEAL::FE_Q<2>(1), P_0, rel. error=0
+DEAL::FE_Q<2>(1), P_1, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=349
+DEAL::FE_Q<2>(2), P_0, rel. error=0
+DEAL::FE_Q<2>(2), P_1, rel. error=0
+DEAL::FE_Q<2>(2), P_2, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=753
+DEAL::FE_Q<2>(3), P_0, rel. error=0
+DEAL::FE_Q<2>(3), P_1, rel. error=0
+DEAL::FE_Q<2>(3), P_2, rel. error=0
+DEAL::FE_Q<2>(3), P_3, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=1309
+DEAL::FE_Q<2>(4), P_0, rel. error=0
+DEAL::FE_Q<2>(4), P_1, rel. error=0
+DEAL::FE_Q<2>(4), P_2, rel. error=0
+DEAL::FE_Q<2>(4), P_3, rel. error=0
+DEAL::FE_Q<2>(4), P_4, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=827
+DEAL::FE_Q<3>(1), P_0, rel. error=0
+DEAL::FE_Q<3>(1), P_1, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=5559
+DEAL::FE_Q<3>(2), P_0, rel. error=0
+DEAL::FE_Q<3>(2), P_1, rel. error=0
+DEAL::FE_Q<3>(2), P_2, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=17587
+DEAL::FE_Q<3>(3), P_0, rel. error=0
+DEAL::FE_Q<3>(3), P_1, rel. error=0
+DEAL::FE_Q<3>(3), P_2, rel. error=0
+DEAL::FE_Q<3>(3), P_3, rel. error=0
diff --git a/tests/numerics/project_parallel_qpmf_common.h b/tests/numerics/project_parallel_qpmf_common.h
new file mode 100644 (file)
index 0000000..cd832d1
--- /dev/null
@@ -0,0 +1,236 @@
+// ---------------------------------------------------------------------
+//
+// 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 whether an element of polynomial order p can
+// represent functions of order q for projection from quadrature points of
+// matrix-free approach.
+
+#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/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+#include <fstream>
+#include <vector>
+
+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 = 0) 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;
+  }
+
+  VectorizedArray<double> value(const Point<dim, VectorizedArray<double>> &p_vec) const
+  {
+    VectorizedArray<double> res = make_vectorized_array (0.);
+    Point<dim> p;
+    for (unsigned int v = 0; v < VectorizedArray<double>::n_array_elements; ++v)
+      {
+        for (unsigned int d = 0; d < dim; d++)
+          p[d] = p_vec[d][v];
+        res[v] = value(p);
+      }
+    return res;
+  }
+
+
+private:
+  const unsigned int q;
+};
+
+
+template <int fe_degree, int n_q_points_1d, int dim>
+void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
+                 const FiniteElement<dim> &fe,
+                 const unsigned int        p)
+{
+  AssertThrow(fe.n_components() ==1,
+              ExcNotImplemented());
+  DoFHandler<dim>        dof_handler(triangulation);
+  dof_handler.distribute_dofs (fe);
+
+  deallog << "n_cells=" << triangulation.n_global_active_cells() << std::endl;
+  deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl;
+
+  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 ();
+
+  QGauss<1> quadrature_formula_1d(n_q_points_1d);
+
+  Table<2, VectorizedArray<double> > qp_data;
+
+  typename MatrixFree<dim,double>::AdditionalData additional_data;
+  additional_data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::partition_color;
+  additional_data.mpi_communicator = MPI_COMM_WORLD;
+  additional_data.mapping_update_flags = update_values | update_JxW_values | update_quadrature_points;
+  MatrixFree<dim,double>  data;
+  data.reinit (dof_handler, constraints, quadrature_formula_1d, additional_data);
+
+  for (unsigned int q=0; q<=p; ++q)
+    {
+      // setup quadrature data:
+      F<dim> function(q, fe.n_components());
+
+      // initialize a quadrature data
+      {
+        FEEvaluation<dim,fe_degree,n_q_points_1d,1,double> fe_eval(data);
+        const unsigned int n_cells = data.n_macro_cells();
+        const unsigned int n_q_points = fe_eval.n_q_points;
+
+        qp_data.reinit(n_cells, n_q_points);
+        for (unsigned int cell=0; cell<n_cells; ++cell)
+          {
+            fe_eval.reinit(cell);
+            for (unsigned int q=0; q<n_q_points; ++q)
+              qp_data(cell,q) = function.value(fe_eval.quadrature_point(q));
+          }
+      }
+
+      LinearAlgebra::distributed::Vector<double> field;
+      data.initialize_dof_vector(field);
+      VectorTools::project
+      (data,
+       constraints,
+       n_q_points_1d,
+       [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray<double> { return qp_data(cell,q); },
+       field);
+
+      field.update_ghost_values();
+
+      const double field_l2_norm = field.l2_norm();
+
+      // L2 norm of the difference between FE field and the function
+      double L2_norm = 0.;
+      {
+        QGauss<dim> quadrature_formula_error(std::max(p,q)+1);
+        FEValues<dim> fe_values (fe, quadrature_formula_error,
+                                 update_values  | update_quadrature_points | update_JxW_values);
+
+        const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+        const unsigned int   n_q_points    = quadrature_formula_error.size();
+        std::vector<double>  values (n_q_points);
+
+        typename DoFHandler<dim>::active_cell_iterator
+        cell = dof_handler.begin_active(),
+        endc = dof_handler.end();
+        for (; cell!=endc; ++cell)
+          if (cell->is_locally_owned())
+            {
+              fe_values.reinit (cell);
+              fe_values.get_function_values(field,
+                                            values);
+
+              for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                L2_norm += Utilities::fixed_power<2>(values[q_point] -  function.value(fe_values.quadrature_point(q_point))) *
+                           fe_values.JxW(q_point);
+
+            }
+      }
+
+      L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD);
+      L2_norm = std::sqrt(L2_norm);
+
+      deallog << fe.get_name() << ", P_" << q
+              << ", rel. error=" << L2_norm / field_l2_norm
+              << std::endl;
+
+      if (q<=p)
+        if (L2_norm > 1e-10*field_l2_norm)
+          deallog << "Projection failed with relative error "
+                  << L2_norm / field_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 fe_degree, int n_q_points_1d, int dim>
+void test_no_hanging_nodes (const FiniteElement<dim> &fe,
+                            const unsigned int        p)
+{
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (3);
+
+  do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+}
+
+
+
+// same test as above, but this time with a mesh that has hanging nodes
+template <int fe_degree, int n_q_points_1d, int dim>
+void test_with_hanging_nodes (const FiniteElement<dim> &fe,
+                              const unsigned int        p)
+{
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (2);
+  if (triangulation.begin_active()->is_locally_owned())
+    triangulation.begin_active()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement ();
+  triangulation.refine_global (1);
+
+  do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+}
+
diff --git a/tests/numerics/project_parallel_qpmf_q_01.cc b/tests/numerics/project_parallel_qpmf_q_01.cc
new file mode 100644 (file)
index 0000000..3979347
--- /dev/null
@@ -0,0 +1,43 @@
+// ---------------------------------------------------------------------
+//
+// 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 VectorTools::project_parallel() for matrix-free quadrature data for
+// FE_Q on a mesh without hanging nodes.
+
+
+#include "project_parallel_qpmf_common.h"
+
+template <int dim>
+void test ()
+{
+  test_no_hanging_nodes<1, 2,dim> (FE_Q<dim>(1), 1);
+  test_no_hanging_nodes<2, 3,dim> (FE_Q<dim>(2), 2);
+  test_no_hanging_nodes<3, 4,dim> (FE_Q<dim>(3), 3);
+}
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv,
+                                                       numbers::invalid_unsigned_int);
+
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/numerics/project_parallel_qpmf_q_01.mpirun=3.with_cxx11=on.output b/tests/numerics/project_parallel_qpmf_q_01.mpirun=3.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..ff5f4b5
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::n_cells=64
+DEAL::n_dofs=81
+DEAL::FE_Q<2>(1), P_0, rel. error=0
+DEAL::FE_Q<2>(1), P_1, rel. error=0
+DEAL::n_cells=64
+DEAL::n_dofs=289
+DEAL::FE_Q<2>(2), P_0, rel. error=0
+DEAL::FE_Q<2>(2), P_1, rel. error=0
+DEAL::FE_Q<2>(2), P_2, rel. error=0
+DEAL::n_cells=64
+DEAL::n_dofs=625
+DEAL::FE_Q<2>(3), P_0, rel. error=0
+DEAL::FE_Q<2>(3), P_1, rel. error=0
+DEAL::FE_Q<2>(3), P_2, rel. error=0
+DEAL::FE_Q<2>(3), P_3, rel. error=0
+DEAL::n_cells=512
+DEAL::n_dofs=729
+DEAL::FE_Q<3>(1), P_0, rel. error=0
+DEAL::FE_Q<3>(1), P_1, rel. error=0
+DEAL::n_cells=512
+DEAL::n_dofs=4913
+DEAL::FE_Q<3>(2), P_0, rel. error=0
+DEAL::FE_Q<3>(2), P_1, rel. error=0
+DEAL::FE_Q<3>(2), P_2, rel. error=0
+DEAL::n_cells=512
+DEAL::n_dofs=15625
+DEAL::FE_Q<3>(3), P_0, rel. error=0
+DEAL::FE_Q<3>(3), P_1, rel. error=0
+DEAL::FE_Q<3>(3), P_2, rel. error=0
+DEAL::FE_Q<3>(3), P_3, rel. error=0
diff --git a/tests/numerics/project_parallel_qpmf_q_02.cc b/tests/numerics/project_parallel_qpmf_q_02.cc
new file mode 100644 (file)
index 0000000..5770980
--- /dev/null
@@ -0,0 +1,42 @@
+// ---------------------------------------------------------------------
+//
+// 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 VectorTools::project_parallel() for matrix-free quadrature data for
+// FE_Q on a mesh with hanging nodes.
+
+#include "project_parallel_qpmf_common.h"
+
+template <int dim>
+void test ()
+{
+  test_with_hanging_nodes<1, 2,dim> (FE_Q<dim>(1), 1);
+  test_with_hanging_nodes<2, 3,dim> (FE_Q<dim>(2), 2);
+  test_with_hanging_nodes<3, 4,dim> (FE_Q<dim>(3), 3);
+}
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv,
+                                                       numbers::invalid_unsigned_int);
+
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/numerics/project_parallel_qpmf_q_02.mpirun=3.with_cxx11=on.output b/tests/numerics/project_parallel_qpmf_q_02.mpirun=3.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..78c248a
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::n_cells=76
+DEAL::n_dofs=97
+DEAL::FE_Q<2>(1), P_0, rel. error=0
+DEAL::FE_Q<2>(1), P_1, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=349
+DEAL::FE_Q<2>(2), P_0, rel. error=0
+DEAL::FE_Q<2>(2), P_1, rel. error=0
+DEAL::FE_Q<2>(2), P_2, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=753
+DEAL::FE_Q<2>(3), P_0, rel. error=0
+DEAL::FE_Q<2>(3), P_1, rel. error=0
+DEAL::FE_Q<2>(3), P_2, rel. error=0
+DEAL::FE_Q<2>(3), P_3, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=827
+DEAL::FE_Q<3>(1), P_0, rel. error=0
+DEAL::FE_Q<3>(1), P_1, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=5559
+DEAL::FE_Q<3>(2), P_0, rel. error=0
+DEAL::FE_Q<3>(2), P_1, rel. error=0
+DEAL::FE_Q<3>(2), P_2, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=17587
+DEAL::FE_Q<3>(3), P_0, rel. error=0
+DEAL::FE_Q<3>(3), P_1, rel. error=0
+DEAL::FE_Q<3>(3), P_2, rel. error=0
+DEAL::FE_Q<3>(3), P_3, rel. error=0

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.