quadraturerules/
lib.rs

1//! Quadrature rules
2#![cfg_attr(feature = "strict", deny(warnings), deny(unused_crate_dependencies))]
3#![warn(missing_docs)]
4
5mod rules;
6
7/// A domain of an integral.
8#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
9#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)]
10#[repr(u8)]
11pub enum Domain {
12    /// interval
13    Interval = 0,
14    /// quadrilateral
15    Quadrilateral = 1,
16    /// edge-adjacent quadrilaterals
17    EdgeAdjacentQuadrilaterals = 2,
18    /// vertex-adjacent quadrilaterals
19    VertexAdjacentQuadrilaterals = 3,
20    /// triangle
21    Triangle = 4,
22    /// edge-adjacent triangle and quadrilateral
23    EdgeAdjacentTriangleAndQuadrilateral = 5,
24    /// vertex-adjacent triangle and quadrilateral
25    VertexAdjacentTriangleAndQuadrilateral = 6,
26    /// edge-adjacent triangles
27    EdgeAdjacentTriangles = 7,
28    /// vertex-adjacent triangles
29    VertexAdjacentTriangles = 8,
30    /// hexahedron
31    Hexahedron = 9,
32    /// tetrahedron
33    Tetrahedron = 10,
34    /// square-based pyramid
35    SquareBasedPyramid = 11,
36    /// triangular prism
37    TriangularPrism = 12,
38}
39
40/// A quadrature rule family.
41#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
42#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)]
43#[repr(u8)]
44pub enum QuadratureRule {
45    /// Gauss--Legendre
46    GaussLegendre = 1,
47    /// Gauss--Lobatto--Legendre
48    GaussLobattoLegendre = 3,
49    /// Hammer--Marlowe--Stroud
50    HammerMarloweStroud = 6,
51    /// Sauter--Schwab
52    SauterSchwab = 7,
53    /// Xiao--Gimbutas
54    XiaoGimbutas = 2,
55    /// centroid quadrature
56    CentroidQuadrature = 5,
57    /// closed Newton--Cotes
58    ClosedNewtonCotes = 8,
59    /// open Newton--Cotes
60    OpenNewtonCotes = 9,
61    /// vertex quadrature
62    VertexQuadrature = 4,
63}
64
65/// Get a quadrature rule for a single integral.
66pub fn single_integral_quadrature(
67    rtype: QuadratureRule,
68    domain: Domain,
69    order: usize,
70) -> Result<(Vec<f64>, Vec<f64>), &'static str> {
71    match rtype {
72        QuadratureRule::GaussLegendre => rules::gauss_legendre(domain, order),
73        QuadratureRule::GaussLobattoLegendre => rules::gauss_lobatto_legendre(domain, order),
74        QuadratureRule::HammerMarloweStroud => rules::hammer_marlowe_stroud(domain, order),
75        QuadratureRule::XiaoGimbutas => rules::xiao_gimbutas(domain, order),
76        QuadratureRule::CentroidQuadrature => rules::centroid_quadrature(domain, order),
77        QuadratureRule::ClosedNewtonCotes => rules::closed_newton_cotes(domain, order),
78        QuadratureRule::OpenNewtonCotes => rules::open_newton_cotes(domain, order),
79        QuadratureRule::VertexQuadrature => rules::vertex_quadrature(domain, order),
80        _ => Err("Unsupported rule for single integral"),
81    }
82}
83
84/// Get a quadrature rule for a double integral.
85pub fn double_integral_quadrature(
86    rtype: QuadratureRule,
87    domain: Domain,
88    order: usize,
89) -> Result<(Vec<f64>, Vec<f64>, Vec<f64>), &'static str> {
90    match rtype {
91        QuadratureRule::SauterSchwab => rules::sauter_schwab(domain, order),
92        _ => Err("Unsupported rule for double integral"),
93    }
94}
95
96#[cfg(test)]
97mod tests {
98    use super::*;
99    use approx::*;
100
101    #[test]
102    fn test_sum_weights_1() {
103        let weights =
104            single_integral_quadrature(QuadratureRule::GaussLegendre, Domain::Interval, 1)
105                .unwrap()
106                .1;
107        assert_relative_eq!(weights.iter().sum::<f64>(), 1.0);
108    }
109
110    #[test]
111    fn test_sum_weights_2() {
112        let weights =
113            single_integral_quadrature(QuadratureRule::GaussLegendre, Domain::Interval, 2)
114                .unwrap()
115                .1;
116        assert_relative_eq!(weights.iter().sum::<f64>(), 1.0);
117    }
118
119    #[test]
120    fn test_sum_weights_3() {
121        let weights =
122            single_integral_quadrature(QuadratureRule::GaussLegendre, Domain::Interval, 3)
123                .unwrap()
124                .1;
125        assert_relative_eq!(weights.iter().sum::<f64>(), 1.0);
126    }
127}
quadraturerules/lib.rs

quadraturerules/
lib.rs
