Added homogeneous 2D mesh; added tests of 1D and homogeneous 2D geo.

.github/workflows/rustbca_compile_check.yml

@@ -23,7 +23,7 @@ jobs:
         sudo apt-get install curl
     - name: Install rust
       run: |
-        curl --proto '=https' --tlsv1.2 -sSf -y https://sh.rustup.rs | sh 
+        curl --proto '=https' --tlsv1.2 -sSf -y https://sh.rustup.rs | sh
         sudo apt-get install rustc cargo
     - name: Install pip for Python-3
       run: |
@@ -70,4 +70,6 @@ jobs:
         ./target/release/RustBCA SPHERE examples/boron_nitride_sphere.toml
         cargo run --release --features parry3d TRIMESH examples/tungsten_twist_trimesh.toml
         ./target/release/RustBCA examples/boron_nitride_wire.toml
-
+        cat boron_nitride_summary.output
+        ./target/release/RustBCA HOMOGENEOUS2D examples/boron_nitride_wire_homogeneous.toml
+        cat boron_nitride_h_summary.output

examples/boron_nitride_wire.toml

@@ -13,13 +13,13 @@ Z = [ 5, 7,]
 m = [ 10.811, 14,]
 
 [particle_parameters]
-length_unit = "MICRON"
+length_unit = "ANGSTROM"
 energy_unit = "EV"
 mass_unit = "AMU"
 N = [ 10000 ]
 m = [ 1.008 ]
 Z = [ 1 ]
-E = [ 1000.0 ]
+E = [ 5000.0 ]
 Ec = [ 1.0 ]
 Es = [ 10.0 ]
 pos = [ [ 0.0, 0.0, 0.0,] ]

examples/boron_nitride_wire_homogeneous.toml

