Optimize sparsity pattern for FE_Q_iso_Q1

[peterrum]

depends on dealii/dealii#14197

[peterrum]

@kronbichler What is a bit annoying that the computation of the element-stiffness matrix introduces rounding errors at entries that should be actually 0 by definition. I am clearing it here by comparing the values to the max value in the matrix. A safer approach would be probably to pass in a mask here. Have ever had the same issue in one of your codes?

[kronbichler]

I tested this on a quick test problem but could not see it:

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/fe/fe_q_iso_q1.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/shape_info.h>

int main()
{
  using namespace dealii;
  const unsigned int dim = 2;
  for (unsigned int i=1; i<4; ++i)
    {
      FE_Q_iso_Q1<dim> fe(i);
      QIterated<1> quad(QGauss<1>(2), i);
      internal::MatrixFreeFunctions::ShapeInfo<double> shape_info;
      shape_info.reinit(quad, fe, 0);
      std::cout << "Values:";
      for (const auto d : shape_info.data[0].shape_values)
        std::cout << " " << d;
      std::cout << std::endl;
      std::cout << "Gradients:";
      for (const auto d : shape_info.data[0].shape_gradients)
        std::cout << " " << d;
      std::cout << std::endl;
      std::cout << "Hessians:";
      for (const auto d : shape_info.data[0].shape_hessians)
        std::cout << " " << d;
      std::cout << std::endl << std::endl;

      Triangulation<dim> tria;
      GridGenerator::hyper_cube(tria, 0., 1.);
      DoFHandler<dim> dofh(tria);
      dofh.distribute_dofs(fe);

      MappingQ<dim> mapping(1);
      MatrixFree<dim> matrix_free;
      AffineConstraints<double> constraints;
      matrix_free.reinit(mapping, dofh, constraints, quad, MatrixFree<dim>::AdditionalData());
      FEEvaluation<dim, -1> eval(matrix_free);
      eval.reinit(0);
      FullMatrix<double> cell_mass(fe.dofs_per_cell, fe.dofs_per_cell), cell_lapl(cell_mass);
      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
        {
          for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
            eval.begin_dof_values()[j] = 0.;
          eval.begin_dof_values()[i] = 1.;
          eval.evaluate(EvaluationFlags::values);
          for (unsigned int q=0; q<eval.n_q_points; ++q)
            eval.submit_value(eval.get_value(q), q);
          eval.integrate(EvaluationFlags::values);
          for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
            cell_mass(j, i) = eval.begin_dof_values()[j][0];
          
          for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
            eval.begin_dof_values()[j] = 0.;
          eval.begin_dof_values()[i] = 1.;
          eval.evaluate(EvaluationFlags::gradients);
          for (unsigned int q=0; q<eval.n_q_points; ++q)
            eval.submit_gradient(eval.get_gradient(q), q);
          eval.integrate(EvaluationFlags::gradients);
          for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
            cell_lapl(j, i) = eval.begin_dof_values()[j][0];
        }
      std::cout << "Mass" << std::endl;
      cell_mass.print_formatted(std::cout);
      std::cout << "Laplace" << std::endl;
      cell_lapl.print_formatted(std::cout);
      std::cout << std::endl;
    }
}

The point is that the numbers should not get zero because two non-zero numbers cancel each other, but because the value or gradient of the test function/solution function are exactly zero, so certain evaluations of the shape functions at quadrature points should be exactly zero. I also do not see if this could go wrong with non-Cartesian geometries or similar, because there should be a structural zero. Do you have a small case that shows it?

[kronbichler]

Or is the origin of the confusion that the Laplacian in 3D on a cube (with constant coefficients) has some additional zero entries (by coincidence). This is for the element Laplace matrix:

 3.333e-01  2.776e-17  2.776e-17 -8.333e-02            -8.333e-02 -8.333e-02 -8.333e-02 
 2.776e-17  3.333e-01 -8.333e-02  2.776e-17 -8.333e-02            -8.333e-02 -8.333e-02 
 2.776e-17 -8.333e-02  3.333e-01  2.776e-17 -8.333e-02 -8.333e-02            -8.333e-02 
-8.333e-02  2.776e-17  2.776e-17  3.333e-01 -8.333e-02 -8.333e-02 -8.333e-02            
           -8.333e-02 -8.333e-02 -8.333e-02  3.333e-01  2.776e-17  2.776e-17 -8.333e-02 
-8.333e-02            -8.333e-02 -8.333e-02  2.776e-17  3.333e-01 -8.333e-02  2.776e-17 
-8.333e-02 -8.333e-02            -8.333e-02  2.776e-17 -8.333e-02  3.333e-01  2.776e-17 
-8.333e-02 -8.333e-02 -8.333e-02            -8.333e-02  2.776e-17  2.776e-17  3.333e-01 

Note that this is coincidence and one should not rely on this fact. (I guess we could argue that the filter to determine the sparsity pattern from a matrix might be too optimistic, and we should rather check if the tensor product indices [i_0, i_1, i_2] for rows are not more than +1 to -1 apart from [j_0, j_1, j_2].

[peterrum]

The problems are the entries of type 2.776e-17 (from your example). Other entries (empty entry in the output of full matrix), are hard zero.

Or is the origin of the confusion that the Laplacian in 3D on a cube (with constant coefficients) has some additional zero entries (by coincidence). This is for the element Laplace matrix:

I don't think we are exploiting any structure on the element-level, but let me check that later.

[kronbichler]

But the entries 2.776e-17 are zero within a cell for a regular FE_Q<3>(1). Also what is empty in the 8x8 matrix I posted above is just coincidence, just that the roundoff errors accumulate differently, I am sure this can be triggered with different numbers to produce 1e-16/1e-17 as well. All these zeros go away as soon as we switch from a square to a rectangle if I remember correctly, and definitely also go away when the elements are rotated or deformed. The 8x8 matrix should be dense in general, as all 8 shape functions overlap.

[peterrum]

Here is an example (2 subdivisions, 2D, 2 comonents, lex numbering):

The picture indicates there is a coupling between dof 6 ("3rd row") of component 0 and dof 0/1 ("1st row") of component 1. However, this should not be the case, since the sparsity pattern should be a tensor product.

[peterrum]

@vovannikov referencing #209

[peterrum]

Merging, since the PR in deal.II is merge. I will revisit the zeroing once we observe problems (Trilinos will kick us out).

@vovannikov Once we have a solution in #204, let's marry number of refinements and number of subdivisions.