Add a hash functions to finite elements (#782)

* add hash function

* fix test

* remove iostream

* Include polyset_type in == and hash for custom elements

* use std::hash and boost's hash_combine

* remove boost

* Change a in hash combine function

cpp/basix/element-families.h

@@ -10,34 +10,34 @@ namespace basix::element
 /// Variants of a Lagrange space that can be created
 enum class lagrange_variant
 {
-  unset = -1,
-  equispaced = 0,
-  gll_warped = 1,
-  gll_isaac = 2,
-  gll_centroid = 3,
-  chebyshev_warped = 4,
-  chebyshev_isaac = 5,
-  chebyshev_centroid = 6,
-  gl_warped = 7,
-  gl_isaac = 8,
-  gl_centroid = 9,
-  legendre = 10,
-  bernstein = 11,
+  unset = 0,
+  equispaced = 1,
+  gll_warped = 2,
+  gll_isaac = 3,
+  gll_centroid = 4,
+  chebyshev_warped = 5,
+  chebyshev_isaac = 6,
+  chebyshev_centroid = 7,
+  gl_warped = 8,
+  gl_isaac = 9,
+  gl_centroid = 10,
+  legendre = 11,
+  bernstein = 12,
 };
 
 /// Variants of a DPC (discontinuous polynomial cubical) space that can
 /// be created. DPC spaces span the same set of polynomials as Lagrange
 /// spaces on simplices but are defined on tensor product cells.
 enum class dpc_variant
 {
-  unset = -1,
-  simplex_equispaced = 0,
-  simplex_gll = 1,
-  horizontal_equispaced = 2,
-  horizontal_gll = 3,
-  diagonal_equispaced = 4,
-  diagonal_gll = 5,
-  legendre = 6,
+  unset = 0,
+  simplex_equispaced = 1,
+  simplex_gll = 2,
+  horizontal_equispaced = 3,
+  horizontal_gll = 4,
+  diagonal_equispaced = 5,
+  diagonal_gll = 6,
+  legendre = 7,
 };
 
 /// Available element families

cpp/basix/finite-element.cpp

@@ -22,7 +22,6 @@
 #include <concepts>
 #include <limits>
 #include <numeric>
-
 #define str_macro(X) #X
 #define str(X) str_macro(X)
 
@@ -182,6 +181,11 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_dual_matrix(
   return {std::move(C), shape};
 }
 //-----------------------------------------------------------------------------
+void combine_hashes(std::size_t& a, std::size_t b)
+{
+  a ^= b + 0x9e3779b9 + (a << 6) + (a >> 2);
+}
+//-----------------------------------------------------------------------------
 } // namespace
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
@@ -1175,7 +1179,8 @@ bool FiniteElement<F>::operator==(const FiniteElement& e) const
            and embedded_superdegree() == e.embedded_superdegree()
            and embedded_subdegree() == e.embedded_subdegree() and coeff_equal
            and entity_dofs() == e.entity_dofs()
