]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Enabel BlockVector for FEEvaluation::gather_evalute() 10865/head
authorPeter Munch <peterrmuench@gmail.com>
Mon, 31 Aug 2020 14:27:05 +0000 (16:27 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 3 Sep 2020 10:48:07 +0000 (12:48 +0200)
doc/news/changes/minor/20200831Munch [new file with mode: 0644]
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
tests/matrix_free/matrix_vector_faces_29.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20200831Munch b/doc/news/changes/minor/20200831Munch
new file mode 100644 (file)
index 0000000..bb36682
--- /dev/null
@@ -0,0 +1,4 @@
+New: The methods FEEvaluation/FEEFacevaluation::gather_evaluate/integrate_scatter() 
+can now also accept block vectors.
+<br>
+(Peter Munch, Magdalena Schreter, Martin Kronbichler, 2020/08/31)
index 4451027de45a7f0636fe133e7648300b382f1441..120f8fe705b7774c358ef21f97621e52e94087fa 100644 (file)
@@ -26,6 +26,7 @@
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/type_traits.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -2041,10 +2042,44 @@ namespace internal
           data, temp1, values_array, integrate_gradients, face_no);
     }
 
-    template <std::size_t n_face_orientations>
+    template <std::size_t n_face_orientations,
+              typename VectorType,
+              typename std::enable_if<
+                !internal::has_begin<VectorType>::value ||
+                  !std::is_same<decltype(std::declval<VectorType>().begin()),
+                                Number2 *>::value,
+                VectorType>::type * = nullptr>
+    static bool
+    gather_evaluate(const VectorType &,
+                    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &,
+                    const MatrixFreeFunctions::DoFInfo &,
+                    VectorizedArrayType *,
+                    VectorizedArrayType *,
+                    VectorizedArrayType *,
+                    const bool,
+                    const bool,
+                    const unsigned int,
+                    const unsigned int,
+                    const std::array<unsigned int, n_face_orientations>,
+                    const std::array<unsigned int, n_face_orientations>,
+                    const unsigned int,
+                    const MatrixFreeFunctions::DoFInfo::DoFAccessIndex,
+                    const std::array<unsigned int, n_face_orientations>,
+                    const Table<2, unsigned int> &)
+    {
+      return false;
+    }
+
+    template <std::size_t n_face_orientations,
+              typename VectorType,
+              typename std::enable_if<
+                internal::has_begin<VectorType>::value &&
+                  std::is_same<decltype(std::declval<VectorType>().begin()),
+                               Number2 *>::value,
+                VectorType>::type * = nullptr>
     static bool
     gather_evaluate(
-      const Number2 *                                            src_ptr,
+      const VectorType &                                         src,
       const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
       const MatrixFreeFunctions::DoFInfo &                       dof_info,
       VectorizedArrayType *                                      values_quad,
@@ -2061,9 +2096,13 @@ namespace internal
       const std::array<unsigned int, n_face_orientations> face_orientations,
       const Table<2, unsigned int> &                      orientation_map)
     {
+      static_assert((std::is_same<decltype(std::declval<VectorType>().begin()),
+                                  Number2 *>::value),
+                    "Types do not match.");
+
       return process_and_io( //
         false /*=evaluate*/,
-        src_ptr,
+        src.begin(),
         data,
         dof_info,
         values_quad,
@@ -2165,10 +2204,63 @@ namespace internal
         });
     }
 
