]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce internal::MatrixFreeFunctions::ConstraintInfo 13519/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 9 Mar 2022 14:50:11 +0000 (15:50 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 10 Mar 2022 06:32:06 +0000 (07:32 +0100)
Add dummy test

include/deal.II/matrix_free/constraint_info.h [new file with mode: 0644]
tests/matrix_free/constraint_info_01.cc [new file with mode: 0644]
tests/matrix_free/constraint_info_01.output [new file with mode: 0644]

diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h
new file mode 100644 (file)
index 0000000..dbb635a
--- /dev/null
@@ -0,0 +1,530 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_constraint_info_h
+#define dealii_matrix_free_constraint_info_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/matrix_free/dof_info.h>
+#include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/hanging_nodes_internal.h>
+#include <deal.II/matrix_free/mapping_info.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    /**
+     * A helper class to apply constraints in matrix-free loops in
+     * user code. It combines constraint related functionalties from
+     * MatrixFree and FEEvaluation.
+     */
+    template <int dim, typename Number>
+    class ConstraintInfo
+    {
+    public:
+      void
+      reinit(const DoFHandler<dim> &dof_handler,
+             const unsigned int     n_cells,
+             const bool             use_fast_hanging_node_algorithm = true);
+
+      void
+      read_dof_indices(
+        const unsigned int                                    cell_no,
+        const unsigned int                                    mg_level,
+        const TriaIterator<DoFCellAccessor<dim, dim, false>> &cell,
+        const dealii::AffineConstraints<typename Number::value_type>
+          &                                                       constraints,
+        const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
+
+      void
+      finalize();
+
+      template <typename T, typename VectorType>
+      void
+      read_write_operation(const T &              operation,
+                           VectorType &           global_vector,
+                           AlignedVector<Number> &local_vector,
+                           const unsigned int     first_cell,
+                           const unsigned int     n_cells,
+                           const unsigned int     n_dofs_per_cell,
+                           const bool             apply_constraints) const;
+
+      void
+      apply_hanging_node_constraints(
+        const unsigned int     first_cell,
+        const unsigned int     n_lanes_filled,
+        const bool             transpose,
+        AlignedVector<Number> &evaluation_data_coarse) const;
+
+    private:
+      // for setup
+      ConstraintValues<double>               constraint_values;
+      std::vector<std::vector<unsigned int>> dof_indices_per_cell;
+      std::vector<std::vector<unsigned int>> plain_dof_indices_per_cell;
+      std::vector<std::vector<std::pair<unsigned short, unsigned short>>>
+        constraint_indicator_per_cell;
+
+      std::unique_ptr<HangingNodes<dim>>     hanging_nodes;
+      std::vector<std::vector<unsigned int>> lexicographic_numbering;
+
+    public:
+      // for read_write_operation()
+      std::vector<unsigned int> dof_indices;
+      std::vector<std::pair<unsigned short, unsigned short>>
+                                                         constraint_indicator;
+      std::vector<std::pair<unsigned int, unsigned int>> row_starts;
+
+      std::vector<unsigned int> plain_dof_indices;
+      std::vector<unsigned int> row_starts_plain_indices;
+
+      // for constraint_pool_begin/end()
+      std::vector<typename Number::value_type> constraint_pool_data;
+      std::vector<unsigned int>                constraint_pool_row_index;
+
+      std::vector<ShapeInfo<Number>> shape_infos;
+      std::vector<ConstraintKinds>   hanging_node_constraint_masks;
+      std::vector<unsigned int>      active_fe_indices;
+
+    private:
+      inline const typename Number::value_type *
+      constraint_pool_begin(const unsigned int row) const;
+
+      inline const typename Number::value_type *
+      constraint_pool_end(const unsigned int row) const;
+    };
+
+
+
+    template <int dim, typename Number>
+    inline void
+    ConstraintInfo<dim, Number>::reinit(
+      const DoFHandler<dim> &dof_handler,
+      const unsigned int     n_cells,
+      const bool             use_fast_hanging_node_algorithm)
+    {
+      this->dof_indices_per_cell.resize(n_cells);
+      this->plain_dof_indices_per_cell.resize(n_cells);
+      this->constraint_indicator_per_cell.resize(n_cells);
+
+      if (use_fast_hanging_node_algorithm &&
+          dof_handler.get_triangulation().has_hanging_nodes())
+        {
+          hanging_nodes = std::make_unique<HangingNodes<dim>>(
+            dof_handler.get_triangulation());
+
+          hanging_node_constraint_masks.resize(n_cells);
+        }
+
+      const auto &fes = dof_handler.get_fe_collection();
+      lexicographic_numbering.resize(fes.size());
+      shape_infos.resize(fes.size());
+
+      for (unsigned int i = 0; i < fes.size(); ++i)
+        {
+          if (fes[i].reference_cell().is_hyper_cube())
+            {
+              const Quadrature<1> dummy_quadrature(
+                std::vector<Point<1>>(1, Point<1>()));
+              shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
+            }
+          else
+            {
+              const auto dummy_quadrature =
+                fes[i].reference_cell().template get_gauss_type_quadrature<dim>(
+                  1);
+              shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
+            }
+
+          lexicographic_numbering[i] = shape_infos[i].lexicographic_numbering;
+        }
+      active_fe_indices.resize(n_cells);
+    }
+
+
+
+    template <int dim, typename Number>
+    inline void
+    ConstraintInfo<dim, Number>::read_dof_indices(
+      const unsigned int                                            cell_no,
+      const unsigned int                                            mg_level,
+      const TriaIterator<DoFCellAccessor<dim, dim, false>> &        cell,
+      const dealii::AffineConstraints<typename Number::value_type> &constraints,
+      const std::shared_ptr<const Utilities::MPI::Partitioner> &    partitioner)
+    {
+      std::vector<types::global_dof_index> local_dof_indices(
+        cell->get_fe().n_dofs_per_cell());
+      std::vector<types::global_dof_index> local_dof_indices_lex(
+        cell->get_fe().n_dofs_per_cell());
+
+      if (mg_level == numbers::invalid_unsigned_int)
+        cell->get_dof_indices(local_dof_indices);
+      else
+        cell->get_mg_dof_indices(local_dof_indices);
+
+      {
+        AssertIndexRange(cell->active_fe_index(), shape_infos.size());
+
+        const auto &lexicographic_numbering =
+          shape_infos[cell->active_fe_index()].lexicographic_numbering;
+
+        AssertDimension(lexicographic_numbering.size(),
+                        local_dof_indices.size());
+
+        for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
+          local_dof_indices_lex[i] =
+            local_dof_indices[lexicographic_numbering[i]];
+      }
+
+      std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
+
+      AssertIndexRange(cell_no, this->constraint_indicator_per_cell.size());
+      AssertIndexRange(cell_no, this->dof_indices_per_cell.size());
+      AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size());
+      AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size());
+
+      auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no];
+      auto &dof_indices          = this->dof_indices_per_cell[cell_no];
+      auto &plain_dof_indices    = this->plain_dof_indices_per_cell[cell_no];
+
+      AssertDimension(constraint_indicator_per_cell[cell_no].size(), 0);
+      AssertDimension(dof_indices_per_cell[cell_no].size(), 0);
+      AssertDimension(plain_dof_indices_per_cell[cell_no].size(), 0);
+
+      const auto global_to_local =
+        [&](const types::global_dof_index global_index) -> unsigned int {
+        if (partitioner)
+          return partitioner->global_to_local(global_index);
+        else
+          return global_index;
+      };
+
+      // plain indices
+      plain_dof_indices.resize(local_dof_indices_lex.size());
+      for (unsigned int i = 0; i < local_dof_indices_lex.size(); ++i)
+        plain_dof_indices[i] = global_to_local(local_dof_indices_lex[i]);
+
+      if (hanging_nodes)
+        {
+          AssertIndexRange(cell_no, this->hanging_node_constraint_masks.size());
+          AssertIndexRange(cell_no, this->active_fe_indices.size());
+
+          std::vector<ConstraintKinds> mask(cell->get_fe().n_components());
+          hanging_nodes->setup_constraints(
+            cell, {}, lexicographic_numbering, local_dof_indices_lex, mask);
+
+          hanging_node_constraint_masks[cell_no] = mask[0];
+          active_fe_indices[cell_no]             = cell->active_fe_index();
+        }
+
+      for (auto current_dof : local_dof_indices_lex)
+        {
+          const auto *entries_ptr =
+            constraints.get_constraint_entries(current_dof);
+
+          // dof is constrained
+          if (entries_ptr != nullptr)
+            {
+              const auto &                  entries   = *entries_ptr;
+              const types::global_dof_index n_entries = entries.size();
+              if (n_entries == 1 &&
+                  std::abs(entries[0].second - 1.) <
+                    100 * std::numeric_limits<double>::epsilon())
+                {
+                  current_dof = entries[0].first;
+                  goto no_constraint;
+                }
+
+              constraint_indicator.push_back(constraint_iterator);
+              constraint_indicator.back().second =
+                constraint_values.insert_entries(entries);
+
+              // reset constraint iterator for next round
+              constraint_iterator.first = 0;
+
+              if (n_entries > 0)
+                {
+                  const std::vector<types::global_dof_index>
+                    &constraint_indices = constraint_values.constraint_indices;
+                  for (unsigned int j = 0; j < n_entries; ++j)
+                    {
+                      dof_indices.push_back(
+                        global_to_local(constraint_indices[j]));
+                    }
+                }
+            }
+          else
+            {
+            no_constraint:
+              dof_indices.push_back(global_to_local(current_dof));
+
+              // make sure constraint_iterator.first is always within the
+              // bounds of unsigned short
+              Assert(constraint_iterator.first <
+                       (1 << (8 * sizeof(unsigned short))) - 1,
+                     ExcInternalError());
+              constraint_iterator.first++;
+            }
+        }
+    }
+
+
+
+    template <int dim, typename Number>
+    inline void
+    ConstraintInfo<dim, Number>::finalize()
+    {
+      this->dof_indices          = {};
+      this->plain_dof_indices    = {};
+      this->constraint_indicator = {};
+
+      this->row_starts = {};
+      this->row_starts.emplace_back(0, 0);
+
+      this->row_starts_plain_indices = {};
+      this->row_starts_plain_indices.emplace_back(0);
+
+      for (unsigned int i = 0; i < dof_indices_per_cell.size(); ++i)
+        {
+          this->dof_indices.insert(this->dof_indices.end(),
+                                   dof_indices_per_cell[i].begin(),
+                                   dof_indices_per_cell[i].end());
+          this->constraint_indicator.insert(
+            this->constraint_indicator.end(),
+            constraint_indicator_per_cell[i].begin(),
+            constraint_indicator_per_cell[i].end());
+
+          this->row_starts.emplace_back(this->dof_indices.size(),
+                                        this->constraint_indicator.size());
+
+          this->plain_dof_indices.insert(this->plain_dof_indices.end(),
+                                         plain_dof_indices_per_cell[i].begin(),
+                                         plain_dof_indices_per_cell[i].end());
+
+          this->row_starts_plain_indices.emplace_back(
+            this->plain_dof_indices.size());
+        }
+
+      std::vector<const std::vector<double> *> constraints(
+        constraint_values.constraints.size());
+      unsigned int length = 0;
+      for (const auto &it : constraint_values.constraints)
+        {
+          AssertIndexRange(it.second, constraints.size());
+          constraints[it.second] = &it.first;
+          length += it.first.size();
+        }
+
+      constraint_pool_data.clear();
+      constraint_pool_data.reserve(length);
+      constraint_pool_row_index.reserve(constraint_values.constraints.size() +
+                                        1);
+      constraint_pool_row_index.resize(1, 0);
+
+      for (const auto &constraint : constraints)
+        {
+          Assert(constraint != nullptr, ExcInternalError());
+          constraint_pool_data.insert(constraint_pool_data.end(),
+                                      constraint->begin(),
+                                      constraint->end());
+          constraint_pool_row_index.push_back(constraint_pool_data.size());
+        }
+
+      AssertDimension(constraint_pool_data.size(), length);
+
+      dof_indices_per_cell.clear();
+      constraint_indicator_per_cell.clear();
+
+      if (hanging_nodes && std::all_of(hanging_node_constraint_masks.begin(),
+                                       hanging_node_constraint_masks.end(),
+                                       [](const auto i) {
+                                         return i ==
+                                                ConstraintKinds::unconstrained;
+                                       }))
+        hanging_node_constraint_masks.clear();
+    }
+
+
+
+    template <int dim, typename Number>
+    template <typename T, typename VectorType>
+    inline void
+    ConstraintInfo<dim, Number>::read_write_operation(
+      const T &              operation,
+      VectorType &           global_vector,
+      AlignedVector<Number> &local_vector,
+      const unsigned int     first_cell,
+      const unsigned int     n_cells,
+      const unsigned int     n_dofs_per_cell,
+      const bool             apply_constraints) const
+    {
+      if (apply_constraints == false)
+        {
+          for (unsigned int v = 0; v < n_cells; ++v)
+            {
+              const unsigned int  cell_index = first_cell + v;
+              const unsigned int *dof_indices =
+                this->plain_dof_indices.data() +
+                this->row_starts_plain_indices[cell_index];
+
+              for (unsigned int i = 0; i < n_dofs_per_cell; ++dof_indices, ++i)
+                operation.process_dof(*dof_indices,
+                                      global_vector,
+                                      local_vector[i][v]);
+            }
+
+          return;
+        }
+
+      for (unsigned int v = 0; v < n_cells; ++v)
+        {
+          const unsigned int  cell_index = first_cell + v;
+          const unsigned int *dof_indices =
+            this->dof_indices.data() + this->row_starts[cell_index].first;
+          unsigned int index_indicators = this->row_starts[cell_index].second;
+          unsigned int next_index_indicators =
+            this->row_starts[cell_index + 1].second;
+
+          unsigned int ind_local = 0;
+          for (; index_indicators != next_index_indicators; ++index_indicators)
+            {
+              const std::pair<unsigned short, unsigned short> indicator =
+                this->constraint_indicator[index_indicators];
+
+              // run through values up to next constraint
+              for (unsigned int j = 0; j < indicator.first; ++j)
+                operation.process_dof(dof_indices[j],
+                                      global_vector,
+                                      local_vector[ind_local + j][v]);
+
+              ind_local += indicator.first;
+              dof_indices += indicator.first;
+
+              // constrained case: build the local value as a linear
+              // combination of the global value according to constraints
+              typename Number::value_type value;
+              operation.pre_constraints(local_vector[ind_local][v], value);
+
+              const typename Number::value_type *data_val =
+                this->constraint_pool_begin(indicator.second);
+              const typename Number::value_type *end_pool =
+                this->constraint_pool_end(indicator.second);
+              for (; data_val != end_pool; ++data_val, ++dof_indices)
+                operation.process_constraint(*dof_indices,
+                                             *data_val,
+                                             global_vector,
+                                             value);
+
+              operation.post_constraints(value, local_vector[ind_local][v]);
+              ind_local++;
+            }
+
+          AssertIndexRange(ind_local, n_dofs_per_cell + 1);
+
+          for (; ind_local < n_dofs_per_cell; ++dof_indices, ++ind_local)
+            operation.process_dof(*dof_indices,
+                                  global_vector,
+                                  local_vector[ind_local][v]);
+        }
+    }
+
+
+
+    template <int dim, typename Number>
+    inline void
+    ConstraintInfo<dim, Number>::apply_hanging_node_constraints(
+      const unsigned int     first_cell,
+      const unsigned int     n_lanes_filled,
+      const bool             transpose,
+      AlignedVector<Number> &evaluation_data_coarse) const
+    {
+      if (hanging_node_constraint_masks.size() == 0)
+        return;
+
+      std::array<ConstraintKinds, Number::size()> constraint_mask;
+
+      bool hn_available = false;
+
+      for (unsigned int v = 0; v < n_lanes_filled; ++v)
+        {
+          const auto mask = hanging_node_constraint_masks[first_cell + v];
+
+          constraint_mask[v] = mask;
+
+          hn_available |= (mask != ConstraintKinds::unconstrained);
+        }
+
+      if (hn_available == true)
+        {
+          for (unsigned int v = n_lanes_filled; v < Number::size(); ++v)
+            constraint_mask[v] = ConstraintKinds::unconstrained;
+
+          for (unsigned int i = 1; i < n_lanes_filled; ++i)
+            AssertDimension(active_fe_indices[first_cell],
+                            active_fe_indices[first_cell + i]);
+
+          const auto &shape_info = shape_infos[active_fe_indices[first_cell]];
+
+          dealii::internal::FEEvaluationHangingNodesFactory<
+            dim,
+            typename Number::value_type,
+            Number>::apply(shape_info.n_components,
+                           shape_info.data.front().fe_degree,
+                           shape_info,
+                           transpose,
+                           constraint_mask,
+                           evaluation_data_coarse.begin());
+        }
+    }
+
+
+
+    template <int dim, typename Number>
+    inline const typename Number::value_type *
+    ConstraintInfo<dim, Number>::constraint_pool_begin(
+      const unsigned int row) const
+    {
+      AssertIndexRange(row, constraint_pool_row_index.size() - 1);
+      return constraint_pool_data.empty() ?
+               nullptr :
+               constraint_pool_data.data() + constraint_pool_row_index[row];
+    }
+
+
+
+    template <int dim, typename Number>
+    inline const typename Number::value_type *
+    ConstraintInfo<dim, Number>::constraint_pool_end(
+      const unsigned int row) const
+    {
+      AssertIndexRange(row, constraint_pool_row_index.size() - 1);
+      return constraint_pool_data.empty() ?
+               nullptr :
+               constraint_pool_data.data() + constraint_pool_row_index[row + 1];
+    }
+
+
+  } // namespace MatrixFreeFunctions
+} // namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/matrix_free/constraint_info_01.cc b/tests/matrix_free/constraint_info_01.cc
new file mode 100644 (file)
index 0000000..b283c0d
--- /dev/null
@@ -0,0 +1,223 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test internal::MatrixFreeFunctions::ConstraintInfo and compare the
+// result with FEEvaluation.
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/constraint_info.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+int
+main()
+{
+  initlog();
+
+  const int dim    = 2;
+  const int degree = 1;
+
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+  using VectorType          = LinearAlgebra::distributed::Vector<Number>;
+
+  Triangulation<dim> tria;
+  GridGenerator::subdivided_hyper_cube(tria, 2);
+  tria.begin()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  QGauss<dim>    quad(degree + 1);
+  FE_Q<dim>      fe(degree);
+  MappingQ1<dim> mapping;
+
+  DoFHandler<dim> dof_handler;
+  dof_handler.reinit(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<Number> constraints;
+
+  DoFTools::make_zero_boundary_constraints(dof_handler, constraints);
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  constraints.close();
+
+  typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+    additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients |
+                                         update_JxW_values |
+                                         dealii::update_quadrature_points;
+
+  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraints, quad, additional_data);
+
+  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArrayType>
+    constraint_info;
+  constraint_info.reinit(dof_handler, matrix_free.n_physical_cells(), true);
+
+  for (unsigned int cell = 0, cell_no = 0; cell < matrix_free.n_cell_batches();
+       ++cell)
+    for (unsigned int v = 0;
+         v < matrix_free.n_active_entries_per_cell_batch(cell);
+         ++v, ++cell_no)
+      constraint_info.read_dof_indices(cell_no,
+                                       -1,
+                                       matrix_free.get_cell_iterator(cell, v),
+                                       constraints,
+                                       matrix_free.get_vector_partitioner());
+
+  constraint_info.finalize();
+
+  const auto print_stat = [](const auto &dof_info) {
+    for (const auto i : dof_info.dof_indices)
+      deallog << i << " ";
+    deallog << std::endl;
+
+    for (const auto i : dof_info.constraint_indicator)
+      deallog << "(" << i.first << "," << i.second << ") ";
+    deallog << std::endl;
+
+    for (const auto i : dof_info.row_starts)
+      deallog << "(" << i.first << "," << i.second << ") ";
+    deallog << std::endl;
+
+    deallog << std::endl;
+  };
+
+  print_stat(constraint_info);
+  print_stat(matrix_free.get_dof_info());
+
+  VectorType src, dst, dst_mf;
+
+  const auto initialize_dof_vectors = [&]() {
+    matrix_free.initialize_dof_vector(src);
+    matrix_free.initialize_dof_vector(dst);
+    matrix_free.initialize_dof_vector(dst_mf);
+
+    for (const auto i : src.locally_owned_elements())
+      src[i] = i;
+  };
+
+  initialize_dof_vectors();
+
+  AlignedVector<VectorizedArrayType> local_vector(fe.n_dofs_per_cell());
+
+  const auto read_dof_values = [&](const unsigned int first_cell,
+                                   const unsigned int n_cells,
+                                   const unsigned int n_dofs_per_cell) {
+    internal::VectorReader<Number, VectorizedArrayType> reader;
+    constraint_info.read_write_operation(
+      reader, src, local_vector, first_cell, n_cells, n_dofs_per_cell, true);
+    constraint_info.apply_hanging_node_constraints(first_cell,
+                                                   n_cells,
+                                                   false,
+                                                   local_vector);
+  };
+
+  const auto distribute_local_to_global =
+    [&](const unsigned int first_cell,
+        const unsigned int n_cells,
+        const unsigned int n_dofs_per_cell) {
+      internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
+        writer;
+      constraint_info.apply_hanging_node_constraints(first_cell,
+                                                     n_cells,
+                                                     true,
+                                                     local_vector);
+      constraint_info.read_write_operation(
+        writer, dst, local_vector, first_cell, n_cells, n_dofs_per_cell, true);
+    };
+
+  const auto read_dof_values_plain = [&](const unsigned int first_cell,
+                                         const unsigned int n_cells,
+                                         const unsigned int n_dofs_per_cell) {
+    internal::VectorReader<Number, VectorizedArrayType> reader;
+    constraint_info.read_write_operation(
+      reader, src, local_vector, first_cell, n_cells, n_dofs_per_cell, false);
+  };
+
+  const auto set_dof_values_plain = [&](const unsigned int first_cell,
+                                        const unsigned int n_cells,
+                                        const unsigned int n_dofs_per_cell) {
+    internal::VectorSetter<Number, VectorizedArrayType> writer;
+    constraint_info.read_write_operation(
+      writer, dst, local_vector, first_cell, n_cells, n_dofs_per_cell, false);
+  };
+
+  unsigned int       first_cell      = 0;
+  const unsigned int n_dofs_per_cell = fe.n_dofs_per_cell();
+
+  FEEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType> fe_eval(matrix_free);
+
+  for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell)
+    {
+      {
+        const unsigned int n_cells =
+          matrix_free.n_active_entries_per_cell_batch(cell);
+
+        read_dof_values(first_cell, n_cells, n_dofs_per_cell);
+        distribute_local_to_global(first_cell, n_cells, n_dofs_per_cell);
+
+        first_cell += n_cells;
+      }
+
+      {
+        fe_eval.reinit(cell);
+
+        fe_eval.read_dof_values(src);
+        fe_eval.distribute_local_to_global(dst_mf);
+      }
+    }
+
+  dst.print(deallog.get_file_stream());
+  dst_mf.print(deallog.get_file_stream());
+  deallog << std::endl;
+
+  initialize_dof_vectors();
+  first_cell = 0;
+
+  for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell)
+    {
+      {
+        const unsigned int n_cells =
+          matrix_free.n_active_entries_per_cell_batch(cell);
+
+        read_dof_values_plain(first_cell, n_cells, n_dofs_per_cell);
+        set_dof_values_plain(first_cell, n_cells, n_dofs_per_cell);
+
+        first_cell += n_cells;
+      }
+
+      {
+        fe_eval.reinit(cell);
+
+        fe_eval.read_dof_values_plain(src);
+        fe_eval.set_dof_values_plain(dst_mf);
+      }
+    }
+
+  dst.print(deallog.get_file_stream());
+  dst_mf.print(deallog.get_file_stream());
+}
diff --git a/tests/matrix_free/constraint_info_01.output b/tests/matrix_free/constraint_info_01.output
new file mode 100644 (file)
index 0000000..3658842
--- /dev/null
@@ -0,0 +1,26 @@
+
+DEAL::11 11 2 11 2 11 2 2 2 2 
+DEAL::(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (1,0) (1,0) (0,0) (0,0) (0,0) (1,0) (0,0) (1,0) (0,0) (1,0) (0,0) (0,0) 
+DEAL::(0,0) (1,3) (3,5) (5,7) (7,9) (8,12) (9,15) (10,18) 
+DEAL::
+DEAL::11 11 2 11 2 11 2 2 2 2 
+DEAL::(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (1,0) (1,0) (0,0) (0,0) (0,0) (1,0) (0,0) (1,0) (0,0) (1,0) (0,0) (0,0) 
+DEAL::(0,0) (1,3) (3,5) (5,7) (7,9) (8,12) (9,15) (10,18) (10,18) 
+DEAL::
+Process #0
+Local range: [0, 14), global size: 14
+Vector data:
+0.000e+00 0.000e+00 1.000e+01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.400e+01 0.000e+00 0.000e+00 
+Process #0
+Local range: [0, 14), global size: 14
+Vector data:
+0.000e+00 0.000e+00 1.000e+01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.400e+01 0.000e+00 0.000e+00 
+DEAL::
+Process #0
+Local range: [0, 14), global size: 14
+Vector data:
+0.000e+00 1.000e+00 2.000e+00 3.000e+00 4.000e+00 5.000e+00 6.000e+00 7.000e+00 8.000e+00 9.000e+00 1.000e+01 1.100e+01 1.200e+01 1.300e+01 
+Process #0
+Local range: [0, 14), global size: 14
+Vector data:
+0.000e+00 1.000e+00 2.000e+00 3.000e+00 4.000e+00 5.000e+00 6.000e+00 7.000e+00 8.000e+00 9.000e+00 1.000e+01 1.100e+01 1.200e+01 1.300e+01 

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.