diff --git a/examples/element_family.rs b/examples/element_family.rs
index 26874d4..b9b8666 100644
--- a/examples/element_family.rs
+++ b/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);
diff --git a/examples/lagrange_element.rs b/examples/lagrange_element.rs
index b3adce6..59f3ed4 100644
--- a/examples/lagrange_element.rs
+++ b/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());
 
diff --git a/examples/test_high_degree.rs b/examples/test_high_degree.rs
index 776a4bf..5ebb0cf 100644
--- a/examples/test_high_degree.rs
+++ b/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);
                 }
             }
         };
diff --git a/python/ndelement/ciarlet.py b/python/ndelement/ciarlet.py
index 5a492de..8f2306f 100644
--- a/python/ndelement/ciarlet.py
+++ b/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):
diff --git a/python/test/test_element.py b/python/test/test_element.py
index aac526f..5f391dd 100644
--- a/python/test/test_element.py
+++ b/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)
diff --git a/src/bindings.rs b/src/bindings.rs
index fb6110c..68c7c61 100644
--- a/src/bindings.rs
+++ b/src/bindings.rs
@@ -15,7 +15,7 @@ pub mod reference_cell {
     pub unsafe extern "C" fn is_simplex(cell: ReferenceCellType) -> bool {
         reference_cell::is_simplex(cell)
     }
-    unsafe fn vertices<T: RlstScalar<Real = T>>(cell: ReferenceCellType, vs: *mut T) {
+    unsafe fn vertices<T: RlstScalar>(cell: ReferenceCellType, vs: *mut T) {
         let mut i = 0;
         for v in reference_cell::vertices::<T>(cell) {
             for c in v {
@@ -38,7 +38,7 @@ pub mod reference_cell {
             vertices(cell, vs);
         }
     }
-    unsafe fn midpoint<T: RlstScalar<Real = T>>(cell: ReferenceCellType, pt: *mut T) {
+    unsafe fn midpoint<T: RlstScalar>(cell: ReferenceCellType, pt: *mut T) {
         for (i, c) in reference_cell::midpoint(cell).iter().enumerate() {
             unsafe {
                 *pt.add(i) = *c;
@@ -199,16 +199,16 @@ pub mod polynomials {
             *shape.add(2) = npts;
         }
     }
-    unsafe fn tabulate_legendre_polynomials<T: RlstScalar>(
+    unsafe fn tabulate_legendre_polynomials<T: RlstScalar, TGeo: RlstScalar>(
         cell: ReferenceCellType,
-        points: *const T::Real,
+        points: *const TGeo,
         npts: usize,
         degree: usize,
         derivatives: usize,
         data: *mut T,
     ) {
         let tdim = reference_cell::dim(cell);
-        let points = SliceArray::<T::Real, 2>::from_shape(
+        let points = SliceArray::<TGeo, 2>::from_shape(
             unsafe { from_raw_parts(points, npts * tdim) },
             [tdim, npts],
         );
@@ -279,8 +279,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn element_value_size<E: FiniteElement>(element: &E) -> usize {
         element.value_size()
@@ -372,7 +373,12 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
-        field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = ["crate::ciarlet::LagrangeElementFamily<{{dtype}}>", "ciarlet::RaviartThomasElementFamily<{{dtype}}>", "ciarlet::NedelecFirstKindElementFamily<{{dtype}}>"])
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = [
+            "ciarlet::LagrangeElementFamily<{{dtype}}, {{geo_dtype}}>",
+            "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{geo_dtype}}>",
+            "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{geo_dtype}}>"
+        ])
     )]
     pub fn element_family_create_element<F: ElementFamily<CellType = ReferenceCellType>>(
         family: &F,
@@ -388,7 +394,12 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
-        field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = ["crate::ciarlet::LagrangeElementFamily<{{dtype}}>", "ciarlet::RaviartThomasElementFamily<{{dtype}}>", "ciarlet::NedelecFirstKindElementFamily<{{dtype}}>"])
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = [
+            "ciarlet::LagrangeElementFamily<{{dtype}}, {{geo_dtype}}>",
+            "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{geo_dtype}}>",
+            "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{geo_dtype}}>"
+        ])
     )]
     pub fn element_family_dtype<
         T: RlstScalar + DTypeIdentifier,
@@ -401,8 +412,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_dtype<T: RlstScalar + DTypeIdentifier, E: FiniteElement<T = T>>(
         _elem: &E,
@@ -412,8 +424,25 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
+    )]
+    pub fn ciarlet_element_geo_dtype<
+        T: RlstScalar + DTypeIdentifier,
+        M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
+    >(
+        _elem: &CiarletElement<T, M, TGeo>,
+    ) -> DType {
+        <TGeo as DTypeIdentifier>::dtype()
+    }
+
+    #[concretise_types(
+        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_tabulate_array_shape<E: FiniteElement>(
         element: &E,
@@ -432,22 +461,19 @@ pub mod ciarlet {
         }
     }
 
-    #[concretise_types(
-        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
-        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
-    )]
-    pub fn ciarlet_element_tabulate<E: FiniteElement<CellType = ReferenceCellType>>(
+    pub unsafe fn ciarlet_element_tabulate<
+        E: FiniteElement<CellType = ReferenceCellType>,
+        TGeo: RlstScalar,
+    >(
         element: &E,
-        points: *const c_void,
+        points: *const TGeo,
         npoints: usize,
         nderivs: usize,
         data: *mut c_void,
     ) {
         let tdim = reference_cell::dim(element.cell_type());
-        let points = points as *mut <E::T as RlstScalar>::Real;
         let data = data as *mut E::T;
-        let points = SliceArray::<<E::T as RlstScalar>::Real, 2>::from_shape(
+        let points = SliceArray::<TGeo, 2>::from_shape(
             unsafe { from_raw_parts(points, npoints * tdim) },
             [tdim, npoints],
         );
@@ -461,8 +487,45 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]),
+    )]
+    pub unsafe fn ciarlet_element_tabulate_f32<E: FiniteElement<CellType = ReferenceCellType>>(
+        element: &E,
+        points: *const f32,
+        npoints: usize,
+        nderivs: usize,
+        data: *mut c_void,
+    ) {
+        unsafe {
+            ciarlet_element_tabulate(element, points, npoints, nderivs, data);
+        }
+    }
+
+    #[concretise_types(
+        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]),
+    )]
+    pub unsafe fn ciarlet_element_tabulate_f64<E: FiniteElement<CellType = ReferenceCellType>>(
+        element: &E,
+        points: *const f64,
+        npoints: usize,
+        nderivs: usize,
+        data: *mut c_void,
+    ) {
+        unsafe {
+            ciarlet_element_tabulate(element, points, npoints, nderivs, data);
+        }
+    }
+
+    #[concretise_types(
+        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_value_size<E: FiniteElement>(element: &E) -> usize {
         element.value_size()
@@ -470,8 +533,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_value_rank<E: FiniteElement>(element: &E) -> usize {
         element.value_shape().len()
@@ -479,8 +543,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_value_shape<E: FiniteElement>(element: &E, shape: *mut usize) {
         for (i, j) in element.value_shape().iter().enumerate() {
@@ -492,8 +557,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_physical_value_size<E: MappedFiniteElement>(
         element: &E,
@@ -504,8 +570,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_physical_value_rank<E: MappedFiniteElement>(
         element: &E,
@@ -516,8 +583,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_physical_value_shape<E: MappedFiniteElement>(
         element: &E,
@@ -532,21 +600,19 @@ pub mod ciarlet {
     }
 
     #[allow(clippy::too_many_arguments)]
-    #[concretise_types(
-        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
-        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
-    )]
-    pub fn ciarlet_element_push_forward<E: MappedFiniteElement<CellType = ReferenceCellType>>(
+    pub unsafe fn ciarlet_element_push_forward<
+        E: MappedFiniteElement<CellType = ReferenceCellType>,
+        TGeo: RlstScalar,
+    >(
         element: &E,
         npoints: usize,
         nfunctions: usize,
         gdim: usize,
-        reference_values: *const c_void,
+        reference_values: *const TGeo,
         nderivs: usize,
-        j: *const c_void,
-        jdet: *const c_void,
-        jinv: *const c_void,
+        j: *const TGeo,
+        jdet: *const TGeo,
+        jinv: *const TGeo,
         physical_values: *mut c_void,
     ) {
         let tdim = reference_cell::dim(element.cell_type());
@@ -562,23 +628,13 @@ pub mod ciarlet {
             },
             [deriv_size, npoints, nfunctions, vs],
         );
-        let j = SliceArray::<<E::T as RlstScalar>::Real, 3>::from_shape(
-            unsafe {
-                from_raw_parts(
-                    j as *const <E::T as RlstScalar>::Real,
-                    npoints * gdim * tdim,
-                )
-            },
+        let j = SliceArray::<TGeo, 3>::from_shape(
+            unsafe { from_raw_parts(j, npoints * gdim * tdim) },
             [npoints, gdim, tdim],
         );
-        let jdet = unsafe { from_raw_parts(jdet as *const <E::T as RlstScalar>::Real, npoints) };
-        let jinv = SliceArray::<<E::T as RlstScalar>::Real, 3>::from_shape(
-            unsafe {
-                from_raw_parts(
-                    jinv as *const <E::T as RlstScalar>::Real,
-                    npoints * tdim * gdim,
-                )
-            },
+        let jdet = unsafe { from_raw_parts(jdet, npoints) };
+        let jinv = SliceArray::<TGeo, 3>::from_shape(
+            unsafe { from_raw_parts(jinv, npoints * tdim * gdim) },
             [npoints, tdim, gdim],
         );
         let mut physical_values = SliceArrayMut::<E::T, 4>::from_shape(
@@ -603,19 +659,91 @@ pub mod ciarlet {
     #[allow(clippy::too_many_arguments)]
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
+    )]
+    pub unsafe fn ciarlet_element_push_forward_f32<
+        E: MappedFiniteElement<CellType = ReferenceCellType>,
+    >(
+        element: &E,
+        npoints: usize,
+        nfunctions: usize,
+        gdim: usize,
+        reference_values: *const f32,
+        nderivs: usize,
+        j: *const f32,
+        jdet: *const f32,
+        jinv: *const f32,
+        physical_values: *mut c_void,
+    ) {
+        unsafe {
+            ciarlet_element_push_forward(
+                element,
+                npoints,
+                nfunctions,
+                gdim,
+                reference_values,
+                nderivs,
+                j,
+                jdet,
+                jinv,
+                physical_values,
+            );
+        }
+    }
+
+    #[allow(clippy::too_many_arguments)]
+    #[concretise_types(
+        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
-    pub fn ciarlet_element_pull_back<E: MappedFiniteElement<CellType = ReferenceCellType>>(
+    pub unsafe fn ciarlet_element_push_forward_f64<
+        E: MappedFiniteElement<CellType = ReferenceCellType>,
+    >(
+        element: &E,
+        npoints: usize,
+        nfunctions: usize,
+        gdim: usize,
+        reference_values: *const f64,
+        nderivs: usize,
+        j: *const f64,
+        jdet: *const f64,
+        jinv: *const f64,
+        physical_values: *mut c_void,
+    ) {
+        unsafe {
+            ciarlet_element_push_forward(
+                element,
+                npoints,
+                nfunctions,
+                gdim,
+                reference_values,
+                nderivs,
+                j,
+                jdet,
+                jinv,
+                physical_values,
+            );
+        }
+    }
+
+    #[allow(clippy::too_many_arguments)]
+    pub unsafe fn ciarlet_element_pull_back<
+        E: MappedFiniteElement<CellType = ReferenceCellType>,
+        TGeo: RlstScalar,
+    >(
         element: &E,
         npoints: usize,
         nfunctions: usize,
         gdim: usize,
-        physical_values: *const c_void,
+        physical_values: *const TGeo,
         nderivs: usize,
-        j: *const c_void,
-        jdet: *const c_void,
-        jinv: *const c_void,
+        j: *const TGeo,
+        jdet: *const TGeo,
+        jinv: *const TGeo,
         reference_values: *mut c_void,
     ) {
         let tdim = reference_cell::dim(element.cell_type());
@@ -631,23 +759,13 @@ pub mod ciarlet {
             },
             [deriv_size, npoints, nfunctions, pvs],
         );
-        let j = SliceArray::<<E::T as RlstScalar>::Real, 3>::from_shape(
-            unsafe {
-                from_raw_parts(
-                    j as *const <E::T as RlstScalar>::Real,
-                    npoints * gdim * tdim,
-                )
-            },
+        let j = SliceArray::<TGeo, 3>::from_shape(
+            unsafe { from_raw_parts(j, npoints * gdim * tdim) },
             [npoints, gdim, tdim],
         );
-        let jdet = unsafe { from_raw_parts(jdet as *const <E::T as RlstScalar>::Real, npoints) };
-        let jinv = SliceArray::<<E::T as RlstScalar>::Real, 3>::from_shape(
-            unsafe {
-                from_raw_parts(
-                    jinv as *const <E::T as RlstScalar>::Real,
-                    npoints * tdim * gdim,
-                )
-            },
+        let jdet = unsafe { from_raw_parts(jdet, npoints) };
+        let jinv = SliceArray::<TGeo, 3>::from_shape(
+            unsafe { from_raw_parts(jinv, npoints * tdim * gdim) },
             [npoints, tdim, gdim],
         );
         let mut reference_values = SliceArrayMut::<E::T, 4>::from_shape(
@@ -669,21 +787,99 @@ pub mod ciarlet {
         );
     }
 
+    #[allow(clippy::too_many_arguments)]
+    #[concretise_types(
+        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
+    )]
+    pub unsafe fn ciarlet_element_pull_back_f32<
+        E: MappedFiniteElement<CellType = ReferenceCellType>,
+    >(
+        element: &E,
+        npoints: usize,
+        nfunctions: usize,
+        gdim: usize,
+        physical_values: *const f32,
+        nderivs: usize,
+        j: *const f32,
+        jdet: *const f32,
+        jinv: *const f32,
+        reference_values: *mut c_void,
+    ) {
+        unsafe {
+            ciarlet_element_pull_back(
+                element,
+                npoints,
+                nfunctions,
+                gdim,
+                physical_values,
+                nderivs,
+                j,
+                jdet,
+                jinv,
+                reference_values,
+            );
+        }
+    }
+
+    #[allow(clippy::too_many_arguments)]
+    #[concretise_types(
+        gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
+    )]
+    pub fn ciarlet_element_pull_back_f64<E: MappedFiniteElement<CellType = ReferenceCellType>>(
+        element: &E,
+        npoints: usize,
+        nfunctions: usize,
+        gdim: usize,
+        physical_values: *const f64,
+        nderivs: usize,
+        j: *const f64,
+        jdet: *const f64,
+        jinv: *const f64,
+        reference_values: *mut c_void,
+    ) {
+        unsafe {
+            ciarlet_element_pull_back(
+                element,
+                npoints,
+                nfunctions,
+                gdim,
+                physical_values,
+                nderivs,
+                j,
+                jdet,
+                jinv,
+                reference_values,
+            );
+        }
+    }
+
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
-    pub fn ciarlet_element_degree<T: RlstScalar + DTypeIdentifier + Getrf + Getri, M: Map>(
-        element: &CiarletElement<T, M>,
+    pub fn ciarlet_element_degree<
+        T: RlstScalar + DTypeIdentifier + Getrf + Getri,
+        M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
+    >(
+        element: &CiarletElement<T, M, TGeo>,
     ) -> usize {
         element.degree()
     }
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_embedded_superdegree<E: MappedFiniteElement>(element: &E) -> usize {
         element.lagrange_superdegree()
@@ -691,8 +887,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_dim<E: FiniteElement>(element: &E) -> usize {
         element.dim()
@@ -700,19 +897,25 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
-    pub fn ciarlet_element_continuity<T: RlstScalar + DTypeIdentifier + Getrf + Getri, M: Map>(
-        element: &CiarletElement<T, M>,
+    pub fn ciarlet_element_continuity<
+        T: RlstScalar + DTypeIdentifier + Getrf + Getri,
+        M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
+    >(
+        element: &CiarletElement<T, M, TGeo>,
     ) -> Continuity {
         element.continuity()
     }
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_cell_type<E: FiniteElement<CellType = ReferenceCellType>>(
         element: &E,
@@ -722,8 +925,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_entity_dofs_size<E: FiniteElement>(
         element: &E,
@@ -735,8 +939,9 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_entity_dofs<E: FiniteElement>(
         element: &E,
@@ -758,11 +963,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
-    pub fn ciarlet_element_entity_closure_dofs_size<T: RlstScalar + DTypeIdentifier, M: Map>(
-        element: &CiarletElement<T, M>,
+    pub fn ciarlet_element_entity_closure_dofs_size<
+        T: RlstScalar + DTypeIdentifier,
+        M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
+    >(
+        element: &CiarletElement<T, M, TGeo>,
         entity_dim: usize,
         entity_index: usize,
     ) -> usize {
@@ -774,11 +984,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
-    pub fn ciarlet_element_entity_closure_dofs<T: RlstScalar + DTypeIdentifier, M: Map>(
-        element: &CiarletElement<T, M>,
+    pub fn ciarlet_element_entity_closure_dofs<
+        T: RlstScalar + DTypeIdentifier,
+        M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
+    >(
+        element: &CiarletElement<T, M, TGeo>,
         entity_dim: usize,
         entity_index: usize,
         entity_closure_dofs: *mut usize,
@@ -797,14 +1012,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_interpolation_npoints<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         entity_dim: usize,
         entity_index: usize,
     ) -> usize {
@@ -813,14 +1030,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_interpolation_ndofs<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         entity_dim: usize,
         entity_index: usize,
     ) -> usize {
@@ -829,19 +1048,21 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_interpolation_points<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
+        TGeo: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         entity_dim: usize,
         entity_index: usize,
         points: *mut c_void,
     ) {
-        let points = points as *mut T::Real;
+        let points = points as *mut TGeo;
         for (i, j) in element.interpolation_points()[entity_dim][entity_index]
             .data()
             .unwrap()
@@ -856,14 +1077,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub fn ciarlet_element_interpolation_weights<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         entity_dim: usize,
         entity_index: usize,
         weights: *mut c_void,
@@ -883,14 +1106,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub unsafe fn ciarlet_element_apply_dof_permutations_usize<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         data: *mut usize,
         data_size: usize,
         orientation: i32,
@@ -902,14 +1127,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub unsafe fn ciarlet_element_apply_dof_permutations<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         data: *mut c_void,
         data_size: usize,
         orientation: i32,
@@ -922,14 +1149,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub unsafe fn ciarlet_element_apply_dof_transformations<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         data: *mut c_void,
         data_size: usize,
         orientation: i32,
@@ -944,14 +1173,16 @@ pub mod ciarlet {
 
     #[concretise_types(
         gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]),
+        gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]),
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
-        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
+        field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"])
     )]
     pub unsafe fn ciarlet_element_apply_dof_permutations_and_transformations<
         T: RlstScalar + DTypeIdentifier + Getrf + Getri,
         M: Map,
+        TGeo: RlstScalar + DTypeIdentifier,
     >(
-        element: &CiarletElement<T, M>,
+        element: &CiarletElement<T, M, TGeo>,
         data: *mut c_void,
         data_size: usize,
         orientation: i32,
diff --git a/src/ciarlet.rs b/src/ciarlet.rs
index fcb702f..38e2701 100644
--- a/src/ciarlet.rs
+++ b/src/ciarlet.rs
@@ -35,7 +35,6 @@ use crate::reference_cell;
 use crate::traits::{FiniteElement, Map, MappedFiniteElement};
 use crate::types::{Continuity, DofTransformation, ReferenceCellType, Transformation};
 use itertools::izip;
-use num::{One, Zero};
 use rlst::RlstScalar;
 use rlst::{Array, DynArray, Inverse, MutableArrayImpl, ValueArrayImpl, rlst_dynamic_array};
 use std::collections::HashMap;
@@ -66,7 +65,7 @@ fn compute_derivative_count(nderivs: usize, cell_type: ReferenceCellType) -> usi
 }
 
 /// A Ciarlet element
-pub struct CiarletElement<T: RlstScalar, M: Map> {
+pub struct CiarletElement<T: RlstScalar, M: Map, TGeo: RlstScalar = <T as RlstScalar>::Real> {
     family_name: String,
     cell_type: ReferenceCellType,
     degree: usize,
@@ -78,13 +77,13 @@ pub struct CiarletElement<T: RlstScalar, M: Map> {
     coefficients: DynArray<T, 3>,
     entity_dofs: [Vec<Vec<usize>>; 4],
     entity_closure_dofs: [Vec<Vec<usize>>; 4],
-    interpolation_points: EntityPoints<T::Real>,
+    interpolation_points: EntityPoints<TGeo>,
     interpolation_weights: EntityWeights<T>,
     map: M,
     dof_transformations: HashMap<(ReferenceCellType, Transformation), DofTransformation<T>>,
 }
 
-impl<T: RlstScalar, M: Map> Debug for CiarletElement<T, M> {
+impl<T: RlstScalar, M: Map, TGeo: RlstScalar> Debug for CiarletElement<T, M, TGeo> {
     fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), std::fmt::Error> {
         f.debug_tuple("CiarletElement")
             .field(&self.family_name)
@@ -94,7 +93,7 @@ impl<T: RlstScalar, M: Map> Debug for CiarletElement<T, M> {
     }
 }
 
-impl<T: RlstScalar, M: Map> CiarletElement<T, M>
+impl<T: RlstScalar, M: Map, TGeo: RlstScalar> CiarletElement<T, M, TGeo>
 where
     DynArray<T, 2>: Inverse<Output = DynArray<T, 2>>,
 {
@@ -112,7 +111,7 @@ where
         degree: usize,
         value_shape: Vec<usize>,
         polynomial_coeffs: DynArray<T, 3>,
-        interpolation_points: EntityPoints<T::Real>,
+        interpolation_points: EntityPoints<TGeo>,
         interpolation_weights: EntityWeights<T>,
         continuity: Continuity,
         embedded_superdegree: usize,
@@ -144,11 +143,11 @@ where
 
         // Format the interpolation points and weights
         let new_pts = if continuity == Continuity::Discontinuous {
-            let mut new_pts: EntityPoints<T::Real> = [vec![], vec![], vec![], vec![]];
-            let mut all_pts = rlst_dynamic_array!(T::Real, [tdim, npts]);
+            let mut new_pts: EntityPoints<TGeo> = [vec![], vec![], vec![], vec![]];
+            let mut all_pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
             for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() {
                 for _pts in pts_i {
-                    new_pts[i].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+                    new_pts[i].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
                 }
             }
             let mut col = 0;
@@ -303,32 +302,32 @@ where
                 ReferenceCellType::Interval,
                 0,
                 Transformation::Reflection,
-                (|x| vec![x[1], x[0]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                (|x| vec![x[1], x[0]]) as fn(&[TGeo]) -> Vec<TGeo>,
             )],
             ReferenceCellType::Quadrilateral => vec![(
                 ReferenceCellType::Interval,
                 0,
                 Transformation::Reflection,
-                (|x| vec![T::Real::one() - x[0], x[1]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                (|x| vec![TGeo::one() - x[0], x[1]]) as fn(&[TGeo]) -> Vec<TGeo>,
             )],
             ReferenceCellType::Tetrahedron => vec![
                 (
                     ReferenceCellType::Interval,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![x[0], x[2], x[1]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[0], x[2], x[1]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Triangle,
                     0,
                     Transformation::Rotation,
-                    (|x| vec![x[1], x[2], x[0]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], x[2], x[0]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Triangle,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![x[0], x[2], x[1]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[0], x[2], x[1]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
             ],
             ReferenceCellType::Hexahedron => vec![
@@ -336,19 +335,19 @@ where
                     ReferenceCellType::Interval,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![T::Real::one() - x[0], x[1], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![TGeo::one() - x[0], x[1], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Quadrilateral,
                     0,
                     Transformation::Rotation,
-                    (|x| vec![x[1], T::Real::one() - x[0], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], TGeo::one() - x[0], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Quadrilateral,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![x[1], x[0], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], x[0], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
             ],
             ReferenceCellType::Prism => vec![
@@ -356,32 +355,31 @@ where
                     ReferenceCellType::Interval,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![T::Real::one() - x[0], x[1], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![TGeo::one() - x[0], x[1], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Triangle,
                     0,
                     Transformation::Rotation,
-                    (|x| vec![x[1], T::Real::one() - x[1] - x[0], x[2]])
-                        as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], TGeo::one() - x[1] - x[0], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Triangle,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![x[1], x[0], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], x[0], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Quadrilateral,
                     1,
                     Transformation::Rotation,
-                    (|x| vec![x[2], T::Real::one() - x[1], x[0]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[2], TGeo::one() - x[1], x[0]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Quadrilateral,
                     1,
                     Transformation::Reflection,
-                    (|x| vec![x[2], x[1], x[0]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[2], x[1], x[0]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
             ],
             ReferenceCellType::Pyramid => vec![
@@ -389,32 +387,31 @@ where
                     ReferenceCellType::Interval,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![T::Real::one() - x[0], x[1], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![TGeo::one() - x[0], x[1], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Triangle,
                     1,
                     Transformation::Rotation,
-                    (|x| vec![x[2], x[1], T::Real::one() - x[2] - x[0]])
-                        as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[2], x[1], TGeo::one() - x[2] - x[0]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Triangle,
                     1,
                     Transformation::Reflection,
-                    (|x| vec![x[2], x[1], x[0]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[2], x[1], x[0]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Quadrilateral,
                     0,
                     Transformation::Rotation,
-                    (|x| vec![x[1], T::Real::one() - x[0], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], TGeo::one() - x[0], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
                 (
                     ReferenceCellType::Quadrilateral,
                     0,
                     Transformation::Reflection,
-                    (|x| vec![x[1], x[0], x[2]]) as fn(&[T::Real]) -> Vec<T::Real>,
+                    (|x| vec![x[1], x[0], x[2]]) as fn(&[TGeo]) -> Vec<TGeo>,
                 ),
             ],
         } {
@@ -424,26 +421,23 @@ where
 
             let finv = match transform {
                 Transformation::Reflection => {
-                    (|x, f| f(x)) as fn(&[T::Real], fn(&[T::Real]) -> Vec<T::Real>) -> Vec<T::Real>
+                    (|x, f| f(x)) as fn(&[TGeo], fn(&[TGeo]) -> Vec<TGeo>) -> Vec<TGeo>
                 }
                 Transformation::Rotation => match entity {
                     ReferenceCellType::Interval => {
-                        (|x, f| f(x))
-                            as fn(&[T::Real], fn(&[T::Real]) -> Vec<T::Real>) -> Vec<T::Real>
+                        (|x, f| f(x)) as fn(&[TGeo], fn(&[TGeo]) -> Vec<TGeo>) -> Vec<TGeo>
                     }
                     ReferenceCellType::Triangle => {
-                        (|x, f| f(&f(x)))
-                            as fn(&[T::Real], fn(&[T::Real]) -> Vec<T::Real>) -> Vec<T::Real>
+                        (|x, f| f(&f(x))) as fn(&[TGeo], fn(&[TGeo]) -> Vec<TGeo>) -> Vec<TGeo>
                     }
                     ReferenceCellType::Quadrilateral => {
-                        (|x, f| f(&f(&f(x))))
-                            as fn(&[T::Real], fn(&[T::Real]) -> Vec<T::Real>) -> Vec<T::Real>
+                        (|x, f| f(&f(&f(x)))) as fn(&[TGeo], fn(&[TGeo]) -> Vec<TGeo>) -> Vec<TGeo>
                     }
                     _ => panic!("Unsupported entity: {entity:?}"),
                 },
             };
 
-            let mut pts = DynArray::<T::Real, 2>::from_shape(ref_pts.shape());
+            let mut pts = DynArray::<TGeo, 2>::from_shape(ref_pts.shape());
             for p in 0..npts {
                 for (i, c) in finv(
                     &ref_pts.data().unwrap()[ref_pts.shape()[0] * p..ref_pts.shape()[0] * (p + 1)],
@@ -459,13 +453,13 @@ where
             let wts = &new_wts[edim][entity_id];
             let edofs = &entity_dofs[edim][entity_id];
             let endofs = edofs.len();
-            let mut j = rlst_dynamic_array![T::Real, [npts, tdim, tdim]];
+            let mut j = rlst_dynamic_array![TGeo, [npts, tdim, tdim]];
             for t_in in 0..tdim {
                 for (t_out, (a, b)) in izip!(
-                    f(&vec![T::Real::zero(); tdim]),
+                    f(&vec![TGeo::zero(); tdim]),
                     f(&{
-                        let mut axis = vec![T::Real::zero(); tdim];
-                        axis[t_in] = T::Real::one();
+                        let mut axis = vec![TGeo::zero(); tdim];
+                        axis[t_in] = TGeo::one();
                         axis
                     })
                 )
@@ -479,19 +473,19 @@ where
             // f is linear. Use this to compute determinants
             let jdet = vec![
                 match transform {
-                    Transformation::Reflection => -T::Real::one(),
-                    Transformation::Rotation => T::Real::one(),
+                    Transformation::Reflection => -TGeo::one(),
+                    Transformation::Rotation => TGeo::one(),
                 };
                 npts
             ];
-            let mut jinv = rlst_dynamic_array![T::Real, [npts, tdim, tdim]];
+            let mut jinv = rlst_dynamic_array![TGeo, [npts, tdim, tdim]];
             for t_in in 0..tdim {
                 for (t_out, (a, b)) in izip!(
-                    finv(&vec![T::Real::zero(); tdim], f),
+                    finv(&vec![TGeo::zero(); tdim], f),
                     finv(
                         &{
-                            let mut axis = vec![T::Real::zero(); tdim];
-                            axis[t_in] = T::Real::one();
+                            let mut axis = vec![TGeo::zero(); tdim];
+                            axis[t_in] = TGeo::one();
                             axis
                         },
                         f
@@ -500,7 +494,7 @@ where
                 .enumerate()
                 {
                     for p in 0..npts {
-                        *jinv.get_mut([p, t_out, t_in]).unwrap() = (b - a).re();
+                        *jinv.get_mut([p, t_out, t_in]).unwrap() = TGeo::from(b - a).unwrap();
                     }
                 }
             }
@@ -585,7 +579,7 @@ where
                 );
             }
         }
-        CiarletElement::<T, M> {
+        CiarletElement::<T, M, TGeo> {
             family_name,
             cell_type,
             degree,
@@ -613,7 +607,7 @@ where
         self.continuity
     }
     /// The interpolation points
-    pub fn interpolation_points(&self) -> &EntityPoints<T::Real> {
+    pub fn interpolation_points(&self) -> &EntityPoints<TGeo> {
         &self.interpolation_points
     }
     /// The interpolation weights
@@ -622,7 +616,9 @@ where
     }
 }
 
-impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
+impl<T: RlstScalar, M: Map, TGeoInternal: RlstScalar> FiniteElement
+    for CiarletElement<T, M, TGeoInternal>
+{
     type CellType = ReferenceCellType;
     type T = T;
 
@@ -643,8 +639,9 @@ impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
     }
 
     fn tabulate<
-        Array2Impl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 2>,
-        Array4MutImpl: MutableArrayImpl<Self::T, 4>,
+        TGeo: RlstScalar,
+        Array2Impl: ValueArrayImpl<TGeo, 2>,
+        Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         points: &Array<Array2Impl, 2>,
@@ -708,7 +705,9 @@ impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
     }
 }
 
-impl<T: RlstScalar, M: Map> MappedFiniteElement for CiarletElement<T, M> {
+impl<T: RlstScalar, M: Map, TGeoInternal: RlstScalar> MappedFiniteElement
+    for CiarletElement<T, M, TGeoInternal>
+{
     type TransformationType = Transformation;
 
     fn lagrange_superdegree(&self) -> usize {
@@ -716,16 +715,17 @@ impl<T: RlstScalar, M: Map> MappedFiniteElement for CiarletElement<T, M> {
     }
 
     fn push_forward<
-        Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
-        Array4Impl: ValueArrayImpl<Self::T, 4>,
-        Array4MutImpl: MutableArrayImpl<Self::T, 4>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
+        Array4Impl: ValueArrayImpl<T, 4>,
+        Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         reference_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[<Self::T as RlstScalar>::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         physical_values: &mut Array<Array4MutImpl, 4>,
     ) {
         self.map.push_forward(
@@ -739,16 +739,17 @@ impl<T: RlstScalar, M: Map> MappedFiniteElement for CiarletElement<T, M> {
     }
 
     fn pull_back<
-        Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
-        Array4Impl: ValueArrayImpl<Self::T, 4>,
-        Array4MutImpl: MutableArrayImpl<Self::T, 4>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
+        Array4Impl: ValueArrayImpl<T, 4>,
+        Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         physical_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[<Self::T as RlstScalar>::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         reference_values: &mut Array<Array4MutImpl, 4>,
     ) {
         self.map.pull_back(
@@ -847,7 +848,7 @@ impl<T: RlstScalar, M: Map> MappedFiniteElement for CiarletElement<T, M> {
         }
     }
 
-    fn apply_dof_transformations(&self, data: &mut [Self::T], cell_orientation: i32) {
+    fn apply_dof_transformations(&self, data: &mut [T], cell_orientation: i32) {
         debug_assert_eq!(data.len() % self.dim, 0);
         let block_size = data.len() / self.dim;
         let tdim = reference_cell::dim(self.cell_type);
@@ -950,13 +951,14 @@ mod test {
 
     #[test]
     fn test_lagrange_1() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
+        let e = lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
     }
 
     #[test]
     fn test_lagrange_0_interval() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Interval, 0, Continuity::Discontinuous);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Interval, 0, Continuity::Discontinuous);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
         let mut points = rlst_dynamic_array!(f64, [1, 4]);
@@ -974,7 +976,7 @@ mod test {
 
     #[test]
     fn test_lagrange_1_interval() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Interval, 1, Continuity::Standard);
+        let e = lagrange::create::<f64, f64>(ReferenceCellType::Interval, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
         let mut points = rlst_dynamic_array!(f64, [1, 4]);
@@ -994,7 +996,8 @@ mod test {
 
     #[test]
     fn test_lagrange_0_triangle() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Triangle, 0, Continuity::Discontinuous);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 0, Continuity::Discontinuous);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
 
@@ -1022,7 +1025,7 @@ mod test {
 
     #[test]
     fn test_lagrange_1_triangle() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
+        let e = lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1052,7 +1055,7 @@ mod test {
 
     #[test]
     fn test_lagrange_0_quadrilateral() {
-        let e = lagrange::create::<f64>(
+        let e = lagrange::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             0,
             Continuity::Discontinuous,
@@ -1082,7 +1085,8 @@ mod test {
 
     #[test]
     fn test_lagrange_1_quadrilateral() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1114,7 +1118,8 @@ mod test {
 
     #[test]
     fn test_lagrange_2_quadrilateral() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1187,8 +1192,11 @@ mod test {
 
     #[test]
     fn test_lagrange_0_tetrahedron() {
-        let e =
-            lagrange::create::<f64>(ReferenceCellType::Tetrahedron, 0, Continuity::Discontinuous);
+        let e = lagrange::create::<f64, f64>(
+            ReferenceCellType::Tetrahedron,
+            0,
+            Continuity::Discontinuous,
+        );
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [3, 6]);
@@ -1220,7 +1228,8 @@ mod test {
 
     #[test]
     fn test_lagrange_1_tetrahedron() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [3, 6]);
@@ -1263,8 +1272,11 @@ mod test {
 
     #[test]
     fn test_lagrange_0_hexahedron() {
-        let e =
-            lagrange::create::<f64>(ReferenceCellType::Hexahedron, 0, Continuity::Discontinuous);
+        let e = lagrange::create::<f64, f64>(
+            ReferenceCellType::Hexahedron,
+            0,
+            Continuity::Discontinuous,
+        );
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [3, 6]);
@@ -1296,7 +1308,8 @@ mod test {
 
     #[test]
     fn test_lagrange_1_hexahedron() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [3, 6]);
@@ -1371,53 +1384,53 @@ mod test {
 
     #[test]
     fn test_lagrange_higher_degree_triangle() {
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 5, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 5, Continuity::Standard);
 
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Discontinuous);
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Discontinuous);
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 4, Continuity::Discontinuous);
-        lagrange::create::<f64>(ReferenceCellType::Triangle, 5, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 2, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 3, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 5, Continuity::Discontinuous);
     }
 
     #[test]
     fn test_lagrange_higher_degree_interval() {
-        lagrange::create::<f64>(ReferenceCellType::Interval, 2, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Interval, 3, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Interval, 4, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Interval, 5, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 2, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 3, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 4, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 5, Continuity::Standard);
 
-        lagrange::create::<f64>(ReferenceCellType::Interval, 2, Continuity::Discontinuous);
-        lagrange::create::<f64>(ReferenceCellType::Interval, 3, Continuity::Discontinuous);
-        lagrange::create::<f64>(ReferenceCellType::Interval, 4, Continuity::Discontinuous);
-        lagrange::create::<f64>(ReferenceCellType::Interval, 5, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 2, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 3, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 4, Continuity::Discontinuous);
+        lagrange::create::<f64, f64>(ReferenceCellType::Interval, 5, Continuity::Discontinuous);
     }
 
     #[test]
     fn test_lagrange_higher_degree_quadrilateral() {
-        lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 3, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 4, Continuity::Standard);
-        lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 5, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Quadrilateral, 3, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Quadrilateral, 4, Continuity::Standard);
+        lagrange::create::<f64, f64>(ReferenceCellType::Quadrilateral, 5, Continuity::Standard);
 
-        lagrange::create::<f64>(
+        lagrange::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             2,
             Continuity::Discontinuous,
         );
-        lagrange::create::<f64>(
+        lagrange::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             3,
             Continuity::Discontinuous,
         );
-        lagrange::create::<f64>(
+        lagrange::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             4,
             Continuity::Discontinuous,
         );
-        lagrange::create::<f64>(
+        lagrange::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             5,
             Continuity::Discontinuous,
@@ -1426,7 +1439,11 @@ mod test {
 
     #[test]
     fn test_raviart_thomas_1_triangle() {
-        let e = raviart_thomas::create(ReferenceCellType::Triangle, 1, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Triangle,
+            1,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 2);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1458,21 +1475,33 @@ mod test {
 
     #[test]
     fn test_raviart_thomas_2_triangle() {
-        let e = raviart_thomas::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Triangle,
+            2,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 2);
         check_dofs(e);
     }
 
     #[test]
     fn test_raviart_thomas_3_triangle() {
-        let e = raviart_thomas::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Triangle,
+            3,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 2);
         check_dofs(e);
     }
 
     #[test]
     fn test_raviart_thomas_1_quadrilateral() {
-        let e = raviart_thomas::create(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Quadrilateral,
+            1,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 2);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1507,7 +1536,7 @@ mod test {
 
     #[test]
     fn test_raviart_thomas_2_quadrilateral() {
-        let e = raviart_thomas::create::<f64>(
+        let e = raviart_thomas::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             2,
             Continuity::Standard,
@@ -1518,7 +1547,7 @@ mod test {
 
     #[test]
     fn test_raviart_thomas_3_quadrilateral() {
-        let e = raviart_thomas::create::<f64>(
+        let e = raviart_thomas::create::<f64, f64>(
             ReferenceCellType::Quadrilateral,
             3,
             Continuity::Standard,
@@ -1529,31 +1558,44 @@ mod test {
 
     #[test]
     fn test_raviart_thomas_1_tetrahedron() {
-        let e =
-            raviart_thomas::create::<f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Tetrahedron,
+            1,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_raviart_thomas_2_tetrahedron() {
-        let e =
-            raviart_thomas::create::<f64>(ReferenceCellType::Tetrahedron, 2, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Tetrahedron,
+            2,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_raviart_thomas_3_tetrahedron() {
-        let e =
-            raviart_thomas::create::<f64>(ReferenceCellType::Tetrahedron, 3, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Tetrahedron,
+            3,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_raviart_thomas_1_hexahedron() {
-        let e = raviart_thomas::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Hexahedron,
+            1,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 3);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [3, 6]);
@@ -1602,15 +1644,18 @@ mod test {
 
     #[test]
     fn test_raviart_thomas_2_hexahedron() {
-        let e =
-            raviart_thomas::create::<f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
+        let e = raviart_thomas::create::<f64, f64>(
+            ReferenceCellType::Hexahedron,
+            2,
+            Continuity::Standard,
+        );
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_1_triangle() {
-        let e = nedelec::create(ReferenceCellType::Triangle, 1, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1642,21 +1687,22 @@ mod test {
 
     #[test]
     fn test_nedelec_2_triangle() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_3_triangle() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_1_quadrilateral() {
-        let e = nedelec::create(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [2, 6]);
@@ -1691,42 +1737,47 @@ mod test {
 
     #[test]
     fn test_nedelec_2_quadrilateral() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_3_quadrilateral() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Quadrilateral, 3, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Quadrilateral, 3, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_1_tetrahedron() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_2_tetrahedron() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Tetrahedron, 2, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Tetrahedron, 2, Continuity::Standard);
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_3_tetrahedron() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Tetrahedron, 3, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Tetrahedron, 3, Continuity::Standard);
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
 
     #[test]
     fn test_nedelec_1_hexahedron() {
-        let e = nedelec::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 3);
         let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
         let mut points = rlst_dynamic_array!(f64, [3, 6]);
@@ -1781,7 +1832,7 @@ mod test {
 
     #[test]
     fn test_nedelec_2_hexahedron() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
         assert_eq!(e.value_size(), 3);
         check_dofs(e);
     }
@@ -1791,7 +1842,7 @@ mod test {
             paste! {
                 #[test]
                 fn [<test_entity_closure_dofs_ $cell:lower _ $degree>]() {
-                    let e = lagrange::create::<f64>(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard);
+                    let e = lagrange::create::<f64, f64>(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard);
                     let c = reference_cell::connectivity(ReferenceCellType::[<$cell>]);
 
                     for (dim, entities) in c.iter().enumerate() {
@@ -1840,7 +1891,7 @@ mod test {
 
                 #[test]
                 fn [<test_dof_transformations_ $cell:lower _ $element:lower _ $degree>]() {
-                    let e = [<$element>]::create::<f64>(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard);
+                    let e = [<$element>]::create::<f64, f64>(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard);
                     let tdim = reference_cell::dim(ReferenceCellType::[<$cell>]);
                     for edim in 1..tdim {
                         for entity in &reference_cell::entity_types(ReferenceCellType::[<$cell>])[edim] {
@@ -1944,7 +1995,7 @@ mod test {
 
     #[test]
     fn test_nedelec1_triangle_dof_transformations() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
         if let DofTransformation::Transformation(m, p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -1960,7 +2011,7 @@ mod test {
 
     #[test]
     fn test_nedelec2_triangle_dof_transformations() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
         if let DofTransformation::Transformation(m, p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -1976,7 +2027,8 @@ mod test {
 
     #[test]
     fn test_nedelec1_tetrahedron_dof_transformations() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
         if let DofTransformation::Transformation(m, p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -2006,7 +2058,8 @@ mod test {
 
     #[test]
     fn test_nedelec2_tetrahedron_dof_transformations() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Tetrahedron, 2, Continuity::Standard);
+        let e =
+            nedelec::create::<f64, f64>(ReferenceCellType::Tetrahedron, 2, Continuity::Standard);
         if let DofTransformation::Transformation(m, p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -2045,7 +2098,7 @@ mod test {
 
     #[test]
     fn test_nedelec2_quadrilateral_dof_transformations() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
         if let DofTransformation::Transformation(m, p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -2061,7 +2114,7 @@ mod test {
 
     #[test]
     fn test_nedelec2_hexahedron_dof_transformations() {
-        let e = nedelec::create::<f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
+        let e = nedelec::create::<f64, f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
         if let DofTransformation::Transformation(m, p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -2111,7 +2164,8 @@ mod test {
 
     #[test]
     fn test_lagrange4_tetrahedron_dof_transformations() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Tetrahedron, 4, Continuity::Standard);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Tetrahedron, 4, Continuity::Standard);
         if let DofTransformation::Permutation(p) = e
             .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection)
             .unwrap()
@@ -2155,7 +2209,7 @@ mod test {
 
     #[test]
     fn test_dof_permuting_triangle() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
+        let e = lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
 
         let mut n = (0..10).collect::<Vec<_>>();
         e.apply_dof_permutations(&mut n, 7);
@@ -2177,7 +2231,8 @@ mod test {
 
     #[test]
     fn test_dof_permuting_tetrahedron() {
-        let e = lagrange::create::<f64>(ReferenceCellType::Tetrahedron, 3, Continuity::Standard);
+        let e =
+            lagrange::create::<f64, f64>(ReferenceCellType::Tetrahedron, 3, Continuity::Standard);
 
         let mut n = (0..20).collect::<Vec<_>>();
         e.apply_dof_permutations(&mut n, 63);
diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs
index af8dfcc..dc9eed1 100644
--- a/src/ciarlet/lagrange.rs
+++ b/src/ciarlet/lagrange.rs
@@ -11,11 +11,11 @@ use rlst::{RlstScalar, rlst_dynamic_array};
 use std::marker::PhantomData;
 
 /// Create a Lagrange element.
-pub fn create<T: RlstScalar + Getrf + Getri>(
+pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, IdentityMap> {
+) -> CiarletElement<T, IdentityMap, TGeo> {
     let dim = polynomial_count(cell_type, degree);
     let tdim = reference_cell::dim(cell_type);
     let mut wcoeffs = rlst_dynamic_array!(T, [dim, 1, dim]);
@@ -26,24 +26,24 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
     let mut x = [vec![], vec![], vec![], vec![]];
     let mut m = [vec![], vec![], vec![], vec![]];
     let entity_counts = reference_cell::entity_counts(cell_type);
-    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+    let vertices = reference_cell::vertices::<TGeo>(cell_type);
     if degree == 0 {
         if continuity == Continuity::Standard {
             panic!("Cannot create continuous degree 0 Lagrange element");
         }
         for (d, counts) in entity_counts.iter().enumerate() {
             for _e in 0..*counts {
-                x[d].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+                x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
                 m[d].push(rlst_dynamic_array!(T, [0, 1, 0]));
             }
         }
-        let mut midp = rlst_dynamic_array!(T::Real, [tdim, 1]);
+        let mut midp = rlst_dynamic_array!(TGeo, [tdim, 1]);
         let nvertices = entity_counts[0];
         for i in 0..tdim {
             for vertex in &vertices {
-                *midp.get_mut([i, 0]).unwrap() += num::cast::<_, T::Real>(vertex[i]).unwrap();
+                *midp.get_mut([i, 0]).unwrap() += num::cast::<_, TGeo>(vertex[i]).unwrap();
             }
-            *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, T::Real>(nvertices).unwrap();
+            *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, TGeo>(nvertices).unwrap();
         }
         x[tdim].push(midp);
         let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]);
@@ -54,9 +54,9 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
         let faces = reference_cell::faces(cell_type);
         let volumes = reference_cell::volumes(cell_type);
         for vertex in &vertices {
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, 1]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, 1]);
             for (i, v) in vertex.iter().enumerate() {
-                *pts.get_mut([i, 0]).unwrap() = num::cast::<_, T::Real>(*v).unwrap();
+                *pts.get_mut([i, 0]).unwrap() = num::cast::<_, TGeo>(*v).unwrap();
             }
             x[0].push(pts);
             let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]);
@@ -64,7 +64,7 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
             m[0].push(mentry);
         }
         for e in &edges {
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, degree - 1]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, degree - 1]);
             let [vn0, vn1] = e[..] else {
                 panic!();
             };
@@ -75,10 +75,9 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
             for i in 1..degree {
                 *ident.get_mut([i - 1, 0, i - 1]).unwrap() = T::from(1.0).unwrap();
                 for j in 0..tdim {
-                    *pts.get_mut([j, i - 1]).unwrap() = num::cast::<_, T::Real>(v0[j]).unwrap()
-                        + num::cast::<_, T::Real>(i).unwrap()
-                            / num::cast::<_, T::Real>(degree).unwrap()
-                            * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap();
+                    *pts.get_mut([j, i - 1]).unwrap() = num::cast::<_, TGeo>(v0[j]).unwrap()
+                        + num::cast::<_, TGeo>(i).unwrap() / num::cast::<_, TGeo>(degree).unwrap()
+                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap();
                 }
             }
             x[1].push(pts);
@@ -101,7 +100,7 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
                     panic!("Unsupported face type");
                 }
             };
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, npts]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
 
             let [vn0, vn1, vn2] = faces[e][..3] else {
                 panic!();
@@ -116,14 +115,14 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
                     for i0 in 1..degree {
                         for i1 in 1..degree - i0 {
                             for j in 0..tdim {
-                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j])
+                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
                                     .unwrap()
-                                    + num::cast::<_, T::Real>(i0).unwrap()
-                                        / num::cast::<_, T::Real>(degree).unwrap()
-                                        * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap()
-                                    + num::cast::<_, T::Real>(i1).unwrap()
-                                        / num::cast::<_, T::Real>(degree).unwrap()
-                                        * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap();
+                                    + num::cast::<_, TGeo>(i0).unwrap()
+                                        / num::cast::<_, TGeo>(degree).unwrap()
+                                        * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
+                                    + num::cast::<_, TGeo>(i1).unwrap()
+                                        / num::cast::<_, TGeo>(degree).unwrap()
+                                        * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap();
                             }
                             n += 1;
                         }
@@ -134,14 +133,14 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
                     for i0 in 1..degree {
                         for i1 in 1..degree {
                             for j in 0..tdim {
-                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j])
+                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
                                     .unwrap()
-                                    + num::cast::<_, T::Real>(i0).unwrap()
-                                        / num::cast::<_, T::Real>(degree).unwrap()
-                                        * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap()
-                                    + num::cast::<_, T::Real>(i1).unwrap()
-                                        / num::cast::<_, T::Real>(degree).unwrap()
-                                        * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap();
+                                    + num::cast::<_, TGeo>(i0).unwrap()
+                                        / num::cast::<_, TGeo>(degree).unwrap()
+                                        * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
+                                    + num::cast::<_, TGeo>(i1).unwrap()
+                                        / num::cast::<_, TGeo>(degree).unwrap()
+                                        * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap();
                             }
                             n += 1;
                         }
@@ -176,7 +175,7 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
                     panic!("Unsupported face type");
                 }
             };
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, npts]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
 
             match volume_type {
                 ReferenceCellType::Tetrahedron => {
@@ -193,17 +192,17 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
                         for i1 in 1..degree - i0 {
                             for i2 in 1..degree - i0 - i1 {
                                 for j in 0..tdim {
-                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j])
+                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
                                         .unwrap()
-                                        + num::cast::<_, T::Real>(i0).unwrap()
-                                            / num::cast::<_, T::Real>(degree).unwrap()
-                                            * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap()
-                                        + num::cast::<_, T::Real>(i1).unwrap()
-                                            / num::cast::<_, T::Real>(degree).unwrap()
-                                            * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap()
-                                        + num::cast::<_, T::Real>(i2).unwrap()
-                                            / num::cast::<_, T::Real>(degree).unwrap()
-                                            * num::cast::<_, T::Real>(v3[j] - v0[j]).unwrap();
+                                        + num::cast::<_, TGeo>(i0).unwrap()
+                                            / num::cast::<_, TGeo>(degree).unwrap()
+                                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
+                                        + num::cast::<_, TGeo>(i1).unwrap()
+                                            / num::cast::<_, TGeo>(degree).unwrap()
+                                            * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap()
+                                        + num::cast::<_, TGeo>(i2).unwrap()
+                                            / num::cast::<_, TGeo>(degree).unwrap()
+                                            * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap();
                                 }
                                 n += 1;
                             }
@@ -224,17 +223,17 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
                         for i1 in 1..degree {
                             for i2 in 1..degree {
                                 for j in 0..tdim {
-                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j])
+                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
                                         .unwrap()
-                                        + num::cast::<_, T::Real>(i0).unwrap()
-                                            / num::cast::<_, T::Real>(degree).unwrap()
-                                            * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap()
-                                        + num::cast::<_, T::Real>(i1).unwrap()
-                                            / num::cast::<_, T::Real>(degree).unwrap()
-                                            * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap()
-                                        + num::cast::<_, T::Real>(i2).unwrap()
-                                            / num::cast::<_, T::Real>(degree).unwrap()
-                                            * num::cast::<_, T::Real>(v3[j] - v0[j]).unwrap();
+                                        + num::cast::<_, TGeo>(i0).unwrap()
+                                            / num::cast::<_, TGeo>(degree).unwrap()
+                                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
+                                        + num::cast::<_, TGeo>(i1).unwrap()
+                                            / num::cast::<_, TGeo>(degree).unwrap()
+                                            * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap()
+                                        + num::cast::<_, TGeo>(i2).unwrap()
+                                            / num::cast::<_, TGeo>(degree).unwrap()
+                                            * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap();
                                 }
                                 n += 1;
                             }
@@ -254,7 +253,7 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
             m[3].push(ident);
         }
     }
-    CiarletElement::<T, IdentityMap>::create(
+    CiarletElement::<T, IdentityMap, TGeo>::create(
         "Lagrange".to_string(),
         cell_type,
         degree,
@@ -272,28 +271,32 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
 ///
 /// A family of Lagrange elements on multiple cell types with appropriate
 /// continuity across different cell types.
-pub struct LagrangeElementFamily<T: RlstScalar + Getrf + Getri> {
+pub struct LagrangeElementFamily<T: RlstScalar + Getrf + Getri = f64, TGeo: RlstScalar = f64> {
     degree: usize,
     continuity: Continuity,
     _t: PhantomData<T>,
+    _tgeo: PhantomData<TGeo>,
 }
 
-impl<T: RlstScalar + Getrf + Getri> LagrangeElementFamily<T> {
+impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> LagrangeElementFamily<T, TGeo> {
     /// Create new family with given `degree` and `continuity`.
     pub fn new(degree: usize, continuity: Continuity) -> Self {
         Self {
             degree,
             continuity,
             _t: PhantomData,
+            _tgeo: PhantomData,
         }
     }
 }
 
-impl<T: RlstScalar + Getrf + Getri> ElementFamily for LagrangeElementFamily<T> {
+impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily
+    for LagrangeElementFamily<T, TGeo>
+{
     type T = T;
-    type FiniteElement = CiarletElement<T, IdentityMap>;
+    type FiniteElement = CiarletElement<T, IdentityMap, TGeo>;
     type CellType = ReferenceCellType;
-    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, IdentityMap> {
-        create::<T>(cell_type, self.degree, self.continuity)
+    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, IdentityMap, TGeo> {
+        create::<T, TGeo>(cell_type, self.degree, self.continuity)
     }
 }
diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs
index 23bd707..2c31a54 100644
--- a/src/ciarlet/nedelec.rs
+++ b/src/ciarlet/nedelec.rs
@@ -13,11 +13,11 @@ use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri};
 use rlst::{DynArray, RlstScalar, SliceArray, rlst_dynamic_array};
 use std::marker::PhantomData;
 
-fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getrf + Getri>(
+fn create_simplex<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, CovariantPiolaMap> {
+) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
     if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Tetrahedron {
         panic!("Invalid cell: {cell_type:?}");
     }
@@ -31,9 +31,9 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     let pts_t = cell_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let pts = SliceArray::<TReal, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
+    let pts = SliceArray::<TGeo, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
 
     let mut phi = DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree, 0));
     tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
@@ -178,10 +178,10 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     let mut m = [vec![], vec![], vec![], vec![]];
 
     let entity_counts = reference_cell::entity_counts(cell_type);
-    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+    let vertices = reference_cell::vertices::<TGeo>(cell_type);
 
     for _ in 0..entity_counts[0] {
-        x[0].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+        x[0].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
         m[0].push(rlst_dynamic_array!(T, [0, tdim, 0]));
     }
 
@@ -190,9 +190,9 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     let edge_pts_t = edge_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let edge_pts = SliceArray::<TReal, 2>::from_shape(&edge_pts_t, [1, edge_q.npoints]);
+    let edge_pts = SliceArray::<TGeo, 2>::from_shape(&edge_pts_t, [1, edge_q.npoints]);
 
     let mut edge_phi = DynArray::<T, 3>::from_shape(legendre_shape(
         ReferenceCellType::Interval,
@@ -209,7 +209,7 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     );
 
     for edge in reference_cell::edges(cell_type) {
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, edge_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, edge_q.npoints]);
         let mut mat = rlst_dynamic_array!(T, [pdim_edge_minus1, tdim, edge_q.npoints]);
 
         for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() {
@@ -232,7 +232,7 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     // DOFs on faces
     if degree == 1 {
         for _ in 0..entity_counts[2] {
-            x[2].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[2].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[2].push(rlst_dynamic_array!(T, [0, tdim, 0]))
         }
     } else {
@@ -240,9 +240,9 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
         let face_pts_t = face_q
             .points
             .iter()
-            .map(|i| TReal::from(*i).unwrap())
+            .map(|i| TGeo::from(*i).unwrap())
             .collect::<Vec<_>>();
-        let face_pts = SliceArray::<TReal, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
+        let face_pts = SliceArray::<TGeo, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
 
         let mut face_phi = DynArray::<T, 3>::from_shape(legendre_shape(
             ReferenceCellType::Triangle,
@@ -259,7 +259,7 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
         );
 
         for face in reference_cell::faces(cell_type) {
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, face_q.npoints]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]);
             let mut mat = rlst_dynamic_array!(T, [2 * pdim_face_minus2, tdim, face_q.npoints]);
 
             for (w_i, wt) in face_q.weights.iter().enumerate() {
@@ -284,14 +284,14 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     }
     if tdim == 3 {
         if degree <= 2 {
-            x[3].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[3].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[3].push(rlst_dynamic_array!(T, [0, tdim, 0]))
         } else {
             let internal_q = gauss_jacobi_rule(cell_type, 2 * degree - 3).unwrap();
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, internal_q.npoints]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, internal_q.npoints]);
             for p in 0..internal_q.npoints {
                 for d in 0..3 {
-                    pts[[d, p]] = TReal::from(internal_q.points[3 * p + d]).unwrap()
+                    pts[[d, p]] = TGeo::from(internal_q.points[3 * p + d]).unwrap()
                 }
             }
 
@@ -328,11 +328,11 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     )
 }
 
-fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getrf + Getri>(
+fn create_tp<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, CovariantPiolaMap> {
+) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
     if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron {
         panic!("Invalid cell: {cell_type:?}");
     }
@@ -346,9 +346,9 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     let pts_t = cell_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let pts = SliceArray::<TReal, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
+    let pts = SliceArray::<TGeo, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
 
     let mut phi = DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree, 0));
     tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
@@ -425,10 +425,10 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     let mut x = [vec![], vec![], vec![], vec![]];
     let mut m = [vec![], vec![], vec![], vec![]];
 
-    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+    let vertices = reference_cell::vertices::<TGeo>(cell_type);
 
     for _ in 0..entity_counts[0] {
-        x[0].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+        x[0].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
         m[0].push(rlst_dynamic_array!(T, [0, tdim, 0]));
     }
 
@@ -437,9 +437,9 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     let edge_pts_t = edge_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let edge_pts = SliceArray::<TReal, 2>::from_shape(&edge_pts_t, [1, edge_q.npoints]);
+    let edge_pts = SliceArray::<TGeo, 2>::from_shape(&edge_pts_t, [1, edge_q.npoints]);
 
     let mut edge_phi = DynArray::<T, 3>::from_shape(legendre_shape(
         ReferenceCellType::Interval,
@@ -456,7 +456,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     );
 
     for edge in reference_cell::edges(cell_type) {
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, edge_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, edge_q.npoints]);
         let mut mat = rlst_dynamic_array!(T, [pdim_edge_minus1, tdim, edge_q.npoints]);
 
         for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() {
@@ -479,7 +479,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     // DOFs on faces
     if degree == 1 {
         for _ in 0..entity_counts[2] {
-            x[2].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[2].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[2].push(rlst_dynamic_array!(T, [0, tdim, 0]))
         }
     } else {
@@ -487,9 +487,9 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
         let face_pts_t = face_q
             .points
             .iter()
-            .map(|i| TReal::from(*i).unwrap())
+            .map(|i| TGeo::from(*i).unwrap())
             .collect::<Vec<_>>();
-        let face_pts = SliceArray::<TReal, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
+        let face_pts = SliceArray::<TGeo, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
 
         let mut face_phi = DynArray::<T, 3>::from_shape(legendre_shape(
             ReferenceCellType::Quadrilateral,
@@ -506,7 +506,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
         );
 
         for face in reference_cell::faces(cell_type) {
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, face_q.npoints]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]);
             let mdim = 2 * pdim_edge_minus2 * pdim_edge_minus1;
             let mut mat = rlst_dynamic_array!(T, [mdim, tdim, face_q.npoints]);
 
@@ -539,7 +539,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     // DOFs on volume
     if tdim == 3 {
         if degree == 1 {
-            x[3].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[3].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[3].push(rlst_dynamic_array!(T, [0, tdim, 0]))
         } else {
             let interior_q =
@@ -547,10 +547,10 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
             let interior_pts_t = interior_q
                 .points
                 .iter()
-                .map(|i| TReal::from(*i).unwrap())
+                .map(|i| TGeo::from(*i).unwrap())
                 .collect::<Vec<_>>();
             let interior_pts =
-                SliceArray::<TReal, 2>::from_shape(&interior_pts_t, [3, interior_q.npoints]);
+                SliceArray::<TGeo, 2>::from_shape(&interior_pts_t, [3, interior_q.npoints]);
 
             let mut interior_phi = DynArray::<T, 3>::from_shape(legendre_shape(
                 ReferenceCellType::Hexahedron,
@@ -566,7 +566,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
                 &mut interior_phi,
             );
 
-            let mut pts = rlst_dynamic_array!(T::Real, [tdim, interior_q.npoints]);
+            let mut pts = rlst_dynamic_array!(TGeo, [tdim, interior_q.npoints]);
             let mdim = 3 * pdim_edge_minus2.pow(2) * pdim_edge_minus1;
             let mut mat = rlst_dynamic_array!(T, [mdim, tdim, interior_q.npoints]);
 
@@ -623,17 +623,17 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
 }
 
 /// Create a Nedelec (first kind) element
-pub fn create<T: RlstScalar + Getrf + Getri>(
+pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, CovariantPiolaMap> {
+) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
     if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron {
-        create_simplex(cell_type, degree, continuity)
+        create_simplex::<T, TGeo>(cell_type, degree, continuity)
     } else if cell_type == ReferenceCellType::Quadrilateral
         || cell_type == ReferenceCellType::Hexahedron
     {
-        create_tp(cell_type, degree, continuity)
+        create_tp::<T, TGeo>(cell_type, degree, continuity)
     } else {
         panic!("Invalid cell: {cell_type:?}");
     }
@@ -643,28 +643,35 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
 ///
 /// A family of Nedelec elements on multiple cell types with appropriate
 /// continuity across different cell types.
-pub struct NedelecFirstKindElementFamily<T: RlstScalar + Getrf + Getri> {
+pub struct NedelecFirstKindElementFamily<
+    T: RlstScalar + Getrf + Getri = f64,
+    TGeo: RlstScalar = f64,
+> {
     degree: usize,
     continuity: Continuity,
     _t: PhantomData<T>,
+    _tgeo: PhantomData<TGeo>,
 }
 
-impl<T: RlstScalar + Getrf + Getri> NedelecFirstKindElementFamily<T> {
+impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> NedelecFirstKindElementFamily<T, TGeo> {
     /// Create new family with given `degree` and `continuity`.
     pub fn new(degree: usize, continuity: Continuity) -> Self {
         Self {
             degree,
             continuity,
             _t: PhantomData,
+            _tgeo: PhantomData,
         }
     }
 }
 
-impl<T: RlstScalar + Getrf + Getri> ElementFamily for NedelecFirstKindElementFamily<T> {
+impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily
+    for NedelecFirstKindElementFamily<T, TGeo>
+{
     type T = T;
     type CellType = ReferenceCellType;
-    type FiniteElement = CiarletElement<T, CovariantPiolaMap>;
-    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, CovariantPiolaMap> {
-        create::<T>(cell_type, self.degree, self.continuity)
+    type FiniteElement = CiarletElement<T, CovariantPiolaMap, TGeo>;
+    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
+        create::<T, TGeo>(cell_type, self.degree, self.continuity)
     }
 }
diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs
index a56b820..b31d8ad 100644
--- a/src/ciarlet/raviart_thomas.rs
+++ b/src/ciarlet/raviart_thomas.rs
@@ -13,11 +13,11 @@ use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri};
 use rlst::{DynArray, RlstScalar, SliceArray, rlst_dynamic_array};
 use std::marker::PhantomData;
 
-fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getrf + Getri>(
+fn create_simplex<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, ContravariantPiolaMap> {
+) -> CiarletElement<T, ContravariantPiolaMap, TGeo> {
     if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Tetrahedron {
         panic!("Invalid cell: {cell_type:?}");
     }
@@ -44,9 +44,9 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     let pts_t = cell_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let pts = SliceArray::<TReal, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
+    let pts = SliceArray::<TGeo, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
 
     let mut phi = DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree, 0));
     tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
@@ -87,11 +87,11 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     let mut m = [vec![], vec![], vec![], vec![]];
 
     let entity_counts = reference_cell::entity_counts(cell_type);
-    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+    let vertices = reference_cell::vertices::<TGeo>(cell_type);
 
     for d in 0..tdim - 1 {
         for _ in 0..entity_counts[d] {
-            x[d].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[d].push(rlst_dynamic_array!(T, [0, tdim, 0]));
         }
     }
@@ -101,9 +101,9 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     let facet_pts_t = facet_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let facet_pts = SliceArray::<TReal, 2>::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]);
+    let facet_pts = SliceArray::<TGeo, 2>::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]);
 
     let mut facet_phi =
         DynArray::<T, 3>::from_shape(legendre_shape(facet_type, &facet_pts, degree - 1, 0));
@@ -115,9 +115,9 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
         } else {
             reference_cell::faces(cell_type)
         },
-        reference_cell::facet_normals::<TReal>(cell_type),
+        reference_cell::facet_normals::<TGeo>(cell_type),
     ) {
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, facet_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, facet_q.npoints]);
         let mut mat = rlst_dynamic_array!(T, [pdim_facet_minus1, tdim, facet_q.npoints]);
 
         for (w_i, wt) in facet_q.weights.iter().enumerate() {
@@ -128,7 +128,7 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
                         &facet_q.points[w_i * (tdim - 1)..(w_i + 1) * (tdim - 1)]
                     )
                     .map(|(v_i, pt)| {
-                        (vertices[*v_i][i] - vertices[facet[0]][i]) * TReal::from(*pt).unwrap()
+                        (vertices[*v_i][i] - vertices[facet[0]][i]) * TGeo::from(*pt).unwrap()
                     })
                     .sum();
 
@@ -146,15 +146,15 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
 
     if degree == 1 {
         for _ in 0..entity_counts[tdim] {
-            x[tdim].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[tdim].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[tdim].push(rlst_dynamic_array!(T, [0, tdim, 0]))
         }
     } else {
         let internal_q = gauss_jacobi_rule(cell_type, 2 * degree - 2).unwrap();
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, internal_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, internal_q.npoints]);
         for p in 0..internal_q.npoints {
             for d in 0..tdim {
-                pts[[d, p]] = TReal::from(internal_q.points[tdim * p + d]).unwrap()
+                pts[[d, p]] = TGeo::from(internal_q.points[tdim * p + d]).unwrap()
             }
         }
 
@@ -190,11 +190,11 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     )
 }
 
-fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getrf + Getri>(
+fn create_tp<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, ContravariantPiolaMap> {
+) -> CiarletElement<T, ContravariantPiolaMap, TGeo> {
     if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron {
         panic!("Invalid cell: {cell_type:?}");
     }
@@ -208,9 +208,9 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     let pts_t = cell_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let pts = SliceArray::<TReal, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
+    let pts = SliceArray::<TGeo, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
 
     let mut phi = DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree, 0));
     tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
@@ -293,11 +293,11 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     let mut x = [vec![], vec![], vec![], vec![]];
     let mut m = [vec![], vec![], vec![], vec![]];
 
-    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+    let vertices = reference_cell::vertices::<TGeo>(cell_type);
 
     for d in 0..tdim - 1 {
         for _ in 0..entity_counts[d] {
-            x[d].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+            x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
             m[d].push(rlst_dynamic_array!(T, [0, tdim, 0]));
         }
     }
@@ -306,9 +306,9 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
     let facet_pts_t = facet_q
         .points
         .iter()
-        .map(|i| TReal::from(*i).unwrap())
+        .map(|i| TGeo::from(*i).unwrap())
         .collect::<Vec<_>>();
-    let facet_pts = SliceArray::<TReal, 2>::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]);
+    let facet_pts = SliceArray::<TGeo, 2>::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]);
 
     let mut facet_phi =
         DynArray::<T, 3>::from_shape(legendre_shape(facet_type, &facet_pts, degree - 1, 0));
@@ -316,9 +316,9 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
 
     for (facet, normal) in izip!(
         reference_cell::facets(cell_type),
-        reference_cell::facet_normals::<T::Real>(cell_type),
+        reference_cell::facet_normals::<TGeo>(cell_type),
     ) {
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, facet_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, facet_q.npoints]);
         let mut mat = rlst_dynamic_array!(T, [pdim_facet_minus1, tdim, facet_q.npoints]);
 
         for (w_i, wt) in facet_q.weights.iter().enumerate() {
@@ -329,7 +329,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
                             (vertices[facet[usize::pow(2, x as u32)]][d] - vertices[facet[0]][d])
                                 * facet_pts[[x, w_i]]
                         })
-                        .sum::<T::Real>();
+                        .sum::<TGeo>();
                 for j in 0..facet_phi.shape()[1] {
                     mat[[j, d, w_i]] = T::from(*wt).unwrap()
                         * facet_phi[[0, j, w_i]]
@@ -344,16 +344,16 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
 
     // DOFs on interior
     if degree == 1 {
-        x[tdim].push(rlst_dynamic_array!(T::Real, [tdim, 0]));
+        x[tdim].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
         m[tdim].push(rlst_dynamic_array!(T, [0, tdim, 0]))
     } else if tdim == 2 {
         let face_q = gauss_jacobi_rule(ReferenceCellType::Quadrilateral, 2 * degree - 1).unwrap();
         let face_pts_t = face_q
             .points
             .iter()
-            .map(|i| TReal::from(*i).unwrap())
+            .map(|i| TGeo::from(*i).unwrap())
             .collect::<Vec<_>>();
-        let face_pts = SliceArray::<TReal, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
+        let face_pts = SliceArray::<TGeo, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
 
         let mut face_phi = DynArray::<T, 3>::from_shape(legendre_shape(
             ReferenceCellType::Quadrilateral,
@@ -369,7 +369,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
             &mut face_phi,
         );
 
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, face_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]);
         let mdim = 2 * pdim_edge_minus2 * pdim_edge_minus1;
         let mut mat = rlst_dynamic_array!(T, [mdim, tdim, face_q.npoints]);
 
@@ -394,10 +394,10 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
         let interior_pts_t = interior_q
             .points
             .iter()
-            .map(|i| TReal::from(*i).unwrap())
+            .map(|i| TGeo::from(*i).unwrap())
             .collect::<Vec<_>>();
         let interior_pts =
-            SliceArray::<TReal, 2>::from_shape(&interior_pts_t, [3, interior_q.npoints]);
+            SliceArray::<TGeo, 2>::from_shape(&interior_pts_t, [3, interior_q.npoints]);
 
         let mut interior_phi = DynArray::<T, 3>::from_shape(legendre_shape(
             ReferenceCellType::Hexahedron,
@@ -413,7 +413,7 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
             &mut interior_phi,
         );
 
-        let mut pts = rlst_dynamic_array!(T::Real, [tdim, interior_q.npoints]);
+        let mut pts = rlst_dynamic_array!(TGeo, [tdim, interior_q.npoints]);
         let mdim = 3 * pdim_edge_minus1.pow(2) * pdim_edge_minus2;
         let mut mat = rlst_dynamic_array!(T, [mdim, tdim, interior_q.npoints]);
 
@@ -466,11 +466,11 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Getr
 }
 
 /// Create a Raviart-Thomas element
-pub fn create<T: RlstScalar + Getrf + Getri>(
+pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
     cell_type: ReferenceCellType,
     degree: usize,
     continuity: Continuity,
-) -> CiarletElement<T, ContravariantPiolaMap> {
+) -> CiarletElement<T, ContravariantPiolaMap, TGeo> {
     if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron {
         create_simplex(cell_type, degree, continuity)
     } else if cell_type == ReferenceCellType::Quadrilateral
@@ -486,28 +486,35 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
 ///
 /// A family of Raviart-Thomas elements on multiple cell types with appropriate
 /// continuity across different cell types.
-pub struct RaviartThomasElementFamily<T: RlstScalar + Getrf + Getri> {
+pub struct RaviartThomasElementFamily<T: RlstScalar + Getrf + Getri = f64, TGeo: RlstScalar = f64> {
     degree: usize,
     continuity: Continuity,
     _t: PhantomData<T>,
+    _tgeo: PhantomData<TGeo>,
 }
 
-impl<T: RlstScalar + Getrf + Getri> RaviartThomasElementFamily<T> {
+impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> RaviartThomasElementFamily<T, TGeo> {
     /// Create new family with given `degree` and `continuity`.
     pub fn new(degree: usize, continuity: Continuity) -> Self {
         Self {
             degree,
             continuity,
             _t: PhantomData,
+            _tgeo: PhantomData,
         }
     }
 }
 
-impl<T: RlstScalar + Getrf + Getri> ElementFamily for RaviartThomasElementFamily<T> {
+impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily
+    for RaviartThomasElementFamily<T, TGeo>
+{
     type T = T;
     type CellType = ReferenceCellType;
-    type FiniteElement = CiarletElement<T, ContravariantPiolaMap>;
-    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, ContravariantPiolaMap> {
-        create::<T>(cell_type, self.degree, self.continuity)
+    type FiniteElement = CiarletElement<T, ContravariantPiolaMap, TGeo>;
+    fn element(
+        &self,
+        cell_type: ReferenceCellType,
+    ) -> CiarletElement<T, ContravariantPiolaMap, TGeo> {
+        create::<T, TGeo>(cell_type, self.degree, self.continuity)
     }
 }
diff --git a/src/map.rs b/src/map.rs
index 12aa282..5ae2ddb 100644
--- a/src/map.rs
+++ b/src/map.rs
@@ -11,16 +11,17 @@ pub struct IdentityMap {}
 impl Map for IdentityMap {
     fn push_forward<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         reference_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        _jacobians: &Array<Array3RealImpl, 3>,
-        _jacobian_determinants: &[T::Real],
-        _inverse_jacobians: &Array<Array3RealImpl, 3>,
+        _jacobians: &Array<Array3GeoImpl, 3>,
+        _jacobian_determinants: &[TGeo],
+        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
         physical_values: &mut Array<Array4MutImpl, 4>,
     ) {
         if nderivs > 0 {
@@ -31,16 +32,17 @@ impl Map for IdentityMap {
     }
     fn pull_back<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         physical_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        _jacobians: &Array<Array3RealImpl, 3>,
-        _jacobian_determinants: &[T::Real],
-        _inverse_jacobians: &Array<Array3RealImpl, 3>,
+        _jacobians: &Array<Array3GeoImpl, 3>,
+        _jacobian_determinants: &[TGeo],
+        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
         reference_values: &mut Array<Array4MutImpl, 4>,
     ) {
         if nderivs > 0 {
@@ -70,16 +72,17 @@ pub struct CovariantPiolaMap {}
 impl Map for CovariantPiolaMap {
     fn push_forward<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         reference_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        _jacobians: &Array<Array3RealImpl, 3>,
-        _jacobian_determinants: &[T::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        _jacobians: &Array<Array3GeoImpl, 3>,
+        _jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         physical_values: &mut Array<Array4MutImpl, 4>,
     ) {
         let tdim = inverse_jacobians.shape()[1];
@@ -109,16 +112,17 @@ impl Map for CovariantPiolaMap {
     }
     fn pull_back<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         physical_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        _jacobian_determinants: &[T::Real],
-        _inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        _jacobian_determinants: &[TGeo],
+        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
         reference_values: &mut Array<Array4MutImpl, 4>,
     ) {
         let gdim = jacobians.shape()[1];
@@ -167,16 +171,17 @@ pub struct ContravariantPiolaMap {}
 impl Map for ContravariantPiolaMap {
     fn push_forward<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         reference_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[T::Real],
-        _inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
         physical_values: &mut Array<Array4MutImpl, 4>,
     ) {
         let gdim = jacobians.shape()[1];
@@ -207,16 +212,17 @@ impl Map for ContravariantPiolaMap {
     }
     fn pull_back<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         physical_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        _jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[T::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        _jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         reference_values: &mut Array<Array4MutImpl, 4>,
     ) {
         let tdim = inverse_jacobians.shape()[1];
@@ -256,7 +262,7 @@ mod test {
     use approx::*;
     use rlst::{Array, rlst_dynamic_array};
 
-    fn fill_jacobians<T: RlstScalar<Real = T>, Array3MutImpl: MutableArrayImpl<T, 3>>(
+    fn fill_jacobians<T: RlstScalar, Array3MutImpl: MutableArrayImpl<T, 3>>(
         j: &mut Array<Array3MutImpl, 3>,
         jdet: &mut [T],
         jinv: &mut Array<Array3MutImpl, 3>,
diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs
index ae4449b..3865d62 100644
--- a/src/polynomials/legendre.rs
+++ b/src/polynomials/legendre.rs
@@ -37,7 +37,8 @@ fn jrc<T: RlstScalar>(a: usize, n: usize) -> (T, T, T) {
 /// Tabulate orthonormal polynomials on an interval
 fn tabulate_interval<
     T: RlstScalar,
-    Array2Impl: ValueArrayImpl<T::Real, 2>,
+    TGeo: RlstScalar,
+    Array2Impl: ValueArrayImpl<TGeo, 2>,
     Array3MutImpl: MutableArrayImpl<T, 3>,
 >(
     points: &Array<Array2Impl, 2>,
@@ -98,7 +99,8 @@ fn tabulate_interval<
 /// Tabulate orthonormal polynomials on a quadrilateral
 fn tabulate_quadrilateral<
     T: RlstScalar,
-    Array2Impl: ValueArrayImpl<T::Real, 2>,
+    TGeo: RlstScalar,
+    Array2Impl: ValueArrayImpl<TGeo, 2>,
     Array3MutImpl: MutableArrayImpl<T, 3>,
 >(
     points: &Array<Array2Impl, 2>,
@@ -251,7 +253,8 @@ fn tabulate_quadrilateral<
 /// Tabulate orthonormal polynomials on a triangle
 fn tabulate_triangle<
     T: RlstScalar,
-    Array2Impl: ValueArrayImpl<T::Real, 2>,
+    TGeo: RlstScalar,
+    Array2Impl: ValueArrayImpl<TGeo, 2>,
     Array3MutImpl: MutableArrayImpl<T, 3>,
 >(
     points: &Array<Array2Impl, 2>,
@@ -452,7 +455,8 @@ fn tabulate_triangle<
 /// Tabulate orthonormal polynomials on a tetrahedron
 fn tabulate_tetrahedron<
     T: RlstScalar,
-    Array2Impl: ValueArrayImpl<T::Real, 2>,
+    TGeo: RlstScalar,
+    Array2Impl: ValueArrayImpl<TGeo, 2>,
     Array3MutImpl: MutableArrayImpl<T, 3>,
 >(
     points: &Array<Array2Impl, 2>,
@@ -815,7 +819,8 @@ fn tabulate_tetrahedron<
 /// Tabulate orthonormal polynomials on a hexahedron
 fn tabulate_hexahedron<
     T: RlstScalar,
-    Array2Impl: ValueArrayImpl<T::Real, 2>,
+    TGeo: RlstScalar,
+    Array2Impl: ValueArrayImpl<TGeo, 2>,
     Array3MutImpl: MutableArrayImpl<T, 3>,
 >(
     points: &Array<Array2Impl, 2>,
@@ -1049,7 +1054,8 @@ pub fn shape<T, Array2Impl: ValueArrayImpl<T, 2>>(
 /// Tabulate orthonormal polynomials
 pub fn tabulate<
     T: RlstScalar,
-    Array2Impl: ValueArrayImpl<T::Real, 2>,
+    TGeo: RlstScalar,
+    Array2Impl: ValueArrayImpl<TGeo, 2>,
     Array3MutImpl: MutableArrayImpl<T, 3>,
 >(
     cell_type: ReferenceCellType,
diff --git a/src/reference_cell.rs b/src/reference_cell.rs
index 905a719..063e3b0 100644
--- a/src/reference_cell.rs
+++ b/src/reference_cell.rs
@@ -40,7 +40,7 @@ pub fn is_simplex(cell: ReferenceCellType) -> bool {
 }
 
 /// The vertices of the reference cell
-pub fn vertices<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<Vec<T>> {
+pub fn vertices<T: RlstScalar>(cell: ReferenceCellType) -> Vec<Vec<T>> {
     let zero = T::from(0.0).unwrap();
     let one = T::from(1.0).unwrap();
     match cell {
@@ -88,7 +88,7 @@ pub fn vertices<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<Vec<T>>
 }
 
 /// The midpoint of the cell
-pub fn midpoint<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<T> {
+pub fn midpoint<T: RlstScalar>(cell: ReferenceCellType) -> Vec<T> {
     let half = T::from(0.5).unwrap();
     let third = T::from(1.0).unwrap() / T::from(3.0).unwrap();
     match cell {
@@ -258,7 +258,7 @@ pub fn peaks(cell: ReferenceCellType) -> Vec<Vec<usize>> {
 }
 
 /// The normals to the facets of the reference cell. The length of each normal is the volume of the facet
-pub fn facet_normals<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<Vec<T>> {
+pub fn facet_normals<T: RlstScalar>(cell: ReferenceCellType) -> Vec<Vec<T>> {
     let zero = T::from(0.0).unwrap();
     let one = T::from(1.0).unwrap();
     match cell {
@@ -303,7 +303,7 @@ pub fn facet_normals<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<Ve
 }
 
 /// The unit normals to the facets of the reference cell
-pub fn facet_unit_normals<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<Vec<T>> {
+pub fn facet_unit_normals<T: RlstScalar>(cell: ReferenceCellType) -> Vec<Vec<T>> {
     let mut normals = facet_normals::<T>(cell);
     for normal in normals.iter_mut() {
         let size = normal.iter().map(|x| x.powi(2)).sum::<T>().sqrt();
diff --git a/src/traits.rs b/src/traits.rs
index 7caae3d..ae21797 100644
--- a/src/traits.rs
+++ b/src/traits.rs
@@ -54,7 +54,8 @@ pub trait FiniteElement {
     ///
     /// For the quadrilaterals and hexahedra, the parameter $n$ denotes the Lagrange superdegree.
     fn tabulate<
-        Array2Impl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 2>,
+        TGeo: RlstScalar,
+        Array2Impl: ValueArrayImpl<TGeo, 2>,
         Array4MutImpl: MutableArrayImpl<Self::T, 4>,
     >(
         &self,
@@ -122,16 +123,17 @@ pub trait MappedFiniteElement: FiniteElement {
     /// - `physical_values`: The output array of the push operation. This shape of this array is the same as the `reference_values`
     ///   input, with the [MappedFiniteElement::physical_value_size] used instead of the reference value size.
     fn push_forward<
-        Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<Self::T, 4>,
         Array4MutImpl: MutableArrayImpl<Self::T, 4>,
     >(
         &self,
         reference_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[<Self::T as RlstScalar>::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         physical_values: &mut Array<Array4MutImpl, 4>,
     );
 
@@ -139,16 +141,17 @@ pub trait MappedFiniteElement: FiniteElement {
     ///
     /// This is the inverse operation to [MappedFiniteElement::push_forward].
     fn pull_back<
-        Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<Self::T, 4>,
         Array4MutImpl: MutableArrayImpl<Self::T, 4>,
     >(
         &self,
         physical_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[<Self::T as RlstScalar>::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         reference_values: &mut Array<Array4MutImpl, 4>,
     );
 
@@ -207,32 +210,34 @@ pub trait Map {
     /// Push values from the reference cell to physical cells.
     fn push_forward<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         reference_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[T::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         physical_values: &mut Array<Array4MutImpl, 4>,
     );
 
     /// Pull values back to the reference cell.
     fn pull_back<
         T: RlstScalar,
-        Array3RealImpl: ValueArrayImpl<T::Real, 3>,
+        TGeo: RlstScalar,
+        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
         Array4Impl: ValueArrayImpl<T, 4>,
         Array4MutImpl: MutableArrayImpl<T, 4>,
     >(
         &self,
         physical_values: &Array<Array4Impl, 4>,
         nderivs: usize,
-        jacobians: &Array<Array3RealImpl, 3>,
-        jacobian_determinants: &[T::Real],
-        inverse_jacobians: &Array<Array3RealImpl, 3>,
+        jacobians: &Array<Array3GeoImpl, 3>,
+        jacobian_determinants: &[TGeo],
+        inverse_jacobians: &Array<Array3GeoImpl, 3>,
         reference_values: &mut Array<Array4MutImpl, 4>,
     );