quadraturerules 0.9.1

quadrature rules library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

//! Quadrature rules
#![cfg_attr(feature = "strict", deny(warnings), deny(unused_crate_dependencies))]
#![warn(missing_docs)]

mod rules;

/// A domain of an integral.
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)]
#[repr(u8)]
pub enum Domain {
    /// interval
    Interval = 0,
    /// quadrilateral
    Quadrilateral = 1,
    /// edge-adjacent quadrilaterals
    EdgeAdjacentQuadrilaterals = 2,
    /// vertex-adjacent quadrilaterals
    VertexAdjacentQuadrilaterals = 3,
    /// triangle
    Triangle = 4,
    /// edge-adjacent triangle and quadrilateral
    EdgeAdjacentTriangleAndQuadrilateral = 5,
    /// vertex-adjacent triangle and quadrilateral
    VertexAdjacentTriangleAndQuadrilateral = 6,
    /// edge-adjacent triangles
    EdgeAdjacentTriangles = 7,
    /// vertex-adjacent triangles
    VertexAdjacentTriangles = 8,
    /// hexahedron
    Hexahedron = 9,
    /// tetrahedron
    Tetrahedron = 10,
    /// square-based pyramid
    SquareBasedPyramid = 11,
    /// triangular prism
    TriangularPrism = 12,
}

/// A quadrature rule family.
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)]
#[repr(u8)]
pub enum QuadratureRule {
    /// Gauss--Legendre
    GaussLegendre = 1,
    /// Gauss--Lobatto--Legendre
    GaussLobattoLegendre = 3,
    /// Hammer--Marlowe--Stroud
    HammerMarloweStroud = 6,
    /// Sauter--Schwab
    SauterSchwab = 7,
    /// Xiao--Gimbutas
    XiaoGimbutas = 2,
    /// centroid quadrature
    CentroidQuadrature = 5,
    /// closed Newton--Cotes
    ClosedNewtonCotes = 8,
    /// open Newton--Cotes
    OpenNewtonCotes = 9,
    /// vertex quadrature
    VertexQuadrature = 4,
}

/// Get a quadrature rule for a single integral.
pub fn single_integral_quadrature(
    rtype: QuadratureRule,
    domain: Domain,
    order: usize,
) -> Result<(Vec<f64>, Vec<f64>), &'static str> {
    match rtype {
        QuadratureRule::GaussLegendre => rules::gauss_legendre(domain, order),
        QuadratureRule::GaussLobattoLegendre => rules::gauss_lobatto_legendre(domain, order),
        QuadratureRule::HammerMarloweStroud => rules::hammer_marlowe_stroud(domain, order),
        QuadratureRule::XiaoGimbutas => rules::xiao_gimbutas(domain, order),
        QuadratureRule::CentroidQuadrature => rules::centroid_quadrature(domain, order),
        QuadratureRule::ClosedNewtonCotes => rules::closed_newton_cotes(domain, order),
        QuadratureRule::OpenNewtonCotes => rules::open_newton_cotes(domain, order),
        QuadratureRule::VertexQuadrature => rules::vertex_quadrature(domain, order),
        _ => Err("Unsupported rule for single integral"),
    }
}

/// Get a quadrature rule for a double integral.
pub fn double_integral_quadrature(
    rtype: QuadratureRule,
    domain: Domain,
    order: usize,
) -> Result<(Vec<f64>, Vec<f64>, Vec<f64>), &'static str> {
    match rtype {
        QuadratureRule::SauterSchwab => rules::sauter_schwab(domain, order),
        _ => Err("Unsupported rule for double integral"),
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::*;

    #[test]
    fn test_sum_weights_1() {
        let weights =
            single_integral_quadrature(QuadratureRule::GaussLegendre, Domain::Interval, 1)
                .unwrap()
                .1;
        assert_relative_eq!(weights.iter().sum::<f64>(), 1.0);
    }

    #[test]
    fn test_sum_weights_2() {
        let weights =
            single_integral_quadrature(QuadratureRule::GaussLegendre, Domain::Interval, 2)
                .unwrap()
                .1;
        assert_relative_eq!(weights.iter().sum::<f64>(), 1.0);
    }

    #[test]
    fn test_sum_weights_3() {
        let weights =
            single_integral_quadrature(QuadratureRule::GaussLegendre, Domain::Interval, 3)
                .unwrap()
                .1;
        assert_relative_eq!(weights.iter().sum::<f64>(), 1.0);
    }
}