Removing from existing AffineConstraints

[Arthur Bawin]

Hello deal.II community,

I'm having an issue with what something I believe to be very simple. I have an existing AffineConstraints "constraints", from which I want to remove a subset of constraints. This subset is the set of fluid velocity dofs on some boundary:

IndexSet velocity_dofs =
    DoFTools::extract_boundary_dofs(dof_handler, velocity_mask, {id});
IndexSet subset = constraints.get_local_lines();
subset.subtract_set(velocity_dofs);

Then I try to re-create the constraints using get_view() and merge() as indicated in affine_constraints.h :

auto tmp_constraints = constraints.get_view(subset);
constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
constraints.close();
constraints.merge(tmp_constraints);

But I get this error from get_view(), in debug and on a single MPI proc:

--------------------------------------------------------
An error occurred in line <2589> of file </home/bawina/Code/dealii/include/deal.II/lac/affine_constraints.h> in function
    dealii::types::global_dof_index dealii::AffineConstraints<number>::calculate_line_index(size_type) const [with number = double; dealii::types::global_dof_index = unsigned int; size_type = unsigned int]
The violated condition was:
    local_lines.is_element(line_n)
Additional information:
    The index set given to this constraints object indicates constraints
    for degree of freedom 10371 should not be stored by this object, but a
    constraint is being added.

Stacktrace:
-----------
#0  ./monolithic_fsi: dealii::AffineConstraints<double>::calculate_line_index(unsigned int) const
#1  ./monolithic_fsi: dealii::AffineConstraints<double>::is_constrained(unsigned int) const
#2  /usr/local/lib/libdeal_II.g.so.9.7.0-pre: dealii::AffineConstraints<double>::add_constraint(unsigned int, dealii::ArrayView<std::pair<unsigned int, double> const, dealii::MemorySpace::Host> const&, double)
#3  /usr/local/lib/libdeal_II.g.so.9.7.0-pre: dealii::AffineConstraints<double>::get_view(dealii::IndexSet const&) const

[...]

It comes from view.add_constraint(mask.index_within_set(line.index), ..., ...) in get_view(), is it a "shift" issue?

Copying and filtering by hand works:

{
    AffineConstraints<double> filtered;
    filtered.reinit(locally_owned_dofs, locally_relevant_dofs);
    for (const auto &line : constraints.get_lines())
    {
      if (velocity_dofs.is_element(line.index))
        continue;

      filtered.add_line(line.index);
      filtered.add_entries(line.index, line.entries);
      filtered.set_inhomogeneity(line.index, line.inhomogeneity);

      for (const auto &entry : line.entries)
        AssertThrow(!velocity_dofs.is_element(entry.first),
                    ExcMessage("Constraint involve a boundary velocity dof"));
    }
    filtered.close();
    constraints = std::move(filtered);
  }

Thank you for your time!
Arthur

[Arthur Bawin]

(The "subset" in the description and in the code snippet are the complement of one another, sorry for the inconsistency, I really mean to say that subset = all constrained dofs \ velocity dofs)

[Wolfgang Bangerth]

On 11/19/25 14:05, Arthur Bawin wrote:
>
> Then I try to re-create the constraints using get_view() and merge() as
> indicated in affine_constraints.h :

You don't show the whole code, so it's difficult to say. But get_view()
doesn't just remove constraints outside the "viewed" subset, but it
re-enumerates. For example, if you have 10 DoFs and your
AffineConstraint stores constraints for indices 1, 2, 5. Then let's say
you call get_view() with an index (sub)set 4...8 inclusive, then the
result is not an AffineConstraint of size 10 that stores a constraint
for index 5. Rather, it is an AffineConstraint of size 5 that stores a
constraint for index 1 (the second index within the range 4...8).

Does that make sense as perhaps the cause for your issue?

Best
W.

[Arthur Bawin]

I'm a bit confused then. The following tests three get_view() with IndexSets : your example (slightly modified), the example from the doc of get_view(), and then a slight modification of it. Your example returns 0-based indices, but the one from the doc returns 11 as first index, which is 1-based. The last example throws an error. Is it a bug or am I missing something obvious?

My output is :

A :{[1,2], 5}
B :{[0,1], [4,9]}
View :{1, 3}
A :{[20,39]}
B :{1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89, 91, 93, 95, 97, 99}
View B:{[11,20]}
C :{22, 26}

--------------------------------------------------------
An error occurred in line <1816> of file </home/bawina/Code/dealii/include/deal.II/base/index_set.h> in function
    void dealii::IndexSet::add_range(size_type, size_type)
The violated condition was:
    (begin < index_space_size) || ((begin == index_space_size) && (end == index_space_size))
Additional information:
    Index 22 is not in the half-open range [0,2).

Stacktrace:
-----------
#0  ./indexset: dealii::IndexSet::add_range(unsigned int, unsigned int)
#1  /usr/local/lib/libdeal_II.g.so.9.7.0-pre: dealii::IndexSet::get_view(dealii::IndexSet const&) const
#2  ./indexset: main
--------------------------------------------------------

(Edit: I accidentally replied to W. Bangerth only, not to the list.)

Thanks,
Arthur

    {
        IndexSet A(10);
        A.add_index(1);
        A.add_index(2);
        A.add_index(5);
        std::cout << "A :";
        A.print(std::cout);

        IndexSet B(10);
        B.add_index(0);
        B.add_index(1);
        B.add_index(4);
        B.add_index(5);
        B.add_index(6);
        B.add_index(7);
        B.add_index(8);
        B.add_index(9);
        std::cout << "B :";
        B.print(std::cout);

        IndexSet view = A.get_view(B);
        std::cout << "View :";
        view.print(std::cout);
      }

      {
        const unsigned int N = 100;

        IndexSet A(N);
        A.add_range(20, 40);
        std::cout << "A :";
        A.print(std::cout);

        // Odd numbers in [0,N)
        IndexSet B(N);
        for (unsigned int i = 1; i < N; i += 2)
          B.add_index(i);
        std::cout << "B :";
        B.print(std::cout);

        {
          IndexSet view = A.get_view(B);
          std::cout << "View B:";
          view.print(std::cout);
        }

        IndexSet C(N);
        C.add_index(22);
        C.add_index(26);
        std::cout << "C :";
        C.print(std::cout);

        {
          IndexSet view = A.get_view(C);
          std::cout << "View C :";
          view.print(std::cout);