Add missing dispatches to _mass_qr

[henrij22]

This adds missing dispatches to _mass_qr() for DiscontinuousLagrange and Serendipity interpolations.

@KnutAM I am not super familiar with the other interpolations, if you give me a hint on how the integrationrules should look like, please let me know, then I can add them

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 94.14%. Comparing base (30185f8) to head (48b30cf).

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1200   +/-   ##
=======================================
  Coverage   94.14%   94.14%           
=======================================
  Files          40       40           
  Lines        6645     6646    +1     
=======================================
+ Hits         6256     6257    +1     
  Misses        389      389           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[KnutAM]

if you give me a hint on how the integrationrules should look like, please let me know, then I can add them

Basically, we have to look at the highest order polynomial term in $N(\xi)$ (let's call that $n$), and since the "mass" matrix is constructed as $M_{ij} = \int_\Omega N_i * N_j \mathrm{d}\Omega$, we need to integrate a polynomial of order $2n$.

I'm not exactly sure how the other cases have been derived (the special-casing of RefSimplex only for interpolation order 2 seems a bit strange to me).

But it could be tested by basically validating that that the local mass matrix becomes the same as for a higher order quadrature rule, i.e.

function _calculate_local_me(cv::CellValues)
    nbasef = getnbasefunctions(cv)
    Me = zeros(nbasef, nbasef)
    for q_point in 1:getnquadpoints(cv)
        dV = getdetJdV(cv, q_point)
        for i in 1:nbasef
        Ni = shape_value(cv, q_point, i)
            for j in 1:nbasef
                Nj = shape_value(cv, q_point, j)
                Me[i, j] += (Ni * Nj) * dV
            end
        end
    end
    return Me
end

We could then reformulate the Ferrite implementation, such that we have

_mass_qr(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip))
_mass_qr_order(::Lagrange{<:AbstractRefShape, order}) = order + 1
_mass_qr_order(::Lagrange{RefSimplex, 2}) = 4
_mass_qr_order(ip::VectorizedInterpolation) = _mass_qr_order(ip.ip)

To mirror the current implementation.

Then we could simply test these orders by defining
_mass_qr_test(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip) + 1) in the test suite, and calculate the Me above for both qr = Ferrite._mass_qr(ip) and qr_test = _mass_qr_test(ip).

[termi-official]

we need to integrate a polynomial of order 2n .

Please also note that this explanation is just a rule of thump, as the actual integral is over the reference element and you get an additional det(J) term which adds more potential nonlinearity to the integral.

I'm not exactly sure how the other cases have been derived

I can't find the issue, but the ,,order'' input to the quadrature rule constructor has different semantics, depending on the reference element.

[KnutAM]

additional det(J) term

Thanks, forgot about that one :D.

From that POV, I would even suggest we should change the docstring,

- a quadrature rule that integrates the mass-matrix exactly for the given interpolation `ip`
+ a quadrature rule that integrates the mass-matrix accurately for the given interpolation `ip`

? (Instead of writing that it is exactly integrated for a cell with the reference geometry / linear geometric interpolation).

And we run the tests using the reference geometry for the linear cell.

[KnutAM]

I can't find the issue, but the ,,order'' input to the quadrature rule constructor has different semantics, depending on the reference element.

Yes, to be clear I didn't mean that order in QuadratureRule should be 2n. I think there are some refs added in #1007 to find the translation between the quadrature rule order and the highest polynomial order that is exactly integrated.

Side-note - would be kindof cool to have this encoded, such that one could construct e.g.
qr = QuadratureRule{RefShape}([type::Symbol]; exact_polyorder) which would look up/calculate the correct order for the given quadrature rule, and call QuadratureRule{RefShape}([type], order)

[henrij22]

Hi, thanks for the input.
I implemented your test as suggested. For now its failing for

The mass matrix calculated with order+1 is not identical to the mass matrix calculated with order for Lagrange{RefTriangle, 3}()
The mass matrix calculated with order+1 is not identical to the mass matrix calculated with order for DiscontinuousLagrange{RefTriangle, 2}()
The mass matrix calculated with order+1 is not identical to the mass matrix calculated with order for DiscontinuousLagrange{RefTriangle, 3}()

I think DiscontinuousLagrange{RefTriangle, 2}() has to be 4 as well, but I am not so sure about order 3

[KnutAM]

I think DiscontinuousLagrange{RefTriangle, 2}() has to be 4 as well, but I am not so sure about order 3

Yes, this was the one I was suspicious of, why just order=2 is hard-coded, but orders > 2 are using the defaults.

For the RefTriangle, it looks like qr_order == polynomial order (ref + code) - can you confirm @termi-official ?

So I guess for RefTriangle, we should have = order^2, and probably for RefTetrahedron also (so for RefSimplex.

[termi-official]

For the RefTriangle, it looks like qr_order == polynomial order (ref + code) - can you confirm @termi-official ?

Yes, this is also what I always thought how it works. But we have this ambiguity in the codebase for a long time, that it is sometimes the quadature order and sometimes the order of the polynomial to be integrated.

[henrij22]

Lagrange{shape, 3}() fails with order^2 because it has 9 integration points, which is not implemented.

  ArgumentError: unsupported order for Dunavant's triangle integration

[henrij22]

From my understanding, shouldn't it be dim * order, so for RefTriangle 2*order and so on

[KnutAM]

The polynomial order is the order of the expression

$$M_{ij} = \int_\Omega N_i(\mathbf{x}) N_j(\mathbf{x}) \mathrm{d}\Omega \approx \sum_{q=1}^{N_\mathrm{quadpts}} N_i(\mathbf{x}_q) N_j(\mathbf{x}_q) \mathrm{det}(\mathbf{J}) w_q$$

So assuming the reference shape (\det(\boldsymbol{J}) is constant), the order is p^2 where p is the polynomial order of N(\boldsymbol{x}). So this is independent of the dimension of the problem. For RefSimplex (at least, probably all interpolations) the maximum polynomial order for Lagrange is the same as p.

Seeing that this fails, and was previously wrong. One option could be to choose to say that we only support automatic detection of the quadrature order for interpolations up to order 2? Thoughts @termi-official and @fredrikekre ?

The other option would be to implement a list of available quadrature (polynomial) orders for each reference shape, such that one can request the lowest order available quadrature rule that integrates a given polynomial order exactly. (Essentially my side-note above)

[KnutAM]

Sorry that there was no further feedback on this @henrij22. Now I've implemented the polyorder idea in #1220, I hope that can help this PR.

[termi-official]

I start to think that we should just implement adaptive quadrature rules for these use-cases here, as I do not see a general way to specify the order (and rule) here.

[termi-official]

Maybe we should add the geometry interpolation to the _mass_qr_order dispatch and at least warn the user if a nonlinear geometry is present that it might fail?

[KnutAM]

What do you mean with fail? Low accuracy or non-solvable system. I think only the former can happen due to nonlinear geometry, and most likely you need quite badly deformed elements?

[termi-official]

If you have a nonlinear geometry, e.g. a quadratic hexahedron, which is deformed with the second order modes, then integrating your detJ can become important. However, as detJ is now a fraction of polynomials the integration rules do not apply anymore. And I do not mean badly deformed elements here, just triggering the nonlinear modes can already cause some trouble. We can partially compensate this by increasing the integration order further.

[KnutAM]

Realized based on the discussion in #1216 (and related #623), that we currently have interpolations that are not polynomial (for pyramids). So we need to have keep a mechanism in the qr-order to deal with that too.

[henrij22]

Hey @KnutAM,
Coming back to this PR after some time. What is the status on #1220?
Otherwise I would suggest that we merge this here but only with some standard options, e.g.

_mass_qr(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip))
_mass_qr_order(::Lagrange{<:AbstractRefShape, order}) where {order} = order + 1
_mass_qr_order(::DiscontinuousLagrange{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}())
_mass_qr_order(::Serendipity{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}())

and lead the rest to the user for now. What do you think?

EDIT:
I think the above is then only valid for hypercubes

[KnutAM]

Sorry about this lingering for so long!
I was waiting for review on #1220, but I will push that soon and just merge it (unless there are protests), but as an internal feature for now before we agree on the public interface.