Merge branch 'main' into dependabot/github_actions/pypa/cibuildwheel-2.22.0

Author: jorgensd

.github/workflows/pythonapp.yml

@@ -17,7 +17,9 @@ on:
     branches:
       - main
   workflow_dispatch:
-
+  # Weekly build on Mondays at 8 am
+  schedule:
+    - cron: "0 8 * * 1"
 jobs:
   lint:
     name: Lint code
@@ -26,7 +28,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: '3.12'
+          python-version: "3.12"
       - uses: actions/checkout@v4
       - name: Ruff check
         run: |
@@ -144,7 +146,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: '3.12'
+          python-version: "3.12"
       - name: Install dependencies
         run: sudo apt-get install -y libopenblas-dev liblapack-dev ninja-build
       - name: Install Python dependencies
@@ -160,7 +162,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: '3.12'
+          python-version: "3.12"
       - name: Install dependencies
         run: sudo apt-get install -y libopenblas-dev liblapack-dev
       - name: Install Basix C++ library

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]: ...
@@ -459,6 +461,8 @@ def restriction(arg0: PolysetType, arg1: CellType, arg2: CellType, /) -> Polyset
 
 def sobolev_space_intersection(arg0: SobolevSpace, arg1: SobolevSpace, /) -> SobolevSpace: ...
 
+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', )]: ...

cpp/basix/cell.cpp

@@ -608,33 +608,62 @@ std::vector<std::vector<cell::type>> cell::subentity_types(cell::type cell_type)
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 3>>
-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<T>(cell_type);
   mdspan_t<const T, 2> x(_x.data(), xshape);
-  std::vector<std::vector<int>> facets = topology(cell_type)[tdim - 1];
+  std::vector<std::vector<int>> entities = topology(cell_type)[e_dim];
 
-  std::array<std::size_t, 3> shape = {facets.size(), tdim, tdim - 1};
+  std::array<std::size_t, 3> shape = {entities.size(), tdim, e_dim};
   std::vector<T> jacobians(shape[0] * shape[1] * shape[2]);
   mdspan_t<T, 3> 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<int>& facet = facets[f];
-    for (std::size_t j = 0; j < tdim - 1; ++j)
+    const std::vector<int>& 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::floating_point T>
+std::pair<std::vector<T>, std::array<std::size_t, 3>>
+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<T>(cell_type, tdim - 1);
+}
+//-----------------------------------------------------------------------------
+template <std::floating_point T>
+std::pair<std::vector<T>, std::array<std::size_t, 3>>
+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<T>(cell_type, 1);
+}
+//-----------------------------------------------------------------------------
 
 /// @cond
 // Explicit instantiation for double and float
@@ -668,6 +697,10 @@ template std::pair<std::vector<float>, std::array<std::size_t, 3>>
     cell::facet_jacobians(cell::type);
 template std::pair<std::vector<double>, std::array<std::size_t, 3>>
     cell::facet_jacobians(cell::type);
+template std::pair<std::vector<float>, std::array<std::size_t, 3>>
+    cell::edge_jacobians(cell::type);
+template std::pair<std::vector<double>, std::array<std::size_t, 3>>
+    cell::edge_jacobians(cell::type);
 /// @endcond
 
 //-----------------------------------------------------------------------------

cpp/basix/cell.h

@@ -121,11 +121,29 @@ std::vector<T> facet_reference_volumes(cell::type cell_type);
 /// @return The subentity types. Indices are (tdim, entity)
 std::vector<std::vector<cell::type>> 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::floating_point T>
+std::pair<std::vector<T>, std::array<std::size_t, 3>>
+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::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 3>>
 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::floating_point T>
+std::pair<std::vector<T>, std::array<std::size_t, 3>>
+edge_jacobians(cell::type cell_type);
+
 } // namespace basix::cell

cpp/basix/mdspan.hpp

@@ -6177,7 +6177,7 @@ layout_stride::mapping<Extents>::submdspan_mapping_impl(
     (defined(__NVCC__) &&                                                      \
      (__CUDACC_VER_MAJOR__ * 100 + __CUDACC_VER_MINOR__ * 10) < 1120)
                       MDSPAN_IMPL_STANDARD_NAMESPACE::detail::tuple<decltype(detail::stride_of(slices))...>(
-                          detail::stride_of(slices)...).values)),
+                          detail::stride_of(slices)...)).values),
 #else
                       MDSPAN_IMPL_STANDARD_NAMESPACE::detail::tuple(detail::stride_of(slices)...)).values),
 #endif

demo/python/demo_create_and_tabulate.py

@@ -46,7 +46,7 @@
 # is 1 as we are only tabulating the function values; it would be higher if we
 # had asked for derivatives too. The second dimension (5) is the number of points.
 # The third dimension (25) is the number of DOFs. The fourth dimension (1) is the
-# value size of the element: this will be greater than 1 for vector-values elements.
+# value size of the element: this will be greater than 1 for vector-valued elements.
 #
 # C++ demo
 # ========