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

//! Hypersingular integrand
use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
use rlst::RlstScalar;

/// Integrand for the curl curl term of a hypersingular boundary operator
pub struct HypersingularCurlCurlBoundaryIntegrand<T: RlstScalar> {
    _t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> HypersingularCurlCurlBoundaryIntegrand<T> {
    /// Create new
    pub fn new() -> Self {
        Self {
            _t: std::marker::PhantomData,
        }
    }
}

impl<T: RlstScalar> Default for HypersingularCurlCurlBoundaryIntegrand<T> {
    fn default() -> Self {
        Self::new()
    }
}

unsafe impl<T: RlstScalar> BoundaryIntegrand for HypersingularCurlCurlBoundaryIntegrand<T> {
    type T = T;

    fn evaluate(
        &self,
        k: &impl Access1D<T = T>,
        test_table: &impl Access2D<T = T>,
        trial_table: &impl Access2D<T = T>,
        test_geometry: &impl GeometryAccess<T = T>,
        trial_geometry: &impl GeometryAccess<T = T>,
    ) -> T {
        unsafe {
            let test0 = test_table.get(1, 0);
            let test1 = test_table.get(2, 0);
            let trial0 = trial_table.get(1, 0);
            let trial1 = trial_table.get(2, 0);

            k.get(0)
                * ((test_geometry.jacobian(3) * test0 - test_geometry.jacobian(0) * test1)
                    * (trial_geometry.jacobian(3) * trial0 - trial_geometry.jacobian(0) * trial1)
                    + (test_geometry.jacobian(4) * test0 - test_geometry.jacobian(1) * test1)
                        * (trial_geometry.jacobian(4) * trial0
                            - trial_geometry.jacobian(1) * trial1)
                    + (test_geometry.jacobian(5) * test0 - test_geometry.jacobian(2) * test1)
                        * (trial_geometry.jacobian(5) * trial0
                            - trial_geometry.jacobian(2) * trial1))
                / test_geometry.jdet()
                / trial_geometry.jdet()
        }
    }
}

/// Integrand for the normal normal term of a hypersingular boundary operator
pub struct HypersingularNormalNormalBoundaryIntegrand<T: RlstScalar> {
    _t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> HypersingularNormalNormalBoundaryIntegrand<T> {
    /// Create new
    pub fn new() -> Self {
        Self {
            _t: std::marker::PhantomData,
        }
    }
}

impl<T: RlstScalar> Default for HypersingularNormalNormalBoundaryIntegrand<T> {
    fn default() -> Self {
        Self::new()
    }
}

unsafe impl<T: RlstScalar> BoundaryIntegrand for HypersingularNormalNormalBoundaryIntegrand<T> {
    type T = T;

    fn evaluate(
        &self,
        k: &impl Access1D<T = T>,
        test_table: &impl Access2D<T = T>,
        trial_table: &impl Access2D<T = T>,
        test_geometry: &impl GeometryAccess<T = T>,
        trial_geometry: &impl GeometryAccess<T = T>,
    ) -> T {
        unsafe {
            k.get(0)
                * (trial_geometry.normal(0) * test_geometry.normal(0)
                    + trial_geometry.normal(1) * test_geometry.normal(1)
                    + trial_geometry.normal(2) * test_geometry.normal(2))
                * test_table.get(0, 0)
                * trial_table.get(0, 0)
        }
    }
}