From 98ad5c28dab7df9fa777f40bd0b5a83d18334683 Mon Sep 17 00:00:00 2001 From: Garth Wells Date: Fri, 14 Feb 2025 12:01:55 +0000 Subject: [PATCH 1/4] Return Basix hash for blocked elements --- python/basix/ufl.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/python/basix/ufl.py b/python/basix/ufl.py index 908cf456f..dfcc6293b 100644 --- a/python/basix/ufl.py +++ b/python/basix/ufl.py @@ -124,7 +124,12 @@ def __hash__(self) -> int: return hash("basix" + self._repr) def basix_hash(self) -> _typing.Optional[int]: - """Return the hash of the Basix element if this is a standard Basix element.""" + """Hash of the Basix element (if this is a standard Basix element). + + Returns: + Hash of the Basix element if this is a Basix element, + otherwise `None`. + """ return None @_abstractmethod @@ -319,7 +324,7 @@ def has_tensor_product_factorisation(self) -> bool: If this value is true, this element's basis functions can be computed as a tensor product of the basis elements of the - elements in the factoriaation. + elements in the factorisation. """ return False @@ -384,7 +389,7 @@ class _BasixElement(_ElementBase): This class allows elements created with `basix.create_element` to be wrapped as UFL compatible elements. Users should not directly call - this class's initiliser, but should use the `element` function + this class's initialiser, but should use the `element` function instead. """ @@ -1230,6 +1235,10 @@ def __hash__(self) -> int: """Return a hash.""" return super().__hash__() + def basix_hash(self) -> _typing.Optional[int]: + """Return the hash of the Basix element.""" + return self._sub_element.basix_hash() + @property def dtype(self) -> _npt.DTypeLike: """Element float type.""" From ab56c21f77d500436146a34915ca449dbe2a1045 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Fri, 14 Feb 2025 12:09:32 +0000 Subject: [PATCH 2/4] Add test --- test/test_ufl_wrapper.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_ufl_wrapper.py b/test/test_ufl_wrapper.py index 428b1789c..647269055 100644 --- a/test/test_ufl_wrapper.py +++ b/test/test_ufl_wrapper.py @@ -170,6 +170,21 @@ def test_finite_element_eq_hash(family, cell, degree, shape): assert (e1 == e2) == (hash(e1) == hash(e2)) +@pytest.mark.parametrize( + "family,cell,degree,shape", + [ + ("Lagrange", "triangle", 1, None), + ("Discontinuous Lagrange", "triangle", 1, None), + ("Lagrange", "quadrilateral", 1, None), + ("Lagrange", "triangle", 1, (2,)), + ("Lagrange", "triangle", 1, None), + ], +) +def test_finite_element_block_hash(family, cell, degree, shape): + e = basix.ufl.element(family, cell, degree, shape=shape) + assert e.basix_hash() is not None + + @pytest.mark.parametrize("component", [0, 1, 0]) def test_component_element_eq_hash(component): base_el = basix.ufl.element("Lagrange", "triangle", 1) From bc8b3e897b8025bfd23a29a405ac463ba6fa3cc3 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Fri, 14 Feb 2025 13:46:59 +0000 Subject: [PATCH 3/4] Simplify use of namespace --- cpp/basix/e-lagrange.cpp | 26 +- cpp/basix/e-regge.cpp | 9 +- cpp/basix/finite-element.cpp | 48 +-- cpp/basix/interpolation.cpp | 7 +- cpp/basix/lattice.cpp | 125 ++---- cpp/basix/moments.cpp | 11 +- cpp/basix/polynomials.cpp | 32 +- cpp/basix/polyset.cpp | 785 +++++++++++++---------------------- cpp/basix/quadrature.cpp | 43 +- 9 files changed, 393 insertions(+), 693 deletions(-) diff --git a/cpp/basix/e-lagrange.cpp b/cpp/basix/e-lagrange.cpp index 1ba9bdcaa..d1959081a 100644 --- a/cpp/basix/e-lagrange.cpp +++ b/cpp/basix/e-lagrange.cpp @@ -15,8 +15,8 @@ #include using namespace basix; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; namespace { @@ -448,9 +448,7 @@ basix::element::create_lagrange(cell::type celltype, int degree, std::array>, 4> x; std::array>, 4> M; x[0].emplace_back(1, 0); - M[0].emplace_back( - MDSPAN_IMPL_STANDARD_NAMESPACE::dextents{1, 1, 1, 1}, - 1); + M[0].emplace_back(md::dextents{1, 1, 1, 1}, 1); return FiniteElement( family::P, cell::type::point, polyset::type::standard, 0, {}, impl::mdspan_t(math::eye(1).data(), 1, 1), impl::to_mdspan(x), @@ -522,10 +520,7 @@ basix::element::create_lagrange(cell::type celltype, int degree, const auto [pt, shape] = lattice::create(celltype, 0, lattice_type, true, simplex_method); - x[tdim].emplace_back( - MDSPAN_IMPL_STANDARD_NAMESPACE::dextents{shape[0], - shape[1]}, - pt); + x[tdim].emplace_back(md::dextents{shape[0], shape[1]}, pt); auto& _M = M[tdim].emplace_back(shape[0], 1, shape[0], 1); std::fill(_M.data(), _M.data() + _M.size(), 0); for (std::size_t i = 0; i < shape[0]; ++i) @@ -543,10 +538,9 @@ basix::element::create_lagrange(cell::type celltype, int degree, = cell::sub_entity_geometry(celltype, dim, e); if (dim == 0) { - x[dim].emplace_back( - MDSPAN_IMPL_STANDARD_NAMESPACE::dextents{ - entity_x_shape[0], entity_x_shape[1]}, - entity_x); + x[dim].emplace_back(md::dextents{entity_x_shape[0], + entity_x_shape[1]}, + entity_x); auto& _M = M[dim].emplace_back(entity_x_shape[0], 1, entity_x_shape[0], 1); std::fill(_M.data(), _M.data() + _M.size(), 0); @@ -557,10 +551,8 @@ basix::element::create_lagrange(cell::type celltype, int degree, { const auto [pt, shape] = lattice::create( celltype, degree, lattice_type, false, simplex_method); - x[dim].emplace_back( - MDSPAN_IMPL_STANDARD_NAMESPACE::dextents{ - shape[0], shape[1]}, - pt); + x[dim].emplace_back(md::dextents{shape[0], shape[1]}, + pt); auto& _M = M[dim].emplace_back(shape[0], 1, shape[0], 1); std::fill(_M.data(), _M.data() + _M.size(), 0); for (std::size_t i = 0; i < shape[0]; ++i) diff --git a/cpp/basix/e-regge.cpp b/cpp/basix/e-regge.cpp index 61a94eb1d..56c1a34a2 100644 --- a/cpp/basix/e-regge.cpp +++ b/cpp/basix/e-regge.cpp @@ -13,8 +13,8 @@ #include using namespace basix; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; //----------------------------------------------------------------------------- template @@ -116,9 +116,8 @@ FiniteElement element::create_regge(cell::type celltype, int degree, // Store up outer(t, t) for all tangents const std::vector& vert_ids = topology[d][e]; const std::size_t ntangents = d * (d + 1) / 2; - stdex::mdarray> - vvt(ntangents, geometry.extent(1), geometry.extent(1)); + stdex::mdarray> vvt( + ntangents, geometry.extent(1), geometry.extent(1)); std::vector edge(geometry.extent(1)); int c = 0; diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp index 411a61fd2..225e608b8 100644 --- a/cpp/basix/finite-element.cpp +++ b/cpp/basix/finite-element.cpp @@ -28,16 +28,13 @@ #define str(X) str_macro(X) using namespace basix; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; template -using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; +using mdspan_t = md::mdspan>; template -using mdarray_t - = stdex::mdarray>; +using mdarray_t = stdex::mdarray>; namespace { @@ -551,14 +548,10 @@ element::make_discontinuous( std::array xshape = {npoints, tdim}; std::vector xb(xshape[0] * xshape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - new_x(xb.data(), xshape); + md::mdspan> new_x(xb.data(), xshape); std::array Mshape = {Mshape0, value_size, npoints, nderivs}; std::vector Mb(Mshape[0] * Mshape[1] * Mshape[2] * Mshape[3]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - new_M(Mb.data(), Mshape); + md::mdspan> new_M(Mb.data(), Mshape); int x_n = 0; int M_n = 0; for (int i = 0; i < 4; ++i) @@ -1520,9 +1513,7 @@ std::size_t FiniteElement::hash() const } std::size_t vs_hash = 0; for (std::size_t i = 0; i < value_shape().size(); ++i) - { combine_hashes(vs_hash, std::hash{}(value_shape()[i])); - } combine_hashes(h, coeff_hash); combine_hashes(h, std::hash{}(embedded_superdegree())); combine_hashes(h, std::hash{}(embedded_subdegree())); @@ -1530,9 +1521,8 @@ std::size_t FiniteElement::hash() const combine_hashes(h, vs_hash); } else - { combine_hashes(h, std::hash{}(degree())); - } + return h; } //----------------------------------------------------------------------------- @@ -1712,14 +1702,10 @@ FiniteElement::push_forward(impl::mdspan_t U, std::vector ub(shape[0] * shape[1] * shape[2]); mdspan_t u(ub.data(), shape); - using u_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; - using U_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; - using J_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; - using K_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using u_t = md::mdspan>; + using U_t = md::mdspan>; + using J_t = md::mdspan>; + using K_t = md::mdspan>; auto map = this->map_fn(); for (std::size_t i = 0; i < u.extent(0); ++i) { @@ -1752,14 +1738,10 @@ FiniteElement::pull_back(impl::mdspan_t u, std::vector Ub(shape[0] * shape[1] * shape[2]); mdspan_t U(Ub.data(), shape); - using u_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; - using U_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; - using J_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; - using K_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const F, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using u_t = md::mdspan>; + using U_t = md::mdspan>; + using J_t = md::mdspan>; + using K_t = md::mdspan>; auto map = this->map_fn(); for (std::size_t i = 0; i < u.extent(0); ++i) { diff --git a/cpp/basix/interpolation.cpp b/cpp/basix/interpolation.cpp index 9c8977376..85fa191f0 100644 --- a/cpp/basix/interpolation.cpp +++ b/cpp/basix/interpolation.cpp @@ -8,12 +8,11 @@ #include using namespace basix; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; template -using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; +using mdspan_t = md::mdspan>; //---------------------------------------------------------------------------- template diff --git a/cpp/basix/lattice.cpp b/cpp/basix/lattice.cpp index d66d9f995..e4329a856 100644 --- a/cpp/basix/lattice.cpp +++ b/cpp/basix/lattice.cpp @@ -15,8 +15,8 @@ #include using namespace basix; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; namespace { @@ -167,56 +167,48 @@ std::vector create_interval(std::size_t n, lattice::type lattice_type, } //----------------------------------------------------------------------------- template -stdex::mdarray> +stdex::mdarray> tabulate_dlagrange(std::size_t n, std::span x) { std::vector equi_pts(n + 1); for (std::size_t i = 0; i < equi_pts.size(); ++i) equi_pts[i] = static_cast(i) / static_cast(n); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - equi_pts_v(equi_pts.data(), n + 1, 1); + md::mdspan> equi_pts_v(equi_pts.data(), + n + 1, 1); const auto [dual_values_b, dshape] = polyset::tabulate( cell::type::interval, polyset::type::standard, n, 0, equi_pts_v); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - dual_values(dual_values_b.data(), dshape); + md::mdspan> dual_values( + dual_values_b.data(), dshape); std::vector dualmat_b(dual_values.extent(1) * dual_values.extent(2)); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - dualmat(dualmat_b.data(), dual_values.extent(1), dual_values.extent(2)); + md::mdspan> dualmat( + dualmat_b.data(), dual_values.extent(1), dual_values.extent(2)); for (std::size_t i = 0; i < dualmat.extent(0); ++i) for (std::size_t j = 0; j < dualmat.extent(1); ++j) dualmat(i, j) = dual_values(0, i, j); - using cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, - MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 1>>; + using cmdspan2_t + = md::mdspan>; const auto [tabulated_values_b, tshape] = polyset::tabulate(cell::type::interval, polyset::type::standard, n, 0, cmdspan2_t(x.data(), x.size(), 1)); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - tabulated_values(tabulated_values_b.data(), tshape); + md::mdspan> tabulated_values( + tabulated_values_b.data(), tshape); std::vector tabulated_b(tabulated_values.extent(1) * tabulated_values.extent(2)); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - tabulated(tabulated_b.data(), tabulated_values.extent(1), - tabulated_values.extent(2)); + md::mdspan> tabulated( + tabulated_b.data(), tabulated_values.extent(1), + tabulated_values.extent(2)); for (std::size_t i = 0; i < tabulated.extent(0); ++i) for (std::size_t j = 0; j < tabulated.extent(1); ++j) tabulated(i, j) = tabulated_values(0, i, j); std::vector c = math::solve(dualmat, tabulated); - return stdex::mdarray< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>( - tabulated.extents(), std::move(c)); + return stdex::mdarray>(tabulated.extents(), + std::move(c)); } //----------------------------------------------------------------------------- template @@ -227,8 +219,7 @@ std::vector warp_function(lattice::type lattice_type, int n, for (int i = 0; i < n + 1; ++i) pts[i] -= static_cast(i) / static_cast(n); - stdex::mdarray> v - = tabulate_dlagrange(n, x); + stdex::mdarray> v = tabulate_dlagrange(n, x); std::vector w(v.extent(1), 0); for (std::size_t i = 0; i < v.extent(0); ++i) for (std::size_t j = 0; j < v.extent(1); ++j) @@ -249,10 +240,8 @@ create_quad(std::size_t n, lattice::type lattice_type, bool exterior) const std::size_t m = r.size(); std::array shape = {m * m, 2}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 2>> - x(xb.data(), shape); + md::mdspan> x(xb.data(), + shape); std::size_t c = 0; for (std::size_t j = 0; j < m; ++j) { @@ -282,11 +271,8 @@ create_hex(int n, lattice::type lattice_type, bool exterior) std::array shape = {m * m * m, 3}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> + md::mdspan> x(xb.data(), m, m, m, 3); for (std::size_t k = 0; k < m; ++k) @@ -314,10 +300,8 @@ create_tri_equispaced(std::size_t n, bool exterior) std::array shape = {(n - 3 * b + 1) * (n - 3 * b + 2) / 2, 2}; std::vector _p(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 2>> - p(_p.data(), shape); + md::mdspan> p(_p.data(), + shape); // Displacement from GLL points in 1D, scaled by 1 /(r * (1 - r)) std::vector r = linspace(0.0, 1.0, 2 * n + 1); @@ -348,10 +332,8 @@ create_tri_warped(std::size_t n, lattice::type lattice_type, bool exterior) // Points std::array shape = {(n - 3 * b + 1) * (n - 3 * b + 2) / 2, 2}; std::vector _p(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 2>> - p(_p.data(), shape); + md::mdspan> p(_p.data(), + shape); // Displacement from GLL points in 1D, scaled by 1 /(r * (1 - r)) const std::vector r = linspace(0.0, 1.0, 2 * n + 1); @@ -421,10 +403,8 @@ create_tri_isaac(std::size_t n, lattice::type lattice_type, bool exterior) const std::size_t b = exterior ? 0 : 1; std::array shape = {(n - 3 * b + 1) * (n - 3 * b + 2) / 2, 2}; std::vector _p(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 2>> - p(_p.data(), shape); + md::mdspan> p(_p.data(), + shape); int c = 0; for (std::size_t j = b; j < (n - b + 1); ++j) @@ -457,10 +437,8 @@ create_tri_centroid(std::size_t n, lattice::type lattice_type, bool exterior) std::array shape = {(n - 2) * (n - 1) / 2, 2}; std::vector _p(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 2>> - p(_p.data(), shape); + md::mdspan> p(_p.data(), + shape); const std::vector x = create_interval(n, lattice_type, false); int c = 0; @@ -526,10 +504,8 @@ create_tet_equispaced(std::size_t n, bool exterior) std::array shape = {(n - 4 * b + 1) * (n - 4 * b + 2) * (n - 4 * b + 3) / 6, 3}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> - x(xb.data(), shape); + md::mdspan> x(xb.data(), + shape); std::size_t c = 0; for (std::size_t k = b; k < (n - b + 1); ++k) @@ -561,10 +537,8 @@ create_tet_isaac(std::size_t n, lattice::type lattice_type, bool exterior) std::array shape = {(n - 4 * b + 1) * (n - 4 * b + 2) * (n - 4 * b + 3) / 6, 3}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> - x(xb.data(), shape); + md::mdspan> x(xb.data(), + shape); int c = 0; for (std::size_t k = b; k < (n - b + 1); ++k) @@ -600,10 +574,8 @@ create_tet_warped(std::size_t n, lattice::type lattice_type, bool exterior) std::array shape = {(n - 4 * b + 1) * (n - 4 * b + 2) * (n - 4 * b + 3) / 6, 3}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> - p(xb.data(), shape); + md::mdspan> p(xb.data(), + shape); std::size_t c = 0; for (std::size_t k = b; k < (n - b + 1); ++k) @@ -658,10 +630,8 @@ create_tet_centroid(std::size_t n, lattice::type lattice_type, bool exterior) std::array shape = {(n - 3) * (n - 2) * (n - 1) / 6, 3}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> - p(xb.data(), shape); + md::mdspan> p(xb.data(), + shape); int c = 0; for (std::size_t i = 0; i + 2 < x.size(); ++i) @@ -726,13 +696,10 @@ std::pair, std::array> create_prism(std::size_t n, lattice::type lattice_type, bool exterior, lattice::simplex_method simplex_method) { - using cmdspan22_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, - MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 2>>; - using mdspan23_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>>; + using cmdspan22_t + = md::mdspan>; + using mdspan23_t + = md::mdspan>; if (n == 0) return {{1.0 / 3.0, 1.0 / 3.0, 0.5}, {1, 3}}; @@ -776,10 +743,8 @@ create_pyramid_equispaced(int n, bool exterior) std::array shape = {m, 3}; std::vector xb(shape[0] * shape[1]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< - std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> - x(xb.data(), shape); + md::mdspan> x(xb.data(), + shape); int c = 0; for (int k = 0; k < n + 1; ++k) { diff --git a/cpp/basix/moments.cpp b/cpp/basix/moments.cpp index 99819c8e3..0bcbcb613 100644 --- a/cpp/basix/moments.cpp +++ b/cpp/basix/moments.cpp @@ -9,18 +9,15 @@ #include "quadrature.h" using namespace basix; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; namespace { -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; template -using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; +using mdspan_t = md::mdspan>; template -using mdarray_t - = stdex::mdarray>; +using mdarray_t = stdex::mdarray>; //---------------------------------------------------------------------------- std::vector axis_points(const cell::type celltype) diff --git a/cpp/basix/polynomials.cpp b/cpp/basix/polynomials.cpp index ece917839..5064c9d00 100644 --- a/cpp/basix/polynomials.cpp +++ b/cpp/basix/polynomials.cpp @@ -11,18 +11,15 @@ #include using namespace basix; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; namespace { template -using mdarray_t - = stdex::mdarray>; +using mdarray_t = stdex::mdarray>; template -using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; +using mdspan_t = md::mdspan>; //----------------------------------------------------------------------------- constexpr int single_choose(int n, int k) @@ -114,11 +111,9 @@ tabulate_bernstein(cell::type celltype, int d, mdspan_t x) //----------------------------------------------------------------------------- template -std::pair, std::array> polynomials::tabulate( - polynomials::type polytype, cell::type celltype, int d, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) +std::pair, std::array> +polynomials::tabulate(polynomials::type polytype, cell::type celltype, int d, + md::mdspan> x) { switch (polytype) { @@ -146,15 +141,10 @@ int polynomials::dim(polynomials::type, cell::type cell, int d) //----------------------------------------------------------------------------- /// @cond template std::pair, std::array> -polynomials::tabulate( - polynomials::type, cell::type, int, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const float, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>); +polynomials::tabulate(polynomials::type, cell::type, int, + md::mdspan>); template std::pair, std::array> -polynomials::tabulate( - polynomials::type, cell::type, int, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const double, - MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>); +polynomials::tabulate(polynomials::type, cell::type, int, + md::mdspan>); /// @endcond //----------------------------------------------------------------------------- diff --git a/cpp/basix/polyset.cpp b/cpp/basix/polyset.cpp index 59d8e5576..6df1f57ec 100644 --- a/cpp/basix/polyset.cpp +++ b/cpp/basix/polyset.cpp @@ -6,16 +6,16 @@ #include "cell.h" #include "indexing.h" #include "mdspan.hpp" +#include #include #include #include -#include using namespace basix; using namespace basix::indexing; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; namespace { @@ -47,13 +47,8 @@ constexpr std::array jrc(int a, int n) // 1 and derivative 0. template void tabulate_polyset_point_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == nderiv + 1); @@ -73,13 +68,8 @@ void tabulate_polyset_point_derivs( /// 1]. The range is rescaled here to [0, 1]. template void tabulate_polyset_line_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == nderiv + 1); @@ -93,8 +83,7 @@ void tabulate_polyset_line_derivs( if (n == 0) return; - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); + auto x0 = md::submdspan(x, md::full_extent, 0); for (std::size_t i = 0; i < P.extent(2); ++i) P(0, 1, i) = (x0[i] * 2.0 - 1.0) * P(0, 0, i); @@ -139,21 +128,15 @@ void tabulate_polyset_line_derivs( /// parts. template void tabulate_polyset_line_macroedge_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == nderiv + 1); assert(P.extent(1) == 2 * n + 1); assert(P.extent(2) == x.extent(0)); - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); + auto x0 = md::submdspan(x, md::full_extent, 0); std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -247,13 +230,8 @@ void tabulate_polyset_line_macroedge_derivs( /// splitting each edge into two parts template void tabulate_polyset_quadrilateral_macroedge_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) / 2); @@ -264,10 +242,8 @@ void tabulate_polyset_quadrilateral_macroedge_derivs( auto quad_idx = [n](std::size_t px, std::size_t py) -> std::size_t { return (2 * n + 1) * px + py; }; - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -466,23 +442,16 @@ void tabulate_polyset_quadrilateral_macroedge_derivs( /// splitting each edge into two parts template void tabulate_polyset_triangle_macroedge_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) / 2); assert(P.extent(1) == (n + 1) * (2 * n + 1)); assert(P.extent(2) == x.extent(0)); - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -963,25 +932,17 @@ void tabulate_polyset_triangle_macroedge_derivs( /// splitting each edge into two parts template void tabulate_polyset_tetrahedron_macroedge_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) * (nderiv + 3) / 6); assert(P.extent(1) == (n + 1) * (2 * n + 1) * (2 * n + 3) / 3); assert(P.extent(2) == x.extent(0)); - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); - auto x2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 2); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); + auto x2 = md::submdspan(x, md::full_extent, 2); std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -1304,13 +1265,8 @@ void tabulate_polyset_tetrahedron_macroedge_derivs( /// splitting each edge into two parts template void tabulate_polyset_hexahedron_macroedge_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(0) > 0); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) * (nderiv + 3) / 6); @@ -1322,12 +1278,9 @@ void tabulate_polyset_hexahedron_macroedge_derivs( = [n](std::size_t px, std::size_t py, std::size_t pz) -> std::size_t { return (2 * n + 1) * (2 * n + 1) * px + (2 * n + 1) * py + pz; }; - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); - auto x2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 2); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); + auto x2 = md::submdspan(x, md::full_extent, 2); std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -1650,20 +1603,13 @@ void tabulate_polyset_hexahedron_macroedge_derivs( /// Karniadakis 1995 (https://doi.org/10.1016/0045-7825(94)00745-9). template void tabulate_polyset_triangle_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(1) == 2); - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) / 2); assert(P.extent(1) == (n + 1) * (n + 2) / 2); @@ -1688,12 +1634,8 @@ void tabulate_polyset_triangle_derivs( { for (std::size_t p = 1; p <= n; ++p) { - auto p0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0 = md::submdspan(P, idx(kx, ky), idx(0, p), md::full_extent); + auto p1 = md::submdspan(P, idx(kx, ky), idx(0, p - 1), md::full_extent); T a = static_cast(2 * p - 1) / static_cast(p); for (std::size_t i = 0; i < p0.extent(0); ++i) { @@ -1703,27 +1645,24 @@ void tabulate_polyset_triangle_derivs( if (kx > 0) { - auto px = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky), idx(0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto px = md::submdspan(P, idx(kx - 1, ky), idx(0, p - 1), + md::full_extent); for (std::size_t i = 0; i < p0.extent(0); ++i) p0[i] += 2 * kx * a * px[i]; } if (ky > 0) { - auto py = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1), idx(0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto py = md::submdspan(P, idx(kx, ky - 1), idx(0, p - 1), + md::full_extent); for (std::size_t i = 0; i < p0.extent(0); ++i) p0[i] += ky * a * py[i]; } if (p > 1) { - auto p2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p2 + = md::submdspan(P, idx(kx, ky), idx(0, p - 2), md::full_extent); // y^2 terms for (std::size_t i = 0; i < p0.extent(0); ++i) @@ -1734,18 +1673,16 @@ void tabulate_polyset_triangle_derivs( if (ky > 0) { - auto p2y = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1), idx(0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p2y = md::submdspan(P, idx(kx, ky - 1), idx(0, p - 2), + md::full_extent); for (std::size_t i = 0; i < p0.extent(0); ++i) p0[i] -= ky * ((x1[i] * 2.0 - 1.0) - 1.0) * p2y[i] * (a - 1.0); } if (ky > 1) { - auto p2y2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 2), idx(0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p2y2 = md::submdspan(P, idx(kx, ky - 2), idx(0, p - 2), + md::full_extent); for (std::size_t i = 0; i < p0.extent(0); ++i) p0[i] -= ky * (ky - 1) * p2y2[i] * (a - 1.0); } @@ -1754,20 +1691,15 @@ void tabulate_polyset_triangle_derivs( for (std::size_t p = 0; p < n; ++p) { - auto p0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0 = md::submdspan(P, idx(kx, ky), idx(0, p), md::full_extent); + auto p1 = md::submdspan(P, idx(kx, ky), idx(1, p), md::full_extent); for (std::size_t i = 0; i < p1.extent(0); ++i) p1[i] = p0[i] * ((x1[i] * 2.0 - 1.0) * (1.5 + p) + 0.5 + p); if (ky > 0) { - auto py = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1), idx(0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto py + = md::submdspan(P, idx(kx, ky - 1), idx(0, p), md::full_extent); for (std::size_t i = 0; i < p1.size(); ++i) p1[i] += 2 * ky * (1.5 + p) * py[i]; } @@ -1775,23 +1707,18 @@ void tabulate_polyset_triangle_derivs( for (std::size_t q = 1; q < n - p; ++q) { const auto [a1, a2, a3] = jrc(2 * p + 1, q); - auto pqp1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(q + 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pqm1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(q - 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pq = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), idx(q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqp1 + = md::submdspan(P, idx(kx, ky), idx(q + 1, p), md::full_extent); + auto pqm1 + = md::submdspan(P, idx(kx, ky), idx(q - 1, p), md::full_extent); + auto pq = md::submdspan(P, idx(kx, ky), idx(q, p), md::full_extent); for (std::size_t i = 0; i < pqp1.extent(0); ++i) pqp1[i] = pq[i] * ((x1[i] * 2.0 - 1.0) * a1 + a2) - pqm1[i] * a3; if (ky > 0) { - auto py = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1), idx(q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto py + = md::submdspan(P, idx(kx, ky - 1), idx(q, p), md::full_extent); for (std::size_t i = 0; i < pqp1.extent(0); ++i) pqp1[i] += 2 * ky * a1 * py[i]; } @@ -1820,25 +1747,17 @@ void tabulate_polyset_triangle_derivs( //----------------------------------------------------------------------------- template void tabulate_polyset_tetrahedron_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(1) == 3); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) * (nderiv + 3) / 6); assert(P.extent(1) == (n + 1) * (n + 2) * (n + 3) / 6); assert(P.extent(2) == x.extent(0)); - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); - auto x2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 2); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); + auto x2 = md::submdspan(x, md::full_extent, 2); // Traverse derivatives in increasing order std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -1860,12 +1779,10 @@ void tabulate_polyset_tetrahedron_derivs( { for (std::size_t p = 1; p <= n; ++p) { - auto p00 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, 0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p0m1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, 0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p00 = md::submdspan(P, idx(kx, ky, kz), idx(0, 0, p), + md::full_extent); + auto p0m1 = md::submdspan(P, idx(kx, ky, kz), idx(0, 0, p - 1), + md::full_extent); T a = static_cast(2 * p - 1) / static_cast(p); for (std::size_t i = 0; i < p00.size(); ++i) { @@ -1876,36 +1793,32 @@ void tabulate_polyset_tetrahedron_derivs( if (kx > 0) { - auto p0m1x = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky, kz), idx(0, 0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m1x = md::submdspan(P, idx(kx - 1, ky, kz), idx(0, 0, p - 1), + md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] += 2 * kx * a * p0m1x[i]; } if (ky > 0) { - auto p0m1y = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, kz), idx(0, 0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m1y = md::submdspan(P, idx(kx, ky - 1, kz), idx(0, 0, p - 1), + md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] += ky * a * p0m1y[i]; } if (kz > 0) { - auto p0m1z = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(0, 0, p - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m1z = md::submdspan(P, idx(kx, ky, kz - 1), idx(0, 0, p - 1), + md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] += kz * a * p0m1z[i]; } if (p > 1) { - auto p0m2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, 0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m2 = md::submdspan(P, idx(kx, ky, kz), idx(0, 0, p - 2), + md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) { T f2 = x1[i] + x2[i] - 1.0; @@ -1913,9 +1826,8 @@ void tabulate_polyset_tetrahedron_derivs( } if (ky > 0) { - auto p0m2y = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, kz), idx(0, 0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m2y = md::submdspan(P, idx(kx, ky - 1, kz), + idx(0, 0, p - 2), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) { p00[i] -= ky * ((x1[i] * 2.0 - 1.0) + (x2[i] * 2.0 - 1.0)) @@ -1925,18 +1837,16 @@ void tabulate_polyset_tetrahedron_derivs( if (ky > 1) { - auto p0m2y2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 2, kz), idx(0, 0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m2y2 = md::submdspan(P, idx(kx, ky - 2, kz), + idx(0, 0, p - 2), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] -= ky * (ky - 1) * p0m2y2[i] * (a - 1.0); } if (kz > 0) { - auto p0m2z = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(0, 0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m2z = md::submdspan(P, idx(kx, ky, kz - 1), + idx(0, 0, p - 2), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] -= kz * ((x1[i] * 2.0 - 1.0) + (x2[i] * 2.0 - 1.0)) * p0m2z[i] * (a - 1.0); @@ -1944,18 +1854,16 @@ void tabulate_polyset_tetrahedron_derivs( if (kz > 1) { - auto p0m2z2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 2), idx(0, 0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m2z2 = md::submdspan(P, idx(kx, ky, kz - 2), + idx(0, 0, p - 2), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] -= kz * (kz - 1) * p0m2z2[i] * (a - 1.0); } if (ky > 0 and kz > 0) { - auto p0m2yz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, kz - 1), idx(0, 0, p - 2), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0m2yz = md::submdspan(P, idx(kx, ky - 1, kz - 1), + idx(0, 0, p - 2), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] -= 2.0 * ky * kz * p0m2yz[i] * (a - 1.0); } @@ -1964,12 +1872,10 @@ void tabulate_polyset_tetrahedron_derivs( for (std::size_t p = 0; p < n; ++p) { - auto p10 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p00 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, 0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p10 = md::submdspan(P, idx(kx, ky, kz), idx(0, 1, p), + md::full_extent); + auto p00 = md::submdspan(P, idx(kx, ky, kz), idx(0, 0, p), + md::full_extent); for (std::size_t i = 0; i < p10.size(); ++i) p10[i] = p00[i] @@ -1978,18 +1884,16 @@ void tabulate_polyset_tetrahedron_derivs( * 0.5); if (ky > 0) { - auto p0y = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, kz), idx(0, 0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0y = md::submdspan(P, idx(kx, ky - 1, kz), idx(0, 0, p), + md::full_extent); for (std::size_t i = 0; i < p10.size(); ++i) p10[i] += 2 * ky * p0y[i] * (1.5 + p); } if (kz > 0) { - auto p0z = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(0, 0, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0z = md::submdspan(P, idx(kx, ky, kz - 1), idx(0, 0, p), + md::full_extent); for (std::size_t i = 0; i < p10.size(); ++i) p10[i] += kz * p0z[i]; } @@ -1997,15 +1901,12 @@ void tabulate_polyset_tetrahedron_derivs( for (std::size_t q = 1; q < n - p; ++q) { auto [aq, bq, cq] = jrc(2 * p + 1, q); - auto pq1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, q + 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pq = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pqm1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, q - 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pq1 = md::submdspan(P, idx(kx, ky, kz), idx(0, q + 1, p), + md::full_extent); + auto pq = md::submdspan(P, idx(kx, ky, kz), idx(0, q, p), + md::full_extent); + auto pqm1 = md::submdspan(P, idx(kx, ky, kz), idx(0, q - 1, p), + md::full_extent); for (std::size_t i = 0; i < pq1.size(); ++i) { T f4 = 1.0 - x2[i]; @@ -2014,21 +1915,18 @@ void tabulate_polyset_tetrahedron_derivs( } if (ky > 0) { - auto pqy = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, kz), idx(0, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqy = md::submdspan(P, idx(kx, ky - 1, kz), idx(0, q, p), + md::full_extent); for (std::size_t i = 0; i < pq1.size(); ++i) pq1[i] += 2 * ky * pqy[i] * aq; } if (kz > 0) { - auto pqz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(0, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pq1z = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(0, q - 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqz = md::submdspan(P, idx(kx, ky, kz - 1), idx(0, q, p), + md::full_extent); + auto pq1z = md::submdspan(P, idx(kx, ky, kz - 1), + idx(0, q - 1, p), md::full_extent); for (std::size_t i = 0; i < pq1.size(); ++i) { pq1[i] += kz * pqz[i] * (aq - bq) @@ -2038,9 +1936,8 @@ void tabulate_polyset_tetrahedron_derivs( if (kz > 1) { - auto pq1z2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 2), idx(0, q - 1, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pq1z2 = md::submdspan(P, idx(kx, ky, kz - 2), + idx(0, q - 1, p), md::full_extent); // Quadratic term in z for (std::size_t i = 0; i < pq1.size(); ++i) pq1[i] -= kz * (kz - 1) * pq1z2[i] * cq; @@ -2052,12 +1949,10 @@ void tabulate_polyset_tetrahedron_derivs( { for (std::size_t q = 0; q < n - p; ++q) { - auto pq = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(1, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pq0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(0, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pq = md::submdspan(P, idx(kx, ky, kz), idx(1, q, p), + md::full_extent); + auto pq0 = md::submdspan(P, idx(kx, ky, kz), idx(0, q, p), + md::full_extent); for (std::size_t i = 0; i < pq.size(); ++i) { pq[i] = pq0[i] @@ -2066,9 +1961,8 @@ void tabulate_polyset_tetrahedron_derivs( if (kz > 0) { - auto pqz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(0, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqz = md::submdspan(P, idx(kx, ky, kz - 1), idx(0, q, p), + md::full_extent); for (std::size_t i = 0; i < pq.size(); ++i) pq[i] += 2 * kz * (2.0 + p + q) * pqz[i]; } @@ -2082,15 +1976,12 @@ void tabulate_polyset_tetrahedron_derivs( for (std::size_t r = 1; r < n - p - q; ++r) { auto [ar, br, cr] = jrc(2 * p + 2 * q + 2, r); - auto pqr1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(r + 1, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pqr = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(r, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pqrm1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), idx(r - 1, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqr1 = md::submdspan(P, idx(kx, ky, kz), idx(r + 1, q, p), + md::full_extent); + auto pqr = md::submdspan(P, idx(kx, ky, kz), idx(r, q, p), + md::full_extent); + auto pqrm1 = md::submdspan(P, idx(kx, ky, kz), idx(r - 1, q, p), + md::full_extent); for (std::size_t i = 0; i < pqr1.size(); ++i) { @@ -2100,9 +1991,8 @@ void tabulate_polyset_tetrahedron_derivs( if (kz > 0) { - auto pqrz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), idx(r, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqrz = md::submdspan(P, idx(kx, ky, kz - 1), idx(r, q, p), + md::full_extent); for (std::size_t i = 0; i < pqr1.size(); ++i) pqr1[i] += 2 * kz * ar * pqrz[i]; } @@ -2120,9 +2010,8 @@ void tabulate_polyset_tetrahedron_derivs( { for (std::size_t r = 0; r <= n - p - q; ++r) { - auto pqr = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, idx(r, q, p), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqr + = md::submdspan(P, md::full_extent, idx(r, q, p), md::full_extent); 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) @@ -2135,13 +2024,8 @@ void tabulate_polyset_tetrahedron_derivs( //----------------------------------------------------------------------------- template void tabulate_polyset_pyramid_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(1) == 3); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) * (nderiv + 3) / 6); @@ -2156,12 +2040,9 @@ void tabulate_polyset_pyramid_derivs( return r0 + p * rv + q; }; - const auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - const auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); - const auto x2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 2); + const auto x0 = md::submdspan(x, md::full_extent, 0); + const auto x1 = md::submdspan(x, md::full_extent, 1); + const auto x2 = md::submdspan(x, md::full_extent, 2); // Traverse derivatives in increasing order std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0); @@ -2191,21 +2072,18 @@ void tabulate_polyset_pyramid_derivs( if (p > 0) { const T a = static_cast(p - 1) / static_cast(p); - auto p00 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p - 1, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p00 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, 0, 0), + md::full_extent); + auto p1 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p - 1, 0, 0), + md::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); if (kx > 0) { - auto p11 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky, kz), pyr_idx(p - 1, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p11 = md::submdspan(P, idx(kx - 1, ky, kz), + pyr_idx(p - 1, 0, 0), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] += 2.0 * kx * p11[i] * (a + 1.0); @@ -2213,9 +2091,8 @@ void tabulate_polyset_pyramid_derivs( 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); + auto pz = md::submdspan(P, idx(kx, ky, kz - 1), + pyr_idx(p - 1, 0, 0), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] += kz * pz[i] * (a + 1.0); @@ -2223,9 +2100,8 @@ void tabulate_polyset_pyramid_derivs( if (p > 1) { - auto p2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p - 2, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p2 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p - 2, 0, 0), + md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) { T f2 = 1.0 - x2[i]; @@ -2234,9 +2110,8 @@ void tabulate_polyset_pyramid_derivs( if (kz > 0) { - auto p2z = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), pyr_idx(p - 2, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p2z = md::submdspan(P, idx(kx, ky, kz - 1), + pyr_idx(p - 2, 0, 0), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] += kz * (1.0 - (x2[i] * 2.0 - 1.0)) * p2z[i] * a; } @@ -2244,9 +2119,8 @@ void tabulate_polyset_pyramid_derivs( if (kz > 1) { // quadratic term in z - auto pz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 2), pyr_idx(p - 2, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pz = md::submdspan(P, idx(kx, ky, kz - 2), + pyr_idx(p - 2, 0, 0), md::full_extent); for (std::size_t i = 0; i < p00.size(); ++i) p00[i] -= kz * (kz - 1) * pz[i] * a; } @@ -2256,13 +2130,11 @@ void tabulate_polyset_pyramid_derivs( for (std::size_t q = 1; q < n + 1; ++q) { const T a = static_cast(q - 1) / static_cast(q); - auto r_pq = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto r_pq = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q, 0), + md::full_extent); - auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q - 1, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _p = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q - 1, 0), + md::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) @@ -2271,27 +2143,24 @@ void tabulate_polyset_pyramid_derivs( if (ky > 0) { - auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, kz), pyr_idx(p, q - 1, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _p = md::submdspan(P, idx(kx, ky - 1, kz), + pyr_idx(p, q - 1, 0), md::full_extent); for (std::size_t i = 0; i < r_pq.size(); ++i) r_pq[i] += 2.0 * ky * _p[i] * (a + 1.0); } if (kz > 0) { - auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), pyr_idx(p, q - 1, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _p = md::submdspan(P, idx(kx, ky, kz - 1), + pyr_idx(p, q - 1, 0), md::full_extent); for (std::size_t i = 0; i < r_pq.size(); ++i) r_pq[i] += kz * _p[i] * (a + 1.0); } if (q > 1) { - auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q - 2, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _p = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q - 2, 0), + md::full_extent); for (std::size_t i = 0; i < r_pq.size(); ++i) { const T f2 = 1.0 - x2[i]; @@ -2300,18 +2169,16 @@ void tabulate_polyset_pyramid_derivs( if (kz > 0) { - auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), pyr_idx(p, q - 2, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _p = md::submdspan(P, idx(kx, ky, kz - 1), + pyr_idx(p, q - 2, 0), md::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; } if (kz > 1) { - auto _p = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 2), pyr_idx(p, q - 2, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _p = md::submdspan(P, idx(kx, ky, kz - 2), + pyr_idx(p, q - 2, 0), md::full_extent); for (std::size_t i = 0; i < r_pq.size(); ++i) r_pq[i] -= kz * (kz - 1) * _p[i] * a; } @@ -2324,13 +2191,11 @@ void tabulate_polyset_pyramid_derivs( { for (std::size_t q = 0; q < n; ++q) { - auto r_pq1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q, 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto r_pq1 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q, 1), + md::full_extent); - auto r_pq0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto r_pq0 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q, 0), + md::full_extent); for (std::size_t i = 0; i < r_pq1.size(); ++i) { r_pq1[i] @@ -2340,9 +2205,8 @@ void tabulate_polyset_pyramid_derivs( if (kz > 0) { - auto r_pq = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), pyr_idx(p, q, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto r_pq = md::submdspan(P, idx(kx, ky, kz - 1), + pyr_idx(p, q, 0), md::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); } @@ -2356,15 +2220,12 @@ void tabulate_polyset_pyramid_derivs( for (std::size_t q = 0; q < n - r; ++q) { auto [ar, br, cr] = jrc(2 * p + 2 * q + 2, 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); - auto _r0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q, r), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto _r1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), pyr_idx(p, q, r - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto r_pqr = md::submdspan(P, idx(kx, ky, kz), + pyr_idx(p, q, r + 1), md::full_extent); + auto _r0 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q, r), + md::full_extent); + auto _r1 = md::submdspan(P, idx(kx, ky, kz), pyr_idx(p, q, r - 1), + md::full_extent); for (std::size_t i = 0; i < r_pqr.size(); ++i) { r_pqr[i] @@ -2373,9 +2234,8 @@ void tabulate_polyset_pyramid_derivs( if (kz > 0) { - auto _r = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), pyr_idx(p, q, r), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _r = md::submdspan(P, idx(kx, ky, kz - 1), + pyr_idx(p, q, r), md::full_extent); for (std::size_t i = 0; i < r_pqr.size(); ++i) r_pqr[i] += ar * 2 * kz * _r[i]; } @@ -2392,9 +2252,8 @@ void tabulate_polyset_pyramid_derivs( { for (std::size_t q = 0; q <= n - r; ++q) { - auto pqr = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, pyr_idx(p, q, r), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqr = md::submdspan(P, md::full_extent, pyr_idx(p, q, r), + md::full_extent); 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) @@ -2406,13 +2265,8 @@ void tabulate_polyset_pyramid_derivs( //----------------------------------------------------------------------------- template void tabulate_polyset_quad_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(1) == 2); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) / 2); @@ -2424,10 +2278,8 @@ void tabulate_polyset_quad_derivs( { return (n + 1) * px + py; }; // Compute 1D basis - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); assert(x0.extent(0) > 0); assert(x1.extent(0) > 0); @@ -2441,9 +2293,7 @@ void tabulate_polyset_quad_derivs( return; { // scope - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, 0), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result = md::submdspan(P, idx(0, 0), md::full_extent, md::full_extent); for (std::size_t i = 0; i < result.extent(1); ++i) { result(quad_idx(0, 1), i) @@ -2464,12 +2314,10 @@ void tabulate_polyset_quad_derivs( for (std::size_t ky = 1; ky <= nderiv; ++ky) { // Get reference to this derivative - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky - 1), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, ky), md::full_extent, md::full_extent); + auto result0 + = md::submdspan(P, idx(0, ky - 1), md::full_extent, md::full_extent); for (std::size_t i = 0; i < result.extent(1); ++i) { result(quad_idx(0, 1), i) @@ -2493,9 +2341,8 @@ void tabulate_polyset_quad_derivs( // Take tensor product with another interval for (std::size_t ky = 0; ky <= nderiv; ++ky) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, ky), md::full_extent, md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2511,9 +2358,8 @@ void tabulate_polyset_quad_derivs( const T a = 1.0 - 1.0 / static_cast(px); for (std::size_t ky = 0; ky <= nderiv; ++ky) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, ky), md::full_extent, md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2531,12 +2377,10 @@ void tabulate_polyset_quad_derivs( { for (std::size_t ky = 0; ky <= nderiv - kx; ++ky) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(kx, ky), md::full_extent, md::full_extent); + auto result0 + = md::submdspan(P, idx(kx - 1, ky), md::full_extent, md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2553,12 +2397,10 @@ void tabulate_polyset_quad_derivs( const T a = 1.0 - 1.0 / static_cast(px); for (std::size_t ky = 0; ky <= nderiv - kx; ++ky) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(kx, ky), md::full_extent, md::full_extent); + auto result0 = md::submdspan(P, idx(kx - 1, ky), md::full_extent, + md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2579,9 +2421,8 @@ void tabulate_polyset_quad_derivs( { for (std::size_t py = 0; py <= n; ++py) { - auto pxy = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, quad_idx(px, py), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pxy = md::submdspan(P, md::full_extent, quad_idx(px, py), + md::full_extent); for (std::size_t i = 0; i < pxy.extent(0); ++i) for (std::size_t j = 0; j < pxy.extent(1); ++j) pxy(i, j) *= std::sqrt((2 * px + 1) * (2 * py + 1)); @@ -2591,13 +2432,8 @@ void tabulate_polyset_quad_derivs( //----------------------------------------------------------------------------- template void tabulate_polyset_hex_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(1) == 3); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) * (nderiv + 3) / 6); @@ -2610,12 +2446,9 @@ void tabulate_polyset_hex_derivs( { return (n + 1) * (n + 1) * px + (n + 1) * py + pz; }; // Compute 1D basis - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); - auto x2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 2); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); + auto x2 = md::submdspan(x, md::full_extent, 2); assert(x0.extent(0) > 0); assert(x1.extent(0) > 0); assert(x2.extent(0) > 0); @@ -2630,9 +2463,8 @@ void tabulate_polyset_hex_derivs( // Tabulate interval for px=py=0 // For kz = 0 { // scope - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, 0, 0), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, 0, 0), md::full_extent, md::full_extent); // for pz = 1 for (std::size_t i = 0; i < result.extent(1); ++i) { @@ -2657,12 +2489,10 @@ void tabulate_polyset_hex_derivs( for (std::size_t kz = 1; kz <= nderiv; ++kz) { // Get reference to this derivative - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, 0, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, 0, kz - 1), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, 0, kz), md::full_extent, md::full_extent); + auto result0 + = md::submdspan(P, idx(0, 0, kz - 1), md::full_extent, md::full_extent); // for pz = 1 for (std::size_t i = 0; i < result.extent(1); ++i) { @@ -2690,9 +2520,8 @@ void tabulate_polyset_hex_derivs( // for py = 1 for (std::size_t kz = 0; kz <= nderiv; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, 0, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, 0, kz), md::full_extent, md::full_extent); for (std::size_t pz = 0; pz <= n; ++pz) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2708,9 +2537,8 @@ void tabulate_polyset_hex_derivs( const T a = 1.0 - 1.0 / static_cast(py); for (std::size_t kz = 0; kz <= nderiv; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, 0, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, 0, kz), md::full_extent, md::full_extent); for (std::size_t pz = 0; pz <= n; ++pz) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2730,12 +2558,10 @@ void tabulate_polyset_hex_derivs( // for py = 1 for (std::size_t kz = 0; kz <= nderiv - ky; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky - 1, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, ky, kz), md::full_extent, md::full_extent); + auto result0 = md::submdspan(P, idx(0, ky - 1, kz), md::full_extent, + md::full_extent); for (std::size_t pz = 0; pz <= n; ++pz) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2752,12 +2578,10 @@ void tabulate_polyset_hex_derivs( const T a = 1.0 - 1.0 / static_cast(py); for (std::size_t kz = 0; kz <= nderiv - ky; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky - 1, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result = md::submdspan(P, idx(0, ky, kz), md::full_extent, + md::full_extent); + auto result0 = md::submdspan(P, idx(0, ky - 1, kz), md::full_extent, + md::full_extent); for (std::size_t pz = 0; pz <= n; ++pz) { for (std::size_t i = 0; i < result.extent(1); ++i) @@ -2780,9 +2604,8 @@ void tabulate_polyset_hex_derivs( { for (std::size_t kz = 0; kz <= nderiv - ky; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result + = md::submdspan(P, idx(0, ky, kz), md::full_extent, md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t pz = 0; pz <= n; ++pz) @@ -2805,9 +2628,8 @@ void tabulate_polyset_hex_derivs( { for (std::size_t kz = 0; kz <= nderiv - ky; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(0, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result = md::submdspan(P, idx(0, ky, kz), md::full_extent, + md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t pz = 0; pz <= n; ++pz) @@ -2834,13 +2656,10 @@ void tabulate_polyset_hex_derivs( { for (std::size_t kz = 0; kz <= nderiv - kx - ky; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky, kz), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result = md::submdspan(P, idx(kx, ky, kz), md::full_extent, + md::full_extent); + auto result0 = md::submdspan(P, idx(kx - 1, ky, kz), md::full_extent, + md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t pz = 0; pz <= n; ++pz) @@ -2865,13 +2684,10 @@ void tabulate_polyset_hex_derivs( { for (std::size_t kz = 0; kz <= nderiv - kx - ky; ++kz) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky, kz), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result = md::submdspan(P, idx(kx, ky, kz), md::full_extent, + md::full_extent); + auto result0 = md::submdspan(P, idx(kx - 1, ky, kz), md::full_extent, + md::full_extent); for (std::size_t py = 0; py <= n; ++py) { for (std::size_t pz = 0; pz <= n; ++pz) @@ -2898,9 +2714,8 @@ void tabulate_polyset_hex_derivs( { for (std::size_t pz = 0; pz <= n; ++pz) { - auto pxyz = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, hex_idx(px, py, pz), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pxyz = md::submdspan(P, md::full_extent, hex_idx(px, py, pz), + md::full_extent); for (std::size_t i = 0; i < pxyz.extent(0); ++i) for (std::size_t j = 0; j < pxyz.extent(1); ++j) pxyz(i, j) *= std::sqrt((2 * px + 1) * (2 * py + 1) * (2 * pz + 1)); @@ -2911,25 +2726,17 @@ void tabulate_polyset_hex_derivs( //----------------------------------------------------------------------------- template void tabulate_polyset_prism_derivs( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - std::size_t n, std::size_t nderiv, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) + md::mdspan> P, std::size_t n, + std::size_t nderiv, md::mdspan> x) { assert(x.extent(1) == 3); assert(P.extent(0) == (nderiv + 1) * (nderiv + 2) * (nderiv + 3) / 6); assert(P.extent(1) == (n + 1) * (n + 1) * (n + 2) / 2); assert(P.extent(2) == x.extent(0)); - auto x0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); - auto x1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 1); - auto x2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - x, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 2); + auto x0 = md::submdspan(x, md::full_extent, 0); + auto x1 = md::submdspan(x, md::full_extent, 1); + auto x2 = md::submdspan(x, md::full_extent, 2); assert(x0.extent(0) > 0); assert(x1.extent(0) > 0); @@ -2958,12 +2765,10 @@ void tabulate_polyset_prism_derivs( { for (std::size_t p = 1; p <= n; ++p) { - auto p0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p - 1, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p, 0, 0), + md::full_extent); + auto p1 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p - 1, 0, 0), + md::full_extent); const T a = static_cast(2 * p - 1) / static_cast(p); for (std::size_t i = 0; i < p0.size(); ++i) { @@ -2973,18 +2778,16 @@ void tabulate_polyset_prism_derivs( if (kx > 0) { - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx - 1, ky, 0), prism_idx(p - 1, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result0 = md::submdspan(P, idx(kx - 1, ky, 0), + prism_idx(p - 1, 0, 0), md::full_extent); for (std::size_t i = 0; i < p0.size(); ++i) p0[i] += 2 * kx * a * result0[i]; } if (ky > 0) { - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, 0), prism_idx(p - 1, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result0 = md::submdspan(P, idx(kx, ky - 1, 0), + prism_idx(p - 1, 0, 0), md::full_extent); for (std::size_t i = 0; i < p0.size(); ++i) p0[i] += ky * a * result0[i]; } @@ -2992,9 +2795,8 @@ void tabulate_polyset_prism_derivs( if (p > 1) { // y^2 terms - auto p2 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p - 2, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p2 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p - 2, 0, 0), + md::full_extent); for (std::size_t i = 0; i < p0.size(); ++i) { T f2 = (1.0 - x1[i]); @@ -3003,9 +2805,8 @@ void tabulate_polyset_prism_derivs( if (ky > 0) { - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, 0), prism_idx(p - 2, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result0 = md::submdspan( + P, idx(kx, ky - 1, 0), prism_idx(p - 2, 0, 0), md::full_extent); for (std::size_t i = 0; i < p0.size(); ++i) { p0[i] @@ -3015,9 +2816,8 @@ void tabulate_polyset_prism_derivs( if (ky > 1) { - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 2, 0), prism_idx(p - 2, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result0 = md::submdspan( + P, idx(kx, ky - 2, 0), prism_idx(p - 2, 0, 0), md::full_extent); for (std::size_t i = 0; i < p0.size(); ++i) p0[i] -= ky * (ky - 1) * result0[i] * (a - 1.0); } @@ -3026,44 +2826,37 @@ void tabulate_polyset_prism_derivs( for (std::size_t p = 0; p < n; ++p) { - auto p0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto p1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p, 1, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto p0 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p, 0, 0), + md::full_extent); + auto p1 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p, 1, 0), + md::full_extent); for (std::size_t i = 0; i < p1.size(); ++i) p1[i] = p0[i] * ((x1[i] * 2.0 - 1.0) * (1.5 + p) + 0.5 + p); if (ky > 0) { - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, 0), prism_idx(p, 0, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result0 = md::submdspan(P, idx(kx, ky - 1, 0), + prism_idx(p, 0, 0), md::full_extent); for (std::size_t i = 0; i < p1.size(); ++i) p1[i] += 2 * ky * (1.5 + p) * result0[i]; } for (std::size_t q = 1; q < n - p; ++q) { - auto pqp1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p, q + 1, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pq = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p, q, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto pqm1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, 0), prism_idx(p, q - 1, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqp1 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p, q + 1, 0), + md::full_extent); + auto pq = md::submdspan(P, idx(kx, ky, 0), prism_idx(p, q, 0), + md::full_extent); + auto pqm1 = md::submdspan(P, idx(kx, ky, 0), prism_idx(p, q - 1, 0), + md::full_extent); const auto [a1, a2, a3] = jrc(2 * p + 1, q); for (std::size_t i = 0; i < p0.size(); ++i) pqp1[i] = pq[i] * ((x1[i] * 2.0 - 1.0) * a1 + a2) - pqm1[i] * a3; if (ky > 0) { - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky - 1, 0), prism_idx(p, q, 0), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result0 = md::submdspan(P, idx(kx, ky - 1, 0), + prism_idx(p, q, 0), md::full_extent); for (std::size_t i = 0; i < pqp1.size(); ++i) pqp1[i] += 2 * ky * a1 * result0[i]; } @@ -3082,13 +2875,10 @@ void tabulate_polyset_prism_derivs( { for (std::size_t ky = 0; ky <= nderiv - kx - kz; ++ky) { - auto result = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); - auto result0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, idx(kx, ky, kz - 1), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto result = md::submdspan(P, idx(kx, ky, kz), md::full_extent, + md::full_extent); + auto result0 = md::submdspan(P, idx(kx, ky, kz - 1), md::full_extent, + md::full_extent); for (std::size_t p = 0; p <= n; ++p) { for (std::size_t q = 0; q <= n - p; ++q) @@ -3132,9 +2922,8 @@ void tabulate_polyset_prism_derivs( { for (std::size_t r = 0; r <= n; ++r) { - auto pqr = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan( - P, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, prism_idx(p, q, r), - MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto pqr = md::submdspan(P, md::full_extent, prism_idx(p, q, r), + md::full_extent); 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((p + 0.5) * (p + q + 1) * (2 * r + 1)) * 2; @@ -3145,14 +2934,9 @@ void tabulate_polyset_prism_derivs( } // namespace //----------------------------------------------------------------------------- template -void polyset::tabulate( - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - P, - cell::type celltype, polyset::type ptype, int d, int n, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) +void polyset::tabulate(md::mdspan> P, + cell::type celltype, polyset::type ptype, int d, int n, + md::mdspan> x) { switch (ptype) { @@ -3216,35 +3000,26 @@ void polyset::tabulate( } //----------------------------------------------------------------------------- template -std::pair, std::array> polyset::tabulate( - cell::type celltype, polyset::type ptype, int d, int n, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - x) +std::pair, std::array> +polyset::tabulate(cell::type celltype, polyset::type ptype, int d, int n, + md::mdspan> x) { std::array shape = {(std::size_t)polyset::nderivs(celltype, n), (std::size_t)polyset::dim(celltype, ptype, d), x.extent(0)}; std::vector P(shape[0] * shape[1] * shape[2]); - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - _P(P.data(), shape); + md::mdspan> _P(P.data(), shape); polyset::tabulate(_P, celltype, ptype, d, n, x); return {std::move(P), std::move(shape)}; } //----------------------------------------------------------------------------- /// @cond template std::pair, std::array> -polyset::tabulate( - cell::type, polyset::type, int, int, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const float, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>); +polyset::tabulate(cell::type, polyset::type, int, int, + md::mdspan>); template std::pair, std::array> -polyset::tabulate( - cell::type, polyset::type, int, int, - MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const double, - MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>); +polyset::tabulate(cell::type, polyset::type, int, int, + md::mdspan>); /// @endcond //----------------------------------------------------------------------------- int polyset::dim(cell::type celltype, polyset::type ptype, int d) diff --git a/cpp/basix/quadrature.cpp b/cpp/basix/quadrature.cpp index 4a433d1a0..b22a1b87d 100644 --- a/cpp/basix/quadrature.cpp +++ b/cpp/basix/quadrature.cpp @@ -5,25 +5,22 @@ #include "quadrature.h" #include "math.h" #include "mdspan.hpp" -#include #include #include #include #include +#include #include #include using namespace basix; +namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; -namespace stdex - = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE; +namespace stdex = md::MDSPAN_IMPL_PROPOSED_NAMESPACE; template -using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; +using mdspan_t = md::mdspan>; template -using mdarray_t - = stdex::mdarray>; +using mdarray_t = stdex::mdarray>; namespace { @@ -54,21 +51,21 @@ std::array, 2> rec_jacobi(int N, T a, T b) std::iota(n.begin(), n.end(), 1.0); std::vector nab(n.size()); std::ranges::transform(n, nab.begin(), - [a, b](auto x) { return 2.0 * x + a + b; }); + [a, b](auto x) { return 2.0 * x + a + b; }); std::vector alpha(N); alpha.front() = nu; - std::ranges::transform(nab, std::next(alpha.begin()), - [a, b](auto nab) - { return (b * b - a * a) / (nab * (nab + 2.0)); }); + std::ranges::transform(nab, std::next(alpha.begin()), [a, b](auto nab) + { return (b * b - a * a) / (nab * (nab + 2.0)); }); std::vector beta(N); beta.front() = mu; - std::ranges::transform(n, nab, std::next(beta.begin()), [a, b](auto n, auto nab) - { - return 4 * (n + a) * (n + b) * n * (n + a + b) - / (nab * nab * (nab + 1.0) * (nab - 1.0)); - }); + std::ranges::transform(n, nab, std::next(beta.begin()), + [a, b](auto n, auto nab) + { + return 4 * (n + a) * (n + b) * n * (n + a + b) + / (nab * nab * (nab + 1.0) * (nab - 1.0)); + }); return {std::move(alpha), std::move(beta)}; } @@ -280,7 +277,8 @@ std::array, 2> make_quadrature_line(int m) { auto [ptx, wx] = compute_gauss_jacobi_rule(0.0, m); std::ranges::transform(wx, wx.begin(), [](auto w) { return 0.5 * w; }); - std::ranges::transform(ptx, ptx.begin(), [](auto x) { return 0.5 * (x + 1.0); }); + std::ranges::transform(ptx, ptx.begin(), + [](auto x) { return 0.5 * (x + 1.0); }); return {std::move(ptx), std::move(wx)}; } //----------------------------------------------------------------------------- @@ -414,7 +412,8 @@ std::array, 2> make_gauss_jacobi_quadrature(cell::type celltype, } case cell::type::pyramid: { - auto [pts, wts] = make_gauss_jacobi_quadrature(cell::type::hexahedron, m + 2); + auto [pts, wts] + = make_gauss_jacobi_quadrature(cell::type::hexahedron, m + 2); mdspan_t x(pts.data(), pts.size() / 3, 3); for (std::size_t i = 0; i < x.extent(0); ++i) { @@ -468,7 +467,8 @@ std::array, 2> make_gll_line(int m) { auto [ptx, wx] = compute_gll_rule(m); std::ranges::transform(wx, wx.begin(), [](auto w) { return 0.5 * w; }); - std::ranges::transform(ptx, ptx.begin(), [](auto x) { return 0.5 * (x + 1.0); }); + std::ranges::transform(ptx, ptx.begin(), + [](auto x) { return 0.5 * (x + 1.0); }); return {ptx, wx}; } //----------------------------------------------------------------------------- @@ -4957,7 +4957,8 @@ template std::vector quadrature::get_gl_points(int m) { std::vector pts = compute_gauss_jacobi_points(0, m); - std::ranges::transform(pts, pts.begin(), [](auto x) { return 0.5 + 0.5 * x; }); + std::ranges::transform(pts, pts.begin(), + [](auto x) { return 0.5 + 0.5 * x; }); return pts; } //----------------------------------------------------------------------------- From 9507310425ae6289623a6ff80f49c229eb5d6638 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Fri, 14 Feb 2025 14:01:33 +0000 Subject: [PATCH 4/4] Doc fixes --- cpp/basix/polynomials.cpp | 8 +++++--- cpp/basix/polyset.cpp | 19 +++++++++++++------ 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/cpp/basix/polynomials.cpp b/cpp/basix/polynomials.cpp index 5064c9d00..b878a0f41 100644 --- a/cpp/basix/polynomials.cpp +++ b/cpp/basix/polynomials.cpp @@ -111,9 +111,11 @@ tabulate_bernstein(cell::type celltype, int d, mdspan_t x) //----------------------------------------------------------------------------- template -std::pair, std::array> -polynomials::tabulate(polynomials::type polytype, cell::type celltype, int d, - md::mdspan> x) +std::pair, std::array> polynomials::tabulate( + polynomials::type polytype, cell::type celltype, int d, + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + x) { switch (polytype) { diff --git a/cpp/basix/polyset.cpp b/cpp/basix/polyset.cpp index 6df1f57ec..0e9ec30b7 100644 --- a/cpp/basix/polyset.cpp +++ b/cpp/basix/polyset.cpp @@ -2934,9 +2934,14 @@ void tabulate_polyset_prism_derivs( } // namespace //----------------------------------------------------------------------------- template -void polyset::tabulate(md::mdspan> P, - cell::type celltype, polyset::type ptype, int d, int n, - md::mdspan> x) +void polyset::tabulate( + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + P, + cell::type celltype, polyset::type ptype, int d, int n, + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + x) { switch (ptype) { @@ -3000,9 +3005,11 @@ void polyset::tabulate(md::mdspan> P, } //----------------------------------------------------------------------------- template -std::pair, std::array> -polyset::tabulate(cell::type celltype, polyset::type ptype, int d, int n, - md::mdspan> x) +std::pair, std::array> polyset::tabulate( + cell::type celltype, polyset::type ptype, int d, int n, + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + x) { std::array shape = {(std::size_t)polyset::nderivs(celltype, n),