use std::fmt::Debug;
use std::fs::File;
use std::io::{BufRead, BufReader, BufWriter, Write};
use std::str::FromStr;
use structopt::StructOpt;
#[derive(StructOpt, Debug)]
#[structopt(name = "dose2gmsh", author = "Max Orok <morok@mevex.com>", about = "Convert dosxyznrc 3ddose files to Gmsh msh files")]
pub struct Cli {
#[structopt(parse(from_os_str), short, long)]
pub input_file: std::path::PathBuf,
#[structopt(parse(from_os_str), short, long)]
pub output_file: Option<std::path::PathBuf>,
}
#[derive(Debug, Clone)]
pub struct DoseBlock {
pub xs: Vec<f64>,
pub ys: Vec<f64>,
pub zs: Vec<f64>,
pub doses: Vec<f64>,
pub uncerts: Vec<f64>,
}
impl DoseBlock {
pub fn from_3d_dose(input_file: &std::path::PathBuf) -> Result<DoseBlock, std::io::Error> {
let dose_input = BufReader::new(File::open(input_file)?);
let mut lines = dose_input.lines().map(|l| l.unwrap());
let (num_x, num_y, num_z) = {
let voxel_nums = lines.next().expect("voxel numbers");
let voxel_nums = parse_simple_line::<usize>(voxel_nums, "voxel number", 3);
(voxel_nums[0], voxel_nums[1], voxel_nums[2])
};
let xs = parse_simple_line::<f64>(
lines.next().expect("x-coordinates"),
"x-coordinate",
num_x + 1,
);
let ys = parse_simple_line::<f64>(
lines.next().expect("y-coordinates"),
"y-coordinate",
num_y + 1,
);
let zs = parse_simple_line::<f64>(
lines.next().expect("z-coordinates"),
"z-coordinate",
num_z + 1,
);
let num_voxels = num_x * num_y * num_z;
let doses = parse_simple_line::<f64>(lines.next().expect("doses"), "dose value", num_voxels);
let uncerts = parse_simple_line::<f64>(
lines.next().expect("uncerts"),
"uncertainty value",
num_voxels,
);
Ok(DoseBlock {
xs,
ys,
zs,
doses,
uncerts,
})
}
pub fn num_x(&self) -> usize {
self.xs.len() - 1
}
pub fn num_y(&self) -> usize {
self.ys.len() - 1
}
pub fn num_z(&self) -> usize {
self.zs.len() - 1
}
pub fn num_voxels(&self) -> usize {
self.num_x() * self.num_y() * self.num_z()
}
pub fn num_nodes(&self) -> usize {
self.xs.len() * self.ys.len() * self.zs.len()
}
pub fn grid_index(&self, i: usize, j: usize, k: usize) -> usize {
i + self.xs.len() * j + self.xs.len() * self.ys.len() * k
}
pub fn write_gmsh(&self, output: &std::path::PathBuf) -> Result<(), std::io::Error> {
use itertools::Itertools;
let mut filestream = BufWriter::new(File::create(output)?);
writeln!(&mut filestream, "$MeshFormat\n2.2 0 8\n$EndMeshFormat")?;
write!(&mut filestream, "$Nodes\n{}\n", self.num_nodes())?;
for (k, z) in self.zs.iter().enumerate() {
for (j, y) in self.ys.iter().enumerate() {
for (i, x) in self.xs.iter().enumerate() {
writeln!(
&mut filestream,
"{} {} {} {}",
self.grid_index(i, j, k) + 1,
x,
y,
z
)?;
}
}
}
writeln!(&mut filestream, "$EndNodes")?;
let mut x_nodes = 1..=self.num_nodes();
let gmsh_y_index = |x_index| x_index + self.xs.len();
let gmsh_z_index = |x_index| x_index + self.xs.len() * self.ys.len();
writeln!(&mut filestream, "$Elements\n{}", self.num_voxels())?;
for index in 0..self.num_voxels() {
if index != 0 && index % self.num_x() == 0 {
x_nodes = x_nodes.dropping(1);
}
if index != 0 && index % (self.num_x() * self.num_y()) == 0 {
x_nodes = x_nodes.dropping(self.xs.len());
}
let xl = x_nodes.next().unwrap();
let xr = xl + 1;
let yl = gmsh_y_index(xl);
let yr = yl + 1;
let zl = gmsh_z_index(xl);
let zr = zl + 1;
let yzl = gmsh_z_index(yl);
let yzr = yzl + 1;
writeln!(
&mut filestream,
"{} 5 2 0 0 {} {} {} {} {} {} {} {}",
index + 1,
xl,
xr,
yr,
yl,
zl,
zr,
yzr,
yzl,
)?;
}
writeln!(&mut filestream, "$EndElements")?;
let mut write_elt_data = |name: &str, data: &Vec<f64>| -> Result<(), std::io::Error> {
writeln!(&mut filestream, "$ElementData")?;
writeln!(&mut filestream, "1\n{}", name)?;
writeln!(&mut filestream, "1\n0.0")?;
writeln!(&mut filestream, "3\n0\n1\n{}", data.len())?;
for (index, val) in data.iter().enumerate() {
writeln!(&mut filestream, "{} {}", index + 1, val)?;
}
writeln!(&mut filestream, "$EndElementData")?;
Ok(())
};
write_elt_data(r#""Dose [Gy·cm2]""#, &self.doses)?;
write_elt_data(r#""Uncertainty fraction""#, &self.uncerts)
}
}
fn parse_simple_line<T>(line: String, title: &'static str, expect_len: usize) -> Vec<T>
where
T: FromStr,
<T as std::str::FromStr>::Err: Debug,
{
let entries: Vec<T> = line
.trim()
.split_whitespace()
.map(|num| num.parse::<T>().expect(title))
.collect();
assert!(entries.len() == expect_len);
entries
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn read_3ddose() {
let mut path = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"));
path.push("props");
path.push("water_block.3ddose");
let data = DoseBlock::from_3d_dose(&path).expect("couldn't parse 3ddose file");
assert_eq!(data.num_voxels(), 40 * 40 * 40);
assert_eq!(data.num_nodes(), 41 * 41 * 41);
assert_eq!(data.grid_index(0, 0, 0), 0);
assert_eq!(data.grid_index(0, 1, 0), 41);
assert_eq!(data.grid_index(0, 0, 1), 1681);
assert_eq!(data.grid_index(0, 0, 1), 1681);
assert_eq!(data.uncerts[21503], 0.37652693977336593);
}
}