support DoF coupling in FESystem

Author: tjhei

source/fe/fe_system.cc

@@ -2009,6 +2009,61 @@ 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);
+
+        // Find shape functions within the same base element. If we do, grab the
+        // coupling from that base element pattern:
+        for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
+          for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
+            {
+              const auto vi = this->system_to_base_index(i);
+              const auto vj = this->system_to_base_index(j);
+
+              const auto base_index_i = vi.first.first;
+              const auto base_index_j = vj.first.first;
+              if (base_index_i == base_index_j)
+                {
+                  const auto shape_index_i = vi.second;
+                  const auto shape_index_j = vj.second;
+
+                  const auto &pattern = this->base_element(base_index_i)
+                                          .get_local_dof_sparsity_pattern();
+
+                  if (!pattern.empty())
+                    this->local_dof_sparsity_pattern(i, j) =
+                      pattern(shape_index_i, shape_index_j);
+                }
+            }
+      }
+  }
+
+
   // wait for all of this to finish
   init_tasks.join_all();
 }

tests/fe/assemble_q_iso_q1.cc tests/fe/assemble_q_iso_q1_01.cc

@@ -0,0 +1,230 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2013 - 2024 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.
+//
+// ------------------------------------------------------------------------
+
+
+// Assemble a 1d,2d,3d vector Poisson problem with FE_Q_iso_Q1 elements to
+// check that the sparsity pattern is computed correctly.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4();
+  void
+  run();
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  assemble_preconditioner();
+
+  Triangulation<dim> triangulation;
+
+  FESystem<dim>   fe_precondition;
+  DoFHandler<dim> dof_handler_precondition;
+
+  AffineConstraints<double> constraints;
+
+  SparsityPattern      prec_pattern;
+  SparseMatrix<double> preconditioner_matrix;
+
+  Vector<double> system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+                          const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 4 * std::pow(p[i], 4);
+
+  return return_value;
+}
+
+
+template <int dim>
+Step4<dim>::Step4()
+  : fe_precondition(FE_Q_iso_Q1<dim>(2), dim)
+  , dof_handler_precondition(triangulation)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+  triangulation.refine_global(1);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+  dof_handler_precondition.distribute_dofs(fe_precondition);
+
+  constraints.clear();
+  std::map<unsigned int, double> boundary_values;
+  VectorTools::interpolate_boundary_values(dof_handler_precondition,
+                                           0,
+                                           Functions::ZeroFunction<dim>(dim),
+                                           constraints);
+  constraints.close();
+
+  {
+    DynamicSparsityPattern c_sparsity(dof_handler_precondition.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler_precondition,
+                                    c_sparsity,
+                                    constraints,
+                                    false);
+    prec_pattern.copy_from(c_sparsity);
+    preconditioner_matrix.reinit(prec_pattern);
+  }
+
+  system_rhs.reinit(dof_handler_precondition.n_dofs());
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::assemble_preconditioner()
+{
+  QIterated<dim> quadrature_formula(QGauss<1>(2), 3);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values(fe_precondition,
+                          quadrature_formula,
+                          update_values | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe_precondition.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler_precondition.begin_active(),
+    endc = dof_handler_precondition.end();
+
+  for (; cell != endc; ++cell)
+    {
+      fe_values.reinit(cell);
+      cell_matrix = 0;
+      cell_rhs    = 0;
+
+      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              cell_matrix(i, j) +=
+                (fe_values.shape_grad(i, q_point) *
+                 fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
+          }
+
+      cell->get_dof_indices(local_dof_indices);
+      constraints.distribute_local_to_global(cell_matrix,
+                                             local_dof_indices,
+                                             preconditioner_matrix);
+    }
+  preconditioner_matrix.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::run()
+{
+  deallog.push("dim " + std::to_string(dim));
+
+  make_grid();
+
+  setup_system();
+  assemble_preconditioner();
+  deallog << "nnz: " << preconditioner_matrix.n_nonzero_elements() << std::endl;
+
+  deallog.pop();
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog(true);
+
+  {
+    Step4<1> test;
+    test.run();
+  }
+  {
+    Step4<2> test;
+    test.run();
+  }
+  {
+    Step4<3> test;
+    test.run();
+  }
+}

tests/fe/assemble_q_iso_q1.output tests/fe/assemble_q_iso_q1_01.output

@@ -0,0 +1,4 @@
+
+DEAL:dim 1::nnz: 11
+DEAL:dim 2::nnz: 228
+DEAL:dim 3::nnz: 3381

tests/fe/assemble_q_iso_q1_02.cc

@@ -365,17 +365,10 @@ create_mass_matrix(const Mapping<dim>        &mapping,
                     }
             }
         }
-      // transfer everything into the
-      // global object. lock the
-      // matrix meanwhile
+      // transfer everything into the global object
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
         for (unsigned int j = 0; j < dofs_per_cell; ++j)
-          if ((n_components == 1) || (cell_matrix(i, j) != 0.0))
-            /*
-              (fe.system_to_component_index(i).first ==
-              fe.system_to_component_index(j).first))
-            */
-            matrix.add(dof_indices[i], dof_indices[j], cell_matrix(i, j));
+          matrix.add(dof_indices[i], dof_indices[j], cell_matrix(i, j));
 
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
         rhs_vector(dof_indices[i]) += cell_vector(i);

tests/fe/assemble_q_iso_q1_02.output

@@ -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);
+}