@@ -0,0 +1,33 @@
+[options]
+name = "boron_nitride_h_"
+weak_collision_order = 0
+num_threads = 4
+
+[material_parameters]
+energy_unit = "EV"
+mass_unit = "AMU"
+Eb = [ 0.0, 0.0,]
+Es = [ 5.76, 0.0,]
+Ec = [ 1.0, 1.0,]
+Z = [ 5, 7,]
+m = [ 10.811, 14,]
+
+[particle_parameters]
+length_unit = "ANGSTROM"
+energy_unit = "EV"
+mass_unit = "AMU"
+N = [ 10000 ]
+m = [ 1.008 ]
+Z = [ 1 ]
+E = [ 5000.0 ]
+Ec = [ 1.0 ]
+Es = [ 10.0 ]
+pos = [ [ 0.0, 0.0, 0.0,] ]
+dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ]
+
+[geometry_input]
+length_unit = "ANGSTROM"
+points = [ [ 0.0, 0.0,], [ 1000.0, 0.0,], [ 998.0267284282716, 62.79051952931337,], [ 992.1147013144779, 125.33323356430427,], [ 982.2872507286887, 187.38131458572462,], [ 968.5831611286311, 248.68988716485478,], [ 951.0565162951535, 309.0169943749474,], [ 929.7764858882514, 368.12455268467795,], [ 904.8270524660195, 425.7792915650727,], [ 876.3066800438636, 481.7536741017153,], [ 844.3279255020151, 535.8267949789966,], [ 809.0169943749474, 587.7852522924732,], [ 770.5132427757892, 637.4239897486898,], [ 728.9686274214115, 684.5471059286888,], [ 684.5471059286887, 728.9686274214115,], [ 637.4239897486897, 770.5132427757893,], [ 587.785252292473, 809.0169943749474,], [ 535.8267949789965, 844.3279255020151,], [ 481.75367410171515, 876.3066800438637,], [ 425.77929156507264, 904.8270524660196,], [ 368.12455268467784, 929.7764858882515,], [ 309.01699437494744, 951.0565162951535,], [ 248.68988716485475, 968.5831611286311,], [ 187.3813145857245, 982.2872507286887,], [ 125.33323356430427, 992.1147013144779,], [ 62.79051952931331, 998.0267284282716,], [ -1.6081226496766366e-13, 1000.0,], [ -62.7905195293134, 998.0267284282716,], [ -125.33323356430436, 992.1147013144778,], [ -187.38131458572482, 982.2872507286886,], [ -248.68988716485484, 968.5831611286311,], [ -309.01699437494756, 951.0565162951535,], [ -368.12455268467795, 929.7764858882515,], [ -425.7792915650727, 904.8270524660195,], [ -481.75367410171543, 876.3066800438635,], [ -535.8267949789969, 844.327925502015,], [ -587.785252292473, 809.0169943749474,], [ -637.4239897486898, 770.5132427757893,], [ -684.5471059286888, 728.9686274214115,], [ -728.9686274214116, 684.5471059286886,], [ -770.5132427757893, 637.4239897486896,], [ -809.0169943749473, 587.7852522924733,], [ -844.3279255020151, 535.8267949789966,], [ -876.3066800438636, 481.7536741017152,], [ -904.8270524660196, 425.77929156507247,], [ -929.7764858882515, 368.1245526846778,], [ -951.0565162951536, 309.0169943749471,], [ -968.5831611286311, 248.6898871648548,], [ -982.2872507286887, 187.38131458572457,], [ -992.1147013144779, 125.3332335643041,], [ -998.0267284282716, 62.79051952931314,], [ -1000.0, -3.216245299353273e-13,], [ -998.0267284282716, -62.79051952931334,], [ -992.1147013144779, -125.33323356430428,], [ -982.2872507286886, -187.38131458572477,], [ -968.5831611286311, -248.689887164855,], [ -951.0565162951535, -309.0169943749477,], [ -929.7764858882513, -368.1245526846783,], [ -904.8270524660195, -425.77929156507264,], [ -876.3066800438636, -481.7536741017154,], [ -844.327925502015, -535.8267949789968,], [ -809.0169943749472, -587.7852522924734,], [ -770.5132427757891, -637.42398974869,], [ -728.9686274214115, -684.5471059286888,], [ -684.5471059286887, -728.9686274214115,], [ -637.4239897486896, -770.5132427757893,], [ -587.7852522924733, -809.0169943749473,], [ -535.8267949789963, -844.3279255020153,], [ -481.75367410171526, -876.3066800438636,], [ -425.7792915650722, -904.8270524660198,], [ -368.1245526846778, -929.7764858882515,], [ -309.01699437494756, -951.0565162951535,], [ -248.68988716485444, -968.5831611286312,], [ -187.38131458572462, -982.2872507286887,], [ -125.33323356430373, -992.1147013144779,], [ -62.79051952931321, -998.0267284282716,], [ -1.8369701987210297e-13, -1000.0,], [ 62.79051952931372, -998.0267284282716,], [ 125.33323356430424, -992.1147013144779,], [ 187.38131458572514, -982.2872507286886,], [ 248.68988716485492, -968.5831611286311,], [ 309.0169943749472, -951.0565162951536,], [ 368.12455268467824, -929.7764858882513,], [ 425.7792915650726, -904.8270524660196,], [ 481.7536741017157, -876.3066800438634,], [ 535.8267949789968, -844.327925502015,], [ 587.7852522924737, -809.016994374947,], [ 637.42398974869, -770.5132427757891,], [ 684.5471059286887, -728.9686274214115,], [ 728.9686274214118, -684.5471059286883,], [ 770.5132427757893, -637.4239897486897,], [ 809.0169943749478, -587.7852522924726,], [ 844.3279255020152, -535.8267949789963,], [ 876.3066800438636, -481.7536741017153,], [ 904.8270524660197, -425.77929156507224,], [ 929.7764858882515, -368.12455268467784,], [ 951.0565162951538, -309.01699437494676,], [ 968.5831611286312, -248.6898871648545,], [ 982.2872507286887, -187.38131458572468,], [ 992.1147013144779, -125.33323356430378,], [ 998.0267284282716, -62.790519529313265,], [ 1000.0, -2.4492935982947065e-13,],]
+simulation_boundary_points = [ [ -1100.0, -1100.0,], [ -1100.0, 1100.0,], [ 1100.0, 1100.0,], [ 1100.0, -1100.0,],]
+densities = [ 0.065, 0.065,]
+electronic_stopping_correction_factor = 1.0

src/enums.rs

@@ -8,7 +8,8 @@ pub enum MaterialType {
     #[cfg(feature = "parry3d")]
     BALL(material::Material<parry::ParryBall>),
     #[cfg(feature = "parry3d")]
-    TRIMESH(material::Material<parry::ParryTriMesh>)
+    TRIMESH(material::Material<parry::ParryTriMesh>),
+    HOMOGENEOUS2D(material::Material<geometry::HomogeneousMesh2D>)
 }
 
 #[derive(Deserialize, PartialEq, Clone, Copy)]
@@ -29,6 +30,7 @@ pub enum GeometryType {
     BALL,
     #[cfg(feature = "parry3d")]
     TRIMESH,
+    HOMOGENEOUS2D,
 }
 
 /// Mode of electronic stopping to use.

src/geometry.rs

@@ -265,6 +265,120 @@ impl Geometry for Mesh1D {
     }
 }
 
