use crate::data::coordinates::Coordinates;
use crate::data::section::Section;
use crate::data::storage::Storage;
use crate::geometry::cell_vertex_ordering::canonical_cell_vertices;
use crate::mesh_error::MeshSieveError;
use crate::topology::cell_type::CellType;
use crate::topology::point::PointId;
use crate::topology::sieve::{OrientedSieve, Sieve};
use std::cmp::Ordering;
const EPS: f64 = 1e-12;
fn lex_coordinate_key_cmp(a: &[f64], b: &[f64]) -> Ordering {
a.iter()
.zip(b)
.map(|(x, y)| x.total_cmp(y))
.find(|ordering| *ordering != Ordering::Equal)
.unwrap_or_else(|| a.len().cmp(&b.len()))
}
#[derive(Clone, Debug)]
pub struct FvmFaceMetrics {
pub face: PointId,
pub owner: PointId,
pub neighbor: Option<PointId>,
pub centroid: [f64; 3],
pub area_vector: [f64; 3],
pub area_magnitude: f64,
pub owner_to_neighbor: Option<[f64; 3]>,
pub orthogonal_distance: f64,
pub non_orthogonality_deg: Option<f64>,
pub skewness_vector: [f64; 3],
pub owner_to_face: [f64; 3],
pub neighbor_to_face: Option<[f64; 3]>,
}
pub fn build_fvm_face_metrics<S, CT, CS>(
sieve: &S,
cell_types: &Section<CellType, CT>,
coordinates: &Coordinates<f64, CS>,
) -> Result<Vec<FvmFaceMetrics>, MeshSieveError>
where
S: OrientedSieve<Point = PointId>,
CT: Storage<CellType>,
CS: Storage<f64>,
{
let faces: Vec<PointId> = sieve
.points()
.filter(|p| is_face(*p, sieve, cell_types).unwrap_or(false))
.collect();
let mut keyed = Vec::with_capacity(faces.len());
for face in faces {
let vertices = face_vertices(sieve, face, coordinates)?;
let key: Vec<f64> = vertices.iter().flat_map(|p| p.iter().copied()).collect();
keyed.push((key, face));
}
keyed.sort_by(|(a, _), (b, _)| lex_coordinate_key_cmp(a, b));
keyed
.into_iter()
.map(|(_, f)| f)
.map(|f| build_face_metric(sieve, cell_types, coordinates, f))
.collect()
}
pub fn build_face_metric<S, CT, CS>(
sieve: &S,
cell_types: &Section<CellType, CT>,
coordinates: &Coordinates<f64, CS>,
face: PointId,
) -> Result<FvmFaceMetrics, MeshSieveError>
where
S: OrientedSieve<Point = PointId>,
CT: Storage<CellType>,
CS: Storage<f64>,
{
let supports: Vec<PointId> = sieve
.support_points(face)
.filter(|p| is_cell(*p, cell_types).unwrap_or(false))
.collect();
if supports.is_empty() {
return Err(MeshSieveError::InvalidGeometry(format!(
"face {face:?} has no cell support"
)));
}
let mut keyed_supports = Vec::with_capacity(supports.len());
for support in supports {
let cell_type = cell_types
.try_restrict(support)
.ok()
.and_then(|v| v.first().copied())
.ok_or_else(|| MeshSieveError::MalformedCellTopology {
cell: support,
cell_type: CellType::Polyhedron,
reason: "missing supported cell type".into(),
})?;
let centroid = cell_centroid(sieve, support, cell_type, coordinates)?;
keyed_supports.push((centroid, support));
}
keyed_supports.sort_by(|(a, _), (b, _)| lex_coordinate_key_cmp(a, b));
let owner = keyed_supports[0].1;
let neighbor = keyed_supports.get(1).map(|(_, point)| *point);
let vertices = face_vertices(sieve, face, coordinates)?;
if vertices.len() < 2 {
return Err(MeshSieveError::InvalidGeometry(
"face with <2 vertices".into(),
));
}
let centroid = polygon_centroid(&vertices);
let owner_type = cell_types
.try_restrict(owner)
.ok()
.and_then(|v| v.first().copied())
.ok_or_else(|| MeshSieveError::MalformedCellTopology {
cell: owner,
cell_type: CellType::Polyhedron,
reason: "missing owner cell type".into(),
})?;
let owner_centroid = cell_centroid(sieve, owner, owner_type, coordinates)?;
let mut area_vec = polygon_area_vector(&vertices);
if let Some(n) = neighbor {
let n_type = cell_types
.try_restrict(n)
.ok()
.and_then(|v| v.first().copied())
.ok_or_else(|| MeshSieveError::MalformedCellTopology {
cell: n,
cell_type: CellType::Polyhedron,
reason: "missing neighbor cell type".into(),
})?;
let nc = cell_centroid(sieve, n, n_type, coordinates)?;
if dot(area_vec, sub(nc, owner_centroid)) < 0.0 {
area_vec = scale(area_vec, -1.0);
}
} else if dot(area_vec, sub(centroid, owner_centroid)) < 0.0 {
area_vec = scale(area_vec, -1.0);
}
let area_mag = norm(area_vec);
let owner_to_face = sub(centroid, owner_centroid);
let (owner_to_neighbor, orth_dist, non_ortho, skew, neigh_to_face) = if let Some(n) = neighbor {
let n_type = cell_types
.try_restrict(n)
.ok()
.and_then(|v| v.first().copied())
.ok_or_else(|| MeshSieveError::MalformedCellTopology {
cell: n,
cell_type: CellType::Polyhedron,
reason: "missing neighbor cell type".into(),
})?;
let nc = cell_centroid(sieve, n, n_type, coordinates)?;
let d = sub(nc, owner_centroid);
let dmag = norm(d);
let ahat = if area_mag > EPS {
scale(area_vec, 1.0 / area_mag)
} else {
[0.0; 3]
};
let orth = dot(d, ahat).abs();
let angle = if dmag > EPS && area_mag > EPS {
let c = (dot(area_vec, d).abs() / (area_mag * dmag)).clamp(-1.0, 1.0);
Some(c.acos().to_degrees())
} else {
None
};
let lambda = if dmag > EPS {
dot(sub(centroid, owner_centroid), d) / (dmag * dmag)
} else {
0.5
};
let closest_on_d = add(owner_centroid, scale(d, lambda));
let skew = sub(centroid, closest_on_d);
(Some(d), orth, angle, skew, Some(sub(centroid, nc)))
} else {
let ahat = if area_mag > EPS {
scale(area_vec, 1.0 / area_mag)
} else {
[0.0; 3]
};
let orth = dot(owner_to_face, ahat).abs();
let skew = sub(
centroid,
add(owner_centroid, scale(ahat, dot(owner_to_face, ahat))),
);
(None, orth, None, skew, None)
};
Ok(FvmFaceMetrics {
face,
owner,
neighbor,
centroid,
area_vector: area_vec,
area_magnitude: area_mag,
owner_to_neighbor,
orthogonal_distance: orth_dist,
non_orthogonality_deg: non_ortho,
skewness_vector: skew,
owner_to_face,
neighbor_to_face: neigh_to_face,
})
}
fn is_cell<C: Storage<CellType>>(
p: PointId,
cell_types: &Section<CellType, C>,
) -> Result<bool, MeshSieveError> {
Ok(cell_types
.try_restrict(p)
.ok()
.and_then(|v| v.first().copied())
.map(|ct| ct.dimension() >= 2)
.unwrap_or(false))
}
fn is_face<S, CT>(
p: PointId,
sieve: &S,
cell_types: &Section<CellType, CT>,
) -> Result<bool, MeshSieveError>
where
S: Sieve<Point = PointId>,
CT: Storage<CellType>,
{
let support_cells = sieve
.support_points(p)
.filter(|q| is_cell(*q, cell_types).unwrap_or(false))
.count();
if support_cells == 0 || support_cells > 2 {
return Ok(false);
}
let cone: Vec<_> = sieve.cone_points(p).collect();
if cone.len() < 2 {
return Ok(false);
}
Ok(cone
.into_iter()
.all(|q| sieve.cone_points(q).next().is_none()))
}
fn face_vertices<S, C: Storage<f64>>(
sieve: &S,
face: PointId,
coordinates: &Coordinates<f64, C>,
) -> Result<Vec<[f64; 3]>, MeshSieveError>
where
S: OrientedSieve<Point = PointId>,
{
let points: Vec<_> = sieve.cone_points(face).collect();
let kind = match points.len() {
2 => {
return points
.into_iter()
.map(|p| coord3(coordinates.try_restrict(p)?))
.collect();
}
3 => CellType::Triangle,
4 => CellType::Quadrilateral,
_ => {
return Err(MeshSieveError::MalformedCellTopology {
cell: face,
cell_type: CellType::Polyhedron,
reason: format!("face has {} vertices", points.len()),
});
}
};
Ok(canonical_cell_vertices(sieve, face, kind, coordinates)?.coordinates)
}
fn cell_centroid<S, C: Storage<f64>>(
sieve: &S,
cell: PointId,
cell_type: CellType,
coordinates: &Coordinates<f64, C>,
) -> Result<[f64; 3], MeshSieveError>
where
S: OrientedSieve<Point = PointId>,
{
match canonical_cell_vertices(sieve, cell, cell_type, coordinates) {
Ok(canonical) => Ok(mean(&canonical.coordinates)),
Err(MeshSieveError::UnsupportedCanonicalCellType { .. }) => {
let mut verts = Vec::new();
for (p, _) in sieve.closure_o([cell]) {
if p != cell && sieve.cone_o(p).next().is_none() && !verts.contains(&p) {
verts.push(p);
}
}
let coords: Result<Vec<_>, _> = verts
.into_iter()
.map(|v| coord3(coordinates.try_restrict(v)?))
.collect();
Ok(mean(&coords?))
}
Err(error) => Err(error),
}
}
fn coord3(values: &[f64]) -> Result<[f64; 3], MeshSieveError> {
match values.len().cmp(&3) {
Ordering::Equal => Ok([values[0], values[1], values[2]]),
Ordering::Less if values.len() == 2 => Ok([values[0], values[1], 0.0]),
_ => Err(MeshSieveError::InvalidGeometry(
"expected 2D/3D coordinates".into(),
)),
}
}
fn polygon_area_vector(vertices: &[[f64; 3]]) -> [f64; 3] {
if vertices.len() == 2 {
return sub(vertices[1], vertices[0]);
}
let c = mean(vertices);
let mut sum = [0.0; 3];
for i in 0..vertices.len() {
let a = sub(vertices[i], c);
let b = sub(vertices[(i + 1) % vertices.len()], c);
sum = add(sum, scale(cross(a, b), 0.5));
}
sum
}
fn polygon_centroid(vertices: &[[f64; 3]]) -> [f64; 3] {
if vertices.len() <= 2 {
return mean(vertices);
}
let c0 = mean(vertices);
let mut weighted = [0.0; 3];
let mut wsum = 0.0;
for i in 0..vertices.len() {
let a = vertices[i];
let b = vertices[(i + 1) % vertices.len()];
let tri_n = cross(sub(a, c0), sub(b, c0));
let w = 0.5 * norm(tri_n);
if w > EPS {
let tri_c = scale(add(add(c0, a), b), 1.0 / 3.0);
weighted = add(weighted, scale(tri_c, w));
wsum += w;
}
}
if wsum > EPS {
scale(weighted, 1.0 / wsum)
} else {
c0
}
}
fn mean(v: &[[f64; 3]]) -> [f64; 3] {
let mut c = [0.0; 3];
for p in v {
c = add(c, *p);
}
scale(c, 1.0 / (v.len() as f64))
}
fn add(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
[a[0] + b[0], a[1] + b[1], a[2] + b[2]]
}
fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
[a[0] - b[0], a[1] - b[1], a[2] - b[2]]
}
fn scale(a: [f64; 3], s: f64) -> [f64; 3] {
[a[0] * s, a[1] * s, a[2] * s]
}
fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}
fn norm(a: [f64; 3]) -> f64 {
dot(a, a).sqrt()
}
fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
[
a[1] * b[2] - a[2] * b[1],
a[2] * b[0] - a[0] * b[2],
a[0] * b[1] - a[1] * b[0],
]
}