honeycomb_kernels/skewness/
mod.rs

1//! Skewness computation routines.
2//!
3//! We use the equiangular skew formula presented
4//! [here](https://en.wikipedia.org/wiki/Types_of_mesh#Skewness)
5
6use honeycomb_core::{
7    cmap::{CMap2, CMap3, DartIdType, FaceIdType},
8    geometry::CoordsFloat,
9    stm::atomically,
10};
11
12/// Return the skewness of a given face.
13///
14/// # Arguments
15///
16/// - `map: &CMap2<T>` -- Input map.
17/// - `fid: FaceIdType` -- Face to compute the skewness of.
18///
19/// # Return
20///
21/// The value returned is comprised in the `[0.0; 1.0]` range. `0.0` corresponds to ideal
22/// (equilateral) cells while `1.0` corresponds to degenerate cells.
23///
24/// # Panics
25///
26/// This function will panic if a topological vertex has no associated coordinates.
27// #[inline] // bench and adjust
28#[must_use = "unused return value"]
29pub fn compute_face_skewness_2d<T: CoordsFloat>(map: &CMap2<T>, fid: FaceIdType) -> T {
30    let (mut d1, mut d2, mut d3) = (
31        fid as DartIdType,
32        map.beta::<1>(fid as DartIdType),
33        map.beta::<1>(map.beta::<1>(fid as DartIdType)),
34    );
35    let (mut vid1, mut vid2, mut vid3) = (map.vertex_id(d1), map.vertex_id(d2), map.vertex_id(d3));
36    let mut cnt = 0;
37    let mut min_theta = T::max_value();
38    let mut max_theta = T::min_value();
39    loop {
40        let theta = atomically(|t| {
41            let v1 = map.read_vertex(t, vid1)?.unwrap();
42            let v2 = map.read_vertex(t, vid2)?.unwrap();
43            let v3 = map.read_vertex(t, vid3)?.unwrap();
44            let vin = v1 - v2;
45            let vout = v3 - v2;
46            Ok(T::acos(vin.dot(&vout) / (vin.norm() * vout.norm())))
47        });
48        min_theta = min_theta.min(theta);
49        max_theta = max_theta.max(theta);
50        // move forward
51        cnt += 1;
52        d1 = d2;
53        d2 = d3;
54        d3 = map.beta::<1>(d3);
55        vid1 = vid2;
56        vid2 = vid3;
57        vid3 = map.vertex_id(d3);
58        if d1 == fid as DartIdType {
59            break;
60        }
61    }
62    let ideal_theta = T::from(f64::from(cnt - 2) * std::f64::consts::PI / f64::from(cnt)).unwrap();
63
64    ((max_theta - ideal_theta) / (T::from(std::f64::consts::PI).unwrap() - ideal_theta))
65        .max((ideal_theta - min_theta) / ideal_theta)
66}
67
68/// Return the skewness of a given face.
69///
70/// # Arguments
71///
72/// - `map: &CMap3<T>` -- Input map.
73/// - `fid: FaceIdType` -- Face to compute the skewness of.
74///
75/// # Return
76///
77/// The value returned is comprised in the `[0.0; 1.0]` range. `0.0` corresponds to ideal
78/// (equilateral) cells while `1.0` corresponds to degenerate cells.
79///
80/// # Panics
81///
82/// This function will panic if a topological vertex has no associated coordinates.
83// #[inline] // bench and adjust
84#[must_use = "unused return value"]
85pub fn compute_face_skewness_3d<T: CoordsFloat>(map: &CMap3<T>, fid: FaceIdType) -> T {
86    let (mut d1, mut d2, mut d3) = (
87        fid as DartIdType,
88        map.beta::<1>(fid as DartIdType),
89        map.beta::<1>(map.beta::<1>(fid as DartIdType)),
90    );
91    let (mut vid1, mut vid2, mut vid3) = (map.vertex_id(d1), map.vertex_id(d2), map.vertex_id(d3));
92    let mut cnt = 0;
93    let mut min_theta = T::max_value();
94    let mut max_theta = T::min_value();
95
96    loop {
97        let theta = atomically(|t| {
98            let v1 = map.read_vertex(t, vid1)?.unwrap();
99            let v2 = map.read_vertex(t, vid2)?.unwrap();
100            let v3 = map.read_vertex(t, vid3)?.unwrap();
101            let vin = v1 - v2;
102            let vout = v3 - v2;
103            Ok(T::acos(vin.dot(&vout) / (vin.norm() * vout.norm())))
104        });
105        min_theta = min_theta.min(theta);
106        max_theta = max_theta.max(theta);
107        // move forward
108        cnt += 1;
109        d1 = d2;
110        d2 = d3;
111        d3 = map.beta::<1>(d3);
112        vid1 = vid2;
113        vid2 = vid3;
114        vid3 = map.vertex_id(d3);
115        if d1 == fid as DartIdType {
116            break;
117        }
118    }
119    let ideal_theta = T::from(f64::from(cnt - 2) * std::f64::consts::PI / f64::from(cnt)).unwrap();
120
121    ((max_theta - ideal_theta) / (T::from(std::f64::consts::PI).unwrap() - ideal_theta))
122        .max((ideal_theta - min_theta) / ideal_theta)
123}
124
125#[cfg(test)]
126mod tests;
honeycomb_kernels/skewness/mod.rs

honeycomb_kernels/skewness/
mod.rs
