aeon-tk 0.3.3

Toolkit for running finite difference simulations with adaptive mesh refinement
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

use crate::element::{Basis, BasisFunction as _, LeastSquares, LinearOperator, Support};
use faer::{
    ColRef, Conj, Mat, MatRef,
    col::generic::Col,
    linalg::{matmul::dot, svd::SvdError},
};
use reborrow::{Reborrow, ReborrowMut};

#[derive(Clone)]
pub struct ApproxOperator {
    lls: LeastSquares,
    /// Cache for a vandermonde matrix.
    vandermonde: Mat<f64>,
    /// Cache for rhs of computation.
    rhs: Mat<f64>,
    /// Shape functions for each operation.
    shape: Mat<f64>,
}

impl ApproxOperator {
    pub fn build<const N: usize>(
        &mut self,
        support: &impl Support<N>,
        basis: &impl Basis<N>,
        ops: &impl LinearOperator<N>,
    ) -> Result<(), SvdError> {
        self.vandermonde
            .resize_with(support.num_points(), basis.order(), |i, j| {
                let point = support.point(i);
                basis.func(j).value(point)
            });

        self.shape
            .resize_with(support.num_points(), ops.num_operations(), |_, _| 0.0);

        self.rhs
            .resize_with(basis.order(), ops.num_operations(), |deg, i| {
                ops.apply(i, &basis.func(deg))
            });

        self.lls.least_squares(
            self.vandermonde.transpose(),
            self.shape.rb_mut(),
            self.rhs.rb(),
        )?;

        Ok(())
    }

    pub fn apply(&self, i: usize, values: &[f64]) -> f64 {
        dot::inner_prod(
            self.shape.row(i),
            Conj::Yes,
            Col::from_slice(values),
            Conj::No,
        )
    }

    pub fn weights(&self, i: usize) -> ColRef<'_, f64> {
        self.shape.col(i)
    }

    pub fn shape(&self) -> MatRef<'_, f64> {
        self.shape.rb()
    }
}

impl Default for ApproxOperator {
    fn default() -> Self {
        Self {
            lls: Default::default(),
            vandermonde: Mat::zeros(0, 0),
            rhs: Mat::zeros(0, 0),
            shape: Mat::zeros(0, 0),
        }
    }
}