autofoam 0.2.5

autofoam and related tools
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

use std::error::Error;

use clap::{ArgGroup, Parser};

use autofoam::{
    histogram::weighted_histogram,
    interpolation::interpolate,
    vtk::{calculate_polygon_areas, VtpProcessor},
};

#[derive(Parser)]
#[command(
    about = "Determines the scalar field value that defines a region with a specified total area."
)]
#[command(group(
    ArgGroup::new("mode")
        .required(true)
        .args(["percentile", "area"]),
))]
pub struct Args {
    #[arg(long, help = "Path to .vtp file")]
    pub file: String,

    #[arg(long, help = "Scalar field name")]
    pub field: String,

    #[arg(long, help = "Percentile threshold (0-100)")]
    pub percentile: Option<f64>,

    #[arg(long, help = "Absolute area threshold")]
    pub area: Option<f64>,
}

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse();

    let vtp = VtpProcessor::from_file(&args.file)?;
    let (points, connectivity, offsets) = vtp.geometry()?;

    let scalar_vec = vtp.field(&args.field)?;

    let area_vec = calculate_polygon_areas(&points, &connectivity, &offsets);

    let bin_width = 0.1;

    let histogram = weighted_histogram(&scalar_vec, &area_vec, &bin_width);

    let area_binned = &histogram.heights;
    let bins = &histogram.bin_edges;

    let mut area_cumsum: Vec<f64> = area_binned
        .iter()
        .scan(0.0, |sum, &val| {
            *sum += val;
            Some(*sum)
        })
        .collect();
    if let Some(&last) = area_cumsum.last() {
        area_cumsum.push(last);
    }

    let area_target = if let Some(percentile) = args.percentile {
        (percentile * 0.01) * area_vec.iter().sum::<f64>()
    } else if let Some(area) = args.area {
        area
    } else {
        panic!("Either percentile or area must be provided");
    };

    let scalar_target = interpolate(&area_cumsum, bins, area_target);

    println!("{:.1}", scalar_target);

    Ok(())
}