bempp_octree/
geometry.rs

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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
//! Geometry information

use mpi::traits::Equivalence;

use crate::constants::DEEPEST_LEVEL;

/// Definition of a point in 3d space.
///
/// A point consists of three coordinates and a global id.
/// The global id makes it easier to identify points as they
/// are distributed across MPI nodes.
#[derive(Clone, Copy, Equivalence)]
pub struct Point {
    coords: [f64; 3],
    global_id: usize,
}

impl Point {
    /// Create a new point from coordinates and global id.
    pub fn new(coords: [f64; 3], global_id: usize) -> Self {
        Self { coords, global_id }
    }

    /// Return the coordintes of a point.
    pub fn coords(&self) -> [f64; 3] {
        self.coords
    }

    /// Return a mutable pointer to the coordinates.
    pub fn coords_mut(&mut self) -> &mut [f64; 3] {
        &mut self.coords
    }

    /// Return the global id of the point.
    pub fn global_id(&self) -> usize {
        self.global_id
    }
}

/// A bounding box describes geometry in which an Octree lives.
pub struct PhysicalBox {
    coords: [f64; 6],
}

impl PhysicalBox {
    /// Create a new bounding box.
    ///
    /// The coordinates are given by `[xmin, ymin, zmin, xmax, ymax, zmax]`.
    pub fn new(coords: [f64; 6]) -> Self {
        Self { coords }
    }

    /// Give a slice of points. Compute an associated bounding box.
    ///
    /// The created bounding box is slightly bigger than the actual point set.
    /// This is to make sure that no point lies exactly on the boundary of the box.
    pub fn from_points(points: &[Point]) -> PhysicalBox {
        let mut xmin = f64::MAX;
        let mut xmax = f64::MIN;

        let mut ymin = f64::MAX;
        let mut ymax = f64::MIN;

        let mut zmin = f64::MAX;
        let mut zmax = f64::MIN;

        for point in points {
            let x = point.coords()[0];
            let y = point.coords()[1];
            let z = point.coords()[2];

            xmin = f64::min(xmin, x);
            xmax = f64::max(xmax, x);

            ymin = f64::min(ymin, y);
            ymax = f64::max(ymax, y);

            zmin = f64::min(zmin, z);
            zmax = f64::max(zmax, z);
        }

        // We want the bounding box to be slightly bigger
        // than the actual point set to avoid issues
        // at the edge of the bounding box.

        let xdiam = xmax - xmin;
        let ydiam = ymax - ymin;
        let zdiam = zmax - zmin;

        let xmean = xmin + 0.5 * xdiam;
        let ymean = ymin + 0.5 * ydiam;
        let zmean = zmin + 0.5 * zdiam;

        // We increase diameters by box size on deepest level
        // and use the maximum diameter to compute a
        // cubic bounding box.

        let deepest_box_diam = 1.0 / (1 << DEEPEST_LEVEL) as f64;

        let max_diam = [xdiam, ydiam, zdiam].into_iter().reduce(f64::max).unwrap();

        let max_diam = max_diam * (1.0 + deepest_box_diam);

        PhysicalBox {
            coords: [
                xmean - 0.5 * max_diam,
                ymean - 0.5 * max_diam,
                zmean - 0.5 * max_diam,
                xmean + 0.5 * max_diam,
                ymean + 0.5 * max_diam,
                zmean + 0.5 * max_diam,
            ],
        }
    }

    /// Return coordinates
    pub fn coordinates(&self) -> [f64; 6] {
        self.coords
    }

    /// Map a point from the reference box [0, 1]^3 to the bounding box.
    pub fn reference_to_physical(&self, point: [f64; 3]) -> [f64; 3] {
        let [xmin, ymin, zmin, xmax, ymax, zmax] = self.coords;

        [
            xmin + (xmax - xmin) * point[0],
            ymin + (ymax - ymin) * point[1],
            zmin + (zmax - zmin) * point[2],
        ]
    }

    /// Map a point from the physical domain to the reference box.
    pub fn physical_to_reference(&self, point: [f64; 3]) -> [f64; 3] {
        let [xmin, ymin, zmin, xmax, ymax, zmax] = self.coords;

        [
            (point[0] - xmin) / (xmax - xmin),
            (point[1] - ymin) / (ymax - ymin),
            (point[2] - zmin) / (zmax - zmin),
        ]
    }

    /// Return an ordered list of corners of the box.
    ///
    /// The ordering of the corners on the unit cube is
    /// [0, 0, 0]
    /// [1, 0, 0]
    /// [1, 1, 0]
    /// [0, 1, 0]
    /// [0, 0, 1]
    /// [1, 0, 1]
    /// [1, 1, 1]
    /// [0, 1, 1]
    pub fn corners(&self) -> [[f64; 3]; 8] {
        let reference_points = [
            [0.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            [1.0, 1.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, 0.0, 1.0],
            [1.0, 0.0, 1.0],
            [1.0, 1.0, 1.0],
            [0.0, 1.0, 1.0],
        ];

        [
            self.reference_to_physical(reference_points[0]),
            self.reference_to_physical(reference_points[1]),
            self.reference_to_physical(reference_points[2]),
            self.reference_to_physical(reference_points[3]),
            self.reference_to_physical(reference_points[4]),
            self.reference_to_physical(reference_points[5]),
            self.reference_to_physical(reference_points[6]),
            self.reference_to_physical(reference_points[7]),
        ]
    }
}

impl std::fmt::Display for PhysicalBox {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        let [xmin, ymin, zmin, xmax, ymax, zmax] = self.coords;

        write!(
            f,
            "(xmin: {}, ymin: {}, zmin: {}, xmax: {}, ymax: {}, zmax: {})",
            xmin, ymin, zmin, xmax, ymax, zmax
        )
    }
}