Implement get_affine_constraints()

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index d8e519ab6d3d3e11788cbc6e97fe91c8ba891fcc..b98f7e32ac72a9316702bb030b5a802191f43206 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1630,6 +1630,15 @@ public:
   const DoFHandler<dim> &
   get_dof_handler(const unsigned int dof_handler_index = 0) const;
 
+  /**
+   * Return the AffineConstraints with the index as given to the
+   * respective `std::vector` argument in the reinit() function. Only available
+   * if the AffineConstraints objects have the same template parameter Number as
+   * MatrixFree. Throws an exception otherwise.
+   */
+  const AffineConstraints<Number> &
+  get_affine_constraints(const unsigned int dof_handler_index = 0) const;
+
   /**
    * Return the cell iterator in deal.II speak to a given cell batch
    * (populating several lanes in a VectorizedArray) and the lane index within
@@ -2002,6 +2011,14 @@ private:
    */
   std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
 
+  /**
+   * Pointers to the AffineConstraints underlying the current problem. Only
+   * filled with an AffineConstraints object if objects of the same `Number`
+   * template parameter as the `Number` template of MatrixFree is passed to
+   * reinit(). Filled with nullptr otherwise.
+   */
+  std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
+
   /**
    * Contains the information about degrees of freedom on the individual cells
    * and constraints.

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 90417b5f0616753f9b03c851a760d89c8b558ca3..9919a0463e00c015c73ad5a9c6833481fd05f960 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -228,6 +228,23 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_dof_handler(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+const AffineConstraints<Number> &
+MatrixFree<dim, Number, VectorizedArrayType>::get_affine_constraints(
+  const unsigned int dof_handler_index) const
+{
+  AssertIndexRange(dof_handler_index, n_components());
+
+  AssertThrow(affine_constraints[dof_handler_index] != nullptr,
+              ExcMessage(
+                "AffineConstraints are only accessible if its Number type "
+                "matches with Number type of MatrixFree!"));
+
+  return *(affine_constraints[dof_handler_index]);
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 typename DoFHandler<dim>::cell_iterator
 MatrixFree<dim, Number, VectorizedArrayType>::get_cell_iterator(
@@ -330,6 +347,29 @@ MatrixFree<dim, Number, VectorizedArrayType>::copy_from(
 
 
 
+namespace internal
+{
+  template <typename Number, typename Number2>
+  void
+  store_affine_constraints(
+    const dealii::AffineConstraints<Number2> *,
+    SmartPointer<const dealii::AffineConstraints<Number>> &stored_constraints)
+  {
+    stored_constraints = nullptr;
+  }
+
+  template <typename Number>
+  void
+  store_affine_constraints(
+    const dealii::AffineConstraints<Number> *              affine_constraints,
+    SmartPointer<const dealii::AffineConstraints<Number>> &stored_constraints)
+  {
+    stored_constraints = affine_constraints;
+  }
+} // namespace internal
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 template <typename number2, int q_dim>
 void
@@ -374,6 +414,11 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
                 .reinit(quad[nq][q_no], dof_handler[no]->get_fe(fe_no), b);
   }
 
+  // Store pointers to AffineConstraints objects if Number type matches
+  affine_constraints.resize(constraints.size());
+  for (unsigned int no = 0; no < constraints.size(); ++no)
+    internal::store_affine_constraints(constraints[no], affine_constraints[no]);
+
   if (additional_data.initialize_indices == true)
     {
       clear();

diff --git a/tests/matrix_free/affine_constraints_store_and_get.cc b/tests/matrix_free/affine_constraints_store_and_get.cc
new file mode 100644 (file)
index 0000000..c0d96d5
--- /dev/null
+++ b/tests/matrix_free/affine_constraints_store_and_get.cc
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test internal::store_affine_constraints() and get_affine_constraints()
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+
+int
+main()
+{
+  initlog();
+
+  const unsigned int fe_degree = 1;
+  constexpr int      dim       = 2;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<1> quadrature(fe_degree + 1);
+
+  AffineConstraints<double> double_ac;
+  double_ac.close();
+
+  {
+    typename MatrixFree<dim, double>::AdditionalData additional_data_double;
+    MatrixFree<dim, double>                          matrix_free_double;
+    matrix_free_double.reinit(MappingQ1<dim>{},
+                              dof_handler,
+                              double_ac,
+                              quadrature,
+                              additional_data_double);
+
+    if (&matrix_free_double.get_affine_constraints() == &double_ac)
+      deallog << "OK" << std::endl;
+  }
+
+  {
+    typename MatrixFree<dim, float>::AdditionalData additional_data_float;
+    MatrixFree<dim, float>                          matrix_free_float;
+    matrix_free_float.reinit(MappingQ1<dim>{},
+                             dof_handler,
+                             double_ac,
+                             quadrature,
+                             additional_data_float);
+
+    try
+      {
+        matrix_free_float.get_affine_constraints();
+      }
+    catch (...)
+      {
+        deallog << "OK" << std::endl;
+      }
+  }
+}

diff --git a/tests/matrix_free/affine_constraints_store_and_get.output b/tests/matrix_free/affine_constraints_store_and_get.output
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
+++ b/tests/matrix_free/affine_constraints_store_and_get.output
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK