basix icon indicating copy to clipboard operation
basix copied to clipboard

Rational pyramidial base functions

Open jmv2009 opened this issue 1 year ago • 0 comments

Currently, the base functions for pyramids are regular polynomials: https://docs.fenicsproject.org/basix/main/cpp/namespacebasix_1_1polyset.html However, they need to be rational polynomials: https://defelement.org/elements/examples/pyramid-lagrange-equispaced-2.html https://www.math.u-bordeaux.fr/~durufle/montjoie/pyramid.php

The following relative minor changes should be performed: min(p,q) should be subtracted from both p+q terms in the polyset definition (both in the power of (1-z) and in the order of the Jacobi Polynomials. Alltogether, such terms then become max(p,q) instead of p+q. In practice, p+(q-min(p,q)) can be implemented, so the x0 polynomial is not affected.

This only changes the z dependency, and keeps othogonality of the rational polynomial. No changes to the quadrature weights are needed. The overall z-integrands of the orthogonality integral are still polynomial after integration over the horizontal directions.

At least in some cases, this is required for compatibility/admissibility at the two slanted faces of the reference pyramid, having the regular functions spaces there. This will enable many other elements as well, beyond the current Lagrange.

Corresponding changes in the derivatives and normalization are required.

We also need to change the documentation.

Suggested for lines 2266 of polyset.cpp, but further adjustments down the code are needed to accomodate.

if (p <= q) {
            for (std::size_t i = 0; i < r_pq.size(); ++i)
            {
              r_pq[i] = (0.5 + (x1[i] * 2.0 - 1.0) + (x2[i] * 2.0 - 1.0) * 0.5) / (1.0 - x2[i])
                        * _p[i] * (a + 1.0);
            }
}
else {
            for (std::size_t i = 0; i < r_pq.size(); ++i)
            {
              r_pq[i] = (0.5 + (x1[i] * 2.0 - 1.0) + (x2[i] * 2.0 - 1.0) * 0.5)
                        * _p[i] * (a + 1.0);
            }
]

or this may be fixed on the z axis recurrence relation.

This is inspired by Cockburn & Fu, https://arxiv.org/pdf/1605.00132, busy implementing the tnt elements. https://github.com/FEniCS/dolfinx/pull/3559

Pull request: https://github.com/FEniCS/basix/pull/895

jmv2009 avatar Jan 06 '25 18:01 jmv2009