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>
<https://dealii.org/current/doxygen/deal.II/classQGauss.html>
(n_q_points_1d),
QGauss<1>
<https://dealii.org/current/doxygen/deal.II/classQGauss.html>(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>
<https://dealii.org/current/doxygen/deal.II/classFEEvaluation.html>
phi(data, 0, 1); in later function calls (the inverse, project...)
Hope this helps,
Abbas
On Thursday, September 11, 2025 at 12:50:18 PM UTC+2 [email protected]
wrote:
> 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
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/707eb626-deae-4efa-ae87-ce70ea0ed602n%40googlegroups.com.
