From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 7 Aug 2022 18:19:27 +0000 (+0200)
Subject: Add second version to initialize MFF::ConstraintInfo
X-Git-Tag: v9.5.0-rc1~1048^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bd8ad048143e53d5ea008752a463637bcd36d6a1;p=dealii.git

Add second version to initialize MFF::ConstraintInfo
---

diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h
index fdf56f4429..fe41903b14 100644
--- a/include/deal.II/matrix_free/constraint_info.h
+++ b/include/deal.II/matrix_free/constraint_info.h
@@ -41,6 +41,10 @@ namespace internal
     class ConstraintInfo
     {
     public:
+      /**
+       * Version 1: indices are extracted from DoFCellAccessor and
+       * constraints are resolved with the help of AffineConstraints.
+       */
       void
       reinit(const DoFHandler<dim> &dof_handler,
              const unsigned int     n_cells,
@@ -55,6 +59,18 @@ namespace internal
           &                                                       constraints,
         const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
 
+      /**
+       * Version 2: no constraints, indices are user-provided.
+       */
+      void
+      reinit(const unsigned int n_cells);
+
+      void
+      read_dof_indices(
+        const unsigned int                                        cell_no,
+        const std::vector<types::global_dof_index> &              dof_indices,
+        const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
+
       void
       finalize();
 
@@ -164,6 +180,17 @@ namespace internal
 
 
 
+    template <int dim, typename Number>
+    inline void
+    ConstraintInfo<dim, Number>::reinit(const unsigned int n_cells)
+    {
+      this->dof_indices_per_cell.resize(n_cells);
+      this->plain_dof_indices_per_cell.resize(0);
+      this->constraint_indicator_per_cell.resize(n_cells);
+    }
+
+
+
     template <int dim, typename Number>
     inline void
     ConstraintInfo<dim, Number>::read_dof_indices(
@@ -291,6 +318,57 @@ namespace internal
 
 
 
+    template <int dim, typename Number>
+    inline void
+    ConstraintInfo<dim, Number>::read_dof_indices(
+      const unsigned int                          cell_no,
+      const std::vector<types::global_dof_index> &local_dof_indices_lex,
+      const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
+    {
+      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;
+      };
+
+      std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
+
+      auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no];
+      auto &dof_indices          = this->dof_indices_per_cell[cell_no];
+
+      for (const auto current_dof : local_dof_indices_lex)
+        {
+          // dof is constrained
+          if (current_dof == numbers::invalid_dof_index)
+            {
+              const std::vector<
+                std::pair<types::global_dof_index, typename Number::value_type>>
+                entries = {};
+
+              constraint_indicator.push_back(constraint_iterator);
+              constraint_indicator.back().second =
+                constraint_values.insert_entries(entries);
+
+              constraint_iterator.first = 0;
+            }
+          else
+            {
+              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()
@@ -302,8 +380,11 @@ namespace internal
       this->row_starts = {};
       this->row_starts.emplace_back(0, 0);
 
-      this->row_starts_plain_indices = {};
-      this->row_starts_plain_indices.emplace_back(0);
+      if (plain_dof_indices_per_cell.empty() == false)
+        {
+          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)
         {
@@ -318,12 +399,16 @@ namespace internal
           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());
+          if (plain_dof_indices_per_cell.empty() == false)
+            {
+              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());
+              this->row_starts_plain_indices.emplace_back(
+                this->plain_dof_indices.size());
+            }
         }
 
       std::vector<const std::vector<double> *> constraints(