Parry3d

.github/workflows/rustbca_compile_check.yml

@@ -37,7 +37,7 @@ jobs:
         sudo apt install libhdf5-dev
     - name: Test RustBCA
       run: |
-        cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions
+        cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions,parry3d
     - name: Run Examples
       run: |
         cargo run --release 0D examples/boron_nitride_0D.toml
@@ -48,4 +48,4 @@ jobs:
         ./target/release/rustBCA examples/layered_geometry.toml
         cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
         ./target/release/rustBCA SPHERE examples/boron_nitride_sphere.toml
-
+        cargo run --release --features parry3d TRIMESH examples/tungsten_twist_trimesh.toml

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "rustBCA"
-version = "0.1.0"
+version = "1.1.0"
 authors = ["Jon Drobny <drobny2@illinois.edu>"]
 edition = "2018"
 
@@ -22,6 +22,7 @@ netlib-src = {version  = "0.8", optional = true}
 intel-mkl-src = {version = "0.6.0", optional = true}
 rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
 ndarray = {version = "0.14.0", features = ["serde"], optional = true}
+parry3d-f64 = {version = "0.2.0", optional = true}
 
 [dev-dependencies]
 float-cmp = "0.8.0"
@@ -39,3 +40,5 @@ cpr_rootfinder_netlib = ["rcpr", "netlib-src"]
 cpr_rootfinder_intel_mkl = ["rcpr", "intel-mkl-src"]
 distributions = ["ndarray"]
 no_list_output = []
+parry3d = ["parry3d-f64"]
+accelerated_ions = []

examples/tungsten_tiles_trimesh.toml

