Merge pull request #16099 from bangerth/view

Add AffineConstraints::get_view().

doc/news/changes/minor/20231005Bangerth

@@ -0,0 +1,5 @@
+New: The function AffineConstraints::get_view() can be used to obtain
+the constraints that correspond to a specific subset of degrees of
+freedom.
+<br>
+(Wolfgang Bangerth, 2023/10/05)

include/deal.II/lac/affine_constraints.h

@@ -942,6 +942,57 @@ class AffineConstraints : public Subscriptor
   void
   shift(const size_type offset);
 
+  /**
+   * This function provides a "view" into a constraint
+   * object. Specifically, given a "mask" index set that describes
+   * which constraints we are interested in, it returns an
+   * AffineConstraints object that contains only those constraints
+   * that correspond to degrees of freedom that are listed in the
+   * mask, with indices shifted so that they correspond to the
+   * position within the mask. This process is the same as how
+   * IndexSet::get_view() computes the shifted indices. The function
+   * is typically used to extract from an AffineConstraints object
+   * corresponding to a DoFHandler only those constraints that
+   * correspond to a specific variable (say, to the velocity in a
+   * Stokes system) so that the resulting AffineConstraints object can
+   * be applied to a single block of a block vector of solutions; in
+   * this case, the mask would be the index set of velocity degrees of
+   * freedom, as a subset of all degrees of freedom.
+   *
+   * This function can only work if the degrees of freedom selected by
+   * the mask are constrained only against other degrees of freedom
+   * that are listed in the mask. In the example above, this means
+   * that constraints for the selected velocity degrees of freedom are
+   * only against other velocity degrees of freedom, but not against
+   * any pressure degrees of freedom. If that is not so, an assertion
+   * will be triggered.
+   *
+   * A typical case where this function is useful is as follows. Say,
+   * you have a block linear system in which you have blocks
+   * corresponding to variables $(u,p,T,c)$ (which you can think of as
+   * velocity, pressure, temperature, and chemical composition -- or
+   * whatever other kind of problem you are currently considering in
+   * your own work). Let's assume we have developed a linear solver or
+   * preconditioner that first solves the coupled $u$-$T$ system, and
+   * once that is done, solves the $p$-$c$ system. In this case, it is
+   * often useful to set up block vectors with only two components
+   * corresponding to the $u$ and $T$ components, and later for only
+   * the $p$-$c$ components of the solution. As part of this, you will
+   * want to apply constraints (using the distribute() function of this
+   * class) to only the 2-block vector, but for this you need to obtain
+   * an AffineConstraints object that represents only those constraints
+   * that correspond to the variables in question, and in the order
+   * in which they appear in the 2-block vector rather than in global
+   * 4-block vectors. This function allows you to extract such an
+   * object corresponding to a subset of constraints by applying a mask
+   * to the global constraints object that corresponds to the
+   * variables we're currently interested in. For the $u$-$T$ system,
+   * we need a mask that contains all indices of $u$ degrees of
+   * freedom as well as all indices of $T$ degrees of freedom.
+   */
+  AffineConstraints
+  get_view(const IndexSet &mask) const;
+
   /**
    * Clear all entries of this matrix. Reset the flag determining whether new
    * entries are accepted or not.

include/deal.II/lac/affine_constraints.templates.h

@@ -1082,6 +1082,68 @@ AffineConstraints<number>::shift(const size_type offset)
 
 
 
+template <typename number>
+AffineConstraints<number>
+AffineConstraints<number>::get_view(const IndexSet &mask) const
+{
+  Assert(sorted == true, ExcMatrixNotClosed());
+
+  AffineConstraints<number> view;
+
+  // If index sets were associated with the current object, take views
+  // of those as well:
+  if (locally_owned_dofs != IndexSet())
+    {
+      Assert(locally_owned_dofs.size() == mask.size(),
+             ExcMessage("This operation only makes sense if the index "
+                        "space described by the mask is the same as "
+                        "the index space associated with the "
+                        "AffineConstraints object into which you take "
+                        "a view."));
+      Assert(local_lines.size() == mask.size(),
+             ExcMessage("This operation only makes sense if the index "
+                        "space described by the mask is the same as "
+                        "the index space associated with the "
+                        "AffineConstraints object into which you take "
+                        "a view."));
+
+      view.reinit(locally_owned_dofs.get_view(mask),
+                  local_lines.get_view(mask));
+    }
+
+  for (const ConstraintLine &line : lines)
+    if (mask.is_element(line.index))
+      {
+        const size_type row = mask.index_within_set(line.index);
+        view.add_line(row);
+        for (const std::pair<size_type, number> &entry : line.entries)
+          {
+            Assert(
+              mask.is_element(entry.first),
+              ExcMessage(
+                "In creating a view of an AffineConstraints "
+                "object, the constraint on degree of freedom " +
+                std::to_string(line.index) + " (which corresponds to the " +
+                std::to_string(row) +
+                "th degree of freedom selected in the mask) "
+                "is constrained against degree of freedom " +
+                std::to_string(entry.first) +
+                ", but this degree of freedom is not listed in the mask and "
+                "consequently cannot be transcribed into the index space "
+                "of the output object."));
+            view.add_entry(row,
+                           mask.index_within_set(entry.first),
+                           entry.second);
+          }
+        view.set_inhomogeneity(row, line.inhomogeneity);
+      }
+
+  view.close();
+  return view;
+}
+
+
+
 template <typename number>
 void
 AffineConstraints<number>::clear()

tests/mpi/affine_constraints_get_view_01.cc

@@ -0,0 +1,210 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check AffineConstraints<double>.get_view()
+//
+// The way this test works is by setting up a Stokes system with a
+// Q2xQ1 element on an adaptively refined mesh for which we compute
+// the corresponding constraints (which will contain constraints for
+// both the velocity and pressure degrees of freedom). Then we create
+// a random vector to which we apply distribute(). This is all
+// existing functionality.
+//
+// Next, we create views into the constraints object for both the
+// velocity variable only and the pressure variable, and apply these
+// views to the individual blocks of the (original) random
+// vector. This must result in the same output vector as above.
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  deallog << "dim=" << dim << std::endl;
+
+  const unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  const unsigned int n_processes =
+    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+  // Set up of triangulation and dof_handler
+
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+
+  const double R0 = 6371000. - 2890000.;
+  const double R1 = 6371000. - 35000.;
+  GridGenerator::hyper_shell(tr, Point<dim>(), R0, R1, 12, true);
+
+  tr.refine_global(1);
+  for (unsigned int step = 0; step < 5; ++step)
+    {
+      typename Triangulation<dim>::active_cell_iterator cell =
+                                                          tr.begin_active(),
+                                                        endc = tr.end();
+
+      for (; cell != endc; ++cell)
+        if (Testing::rand() % 42 == 1)
+          cell->set_refine_flag();
+
+      tr.execute_coarsening_and_refinement();
+    }
+  deallog << "Cells on process #" << myid << ": " << tr.n_active_cells()
+          << std::endl;
+
+  FESystem<dim>   fe(FE_Q<dim>(2), dim, FE_Q<dim>(1), 1);
+  DoFHandler<dim> dof_handler(tr);
+
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
+  stokes_sub_blocks[dim] = 1;
+  DoFRenumbering::component_wise(dof_handler, stokes_sub_blocks);
+
+  deallog << "Total number of DoFs: " << dof_handler.n_dofs() << std::endl;
+
+  const std::vector<types::global_dof_index> stokes_dofs_per_block =
+    DoFTools::count_dofs_per_fe_block(dof_handler, stokes_sub_blocks);
+
+  const types::global_dof_index n_u = stokes_dofs_per_block[0],
+                                n_p = stokes_dofs_per_block[1];
+
+  const IndexSet &stokes_locally_owned_index_set =
+    dof_handler.locally_owned_dofs();
+  const IndexSet stokes_locally_relevant_set =
+    DoFTools::extract_locally_relevant_dofs(dof_handler);
+
+  std::vector<IndexSet> stokes_partitioning;
+  stokes_partitioning.push_back(
+    stokes_locally_owned_index_set.get_view(0, n_u));
+  stokes_partitioning.push_back(
+    stokes_locally_owned_index_set.get_view(n_u, n_u + n_p));
+
+  std::vector<IndexSet> stokes_relevant_partitioning;
+  stokes_relevant_partitioning.push_back(
+    stokes_locally_relevant_set.get_view(0, n_u));
+  stokes_relevant_partitioning.push_back(
+    stokes_locally_relevant_set.get_view(n_u, n_u + n_p));
+
+  // Set up a constraints object from hanging node constraints and
+  // boundary conditions.
+  AffineConstraints<double> stokes_constraints(stokes_locally_owned_index_set,
+                                               stokes_locally_relevant_set);
+
+  DoFTools::make_hanging_node_constraints(dof_handler, stokes_constraints);
+
+  const FEValuesExtractors::Vector velocity_components(0);
+  VectorTools::interpolate_boundary_values(
+    dof_handler,
+    0,
+    Functions::ZeroFunction<dim>(dim + 1),
+    stokes_constraints,
+    fe.component_mask(velocity_components));
+
+  std::set<types::boundary_id> no_normal_flux_boundaries;
+  no_normal_flux_boundaries.insert(1);
+  VectorTools::compute_no_normal_flux_constraints(dof_handler,
+                                                  0,
+                                                  no_normal_flux_boundaries,
+                                                  stokes_constraints);
+  stokes_constraints.close();
+  deallog << "Number of constrains: " << stokes_constraints.n_constraints()
+          << std::endl;
+
+  // Now set up a block vector for this situation. Fill it with random numbers.
+  TrilinosWrappers::MPI::BlockVector random_vector;
+  random_vector.reinit(stokes_partitioning, MPI_COMM_WORLD);
+  for (const unsigned int i : stokes_locally_owned_index_set)
+    random_vector(i) = random_value<double>(-100, 100);
+  random_vector.compress(VectorOperation::insert);
+
+
+  // Check 1: Call distribute() on a vector that contains
+  // everything, with an AffineConstraint object that contains all
+  // constraints.
+  TrilinosWrappers::MPI::BlockVector check_vector_1;
+  {
+    check_vector_1.reinit(stokes_partitioning, MPI_COMM_WORLD);
+    check_vector_1 = random_vector;
+
+    stokes_constraints.distribute(check_vector_1);
+  }
+
+  // Check 2: Call distribute() on the two blocks individually, with
+  // AffineConstraint objects that only contain views.
+  TrilinosWrappers::MPI::BlockVector check_vector_2;
+  {
+    check_vector_2.reinit(stokes_partitioning, MPI_COMM_WORLD);
+    check_vector_2 = random_vector;
+
+
+    IndexSet velocity_dofs(dof_handler.n_dofs());
+    velocity_dofs.add_range(0, n_u);
+    velocity_dofs.compress();
+
+    const AffineConstraints<double> velocity_constraints =
+      stokes_constraints.get_view(velocity_dofs);
+
+
+    IndexSet pressure_dofs(dof_handler.n_dofs());
+    pressure_dofs.add_range(n_u, n_u + n_p);
+    pressure_dofs.compress();
+
+    const AffineConstraints<double> pressure_constraints =
+      stokes_constraints.get_view(pressure_dofs);
+
+    velocity_constraints.distribute(check_vector_2.block(0));
+    pressure_constraints.distribute(check_vector_2.block(1));
+  }
+
+  // Now the assertion part: These vectors are the same:
+  Assert(check_vector_1 == check_vector_2, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog();
+
+  test<2>();
+  test<3>();
+}

tests/mpi/affine_constraints_get_view_01.with_trilinos=true.mpirun=1.output

@@ -0,0 +1,11 @@
+
+DEAL::dim=2
+DEAL::Cells on process #0: 72
+DEAL::Total number of DoFs: 878
+DEAL::Number of constrains: 306
+DEAL::OK
+DEAL::dim=3
+DEAL::Cells on process #0: 509
+DEAL::Total number of DoFs: 18124
+DEAL::Number of constrains: 7698
+DEAL::OK

tests/mpi/affine_constraints_get_view_01.with_trilinos=true.mpirun=2.output

@@ -0,0 +1,11 @@
+
+DEAL::dim=2
+DEAL::Cells on process #0: 51
+DEAL::Total number of DoFs: 853
+DEAL::Number of constrains: 170
+DEAL::OK
+DEAL::dim=3
+DEAL::Cells on process #0: 481
+DEAL::Total number of DoFs: 20074
+DEAL::Number of constrains: 6280
+DEAL::OK

tests/mpi/affine_constraints_get_view_01.with_trilinos=true.mpirun=4.output

@@ -0,0 +1,11 @@
+
+DEAL::dim=2
+DEAL::Cells on process #0: 45
+DEAL::Total number of DoFs: 903
+DEAL::Number of constrains: 123
+DEAL::OK
+DEAL::dim=3
+DEAL::Cells on process #0: 215
+DEAL::Total number of DoFs: 9429
+DEAL::Number of constrains: 2582
+DEAL::OK