Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rational pyramidial base functions #894

Open
jmv2009 opened this issue Jan 6, 2025 · 0 comments
Open

Rational pyramidial base functions #894

jmv2009 opened this issue Jan 6, 2025 · 0 comments
Assignees

Comments

@jmv2009
Copy link

jmv2009 commented Jan 6, 2025

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. FEniCS/dolfinx#3559

Pull request: #895

@jmv2009 jmv2009 changed the title Pyramidial base functions Fractional pyramidial base functions Jan 6, 2025
@jmv2009 jmv2009 changed the title Fractional pyramidial base functions Rational pyramidial base functions Jan 7, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants