Add Nedelec elements on simplices (#58)

* add nedelec on triangles and tets

* update python wrapper

* fmt

* ruff

* remove asserts from debugging

python/ndelement/ciarlet.py

@@ -21,6 +21,7 @@ class Family(Enum):
 
     Lagrange = 0
     RaviartThomas = 1
+    NedelecFirstKind = 2
 
 
 class MapType(Enum):
@@ -223,5 +224,12 @@ def create_family(
             )
         else:
             raise TypeError(f"Unsupported dtype: {dtype}")
+    elif family == Family.NedelecFirstKind:
+        if dtype == np.float64:
+            return ElementFamily(_lib.nedelec_element_family_new_f64(degree, continuity.value))
+        elif dtype == np.float32:
+            return ElementFamily(_lib.nedelec_element_family_new_f64(degree, continuity.value))
+        else:
+            raise TypeError(f"Unsupported dtype: {dtype}")
     else:
         raise ValueError(f"Unsupported family: {family}")

python/test/test_element.py

@@ -196,3 +196,49 @@ def test_raviart_thomas_1_triangle(continuity):
             assert element.entity_closure_dofs(1, i) == []
         assert element.entity_dofs(2, 0) == [0, 1, 2]
         assert element.entity_closure_dofs(2, 0) == [0, 1, 2]
+
+
+@pytest.mark.parametrize("continuity", [Continuity.Standard, Continuity.Discontinuous])
+def test_nedelec_1_triangle(continuity):
+    family = create_family(Family.NedelecFirstKind, 1, continuity=continuity)
+    element = family.element(ReferenceCellType.Triangle)
+
+    assert element.value_size == 2
+    assert element.value_shape == (2,)
+    assert element.degree == 1
+    assert element.embedded_superdegree == 1
+    assert element.dim == 3
+    assert element.continuity == continuity
+    assert element.map_type == MapType.CovariantPiola
+    assert element.cell_type == ReferenceCellType.Triangle
+
+    if continuity == Continuity.Standard:
+        assert element.entity_dofs(0, 0) == []
+        assert element.entity_dofs(0, 1) == []
+        assert element.entity_dofs(0, 2) == []
+        assert element.entity_dofs(1, 0) == [0]
+        assert element.entity_dofs(1, 1) == [1]
+        assert element.entity_dofs(1, 2) == [2]
+        assert element.entity_dofs(2, 0) == []
+
+        assert element.entity_closure_dofs(0, 0) == []
+        assert element.entity_closure_dofs(0, 1) == []
+        assert element.entity_closure_dofs(0, 2) == []
+        assert element.entity_closure_dofs(1, 0) == [0]
+        assert element.entity_closure_dofs(1, 1) == [1]
+        assert element.entity_closure_dofs(1, 2) == [2]
+        assert element.entity_closure_dofs(2, 0) == [0, 1, 2]
+
+        ip = element.interpolation_points()
+        assert ip[0][0].shape == (0, 2)
+        assert ip[0][1].shape == (0, 2)
+        assert ip[0][2].shape == (0, 2)
+        assert ip[2][0].shape == (0, 2)
+    else:
+        for i in range(3):
+            assert element.entity_dofs(0, i) == []
+            assert element.entity_dofs(1, i) == []
+            assert element.entity_closure_dofs(0, i) == []
+            assert element.entity_closure_dofs(1, i) == []
+        assert element.entity_dofs(2, 0) == [0, 1, 2]
+        assert element.entity_closure_dofs(2, 0) == [0, 1, 2]

src/bindings.rs

@@ -262,6 +262,7 @@ pub mod ciarlet {
     pub enum ElementType {
         Lagrange = 0,
         RaviartThomas = 1,
+        NedelecFirstKind = 2,
     }
 
     #[repr(C)]
@@ -339,6 +340,20 @@ pub mod ciarlet {
                         Box::from_raw(*family as *mut ciarlet::RaviartThomasElementFamily<c64>)
                     }),
                 },
+                ElementType::NedelecFirstKind => match dtype {
+                    DType::F32 => drop(unsafe {
+                        Box::from_raw(*family as *mut ciarlet::NedelecFirstKindElementFamily<f32>)
+                    }),
+                    DType::F64 => drop(unsafe {
+                        Box::from_raw(*family as *mut ciarlet::NedelecFirstKindElementFamily<f64>)
+                    }),
+                    DType::C32 => drop(unsafe {
+                        Box::from_raw(*family as *mut ciarlet::NedelecFirstKindElementFamily<c32>)
+                    }),
+                    DType::C64 => drop(unsafe {
+                        Box::from_raw(*family as *mut ciarlet::NedelecFirstKindElementFamily<c64>)
+                    }),
+                },
             }
         }
     }
@@ -832,6 +847,24 @@ pub mod ciarlet {
                     ciarlet::RaviartThomasElementFamily<c64>,
                 >(family, cell),
             },
