conspire/geometry/grid/from/mesh/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::geometry::{grid::Voxels, mesh::Mesh};
5use std::array::from_fn;
6
7const TETS_4: [[usize; 4]; 1] = [[0, 1, 2, 3]];
8const TETS_5: [[usize; 4]; 2] = [[0, 1, 2, 4], [0, 2, 3, 4]];
9const TETS_6: [[usize; 4]; 3] = [[0, 1, 2, 3], [1, 4, 5, 3], [1, 2, 5, 3]];
10const TETS_8: [[usize; 4]; 6] = [
11    [0, 1, 2, 6],
12    [0, 2, 3, 6],
13    [0, 3, 7, 6],
14    [0, 7, 4, 6],
15    [0, 4, 5, 6],
16    [0, 5, 1, 6],
17];
18
19impl Voxels<usize> {
20    pub fn from_finite_elements(mesh: &Mesh<3>, size: f64) -> Self {
21        let coordinates = mesh.coordinates();
22        let mut min = [f64::INFINITY; 3];
23        let mut max = [f64::NEG_INFINITY; 3];
24        for point in coordinates {
25            (0..3).for_each(|ax| {
26                min[ax] = min[ax].min(point[ax]);
27                max[ax] = max[ax].max(point[ax]);
28            });
29        }
30        let nel: [usize; 3] = from_fn(|ax| (((max[ax] - min[ax]) / size).ceil() as usize).max(1));
31        let [nx, ny, nz] = nel;
32        let mut data = vec![0usize; nx * ny * nz];
33        let numbers = mesh.blocks();
34        for (block, connectivity) in mesh.iter().enumerate() {
35            let material = numbers.map_or(block + 1, |numbers| numbers[block]);
36            for nodes in connectivity {
37                let tets: &[[usize; 4]] = match nodes.len() {
38                    4 => &TETS_4,
39                    5 => &TETS_5,
40                    6 => &TETS_6,
41                    8 => &TETS_8,
42                    _ => continue,
43                };
44                let points: Vec<[f64; 3]> = nodes
45                    .iter()
46                    .map(|&node| from_fn(|ax| coordinates[node][ax]))
47                    .collect();
48                let mut lo = [usize::MAX; 3];
49                let mut hi = [0usize; 3];
50                for point in &points {
51                    (0..3).for_each(|ax| {
52                        let l = ((point[ax] - min[ax]) / size).floor().max(0.0) as usize;
53                        let h = (((point[ax] - min[ax]) / size).floor() as usize + 1).min(nel[ax]);
54                        lo[ax] = lo[ax].min(l);
55                        hi[ax] = hi[ax].max(h);
56                    });
57                }
58                for k in lo[2]..hi[2] {
59                    for j in lo[1]..hi[1] {
60                        for i in lo[0]..hi[0] {
61                            let center =
62                                from_fn(|ax| min[ax] + ([i, j, k][ax] as f64 + 0.5) * size);
63                            if tets
64                                .iter()
65                                .any(|tet| inside(center, from_fn(|vertex| points[tet[vertex]])))
66                            {
67                                data[i + nx * j + nx * ny * k] = material;
68                            }
69                        }
70                    }
71                }
72            }
73        }
74        Voxels::new(data, nel)
75    }
76}
77
78fn inside(query: [f64; 3], tet: [[f64; 3]; 4]) -> bool {
79    let volume = orient(tet[0], tet[1], tet[2], tet[3]);
80    if volume == 0.0 {
81        return false;
82    }
83    let sign = volume.signum();
84    let tolerance = -1e-9 * volume.abs();
85    sign * orient(query, tet[1], tet[2], tet[3]) >= tolerance
86        && sign * orient(tet[0], query, tet[2], tet[3]) >= tolerance
87        && sign * orient(tet[0], tet[1], query, tet[3]) >= tolerance
88        && sign * orient(tet[0], tet[1], tet[2], query) >= tolerance
89}
90
91fn orient(a: [f64; 3], b: [f64; 3], c: [f64; 3], d: [f64; 3]) -> f64 {
92    let u = from_fn::<_, 3, _>(|ax| b[ax] - a[ax]);
93    let v = from_fn::<_, 3, _>(|ax| c[ax] - a[ax]);
94    let w = from_fn::<_, 3, _>(|ax| d[ax] - a[ax]);
95    u[0] * (v[1] * w[2] - v[2] * w[1]) - u[1] * (v[0] * w[2] - v[2] * w[0])
96        + u[2] * (v[0] * w[1] - v[1] * w[0])
97}
conspire/geometry/grid/from/mesh/mod.rs

conspire/geometry/grid/from/mesh/
mod.rs
