Resolution of identity constraints during FEEvaluation::read_dof_values_plain()

[peterrum]

The following program:

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>

#include <deal.II/lac/affine_constraints.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

using namespace dealii;

int
main()
{
  const unsigned int dim       = 2;
  const unsigned int fe_degree = 1;
  using Number                 = double;
  using VectorizedArrayType    = VectorizedArray<Number, 1>;

  Triangulation<dim> tria;
  GridGenerator::hyper_cube(tria, 0, 1, true);

  {
    std::vector<
      GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
      periodicity_vector;

    GridTools::collect_periodic_faces(tria, 0, 1, 0, periodicity_vector);
    tria.add_periodicity(periodicity_vector);
  }

  DoFHandler<dim> dof_handler(tria);
  dof_handler.distribute_dofs(FE_Q<dim>{fe_degree});

  AffineConstraints<Number> constraints;

  {
    std::vector<GridTools::PeriodicFacePair<
      typename dealii::DoFHandler<dim>::cell_iterator>>
      periodicity_vector;
    GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vector);

    DoFTools::make_periodicity_constraints<dim, dim>(periodicity_vector,
                                                     constraints);
  }

  constraints.close();

  typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
    additional_data;
  additional_data.mapping_update_flags = update_values;

  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
  matrix_free.reinit(dof_handler,
                     constraints,
                     QGauss<dim>(fe_degree + 1),
                     additional_data);

  Vector<Number> vec;
  matrix_free.initialize_dof_vector(vec);
  vec = 1.0;
  constraints.set_zero(vec);
  // vec.print(std::cout);

  FEEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType> phi(matrix_free);
  phi.reinit(0);

  phi.read_dof_values(vec);
  for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
    std::cout << phi.begin_dof_values()[i] << " ";
  std::cout << std::endl;

  phi.read_dof_values_plain(vec);
  for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
    std::cout << phi.begin_dof_values()[i] << " ";
  std::cout << std::endl;
}

gives the following output:

1 1 1 1 // FEEvaluation::read_dof_values()
1 1 1 1 // FEEvaluation::read_dof_values_plain()

although what I have expected is:

1 1 1 1 // FEEvaluation::read_dof_values()
1 0 1 0 // FEEvaluation::read_dof_values_plain()

This means that FEEvaluation::read_dof_values_plain() applies constraints (here PBC), what is not what I expect from DoFAccessor::get_interpolated_dof_values.

The reason for this behavior is:

dealii/include/deal.II/matrix_free/dof_info.templates.h

Lines 158 to 169 in 37b1cca

	// check whether this dof is identity constrained to another
	// dof. then we can simply insert that dof and there is no
	// need to actually resolve the constraint entries
	const auto & entries = *entries_ptr;
	const types::global_dof_index n_entries = entries.size();
	if (n_entries == 1 &&
	std::abs(entries[0].second - 1.) <
	100 * std::numeric_limits<double>::epsilon())
	{
	current_dof = entries[0].first;
	goto no_constraint;
	}

with the result that

dealii/include/deal.II/matrix_free/dof_info.templates.h

Lines 250 to 253 in 37b1cca

	const bool cell_has_constraints =
	cell_has_hanging_nodes ||
	(row_starts[(cell_number + 1) * n_components].second >
	row_starts[cell_number * n_components].second);

is false.

I would say that replacing DoFs in the case of identity constraints is valid for FEEvaluation::read_dof_values() but not for FEEvaluation::read_dof_values_plain(). @kronbichler What do you think?

The automatic resolution of identity constrains during FEEvaluation::read_dof_values_plain() is often what one wants. But in the case of FEEvaluation::set_values_plain() this leads to unexpected behavior: zero velues at constrained entries.

FYI @vovannikov @mschreter

[kronbichler]

In my opinion we should fix the particular case of periodic boundary conditions with simple identity constraints (no transformation matrix) in DoFHandler/DoFAccessor, as fewer DoFs in the vectors are better and I can't see any application wanting something different. It's an incompatible change, but I would believe that it's less work on the user's side than many other incompatible changes we keep introducing over the years.

The problem is that with the current code we would need to change several places to actually go to the plain dof values, not just the code you are using, as the check here

dealii/include/deal.II/matrix_free/fe_evaluation.h

Lines 3517 to 3527 in 37b1cca

	if (apply_constraints == false &&
	(this->dof_info->row_starts[cell_dof_index].second !=
	this->dof_info->row_starts[cell_dof_index + n_components_read]
	.second ||
	((this->dof_info->hanging_node_constraint_masks.size() > 0 &&
	this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
	this->dof_info->hanging_node_constraint_masks[cell_index] !=
	internal::MatrixFreeFunctions::
	unconstrained_compressed_constraint_kind) &&
	this->dof_info->hanging_node_constraint_masks_comp
	[this->active_fe_index][this->first_selected_component])))

might fail in case there is no other constraint on the current cell. Then, we would also need to touch the logic of setting up the plain dof values.