Behavior at shared faces between p-refined and h-refined cells

[peterrum]

Let me post a test @mschreter and I wrote last year (Nov 9, 2022) and which reminds me of the problematic configuration we have here. Left cell is refined and linear and right cell is not refined and quadratic. The DoFs on the shared face (there one DoF on each side) are constrained, degenerating the function space at the face to a linear function.

#include <deal.II/base/convergence_table.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/parameter_handler.h>

#include <deal.II/distributed/cell_weights.h>
#include <deal.II/distributed/tria.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/grid/grid_generator.h>

#include <deal.II/hp/fe_values.h>

#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>

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

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_creator.h>

using namespace dealii;

const bool manual_constraint = false;
const unsigned int fe_degree = 1;

/**
  Test contraints between a locally refined cell and a cell with
  two subdivisions.

  The two DoFs in the center of the shared face are constrained
  by the other ones:
     2 11:  0.5
     2 15:  0.5
    12 11:  0.5
    12 15:  0.5
 */
int
main(int argc, char *argv[])
{
  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

  const unsigned int dim = 2;

  // 1) create system
  hp::FECollection<dim> fe;
  fe.push_back(FE_Q_iso_Q1<dim>(fe_degree));
  fe.push_back(FE_Q_iso_Q1<dim>(fe_degree + 1));

  hp::QCollection<dim> quadrature_collection;
  quadrature_collection.push_back(QIterated<dim>(QGauss<1>(2), fe_degree));
  quadrature_collection.push_back(QIterated<dim>(QGauss<1>(2), fe_degree + 1));

  hp::QCollection<dim> quadrature;

  for (const auto &f : fe)
    quadrature.push_back(Quadrature<dim>(f.get_unit_support_points()));

  Triangulation<dim> tria;
  GridGenerator::subdivided_hyper_rectangle(tria,
                                            {2, 1},
                                            Point<dim>(0.0, 0.0),
                                            Point<dim>(2.0, 1.0));

  tria.begin()->set_refine_flag();
  tria.execute_coarsening_and_refinement();

  DoFHandler<dim> dof_handler(tria);

  for (const auto &cell : dof_handler.active_cell_iterators())
    if (cell->center()[0] < 1.0)
      cell->set_active_fe_index(0);
    else
      cell->set_active_fe_index(1);

  dof_handler.distribute_dofs(fe);

  MappingQ1<dim> mapping;

  // 2) match support points with DoF indices
  hp::FEValues<dim> fe_values_hp(fe, quadrature, update_quadrature_points);

  std::vector<types::global_dof_index> dof_indices;

  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      fe_values_hp.reinit(cell);
      const auto &fe_values = fe_values_hp.get_present_fe_values();

      dof_indices.resize(fe_values.dofs_per_cell);

      cell->get_dof_indices(dof_indices);

      for (const auto q : fe_values.quadrature_point_indices())
        std::cout << dof_indices[q] << ": " << fe_values.quadrature_point(q)
                  << std::endl;

      std::cout << std::endl;
    }

  std::cout << dof_handler.n_dofs() << std::endl;

  // set up constraints
  AffineConstraints<double> constraints;
  if(true)
    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  else
    {
      constraints.add_line(12);
      constraints.add_entry(12, 2, 1.0);
    }
  constraints.close();

  constraints.print(std::cout);
  std::cout << std::endl;

  DynamicSparsityPattern dsp(dof_handler.n_dofs());
      DoFTools::make_sparsity_pattern(dof_handler,
                                      dsp,
                                      constraints);

  SparsityPattern sparsity_pattern;
      sparsity_pattern.copy_from(dsp);

  SparseMatrix<double> matrix;
      matrix.reinit(sparsity_pattern);
      MatrixCreator::create_laplace_matrix<dim, dim>(
        dof_handler,
        quadrature_collection,
        matrix,
        nullptr,
        constraints);

  FullMatrix<double> full_matrix(dof_handler.n_dofs(),dof_handler.n_dofs());
  full_matrix.copy_from(matrix);

  full_matrix.print_formatted(std::cout, 2, false, 7, "0.00");
  std::cout << std::endl;

  LAPACKFullMatrix<double> lapack_full_matrix;
  lapack_full_matrix.copy_from(full_matrix);
  lapack_full_matrix.compute_eigenvalues();

  std::vector<double> eigenvalues;

  for(unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
    eigenvalues.push_back(lapack_full_matrix.eigenvalue(i).real());

  std::sort(eigenvalues.begin(), eigenvalues.end());

  for(const auto i : eigenvalues)
    std::cout << i << " ";
  std::cout <<std::endl;

}

