Merge pull request #43 from Andlon/compact-quadrature-table

Compact quadrature table

Author: Andlon

examples/poisson2d.rs

@@ -1,6 +1,4 @@
 use eyre::eyre;
-use nalgebra::Vector1;
-
 use fenris::assembly::global::{
     apply_homogeneous_dirichlet_bc_csr, apply_homogeneous_dirichlet_bc_rhs, CsrAssembler, VectorAssembler,
 };
@@ -14,6 +12,7 @@ use fenris::mesh::QuadMesh2d;
 use fenris::nalgebra::{DMatrix, DVector, Point2, U1, U2};
 use fenris::nalgebra_sparse::CsrMatrix;
 use fenris::quadrature;
+use nalgebra::Vector1;
 
 fn main() -> eyre::Result<()> {
     // TODO: Make it easy to construct triangle meshes as well.

src/assembly/local/quadrature_table.rs

@@ -67,21 +67,64 @@ where
     data: NestedVec<Data>,
 }
 
+fn unit_data_table_for_weights<T>(points: &NestedVec<T>) -> NestedVec<()> {
+    let mut data = NestedVec::new();
+    for i in 0..points.len() {
+        data.push(&vec![(); points.get(i).unwrap().len()]);
+    }
+    data
+}
+
 impl<T, GeometryDim> GeneralQuadratureTable<T, GeometryDim>
 where
     T: Scalar,
     GeometryDim: DimName,
     DefaultAllocator: Allocator<T, GeometryDim>,
 {
     pub fn from_points_and_weights(points: NestedVec<OPoint<T, GeometryDim>>, weights: NestedVec<T>) -> Self {
-        let mut data = NestedVec::new();
-        for i in 0..points.len() {
-            data.push(&vec![(); points.get(i).unwrap().len()]);
-        }
+        let data = unit_data_table_for_weights(&weights);
         Self::from_points_weights_and_data(points, weights, data)
     }
 }
 
+/// Checks that the provided quadrature rules are consistent, in the sense that
+/// the number of elements for each table is identical, and that each rule has
+/// consistent numbers of points, weights and data entries.
+fn check_rules_consistency<T, D, Data>(points: &NestedVec<OPoint<T, D>>, weights: &NestedVec<T>, data: &NestedVec<Data>)
+where
+    T: Scalar,
+    D: DimName,
+    DefaultAllocator: Allocator<T, D>,
+{
+    assert_eq!(
+        points.len(),
+        weights.len(),
+        "Quadrature point and weight tables must have the same number of rules."
+    );
+    assert_eq!(
+        points.len(),
+        data.len(),
+        "Quadrature point and data tables must have the same number of rules."
+    );
+
+    // Ensure that each element has a consistent quadrature rule
+    let iter = izip!(points.iter(), weights.iter(), data.iter());
+    for (element_index, (element_points, element_weights, element_data)) in iter.enumerate() {
+        assert_eq!(
+            element_points.len(),
+            element_weights.len(),
+            "Element {} has mismatched number of points and weights.",
+            element_index
+        );
+        assert_eq!(
+            element_points.len(),
+            element_data.len(),
+            "Element {} has mismatched number of points and data.",
+            element_index
+        );
+    }
+}
+
 impl<T, GeometryDim, Data> GeneralQuadratureTable<T, GeometryDim, Data>
 where
     T: Scalar,
@@ -93,26 +136,7 @@ where
         weights: NestedVec<T>,
         data: NestedVec<Data>,
     ) -> Self {
-        assert_eq!(points.len(), weights.len());
-        assert_eq!(points.len(), data.len());
-
-        // Ensure that each element has a consistent quadrature rule
-        let iter = izip!(points.iter(), weights.iter(), data.iter());
-        for (element_index, (element_points, element_weights, element_data)) in iter.enumerate() {
-            assert_eq!(
-                element_points.len(),
-                element_weights.len(),
-                "Element {} has mismatched number of points and weights.",
-                element_index
-            );
-            assert_eq!(
-                element_points.len(),
-                element_data.len(),
-                "Element {} has mismatched number of points and data.",
-                element_index
-            );
-        }
-
+        check_rules_consistency(&points, &weights, &data);
         Self { points, weights, data }
     }
 
@@ -272,3 +296,144 @@ where
         weights.clone_from_slice(&self.weights);
     }
 }
