115 lines
3.6 KiB
Rust
115 lines
3.6 KiB
Rust
// Mesh quality and statistics calculations
|
|
|
|
use super::MeshStatistics;
|
|
use crate::mesh::entities::{Cell, Face, Node};
|
|
use nalgebra::Vector3;
|
|
|
|
/// Mesh quality analyzer
|
|
pub struct MeshQualityAnalyzer;
|
|
|
|
impl MeshQualityAnalyzer {
|
|
/// Calculate comprehensive mesh statistics
|
|
#[must_use]
|
|
pub fn analyze_mesh(
|
|
nodes: &indexmap::IndexMap<usize, Node>,
|
|
cells: &indexmap::IndexMap<usize, Cell>,
|
|
faces: &indexmap::IndexMap<usize, Face>,
|
|
) -> MeshStatistics {
|
|
let mut stats = MeshStatistics::new();
|
|
|
|
stats.total_nodes = nodes.len();
|
|
stats.total_cells = cells.len();
|
|
stats.total_faces = faces.len();
|
|
|
|
if !cells.is_empty() {
|
|
// Volume statistics
|
|
let volumes: Vec<f64> = cells
|
|
.values()
|
|
.map(super::super::traits::MeshEntity::volume)
|
|
.collect();
|
|
stats.min_cell_volume = volumes.iter().fold(f64::INFINITY, |a, &b| a.min(b));
|
|
stats.max_cell_volume = volumes.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
|
|
stats.average_cell_volume = volumes.iter().sum::<f64>() / volumes.len() as f64;
|
|
|
|
// Quality metrics
|
|
stats.max_skewness = Self::calculate_max_skewness(cells);
|
|
}
|
|
|
|
// Boundary face count
|
|
stats.boundary_faces = faces.values().filter(|face| face.is_boundary()).count();
|
|
|
|
// Aspect ratio calculation
|
|
stats.aspect_ratio = Self::calculate_aspect_ratio(nodes);
|
|
|
|
stats
|
|
}
|
|
|
|
/// Calculate maximum skewness in the mesh
|
|
fn calculate_max_skewness(cells: &indexmap::IndexMap<usize, Cell>) -> f64 {
|
|
cells
|
|
.values()
|
|
.map(super::entities::Cell::skewness)
|
|
.fold(0.0, f64::max)
|
|
}
|
|
|
|
/// Calculate mesh aspect ratio
|
|
fn calculate_aspect_ratio(nodes: &indexmap::IndexMap<usize, Node>) -> f64 {
|
|
if nodes.is_empty() {
|
|
return 1.0;
|
|
}
|
|
|
|
let mut min = Vector3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
|
|
let mut max = Vector3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
|
|
|
|
for (_, node) in nodes {
|
|
let pos = node.position();
|
|
min.x = min.x.min(pos.x);
|
|
min.y = min.y.min(pos.y);
|
|
min.z = min.z.min(pos.z);
|
|
max.x = max.x.max(pos.x);
|
|
max.y = max.y.max(pos.y);
|
|
max.z = max.z.max(pos.z);
|
|
}
|
|
|
|
let dimensions = max - min;
|
|
let max_dim = dimensions.x.max(dimensions.y).max(dimensions.z);
|
|
let min_dim = dimensions.x.min(dimensions.y).min(dimensions.z.max(1e-10));
|
|
|
|
max_dim / min_dim
|
|
}
|
|
|
|
/// Calculate mesh orthogonality
|
|
#[must_use]
|
|
pub fn calculate_orthogonality(faces: &indexmap::IndexMap<usize, Face>) -> f64 {
|
|
if faces.is_empty() {
|
|
return 1.0;
|
|
}
|
|
|
|
let mut total_orthogonality = 0.0;
|
|
let mut count = 0;
|
|
|
|
for (_, face) in faces {
|
|
if !face.is_boundary() {
|
|
// For internal faces, calculate orthogonality based on face normal
|
|
// and line connecting cell centers
|
|
// This is a simplified calculation
|
|
let orthogonality = 1.0; // Would need cell center information for real calculation
|
|
total_orthogonality += orthogonality;
|
|
count += 1;
|
|
}
|
|
}
|
|
|
|
if count > 0 {
|
|
total_orthogonality / f64::from(count)
|
|
} else {
|
|
1.0
|
|
}
|
|
}
|
|
|
|
/// Calculate mesh non-orthogonality
|
|
#[must_use]
|
|
pub fn calculate_non_orthogonality(faces: &indexmap::IndexMap<usize, Face>) -> f64 {
|
|
1.0 - Self::calculate_orthogonality(faces)
|
|
}
|
|
}
|