diff --git a/src/enums.rs b/src/enums.rs index 908c040..fbe2b13 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -5,6 +5,7 @@ pub enum MaterialType { MESH1D(material::Material), MESH2D(material::Material), SPHERE(material::Material), + SPHEREINCUBOID(material::Material), #[cfg(feature = "parry3d")] BALL(material::Material), #[cfg(feature = "parry3d")] @@ -26,6 +27,7 @@ pub enum GeometryType { MESH1D, MESH2D, SPHERE, + SPHEREINCUBOID, #[cfg(feature = "parry3d")] BALL, #[cfg(feature = "parry3d")] diff --git a/src/lib.rs b/src/lib.rs index 3a9f53e..40eb996 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -66,6 +66,7 @@ pub mod enums; pub mod consts; pub mod structs; pub mod sphere; +pub mod sphereincuboid; #[cfg(feature = "parry3d")] pub mod parry; @@ -77,6 +78,7 @@ pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, I pub use crate::output::{OutputUnits}; pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; +pub use crate::sphereincuboid::{SphereInCuboid, SphereInCuboidInput, InputSphereInCuboid}; #[cfg(feature = "parry3d")] pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, ParryTriMeshInput, InputParryTriMesh}; diff --git a/src/main.rs b/src/main.rs index 7da24c6..3bc4342 100644 --- a/src/main.rs +++ b/src/main.rs @@ -48,6 +48,7 @@ pub mod enums; pub mod consts; pub mod structs; pub mod sphere; +pub mod sphereincuboid; pub mod physics; #[cfg(feature = "parry3d")] @@ -60,6 +61,7 @@ pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, I pub use crate::output::{OutputUnits}; pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, HomogeneousMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; +pub use crate::sphereincuboid::{SphereInCuboid, SphereInCuboidInput, InputSphereInCuboid}; pub use crate::physics::{physics_loop}; #[cfg(feature = "parry3d")] @@ -77,6 +79,7 @@ fn main() { "1D" => GeometryType::MESH1D, "2D" => GeometryType::MESH2D, "SPHERE" => GeometryType::SPHERE, + "SPHEREINCUBOID" => GeometryType::SPHEREINCUBOID, #[cfg(feature = "parry3d")] "BALL" => GeometryType::BALL, #[cfg(feature = "parry3d")] @@ -104,6 +107,10 @@ fn main() { let (particle_input_array, material, options, output_units) = input::input::(input_file); physics_loop::(particle_input_array, material, options, output_units); }, + GeometryType::SPHEREINCUBOID => { + let (particle_input_array, material, options, output_units) = input::input::(input_file); + physics_loop::(particle_input_array, material, options, output_units); + }, #[cfg(feature = "parry3d")] GeometryType::BALL => { let (particle_input_array, material, options, output_units) = input::input::(input_file); diff --git a/src/sphereincuboid.rs b/src/sphereincuboid.rs new file mode 100644 index 0000000..0add566 --- /dev/null +++ b/src/sphereincuboid.rs @@ -0,0 +1,207 @@ +use super::*; + +#[derive(Deserialize, Clone)] +pub struct InputSphereInCuboid { + pub options: Options, + pub material_parameters: material::MaterialParameters, + pub particle_parameters: particle::ParticleParameters, + pub geometry_input: sphereincuboid::SphereInCuboidInput, +} + +impl InputFile for InputSphereInCuboid { + + fn new(string: &str) -> InputSphereInCuboid { + toml::from_str(string).context("Could not parse TOML file. Be sure you are using the correct input file mode (e.g., ./RustBCA SPHERE sphere.toml or RustBCA.exe 0D mesh_0d.toml).").unwrap() + } + + fn get_options(&self) -> &Options{ + &self.options + } + fn get_material_parameters(&self) -> &material::MaterialParameters{ + &self.material_parameters + } + fn get_particle_parameters(&self) -> &particle::ParticleParameters{ + &self.particle_parameters + } + fn get_geometry_input(&self) -> &Self::GeometryInput{ + &self.geometry_input + } +} + +#[derive(Deserialize, Clone)] +pub struct SphereInCuboidInput { + pub length_unit: String, + + pub sphere_radius: f64, + pub sphere_densities: Vec, + pub sphere_electronic_stopping_correction_factor: f64, + + pub cuboid_corners: Vec<(f64, f64, f64)>, + pub cuboid_densities: Vec, + pub cuboid_electronic_stopping_correction_factor: f64, +} + +#[derive(Clone)] +pub struct SphereInCuboid { + pub sphere_radius: f64, + pub sphere_densities: Vec, + pub sphere_concentrations: Vec, + pub sphere_electronic_stopping_correction_factor: f64, + + pub cuboid_corners: Vec<(f64, f64, f64)>, + pub cuboid_densities: Vec, + pub cuboid_concentrations: Vec, + pub cuboid_electronic_stopping_correction_factor: f64, + + pub energy_barrier_thickness: f64, +} + +impl GeometryInput for InputSphereInCuboid { + type GeometryInput = SphereInCuboidInput; +} + +impl Geometry for SphereInCuboid { + + type InputFileFormat = InputSphereInCuboid; + + fn new(input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> SphereInCuboid { + + let length_unit: f64 = match input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "MM" => MM, + "ANGSTROM" => ANGSTROM, + "NM" => NM, + "M" => 1., + _ => input.length_unit.parse() + .expect(format!( + "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + &input.length_unit.as_str() + ).as_str()), + }; + + let sphere_radius: f64 = input.sphere_radius*length_unit; + let sphere_densities: Vec = input.sphere_densities + .iter().map(|element| element/(length_unit).powi(3)).collect(); + let sphere_total_density: f64 = sphere_densities.iter().sum(); + let sphere_concentrations: Vec = sphere_densities + .iter().map(|&density| density/sphere_total_density).collect::>(); + let sphere_electronic_stopping_correction_factor: f64 = input.sphere_electronic_stopping_correction_factor; + + let cuboid_corners: Vec<(f64, f64, f64)> = input.cuboid_corners + .iter().map( + |&(x, y, z)| (x * length_unit, y*length_unit, z*length_unit) + ).collect(); + let cuboid_densities: Vec = input.cuboid_densities + .iter().map(|element| element/(length_unit).powi(3)).collect(); + let cuboid_total_density: f64 = cuboid_densities.iter().sum(); + let cuboid_concentrations: Vec = cuboid_densities + .iter().map(|&density| density/cuboid_total_density).collect::>(); + let cuboid_electronic_stopping_correction_factor: f64 = input.cuboid_electronic_stopping_correction_factor; + + let energy_barrier_thickness: f64 = cuboid_total_density.powf(-1./3.)/SQRTPI*2.; + + SphereInCuboid { + sphere_radius, + sphere_densities, + sphere_concentrations, + sphere_electronic_stopping_correction_factor, + + cuboid_corners, + cuboid_densities, + cuboid_concentrations, + cuboid_electronic_stopping_correction_factor, + + energy_barrier_thickness, + } + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + let r: f64 = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + if r > self.sphere_radius { + &self.cuboid_densities + } else { + &self.sphere_densities + } + } + + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { + let r: f64 = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + if r > self.sphere_radius { + self.cuboid_electronic_stopping_correction_factor + } else { + self.sphere_electronic_stopping_correction_factor + } + } + + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64{ + let r: f64 = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + if r > self.sphere_radius { + self.cuboid_densities.iter().sum() + } else { + self.sphere_densities.iter().sum() + } + } + + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + let r: f64 = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + if r > self.sphere_radius { + &self.cuboid_concentrations + } else { + &self.sphere_concentrations + } + } + + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + let (xlo, ylo, zlo): (f64, f64, f64) = self.cuboid_corners[0]; + let (xhi, yhi, zhi): (f64, f64, f64) = self.cuboid_corners[1]; + + xlo <= x && x <= xhi && + ylo <= y && y <= yhi && + zlo <= z && z <= zhi + } + + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + let (xlo, ylo, zlo): (f64, f64, f64) = self.cuboid_corners[0]; + let (xhi, yhi, zhi): (f64, f64, f64) = self.cuboid_corners[1]; + let ebt: f64 = self.energy_barrier_thickness; + + (xlo - 10.0 * ebt) <= x && x <= (xhi + 10.0 * ebt) && + (ylo - 10.0 * ebt) <= y && y <= (yhi + 10.0 * ebt) && + (zlo - 10.0 * ebt) <= z && z <= (zhi + 10.0 * ebt) + } + + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + let (xlo, ylo, zlo): (f64, f64, f64) = self.cuboid_corners[0]; + let (xhi, yhi, zhi): (f64, f64, f64) = self.cuboid_corners[1]; + let ebt: f64 = self.energy_barrier_thickness; + + (xlo - ebt) <= x && x <= (xhi + ebt) && + (ylo - ebt) <= y && y <= (yhi + ebt) && + (zlo - ebt) <= z && z <= (zhi + ebt) + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let (xlo, ylo, zlo): (f64, f64, f64) = self.cuboid_corners[0]; + let (xhi, yhi, zhi): (f64, f64, f64) = self.cuboid_corners[1]; + + if self.inside(x, y, z) { + let faces: [((f64, f64, f64), f64); 6] = [ + ((xlo, y, z), x - xlo), + ((xhi, y, z), xhi - x), + ((x, ylo, z), y - ylo), + ((x, yhi, z), yhi - y), + ((x, y, zlo), z - zlo), + ((x, y, zhi), zhi - z), + ]; + + let closest: ((f64, f64, f64), f64) = *faces + .iter() + .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap()) + .unwrap(); + closest.0 + } else { + (x.max(xlo).min(xhi), y.max(ylo).min(yhi), z.max(zlo).min(zhi)) + } + } +}