cmtool_core/utils/
mod.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2
3use crate::coordinates::*;
4use crate::ensight_gold::types::{ElementsType, VolumeElementTypes};
5mod area;
6pub use area::*;
7mod polygon;
8
9#[inline(always)]
10pub fn linear_index_2d_matrix_row_major(i_coord: usize, i_axis: usize, n_col: usize) -> usize {
11    n_col * i_coord + i_axis
12}
13#[inline(always)]
14pub fn linear_index_coordinates_matrix(i_coord: usize, i_axis: usize) -> usize {
15    linear_index_2d_matrix_row_major(i_coord, i_axis, NUMBER_OF_AXIS)
16}
17
18pub type AxisPoints = [usize; NUMBER_OF_AXIS];
19
20/// Computes the signed volume of a tetrahedron defined by four 3D points.
21///
22/// # Arguments
23///
24/// * `a`, `b`, `c`, `d` - The vertices of the tetrahedron (3D coordinates).
25///
26/// # Returns
27///
28/// The absolute value of the mixed product divided by 6, which gives the
29/// volume of the tetrahedron.
30///
31fn tetra_volume(
32    a: CartesianCoordinates,
33    b: CartesianCoordinates,
34    c: CartesianCoordinates,
35    d: CartesianCoordinates,
36) -> f64 {
37    // let ab = b.sub(&a);
38    // let ac = c.sub(&a);
39    // let ad = d.sub(&a);
40    let ac = CartesianVec3::from_point(c, a);
41    let ab = CartesianVec3::from_point(b, a);
42    let ad = CartesianVec3::from_point(d, a);
43
44    ab.dot(&ac.cross(&ad)).abs() / 6.0
45}
46
47impl VolumeElementTypes {
48    /// Returns the list of tetrahedral subdivisions for this volume element type.
49    ///
50    /// Each 4-index array represents one tetrahedron by specifying vertex indices
51    /// within the original element. These are used to compute the total volume
52    /// of the element as a sum of sub-tetrahedron volumes.
53    ///
54    /// The subdivisions follow the same decomposition logic as the original
55    /// `c_numbering_tetra` table from the C++ implementation.
56    ///
57    /// # Returns
58    ///
59    /// A reference to a static array of vertex index groups (`[usize; 4]`),
60    /// where each group defines one tetrahedron.
61    fn tetra_subdivisions(&self) -> &'static [[usize; 4]] {
62        match self {
63            Self::Tetra4 | Self::GTetra4 => &[[0, 1, 2, 3]],
64            Self::Pyramid5 | Self::GPyramid5 => &[[0, 1, 3, 4], [1, 2, 3, 4]],
65            Self::Penta6 | Self::GPenta6 => &[[0, 1, 2, 3], [1, 2, 3, 4], [2, 3, 4, 5]],
66            Self::Hexa8 | Self::GHexa8 => &[
67                [0, 1, 3, 4],
68                [1, 3, 4, 5],
69                [3, 4, 5, 7],
70                [2, 3, 1, 6],
71                [3, 1, 6, 7],
72                [1, 6, 7, 5],
73            ],
74            Self::Tetra10 | Self::GTetra10 => &[[0, 1, 2, 3]],
75            Self::Pyramid13 | Self::GPyramid13 => &[[0, 1, 3, 4], [1, 2, 3, 4]],
76            Self::Penta15 | Self::GPenta15 => &[[0, 1, 2, 3], [1, 2, 3, 4], [2, 3, 4, 5]],
77            Self::Hexa20 | Self::GHexa20 => &[
78                [0, 1, 3, 4],
79                [1, 3, 4, 5],
80                [3, 4, 5, 7],
81                [2, 3, 1, 6],
82                [3, 1, 6, 7],
83                [1, 6, 7, 5],
84            ],
85        }
86    }
87}
88
89pub fn compute_centroid<'a, I>(vertices: I) -> CartesianCoordinates
90where
91    I: Iterator<Item = &'a [f64; 3]>,
92{
93    let mut coords: Coords3 = Default::default();
94    let mut count = 0;
95
96    // for vertex in vertices {
97    //     coords
98    //         .iter_mut()
99    //         .zip(vertex.iter())
100    //         .for_each(|(c, v)| *c += *v);
101    //     count += 1;
102    // }
103    // if count > 0 {
104    //     coords.iter_mut().for_each(|c| *c /= count as f64);
105    // }
106
107    //SIMD friendly version ? Even if it's not the case, loop unroling here is still very readable
108    for vertex in vertices {
109        coords[0] += vertex[0];
110        coords[1] += vertex[1];
111        coords[2] += vertex[2];
112        count += 1;
113    }
114
115    //Same here loop unrolling may improve SIMD and still elegant
116    if count > 0 {
117        let inv_count = 1.0 / count as f64;
118        coords[0] *= inv_count;
119        coords[1] *= inv_count;
120        coords[2] *= inv_count;
121    }
122
123    CartesianCoordinates(coords)
124}
125
126/// Computes the total volume of a given volume element by summing the
127/// volumes of its tetrahedral subdivisions.
128///
129/// # Arguments
130///
131/// * `local_vertices` - A slice of the 3D coordinates of the element's vertices.
132/// * `elem_type` - The type of the volume element (e.g., Hexa8, Penta6).
133///
134/// # Returns
135///
136/// An `Option<f64>` containing the total volume if the input is valid.
137/// Returns `None` if the number of vertices does not match the element type.
138///
139pub fn compute_volume(
140    local_vertices: &[CartesianCoordinates],
141    elem_type: VolumeElementTypes,
142) -> Option<f64> {
143    let tetra_indices = elem_type.tetra_subdivisions();
144
145    if local_vertices.len() != ElementsType::VolumeElementType(elem_type).node_count() as usize {
146        return None;
147    }
148
149    let vol = tetra_indices
150        .iter()
151        .map(|&[i0, i1, i2, i3]| {
152            let a = local_vertices[i0];
153            let b = local_vertices[i1];
154            let c = local_vertices[i2];
155            let d = local_vertices[i3];
156            tetra_volume(a, b, c, d)
157        })
158        .sum();
159
160    Some(vol)
161}
162
163#[cfg(test)]
164mod tests {
165    use super::*;
166
167    #[test]
168    fn unit_tetrahedron_volume() {
169        let a = CartesianCoordinates([0.0, 0.0, 0.0]);
170        let b = CartesianCoordinates([1.0, 0.0, 0.0]);
171        let c = CartesianCoordinates([0.0, 1.0, 0.0]);
172        let d = CartesianCoordinates([0.0, 0.0, 1.0]);
173
174        let vol = tetra_volume(a, b, c, d);
175        assert!((vol - 1.0 / 6.0).abs() < 1e-10);
176    }
177
178    #[test]
179    fn hexa8_unit_cube_volume() {
180        let vertices: Vec<CartesianCoordinates> = vec![
181            [0.0, 0.0, 0.0], // 0
182            [1.0, 0.0, 0.0], // 1
183            [1.0, 1.0, 0.0], // 2
184            [0.0, 1.0, 0.0], // 3
185            [0.0, 0.0, 1.0], // 4
186            [1.0, 0.0, 1.0], // 5
187            [1.0, 1.0, 1.0], // 6
188            [0.0, 1.0, 1.0], // 7
189        ]
190        .into_iter()
191        .map(CartesianCoordinates)
192        .collect();
193
194        let volume = compute_volume(&vertices, VolumeElementTypes::Hexa8);
195        assert!(volume.is_some());
196        let volume = volume.unwrap();
197        assert!(
198            (volume - 1.0).abs() < 1e-10,
199            "Expected volume ≈ 1.0, got {}",
200            volume
201        );
202    }
203
204    #[test]
205    fn translated_hexa8_cube_volume() {
206        let size = 2.0;
207        let offset = 3.0;
208
209        let vertices: Vec<CartesianCoordinates> = vec![
210            [offset, offset, offset],
211            [offset + size, offset, offset],
212            [offset + size, offset + size, offset],
213            [offset, offset + size, offset],
214            [offset, offset, offset + size],
215            [offset + size, offset, offset + size],
216            [offset + size, offset + size, offset + size],
217            [offset, offset + size, offset + size],
218        ]
219        .into_iter()
220        .map(CartesianCoordinates)
221        .collect();
222
223        let volume = compute_volume(&vertices, VolumeElementTypes::Hexa8);
224        assert!(volume.is_some());
225        let volume = volume.unwrap();
226        let expected = size.powi(3);
227        assert!(
228            (volume - expected).abs() < 1e-10,
229            "Expected volume ≈ {}, got {}",
230            expected,
231            volume
232        );
233    }
234
235    #[test]
236    fn test_centroid_triangle_2d() {
237        let v1 = &[0.0, 0.0, 0.0];
238        let v2 = &[1.0, 0.0, 0.0];
239        let v3 = &[0.0, 1.0, 0.0];
240
241        let centroid = compute_centroid([v1, v2, v3].iter().copied());
242
243        assert_eq!(centroid.0, [1.0 / 3.0, 1.0 / 3.0, 0.0]);
244    }
245
246    #[test]
247    fn test_centroid_hexa() {
248        //Vertices are chosen to be -1 0 or 1 to have origin as centroid
249        let cube_vertices = [
250            &[-1.0, -1.0, -1.0],
251            &[1.0, -1.0, -1.0],
252            &[1.0, 1.0, -1.0],
253            &[-1.0, 1.0, -1.0],
254            &[-1.0, -1.0, 1.0],
255            &[1.0, -1.0, 1.0],
256            &[1.0, 1.0, 1.0],
257            &[-1.0, 1.0, 1.0],
258        ];
259
260        let centroid = compute_centroid(cube_vertices.iter().copied());
261
262        assert!(
263            centroid.0.iter().all(|c| c.abs() < 1e-12),
264            "Centroid is not at origin: got {:?}",
265            centroid.0
266        );
267    }
268
269    #[test]
270    fn test_centroid_tetraheadron() {
271        //for tetra: Centroid=1/4(A+B+C+D)
272
273        let tetrahedron_vertices = [
274            &[0.0, 0.0, 0.0],
275            &[2.0, 0.0, 0.0],
276            &[0.0, 2.0, 0.0],
277            &[2.0, 2.0, 4.0],
278        ];
279
280        let centroid = compute_centroid(tetrahedron_vertices.iter().copied());
281
282        let expect_centroid = [1., 1., 1.];
283        centroid
284            .0
285            .iter()
286            .zip(expect_centroid)
287            .for_each(|(c, e)| assert!((c - e).abs() < 1e-12, "expected {e}, got {c}"));
288    }
289}
cmtool_core/utils/mod.rs

cmtool_core/utils/
mod.rs