+
+/// A general quadrature table that avoids duplication of identical rules.
+///
+/// In a nutshell, [`CompactQuadratureTable`] sits in between [`UniformQuadratureTable`]
+/// and [`GeneralQuadratureTable`]. Like [`GeneralQuadratureTable`], it can store an arbitrary
+/// rule per element, but it uses a layer of indirection so that `M` quadrature rules
+/// are associated with `N` elements.
+///
+/// This can be useful in settings where many elements share the same quadrature data, such
+/// as a finite element space with elements of different degrees, or the common case
+/// where the elements are the same but different quadrature data is needed in different
+/// regions of the mesh.
+#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
+pub struct CompactQuadratureTable<T, D, Data = ()>
+where
+    T: Scalar,
+    D: DimName,
+    DefaultAllocator: Allocator<T, D>,
+{
+    #[serde(bound(serialize = "OPoint<T, D>: Serialize"))]
+    #[serde(bound(deserialize = "OPoint<T, D>: Deserialize<'de>"))]
+    points: NestedVec<OPoint<T, D>>,
+    weights: NestedVec<T>,
+    data: NestedVec<Data>,
+    element_to_rule_map: Vec<usize>,
+}
+
+impl<T, D> CompactQuadratureTable<T, D>
+where
+    T: Scalar,
+    D: DimName,
+    DefaultAllocator: Allocator<T, D>,
+{
+    pub fn from_points_weights_and_map(
+        points: NestedVec<OPoint<T, D>>,
+        weights: NestedVec<T>,
+        element_to_rule_map: Vec<usize>,
+    ) -> Self {
+        let data = unit_data_table_for_weights(&weights);
+        Self::from_quadrature_rules_and_map(points, weights, data, element_to_rule_map)
+    }
+}
+
+impl<T, D, Data> CompactQuadratureTable<T, D, Data>
+where
+    T: Scalar,
+    D: DimName,
+    DefaultAllocator: Allocator<T, D>,
+{
+    /// Construct a new table from the given quadrature rules and a map from elements
+    /// to quadrature rules.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `points`, `weights` and `data` are not consistent with each other.
+    ///
+    /// Panics if the mapping from elements to quadrature rules contains indices that are
+    /// out of bounds with respect to the number of quadrature rules.
+    pub fn from_quadrature_rules_and_map(
+        points: NestedVec<OPoint<T, D>>,
+        weights: NestedVec<T>,
+        data: NestedVec<Data>,
+        element_to_rule_map: Vec<usize>,
+    ) -> Self {
+        check_rules_consistency(&points, &weights, &data);
+        let num_rules = points.len();
+        let rule_indices_in_bounds = element_to_rule_map
+            .iter()
+            .all(|rule_index| rule_index < &num_rules);
+        assert!(
+            rule_indices_in_bounds,
+            "Each rule index must correspond to a provided quadrature rule."
+        );
+        Self {
+            element_to_rule_map,
+            points,
+            weights,
+            data,
+        }
+    }
+
+    fn rule_index_for_element(&self, element_index: usize) -> usize {
+        self.element_to_rule_map[element_index]
+    }
+}
+
+impl<T, D, Data> QuadratureTable<T, D> for CompactQuadratureTable<T, D, Data>
+where
+    T: Scalar,
+    D: SmallDim,
+    Data: Default + Clone,
+    DefaultAllocator: Allocator<T, D>,
+{
+    type Data = Data;
+
+    fn element_quadrature_size(&self, element_index: usize) -> usize {
+        let rule_index = self.rule_index_for_element(element_index);
+        self.points
+            .get(rule_index)
+            .expect("Internal error: Rule index out of bounds")
+            .len()
+    }
+
+    fn populate_element_data(&self, element_index: usize, data: &mut [Self::Data]) {
+        let rule_index = self.rule_index_for_element(element_index);
+        let data_array = self
+            .data
+            .get(rule_index)
+            .expect("Internal error: Rule index out of bounds");
+        assert_eq!(
+            data.len(),
+            data_array.len(),
+            "Length mismatch in data array: Stored quadrature data array has different length than output array."
+        );
+        data.clone_from_slice(data_array);
+    }
+
+    fn populate_element_quadrature(&self, element_index: usize, points: &mut [OPoint<T, D>], weights: &mut [T]) {
+        let rule_index = self.rule_index_for_element(element_index);
+        let points_array = self
+            .points
+            .get(rule_index)
+            .expect("Internal error: Rule index out of bounds");
+        let weights_array = self
+            .weights
+            .get(rule_index)
+            .expect("Internal error: Rule index out of bounds");
+        assert_eq!(
+            points.len(),
+            points_array.len(),
+            "Length mismatch in points array: Stored points array has different length than output array."
+        );
+        assert_eq!(
+            weights.len(),
+            weights_array.len(),
+            "Length mismatch in points array: Stored points array has different length than output array."
+        );
+        points.clone_from_slice(points_array);
+        weights.clone_from_slice(weights_array);
+    }
+}

src/element/tetrahedron.rs

@@ -1,6 +1,6 @@
 use numeric_literals::replace_float_literals;
 
-use crate::connectivity::{Tet10Connectivity, Tet4Connectivity};
+use crate::connectivity::{Tet10Connectivity, Tet20Connectivity, Tet4Connectivity};
 use crate::element::{ElementConnectivity, FiniteElement, FixedNodesReferenceFiniteElement};
 use crate::nalgebra::{
     distance, Matrix1x4, Matrix3, Matrix3x4, OMatrix, OPoint, Point3, Scalar, Vector3, U1, U10, U20, U3, U4,
@@ -52,6 +52,30 @@ where
     }
 }
 
+impl<T> ElementConnectivity<T> for Tet20Connectivity
+where
+    T: Real,
+{
+    type Element = Tet20Element<T>;
+    type GeometryDim = U3;
+    type ReferenceDim = U3;
+
+    fn element(&self, vertices: &[OPoint<T, Self::GeometryDim>]) -> Option<Self::Element> {
+        let mut tet20_vertices = [Point3::origin(); 20];
+        for (i, v) in tet20_vertices.iter_mut().enumerate() {
+            *v = vertices.get(self.0[i])?.clone();
+        }
+
+        let mut tet4_vertices = [Point3::origin(); 4];
+        tet4_vertices.copy_from_slice(&tet20_vertices[0..4]);
+
+        Some(Tet20Element {
+            tet4: Tet4Element::from_vertices(tet4_vertices),
+            vertices: tet20_vertices,
+        })
+    }
+}
+
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub struct Tet10Element<T>
 where