From 3e816e0a29a0eff55e2f2e7663a45e1a515890ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 5 Dec 2024 14:30:38 +0100 Subject: [PATCH 1/8] Add edge jacobians for 3D-1D code compilation --- basix/_basixcpp.pyi | 2 ++ cpp/basix/cell.cpp | 51 ++++++++++++++++++++++++++++++++++++-------- cpp/basix/cell.h | 20 ++++++++++++++++- python/basix/cell.py | 14 ++++++++++++ python/wrapper.cpp | 3 +++ 5 files changed, 80 insertions(+), 10 deletions(-) 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) From cd28bae01d1613cd117f17326a14463d841c311a Mon Sep 17 00:00:00 2001 From: jorgensd Date: Thu, 5 Dec 2024 20:23:18 +0000 Subject: [PATCH 2/8] Update stubs --- python/basix/_basixcpp.pyi | 84 +++++++++++++++++++------------------- python/wrapper.cpp | 14 +++---- 2 files changed, 49 insertions(+), 49 deletions(-) diff --git a/python/basix/_basixcpp.pyi b/python/basix/_basixcpp.pyi index e78e5d913..410532410 100644 --- a/python/basix/_basixcpp.pyi +++ b/python/basix/_basixcpp.pyi @@ -75,19 +75,19 @@ class ElementFamily(enum.IntEnum): iso = 13 class FiniteElement_float32: - def tabulate(self, arg0: int, arg1: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float32', )]: ... + def tabulate(self, arg0: int, arg1: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float32')]: ... def __eq__(self, arg: object, /) -> bool: ... def hash(self) -> int: ... - def permute_subentity_closure(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32', )]: ... + def permute_subentity_closure(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32')]: ... - def permute_subentity_closure_inv(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32', )]: ... + def permute_subentity_closure_inv(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32')]: ... - def push_forward(self, arg0: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None), order='C')], arg1: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None), order='C')], arg2: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None), order='C')], arg3: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float32', )]: ... + def push_forward(self, arg0: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None), order='C', writable=False)], arg1: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None), order='C', writable=False)], arg2: Annotated[ArrayLike, dict(dtype='float32', shape=(None), order='C', writable=False)], arg3: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float32')]: ... - def pull_back(self, arg0: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None), order='C')], arg1: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None), order='C')], arg2: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None), order='C')], arg3: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float32', )]: ... + def pull_back(self, arg0: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None), order='C', writable=False)], arg1: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None), order='C', writable=False)], arg2: Annotated[ArrayLike, dict(dtype='float32', shape=(None), order='C', writable=False)], arg3: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float32')]: ... def T_apply(self, arg0: Annotated[ArrayLike, dict(dtype='float32', shape=(None), order='C')], arg1: int, arg2: int, /) -> None: ... @@ -95,7 +95,7 @@ class FiniteElement_float32: def Tt_inv_apply(self, arg0: Annotated[ArrayLike, dict(dtype='float32', shape=(None), order='C')], arg1: int, arg2: int, /) -> None: ... - def base_transformations(self) -> Annotated[ArrayLike, dict(dtype='float32', )]: ... + def base_transformations(self) -> Annotated[ArrayLike, dict(dtype='float32')]: ... def entity_transformations(self) -> dict: ... @@ -165,26 +165,26 @@ class FiniteElement_float32: def sobolev_space(self) -> SobolevSpace: ... @property - def points(self) -> Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), )]: ... + def points(self) -> Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), writable=False)]: ... @property - def interpolation_matrix(self) -> Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), )]: ... + def interpolation_matrix(self) -> Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), writable=False)]: ... @property - def dual_matrix(self) -> Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), )]: ... + def dual_matrix(self) -> Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), writable=False)]: ... @property - def coefficient_matrix(self) -> Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), )]: + def coefficient_matrix(self) -> Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), writable=False)]: """Coefficient matrix.""" @property - def wcoeffs(self) -> Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), )]: ... + def wcoeffs(self) -> Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), writable=False)]: ... @property - def M(self) -> list[list[Annotated[ArrayLike, dict(dtype='float32', writable=False, )]]]: ... + def M(self) -> list[list[Annotated[ArrayLike, dict(dtype='float32', writable=False)]]]: ... @property - def x(self) -> list[list[Annotated[ArrayLike, dict(dtype='float32', writable=False, )]]]: ... + def x(self) -> list[list[Annotated[ArrayLike, dict(dtype='float32', writable=False)]]]: ... @property def has_tensor_product_factorisation(self) -> bool: ... @@ -199,19 +199,19 @@ class FiniteElement_float32: def dtype(self) -> str: ... class FiniteElement_float64: - def tabulate(self, arg0: int, arg1: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... + def tabulate(self, arg0: int, arg1: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def __eq__(self, arg: object, /) -> bool: ... def hash(self) -> int: ... - def permute_subentity_closure(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32', )]: ... + def permute_subentity_closure(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32')]: ... - def permute_subentity_closure_inv(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32', )]: ... + def permute_subentity_closure_inv(self, arg0: Annotated[ArrayLike, dict(dtype='int32', shape=(None), order='C')], arg1: int, arg2: CellType, /) -> Annotated[ArrayLike, dict(dtype='int32')]: ... - def push_forward(self, arg0: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None), order='C')], arg1: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None), order='C')], arg2: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None), order='C')], arg3: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... + def push_forward(self, arg0: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None), order='C', writable=False)], arg1: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None), order='C', writable=False)], arg2: Annotated[ArrayLike, dict(dtype='float64', shape=(None), order='C', writable=False)], arg3: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... - def pull_back(self, arg0: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None), order='C')], arg1: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None), order='C')], arg2: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None), order='C')], arg3: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... + def pull_back(self, arg0: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None), order='C', writable=False)], arg1: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None), order='C', writable=False)], arg2: Annotated[ArrayLike, dict(dtype='float64', shape=(None), order='C', writable=False)], arg3: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def T_apply(self, arg0: Annotated[ArrayLike, dict(dtype='float64', shape=(None), order='C')], arg1: int, arg2: int, /) -> None: ... @@ -219,7 +219,7 @@ class FiniteElement_float64: def Tt_inv_apply(self, arg0: Annotated[ArrayLike, dict(dtype='float64', shape=(None), order='C')], arg1: int, arg2: int, /) -> None: ... - def base_transformations(self) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... + def base_transformations(self) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def entity_transformations(self) -> dict: ... @@ -289,26 +289,26 @@ class FiniteElement_float64: def sobolev_space(self) -> SobolevSpace: ... @property - def points(self) -> Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), )]: ... + def points(self) -> Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), writable=False)]: ... @property - def interpolation_matrix(self) -> Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), )]: ... + def interpolation_matrix(self) -> Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), writable=False)]: ... @property - def dual_matrix(self) -> Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), )]: ... + def dual_matrix(self) -> Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), writable=False)]: ... @property - def coefficient_matrix(self) -> Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), )]: + def coefficient_matrix(self) -> Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), writable=False)]: """Coefficient matrix.""" @property - def wcoeffs(self) -> Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), )]: ... + def wcoeffs(self) -> Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), writable=False)]: ... @property - def M(self) -> list[list[Annotated[ArrayLike, dict(dtype='float64', writable=False, )]]]: ... + def M(self) -> list[list[Annotated[ArrayLike, dict(dtype='float64', writable=False)]]]: ... @property - def x(self) -> list[list[Annotated[ArrayLike, dict(dtype='float64', writable=False, )]]]: ... + def x(self) -> list[list[Annotated[ArrayLike, dict(dtype='float64', writable=False)]]]: ... @property def has_tensor_product_factorisation(self) -> bool: ... @@ -434,35 +434,37 @@ class SobolevSpace(enum.IntEnum): HDivDiv = 13 -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_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]: ... -def cell_facet_outward_normals(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def cell_facet_outward_normals(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... -def cell_facet_reference_volumes(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def cell_facet_reference_volumes(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def cell_volume(arg: CellType, /) -> float: ... @overload -def compute_interpolation_operator(arg0: FiniteElement_float32, arg1: FiniteElement_float32, /) -> Annotated[ArrayLike, dict(dtype='float32', )]: ... +def compute_interpolation_operator(arg0: FiniteElement_float32, arg1: FiniteElement_float32, /) -> Annotated[ArrayLike, dict(dtype='float32')]: ... @overload -def compute_interpolation_operator(arg0: FiniteElement_float64, arg1: FiniteElement_float64, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def compute_interpolation_operator(arg0: FiniteElement_float64, arg1: FiniteElement_float64, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... -def create_custom_element_float32(cell_type: CellType, value_shape: Sequence[int], wcoeffs: Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), order='C')], x: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None), order='C')]]], M: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float32', writable=False, shape=(None, None, None, None), order='C')]]], interpolation_nderivs: int, map_type: MapType, sobolev_space: SobolevSpace, discontinuous: bool, embedded_subdegree: int, embedded_superdegree: int, poly_type: PolysetType) -> FiniteElement_float32: ... +def create_custom_element_float32(cell_type: CellType, value_shape: Sequence[int], wcoeffs: Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), order='C', writable=False)], x: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float32', shape=(None, None), order='C', writable=False)]]], M: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float32', shape=(None, None, None, None), order='C', writable=False)]]], interpolation_nderivs: int, map_type: MapType, sobolev_space: SobolevSpace, discontinuous: bool, embedded_subdegree: int, embedded_superdegree: int, poly_type: PolysetType) -> FiniteElement_float32: ... -def create_custom_element_float64(cell_type: CellType, value_shape: Sequence[int], wcoeffs: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), order='C')], x: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), order='C')]]], M: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None, None, None), order='C')]]], interpolation_nderivs: int, map_type: MapType, sobolev_space: SobolevSpace, discontinuous: bool, embedded_subdegree: int, embedded_superdegree: int, poly_type: PolysetType) -> FiniteElement_float64: ... +def create_custom_element_float64(cell_type: CellType, value_shape: Sequence[int], wcoeffs: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), order='C', writable=False)], x: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), order='C', writable=False)]]], M: Sequence[Sequence[Annotated[ArrayLike, dict(dtype='float64', shape=(None, None, None, None), order='C', writable=False)]]], interpolation_nderivs: int, map_type: MapType, sobolev_space: SobolevSpace, discontinuous: bool, embedded_subdegree: int, embedded_superdegree: int, poly_type: PolysetType) -> FiniteElement_float64: ... def create_element(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, arg6: Sequence[int], arg7: str, /) -> FiniteElement_float32 | FiniteElement_float64: ... -def create_lattice(arg0: CellType, arg1: int, arg2: LatticeType, arg3: bool, arg4: LatticeSimplexMethod, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def create_lattice(arg0: CellType, arg1: int, arg2: LatticeType, arg3: bool, arg4: LatticeSimplexMethod, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def create_tp_element(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, arg6: str, /) -> FiniteElement_float32 | FiniteElement_float64: ... -def geometry(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def geometry(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... @overload def index(arg: int, /) -> int: ... @@ -473,7 +475,7 @@ def index(arg0: int, arg1: int, /) -> int: ... @overload def index(arg0: int, arg1: int, arg2: int, /) -> int: ... -def make_quadrature(arg0: QuadratureType, arg1: CellType, arg2: PolysetType, arg3: int, /) -> tuple[Annotated[ArrayLike, dict(dtype='float64', )], Annotated[ArrayLike, dict(dtype='float64', )]]: ... +def make_quadrature(arg0: QuadratureType, arg1: CellType, arg2: PolysetType, arg3: int, /) -> tuple[Annotated[ArrayLike, dict(dtype='float64')], Annotated[ArrayLike, dict(dtype='float64')]]: ... def polynomials_dim(arg0: PolynomialType, arg1: CellType, arg2: int, /) -> int: ... @@ -481,17 +483,17 @@ def restriction(arg0: PolysetType, arg1: CellType, arg2: CellType, /) -> Polyset def sobolev_space_intersection(arg0: SobolevSpace, arg1: SobolevSpace, /) -> SobolevSpace: ... -def subentity_types(arg: CellType, /) -> list[list[CellType]]: ... - def sub_entity_connectivity(arg: CellType, /) -> list[list[list[list[int]]]]: ... -def sub_entity_geometry(arg0: CellType, arg1: int, arg2: int, /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def sub_entity_geometry(arg0: CellType, arg1: int, arg2: int, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... + +def subentity_types(arg: CellType, /) -> list[list[CellType]]: ... def superset(arg0: CellType, arg1: PolysetType, arg2: PolysetType, /) -> PolysetType: ... def tabulate_polynomial_set(celltype: CellType, polytype: PolysetType, d: int, n: int, x: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), order='C')]) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... -def tabulate_polynomials(arg0: PolynomialType, arg1: CellType, arg2: int, arg3: Annotated[ArrayLike, dict(dtype='float64', writable=False, shape=(None, None), order='C')], /) -> Annotated[ArrayLike, dict(dtype='float64', )]: ... +def tabulate_polynomials(arg0: PolynomialType, arg1: CellType, arg2: int, arg3: Annotated[ArrayLike, dict(dtype='float64', shape=(None, None), order='C', writable=False)], /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def topology(arg: CellType, /) -> list[list[list[int]]]: ... diff --git a/python/wrapper.cpp b/python/wrapper.cpp index 16b8d89b3..c15a95848 100644 --- a/python/wrapper.cpp +++ b/python/wrapper.cpp @@ -510,14 +510,12 @@ NB_MODULE(_basixcpp, m) { return cell::volume(cell_type); }); m.def("cell_facet_normals", [](cell::type cell_type) { return as_nbarrayp(cell::facet_normals(cell_type)); }); - m.def("cell_facet_reference_volumes", - [](cell::type cell_type) { - return as_nbarray(cell::facet_reference_volumes(cell_type)); - }); - m.def("cell_facet_outward_normals", - [](cell::type cell_type) { - return as_nbarrayp(cell::facet_outward_normals(cell_type)); - }); + m.def( + "cell_facet_reference_volumes", [](cell::type cell_type) + { return as_nbarray(cell::facet_reference_volumes(cell_type)); }); + m.def( + "cell_facet_outward_normals", [](cell::type cell_type) + { return as_nbarrayp(cell::facet_outward_normals(cell_type)); }); m.def("cell_facet_orientations", [](cell::type cell_type) { From aaa5b2536a454eacb843ae6100823974c5ac0bc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 14 Jan 2025 15:40:20 +0000 Subject: [PATCH 3/8] Fix stubs --- python/basix/_basixcpp.pyi | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/python/basix/_basixcpp.pyi b/python/basix/_basixcpp.pyi index e16499b50..4736eb584 100644 --- a/python/basix/_basixcpp.pyi +++ b/python/basix/_basixcpp.pyi @@ -434,7 +434,7 @@ class SobolevSpace(enum.IntEnum): HDivDiv = 13 -def cell_ridge_jacobians(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... +def cell_edge_jacobians(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... def cell_facet_jacobians(arg: CellType, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... @@ -483,14 +483,12 @@ def restriction(arg0: PolysetType, arg1: CellType, arg2: CellType, /) -> Polyset def sobolev_space_intersection(arg0: SobolevSpace, arg1: SobolevSpace, /) -> SobolevSpace: ... -def subentity_types(arg: CellType, /) -> list[list[CellType]]: ... - -def sub_entity_type(arg: CellType, dim: int, index: int, /) -> CellType: ... - def sub_entity_connectivity(arg: CellType, /) -> list[list[list[list[int]]]]: ... def sub_entity_geometry(arg0: CellType, arg1: int, arg2: int, /) -> Annotated[ArrayLike, dict(dtype='float64')]: ... +def sub_entity_type(arg0: CellType, arg1: int, arg2: int, /) -> CellType: ... + def subentity_types(arg: CellType, /) -> list[list[CellType]]: ... def superset(arg0: CellType, arg1: PolysetType, arg2: PolysetType, /) -> PolysetType: ... From 3f7ef9adcb8d6543725ac738cd6f347a402f7433 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 16 Jan 2025 15:03:45 +0000 Subject: [PATCH 4/8] Add tests --- test/test_cell.py | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/test/test_cell.py b/test/test_cell.py index 1553da731..582434e29 100644 --- a/test/test_cell.py +++ b/test/test_cell.py @@ -89,3 +89,47 @@ def test_sub_entity_type(): for i in range(4): assert basix.cell.sub_entity_type(cell_type, 2, i) == basix.CellType.triangle assert basix.cell.sub_entity_type(cell_type, 3, 0) == basix.CellType.tetrahedron + + +def test_facet_jacobians_2D_simplex(): + cell = basix.cell.CellType.triangle + facet_jacobian = basix.cell.facet_jacobians(cell) + + reference_vertices = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]) + mask = np.zeros(3, dtype=np.bool_) + for i in range(3): + mask[:] = True + mask[i] = False + facet = reference_vertices[mask] + + reference_facet_jacobian = -facet[0:1, :].T + facet[1:2, :].T + np.testing.assert_allclose(reference_facet_jacobian, facet_jacobian[i]) + + +def test_facet_jacobians_3D_simplex(): + cell = basix.cell.CellType.tetrahedron + facet_jacobian = basix.cell.facet_jacobians(cell) + + reference_vertices = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + ) + mask = np.zeros(4, dtype=np.bool_) + for i in range(4): + mask[:] = True + mask[i] = False + facet = reference_vertices[mask] + reference_facet_jacobian = np.array([-facet[0] + facet[1], -facet[0] + facet[2]]).T + np.testing.assert_allclose(reference_facet_jacobian, facet_jacobian[i]) + + +@pytest.mark.parametrize("cell", [basix.cell.CellType.hexahedron, basix.cell.CellType.tetrahedron]) +def test_edge_jacobian_3D_simplex(cell): + edge_jacobian = basix.cell.edge_jacobians(cell) + geom = basix.geometry(cell) + topology = basix.topology(cell) + edges = topology[1] + + for i, edge in enumerate(edges): + points = geom[edge] + reference_edge_jacobian = (points[1:2, :] - points[0:1, :]).T + np.testing.assert_allclose(reference_edge_jacobian, edge_jacobian[i]) From 833299df5d2b9e28cba6685a6c546a34aa35ec23 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 16 Jan 2025 15:04:48 +0000 Subject: [PATCH 5/8] Update copyright notice --- test/test_cell.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_cell.py b/test/test_cell.py index 582434e29..1dbcc4836 100644 --- a/test/test_cell.py +++ b/test/test_cell.py @@ -1,4 +1,4 @@ -# Copyright (c) 2021 Matthew Scroggs +# Copyright (c) 2021-2025 Matthew Scroggs, Jørgen S. Dokken # FEniCS Project # SPDX-License-Identifier: MIT From dfcfe1a2f2233cd8aced4fb7d62e449d2587a2c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 16 Jan 2025 15:33:50 +0000 Subject: [PATCH 6/8] Remove cell dimension from jacobian --- cpp/basix/cell.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/basix/cell.cpp b/cpp/basix/cell.cpp index 20e0da1d4..9894c13b4 100644 --- a/cpp/basix/cell.cpp +++ b/cpp/basix/cell.cpp @@ -611,7 +611,7 @@ std::pair, std::array> cell::entity_jacobians(cell::type cell_type, std::size_t e_dim) { std::size_t tdim = cell::topological_dimension(cell_type); - if ((e_dim > tdim)) + if ((e_dim >= tdim)) { throw std::runtime_error( "Entity jacobians supported for entity of dimension " From 1838c73a2db8b2e09932df4cf2b1709bb476f139 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 16 Jan 2025 15:34:03 +0000 Subject: [PATCH 7/8] Add prism and pyramid --- test/test_cell.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_cell.py b/test/test_cell.py index 1dbcc4836..4d862a16d 100644 --- a/test/test_cell.py +++ b/test/test_cell.py @@ -122,7 +122,8 @@ def test_facet_jacobians_3D_simplex(): np.testing.assert_allclose(reference_facet_jacobian, facet_jacobian[i]) -@pytest.mark.parametrize("cell", [basix.cell.CellType.hexahedron, basix.cell.CellType.tetrahedron]) +@pytest.mark.parametrize("cell", [basix.cell.CellType.hexahedron, basix.cell.CellType.tetrahedron, + basix.cell.CellType.prism, basix.cell.CellType.pyramid]) def test_edge_jacobian_3D_simplex(cell): edge_jacobian = basix.cell.edge_jacobians(cell) geom = basix.geometry(cell) From ee94913086eb28b75697e36e4616b9df5cdd985f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 16 Jan 2025 15:47:13 +0000 Subject: [PATCH 8/8] Ruff formatting --- test/test_cell.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/test/test_cell.py b/test/test_cell.py index 4d862a16d..3748d4bec 100644 --- a/test/test_cell.py +++ b/test/test_cell.py @@ -122,8 +122,15 @@ def test_facet_jacobians_3D_simplex(): np.testing.assert_allclose(reference_facet_jacobian, facet_jacobian[i]) -@pytest.mark.parametrize("cell", [basix.cell.CellType.hexahedron, basix.cell.CellType.tetrahedron, - basix.cell.CellType.prism, basix.cell.CellType.pyramid]) +@pytest.mark.parametrize( + "cell", + [ + basix.cell.CellType.hexahedron, + basix.cell.CellType.tetrahedron, + basix.cell.CellType.prism, + basix.cell.CellType.pyramid, + ], +) def test_edge_jacobian_3D_simplex(cell): edge_jacobian = basix.cell.edge_jacobians(cell) geom = basix.geometry(cell)