working on RT element

element/src/element/raviart_thomas.rs

@@ -1,5 +1,6 @@
 //! Raviart-Thomas elements
 
+use crate::element::CiarletElement;
 use crate::element::OldCiarletElement;
 use bempp_tools::arrays::{AdjacencyList, Array3D};
 use bempp_traits::arrays::Array3DAccess;
@@ -84,6 +85,98 @@ pub fn create(
     }
 }
 
+/// Create a Raviart-Thomas element
+pub fn create_new(
+    cell_type: ReferenceCellType,
+    degree: usize,
+    continuity: Continuity,
+) -> CiarletElement {
+    if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Quadrilateral {
+        panic!("Unsupported cell type");
+    }
+
+    if cell_type != ReferenceCellType::Triangle {
+        panic!("RT elements on quadrilaterals not implemented yet");
+    }
+    if degree != 1 {
+        panic!("Degree > 1 RT elements not implemented yet");
+    }
+
+    let cell = create_cell(cell_type);
+    let pdim = polynomial_count(cell_type, degree);
+    let tdim = cell.dim();
+    let edim = tdim * polynomial_count(cell_type, degree - 1) + degree;
+
+    let mut wcoeffs = Array3D::<f64>::new((edim, tdim, pdim));
+
+    // [sqrt(2), 6*y - 2, 4*sqrt(3)*(x + y/2 - 1/2)]
+
+    // norm(x**2 + y**2)
+    // sqrt(70)/30
+
+    *wcoeffs.get_mut(0, 0, 0) = 1.0;
+    *wcoeffs.get_mut(1, 1, 0) = 1.0;
+
+    let mut x = [vec![], vec![], vec![], vec![]];
+    let mut m = [vec![], vec![], vec![], vec![]];
+    for _e in 0..cell.entity_count(0) {
+        x[0].push(Array2D::<f64>::new((0, tdim)));
+        m[0].push(Array3D::<f64>::new((0, 2, 0)));
+    }
+
+    for _e in 0..cell.entity_count(1) {
+        let mut pts = vec![0.0; tdim];
+        let mut mat = vec![0.0; 2];
+        let vn0 = cell.edges()[2 * e];
+        let vn1 = cell.edges()[2 * e + 1];
+        let v0 = &cell.vertices()[vn0 * tdim..(vn0 + 1) * tdim];
+        let v1 = &cell.vertices()[vn1 * tdim..(vn1 + 1) * tdim];
+        for i in 0..tdim {
+            pts[i] = (v0[i] + v1[i]) / 2.0;
+        }
+        x[1].push(Array2D::<f64>::from_data(pts, (1, tdim)));
+        m[1].push(Array3D::<f64>::from_data(mat, (1, 2, 1)));
+    }
+
+    for _e in 0..cell.entity_count(2) {
+        x[2].push(Array2D::<f64>::new((0, tdim)));
+        m[2].push(Array3D::<f64>::new((0, 2, 0)));
+    }
+
+
+    let coefficients = match cell_type {
+        ReferenceCellType::Triangle => match degree {
+            // Basis = {(-x, -y), (x-1,y), (-x, 1-y)}
+            1 => Array3D::from_data(
+                vec![
+                    0.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0,
+                    1.0, 0.0, -1.0,
+                ],
+                (3, 2, 3),
+            ),
+            _ => {
+                panic!("Degree not supported");
+            }
+        },
+        _ => {
+            panic!("Cell type not supported");
+        }
+    };
+
+    CiarletElement::create(
+        ElementFamily::RaviartThomas
+        cell_type,
+        degree,
+        vec![2],
+        wcoeffs,
+        x,
+        m,
+        MapType::ContravariantPiola,
+        continuity,
+        degree,
+    }
+}
+
 #[cfg(test)]
 mod test {
     use crate::cell::*;
@@ -124,7 +217,7 @@ mod test {
     }
 
     #[test]
-    fn test_raviart_thomas_1_triangle() {
+    fn test_oldraviart_thomas_1_triangle() {
         let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous);
         assert_eq!(e.value_size(), 2);
         let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
@@ -159,4 +252,41 @@ mod test {
         }
         check_dofs(e);
     }
+
+    #[test]
+    fn test_raviart_thomas_1_triangle() {
+        let e = create_new(ReferenceCellType::Triangle, 1, Continuity::Continuous);
+        assert_eq!(e.value_size(), 2);
+        let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
+        let points = Array2D::from_data(
+            vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
+            (6, 2),
+        );
+        e.tabulate(&points, 0, &mut data);
+
+        for pt in 0..6 {
+            assert_relative_eq!(
+                *data.get(0, pt, 0, 0).unwrap(),
+                -*points.get(pt, 0).unwrap()
+            );
+            assert_relative_eq!(
+                *data.get(0, pt, 0, 1).unwrap(),
+                -*points.get(pt, 1).unwrap()
+            );
+            assert_relative_eq!(
+                *data.get(0, pt, 1, 0).unwrap(),
+                *points.get(pt, 0).unwrap() - 1.0
+            );
+            assert_relative_eq!(*data.get(0, pt, 1, 1).unwrap(), *points.get(pt, 1).unwrap());
+            assert_relative_eq!(
+                *data.get(0, pt, 2, 0).unwrap(),
+                -*points.get(pt, 0).unwrap()
+            );
+            assert_relative_eq!(
+                *data.get(0, pt, 2, 1).unwrap(),
+                1.0 - *points.get(pt, 1).unwrap()
+            );
+        }
+        check_dofs(e);
+    }
 }

tools/src/arrays.rs

@@ -6,7 +6,7 @@ use std::clone::Clone;
 /// A two-dimensional rectangular array
 pub struct Array2D<T: Num> {
     /// The data in the array, in row-major order
-    pub data: Vec<T>,
+    data: Vec<T>,
     /// The shape of the array
     shape: (usize, usize),
 }
@@ -226,7 +226,7 @@ impl<T: Num> Array4DAccess<T> for Array4D<T> {
 /// An adjacency list stores two-dimensional data where each row may have a different number of items
 pub struct AdjacencyList<T: Num> {
     /// The data in the array, in row-major order
-    pub data: Vec<T>,
+    data: Vec<T>,
     /// The starting index of each row, plus a final entry that is the length of data
     offsets: Vec<usize>,
 }