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
109
110
111
112
113
114

//! Assemblers that assemble the contributions to the global matrix due to a single cell

use crate::assembly::common::{AssemblerGeometry, RlstArray};
use crate::traits::{CellAssembler, KernelEvaluator, PotentialIntegrand};
use ndgrid::traits::GeometryMap;
use num::Zero;
use rlst::{
    rlst_dynamic_array2, rlst_dynamic_array3, DefaultIteratorMut, RandomAccessByRef, RawAccess,
    RawAccessMut, RlstScalar, Shape,
};

pub struct PotentialCellAssembler<
    'a,
    T: RlstScalar,
    I: PotentialIntegrand<T = T>,
    G: GeometryMap<T = T::Real>,
    K: KernelEvaluator<T = T>,
    Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
> {
    integrand: &'a I,
    kernel: &'a K,
    evaluator: G,
    table: &'a RlstArray<T, 4>,
    k: RlstArray<T, 3>,
    mapped_pts: RlstArray<T::Real, 2>,
    normals: RlstArray<T::Real, 2>,
    jacobians: RlstArray<T::Real, 2>,
    jdet: Vec<T::Real>,
    evaluation_points: &'a Array2,
    weights: &'a [T::Real],
    cell: usize,
}

impl<
        'a,
        T: RlstScalar,
        I: PotentialIntegrand<T = T>,
        G: GeometryMap<T = T::Real>,
        K: KernelEvaluator<T = T>,
        Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
    > PotentialCellAssembler<'a, T, I, G, K, Array2>
{
    #[allow(clippy::too_many_arguments)]
    pub fn new(
        npts: usize,
        nevalpts: usize,
        deriv_size: usize,
        integrand: &'a I,
        kernel: &'a K,
        evaluator: G,
        table: &'a RlstArray<T, 4>,
        evaluation_points: &'a Array2,
        weights: &'a [T::Real],
    ) -> Self {
        Self {
            integrand,
            kernel,
            evaluator,
            table,
            k: rlst_dynamic_array3!(T, [deriv_size, npts, nevalpts]),
            mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]),
            normals: rlst_dynamic_array2!(T::Real, [3, npts]),
            jacobians: rlst_dynamic_array2!(T::Real, [6, npts]),
            jdet: vec![T::Real::zero(); npts],
            evaluation_points,
            weights,
            cell: 0,
        }
    }
}
impl<
        'a,
        T: RlstScalar,
        I: PotentialIntegrand<T = T>,
        G: GeometryMap<T = T::Real>,
        K: KernelEvaluator<T = T>,
        Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
    > CellAssembler for PotentialCellAssembler<'a, T, I, G, K, Array2>
{
    type T = T;
    fn set_cell(&mut self, cell: usize) {
        self.cell = cell;
        self.evaluator.points(cell, self.mapped_pts.data_mut());
        self.evaluator.jacobians_dets_normals(
            cell,
            self.jacobians.data_mut(),
            &mut self.jdet,
            self.normals.data_mut(),
        );
    }

    fn assemble(&mut self, local_mat: &mut RlstArray<T, 2>) {
        self.kernel.assemble_st(
            self.mapped_pts.data(),
            self.evaluation_points.data(),
            self.k.data_mut(),
        );

        let geometry =
            AssemblerGeometry::new(&self.mapped_pts, &self.normals, &self.jacobians, &self.jdet);

        for (i, mut col) in local_mat.col_iter_mut().enumerate() {
            for (eval_index, entry) in col.iter_mut().enumerate() {
                *entry = T::zero();
                for (index, wt) in self.weights.iter().enumerate() {
                    *entry += self
                        .integrand
                        .evaluate(self.table, index, eval_index, i, &self.k, &geometry)
                        * num::cast::<T::Real, T>(*wt * self.jdet[index]).unwrap();
                }
            }
        }
    }
}