+#[derive(Deserialize, Clone)]
+pub struct HomogeneousMesh2DInput {
+    pub length_unit: String,
+    pub points: Vec<(f64, f64)>,
+    pub simulation_boundary_points: Vec<(f64, f64)>,
+    pub densities: Vec<f64>,
+    pub electronic_stopping_correction_factor: f64
+}
+
+#[derive(Clone)]
+pub struct HomogeneousMesh2D {
+    pub boundary: Polygon<f64>,
+    pub simulation_boundary: Polygon<f64>,
+    pub energy_barrier_thickness: f64,
+    pub densities: Vec<f64>,
+    pub electronic_stopping_correction_factor: f64,
+    pub concentrations: Vec<f64>
+}
+
+impl Geometry for HomogeneousMesh2D {
+
+    type InputFileFormat = InputHomogeneous2D;
+
+    fn new(input: &<<Self as Geometry>::InputFileFormat as GeometryInput>::GeometryInput) -> Self {
+        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 boundary_points_converted: Vec<(f64, f64)> = input.points.iter().map(|(x, y)| (x*length_unit, y*length_unit)).collect();
+
+        let simulation_boundary_points_converted: Vec<(f64, f64)> = input.simulation_boundary_points.iter().map(|(x, y)| (x*length_unit, y*length_unit)).collect();
+
+        let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor;
+
+        let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
+
+        let total_density: f64 = densities.iter().sum();
+
+        let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
+
+        let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
+
+        let boundary = Polygon::new(LineString::from(boundary_points_converted), vec![]);
+
+        let simulation_boundary =  Polygon::new(LineString::from(simulation_boundary_points_converted), vec![]);
+
+        HomogeneousMesh2D {
+            densities,
+            simulation_boundary,
+            boundary,
+            electronic_stopping_correction_factor,
+            energy_barrier_thickness,
+            concentrations
+        }
+    }
+
+    fn get_densities(&self,  x: f64, y: f64, z: f64) -> &Vec<f64> {
+        &self.densities
+    }
+
+    fn get_ck(&self,  x: f64, y: f64, z: f64) -> f64 {
+        self.electronic_stopping_correction_factor
+    }
+
+    fn get_total_density(&self,  x: f64, y: f64, z: f64) -> f64 {
+        self.densities.iter().sum::<f64>()
+    }
+
+    fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec<f64> {
+        &self.concentrations
+    }
+
+    fn inside(&self, x: f64, y: f64, z: f64) -> bool {
+        self.boundary.contains(&point!(x: x, y: y))
+    }
+
+    fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool {
+        self.simulation_boundary.contains(&point!(x: x, y: y))
+    }
+    fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool {
+        if self.inside(x, y, z) {
+            true
+        } else {
+            if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+                let distance = ((x - p.x()).powf(2.) +  (y - p.y()).powf(2.)).sqrt();
+                distance < self.energy_barrier_thickness
+            } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+                true
+            } else {
+                panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
+            }
+        }
+    }
+    fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+            (p.x(), p.y(), z)
+        } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+            return (p.x(), p.y(), z)
+        } else {
+            panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
+        }
+    }
+}
+
 /// Object that contains raw mesh input data.
 #[derive(Deserialize, Clone)]
 pub struct Mesh2DInput {
@@ -434,6 +548,8 @@ impl Geometry for Mesh2D {
     fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
         if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
             (p.x(), p.y(), z)
+        } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+            (p.x(), p.y(), z)
         } else {
             panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
         }

src/input.rs

@@ -49,6 +49,42 @@ impl InputFile for Input2D {
     }
 }
 
+/// Rustbca's internal representation of an input file for homogeneous 2D meshes.
+#[derive(Deserialize, Clone)]
+pub struct InputHomogeneous2D {
+    pub options: Options,
+    pub material_parameters: material::MaterialParameters,
+    pub particle_parameters: particle::ParticleParameters,
+    pub geometry_input: geometry::HomogeneousMesh2DInput,
+}
+
+impl GeometryInput for InputHomogeneous2D {
+    type GeometryInput = geometry::HomogeneousMesh2DInput;
+}
+
+impl InputFile for InputHomogeneous2D {
+
+    fn new(string: &str) -> InputHomogeneous2D {
+        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 Input0D {
     pub options: Options,

src/lib.rs

@@ -77,7 +77,7 @@ pub mod parry;
 pub use crate::enums::*;
 pub use crate::consts::*;
 pub use crate::structs::*;
-pub use crate::input::{Input2D, Input1D, Input0D, Options, InputFile, GeometryInput};
+pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput};
 pub use crate::output::{OutputUnits};
 pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D};
 pub use crate::sphere::{Sphere, SphereInput, InputSphere};

src/main.rs

@@ -63,9 +63,9 @@ pub mod parry;
 pub use crate::enums::*;
 pub use crate::consts::*;
 pub use crate::structs::*;
-pub use crate::input::{Input2D, Input1D, Input0D, Options, InputFile, GeometryInput};
+pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput};
 pub use crate::output::{OutputUnits};
-pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D};
+pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, HomogeneousMesh2D};
 pub use crate::sphere::{Sphere, SphereInput, InputSphere};
 pub use crate::physics::{physics_loop};
 
@@ -88,6 +88,7 @@ fn main() {
             "BALL" => GeometryType::BALL,
             #[cfg(feature = "parry3d")]
             "TRIMESH" => GeometryType::TRIMESH,
+            "HOMOGENEOUS2D" => GeometryType::HOMOGENEOUS2D,
             _ => panic!("Unimplemented geometry {}.", args[1].clone())
         }),
         _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 (<input file name>) or 2 (<geometry type> <input file name>)"),
@@ -120,5 +121,9 @@ fn main() {
             let (particle_input_array, material, options, output_units) = input::input::<ParryTriMesh>(input_file);
             physics_loop::<ParryTriMesh>(particle_input_array, material, options, output_units);
         }
+        GeometryType::HOMOGENEOUS2D => {
+            let (particle_input_array, material, options, output_units) = input::input::<geometry::HomogeneousMesh2D>(input_file);
+            physics_loop::<HomogeneousMesh2D>(particle_input_array, material, options, output_units);
+        }
     }
 }

src/tests.rs

@@ -468,13 +468,31 @@ fn test_geometry() {
         electronic_stopping_correction_factors: vec![1.0, 1.0],
     };
 
+    let geometry_input_homogeneous_2D = geometry::HomogeneousMesh2DInput {
+        length_unit: "ANGSTROM".to_string(),
+        points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)],
+        densities: vec![0.03, 0.03],
+        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)],
+        electronic_stopping_correction_factor: 1.0,
+    };
+
+    let geometry_input_1D = geometry::Mesh1DInput {
+        length_unit: "ANGSTROM".to_string(),
+        densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]],
+        layer_thicknesses: vec![thickness/2., thickness/2.],
+        electronic_stopping_correction_factors: vec![1.0, 1.0, 1.0]
+
+    };
+
     let geometry_input_0D = geometry::Mesh0DInput {
         length_unit: "ANGSTROM".to_string(),
         densities: vec![0.03, 0.03],
         electronic_stopping_correction_factor: 1.0
     };
 
     let material_2D: material::Material<geometry::Mesh2D> = material::Material::<geometry::Mesh2D>::new(&material_parameters, &geometry_input_2D);
+    let material_homogeneous_2D: material::Material<geometry::HomogeneousMesh2D> = material::Material::<geometry::HomogeneousMesh2D>::new(&material_parameters, &geometry_input_homogeneous_2D);
+    let material_1D: material::Material<geometry::Mesh1D> = material::Material::<geometry::Mesh1D>::new(&material_parameters, &geometry_input_1D);
     let mut material_0D: material::Material<geometry::Mesh0D> = material::Material::<geometry::Mesh0D>::new(&material_parameters, &geometry_input_0D);
     material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM;
 
@@ -484,13 +502,30 @@ fn test_geometry() {
     particle_1.pos_old.x = -500.*ANGSTROM;
     particle_1.pos_old.y = 0.;
 
+    assert!(material_2D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z));
+    assert!(!material_2D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z));
+
+    assert!(material_homogeneous_2D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z));
+    assert!(!material_homogeneous_2D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z));
+
+    assert!(material_0D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z));
+    assert!(!material_0D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z));
+
     assert_eq!(material_2D.geometry.get_ck(0., 0., 0.), material_0D.geometry.get_ck(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_0D.geometry.get_densities(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_0D.geometry.get_total_density(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_0D.geometry.get_concentrations(0., 0., 0.));
     assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_0D.geometry.closest_point(-10., 0., 5.));
     assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_0D.geometry.get_densities(-10., 0., 5.));
     assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_0D.geometry.get_ck(-10., 0., 5.));
+
+    assert_eq!(material_2D.geometry.get_ck(0., 0., 0.), material_homogeneous_2D.geometry.get_ck(0., 0., 0.));
+    assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_homogeneous_2D.geometry.get_densities(0., 0., 0.));
+    assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_homogeneous_2D.geometry.get_total_density(0., 0., 0.));
+    assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_homogeneous_2D.geometry.get_concentrations(0., 0., 0.));
+    assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_homogeneous_2D.geometry.closest_point(-10., 0., 5.));
+    assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_homogeneous_2D.geometry.get_densities(-10., 0., 5.));
+    assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_homogeneous_2D.geometry.get_ck(-10., 0., 5.));
 }
 
 #[test]