-           and dof_ordering() == e.dof_ordering();
+           and dof_ordering() == e.dof_ordering()
+           and polyset_type() == e.polyset_type();
   }
   else
   {
@@ -1189,6 +1194,55 @@ bool FiniteElement<F>::operator==(const FiniteElement& e) const
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point F>
+std::size_t FiniteElement<F>::hash() const
+{
+  std::size_t dof_ordering_hash = 0;
+  for (std::size_t i = 0; i < dof_ordering().size(); ++i)
+  {
+    if (dof_ordering()[i] != static_cast<int>(i))
+    {
+      combine_hashes(dof_ordering_hash,
+                     std::hash<int>{}(dof_ordering()[i] - i));
+    }
+  }
+
+  std::size_t h = std::hash<int>{}(static_cast<int>(family()));
+  combine_hashes(h, dof_ordering_hash);
+  combine_hashes(h, dof_ordering_hash);
+  combine_hashes(h, std::hash<int>{}(static_cast<int>(cell_type())));
+  combine_hashes(h, std::hash<int>{}(static_cast<int>(lagrange_variant())));
+  combine_hashes(h, std::hash<int>{}(static_cast<int>(dpc_variant())));
+  combine_hashes(h, std::hash<int>{}(static_cast<int>(sobolev_space())));
+  combine_hashes(h, std::hash<int>{}(static_cast<int>(map_type())));
+
+  if (family() == element::family::custom)
+  {
+    std::size_t coeff_hash = 0;
+    for (auto i : _coeffs.first)
+    {
+      // This takes five decimal places of each matrix entry. We should revisit
+      // this
+      combine_hashes(coeff_hash, int(i * 100000));
+    }
+    std::size_t vs_hash = 0;
+    for (std::size_t i = 0; i < value_shape().size(); ++i)
+    {
+      combine_hashes(vs_hash, std::hash<int>{}(value_shape()[i]));
+    }
+    combine_hashes(h, coeff_hash);
+    combine_hashes(h, std::hash<int>{}(embedded_superdegree()));
+    combine_hashes(h, std::hash<int>{}(embedded_subdegree()));
+    combine_hashes(h, std::hash<int>{}(static_cast<int>(polyset_type())));
+    combine_hashes(h, vs_hash);
+  }
+  else
+  {
+    combine_hashes(h, std::hash<int>{}(degree()));
+  }
+  return h;
+}
+//-----------------------------------------------------------------------------
+template <std::floating_point F>
 std::pair<std::vector<F>, std::array<std::size_t, 4>>
 FiniteElement<F>::tabulate(int nd, impl::mdspan_t<const F, 2> x) const
 {

cpp/basix/finite-element.h

@@ -348,6 +348,9 @@ class FiniteElement
   /// @return True if elements are the same
   bool operator==(const FiniteElement& e) const;
 
+  /// Get a unique hash of this element
+  std::size_t hash() const;
+
   /// Array shape for tabulate basis values and derivatives at set of
   /// points.
   ///

python/basix/_basixcpp.pyi

@@ -116,6 +116,7 @@ class FiniteElement_float32:
     def push_forward(self, *args, **kwargs) -> Any: ...
     def tabulate(self, *args, **kwargs) -> Any: ...
     def __eq__(self, other) -> Any: ...
+    def __hash__(self) -> int: ...
     @property
     def M(self) -> Any: ...
     @property
@@ -193,6 +194,7 @@ class FiniteElement_float64:
     def push_forward(self, *args, **kwargs) -> Any: ...
     def tabulate(self, *args, **kwargs) -> Any: ...
     def __eq__(self, other) -> Any: ...
+    def __hash__(self) -> int: ...
     @property
     def M(self) -> Any: ...
     @property

python/basix/finite_element.py

@@ -143,6 +143,10 @@ def __eq__(self, other) -> bool:
         except TypeError:
             return False
 
+    def __hash__(self) -> int:
+        """Hash."""
+        return hash(self._e)
+
     def push_forward(self, U, J, detJ, K) -> npt.NDArray[np.floating]:
         """Map function values from the reference to a physical cell.
 

python/wrapper.cpp

@@ -97,6 +97,7 @@ void declare_float(nb::module_& m, std::string type)
              return as_nbarrayp(self.tabulate(n, _x));
            })
       .def("__eq__", &FiniteElement<T>::operator==)
+      .def("__hash__", &FiniteElement<T>::hash)
       .def("push_forward",
            [](const FiniteElement<T>& self,
               nb::ndarray<const T, nb::ndim<3>, nb::c_contig> U,

test/test_elements.py

@@ -66,12 +66,22 @@ def test_dimension(family, cell, degree, dim):
     "family, cell, degree, functions",
     [
         ("Lagrange", "interval", 1, [lambda x: 1 - x[0], lambda x: x[0]]),
-        ("Lagrange", "triangle", 1, [lambda x: 1 - x[0] - x[1], lambda x: x[0], lambda x: x[1]]),
+        (
+            "Lagrange",
+            "triangle",
+            1,
+            [lambda x: 1 - x[0] - x[1], lambda x: x[0], lambda x: x[1]],
+        ),
         (
             "Lagrange",
             "tetrahedron",
             1,
-            [lambda x: 1 - x[0] - x[1] - x[2], lambda x: x[0], lambda x: x[1], lambda x: x[2]],
+            [
+                lambda x: 1 - x[0] - x[1] - x[2],
+                lambda x: x[0],
+                lambda x: x[1],
+                lambda x: x[2],
+            ],
         ),
         (
             "Lagrange",
@@ -116,7 +126,11 @@ def test_dimension(family, cell, degree, dim):
             "Raviart-Thomas",
             "triangle",
             1,
-            [lambda x: (-x[0], -x[1]), lambda x: (x[0] - 1, x[1]), lambda x: (-x[0], 1 - x[1])],
+            [
+                lambda x: (-x[0], -x[1]),
+                lambda x: (x[0] - 1, x[1]),
+                lambda x: (-x[0], 1 - x[1]),
+            ],
         ),
         (
             "Raviart-Thomas",
@@ -133,7 +147,11 @@ def test_dimension(family, cell, degree, dim):
             "N1curl",
             "triangle",
             1,
-            [lambda x: (-x[1], x[0]), lambda x: (x[1], 1 - x[0]), lambda x: (1.0 - x[1], x[0])],
+            [
+                lambda x: (-x[1], x[0]),
+                lambda x: (x[1], 1 - x[0]),
+                lambda x: (1.0 - x[1], x[0]),
+            ],
         ),
         (
             "N1curl",
@@ -160,3 +178,53 @@ def test_values(family, cell, degree, functions):
     for x, t in zip(points, tables):
         for i, f in enumerate(functions):
             assert np.allclose(t[i :: len(functions)], f(x))
+
+
+def test_hash():
+    e0 = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, 1)
+    e1 = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, 1)
+    e2 = basix.create_element(
+        basix.ElementFamily.P, basix.CellType.interval, 1, dof_ordering=[0, 1]
+    )
+    e3 = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, 2)
+    e4 = basix.create_element(basix.ElementFamily.P, basix.CellType.triangle, 1)
+
+    wcoeffs = np.eye(3)
+    z = np.zeros((0, 2))
+    x = [
+        [np.array([[0.0, 0.0]]), np.array([[1.0, 0.0]]), np.array([[0.0, 1.0]])],
+        [z, z, z],
+        [z],
+        [],
+    ]
+    z = np.zeros((0, 1, 0, 1))
+    M = [
+        [np.array([[[[1.0]]]]), np.array([[[[1.0]]]]), np.array([[[[1.0]]]])],
+        [z, z, z],
+        [z],
+        [],
+    ]
+    e5 = basix.create_custom_element(
+        basix.CellType.triangle,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.L2Piola,
+        basix.SobolevSpace.L2,
+        False,
+        1,
+        1,
+        basix.PolysetType.standard,
+    )
+
+    e6 = basix.create_element(basix.ElementFamily.P, basix.CellType.quadrilateral, 2)
+    e7 = basix.create_tp_element(basix.ElementFamily.P, basix.CellType.quadrilateral, 2)
+
+    assert hash(e0) == hash(e1) == hash(e2)
+
+    different_elements = [e2, e3, e4, e5, e6, e7]
+    for i, d0 in enumerate(different_elements):
+        for d1 in different_elements[:i]:
+            assert hash(d0) != hash(d1)