@@ -0,0 +1,105 @@
+[options]
+name = "tungsten_tiles_"
+track_trajectories = true
+track_recoils = true
+track_recoil_trajectories = true
+write_buffer_size = 8000
+weak_collision_order = 0
+suppress_deep_recoils = false
+high_energy_free_flight_paths = false
+num_threads = 4
+num_chunks = 5
+use_hdf5 = false
+electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
+mean_free_path_model = "LIQUID"
+interaction_potential = [["KR_C"]]
+scattering_integral = [["MENDENHALL_WELLER"]]
+root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-3}}]]
+track_displacements = false
+track_energy_losses = false
+
+[material_parameters]
+energy_unit = "EV"
+mass_unit = "AMU"
+Eb = [ 3.0,]
+Es = [ 8.79,]
+Ec = [ 1.0,]
+Z = [ 74 ]
+m = [ 183.84 ]
+interaction_index = [0]
+surface_binding_model = {"PLANAR"={calculation="TARGET"}}
+bulk_binding_model = "AVERAGE"
+
+[particle_parameters]
+length_unit = "1e-8"
+energy_unit = "EV"
+mass_unit = "AMU"
+N = [ 10 ]
+m = [ 4.008 ]
+Z = [ 2 ]
+E = [ 250.0 ]
+Ec = [ 1.0 ]
+Es = [ 0.0 ]
+pos = [ [ 0.35, 1.0, 0.5,] ]
+dir = [ [ 0.98, -0.3, 0.0,] ]
+interaction_index = [ 0 ]
+
+[geometry_input]
+length_unit = "1e-8"
+electronic_stopping_correction_factor = 1.0
+densities = [60000]
+vertices = [
+  [0.00000000, 0.00000000, 0.00000000],
+  [0.00000000, 0.00000000, 0.20000000],
+  [0.00000000, 0.00000000, 0.40000000],
+  [0.00000000, 0.00000000, 0.60000000],
+  [0.00000000, 0.00000000, 0.80000000],
+  [0.00000000, 0.00000000, 1.00000000],
+  [0.00000000, 1.00000000, 0.00000000],
+  [0.00000000, 1.00000000, 0.20000000],
+  [0.00000000, 1.00000000, 0.40000000],
+  [0.00000000, 1.00000000, 0.60000000],
+  [0.00000000, 1.00000000, 0.80000000],
+  [0.00000000, 1.00000000, 1.00000000],
+  [0.20000000, 0.55000000, 0.00000000],
+  [0.20000000, 0.55000000, 1.00000000],
+  [0.40000000, 0.20000000, 0.00000000],
+  [0.40000000, 0.20000000, 0.20000000],
+  [0.40000000, 0.20000000, 0.40000000],
+  [0.40000000, 0.20000000, 0.60000000],
+  [0.40000000, 0.20000000, 0.80000000],
+  [0.40000000, 0.20000000, 1.00000000],
+  [0.40000000, 1.00000000, 0.00000000],
+  [0.40000000, 1.00000000, 0.20000000],
+  [0.40000000, 1.00000000, 0.40000000],
+  [0.40000000, 1.00000000, 0.60000000],
+  [0.40000000, 1.00000000, 0.80000000],
+  [0.40000000, 1.00000000, 1.00000000],
+  [0.60000000, 0.20000000, 0.00000000],
+  [0.60000000, 0.20000000, 0.20000000],
+  [0.60000000, 0.20000000, 0.40000000],
+  [0.60000000, 0.20000000, 0.60000000],
+  [0.60000000, 0.20000000, 0.80000000],
+  [0.60000000, 0.20000000, 1.00000000],
+  [0.60000000, 1.00000000, 0.00000000],
+  [0.60000000, 1.00000000, 0.20000000],
+  [0.60000000, 1.00000000, 0.40000000],
+  [0.60000000, 1.00000000, 0.60000000],
+  [0.60000000, 1.00000000, 0.80000000],
+  [0.60000000, 1.00000000, 1.00000000],
+  [0.80000000, 0.55000000, 0.00000000],
+  [0.80000000, 0.55000000, 1.00000000],
+  [1.00000000, 0.00000000, 0.00000000],
+  [1.00000000, 0.00000000, 0.20000000],
+  [1.00000000, 0.00000000, 0.40000000],
+  [1.00000000, 0.00000000, 0.60000000],
+  [1.00000000, 0.00000000, 0.80000000],
+  [1.00000000, 0.00000000, 1.00000000],
+  [1.00000000, 1.00000000, 0.00000000],
+  [1.00000000, 1.00000000, 0.20000000],
+  [1.00000000, 1.00000000, 0.40000000],
+  [1.00000000, 1.00000000, 0.60000000],
+  [1.00000000, 1.00000000, 0.80000000],
+  [1.00000000, 1.00000000, 1.00000000],
+]
+indices = [[0, 14, 12], [26, 40, 38], [46, 32, 38], [20, 6, 12], [0, 40, 14], [6, 0, 12], [40, 46, 38], [40, 26, 14], [32, 26, 38], [14, 20, 12], [1, 0, 40], [1, 40, 41], [2, 1, 41], [2, 41, 42], [3, 2, 42], [3, 42, 43], [4, 3, 43], [4, 43, 44], [5, 4, 44], [5, 44, 45], [41, 40, 46], [41, 46, 47], [42, 41, 47], [42, 47, 48], [43, 42, 48], [43, 48, 49], [44, 43, 49], [44, 49, 50], [45, 44, 50], [45, 50, 51], [47, 46, 32], [47, 32, 33], [48, 47, 33], [48, 33, 34], [49, 48, 34], [49, 34, 35], [50, 49, 35], [50, 35, 36], [51, 50, 36], [51, 36, 37], [33, 32, 26], [33, 26, 27], [34, 33, 27], [34, 27, 28], [35, 34, 28], [35, 28, 29], [36, 35, 29], [36, 29, 30], [37, 36, 30], [37, 30, 31], [27, 26, 15], [26, 14, 15], [28, 27, 16], [27, 15, 16], [29, 28, 17], [28, 16, 17], [30, 29, 18], [29, 17, 18], [31, 30, 19], [30, 18, 19], [15, 14, 20], [15, 20, 21], [16, 15, 21], [16, 21, 22], [17, 16, 22], [17, 22, 23], [18, 17, 23], [18, 23, 24], [19, 18, 24], [19, 24, 25], [21, 20, 6], [21, 6, 7], [22, 21, 7], [22, 7, 8], [23, 22, 8], [23, 8, 9], [24, 23, 9], [24, 9, 10], [25, 24, 10], [25, 10, 11], [7, 6, 0], [7, 0, 1], [8, 7, 1], [8, 1, 2], [9, 8, 2], [9, 2, 3], [10, 9, 3], [10, 3, 4], [11, 10, 4], [11, 4, 5], [5, 19, 13], [31, 45, 39], [51, 37, 39], [25, 11, 13], [5, 45, 19], [11, 5, 13], [45, 51, 39], [45, 31, 19], [37, 31, 39], [19, 25, 13]]

scripts/rustbca.py

