Remove assumption that geometry uses real scalars (#141)

* remove assumption that everything is real

* use branch of ndelement

* main branch

Author: mscroggs

src/geometry/geometry_map.rs

@@ -1,5 +1,5 @@
 //! Geometry map
-use crate::{traits::GeometryMap as GeometryMapTrait, types::RealScalar};
+use crate::{traits::GeometryMap as GeometryMapTrait, types::Scalar};
 use ndelement::{reference_cell, traits::MappedFiniteElement, types::ReferenceCellType};
 use rlst::{Array, DynArray, RlstScalar, ValueArrayImpl};
 
@@ -43,7 +43,7 @@ fn cross<T: RlstScalar>(mat: &[T], result: &mut [T]) {
     }
 }
 
-impl<'a, T: RealScalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>>
+impl<'a, T: Scalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>>
     GeometryMap<'a, T, B2D, C2D>
 {
     /// Create new
@@ -71,7 +71,7 @@ impl<'a, T: RealScalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>
     }
 }
 
-impl<T: RealScalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>> GeometryMapTrait
+impl<T: Scalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>> GeometryMapTrait
     for GeometryMap<'_, T, B2D, C2D>
 {
     type T = T;

src/geometry/mixed/entity_geometry.rs

@@ -3,21 +3,21 @@ use super::MixedGeometry;
 use crate::{
     geometry::{Point, PointIter},
     traits::Geometry,
-    types::RealScalar,
+    types::Scalar,
 };
 use ndelement::traits::MappedFiniteElement;
 
 /// Geometry of an entity
 #[derive(Debug)]
-pub struct MixedEntityGeometry<'a, T: RealScalar, E: MappedFiniteElement> {
+pub struct MixedEntityGeometry<'a, T: Scalar, E: MappedFiniteElement> {
     geometry: &'a MixedGeometry<T, E>,
     element_index: usize,
     cell_index: usize,
     sub_entity_dim: usize,
     sub_entity_index: usize,
 }
 
