ndelement 0.4.0

Finite element definition 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

//! Quadrature.
use crate::types::ReferenceCellType;
use bempp_quadrature::{
    gauss_jacobi, simplex_rules,
    types::{NumericalQuadratureDefinition, QuadratureError},
};

/// Return a simplex rule for a given number of points.
///
/// If the rule does not exist `Err(())` is returned.
pub fn simplex_rule(
    cell: ReferenceCellType,
    npoints: usize,
) -> Result<NumericalQuadratureDefinition, QuadratureError> {
    match cell {
        ReferenceCellType::Point => {
            if npoints == 1 {
                Ok(NumericalQuadratureDefinition {
                    dim: 0,
                    order: 1,
                    npoints,
                    weights: vec![1.0],
                    points: vec![],
                })
            } else {
                Err(QuadratureError::RuleNotFound)
            }
        }
        ReferenceCellType::Interval => simplex_rules::simplex_rule_interval(npoints),
        ReferenceCellType::Triangle => simplex_rules::simplex_rule_triangle(npoints),
        ReferenceCellType::Quadrilateral => simplex_rules::simplex_rule_quadrilateral(npoints),
        ReferenceCellType::Tetrahedron => simplex_rules::simplex_rule_tetrahedron(npoints),
        ReferenceCellType::Hexahedron => simplex_rules::simplex_rule_hexahedron(npoints),
        ReferenceCellType::Prism => simplex_rules::simplex_rule_prism(npoints),
        ReferenceCellType::Pyramid => simplex_rules::simplex_rule_pyramid(npoints),
    }
}

/// For a given cell type return a vector with the numbers of points for which simplex rules are available.
pub fn available_simplex_rules(cell: ReferenceCellType) -> Vec<usize> {
    match cell {
        ReferenceCellType::Point => vec![1],
        ReferenceCellType::Interval => simplex_rules::available_rules_interval(),
        ReferenceCellType::Triangle => simplex_rules::available_rules_triangle(),
        ReferenceCellType::Quadrilateral => simplex_rules::available_rules_quadrilateral(),
        ReferenceCellType::Tetrahedron => simplex_rules::available_rules_tetrahedron(),
        ReferenceCellType::Hexahedron => simplex_rules::available_rules_hexahedron(),
        ReferenceCellType::Prism => simplex_rules::available_rules_prism(),
        ReferenceCellType::Pyramid => simplex_rules::available_rules_pyramid(),
    }
}

/// Get the points and weights for a Gauss-Jacobi quadrature rule of degree `m` on the given cell type.
pub fn gauss_jacobi_rule(
    celltype: ReferenceCellType,
    m: usize,
) -> Result<NumericalQuadratureDefinition, QuadratureError> {
    let np = (m + 2) / 2;
    match celltype {
        ReferenceCellType::Interval => Ok(gauss_jacobi::gauss_jacobi_interval(np)),
        ReferenceCellType::Triangle => Ok(gauss_jacobi::gauss_jacobi_triangle(np)),
        ReferenceCellType::Quadrilateral => Ok(gauss_jacobi::gauss_jacobi_quadrilateral(np)),
        ReferenceCellType::Tetrahedron => Ok(gauss_jacobi::gauss_jacobi_tetrahedron(np)),
        ReferenceCellType::Hexahedron => Ok(gauss_jacobi::gauss_jacobi_hexahedron(np)),
        _ => Err(QuadratureError::RuleNotFound),
    }
}

/// Get the number of quadrature points for a Gauss-Jacobi rule of degree `m` on a given cell type.
pub fn gauss_jacobi_npoints(celltype: ReferenceCellType, m: usize) -> usize {
    let np = (m + 2) / 2;
    match celltype {
        ReferenceCellType::Interval => np,
        ReferenceCellType::Quadrilateral => np.pow(2),
        ReferenceCellType::Hexahedron => np.pow(3),
        ReferenceCellType::Triangle => np.pow(2),
        ReferenceCellType::Tetrahedron => np.pow(3),
        _ => {
            panic!("Unsupported cell type");
        }
    }
}