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();
}
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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();
+ }
+}
--- /dev/null
+
+DEAL:dim 1::nnz: 11
+DEAL:dim 2::nnz: 228
+DEAL:dim 3::nnz: 3381
}
}
}
- // 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);
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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]:
+XX..XXXX
+XX..XXXX
+..XXXXXX
+..XXXXXX
+XXXXXXXX
+XXXXXXXX
+XXXXXXXX
+XXXXXXXX
+
+DEAL::FESystem<1>[FE_Q_iso_Q1<1>(2)^2-FE_Q_iso_Q1<1>(2)]:
+XXX..XXXX
+XXX..XXXX
+XXXXX.XXX
+..XXXXXXX
+..XXXXXXX
+XX.XXXXXX
+XXXXXXXXX
+XXXXXXXXX
+XXXXXXXXX
+
+DEAL::FESystem<1>[FE_Q_iso_Q1<1>(2)-FE_Q_iso_Q1<1>(3)]:
+XX.XXXX
+XXX.XX.
+.XXXXXX
+X.XXX.X
+XXXXXXX
+XXX.XXX
+X.XXXXX
+
//
// ------------------------------------------------------------------------
-// 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.