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/12a5196b-ebb2-4104-a29b-ddb5b553473cn%40googlegroups.com.