-    template <std::size_t n_face_orientations>
+    template <std::size_t n_face_orientations,
+              typename VectorType,
+              typename std::enable_if<
+                !internal::has_begin<VectorType>::value ||
+                  !std::is_same<decltype(std::declval<VectorType>().begin()),
+                                Number2 *>::value,
+                VectorType>::type * = nullptr>
+    static bool
+    integrate_scatter(
+      VectorType &,
+      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+      const MatrixFreeFunctions::DoFInfo &,
+      VectorizedArrayType *values_array,
+      VectorizedArrayType *values_quad,
+      VectorizedArrayType *gradients_quad,
+      VectorizedArrayType *scratch_data,
+      const bool           integrate_values,
+      const bool           integrate_gradients,
+      const unsigned int,
+      const unsigned int,
+      const std::array<unsigned int, n_face_orientations>,
+      const std::array<unsigned int, n_face_orientations> face_nos,
+      const unsigned int                                  subface_index,
+      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex,
+      const std::array<unsigned int, n_face_orientations> face_orientations,
+      const Table<2, unsigned int> &                      orientation_map)
+    {
+      AssertDimension(face_nos.size(), 1);
+      AssertDimension(face_orientations.size(), 1);
+
+      // for block vectors simply integrate
+      integrate(data,
+                values_array,
+                values_quad,
+                gradients_quad,
+                scratch_data,
+                integrate_values,
+                integrate_gradients,
+                face_nos[0],
+                subface_index,
+                face_orientations[0],
+                orientation_map);
+
+      // default vector access
+      return false;
+    }
+
+    template <std::size_t n_face_orientations,
+              typename VectorType,
+              typename std::enable_if<
+                internal::has_begin<VectorType>::value &&
+                  std::is_same<decltype(std::declval<VectorType>().begin()),
+                               Number2 *>::value,
+                VectorType>::type * = nullptr>
     static bool
     integrate_scatter(
-      Number2 *                                                  dst_ptr,
+      VectorType &                                               dst,
       const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
       const MatrixFreeFunctions::DoFInfo &                       dof_info,
       VectorizedArrayType *                                      values_array,
@@ -2186,9 +2278,13 @@ namespace internal
       const std::array<unsigned int, n_face_orientations> face_orientations,
       const Table<2, unsigned int> &                      orientation_map)
     {
+      static_assert((std::is_same<decltype(std::declval<VectorType>().begin()),
+                                  Number2 *>::value),
+                    "Types do not match.");
+
       return process_and_io( //
         true /*=integrate*/,
-        dst_ptr,
+        dst.begin(),
         data,
         dof_info,
         values_quad,
index 560d7954effeae34d4de1099bb0fcbc208123c28..b91ba87259d86071a206d47377bf340b963b7f7c 100644 (file)
@@ -7315,6 +7315,192 @@ FEEvaluation<
   gather_evaluate(input_vector, flag);
 }
 
+
+namespace internal
+{
+  /**
+   * Implementation for vectors that have the begin() methods.
+   */
+  template <typename Number,
+            typename VectorizedArrayType,
+            typename VectorType,
+            typename T,
+            typename std::enable_if<
+              internal::has_begin<VectorType>::value &&
+                std::is_same<decltype(std::declval<VectorType>().begin()),
+                             Number *>::value,
+              VectorType>::type * = nullptr>
+  bool
+  try_gather_evaluate_inplace(
+    T                                             phi,
+    const VectorType &                            input_vector,
+    const unsigned int                            cell,
+    const unsigned int                            active_fe_index,
+    const unsigned int                            first_selected_component,
+    const internal::MatrixFreeFunctions::DoFInfo *dof_info,
+    const EvaluationFlags::EvaluationFlags        evaluation_flag)
+  {
+    // If the index storage is interleaved and contiguous and the vector storage
+    // has the correct alignment, we can directly pass the pointer into the
+    // vector to the evaluate() call, without reading the vector entries into a
+    // separate data field. This saves some operations.
+    if (std::is_same<typename VectorType::value_type, Number>::value &&
+        dof_info->index_storage_variants
+            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
+          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+            interleaved_contiguous &&
+        reinterpret_cast<std::size_t>(
+          input_vector.begin() +
+          dof_info->dof_indices_contiguous
+            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+            [cell * VectorizedArrayType::size()]) %
+            sizeof(VectorizedArrayType) ==
+          0)
+      {
+        const VectorizedArrayType *vec_values =
+          reinterpret_cast<const VectorizedArrayType *>(
+            input_vector.begin() +
+            dof_info->dof_indices_contiguous
+              [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+              [cell * VectorizedArrayType::size()] +
+            dof_info->component_dof_indices_offset[active_fe_index]
+                                                  [first_selected_component] *
+              VectorizedArrayType::size());
+
+        phi->evaluate(vec_values, evaluation_flag);
+
+        return true;
+      }
+
+    return false;
+  }
+
+  /**
+   * Implementation for all other vectors like block vectors.
+   */
+  template <typename Number,
+            typename VectorizedArrayType,
+            typename VectorType,
+            typename T,
+            typename std::enable_if<
+              !internal::has_begin<VectorType>::value ||
+                !std::is_same<decltype(std::declval<VectorType>().begin()),
+                              Number *>::value,
+              VectorType>::type * = nullptr>
+  bool
+  try_gather_evaluate_inplace(T,
+                              const VectorType &,
+                              const unsigned int,
+                              const unsigned int,
+                              const unsigned int,
+                              const internal::MatrixFreeFunctions::DoFInfo *,
+                              const EvaluationFlags::EvaluationFlags)
+  {
+    return false;
+  }
+
+  /**
+   * Implementation for vectors that have the begin() methods.
+   */
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType,
+            typename std::enable_if<
+              internal::has_begin<VectorType>::value &&
+                std::is_same<decltype(std::declval<VectorType>().begin()),
+                             Number *>::value,
+              VectorType>::type * = nullptr>
+  bool
+  try_integrate_scatter_inplace(
+    VectorType &                                  destination,
+    const unsigned int                            cell,
+    const unsigned int                            n_components,
+    const unsigned int                            active_fe_index,
+    const unsigned int                            first_selected_component,
+    const internal::MatrixFreeFunctions::DoFInfo *dof_info,
+    VectorizedArrayType *                         values_quad,
+    VectorizedArrayType *                         gradients_quad,
+    VectorizedArrayType *                         scratch_data,
+    const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data,
+    const EvaluationFlags::EvaluationFlags integration_flag)
+  {
+    // If the index storage is interleaved and contiguous and the vector storage
+    // has the correct alignment, we can directly pass the pointer into the
+    // vector to the integrate() call, without writing temporary results into a
+    // separate data field that will later be added into the vector. This saves
+    // some operations.
+    if (std::is_same<typename VectorType::value_type, Number>::value &&
+        dof_info->index_storage_variants
+            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
+          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+            interleaved_contiguous &&
+        reinterpret_cast<std::size_t>(
+          destination.begin() +
+          dof_info->dof_indices_contiguous
+            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+            [cell * VectorizedArrayType::size()]) %
+            sizeof(VectorizedArrayType) ==
+          0)
+      {
+        VectorizedArrayType *vec_values =
+          reinterpret_cast<VectorizedArrayType *>(
+            destination.begin() +
+            dof_info->dof_indices_contiguous
+              [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+              [cell * VectorizedArrayType::size()] +
+            dof_info->component_dof_indices_offset[active_fe_index]
+                                                  [first_selected_component] *
+              VectorizedArrayType::size());
+        SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
+          integrate(n_components,
+                    integration_flag,
+                    *data,
+                    vec_values,
+                    values_quad,
+                    gradients_quad,
+                    scratch_data,
+                    true);
+        return true;
+      }
+
+    return false;
+  }
+
+  /**
+   * Implementation for all other vectors like block vectors.
+   */
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType,
+            typename std::enable_if<
+              !internal::has_begin<VectorType>::value ||
+                !std::is_same<decltype(std::declval<VectorType>().begin()),
+                              Number *>::value,
+              VectorType>::type * = nullptr>
+  bool
+  try_integrate_scatter_inplace(
+    VectorType &,
+    const unsigned int,
+    const unsigned int,
+    const unsigned int,
+    const unsigned int,
+    const internal::MatrixFreeFunctions::DoFInfo *,
+    const VectorizedArrayType *,
+    const VectorizedArrayType *,
+    const VectorizedArrayType *,
+    const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *,
+    const EvaluationFlags::EvaluationFlags)
+  {
+    return false;
+  }
+} // namespace internal
+
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
@@ -7332,37 +7518,14 @@ FEEvaluation<dim,
   gather_evaluate(const VectorType &                     input_vector,
                   const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
-  // If the index storage is interleaved and contiguous and the vector storage
-  // has the correct alignment, we can directly pass the pointer into the
-  // vector to the evaluate() call, without reading the vector entries into a
-  // separate data field. This saves some operations.
-  if (std::is_same<typename VectorType::value_type, Number>::value &&
-      this->dof_info->index_storage_variants
-          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
-                            IndexStorageVariants::interleaved_contiguous &&
-      reinterpret_cast<std::size_t>(
-        input_vector.begin() +
-        this->dof_info->dof_indices_contiguous
-          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-          [this->cell * VectorizedArrayType::size()]) %
-          sizeof(VectorizedArrayType) ==
-        0)
-    {
-      const VectorizedArrayType *vec_values =
-        reinterpret_cast<const VectorizedArrayType *>(
-          input_vector.begin() +
-          this->dof_info->dof_indices_contiguous
-            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-            [this->cell * VectorizedArrayType::size()] +
-          this->dof_info
-              ->component_dof_indices_offset[this->active_fe_index]
-                                            [this->first_selected_component] *
-            VectorizedArrayType::size());
-
-      evaluate(vec_values, evaluation_flag);
-    }
-  else
+  if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
+        this,
+        input_vector,
+        this->cell,
+        this->active_fe_index,
+        this->first_selected_component,
+        this->dof_info,
+        evaluation_flag) == false)
     {
       this->read_dof_values(input_vector);
       evaluate(this->begin_dof_values(), evaluation_flag);
@@ -7540,44 +7703,22 @@ FEEvaluation<dim,
   integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                     VectorType &                           destination)
 {
-  // If the index storage is interleaved and contiguous and the vector storage
-  // has the correct alignment, we can directly pass the pointer into the
-  // vector to the integrate() call, without writing temporary results into a
-  // separate data field that will later be added into the vector. This saves
-  // some operations.
-  if (std::is_same<typename VectorType::value_type, Number>::value &&
-      this->dof_info->index_storage_variants
-          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
-                            IndexStorageVariants::interleaved_contiguous &&
-      reinterpret_cast<std::size_t>(
-        destination.begin() +
-        this->dof_info->dof_indices_contiguous
-          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-          [this->cell * VectorizedArrayType::size()]) %
-          sizeof(VectorizedArrayType) ==
-        0)
-    {
-      VectorizedArrayType *vec_values = reinterpret_cast<VectorizedArrayType *>(
-        destination.begin() +
-        this->dof_info->dof_indices_contiguous
-          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-          [this->cell * VectorizedArrayType::size()] +
-        this->dof_info
-            ->component_dof_indices_offset[this->active_fe_index]
-                                          [this->first_selected_component] *
-          VectorizedArrayType::size());
-      SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
-        integrate(n_components,
-                  integration_flag,
-                  *this->data,
-                  vec_values,
-                  this->values_quad,
-                  this->gradients_quad,
-                  this->scratch_data,
-                  true);
-    }
-  else
+  if (internal::try_integrate_scatter_inplace<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              Number,
+                                              VectorizedArrayType>(
+        destination,
+        this->cell,
+        n_components,
+        this->active_fe_index,
+        this->first_selected_component,
+        this->dof_info,
+        this->values_quad,
+        this->gradients_quad,
+        this->scratch_data,
+        this->data,
+        integration_flag) == false)
     {
       integrate(integration_flag, this->begin_dof_values());
       this->distribute_local_to_global(destination);
@@ -8065,15 +8206,6 @@ FEFaceEvaluation<dim,
   gather_evaluate(const VectorType &                     input_vector,
                   const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
-  static_assert(internal::has_begin<VectorType>::value &&
-                  (std::is_same<decltype(std::declval<VectorType>().begin()),
-                                double *>::value ||
-                   std::is_same<decltype(std::declval<VectorType>().begin()),
-                                float *>::value),
-                "This function requires a vector type with begin() function "
-                "evaluating to a pointer to basic number (float,double). "
-                "Use read_dof_values() followed by evaluate() instead.");
-
   Assert(
     (evaluation_flag &
      ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
@@ -8095,7 +8227,7 @@ FEFaceEvaluation<dim,
                                                   n_components,
                                                   Number,
                                                   VectorizedArrayType>::
-          gather_evaluate(input_vector.begin(),
+          gather_evaluate(input_vector,
                           *this->data,
                           *this->dof_info,
                           this->begin_values(),
@@ -8121,7 +8253,7 @@ FEFaceEvaluation<dim,
                                                   n_components,
                                                   Number,
                                                   VectorizedArrayType>::
-          gather_evaluate(input_vector.begin(),
+          gather_evaluate(input_vector,
                           *this->data,
                           *this->dof_info,
                           this->begin_values(),
@@ -8202,16 +8334,6 @@ FEFaceEvaluation<dim,
   integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
                     VectorType &                           destination)
 {
-  static_assert(internal::has_begin<VectorType>::value &&
-                  (std::is_same<decltype(std::declval<VectorType>().begin()),
-                                double *>::value ||
-                   std::is_same<decltype(std::declval<VectorType>().begin()),
-                                float *>::value),
-                "This function requires a vector type with begin() function "
-                "evaluating to a pointer to basic number (float,double). "
-                "Use integrate() followed by distribute_local_to_global() "
-                "instead.");
-
   Assert((this->dof_access_index ==
             internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
           this->is_interior_face == false) == false,
@@ -8223,7 +8345,7 @@ FEFaceEvaluation<dim,
                                           n_components,
                                           Number,
                                           VectorizedArrayType>::
-        integrate_scatter(destination.begin(),
+        integrate_scatter(destination,
                           *this->data,
                           *this->dof_info,
                           this->begin_dof_values(),
diff --git a/tests/matrix_free/matrix_vector_faces_29.cc b/tests/matrix_free/matrix_vector_faces_29.cc
new file mode 100644 (file)
index 0000000..e1706bb
--- /dev/null
@@ -0,0 +1,148 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// similar to matrix_vector_faces_13 (renumbering of degrees of freedom to
+// better vectorized access, gather_evaluate and integrate_scatter) but with
+// block vectors
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/la_parallel_block_vector.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_faces_common.h"
+
+
+
+template <int dim, int fe_degree_>
+void
+test()
+{
+  // raise element degree by one to test quadratic and cubic shape functions
+  // rather than linears and quadratics according to the main function in
+  // matrix_vector_faces_common.h
+
+  const unsigned int                        fe_degree = fe_degree_ + 1;
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+  for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
+       cell != tria.end();
+       ++cell)
+    for (const unsigned int f : GeometryInfo<dim>::face_indices())
+      if (cell->at_boundary(f))
+        cell->face(f)->set_all_boundary_ids(f);
+  std::vector<
+    GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
+    periodic_faces;
+  for (unsigned int d = 0; d < dim; ++d)
+    GridTools::collect_periodic_faces(
+      tria, 2 * d, 2 * d + 1, d, periodic_faces);
+  tria.add_periodicity(periodic_faces);
+
+  tria.refine_global(3 - dim);
+
+  FE_DGQHermite<dim>        fe(fe_degree);
+  DoFHandler<dim>           dof(tria);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  for (unsigned int test = 0; test < 2; ++test)
+    {
+      tria.refine_global(1);
+
+      dof.distribute_dofs(fe);
+
+      deallog << "Testing " << dof.get_fe().get_name();
+      deallog << " on " << dof.n_dofs() << " DoFs";
+      deallog << std::endl;
+
+      MappingQGeneric<dim> mapping(dof.get_fe().degree + 1);
+
+      LinearAlgebra::distributed::BlockVector<double> in(1);
+      LinearAlgebra::distributed::BlockVector<double> out(1);
+      LinearAlgebra::distributed::BlockVector<double> out_dist(1);
+
+      MatrixFree<dim, double>                          mf_data;
+      const QGauss<1>                                  quad(fe_degree + 1);
+      typename MatrixFree<dim, double>::AdditionalData data;
+      data.tasks_parallel_scheme =
+        MatrixFree<dim, double>::AdditionalData::none;
+      data.mapping_update_flags_inner_faces =
+        (update_gradients | update_JxW_values);
+      data.mapping_update_flags_boundary_faces =
+        (update_gradients | update_JxW_values);
+      data.initialize_mapping = false;
+
+      mf_data.reinit(mapping, dof, constraints, quad, data);
+
+      std::vector<types::global_dof_index> renumbering;
+      mf_data.renumber_dofs(renumbering);
+      dof.renumber_dofs(renumbering);
+
+      data.initialize_mapping = true;
+      mf_data.reinit(mapping, dof, constraints, quad, data);
+
+      mf_data.initialize_dof_vector(in.block(0));
+      mf_data.initialize_dof_vector(out.block(0));
+      mf_data.initialize_dof_vector(out_dist.block(0));
+
+      // Set random seed for reproducibility
+      Testing::srand(42);
+      for (unsigned int i = 0; i < in.block(0).local_size(); ++i)
+        {
+          const double entry           = Testing::rand() / (double)RAND_MAX;
+          in.block(0).local_element(i) = entry;
+        }
+
+      MatrixFreeTest<dim,
+                     fe_degree,
+                     fe_degree + 1,
+                     double,
+                     LinearAlgebra::distributed::BlockVector<double>>
+        mf(mf_data);
+      mf.vmult(out, in);
+
+      MatrixFreeVariant<dim,
+                        fe_degree,
+                        fe_degree + 1,
+                        double,
+                        LinearAlgebra::distributed::BlockVector<double>>
+        mf2(mf_data);
+      mf2.vmult(out_dist, in);
+
+      out_dist -= out;
+
+      double diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+      deallog << "Norm of difference:          " << diff_norm << " ";
+
+      // test again, now doing matrix-vector product twice
+      mf2.vmult(out_dist, in);
+      mf2.vmult(out_dist, in);
+      out_dist -= out;
+      diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+      deallog << diff_norm << std::endl;
+    }
+}
diff --git a/tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..bf136d0
--- /dev/null
@@ -0,0 +1,17 @@
+
+DEAL:2d::Testing FE_DGQHermite<2>(2) on 144 DoFs
+DEAL:2d::Norm of difference:          0 0
+DEAL:2d::Testing FE_DGQHermite<2>(2) on 576 DoFs
+DEAL:2d::Norm of difference:          0 0
+DEAL:2d::Testing FE_DGQHermite<2>(3) on 256 DoFs
+DEAL:2d::Norm of difference:          0 0
+DEAL:2d::Testing FE_DGQHermite<2>(3) on 1024 DoFs
+DEAL:2d::Norm of difference:          0 0
+DEAL:3d::Testing FE_DGQHermite<3>(2) on 216 DoFs
+DEAL:3d::Norm of difference:          0 0
+DEAL:3d::Testing FE_DGQHermite<3>(2) on 1728 DoFs
+DEAL:3d::Norm of difference:          0 0
+DEAL:3d::Testing FE_DGQHermite<3>(3) on 512 DoFs
+DEAL:3d::Norm of difference:          0 0
+DEAL:3d::Testing FE_DGQHermite<3>(3) on 4096 DoFs
+DEAL:3d::Norm of difference:          0 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.