fenris 0.0.15

A library for advanced finite element computations in Rust
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
128
129
130
131
132
133
134
135
136
137
138
139

use std::ops::{Add, AddAssign, Deref, Mul};

use nalgebra::allocator::Allocator;
use nalgebra::{DefaultAllocator, DimName, OPoint, Point1, Scalar, U2, U3};
use num::Zero;

/// Errors returned by quadrature methods.
///
/// TODO: How to prevent collapse?
pub use fenris_quadrature::Error as QuadratureError;

use crate::nalgebra::{convert, Point2, Point3, RealField, U1};

pub mod tensor;
pub mod total_order;
pub mod univariate;

pub type QuadraturePair<T, D> = (Vec<T>, Vec<OPoint<T, D>>);
pub type QuadraturePair1d<T> = QuadraturePair<T, U1>;
pub type QuadraturePair2d<T> = QuadraturePair<T, U2>;
pub type QuadraturePair3d<T> = QuadraturePair<T, U3>;

pub trait Quadrature<T, D>
where
    T: Scalar,
    D: DimName,
    DefaultAllocator: Allocator<T, D>,
{
    type Data;

    fn weights(&self) -> &[T];
    fn points(&self) -> &[OPoint<T, D>];
    fn data(&self) -> &[Self::Data];

    /// Approximates the integral of the given function using this quadrature rule.
    fn integrate<U, Function>(&self, f: Function) -> U
    where
        Function: Fn(&OPoint<T, D>) -> U,
        U: Zero + Mul<T, Output = U> + Add<T, Output = U> + AddAssign<U>,
    {
        let mut integral = U::zero();
        for (w, p) in self.weights().iter().zip(self.points()) {
            integral += f(p) * w.clone();
        }
        integral
    }
}

/// Helper trait for 2D quadratures.
pub trait Quadrature2d<T>: Quadrature<T, U2>
where
    T: Scalar,
{
}

impl<T, X> Quadrature2d<T> for X
where
    T: Scalar,
    X: Quadrature<T, U2>,
{
}

impl<T, D, A, B> Quadrature<T, D> for (A, B)
where
    T: Scalar,
    D: DimName,
    A: Deref<Target = [T]>,
    B: Deref<Target = [OPoint<T, D>]>,
    DefaultAllocator: Allocator<T, D>,
{
    type Data = ();

    fn weights(&self) -> &[T] {
        self.0.deref()
    }

    fn points(&self) -> &[OPoint<T, D>] {
        self.1.deref()
    }

    fn data(&self) -> &[()] {
        // This may look absurd, but since we're just returning a slice to a zero-sized type (the unit type),
        // the (global) allocator never allocates anything and most likely the whole thing will get completely
        // optimized out
        vec![(); self.0.len()].leak()
    }
}

impl<T, D, X> Quadrature<T, D> for &X
where
    T: Scalar,
    D: DimName,
    X: Quadrature<T, D>,
    DefaultAllocator: Allocator<T, D>,
{
    type Data = X::Data;

    fn weights(&self) -> &[T] {
        X::weights(self)
    }

    fn points(&self) -> &[OPoint<T, D>] {
        X::points(self)
    }

    fn data(&self) -> &[Self::Data] {
        X::data(self)
    }
}

fn convert_quadrature_rule_from_1d_f64<T>(quadrature: fenris_quadrature::Rule<1>) -> QuadraturePair1d<T>
where
    T: RealField,
{
    let (weights, points) = quadrature;
    let weights = weights.into_iter().map(convert).collect();
    let points = points.into_iter().map(Point1::from).map(convert).collect();
    (weights, points)
}

fn convert_quadrature_rule_from_2d_f64<T>(quadrature: fenris_quadrature::Rule<2>) -> QuadraturePair2d<T>
where
    T: RealField,
{
    let (weights, points) = quadrature;
    let weights = weights.into_iter().map(convert).collect();
    let points = points.into_iter().map(Point2::from).map(convert).collect();
    (weights, points)
}

fn convert_quadrature_rule_from_3d_f64<T>(quadrature: fenris_quadrature::Rule<3>) -> QuadraturePair3d<T>
where
    T: RealField,
{
    let (weights, points) = quadrature;
    let weights = weights.into_iter().map(convert).collect();
    let points = points.into_iter().map(Point3::from).map(convert).collect();
    (weights, points)
}