@@ -136,7 +136,7 @@ def do_trajectory_plot(name, thickness=None, depth=None, boundary=None, plot_fin
         plt.savefig(name+'trajectories_.png')
         plt.close()
 
-def do_trajectory_plot_3d(name, thickness=None, depth=None, boundary=None, plot_final_positions=True, plot_origins=True, radius=None):
+def do_trajectory_plot_3d(name, thickness=None, depth=None, boundary=None, plot_final_positions=True, plot_origins=True, radius=None, cube_length=None, input_file=None):
     '''
     Plots trajectories of ions and recoils from [name]trajectories.output.
     Optionally marks final positions/origins and draws material geometry.
@@ -154,7 +154,7 @@ def do_trajectory_plot_3d(name, thickness=None, depth=None, boundary=None, plot_
 
     '''
 
-    from mayavi.mlab import points3d, plot3d, mesh
+    from mayavi.mlab import points3d, plot3d, mesh, triangular_mesh
 
     reflected = np.atleast_2d(np.genfromtxt(name+'reflected.output', delimiter=','))
     sputtered = np.atleast_2d(np.genfromtxt(name+'sputtered.output', delimiter=','))
@@ -169,16 +169,17 @@ def do_trajectory_plot_3d(name, thickness=None, depth=None, boundary=None, plot_
 
     index = 0
     x_max = 0
-    min_length = 3
+    min_length = 1
+    scale_factor = 200.0
     if np.size(trajectories) > 0:
         for trajectory_length in trajectory_data:
 
             M = trajectories[index, 0]
             Z = trajectories[index, 1]
             E = trajectories[index:(trajectory_length + index), 2]
-            x = trajectories[index:(trajectory_length + index), 3]
-            y = trajectories[index:(trajectory_length + index), 4]
-            z = trajectories[index:(trajectory_length + index), 5]
+            x = trajectories[index:(trajectory_length + index), 3]*scale_factor
+            y = trajectories[index:(trajectory_length + index), 4]*scale_factor
+            z = trajectories[index:(trajectory_length + index), 5]*scale_factor
 
             if np.max(x) > x_max:
                 x_max = np.max(x)
@@ -193,17 +194,17 @@ def do_trajectory_plot_3d(name, thickness=None, depth=None, boundary=None, plot_
             if np.size(sputtered) > 0:
                 sputtered_colors = [colormap.to_rgba(Z)[:3] for Z in sputtered[:,1]]
                 for x, y, z, c in zip(sputtered[:,3], sputtered[:,4], sputtered[:,5], sputtered_colors):
-                    points3d(x, y, z, color=c, scale_factor=2)
+                    points3d(x*scale_factor, y*scale_factor, z*scale_factor, color=c, scale_factor=2)
 
             if np.size(reflected) > 0:
                 reflected_colors = [colormap.to_rgba(Z)[:3] for Z in reflected[:,1]]
                 for x, y, z, c in zip(reflected[:,3], reflected[:,4], reflected[:,5], reflected_colors):
-                    points3d(x, y, z, color=c, scale_factor=2)
+                    points3d(x*scale_factor, y*scale_factor, z*scale_factor, color=c, scale_factor=2)
 
             if np.size(deposited) > 0:
                 deposited_colors = [colormap.to_rgba(Z)[:3] for Z in deposited[:,1]]
                 for x, y, z, c in zip(deposited[:,2], deposited[:,3], deposited[:,4], deposited_colors):
-                    points3d(x, y, z, color=c, scale_factor=2)
+                    points3d(x*scale_factor, y*scale_factor, z*scale_factor, color=c, scale_factor=2)
 
         if boundary:
             x = [x_ for (x_, y_) in boundary]
@@ -221,6 +222,53 @@ def do_trajectory_plot_3d(name, thickness=None, depth=None, boundary=None, plot_
             z = np.cos(theta)
             mesh(radius*x, radius*y, radius*z, color=(0.1,0.7,0.3), opacity=0.2)
 
+        if cube_length:
+            faces = []
+
+            xmin = -cube_length/2.*scale_factor
+            xmax = cube_length/2.*scale_factor
+            ymin = -cube_length/2.*scale_factor
+            ymax = cube_length/2.*scale_factor
+            zmin = -cube_length/2.*scale_factor
+            zmax = cube_length/2.*scale_factor
+
+            x,y = np.mgrid[xmin:xmax:3j,ymin:ymax:3j]
+            z = np.ones(y.shape)*zmin
+            faces.append((x,y,z))
+
+            x,y = np.mgrid[xmin:xmax:3j,ymin:ymax:3j]
+            z = np.ones(y.shape)*zmax
+            faces.append((x,y,z))
+
+            x,z = np.mgrid[xmin:xmax:3j,zmin:zmax:3j]
+            y = np.ones(z.shape)*ymin
+            faces.append((x,y,z))
+
+            x,z = np.mgrid[xmin:xmax:3j,zmin:zmax:3j]
+            y = np.ones(z.shape)*ymax
+            faces.append((x,y,z))
+
+            y,z = np.mgrid[ymin:ymax:3j,zmin:zmax:3j]
+            x = np.ones(z.shape)*xmin
+            faces.append((x,y,z))
+
+            y,z = np.mgrid[ymin:ymax:3j,zmin:zmax:3j]
+            x = np.ones(z.shape)*xmax
+            faces.append((x,y,z))
+
+            for grid in faces:
+                x,y,z = grid
+                mesh(x, y, z, opacity=0.4, color=(0.1,0.7,0.3))
+
+        if input_file:
+            input = toml.load(input_file)
+            vertices = input['geometry_input']['vertices']
+            triangles = input['geometry_input']['indices']
+            x = [vertex[0]*scale_factor for vertex in vertices]
+            y = [vertex[1]*scale_factor for vertex in vertices]
+            z = [vertex[2]*scale_factor for vertex in vertices]
+            triangular_mesh(x, y, z, triangles, opacity=0.3, color=(0.1, 0.7, 0.3), representation='surface')
+
 
 def generate_rustbca_input(Zb, Mb, n, Eca, Ecb, Esa, Esb, Eb, Ma, Za, E0, N, N_, theta,
     thickness, depth, track_trajectories=True, track_recoils=True,

src/bca.rs

@@ -84,6 +84,13 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
             //Choose impact parameters and azimuthal angles for all collisions, and determine mean free path
             let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material, &options);
 
+            #[cfg(feature = "accelerated_ions")]
+            if !material.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z) {
+                let (x, y, z) = material.geometry.closest_point(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z);
+                let distance_to = ((x - particle_1.pos.x).powi(2) + (y - particle_1.pos.y).powi(2) + (z - particle_1.pos.z).powi(2)).sqrt();
+                mfp += distance_to;
+            }
+
             let mut total_energy_loss = 0.;
             let mut total_asymptotic_deflection = 0.;
             let mut normalized_distance_of_closest_approach = 0.;
@@ -330,9 +337,11 @@ pub fn determine_mfp_phi_impact_parameter<T: Geometry>(particle_1: &mut particle
             mfp *= -rand::random::<f64>().ln();
         }
 
+
         for k in 0..(options.weak_collision_order + 1) {
             binary_collision_geometries.push(BinaryCollisionGeometry::new(phis_azimuthal[k], impact_parameters[k], mfp))
         }
+
         return binary_collision_geometries;
     }
 }

src/enums.rs

@@ -5,6 +5,10 @@ pub enum MaterialType {
     MESH1D(material::Material<geometry::Mesh1D>),
     MESH2D(material::Material<geometry::Mesh2D>),
     SPHERE(material::Material<sphere::Sphere>),
+    #[cfg(feature = "parry3d")]
+    BALL(material::Material<parry::ParryBall>),
+    #[cfg(feature = "parry3d")]
+    TRIMESH(material::Material<parry::ParryTriMesh>)
 }
 
 #[derive(Deserialize)]
@@ -13,6 +17,10 @@ pub enum GeometryType {
     MESH1D,
     MESH2D,
     SPHERE,
+    #[cfg(feature = "parry3d")]
+    BALL,
+    #[cfg(feature = "parry3d")]
+    TRIMESH,
 }
 
 /// Mode of electronic stopping to use.

src/geometry.rs

@@ -427,7 +427,7 @@ impl Geometry for Mesh2D {
     fn get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec<f64> {
         self.nearest_to(x, y, z).get_densities()
     }
-    
+
     fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool {
         self.simulation_boundary.contains(&point!(x: x, y: y))
     }

src/main.rs

@@ -60,6 +60,9 @@ pub mod consts;
 pub mod structs;
 pub mod sphere;
 
+#[cfg(feature = "parry3d")]
+pub mod parry;
+
 pub use crate::enums::*;
 pub use crate::consts::*;
 pub use crate::structs::*;
@@ -68,6 +71,9 @@ pub use crate::output::{OutputUnits};
 pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D};
 pub use crate::sphere::{Sphere, SphereInput, InputSphere};
 
+#[cfg(feature = "parry3d")]
+pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, ParryTriMeshInput, InputParryTriMesh};
+
 fn physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) {
 
         println!("Processing {} ions...", particle_input_array.len());
@@ -158,6 +164,10 @@ fn main() {
             "1D" => GeometryType::MESH1D,
             "2D" => GeometryType::MESH2D,
             "SPHERE" => GeometryType::SPHERE,
+            #[cfg(feature = "parry3d")]
+            "BALL" => GeometryType::BALL,
+            #[cfg(feature = "parry3d")]
+            "TRIMESH" => GeometryType::TRIMESH,
             _ => 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>)"),
@@ -179,6 +189,16 @@ fn main() {
         GeometryType::SPHERE => {
             let (particle_input_array, material, options, output_units) = input::input::<Sphere>(input_file);
             physics_loop::<Sphere>(particle_input_array, material, options, output_units);
+        },
+        #[cfg(feature = "parry3d")]
+        GeometryType::BALL => {
+            let (particle_input_array, material, options, output_units) = input::input::<ParryBall>(input_file);
+            physics_loop::<ParryBall>(particle_input_array, material, options, output_units);
+        }
+        #[cfg(feature = "parry3d")]
+        GeometryType::TRIMESH => {
+            let (particle_input_array, material, options, output_units) = input::input::<ParryTriMesh>(input_file);
+            physics_loop::<ParryTriMesh>(particle_input_array, material, options, output_units);
         }
     }
 }