modification of the mass-matrix in tutorial-67

[Luca Arpaia]

Dear deal.ii user group,

I am developing a discontinuous Galerkin code. I have started from the tutorial-67 which is great but now I need to modify the quadrature formula to integrate the mass matrix. The integration of the tutorial rely on the class "CellwiseInverseMassMatrix" where a very specific quadrature, with the number of points equal to the number of degrees of freedom, is used. I need something more general with an arbitrary quadrature formula.

I have started implemented the inversion of the cell-wise matrices with the class "FullMatrix", something like:

FEValues<dim> fe_values(fe, quadrature_formula
  update_values | update_JxW_values);

for (unsigned int q=0; q<n_q_points; ++q)

  for (unsigned int i=0; i<dofs_per_cell; ++i)

    for (unsigned int j=0; j<dofs_per_cell; ++j)

      cell_matrix(i,j) += fe_values.shape_values(i,q) *

                            fe_values.shape_values(j,q) *

                              fe_values.JxW(q);  // (1)

cell_matrix.gauss_jordan();
cell_matrix.vmult(cell_dst, cell_src);

However I wonder if I could tensorize this operation making it more efficient: computing smaller 1d mass-matrices and later computing the multidimensional mass matrix by tensorization. I was able to recover the 1d shape functions at the quadrature points:

FEEvaluation<dim, degree, n_q_points_1d, 1, Number> fe_eval(data, 0, 1);

AlignedVector<VectorizedArrayType> shape_values =

  fe_eval.get_shape_info().data.front().shape_values;

to later build the 1d mass-matrix as:

for (int i = 0; i < degree + 1; ++i)

  for (int j = 0; j < degree + 1; ++j)

    for (int q = 0; q < n_q_points_1d; ++q)

      {

        double phi_i = shape_values[i * n_q_points_1d + q][cell_in_lane];

        double phi_j = shape_values[j * n_q_points_1d + q][cell_in_lane];

        cell_matrix_1d(i,j) += phi_i * phi_j *JxW_1d[q] ;

      }

The problem that I face is where to pick the correct JxW_1d. In general the Jacobian determinant for the quadrilateral is non-constant over the cell, so do I have to abandon the tensorization strategy and go back to option (1)?

As an aside, if I have to go back to option (1) is there a way to get the shape values without defining the object:

FEValues<dim> fe_values(fe, quadrature_formula
  update_values | update_JxW_values);

since I have already defined an object to evaluate functions at quadrature points:

FEEvaluation<dim, degree, n_q_points_1d, 1, Number> fe_eval(data, 0, 1);

thank you very much,

Luca

[Abbas Ballout]

Luca, 

I believe that you might be missing something.
In the reinit function you have two quadratures:

const std::vector<Quadrature<1>> quadratures = {QGauss<1>(n_q_points_1d),

  QGauss<1>(fe_degree + 1)};

The mass matrix is using the second quadrature which is selected through the 1 passed to: 
FEEvaluation<dim, degree, degree + 1, dim + 2, Number> phi(data, 0, 1); in later function calls (the inverse, project...) 
Hope this helps, 

Abbas