Loop-invariant reconstruction of the consistent mass matrix in QSConvectionDiffusionExplicit

[rzhang20]

Description

While inspecting the explicit convection–diffusion element implementation, I noticed that the consistent mass matrix is manually reconstructed on every call, even though its normalized structure is invariant for a given element type.

In the current implementation, the mass matrix is assembled by repeatedly defining constants and assigning fixed coefficients:

Kratos/applications/ConvectionDiffusionApplication/custom_elements/qs_convection_diffusion_explicit.cpp

Lines 177 to 181 in 96eeda3

	const double one_six = 1.0 / 6.0;
	const double one_twelve = 1.0 / 12.0;
	rMassMatrix(0,0) = one_six; rMassMatrix(0,1) = one_twelve; rMassMatrix(0,2) = one_twelve;
	rMassMatrix(1,0) = one_twelve; rMassMatrix(1,1) = one_six; rMassMatrix(1,2) = one_twelve;
	rMassMatrix(2,0) = one_twelve; rMassMatrix(2,1) = one_twelve; rMassMatrix(2,2) = one_six;

This pattern appears in the section where the mass matrix entries are filled manually.

However, for linear simplex elements, the structure of the consistent mass matrix is mathematically invariant:

• 2D linear triangle
  • diagonal entries: 1/6
  • off-diagonal entries: 1/12
• 3D linear tetrahedron
  • diagonal entries: 1/10
  • off-diagonal entries: 1/20

Only the final scaling factor (area or volume) depends on the actual geometry.

Despite this, the current code:

• redeclares the same double constants,
• and assigns the same fixed ratios entry by entry,
  on every invocation of the mass matrix calculation.

This is a classic case of loop-invariant static data being recomputed repeatedly.

---

Impact

From a performance and maintainability perspective, this has several drawbacks:

• Unnecessary repeated work: the same invariant coefficients are rebuilt on every call.
• Reduced clarity: the mathematical invariance of the mass matrix structure is not explicitly reflected in the code.
• Missed optimization opportunities:
  • prevents reuse of a precomputed matrix,
  • prevents taking advantage of optimized matrix copy / scaling operations.

While the cost of a single reconstruction is small, this function is typically called many times in explicit schemes and time-stepping loops, so the overhead can accumulate.

Conceptually, the current logic is equivalent to “rebuilding the same mold every time” instead of “reusing the mold and scaling it by the element measure”.

---

Suggested fix

The normalized (unit) mass matrix can be extracted as a static const matrix, initialized once, and reused by scaling it with the element area (or volume).

For example, for the 2D triangular case:

template <>
void QSConvectionDiffusionExplicit<2,3>::CalculateMassMatrix(
    MatrixType& rMassMatrix,
    const ProcessInfo& rCurrentProcessInfo)
{
    KRATOS_TRY;

    const unsigned int local_size = 3;

    // Ensure correct size (can be omitted or replaced by an assert if guaranteed by Kratos)
    if (rMassMatrix.size1() != local_size)
        rMassMatrix.resize(local_size, local_size, false);

    // Loop-invariant: unit consistent mass matrix for a linear triangle
    // Initialized once
    static const MatrixType unit_mass_matrix_2d = []{
        MatrixType m(3, 3);
        const double d = 1.0 / 6.0;   // diagonal
        const double n = 1.0 / 12.0;  // off-diagonal
        m(0,0)=d; m(0,1)=n; m(0,2)=n;
        m(1,0)=n; m(1,1)=d; m(1,2)=n;
        m(2,0)=n; m(2,1)=n; m(2,2)=d;
        return m;
    }();

    // Copy + scale in one step
    const double area = GetGeometry().Area();
    noalias(rMassMatrix) = unit_mass_matrix_2d * area;

    KRATOS_CATCH("");
}

[rubenzorrilla]

Fair enough. We will take care. Thanks for reporting 👍 👍 .