Note: one should be able to introduce identity constraints.

[peterrum]

with "hanging-node" constraints (system matrix and eigenvalues):

    2 11:  0.5
    2 15:  0.5
    12 11:  0.5
    12 15:  0.5

   0.67    0.00    0.00   -0.17   -0.17    0.00   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
   0.00    0.67    0.00   -0.17    0.00   -0.17   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
   0.00    0.00    1.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
  -0.17   -0.17    0.00    1.33   -0.33   -0.33   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
  -0.17    0.00    0.00   -0.33    1.33    0.00   -0.33    0.00    0.00    0.00    0.00   -0.33    0.00    0.00    0.00   -0.17 
   0.00   -0.17    0.00   -0.33    0.00    1.33   -0.33    0.00    0.00    0.00    0.00   -0.17    0.00    0.00    0.00   -0.33 
  -0.33   -0.33    0.00   -0.33   -0.33   -0.33    2.67    0.00    0.00    0.00    0.00   -0.50    0.00    0.00    0.00   -0.50 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.67   -0.17   -0.17   -0.33    0.00    0.00    0.00    0.00    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.17    1.33   -0.33   -0.33   -0.33    0.00    0.00    0.00   -0.17 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.17   -0.33    1.33   -0.33    0.00    0.00   -0.17   -0.33    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.33   -0.33   -0.33    2.67   -0.50    0.00   -0.33   -0.33   -0.50 
   0.00    0.00    0.00    0.00   -0.33   -0.17   -0.50    0.00   -0.33    0.00   -0.50    1.67    0.00    0.00   -0.17    0.33 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    1.33    0.00    0.00    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.17   -0.33    0.00    0.00    0.67   -0.17    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.33   -0.33   -0.17    0.00   -0.17    1.33   -0.33 
   0.00    0.00    0.00    0.00   -0.17   -0.33   -0.50    0.00   -0.17    0.00   -0.50    0.33    0.00    0.00   -0.33    1.67 

3.05222e-16 0.248545 0.622743 0.627322 0.842638 1.04866 1.12594 1.33333 1.33333 1.37268 1.47005 1.58465 1.79745 2.21538 2.90535 3.47194 

with identity constraints:

    12 2:  1

   0.67    0.00    0.00   -0.17   -0.17    0.00   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
   0.00    0.67    0.00   -0.17    0.00   -0.17   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
   0.00    0.00    2.67    0.00   -0.33   -0.33   -0.33    0.00   -0.33    0.00   -0.33   -0.33    0.00    0.00   -0.33   -0.33 
  -0.17   -0.17    0.00    1.33   -0.33   -0.33   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00 
  -0.17    0.00   -0.33   -0.33    1.33    0.00   -0.33    0.00    0.00    0.00    0.00   -0.17    0.00    0.00    0.00    0.00 
   0.00   -0.17   -0.33   -0.33    0.00    1.33   -0.33    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.17 
  -0.33   -0.33   -0.33   -0.33   -0.33   -0.33    2.67    0.00    0.00    0.00    0.00   -0.33    0.00    0.00    0.00   -0.33 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.67   -0.17   -0.17   -0.33    0.00    0.00    0.00    0.00    0.00 
   0.00    0.00   -0.33    0.00    0.00    0.00    0.00   -0.17    1.33   -0.33   -0.33   -0.17    0.00    0.00    0.00    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.17   -0.33    1.33   -0.33    0.00    0.00   -0.17   -0.33    0.00 
   0.00    0.00   -0.33    0.00    0.00    0.00    0.00   -0.33   -0.33   -0.33    2.67   -0.33    0.00   -0.33   -0.33   -0.33 
   0.00    0.00   -0.33    0.00   -0.17    0.00   -0.33    0.00   -0.17    0.00   -0.33    1.33    0.00    0.00    0.00    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    1.33    0.00    0.00    0.00 
   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.17   -0.33    0.00    0.00    0.67   -0.17    0.00 
   0.00    0.00   -0.33    0.00    0.00    0.00    0.00    0.00    0.00   -0.33   -0.33    0.00    0.00   -0.17    1.33   -0.17 
   0.00    0.00   -0.33    0.00    0.00   -0.17   -0.33    0.00    0.00    0.00   -0.33    0.00    0.00    0.00   -0.17    1.33 

-3.12228e-16 0.248545 0.622743 0.627322 0.779393 1.04866 1.12594 1.25832 1.33333 1.37268 1.58465 1.78779 1.79745 2.90535 2.92792 3.24658