Remove assumption that geometry uses real scalars (#110)

* remove real from reference_cell

* default f64

* replace more real with TGeo scalar type

* remove another real

* remove real from bindings
Also cargo fmt on other files

* remove more real

* clippy

* fix dtypes

* geo_dtype

Author: mscroggs

examples/element_family.rs

@@ -5,7 +5,7 @@ use ndelement::types::{Continuity, ReferenceCellType};
 fn main() {
     // Create the degree 2 Lagrange element family. A family is a set of finite elements with the
     // same family type, degree, and continuity across a set of cells
-    let family = LagrangeElementFamily::<f64>::new(2, Continuity::Standard);
+    let family = LagrangeElementFamily::<f64, f64>::new(2, Continuity::Standard);
 
     // Get the element in the family on a triangle
     let element = family.element(ReferenceCellType::Triangle);

examples/lagrange_element.rs

@@ -7,7 +7,8 @@ use rlst::{DynArray, rlst_dynamic_array};
 
 fn main() {
     // Create a P2 element on a triangle
-    let element = lagrange::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
+    let element =
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
 
     println!("This element has {} basis functions.", element.dim());
 

examples/test_high_degree.rs

@@ -10,7 +10,7 @@ fn main() {
             paste! {
                 for d in 1..[<$max_degree>] {
                     println!("Constructing Lagrange(degree={d}, cell={:?})", ReferenceCellType::[<$cell>]);
-                    let _e = lagrange::create::<f64>(ReferenceCellType::[<$cell>], d, Continuity::Standard);
+                    let _e = lagrange::create::<f64, f64>(ReferenceCellType::[<$cell>], d, Continuity::Standard);
                 }
             }
         };
@@ -21,7 +21,7 @@ fn main() {
             paste! {
                 for d in 1..[<$max_degree>] {
                     println!("Constructing RaviartThomas(degree={d}, cell={:?})", ReferenceCellType::[<$cell>]);
-                    let _e = raviart_thomas::create::<f64>(ReferenceCellType::[<$cell>], d, Continuity::Standard);
+                    let _e = raviart_thomas::create::<f64, f64>(ReferenceCellType::[<$cell>], d, Continuity::Standard);
                 }
             }
         };
@@ -32,7 +32,7 @@ fn main() {
             paste! {
                 for d in 1..[<$max_degree>] {
                     println!("Constructing Nedelec(degree={d}, cell={:?})", ReferenceCellType::[<$cell>]);
-                    let _e = nedelec::create::<f64>(ReferenceCellType::[<$cell>], d, Continuity::Standard);
+                    let _e = nedelec::create::<f64, f64>(ReferenceCellType::[<$cell>], d, Continuity::Standard);
                 }
             }
         };

python/ndelement/ciarlet.py

@@ -52,6 +52,11 @@ def dtype(self) -> typing.Type[np.floating]:
         """Data type."""
         return _dtypes[_lib.ciarlet_element_dtype(self._rs_element)]
 
+    @property
+    def geo_dtype(self) -> typing.Type[np.floating]:
+        """Data type."""
+        return _dtypes[_lib.ciarlet_element_geo_dtype(self._rs_element)]
+
     @property
     def value_size(self) -> int:
         """Value size of the element."""
@@ -123,7 +128,7 @@ def interpolation_points(self) -> list[list[npt.NDArray]]:
             points_d = []
             for i in range(n):
                 shape = (_lib.ciarlet_element_interpolation_npoints(self._rs_element, d, i), tdim)
-                points_di = np.empty(shape, dtype=self.dtype(0).real.dtype)
+                points_di = np.empty(shape, dtype=self.geo_dtype)
                 _lib.ciarlet_element_interpolation_points(
                     self._rs_element, d, i, _ffi.cast("void*", points_di.ctypes.data)
                 )
@@ -152,20 +157,31 @@ def interpolation_weights(self) -> list[list[npt.NDArray]]:
 
     def tabulate(self, points: npt.NDArray[np.floating], nderivs: int) -> npt.NDArray:
         """Tabulate the basis functions at a set of points."""
-        if points.dtype != self.dtype(0).real.dtype:
-            raise TypeError("points has incorrect type")
         shape = np.empty(4, dtype=np.uintp)
         _lib.ciarlet_element_tabulate_array_shape(
             self._rs_element, nderivs, points.shape[0], _ffi.cast("uintptr_t*", shape.ctypes.data)
         )
         data = np.empty(shape[::-1], dtype=self.dtype)
-        _lib.ciarlet_element_tabulate(
-            self._rs_element,
-            _ffi.cast("void*", points.ctypes.data),
-            points.shape[0],
-            nderivs,
-            _ffi.cast("void*", data.ctypes.data),
-        )
+
+        if points.dtype == np.float64:
+            _lib.ciarlet_element_tabulate_f64(
+                self._rs_element,
+                _ffi.cast("double*", points.ctypes.data),
+                points.shape[0],
+                nderivs,
+                _ffi.cast("void*", data.ctypes.data),
+            )
+        elif points.dtype == np.float32:
+            _lib.ciarlet_element_tabulate_f32(
+                self._rs_element,
+                _ffi.cast("float*", points.ctypes.data),
+                points.shape[0],
+                nderivs,
+                _ffi.cast("void*", data.ctypes.data),
+            )
+        else:
+            raise TypeError(f"Unsupported dtype: {points.dtype}")
+
         return data
 
     def physical_value_size(self, gdim: int) -> int:
@@ -189,13 +205,12 @@ def push_forward(
         inverse_jacobians: npt.NDArray[np.floating],
     ) -> npt.NDArray[np.floating]:
         """Push values forward to a physical cell."""
-        if reference_values.dtype != self.dtype(0):
-            raise TypeError("reference_values has incorrect type")
-        if jacobians.dtype != self.dtype(0).real.dtype:
+        geo_dtype = reference_values.dtype
+        if jacobians.dtype != geo_dtype:
             raise TypeError("jacobians has incorrect type")
-        if jacobian_determinants.dtype != self.dtype(0).real.dtype:
+        if jacobian_determinants.dtype != geo_dtype:
             raise TypeError("jacobian_determinants has incorrect type")
-        if inverse_jacobians.dtype != self.dtype(0).real.dtype:
+        if inverse_jacobians.dtype != geo_dtype:
             raise TypeError("inverse_jacobians has incorrect type")
 
         gdim = jacobians.shape[1]
@@ -205,18 +220,29 @@ def push_forward(
 
         shape = (pvs,) + reference_values.shape[1:]
         data = np.empty(shape, dtype=self.dtype)
-        _lib.ciarlet_element_push_forward(
+
+        if reference_values.dtype == np.float64:
+            push_function = _lib.ciarlet_element_push_forward_f64
+            geo_type = "double*"
+        elif reference_values.dtype == np.float32:
+            push_function = _lib.ciarlet_element_push_forward_f32
+            geo_type = "float*"
+        else:
+            raise TypeError(f"Unsupported dtype: {reference_values.dtype}")
+
+        push_function(
             self._rs_element,
             npts,
             nfuncs,
             gdim,
-            _ffi.cast("void*", reference_values.ctypes.data),
+            _ffi.cast(geo_type, reference_values.ctypes.data),
             nderivs,
-            _ffi.cast("void*", jacobians.ctypes.data),
-            _ffi.cast("void*", jacobian_determinants.ctypes.data),
-            _ffi.cast("void*", inverse_jacobians.ctypes.data),
+            _ffi.cast(geo_type, jacobians.ctypes.data),
+            _ffi.cast(geo_type, jacobian_determinants.ctypes.data),
+            _ffi.cast(geo_type, inverse_jacobians.ctypes.data),
             _ffi.cast("void*", data.ctypes.data),
         )
+
         return data
 
     def pull_back(
@@ -228,13 +254,12 @@ def pull_back(
         inverse_jacobians: npt.NDArray[np.floating],
     ) -> npt.NDArray[np.floating]:
         """Push values back from a physical cell."""
-        if physical_values.dtype != self.dtype(0):
-            raise TypeError("physical_values has incorrect type")
-        if jacobians.dtype != self.dtype(0).real.dtype:
+        geo_dtype = physical_values.dtype
+        if jacobians.dtype != geo_dtype:
             raise TypeError("jacobians has incorrect type")
-        if jacobian_determinants.dtype != self.dtype(0).real.dtype:
+        if jacobian_determinants.dtype != geo_dtype:
             raise TypeError("jacobian_determinants has incorrect type")
-        if inverse_jacobians.dtype != self.dtype(0).real.dtype:
+        if inverse_jacobians.dtype != geo_dtype:
             raise TypeError("inverse_jacobians has incorrect type")
 
         gdim = jacobians.shape[1]
@@ -244,18 +269,29 @@ def pull_back(
 
         shape = (vs,) + physical_values.shape[1:]
         data = np.empty(shape, dtype=self.dtype)
-        _lib.ciarlet_element_pull_back(
+
+        if physical_values.dtype == np.float64:
+            pull_function = _lib.ciarlet_element_pull_back_f64
+            geo_type = "double*"
+        elif physical_values.dtype == np.float32:
+            pull_function = _lib.ciarlet_element_pull_back_f32
+            geo_type = "float*"
+        else:
+            raise TypeError(f"Unsupported dtype: {physical_values.dtype}")
+
+        pull_function(
             self._rs_element,
             npts,
             nfuncs,
             gdim,
-            _ffi.cast("void*", physical_values.ctypes.data),
+            _ffi.cast(geo_type, physical_values.ctypes.data),
             nderivs,
-            _ffi.cast("void*", jacobians.ctypes.data),
-            _ffi.cast("void*", jacobian_determinants.ctypes.data),
-            _ffi.cast("void*", inverse_jacobians.ctypes.data),
+            _ffi.cast(geo_type, jacobians.ctypes.data),
+            _ffi.cast(geo_type, jacobian_determinants.ctypes.data),
+            _ffi.cast(geo_type, inverse_jacobians.ctypes.data),
             _ffi.cast("void*", data.ctypes.data),
         )
+
         return data
 
     def apply_dof_permutations(self, data: npt.NDArray, orientation: np.int32):

python/test/test_element.py

@@ -27,25 +27,6 @@ def test_value_size(cell, degree):
     assert element.value_size == 1
 
 
-@pytest.mark.parametrize(
-    "dt0,dt1",
-    [
-        (np.float32, np.float64),
-        (np.float64, np.float32),
-        (np.complex64, np.float64),
-        (np.complex128, np.float32),
-    ]
-    + [(dt0, dt1) for dt0 in dtypes for dt1 in [np.complex64, np.complex128]],
-)
-def test_incompatible_types(dt0, dt1):
-    family = create_family(Family.Lagrange, 2, dtype=dt0)
-    element = family.element(ReferenceCellType.Triangle)
-    points = np.array([[0.0, 0.0], [0.2, 0.1], [0.8, 0.05]], dtype=dt1)
-
-    with pytest.raises(TypeError):
-        element.tabulate(points, 1)
-
-
 @pytest.mark.parametrize("dtype", dtypes)
 def test_lagrange_2_triangle_tabulate(dtype):
     family = create_family(Family.Lagrange, 2, dtype=dtype)