diff --git a/README.md b/README.md index 502974ce3..343fc3b9e 100644 --- a/README.md +++ b/README.md @@ -92,6 +92,7 @@ The following elements are supported on a triangle: - [Crouzeix-Raviart](https://defelement.com/elements/crouzeix-raviart.html) - [Bubble](https://defelement.com/elements/bubble.html) - [Hermite](https://defelement.com/elements/hermite.html) + - [iso](https://defelement.com/elements/p1-iso-p2.html) ### Quadrilateral @@ -132,6 +133,7 @@ The following elements are supported on a tetrahedron: - [Crouzeix-Raviart](https://defelement.com/elements/crouzeix-raviart.html) - [Bubble](https://defelement.com/elements/bubble.html) - [Hermite](https://defelement.com/elements/hermite.html) + - [iso](https://defelement.com/elements/p1-iso-p2.html) ### Hexahedron diff --git a/cpp/basix/polyset.cpp b/cpp/basix/polyset.cpp index 4a7d55851..a26f10ad5 100644 --- a/cpp/basix/polyset.cpp +++ b/cpp/basix/polyset.cpp @@ -931,6 +931,340 @@ void tabulate_polyset_triangle_macroedge_derivs( } //----------------------------------------------------------------------------- /// Compute the complete set of derivatives from 0 to nderiv, for all +/// the piecewise polynomials up to order n on a tetrahedron split into 8 by +/// splitting each edge into two parts +template +void tabulate_polyset_tetrahedron_macroedge_derivs( + stdex::mdspan> P, std::size_t n, + std::size_t nderiv, + stdex::mdspan> x) +{ + assert(x.extent(0) > 0); + assert(P.extent(0) == nderiv + 1); + assert(P.extent(1) == (n + 1) * (2 * n + 1) * (2 * n + 3) / 3); + assert(P.extent(2) == x.extent(0)); + + auto x0 = stdex::submdspan(x, stdex::full_extent, 0); + auto x1 = stdex::submdspan(x, stdex::full_extent, 1); + auto x2 = stdex::submdspan(x, stdex::full_extent, 2); + + std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); + + if (n == 0) + { + for (std::size_t p = 0; p < P.extent(2); ++p) + P(idx(0, 0), 0, p) = std::sqrt(6); + } + else if (n == 1) + { + for (std::size_t p = 0; p < P.extent(2); ++p) + { + if (x0[p] + x1[p] + x2[p] < 0.5) + { + P(idx(0, 0), 0, p) = 21.908902300206645 - 43.81780460041329 * x2[p] + - 43.81780460041329 * x1[p] + - 43.81780460041329 * x0[p]; + P(idx(0, 0), 4, p) = -5.855400437691199 + 11.710800875382398 * x2[p] + + 11.710800875382398 * x1[p] + + 35.1324026261472 * x0[p]; + P(idx(0, 0), 5, p) = -3.1326068447244277 + 6.2652136894488555 * x2[p] + + 25.75698961217863 * x1[p] + - 0.6961348543832062 * x0[p]; + P(idx(0, 0), 6, p) = -4.079436335011508 + 32.853642355761124 * x2[p] + + 3.359535805303594 * x1[p] + + 4.581185189050355 * x0[p]; + P(idx(0, 0), 7, p) = 1.6490105505038288 - 7.300270809207225 * x2[p] + - 8.666892660787566 * x1[p] + - 0.5229420350434975 * x0[p]; + P(idx(0, 0), 8, p) = 3.365373344712382 - 10.492441445989503 * x2[p] + - 11.258215021433038 * x1[p] + - 11.903076979701279 * x0[p]; + P(idx(0, 0), 9, p) = 0.8017837257372732 + 2.6726124191242437 * x2[p] + - 5.3452248382484875 * x1[p] + - 5.3452248382484875 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 0, p) = -43.81780460041329; + P(idx(0, 0, 1), 4, p) = 11.710800875382398; + P(idx(0, 0, 1), 5, p) = 6.2652136894488555; + P(idx(0, 0, 1), 6, p) = 32.853642355761124; + P(idx(0, 0, 1), 7, p) = -7.300270809207225; + P(idx(0, 0, 1), 8, p) = -10.492441445989503; + P(idx(0, 0, 1), 9, p) = 2.6726124191242437; + P(idx(0, 1, 0), 0, p) = -43.81780460041329; + P(idx(0, 1, 0), 4, p) = 11.710800875382398; + P(idx(0, 1, 0), 5, p) = 25.75698961217863; + P(idx(0, 1, 0), 6, p) = 3.359535805303594; + P(idx(0, 1, 0), 7, p) = -8.666892660787566; + P(idx(0, 1, 0), 8, p) = -11.258215021433038; + P(idx(0, 1, 0), 9, p) = -5.3452248382484875; + P(idx(1, 0, 0), 0, p) = -43.81780460041329; + P(idx(1, 0, 0), 4, p) = 35.1324026261472; + P(idx(1, 0, 0), 5, p) = -0.6961348543832062; + P(idx(1, 0, 0), 6, p) = 4.581185189050355; + P(idx(1, 0, 0), 7, p) = -0.5229420350434975; + P(idx(1, 0, 0), 8, p) = -11.903076979701279; + P(idx(1, 0, 0), 9, p) = -5.3452248382484875; + } + } + else if (x0[p] > 0.5) + { + P(idx(0, 0), 1, p) = -21.908902300206645 + 43.81780460041329 * x0[p]; + P(idx(0, 0), 4, p) = 29.277002188455995 - 23.421601750764797 * x2[p] + - 23.421601750764797 * x1[p] + - 35.1324026261472 * x0[p]; + P(idx(0, 0), 5, p) = -8.701685679790078 + 6.961348543832062 * x2[p] + + 6.961348543832062 * x1[p] + + 10.442022815748093 * x0[p]; + P(idx(0, 0), 6, p) = -4.472109351215823 + 3.5776874809726587 * x2[p] + + 3.5776874809726587 * x1[p] + + 5.366531221458988 * x0[p]; + P(idx(0, 0), 7, p) = 3.4688488324552 - 2.77507906596416 * x2[p] + - 2.77507906596416 * x1[p] + - 4.16261859894624 * x0[p]; + P(idx(0, 0), 8, p) = -1.148660363165304 + 26.439340288997876 * x2[p] + + 5.172330290276515 * x1[p] + - 2.8750095639459072 * x0[p]; + P(idx(0, 0), 9, p) = 0.8017837257372732 + 29.398736610366683 * x1[p] + - 5.3452248382484875 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 4, p) = -23.421601750764797; + P(idx(0, 0, 1), 5, p) = 6.961348543832062; + P(idx(0, 0, 1), 6, p) = 3.5776874809726587; + P(idx(0, 0, 1), 7, p) = -2.77507906596416; + P(idx(0, 0, 1), 8, p) = 26.439340288997876; + P(idx(0, 1, 0), 4, p) = -23.421601750764797; + P(idx(0, 1, 0), 5, p) = 6.961348543832062; + P(idx(0, 1, 0), 6, p) = 3.5776874809726587; + P(idx(0, 1, 0), 7, p) = -2.77507906596416; + P(idx(0, 1, 0), 8, p) = 5.172330290276515; + P(idx(0, 1, 0), 9, p) = 29.398736610366683; + P(idx(1, 0, 0), 1, p) = 43.81780460041329; + P(idx(1, 0, 0), 4, p) = -35.1324026261472; + P(idx(1, 0, 0), 5, p) = 10.442022815748093; + P(idx(1, 0, 0), 6, p) = 5.366531221458988; + P(idx(1, 0, 0), 7, p) = -4.16261859894624; + P(idx(1, 0, 0), 8, p) = -2.8750095639459072; + P(idx(1, 0, 0), 9, p) = -5.3452248382484875; + } + } + else if (x1[p] > 0.5) + { + P(idx(0, 0), 2, p) = -21.908902300206645 + 43.81780460041329 * x1[p]; + P(idx(0, 0), 5, p) = 24.36471990341222 - 19.491775922729772 * x2[p] + - 29.237663884094662 * x1[p] + - 19.491775922729772 * x0[p]; + P(idx(0, 0), 6, p) = -5.999171080899275 + 4.79933686471942 * x2[p] + + 7.199005297079131 * x1[p] + + 4.79933686471942 * x0[p]; + P(idx(0, 0), 7, p) = -0.4985380734081343 + 30.219077065046907 * x2[p] + - 4.371795412963639 * x1[p] + + 5.368871559779907 * x0[p]; + P(idx(0, 0), 8, p) = -6.952417987579472 - 0.6448619582682409 * x2[p] + + 9.377367643150668 * x1[p] + + 4.527468332008274 * x0[p]; + P(idx(0, 0), 9, p) = 0.8017837257372732 - 5.3452248382484875 * x1[p] + + 29.398736610366683 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 5, p) = -19.491775922729772; + P(idx(0, 0, 1), 6, p) = 4.79933686471942; + P(idx(0, 0, 1), 7, p) = 30.219077065046907; + P(idx(0, 0, 1), 8, p) = -0.6448619582682409; + P(idx(0, 1, 0), 2, p) = 43.81780460041329; + P(idx(0, 1, 0), 5, p) = -29.237663884094662; + P(idx(0, 1, 0), 6, p) = 7.199005297079131; + P(idx(0, 1, 0), 7, p) = -4.371795412963639; + P(idx(0, 1, 0), 8, p) = 9.377367643150668; + P(idx(0, 1, 0), 9, p) = -5.3452248382484875; + P(idx(1, 0, 0), 5, p) = -19.491775922729772; + P(idx(1, 0, 0), 6, p) = 4.79933686471942; + P(idx(1, 0, 0), 7, p) = 5.368871559779907; + P(idx(1, 0, 0), 8, p) = 4.527468332008274; + P(idx(1, 0, 0), 9, p) = 29.398736610366683; + } + } + else if (x2[p] > 0.5) + { + P(idx(0, 0), 3, p) = -21.908902300206645 + 43.81780460041329 * x2[p]; + P(idx(0, 0), 6, p) = 30.868462107172636 - 37.04215452860716 * x2[p] + - 24.69476968573811 * x1[p] + - 24.69476968573811 * x0[p]; + P(idx(0, 0), 7, p) = 1.209739241067291 - 6.421728190334149 * x2[p] + + 28.85245521346657 * x1[p] + + 4.002249708199567 * x0[p]; + P(idx(0, 0), 8, p) = -0.6784485185947118 - 2.4047977193753147 * x2[p] + - 1.4106355337117769 * x1[p] + + 25.0287047552861 * x0[p]; + P(idx(0, 0), 9, p) = 3.474396144861517 - 2.6726124191242437 * x2[p] + - 8.017837257372731 * x1[p] + - 8.017837257372731 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 3, p) = 43.81780460041329; + P(idx(0, 0, 1), 6, p) = -37.04215452860716; + P(idx(0, 0, 1), 7, p) = -6.421728190334149; + P(idx(0, 0, 1), 8, p) = -2.4047977193753147; + P(idx(0, 0, 1), 9, p) = -2.6726124191242437; + P(idx(0, 1, 0), 6, p) = -24.69476968573811; + P(idx(0, 1, 0), 7, p) = 28.85245521346657; + P(idx(0, 1, 0), 8, p) = -1.4106355337117769; + P(idx(0, 1, 0), 9, p) = -8.017837257372731; + P(idx(1, 0, 0), 6, p) = -24.69476968573811; + P(idx(1, 0, 0), 7, p) = 4.002249708199567; + P(idx(1, 0, 0), 8, p) = 25.0287047552861; + P(idx(1, 0, 0), 9, p) = -8.017837257372731; + } + } + else if (x1[p] + x2[p] < 0.5 && x0[p] + x1[p] < 0.5) + { + P(idx(0, 0), 4, p) = 11.710800875382398 - 23.421601750764797 * x2[p] + - 23.421601750764797 * x1[p]; + P(idx(0, 0), 5, p) = -3.480674271916031 + 6.961348543832062 * x2[p] + + 26.453124466561835 * x1[p]; + P(idx(0, 0), 6, p) = 10.558541102382724 + 3.5776874809726587 * x2[p] + - 25.91641906948487 * x1[p] + - 24.69476968573811 * x0[p]; + P(idx(0, 0), 7, p) = -0.6135853211177037 - 2.77507906596416 * x2[p] + - 4.1417009175445 * x1[p] + + 4.002249708199567 * x0[p]; + P(idx(0, 0), 8, p) = -15.100517522781306 + 26.439340288997876 * x2[p] + + 25.67356671355434 * x1[p] + + 25.0287047552861 * x0[p]; + P(idx(0, 0), 9, p) = 2.138089935299395 - 8.017837257372731 * x1[p] + - 8.017837257372731 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 4, p) = -23.421601750764797; + P(idx(0, 0, 1), 5, p) = 6.961348543832062; + P(idx(0, 0, 1), 6, p) = 3.5776874809726587; + P(idx(0, 0, 1), 7, p) = -2.77507906596416; + P(idx(0, 0, 1), 8, p) = 26.439340288997876; + P(idx(0, 1, 0), 4, p) = -23.421601750764797; + P(idx(0, 1, 0), 5, p) = 26.453124466561835; + P(idx(0, 1, 0), 6, p) = -25.91641906948487; + P(idx(0, 1, 0), 7, p) = -4.1417009175445; + P(idx(0, 1, 0), 8, p) = 25.67356671355434; + P(idx(0, 1, 0), 9, p) = -8.017837257372731; + P(idx(1, 0, 0), 6, p) = -24.69476968573811; + P(idx(1, 0, 0), 7, p) = 4.002249708199567; + P(idx(1, 0, 0), 8, p) = 25.0287047552861; + P(idx(1, 0, 0), 9, p) = -8.017837257372731; + } + } + else if (x1[p] + x2[p] < 0.5) + { + P(idx(0, 0), 4, p) = 11.710800875382398 - 23.421601750764797 * x2[p] + - 23.421601750764797 * x1[p]; + P(idx(0, 0), 5, p) = 6.2652136894488555 + 6.961348543832062 * x2[p] + + 6.961348543832062 * x1[p] + - 19.491775922729772 * x0[p]; + P(idx(0, 0), 6, p) = -4.18851217284604 + 3.5776874809726587 * x2[p] + + 3.5776874809726587 * x1[p] + + 4.79933686471942 * x0[p]; + P(idx(0, 0), 7, p) = -1.2968962469078738 - 2.77507906596416 * x2[p] + - 2.77507906596416 * x1[p] + + 5.368871559779907 * x0[p]; + P(idx(0, 0), 8, p) = -4.849899311142394 + 26.439340288997876 * x2[p] + + 5.172330290276515 * x1[p] + + 4.527468332008274 * x0[p]; + P(idx(0, 0), 9, p) = -16.57019699857031 + 29.398736610366683 * x1[p] + + 29.398736610366683 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 4, p) = -23.421601750764797; + P(idx(0, 0, 1), 5, p) = 6.961348543832062; + P(idx(0, 0, 1), 6, p) = 3.5776874809726587; + P(idx(0, 0, 1), 7, p) = -2.77507906596416; + P(idx(0, 0, 1), 8, p) = 26.439340288997876; + P(idx(0, 1, 0), 4, p) = -23.421601750764797; + P(idx(0, 1, 0), 5, p) = 6.961348543832062; + P(idx(0, 1, 0), 6, p) = 3.5776874809726587; + P(idx(0, 1, 0), 7, p) = -2.77507906596416; + P(idx(0, 1, 0), 8, p) = 5.172330290276515; + P(idx(0, 1, 0), 9, p) = 29.398736610366683; + P(idx(1, 0, 0), 5, p) = -19.491775922729772; + P(idx(1, 0, 0), 6, p) = 4.79933686471942; + P(idx(1, 0, 0), 7, p) = 5.368871559779907; + P(idx(1, 0, 0), 8, p) = 4.527468332008274; + P(idx(1, 0, 0), 9, p) = 29.398736610366683; + } + } + else if (x0[p] + x1[p] > 0.5) + { + P(idx(0, 0), 5, p) = 19.491775922729772 - 19.491775922729772 * x2[p] + - 19.491775922729772 * x1[p] + - 19.491775922729772 * x0[p]; + P(idx(0, 0), 6, p) = -4.79933686471942 + 4.79933686471942 * x2[p] + + 4.79933686471942 * x1[p] + + 4.79933686471942 * x0[p]; + P(idx(0, 0), 7, p) = -17.793974312413408 + 30.219077065046907 * x2[p] + + 30.219077065046907 * x1[p] + + 5.368871559779907 * x0[p]; + P(idx(0, 0), 8, p) = 8.692201812490664 - 0.6448619582682409 * x2[p] + - 21.9118719569896 * x1[p] + + 4.527468332008274 * x0[p]; + P(idx(0, 0), 9, p) = -16.57019699857031 + 29.398736610366683 * x1[p] + + 29.398736610366683 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 5, p) = -19.491775922729772; + P(idx(0, 0, 1), 6, p) = 4.79933686471942; + P(idx(0, 0, 1), 7, p) = 30.219077065046907; + P(idx(0, 0, 1), 8, p) = -0.6448619582682409; + P(idx(0, 1, 0), 5, p) = -19.491775922729772; + P(idx(0, 1, 0), 6, p) = 4.79933686471942; + P(idx(0, 1, 0), 7, p) = 30.219077065046907; + P(idx(0, 1, 0), 8, p) = -21.9118719569896; + P(idx(0, 1, 0), 9, p) = 29.398736610366683; + P(idx(1, 0, 0), 5, p) = -19.491775922729772; + P(idx(1, 0, 0), 6, p) = 4.79933686471942; + P(idx(1, 0, 0), 7, p) = 5.368871559779907; + P(idx(1, 0, 0), 8, p) = 4.527468332008274; + P(idx(1, 0, 0), 9, p) = 29.398736610366683; + } + } + else + { + P(idx(0, 0), 5, p) = 9.745887961364886 - 19.491775922729772 * x2[p]; + P(idx(0, 0), 6, p) = 9.947716410509344 + 4.79933686471942 * x2[p] + - 24.69476968573811 * x1[p] + - 24.69476968573811 * x0[p]; + P(idx(0, 0), 7, p) = -17.110663386623237 + 30.219077065046907 * x2[p] + + 28.85245521346657 * x1[p] + + 4.002249708199567 * x0[p]; + P(idx(0, 0), 8, p) = -1.5584163991482487 - 0.6448619582682409 * x2[p] + - 1.4106355337117769 * x1[p] + + 25.0287047552861 * x0[p]; + P(idx(0, 0), 9, p) = 2.138089935299395 - 8.017837257372731 * x1[p] + - 8.017837257372731 * x0[p]; + if (nderiv >= 1) + { + P(idx(0, 0, 1), 5, p) = -19.491775922729772; + P(idx(0, 0, 1), 6, p) = 4.79933686471942; + P(idx(0, 0, 1), 7, p) = 30.219077065046907; + P(idx(0, 0, 1), 8, p) = -0.6448619582682409; + P(idx(0, 1, 0), 6, p) = -24.69476968573811; + P(idx(0, 1, 0), 7, p) = 28.85245521346657; + P(idx(0, 1, 0), 8, p) = -1.4106355337117769; + P(idx(0, 1, 0), 9, p) = -8.017837257372731; + P(idx(1, 0, 0), 6, p) = -24.69476968573811; + P(idx(1, 0, 0), 7, p) = 4.002249708199567; + P(idx(1, 0, 0), 8, p) = 25.0287047552861; + P(idx(1, 0, 0), 9, p) = -8.017837257372731; + } + } + } + } + else + { + throw std::runtime_error("Only degree 0 and 1 macro polysets are currently " + "implemented on a tetrahedron."); + } +} +//----------------------------------------------------------------------------- +/// Compute the complete set of derivatives from 0 to nderiv, for all /// the piecewise polynomials up to order n on a hexahedron split into 4 by /// splitting each edge into two parts template @@ -2682,6 +3016,9 @@ void polyset::tabulate( case cell::type::triangle: tabulate_polyset_triangle_macroedge_derivs(P, d, n, x); return; + case cell::type::tetrahedron: + tabulate_polyset_tetrahedron_macroedge_derivs(P, d, n, x); + return; case cell::type::quadrilateral: tabulate_polyset_quadrilateral_macroedge_derivs(P, d, n, x); return; @@ -2760,6 +3097,8 @@ int polyset::dim(cell::type celltype, polyset::type ptype, int d) return 2 * d + 1; case cell::type::triangle: return (d + 1) * (2 * d + 1); + case cell::type::tetrahedron: + return (d + 1) * (2 * d + 1) * (2 * d + 3) / 3; case cell::type::quadrilateral: return (2 * d + 1) * (2 * d + 1); case cell::type::hexahedron: diff --git a/cpp/basix/quadrature.cpp b/cpp/basix/quadrature.cpp index 3b66745b5..5652c5aa5 100644 --- a/cpp/basix/quadrature.cpp +++ b/cpp/basix/quadrature.cpp @@ -4714,14 +4714,14 @@ make_macroedge_quadrature(quadrature::type rule, cell::type celltype, int m) { auto standard_q = quadrature::make_quadrature(rule, celltype, polyset::type::standard, m); + if (m == 0) + { + return standard_q; + } switch (celltype) { case cell::type::interval: { - if (m == 0) - { - return standard_q; - } const std::size_t npts = standard_q[0].size(); std::vector x(npts * 2); std::vector w(npts * 2); @@ -4736,10 +4736,6 @@ make_macroedge_quadrature(quadrature::type rule, cell::type celltype, int m) } case cell::type::triangle: { - if (m == 0) - { - return standard_q; - } const std::size_t npts = standard_q[0].size() / 2; std::vector x(npts * 8); std::vector w(npts * 4); @@ -4760,12 +4756,68 @@ make_macroedge_quadrature(quadrature::type rule, cell::type celltype, int m) } return {std::move(x), std::move(w)}; } - case cell::type::quadrilateral: + case cell::type::tetrahedron: { - if (m == 0) + const std::size_t npts = standard_q[0].size() / 3; + std::vector x(npts * 24); + std::vector w(npts * 8); + for (std::size_t i = 0; i < npts; ++i) { - return standard_q; + x[3 * i] = 0.5 * standard_q[0][3 * i]; + x[3 * i + 1] = 0.5 * standard_q[0][3 * i + 1]; + x[3 * i + 2] = 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + npts)] = 1.0 - 0.5 * standard_q[0][3 * i] + - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + npts) + 1] = 0.5 * standard_q[0][3 * i]; + x[3 * (i + npts) + 2] = 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 2 * npts)] = 0.5 * standard_q[0][3 * i]; + x[3 * (i + 2 * npts) + 1] = 1.0 - 0.5 * standard_q[0][3 * i] + - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 2 * npts) + 2] = 0.5 * standard_q[0][3 * i + 1]; + x[3 * (i + 3 * npts)] = 0.5 * standard_q[0][3 * i + 1]; + x[3 * (i + 3 * npts) + 1] = 0.5 * standard_q[0][3 * i]; + x[3 * (i + 3 * npts) + 2] = 1.0 - 0.5 * standard_q[0][3 * i] + - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 4 * npts)] = 0.5 - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 4 * npts) + 1] = 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 4 * npts) + 2] + = 0.5 - 0.5 * standard_q[0][3 * i] - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 5 * npts)] = 0.5 * standard_q[0][3 * i] + + 0.5 * standard_q[0][3 * i + 1] + + 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 5 * npts) + 1] = 0.5 - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 5 * npts) + 2] = 0.5 * standard_q[0][3 * i + 1]; + x[3 * (i + 6 * npts)] = 0.5 - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 6 * npts) + 1] = 0.5 * standard_q[0][3 * i] + + 0.5 * standard_q[0][3 * i + 1] + + 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 6 * npts) + 2] + = 0.5 - 0.5 * standard_q[0][3 * i] - 0.5 * standard_q[0][3 * i + 1]; + x[3 * (i + 7 * npts)] = 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 7 * npts) + 1] = 0.5 - 0.5 * standard_q[0][3 * i + 1] + - 0.5 * standard_q[0][3 * i + 2]; + x[3 * (i + 7 * npts) + 2] = 0.5 * standard_q[0][3 * i] + + 0.5 * standard_q[0][3 * i + 1] + + 0.5 * standard_q[0][3 * i + 2]; + w[i] = 0.125 * standard_q[1][i]; + w[npts + i] = 0.125 * standard_q[1][i]; + w[2 * npts + i] = 0.125 * standard_q[1][i]; + w[3 * npts + i] = 0.125 * standard_q[1][i]; + w[4 * npts + i] = 0.125 * standard_q[1][i]; + w[5 * npts + i] = 0.125 * standard_q[1][i]; + w[6 * npts + i] = 0.125 * standard_q[1][i]; + w[7 * npts + i] = 0.125 * standard_q[1][i]; } + return {std::move(x), std::move(w)}; + } + case cell::type::quadrilateral: + { const std::size_t npts = standard_q[0].size() / 2; std::vector x(npts * 8); std::vector w(npts * 4); @@ -4788,10 +4840,6 @@ make_macroedge_quadrature(quadrature::type rule, cell::type celltype, int m) } case cell::type::hexahedron: { - if (m == 0) - { - return standard_q; - } const std::size_t npts = standard_q[0].size() / 3; std::vector x(npts * 24); std::vector w(npts * 8); diff --git a/python/basix/finite_element.py b/python/basix/finite_element.py index 0f42109f3..8f1bfab2e 100644 --- a/python/basix/finite_element.py +++ b/python/basix/finite_element.py @@ -20,6 +20,7 @@ def string_to_family(family: str, cell: str) -> _EF: "P": _EF.P, "Bubble": _EF.bubble, "bubble": _EF.bubble, + "iso": _EF.iso, } # Family names that are valid on non-interval cells if cell != "interval": @@ -46,11 +47,6 @@ def string_to_family(family: str, cell: str) -> _EF: families.update({ "DPC": _EF.P, }) - # Cells that iso elements are implemented for - if cell in ["interval", "quadrilateral", "hexahedron", "triangle"]: - families.update({ - "iso": _EF.iso, - }) # Family names that are valid for tensor product cells if cell in ["interval", "quadrilateral", "hexahedron"]: families.update({ diff --git a/test/test_orthonormal_basis.py b/test/test_orthonormal_basis.py index 704e41c5b..db47f024c 100644 --- a/test/test_orthonormal_basis.py +++ b/test/test_orthonormal_basis.py @@ -121,6 +121,7 @@ def test_standard(cell_type, order): @pytest.mark.parametrize("cell_type", [ basix.CellType.interval, basix.CellType.triangle, + basix.CellType.tetrahedron, basix.CellType.quadrilateral, basix.CellType.hexahedron, ]) @@ -128,6 +129,8 @@ def test_standard(cell_type, order): def test_macroedge(cell_type, order): if cell_type == basix.CellType.triangle and order > 2: pytest.xfail("Degree > 2 edge macro polysets not implemented on triangles.") + if cell_type == basix.CellType.tetrahedron and order > 1: + pytest.xfail("Degree > 1 edge macro polysets not implemented on tetrahedra.") Qpts, Qwts = basix.make_quadrature(cell_type, 2*order + 1, polyset_type=basix.PolysetType.macroedge) basis = basix._basixcpp.tabulate_polynomial_set(