diff --git a/basix/_basixcpp.pyi b/basix/_basixcpp.pyi index caa601e26..cb94e11a7 100644 --- a/basix/_basixcpp.pyi +++ b/basix/_basixcpp.pyi @@ -414,6 +414,8 @@ class SobolevSpace(enum.IntEnum): def cell_facet_jacobians(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def cell_edge_jacobians(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... + def cell_facet_normals(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... def cell_facet_orientations(arg: CellType, /) -> list[int]: ... diff --git a/cpp/basix/cell.cpp b/cpp/basix/cell.cpp index f3a53b6cc..20e0da1d4 100644 --- a/cpp/basix/cell.cpp +++ b/cpp/basix/cell.cpp @@ -608,33 +608,62 @@ std::vector> cell::subentity_types(cell::type cell_type) //----------------------------------------------------------------------------- template std::pair, std::array> -cell::facet_jacobians(cell::type cell_type) +cell::entity_jacobians(cell::type cell_type, std::size_t e_dim) { std::size_t tdim = cell::topological_dimension(cell_type); - if (tdim != 2 and tdim != 3) + if ((e_dim > tdim)) { throw std::runtime_error( - "Facet jacobians not supported for this cell type."); + "Entity jacobians supported for entity of dimension " + + std::to_string(e_dim) + "for cell of dimension " + + std::to_string(tdim)); } const auto [_x, xshape] = cell::geometry(cell_type); mdspan_t x(_x.data(), xshape); - std::vector> facets = topology(cell_type)[tdim - 1]; + std::vector> entities = topology(cell_type)[e_dim]; - std::array shape = {facets.size(), tdim, tdim - 1}; + std::array shape = {entities.size(), tdim, e_dim}; std::vector jacobians(shape[0] * shape[1] * shape[2]); mdspan_t J(jacobians.data(), shape); - for (std::size_t f = 0; f < facets.size(); ++f) + for (std::size_t e = 0; e < entities.size(); ++e) { - const std::vector& facet = facets[f]; - for (std::size_t j = 0; j < tdim - 1; ++j) + const std::vector& entity = entities[e]; + for (std::size_t j = 0; j < e_dim; ++j) for (std::size_t k = 0; k < J.extent(1); ++k) - J(f, k, j) = x(facet[1 + j], k) - x(facet[0], k); + J(e, k, j) = x(entity[1 + j], k) - x(entity[0], k); } return {std::move(jacobians), std::move(shape)}; } //----------------------------------------------------------------------------- +template +std::pair, std::array> +cell::facet_jacobians(cell::type cell_type) +{ + std::size_t tdim = cell::topological_dimension(cell_type); + if (tdim != 2 and tdim != 3) + { + throw std::runtime_error( + "Facet jacobians not supported for this cell type."); + } + + return entity_jacobians(cell_type, tdim - 1); +} +//----------------------------------------------------------------------------- +template +std::pair, std::array> +cell::edge_jacobians(cell::type cell_type) +{ + std::size_t tdim = cell::topological_dimension(cell_type); + if (tdim != 3) + { + throw std::runtime_error( + "Edge jacobians not supported for this cell type."); + } + return entity_jacobians(cell_type, 1); +} +//----------------------------------------------------------------------------- /// @cond // Explicit instantiation for double and float @@ -668,6 +697,10 @@ template std::pair, std::array> cell::facet_jacobians(cell::type); template std::pair, std::array> cell::facet_jacobians(cell::type); +template std::pair, std::array> + cell::edge_jacobians(cell::type); +template std::pair, std::array> + cell::edge_jacobians(cell::type); /// @endcond //----------------------------------------------------------------------------- diff --git a/cpp/basix/cell.h b/cpp/basix/cell.h index 6fcce0089..1f8802033 100644 --- a/cpp/basix/cell.h +++ b/cpp/basix/cell.h @@ -121,11 +121,29 @@ std::vector facet_reference_volumes(cell::type cell_type); /// @return The subentity types. Indices are (tdim, entity) std::vector> subentity_types(cell::type cell_type); +/// Get the jacobians of a set of entities of a reference cell +/// @param cell_type Type of cell +/// @param e_dim Dimension of entities +/// @return The jacobians of the facets (flattened) and the shape (nfacets, +/// tdim, tdim - 1). +template +std::pair, std::array> +entity_jacobians(cell::type cell_type, std::size_t e_dim); + /// Get the jacobians of the facets of a reference cell /// @param cell_type Type of cell -/// @return The jacobians of the facets. Shape is (nfacets, gdim, gdim - 1) +/// @return The jacobians of the facets (flattened) and the shape (nfacets, +/// tdim, tdim - 1). template std::pair, std::array> facet_jacobians(cell::type cell_type); +/// Get the jacobians of the edegs of a reference cell +/// @param cell_type Type of cell +/// @return The jacobians of the edges (flattened) and the shape (nedges, tdim, +/// tdim-2). +template +std::pair, std::array> +edge_jacobians(cell::type cell_type); + } // namespace basix::cell diff --git a/python/basix/cell.py b/python/basix/cell.py index 4fd13f300..5ad807b30 100644 --- a/python/basix/cell.py +++ b/python/basix/cell.py @@ -9,6 +9,7 @@ from basix._basixcpp import CellType from basix._basixcpp import cell_facet_jacobians as _fj +from basix._basixcpp import cell_edge_jacobians as _ej from basix._basixcpp import cell_facet_normals as _fn from basix._basixcpp import cell_facet_orientations as _fo from basix._basixcpp import cell_facet_outward_normals as _fon @@ -23,6 +24,7 @@ "sub_entity_connectivity", "subentity_types", "volume", + "edge_jacobians", "facet_jacobians", "facet_normals", "facet_orientations", @@ -67,6 +69,18 @@ def facet_jacobians(celltype: CellType) -> npt.ArrayLike: return _fj(celltype) +def edge_jacobians(celltype: CellType) -> npt.ArrayLike: + """Jacobians of the edges of a reference cell. + + Args: + celltype: cell type. + + Returns: + Jacobians of the edges. + """ + return _ej(celltype) + + def facet_normals(celltype: CellType) -> npt.ArrayLike: """Normals to the facets of a reference cell. diff --git a/python/wrapper.cpp b/python/wrapper.cpp index 075ad7ab6..16b8d89b3 100644 --- a/python/wrapper.cpp +++ b/python/wrapper.cpp @@ -528,6 +528,9 @@ NB_MODULE(_basixcpp, m) m.def("cell_facet_jacobians", [](cell::type cell_type) { return as_nbarrayp(cell::facet_jacobians(cell_type)); }); + m.def("cell_edge_jacobians", [](cell::type cell_type) + { return as_nbarrayp(cell::edge_jacobians(cell_type)); }); + nb::enum_(m, "ElementFamily", nb::is_arithmetic(), "Finite element family.") .value("custom", element::family::custom)