support DoF coupling in FESystem

Author: tjhei

Diff for: source/fe/fe_system.cc

@@ -2009,6 +2009,69 @@ FESystem<dim, spacedim>::initialize(
       Assert(index == this->n_dofs_per_line(), ExcInternalError());
     });
 
+
+  // Compute local_dof_sparsity_pattern if any of our base elements contains a
+  // non-empty one (empty denotes the default of all DoFs coupling within a
+  // cell). Note the we currently only handle coupling within a base element and
+  // not between two different base elements. Handling the latter could be
+  // doable if the underlying element happens to be identical, but we currently
+  // have no functionality to compute the coupling between different elements
+  // with a pattern (for example FE_Q_iso_Q1 with different degrees).
+  {
+    // Does any of our base elements not couple all DoFs?
+    const bool have_nonempty = [&]() -> bool {
+      for (unsigned int b = 0; b < this->n_base_elements(); ++b)
+        {
+          if (!this->base_element(b).get_local_dof_sparsity_pattern().empty() &&
+              (this->element_multiplicity(b) > 0))
+            return true;
+        }
+      return false;
+    }();
+
+    if (have_nonempty)
+      {
+        this->local_dof_sparsity_pattern.reinit(this->n_dofs_per_cell(),
+                                                this->n_dofs_per_cell());
+        // by default, everything couples:
+        this->local_dof_sparsity_pattern.fill(true);
+
+        unsigned int current_dof_index = 0;
+
+        for (unsigned int b = 0; b < this->n_base_elements(); ++b)
+          {
+            const auto &pattern =
+              this->base_element(b).get_local_dof_sparsity_pattern();
+            const unsigned int m = this->element_multiplicity(b);
+
+            if (!pattern.empty())
+              {
+                // the base element b with multiplicity m has a sparsity pattern
+                // where not all DoFs couple. Within this base element, each of
+                // them couples to each other one only based on this pattern.
+                // This means we can copy m times m copies of the pattern into
+                // ours:
+                const unsigned int n_dofs = pattern.size(0);
+                for (unsigned int bi = 0; bi < m; ++bi)
+                  for (unsigned int bj = 0; bj < m; ++bj)
+                    for (unsigned int i = 0; i < n_dofs; ++i)
+                      for (unsigned int j = 0; j < n_dofs; ++j)
+                        this->local_dof_sparsity_pattern(
+                          current_dof_index + bi * n_dofs + i,
+                          current_dof_index + bj * n_dofs + j) = pattern(i, j);
+              }
+
+            current_dof_index += this->base_element(b).n_dofs_per_cell() * m;
+          }
+      }
+    else
+      {
+        // None of our base elements has a sparsity pattern set, meaning we
+        // don't have to do anything (the default empty
+        // local_dof_sparsity_pattern is correct).
+      }
+  }
+
   // wait for all of this to finish
   init_tasks.join_all();
 }

Diff for: tests/fe/local_sparsity_01.cc

@@ -0,0 +1,67 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2024 - 2025 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+// Check FiniteElement::get_local_dof_sparsity_pattern() for FESystem.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include "deal.II/fe/fe_system.h"
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const FiniteElement<dim> &fe, const unsigned int n_subdivisions = 1)
+{
+  const auto &pattern = fe.get_local_dof_sparsity_pattern();
+  deallog << fe.get_name() << ":\n";
+
+  if (!pattern.empty())
+    {
+      for (unsigned int i = 0; i < pattern.size(0); ++i)
+        {
+          for (unsigned int j = 0; j < pattern.size(1); ++j)
+            deallog << (pattern(i, j) == true ? "X" : ".");
+          deallog << "\n";
+        }
+      deallog << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  // The simplest case is a single iso Q1:
+  test<1>(FE_Q_iso_Q1<1>(2), 1);
+  // The 2 iso Q1 elements couple using their pattern:
+  test<1>(FESystem<1, 1>(FE_DGQ<1>(1), 1, FE_Q_iso_Q1<1>(2), 2), 1);
+  // The coupling between the first two to the third copy is a full coupling
+  // currently, because we don't detect this yet:
+  test<1>(FESystem<1, 1>(FE_Q_iso_Q1<1>(2), 2, FE_Q_iso_Q1<1>(2), 1), 1);
+  // Different iso_Q1 degrees always couple fully (off diagonal blocks):
+  test<1>(FESystem<1, 1>(FE_Q_iso_Q1<1>(2), 1, FE_Q_iso_Q1<1>(3), 1), 1);
+}

Diff for: tests/fe/local_sparsity_01.output

@@ -0,0 +1,36 @@
+
+DEAL::FE_Q_iso_Q1<1>(2):
+X.X
+.XX
+XXX
+
+DEAL::FESystem<1>[FE_DGQ<1>(1)-FE_Q_iso_Q1<1>(2)^2]:
+XXXXXXXX
+XXXXXXXX
+XXX.XX.X
+XX.XX.XX
+XXXXXXXX
+XXX.XX.X
+XX.XX.XX
+XXXXXXXX
+
+DEAL::FESystem<1>[FE_Q_iso_Q1<1>(2)^2-FE_Q_iso_Q1<1>(2)]:
+X.XX.XXXX
+.XX.XXXXX
+XXXXXXXXX
+X.XX.XXXX
+.XX.XXXXX
+XXXXXXXXX
+XXXXXXX.X
+XXXXXX.XX
+XXXXXXXXX
+
+DEAL::FESystem<1>[FE_Q_iso_Q1<1>(2)-FE_Q_iso_Q1<1>(3)]:
+X.XXXXX
+.XXXXXX
+XXXXXXX
+XXXX.X.
+XXX.X.X
+XXXX.XX
+XXX.XXX
+

Diff for: tests/fe/local_sparsity_q_iso_q1.cc

@@ -12,7 +12,7 @@
 //
 // ------------------------------------------------------------------------
 
-// Assemble the sparsity pattern for a Q1 finite element space and a an iso Q
+// Assemble the sparsity pattern for a Q1 finite element space and an iso Q
 // element space with the same number of subdivisions, and check we get the same
 // sparsity patterns.