implement dof fe coupling support in make_sparsity_pattern

include/deal.II/fe/fe.h

@@ -1610,6 +1610,12 @@ class FiniteElement : public EnableObserverPointer,
   bool
   is_primitive(const unsigned int i) const;
 
+  /**
+   *
+   */
+  virtual const Table<2, bool> &
+  get_local_dof_sparsity_pattern() const;
+
   /**
    * Number of base elements in a mixed discretization.
    *
@@ -3345,6 +3351,16 @@ FiniteElement<dim, spacedim>::is_primitive(const unsigned int i) const
 
 
 
+template <int dim, int spacedim>
+inline const Table<2, bool> &
+FiniteElement<dim, spacedim>::get_local_dof_sparsity_pattern() const
+{
+  static Table<2, bool> return_value;
+  return return_value;
+}
+
+
+
 template <int dim, int spacedim>
 inline GeometryPrimitive
 FiniteElement<dim, spacedim>::get_associated_geometry_primitive(

include/deal.II/fe/fe_q_iso_q1.h

@@ -159,6 +159,13 @@ class FE_Q_iso_Q1 : public FE_Q_Base<dim, spacedim>
     const std::vector<Vector<double>> &support_point_values,
     std::vector<double>               &nodal_values) const override;
 
+  /**
+   * Override the function to specify the subset of dofs that are really
+   * coupling for each "sub" cell.
+   */
+  virtual const Table<2, bool> &
+  get_local_dof_sparsity_pattern() const override;
+
   /**
    * @name Functions to support hp
    * @{
@@ -172,10 +179,20 @@ class FE_Q_iso_Q1 : public FE_Q_Base<dim, spacedim>
                          const unsigned int codim = 0) const override final;
 
   /** @} */
+
+private:
+  Table<2, bool> local_dof_sparsity_pattern;
 };
 
 
 
+template <int dim, int spacedim>
+inline const Table<2, bool> &
+FE_Q_iso_Q1<dim, spacedim>::get_local_dof_sparsity_pattern() const
+{
+  return local_dof_sparsity_pattern;
+}
+
 /** @} */
 
 DEAL_II_NAMESPACE_CLOSE

source/dofs/dof_tools_sparsity.cc

@@ -82,6 +82,13 @@ namespace DoFTools
                  "locally owned one does not make sense."));
       }
 
+    const auto                 &fe_collection = dof.get_fe_collection();
+    std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
+    for (unsigned int f = 0; f < fe_collection.size(); ++f)
+      {
+        fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
+      }
+
     std::vector<types::global_dof_index> dofs_on_this_cell;
     dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
 
@@ -100,9 +107,16 @@ namespace DoFTools
           // make sparsity pattern for this cell. if no constraints pattern
           // was given, then the following call acts as if simply no
           // constraints existed
-          constraints.add_entries_local_to_global(dofs_on_this_cell,
-                                                  sparsity,
-                                                  keep_constrained_dofs);
+          const types::fe_index fe_index = cell->active_fe_index();
+          if (fe_dof_mask[fe_index].empty())
+            constraints.add_entries_local_to_global(dofs_on_this_cell,
+                                                    sparsity,
+                                                    keep_constrained_dofs);
+          else
+            constraints.add_entries_local_to_global(dofs_on_this_cell,
+                                                    sparsity,
+                                                    keep_constrained_dofs,
+                                                    fe_dof_mask[fe_index]);
         }
   }
 
@@ -152,6 +166,12 @@ namespace DoFTools
     const std::vector<Table<2, Coupling>> dof_mask //(fe_collection.size())
       = dof_couplings_from_component_couplings(fe_collection, couplings);
 
+    std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
+    for (unsigned int f = 0; f < fe_collection.size(); ++f)
+      {
+        fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
+      }
+
     // Convert the dof_mask to bool_dof_mask so we can pass it
     // to constraints.add_entries_local_to_global()
     std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
@@ -163,7 +183,8 @@ namespace DoFTools
         bool_dof_mask[f].fill(false);
         for (unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
           for (unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
-            if (dof_mask[f](i, j) != none)
+            if (dof_mask[f](i, j) != none &&
+                (fe_dof_mask[f].empty() || fe_dof_mask[f](i, j)))
               bool_dof_mask[f](i, j) = true;
       }
 

source/fe/fe_q_iso_q1.cc

@@ -14,10 +14,12 @@
 
 
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
 
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_tools.h>
 
 #include <deal.II/lac/vector.h>
 
@@ -50,6 +52,42 @@ FE_Q_iso_Q1<dim, spacedim>::FE_Q_iso_Q1(const unsigned int subdivisions)
   const QIterated<1>  points(trapez, subdivisions);
 
   this->initialize(points.get_points());
+
+  {
+    const std::vector<unsigned int> R =
+      FETools::lexicographic_to_hierarchic_numbering<dim>(subdivisions);
+
+    local_dof_sparsity_pattern.reinit(this->n_dofs_per_cell(),
+                                      this->n_dofs_per_cell());
+
+    {
+      local_dof_sparsity_pattern.fill(false);
+
+      const unsigned int N   = Utilities::pow(points.size(), dim);
+      const int          N1d = points.size();
+      for (unsigned int i = 0; i < N; ++i)
+        for (unsigned int j = 0; j < N; ++j)
+          {
+            // compute l1 distance:
+            int distance = 0;
+
+            int xi = i;
+            int xj = j;
+            for (unsigned int d = 0; d < dim; ++d)
+              {
+                int current_distance = std::abs((xi % N1d) - (xj % N1d));
+                xi /= N1d;
+                xj /= N1d;
+                distance = std::max(distance, current_distance);
+                if (distance > 1)
+                  break;
+              }
+
+            if (distance <= 1)
+              local_dof_sparsity_pattern(R[i], R[j]) = true;
+          }
+    }
+  }
 }
 
 
@@ -178,6 +216,7 @@ FE_Q_iso_Q1<dim, spacedim>::compare_for_domination(
 }
 
 
+
 // explicit instantiations
 #include "fe_q_iso_q1.inst"