-impl<'a, T: RealScalar, E: MappedFiniteElement> MixedEntityGeometry<'a, T, E> {
+impl<'a, T: Scalar, E: MappedFiniteElement> MixedEntityGeometry<'a, T, E> {
     /// Create new
     pub fn new(
         geometry: &'a MixedGeometry<T, E>,
@@ -36,7 +36,7 @@ impl<'a, T: RealScalar, E: MappedFiniteElement> MixedEntityGeometry<'a, T, E> {
     }
 }
 
-impl<T: RealScalar, E: MappedFiniteElement> Geometry for MixedEntityGeometry<'_, T, E> {
+impl<T: Scalar, E: MappedFiniteElement> Geometry for MixedEntityGeometry<'_, T, E> {
     type T = T;
     type Point<'a>
         = Point<'a, T>

src/geometry/mixed/geometry.rs

@@ -1,7 +1,7 @@
 //! Geometry where each entity of a given dimension is represented by the same element
 #[cfg(feature = "serde")]
 use crate::traits::ConvertToSerializable;
-use crate::types::RealScalar;
+use crate::types::Scalar;
 use itertools::izip;
 #[cfg(feature = "serde")]
 use ndelement::{
@@ -19,15 +19,15 @@ use std::collections::HashMap;
 use std::fmt::{Debug, Formatter};
 
 /// Single element geometry
-pub struct MixedGeometry<T: RealScalar, E: MappedFiniteElement> {
+pub struct MixedGeometry<T: Scalar, E: MappedFiniteElement> {
     points: DynArray<T, 2>,
     cells: Vec<DynArray<usize, 2>>,
     elements: Vec<HashMap<ReferenceCellType, E>>,
     pub(crate) insertion_indices_to_element_indices: Vec<usize>,
     pub(crate) insertion_indices_to_cell_indices: Vec<usize>,
 }
 
-impl<T: RealScalar, E: MappedFiniteElement> Debug for MixedGeometry<T, E> {
+impl<T: Scalar, E: MappedFiniteElement> Debug for MixedGeometry<T, E> {
     fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), std::fmt::Error> {
         f.debug_struct("MixedGeometry")
             .field("points", &self.points)
@@ -39,7 +39,7 @@ impl<T: RealScalar, E: MappedFiniteElement> Debug for MixedGeometry<T, E> {
 #[cfg(feature = "serde")]
 #[derive(serde::Serialize, Debug, serde::Deserialize)]
 #[serde(bound = "for<'de2> T: serde::Deserialize<'de2>")]
-pub struct SerializableGeometry<T: RealScalar + serde::Serialize>
+pub struct SerializableGeometry<T: Scalar + serde::Serialize>
 where
     for<'de2> T: serde::Deserialize<'de2>,
 {
@@ -51,8 +51,8 @@ where
 }
 
 #[cfg(feature = "serde")]
-impl<T: RealScalar + serde::Serialize> ConvertToSerializable
-    for MixedGeometry<T, CiarletElement<T, IdentityMap>>
+impl<T: Scalar + serde::Serialize> ConvertToSerializable
+    for MixedGeometry<T, CiarletElement<T, IdentityMap, T>>
 where
     for<'de2> T: serde::Deserialize<'de2>,
 {
@@ -115,7 +115,7 @@ where
     }
 }
 
-impl<T: RealScalar, E: MappedFiniteElement> MixedGeometry<T, E> {
+impl<T: Scalar, E: MappedFiniteElement> MixedGeometry<T, E> {
     /// Create single element geometry
     pub fn new(
         cell_types_in: &[ReferenceCellType],

src/geometry/point.rs

@@ -1,20 +1,20 @@
 //! Point
-use crate::{traits::Point as PointTrait, types::RealScalar};
+use crate::{traits::Point as PointTrait, types::Scalar};
 
 /// A points
 #[derive(Debug, Clone, Copy)]
-pub struct Point<'a, T: RealScalar> {
+pub struct Point<'a, T: Scalar> {
     index: usize,
     coordinates: &'a [T],
 }
 
-impl<'a, T: RealScalar> Point<'a, T> {
+impl<'a, T: Scalar> Point<'a, T> {
     /// Create new
     pub fn new(index: usize, coordinates: &'a [T]) -> Self {
         Self { index, coordinates }
     }
 }
-impl<T: RealScalar> PointTrait for Point<'_, T> {
+impl<T: Scalar> PointTrait for Point<'_, T> {
     type T = T;
 
     fn index(&self) -> usize {
@@ -32,17 +32,17 @@ impl<T: RealScalar> PointTrait for Point<'_, T> {
 
 /// Iterator over points
 #[derive(Debug)]
-pub struct PointIter<'a, T: RealScalar> {
+pub struct PointIter<'a, T: Scalar> {
     points: Vec<(usize, &'a [T])>,
     index: usize,
 }
-impl<'a, T: RealScalar> PointIter<'a, T> {
+impl<'a, T: Scalar> PointIter<'a, T> {
     /// Create new
     pub fn new(points: Vec<(usize, &'a [T])>) -> Self {
         Self { points, index: 0 }
     }
 }
-impl<'a, T: RealScalar> Iterator for PointIter<'a, T> {
+impl<'a, T: Scalar> Iterator for PointIter<'a, T> {
     type Item = Point<'a, T>;
 
     fn next(&mut self) -> Option<Point<'a, T>> {

src/geometry/single_element/entity_geometry.rs

@@ -3,20 +3,20 @@ use super::SingleElementGeometry;
 use crate::{
     geometry::{Point, PointIter},
     traits::Geometry,
-    types::RealScalar,
+    types::Scalar,
 };
 use ndelement::traits::MappedFiniteElement;
 
 /// Geometry of an entity
 #[derive(Debug)]
-pub struct SingleElementEntityGeometry<'a, T: RealScalar, E: MappedFiniteElement> {
+pub struct SingleElementEntityGeometry<'a, T: Scalar, E: MappedFiniteElement> {
     geometry: &'a SingleElementGeometry<T, E>,
     cell_index: usize,
     sub_entity_dimension: usize,
     sub_entity_index: usize,
 }
 
-impl<'a, T: RealScalar, E: MappedFiniteElement> SingleElementEntityGeometry<'a, T, E> {
+impl<'a, T: Scalar, E: MappedFiniteElement> SingleElementEntityGeometry<'a, T, E> {
     /// Create new
     pub fn new(
         geometry: &'a SingleElementGeometry<T, E>,
@@ -33,7 +33,7 @@ impl<'a, T: RealScalar, E: MappedFiniteElement> SingleElementEntityGeometry<'a,
     }
 }
 
-impl<T: RealScalar, E: MappedFiniteElement> Geometry for SingleElementEntityGeometry<'_, T, E> {
+impl<T: Scalar, E: MappedFiniteElement> Geometry for SingleElementEntityGeometry<'_, T, E> {
     type T = T;
     type Point<'a>
         = Point<'a, T>

src/geometry/single_element/geometry.rs

@@ -1,7 +1,7 @@
 //! Geometry where each entity of a given dimension is represented by the same element
 #[cfg(feature = "serde")]
 use crate::traits::ConvertToSerializable;
-use crate::types::RealScalar;
+use crate::types::Scalar;
 #[cfg(feature = "serde")]
 use ndelement::{
     ciarlet::{CiarletElement, lagrange},
@@ -18,13 +18,13 @@ use rlst::{DynArray, rlst_dynamic_array};
 use std::fmt::{Debug, Formatter};
 
 /// Single element geometry
-pub struct SingleElementGeometry<T: RealScalar, E: MappedFiniteElement> {
+pub struct SingleElementGeometry<T: Scalar, E: MappedFiniteElement> {
     points: DynArray<T, 2>,
     cells: DynArray<usize, 2>,
     elements: Vec<E>,
 }
 
-impl<T: RealScalar, E: MappedFiniteElement> Debug for SingleElementGeometry<T, E> {
+impl<T: Scalar, E: MappedFiniteElement> Debug for SingleElementGeometry<T, E> {
     fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), std::fmt::Error> {
         f.debug_struct("SingleElementGeometry")
             .field("points", &self.points)
@@ -36,7 +36,7 @@ impl<T: RealScalar, E: MappedFiniteElement> Debug for SingleElementGeometry<T, E
 #[cfg(feature = "serde")]
 #[derive(serde::Serialize, Debug, serde::Deserialize)]
 #[serde(bound = "for<'de2> T: serde::Deserialize<'de2>")]
-pub struct SerializableGeometry<T: RealScalar + serde::Serialize>
+pub struct SerializableGeometry<T: Scalar + serde::Serialize>
 where
     for<'de2> T: serde::Deserialize<'de2>,
 {
@@ -46,8 +46,8 @@ where
 }
 
 #[cfg(feature = "serde")]
-impl<T: RealScalar + serde::Serialize> ConvertToSerializable
-    for SingleElementGeometry<T, CiarletElement<T, IdentityMap>>
+impl<T: Scalar + serde::Serialize> ConvertToSerializable
+    for SingleElementGeometry<T, CiarletElement<T, IdentityMap, T>>
 where
     for<'de2> T: serde::Deserialize<'de2>,
 {
@@ -86,7 +86,7 @@ where
     }
 }
 
-impl<T: RealScalar, E: MappedFiniteElement> SingleElementGeometry<T, E> {
+impl<T: Scalar, E: MappedFiniteElement> SingleElementGeometry<T, E> {
     /// Create single element geometry
     pub fn new(
         cell_type: ReferenceCellType,

src/grid/local_grid/mixed/builder.rs

@@ -5,7 +5,7 @@ use crate::{
     geometry::MixedGeometry,
     topology::mixed::MixedTopology,
     traits::{Builder, GeometryBuilder, GridBuilder, TopologyBuilder},
-    types::RealScalar,
+    types::Scalar,
 };
 use itertools::izip;
 use ndelement::{
@@ -20,10 +20,10 @@ use std::collections::{HashMap, HashSet};
 
 /// Grid builder for a grid with a mixture of element types
 #[derive(Debug)]
-pub struct MixedGridBuilder<T: RealScalar> {
+pub struct MixedGridBuilder<T: Scalar> {
     gdim: usize,
     element_indices: HashMap<(ReferenceCellType, usize), usize>,
-    elements: Vec<CiarletElement<T, IdentityMap>>,
+    elements: Vec<CiarletElement<T, IdentityMap, T>>,
     points_per_cell: Vec<usize>,
     pub(crate) points: Vec<T>,
     cells: Vec<usize>,
@@ -39,7 +39,7 @@ pub struct MixedGridBuilder<T: RealScalar> {
     cell_indices: HashSet<usize>,
 }
 
-impl<T: RealScalar> MixedGridBuilder<T> {
+impl<T: Scalar> MixedGridBuilder<T> {
     /// Create a new grid builder
     pub fn new(gdim: usize) -> Self {
         Self {
@@ -63,8 +63,8 @@ impl<T: RealScalar> MixedGridBuilder<T> {
     }
 }
 
-impl<T: RealScalar> Builder for MixedGridBuilder<T> {
-    type Grid = MixedGrid<T, CiarletElement<T, IdentityMap>>;
+impl<T: Scalar> Builder for MixedGridBuilder<T> {
+    type Grid = MixedGrid<T, CiarletElement<T, IdentityMap, T>>;
     type T = T;
     type CellData<'a> = (ReferenceCellType, usize, &'a [usize]);
     type EntityDescriptor = ReferenceCellType;
@@ -89,7 +89,7 @@ impl<T: RealScalar> Builder for MixedGridBuilder<T> {
             .entry((cell_data.0, cell_data.1))
             .or_insert_with(|| {
                 let n = self.cell_indices_to_ids.len();
-                self.elements.push(lagrange::create::<T>(
+                self.elements.push(lagrange::create::<T, T>(
                     cell_data.0,
                     cell_data.1,
                     Continuity::Standard,
@@ -127,7 +127,7 @@ impl<T: RealScalar> Builder for MixedGridBuilder<T> {
         self.add_cell(id, (cell_type, cell_degree, nodes));
     }
 
-    fn create_grid(&self) -> MixedGrid<T, CiarletElement<T, IdentityMap>> {
+    fn create_grid(&self) -> MixedGrid<T, CiarletElement<T, IdentityMap, T>> {
         let cell_vertices =
             self.extract_vertices(&self.cells, &self.cell_types, &self.cell_degrees);
 
@@ -209,16 +209,16 @@ impl<T: RealScalar> Builder for MixedGridBuilder<T> {
     }
 }
 
-impl<T: RealScalar> GeometryBuilder for MixedGridBuilder<T> {
-    type GridGeometry = MixedGeometry<T, CiarletElement<T, IdentityMap>>;
+impl<T: Scalar> GeometryBuilder for MixedGridBuilder<T> {
+    type GridGeometry = MixedGeometry<T, CiarletElement<T, IdentityMap, T>>;
     fn create_geometry(
         &self,
         point_ids: &[usize],
         coordinates: &[Self::T],
         cell_points: &[usize],
         cell_types: &[ReferenceCellType],
         cell_degrees: &[usize],
-    ) -> MixedGeometry<T, CiarletElement<T, IdentityMap>> {
+    ) -> MixedGeometry<T, CiarletElement<T, IdentityMap, T>> {
         let npts = point_ids.len();
         let mut points = rlst_dynamic_array!(T, [self.gdim(), npts]);
         points.data_mut().unwrap().copy_from_slice(coordinates);
@@ -229,15 +229,15 @@ impl<T: RealScalar> GeometryBuilder for MixedGridBuilder<T> {
         for degree in cell_degrees {
             cell_families.push(*ef_indices.entry(*degree).or_insert_with(|| {
                 let n = element_families.len();
-                element_families.push(LagrangeElementFamily::<T>::new(
+                element_families.push(LagrangeElementFamily::<T, T>::new(
                     *degree,
                     Continuity::Standard,
                 ));
                 n
             }))
         }
 
-        MixedGeometry::<T, CiarletElement<T, IdentityMap>>::new(
+        MixedGeometry::<T, CiarletElement<T, IdentityMap, T>>::new(
             cell_types,
             points,
             cell_points,
@@ -247,7 +247,7 @@ impl<T: RealScalar> GeometryBuilder for MixedGridBuilder<T> {
     }
 }
 
-impl<T: RealScalar> TopologyBuilder for MixedGridBuilder<T> {
+impl<T: Scalar> TopologyBuilder for MixedGridBuilder<T> {
     type GridTopology = MixedTopology;
     fn create_topology(
         &self,
@@ -297,7 +297,7 @@ impl<T: RealScalar> TopologyBuilder for MixedGridBuilder<T> {
     }
 }
 
-impl<T: RealScalar> GridBuilder for MixedGridBuilder<T> {
+impl<T: Scalar> GridBuilder for MixedGridBuilder<T> {
     fn create_grid_from_topology_geometry(
         &self,
         topology: <Self as TopologyBuilder>::GridTopology,