+            ElementType::NedelecFirstKind => match (*family).dtype {
+                DType::F32 => element_family_element_internal::<
+                    f32,
+                    ciarlet::NedelecFirstKindElementFamily<f32>,
+                >(family, cell),
+                DType::F64 => element_family_element_internal::<
+                    f64,
+                    ciarlet::NedelecFirstKindElementFamily<f64>,
+                >(family, cell),
+                DType::C32 => element_family_element_internal::<
+                    c32,
+                    ciarlet::NedelecFirstKindElementFamily<c32>,
+                >(family, cell),
+                DType::C64 => element_family_element_internal::<
+                    c64,
+                    ciarlet::NedelecFirstKindElementFamily<c64>,
+                >(family, cell),
+            },
         }
     }
 
@@ -902,4 +935,42 @@ pub mod ciarlet {
             dtype: DType::F64,
         }))
     }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn nedelec_element_family_new_f32(
+        degree: usize,
+        continuity: u8,
+    ) -> *const ElementFamilyWrapper {
+        let family = Box::into_raw(Box::new(
+            ciarlet::NedelecFirstKindElementFamily::<f32>::new(
+                degree,
+                Continuity::from(continuity).expect("Invalid continuity"),
+            ),
+        )) as *mut c_void;
+
+        Box::into_raw(Box::new(ElementFamilyWrapper {
+            family,
+            etype: ElementType::NedelecFirstKind,
+            dtype: DType::F32,
+        }))
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn nedelec_element_family_new_f64(
+        degree: usize,
+        continuity: u8,
+    ) -> *const ElementFamilyWrapper {
+        let family = Box::into_raw(Box::new(
+            ciarlet::NedelecFirstKindElementFamily::<f64>::new(
+                degree,
+                Continuity::from(continuity).expect("Invalid continuity"),
+            ),
+        )) as *mut c_void;
+
+        Box::into_raw(Box::new(ElementFamilyWrapper {
+            family,
+            etype: ElementType::NedelecFirstKind,
+            dtype: DType::F64,
+        }))
+    }
 }

src/ciarlet.rs

@@ -12,8 +12,10 @@ use rlst::{
 use std::fmt::{Debug, Formatter};
 
 pub mod lagrange;
+pub mod nedelec;
 pub mod raviart_thomas;
 pub use lagrange::LagrangeElementFamily;
+pub use nedelec::NedelecFirstKindElementFamily;
 pub use raviart_thomas::RaviartThomasElementFamily;
 
 type EntityPoints<T> = [Vec<Array<T, BaseArray<T, VectorContainer<T>, 2>, 2>>; 4];
@@ -972,6 +974,13 @@ mod test {
         check_dofs(e);
     }
 
+    #[test]
+    fn test_raviart_thomas_3_triangle() {
+        let e = raviart_thomas::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
+        assert_eq!(e.value_size(), 2);
+        check_dofs(e);
+    }
+
     #[test]
     fn test_raviart_thomas_1_tetrahedron() {
         let e =
@@ -988,6 +997,104 @@ mod test {
         check_dofs(e);
     }
 
+    #[test]
+    fn test_raviart_thomas_3_tetrahedron() {
+        let e =
+            raviart_thomas::create::<f64>(ReferenceCellType::Tetrahedron, 3, 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);
+        assert_eq!(e.value_size(), 2);
+        let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
+        *points.get_mut([0, 0]).unwrap() = 0.0;
+        *points.get_mut([1, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).unwrap() = 1.0;
+        *points.get_mut([1, 1]).unwrap() = 0.0;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([0, 4]).unwrap() = 0.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
+        e.tabulate(&points, 0, &mut data);
+
+        for pt in 0..6 {
+            assert_relative_eq!(
+                *data.get([0, pt, 0, 0]).unwrap(),
+                -*points.get([1, pt]).unwrap(),
+                epsilon = 1e-14
+            );
+            assert_relative_eq!(
+                *data.get([0, pt, 0, 1]).unwrap(),
+                *points.get([0, pt]).unwrap(),
+                epsilon = 1e-14
+            );
+            assert_relative_eq!(
+                *data.get([0, pt, 1, 0]).unwrap(),
+                *points.get([1, pt]).unwrap(),
+                epsilon = 1e-14
+            );
+            assert_relative_eq!(
+                *data.get([0, pt, 1, 1]).unwrap(),
+                1.0 - *points.get([0, pt]).unwrap(),
+                epsilon = 1e-14
+            );
+            assert_relative_eq!(
+                *data.get([0, pt, 2, 0]).unwrap(),
+                1.0 - *points.get([1, pt]).unwrap(),
+                epsilon = 1e-14
+            );
+            assert_relative_eq!(
+                *data.get([0, pt, 2, 1]).unwrap(),
+                *points.get([0, pt]).unwrap(),
+                epsilon = 1e-14
+            );
+        }
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_nedelec_2_triangle() {
+        let e = nedelec::create::<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);
+        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);
+        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);
+        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);
+        assert_eq!(e.value_size(), 3);
+        check_dofs(e);
+    }
+
     macro_rules! test_entity_closure_dofs_lagrange {
         ($cell:ident, $degree:expr) => {
             paste! {