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

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

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