ndelement/
quadrature.rs

1//! Quadrature.
2use crate::types::ReferenceCellType;
3use bempp_quadrature::{
4    gauss_jacobi, simplex_rules,
5    types::{NumericalQuadratureDefinition, QuadratureError},
6};
7
8/// Return a simplex rule for a given number of points.
9///
10/// If the rule does not exist `Err(())` is returned.
11pub fn simplex_rule(
12    cell: ReferenceCellType,
13    npoints: usize,
14) -> Result<NumericalQuadratureDefinition, QuadratureError> {
15    match cell {
16        ReferenceCellType::Point => {
17            if npoints == 1 {
18                Ok(NumericalQuadratureDefinition {
19                    dim: 0,
20                    order: 1,
21                    npoints,
22                    weights: vec![1.0],
23                    points: vec![],
24                })
25            } else {
26                Err(QuadratureError::RuleNotFound)
27            }
28        }
29        ReferenceCellType::Interval => simplex_rules::simplex_rule_interval(npoints),
30        ReferenceCellType::Triangle => simplex_rules::simplex_rule_triangle(npoints),
31        ReferenceCellType::Quadrilateral => simplex_rules::simplex_rule_quadrilateral(npoints),
32        ReferenceCellType::Tetrahedron => simplex_rules::simplex_rule_tetrahedron(npoints),
33        ReferenceCellType::Hexahedron => simplex_rules::simplex_rule_hexahedron(npoints),
34        ReferenceCellType::Prism => simplex_rules::simplex_rule_prism(npoints),
35        ReferenceCellType::Pyramid => simplex_rules::simplex_rule_pyramid(npoints),
36    }
37}
38
39/// For a given cell type return a vector with the numbers of points for which simplex rules are available.
40pub fn available_simplex_rules(cell: ReferenceCellType) -> Vec<usize> {
41    match cell {
42        ReferenceCellType::Point => vec![1],
43        ReferenceCellType::Interval => simplex_rules::available_rules_interval(),
44        ReferenceCellType::Triangle => simplex_rules::available_rules_triangle(),
45        ReferenceCellType::Quadrilateral => simplex_rules::available_rules_quadrilateral(),
46        ReferenceCellType::Tetrahedron => simplex_rules::available_rules_tetrahedron(),
47        ReferenceCellType::Hexahedron => simplex_rules::available_rules_hexahedron(),
48        ReferenceCellType::Prism => simplex_rules::available_rules_prism(),
49        ReferenceCellType::Pyramid => simplex_rules::available_rules_pyramid(),
50    }
51}
52
53/// Get the points and weights for a Gauss-Jacobi quadrature rule of degree `m` on the given cell type.
54pub fn gauss_jacobi_rule(
55    celltype: ReferenceCellType,
56    m: usize,
57) -> Result<NumericalQuadratureDefinition, QuadratureError> {
58    let np = (m + 2) / 2;
59    match celltype {
60        ReferenceCellType::Interval => Ok(gauss_jacobi::gauss_jacobi_interval(np)),
61        ReferenceCellType::Triangle => Ok(gauss_jacobi::gauss_jacobi_triangle(np)),
62        ReferenceCellType::Quadrilateral => Ok(gauss_jacobi::gauss_jacobi_quadrilateral(np)),
63        ReferenceCellType::Tetrahedron => Ok(gauss_jacobi::gauss_jacobi_tetrahedron(np)),
64        ReferenceCellType::Hexahedron => Ok(gauss_jacobi::gauss_jacobi_hexahedron(np)),
65        _ => Err(QuadratureError::RuleNotFound),
66    }
67}
68
69/// Get the number of quadrature points for a Gauss-Jacobi rule of degree `m` on a given cell type.
70pub fn gauss_jacobi_npoints(celltype: ReferenceCellType, m: usize) -> usize {
71    let np = (m + 2) / 2;
72    match celltype {
73        ReferenceCellType::Interval => np,
74        ReferenceCellType::Quadrilateral => np.pow(2),
75        ReferenceCellType::Hexahedron => np.pow(3),
76        ReferenceCellType::Triangle => np.pow(2),
77        ReferenceCellType::Tetrahedron => np.pow(3),
78        _ => {
79            panic!("Unsupported cell type");
80        }
81    }
82}
ndelement/quadrature.rs

ndelement/
quadrature.rs
