Skip to content

Rational functions for pyramid #895

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.

Sign up for GitHub

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

Closed
wants to merge 10 commits into from
Closed

Rational functions for pyramid #895

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
0294f06
Rational functions for pyramid
jmv2009 Jan 9, 2025
fe781dd
Fractional functions with derivatives
jmv2009 Jan 11, 2025
5bd608a
Update polyset.cpp
jmv2009 Jan 11, 2025
b0d13f6
Update polyset.h doc update
jmv2009 Jan 11, 2025
8b7f6e6
Update polyset.cpp
jmv2009 Jan 12, 2025
04a647f
Update polyset.h
jmv2009 Jan 12, 2025
6ad5e9c
Update polyset.cpp
jmv2009 Feb 1, 2025
4385c7a
Update polyset.h
jmv2009 Feb 1, 2025
65dec06
Update cell.cpp
jmv2009 Feb 1, 2025
0b57159
Update quadrature.cpp
jmv2009 Feb 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion cpp/basix/cell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ cell::geometry(cell::type celltype)
0.0, 1.0, 0.0, 1.0, 1.0},
{6, 3}};
case cell::type::pyramid:
return {{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0,
return {{-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0,
0.0, 1.0},
{5, 3}};
case cell::type::hexahedron:
Expand Down
90 changes: 54 additions & 36 deletions cpp/basix/polyset.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2166,15 +2166,22 @@ void tabulate_polyset_pyramid_derivs(
// Traverse derivatives in increasing order
std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0);

for (std::size_t j = 0; j < P.extent(2); ++j)
P(idx(0, 0, 0), pyr_idx(0, 0, 0), j) = 1.0;

if (n == 0)
{
for (std::size_t j = 0; j < P.extent(2); ++j)
P(idx(0, 0, 0), pyr_idx(0, 0, 0), j) = std::sqrt(3);
return;
}
for (std::size_t j = 0; j < P.extent(2); ++j) {
//# P(idx(0, 0, 0), pyr_idx(0, 0, 0), j) = 1.0/(1.0-x2);
auto p00 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
P, idx(0, 0, 0), pyr_idx(0, 0, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);

for (std::size_t i = 0; i < p00.size(); ++i)
p00[i] = 1.0/(1.0-x2[i]);
}

// if (n == 0)
// {
// for (std::size_t j = 0; j < P.extent(2); ++j)
// P(idx(0, 0, 0), pyr_idx(0, 0, 0), j) = std::sqrt(3) / 2.0;
// return;
// }

for (std::size_t k = 0; k <= nderiv; ++k)
{
Expand All @@ -2198,8 +2205,7 @@ void tabulate_polyset_pyramid_derivs(
P, idx(kx, ky, kz), pyr_idx(p - 1, 0, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
for (std::size_t i = 0; i < p00.size(); ++i)
p00[i] = (0.5 + (x0[i] * 2.0 - 1.0) + (x2[i] * 2.0 - 1.0) * 0.5)
* p1[i] * (a + 1.0);
p00[i] = x0[i] * p1[i] * (a + 1.0);

if (kx > 0)
{
Expand All @@ -2208,17 +2214,7 @@ void tabulate_polyset_pyramid_derivs(
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);

for (std::size_t i = 0; i < p00.size(); ++i)
p00[i] += 2.0 * kx * p11[i] * (a + 1.0);
}

if (kz > 0)
{
auto pz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
P, idx(kx, ky, kz - 1), pyr_idx(p - 1, 0, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);

for (std::size_t i = 0; i < p00.size(); ++i)
p00[i] += kz * pz[i] * (a + 1.0);
p00[i] += kx * p11[i] * (a + 1.0);
}

if (p > 1)
Expand Down Expand Up @@ -2265,8 +2261,7 @@ void tabulate_polyset_pyramid_derivs(
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
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);
r_pq[i] = x1[i] * _p[i] * (a + 1.0);
}

if (ky > 0)
Expand All @@ -2275,18 +2270,24 @@ void tabulate_polyset_pyramid_derivs(
P, idx(kx, ky - 1, kz), pyr_idx(p, q - 1, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
for (std::size_t i = 0; i < r_pq.size(); ++i)
r_pq[i] += 2.0 * ky * _p[i] * (a + 1.0);
r_pq[i] += ky * _p[i] * (a + 1.0);
}

if (kz > 0)
if (kz > 0 and q <= p)
{
auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
P, idx(kx, ky, kz - 1), pyr_idx(p, q - 1, 0),
P, idx(kx, ky, kz - 1), pyr_idx(p, q, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
for (std::size_t i = 0; i < r_pq.size(); ++i)
r_pq[i] += kz * _p[i] * (a + 1.0);
r_pq[i] += kz * _p[i];
}


if (q <= p)
for (std::size_t i = 0; i < r_pq.size(); ++i)
if (x2[i] != 1.0)
{
r_pq[i] /= (1 - x2[i]);
}
if (q > 1)
{
auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
Expand All @@ -2295,7 +2296,12 @@ void tabulate_polyset_pyramid_derivs(
for (std::size_t i = 0; i < r_pq.size(); ++i)
{
const T f2 = 1.0 - x2[i];
r_pq[i] -= f2 * f2 * _p[i] * a;
if (q <= p)
r_pq[i] -=_p[i] * a;
else if (q <= p + 1)
r_pq[i] -= f2 * _p[i] * a;
else
r_pq[i] -= f2 * f2 * _p[i] * a;
}

if (kz > 0)
Expand All @@ -2304,10 +2310,22 @@ void tabulate_polyset_pyramid_derivs(
P, idx(kx, ky, kz - 1), pyr_idx(p, q - 2, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
for (std::size_t i = 0; i < r_pq.size(); ++i)
r_pq[i] += kz * (1.0 - (x2[i] * 2.0 - 1.0)) * _p[i] * a;
{
const T f2 = 1.0 - x2[i];
if (q <= p)
{
if (x2[i] != 1.0) r_pq[i] += 0.5 * kz * (1.0 - (x2[i] * 2.0 - 1.0)) * _p[i] * a /f2 /f2;
}
else if (q <= p + 1)
{
if (x2[i] != 1.0) r_pq[i] += 0.5 * kz * (1.0 - (x2[i] * 2.0 - 1.0)) * _p[i] * a /f2;
}
else
r_pq[i] += kz * (1.0 - (x2[i] * 2.0 - 1.0)) * _p[i] * a;
}
}

if (kz > 1)
if (kz > 1 and q > p + 1)
{
auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
P, idx(kx, ky, kz - 2), pyr_idx(p, q - 2, 0),
Expand Down Expand Up @@ -2335,7 +2353,7 @@ void tabulate_polyset_pyramid_derivs(
{
r_pq1[i]
= r_pq0[i]
* ((1.0 + p + q) + (x2[i] * 2.0 - 1.0) * (2.0 + p + q));
* (std::max(p, q) + (x2[i] * 2.0 - 1.0) * (1.0 + std::max(p, q)));
}

if (kz > 0)
Expand All @@ -2344,7 +2362,7 @@ void tabulate_polyset_pyramid_derivs(
P, idx(kx, ky, kz - 1), pyr_idx(p, q, 0),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
for (std::size_t i = 0; i < r_pq1.size(); ++i)
r_pq1[i] += 2 * kz * r_pq[i] * (2.0 + p + q);
r_pq1[i] += 2 * kz * r_pq[i] * (1.0 + std::max(p, q));
}
}
}
Expand All @@ -2355,7 +2373,7 @@ void tabulate_polyset_pyramid_derivs(
{
for (std::size_t q = 0; q < n - r; ++q)
{
auto [ar, br, cr] = jrc<T>(2 * p + 2 * q + 2, r);
auto [ar, br, cr] = jrc<T>(2 * std::max(p, q), r);
auto r_pqr = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
P, idx(kx, ky, kz), pyr_idx(p, q, r + 1),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
Expand Down Expand Up @@ -2398,7 +2416,7 @@ void tabulate_polyset_pyramid_derivs(
for (std::size_t i = 0; i < pqr.extent(0); ++i)
for (std::size_t j = 0; j < pqr.extent(1); ++j)
pqr(i, j)
*= std::sqrt(2 * (q + 0.5) * (p + 0.5) * (p + q + r + 1.5)) * 2;
*= std::sqrt(2 * (q + 0.5) * (p + 0.5) * (std::max(p, q) + r + 0.5));
}
}
}
Expand Down
14 changes: 11 additions & 3 deletions cpp/basix/polyset.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,13 +99,21 @@
/// \f$z\f$s.
///
/// ### Pyramid
/// Orthogonal polynomials on the pyramid element are best calculated in
/// Orthogonal rational functions on the pyramid element are best calculated in
/// the same way as the tetrahedron, using recurrence relations on each
/// axis. Let \f$\zeta_x = 2\frac{1+x}{1-z} - 1\f$, \f$\zeta_y =
/// 2\frac{1+y}{1-z} - 1\f$. The recurrence relation is then
/// 2\frac{1+y}{1-z} - 1\f$. The functions are then constituted as
///
/// \f[Q_{p, q, r} = P^{0,0}_p(\zeta_x) P^{0,0}_q(\zeta_y)
/// \left(\frac{1-z}{2}\right)^{(p+q)} P_r^{2(p+q+1), 0}(z).\f]
/// \left(1-z\right)^{max(p,q)-1} P_r^{2(max(p,q)+1), 0}(2z-1).\f]
///
/// Some nuance: An extra 1/(1-z) factor is included. To get the polynomials to span
/// the Lagrange spaces (0-Form), a factor of (1-z) need to be applied.
/// The orthogonal polynomials as defined here then include the derivatives of those polynomials,
/// as required to span the 1-Forms (Curl-spaces).
/// The derivatives provided are actually D[(1-z)Q, z]/(1-z),
/// which are the to be used anyway, after multiplication by (1-z),
/// so this was purposefully/lazily not corrected.
///
/// ### Normalisation
/// For each cell type, we obtain an orthonormal set of polynomials by
Expand Down
8 changes: 5 additions & 3 deletions cpp/basix/quadrature.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -419,9 +419,11 @@ std::array<std::vector<T>, 2> make_gauss_jacobi_quadrature(cell::type celltype,
for (std::size_t i = 0; i < x.extent(0); ++i)
{
const auto z = x(i, 2);
x(i, 0) *= (1 - z);
x(i, 1) *= (1 - z);
wts[i] *= (1 - z) * (1 - z);
x(i, 0) -= 0.5;
x(i, 1) -= 0.5;
x(i, 0) *= 2.0*(1 - z);
x(i, 1) *= 2.0*(1 - z);
wts[i] *= 4.0*(1 - z) * (1 - z);
}
return {std::move(pts), std::move(wts)};
}
Expand Down