Implement CompactQuadratureTable

Author: Andlon

src/assembly/local/quadrature_table.rs

@@ -82,6 +82,42 @@ where
     }
 }
 
+/// 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 +129,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 +289,105 @@ 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, 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);
+    }
+}