bempp 0.2.0

Boundary element method library.
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

//! Get rules on simplices.

use crate::quadrature::simplex_rule_definitions::SIMPLEX_RULE_DEFINITIONS;
use crate::quadrature::types::NumericalQuadratureDefinition;
use crate::quadrature::types::QuadratureError;
use ndelement::types::ReferenceCellType;

/// Return a simplex rule for a given number of points.
///
/// If the rule does not exist `Err(())` is returned.
pub fn simplex_rule(
    cell_type: ReferenceCellType,
    npoints: usize,
) -> Result<NumericalQuadratureDefinition, QuadratureError> {
    let dim: usize = match cell_type {
        ReferenceCellType::Point => 0,
        ReferenceCellType::Interval => 1,
        ReferenceCellType::Triangle => 2,
        ReferenceCellType::Quadrilateral => 2,
        ReferenceCellType::Tetrahedron => 3,
        ReferenceCellType::Hexahedron => 3,
        ReferenceCellType::Prism => 3,
        ReferenceCellType::Pyramid => 3,
    };

    if let Some((order, points, weights)) = SIMPLEX_RULE_DEFINITIONS
        .get(&cell_type)
        .unwrap()
        .get(&npoints)
    {
        Ok(NumericalQuadratureDefinition {
            dim,
            order: *order,
            npoints,
            weights: weights.to_vec(),
            points: points.to_vec(),
        })
    } else {
        Err(QuadratureError::RuleNotFound)
    }
}

/// For a given cell type return a vector with the numbers of points for which simplex rules are available.
pub fn available_rules(cell_type: ReferenceCellType) -> Vec<usize> {
    SIMPLEX_RULE_DEFINITIONS
        .get(&cell_type)
        .unwrap()
        .iter()
        .map(|(npoints, _)| *npoints)
        .collect()
}

#[cfg(test)]
mod test {

    use super::*;
    use paste::paste;

    use approx::*;

    fn get_volume(cell_type: ReferenceCellType) -> f64 {
        match cell_type {
            ReferenceCellType::Point => 0.0,
            ReferenceCellType::Interval => 1.0,
            ReferenceCellType::Triangle => 0.5,
            ReferenceCellType::Quadrilateral => 1.0,
            ReferenceCellType::Tetrahedron => 1.0 / 6.0,
            ReferenceCellType::Hexahedron => 1.0,
            ReferenceCellType::Prism => 0.5,
            ReferenceCellType::Pyramid => 1.0 / 3.0,
        }
    }

    macro_rules! test_cell {

        ($($cell:ident),+) => {

        $(
            paste! {

                #[test]
                fn [<test_volume_ $cell:lower>]() {
                    let cell_type = ReferenceCellType::[<$cell>];
                    let rules = available_rules(cell_type);
                    for npoints in rules {
                        let rule = simplex_rule(cell_type, npoints).unwrap();
                        let volume_actual: f64 = rule.weights.iter().sum();
                        let volume_expected = get_volume(cell_type);
                        assert_relative_eq!(volume_actual, volume_expected, max_relative=1E-14);
                        }

                }

            }
        )*
        };
    }

    test_cell!(
        Interval,
        Triangle,
        Quadrilateral,
        Tetrahedron,
        Hexahedron,
        Prism,
        Pyramid
    );
}