From aac64373d198413bc0e6457cf12a7f9ec1e27a39 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 17 Jul 2024 10:03:06 -0700 Subject: [PATCH 01/13] Draft version of parry2D. --- src/geometry.rs | 386 +++++++++++++++++++++++++++++++++++++++++++++--- src/main.rs | 6 +- src/material.rs | 33 ++++- src/parry.rs | 25 ++++ src/particle.rs | 17 ++- src/sphere.rs | 9 ++ src/structs.rs | 6 +- src/tests.rs | 186 ++++++++++++++++++++++- 8 files changed, 633 insertions(+), 35 deletions(-) diff --git a/src/geometry.rs b/src/geometry.rs index 8e3d31a1..8f41ebe3 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -1,5 +1,14 @@ use super::*; +use itertools::Itertools; + +use parry2d_f64::shape::{Polyline, TriMesh, FeatureId::*}; +use parry2d_f64::shape::SegmentPointLocation::{OnEdge, OnVertex}; +use parry2d_f64::query::{PointQuery, Ray, RayCast}; +use parry2d_f64::math::{Isometry, Vector}; +use parry2d_f64::math::Point as Point2; +use parry2d_f64::bounding_volume::aabb; + use geo::algorithm::contains::Contains; use geo::{Polygon, LineString, Point, point, Closest}; use geo::algorithm::closest_point::ClosestPoint; @@ -18,6 +27,7 @@ pub trait Geometry { fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool; fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool; fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); } pub trait GeometryElement: Clone { @@ -108,6 +118,10 @@ impl Geometry for Mesh0D { fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { (0., y, z) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + (-1., 0., 0.) + } } #[derive(Deserialize, Clone)] @@ -263,6 +277,164 @@ impl Geometry for Mesh1D { (self.top, y, z) } } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + (-1., 0., 0.) + } +} + +#[derive(Deserialize, Clone)] +pub struct ParryHomogeneousMesh2DInput { + pub length_unit: String, + pub points: Vec<(f64, f64)>, + pub simulation_boundary_points: Vec<(f64, f64)>, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64 +} + +#[derive(Clone)] +pub struct ParryHomogeneousMesh2D { + pub boundary: Polyline, + pub simulation_boundary: Polyline, + pub energy_barrier_thickness: f64, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64, + pub concentrations: Vec +} + + +impl Geometry for ParryHomogeneousMesh2D { + + type InputFileFormat = InputHomogeneous2D; + + fn new(input: &<::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_converted2: Vec> = input.points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + let number_boundary_points = boundary_points_converted2.len() as u32; + + for p in boundary_points_converted2.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let simulation_boundary_points_converted2: Vec> = input.simulation_boundary_points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + + + for p in simulation_boundary_points_converted2.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let test_ccw = (0..number_boundary_points as usize) + .map(|i| (boundary_points_converted2[(i + 1) % number_boundary_points as usize].x - boundary_points_converted2[i].x)*(boundary_points_converted2[i].y + boundary_points_converted2[(i + 1) % number_boundary_points as usize].y)) + .sum::(); + + assert!(test_ccw <= 0.0, "Input error: signed area of boundary is less than or equal to zero"); + + let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; + + let densities: Vec = 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 = densities.iter().map(|&density| density/total_density).collect::>(); + + + + let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_boundary_points.push([number_boundary_points - 1, 0]); + let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_boundary_points)); + + let number_simulation_boundary_points = simulation_boundary_points_converted2.clone().len() as u32; + let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); + let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted2, Some(linked_simulation_boundary_points)); + + ParryHomogeneousMesh2D { + densities, + simulation_boundary: simulation_boundary2, + boundary: boundary2, + electronic_stopping_correction_factor, + energy_barrier_thickness, + concentrations + } + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + &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::() + } + + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.concentrations + } + + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2::new(x, y); + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + } + + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2::new(x, y); + + let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + + } + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + let p = Point2::new(x, y); + (self.boundary.distance_to_local_point(&p, true) as f64) < self.energy_barrier_thickness + } + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2::new(x, y); + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_, y_) = (point_projection.point.x, point_projection.point.y); + (x_ as f64, y_ as f64, z) + } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2::new(x, y); + + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z); + + let dx = x_intersect - x; + let dy = y_intersect - y; + let dz = z_intersect - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + + if point_projection.is_inside { + (dx/mag, dy/mag, dz/mag) + } else { + (-dx/mag, -dy/mag, -dz/mag) + } + } } #[derive(Deserialize, Clone)] @@ -377,6 +549,10 @@ impl Geometry for HomogeneousMesh2D { panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); } } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + panic!("Not implemented.") + } } /// Object that contains raw mesh input data. @@ -392,6 +568,182 @@ pub struct Mesh2DInput { pub energy_barrier_thickness: f64, } +#[derive(Clone)] +pub struct ParryMesh2D { + trimesh: TriMesh, + densities: Vec>, + electronic_stopping_correction_factors: Vec, + concentrations: Vec>, + pub boundary: Polyline, + pub simulation_boundary: Polyline, + pub energy_barrier_thickness: f64 +} + +impl ParryMesh2D { + fn get_id(&self, x: f64, y: f64, z: f64) -> usize { + let p = Point2::new(x, y); + let (point_projection, feature_id) = self.trimesh.project_local_point_and_get_feature(&p); + match feature_id { + Vertex(vertex_id) => { + for (triangle_id, triangle) in self.trimesh.triangles().enumerate() { + if triangle.contains_local_point(&self.trimesh.vertices()[vertex_id as usize]) { + return triangle_id as usize + } + } + panic!("Geometry error: point ({}, {}) located on a vertex that is not associated with any triangle. Check mesh input.", x, y); + }, + Face(triangle_id) => { + triangle_id as usize + } + Unknown => { + panic!("Geometry error: unknown feature detected near ({}, {}). Check mesh input.", x, y); + } + } + } +} + +impl Geometry for ParryMesh2D { + type InputFileFormat = Input2D; + + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + let p = Point2::new(x, y); + if self.trimesh.distance_to_local_point(&p, true) < self.energy_barrier_thickness { + true + } else { + false + } + } + } + + fn inside_simulation_boundary (&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2::new(x, y); + + let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p); + + point_projection.is_inside + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2::new(x, y); + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_, y_) = (point_projection.point.x, point_projection.point.y); + (x_ as f64, y_ as f64, z) + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities[self.get_id(x, y, z)] + } + + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factors[self.get_id(x, y, z)] + } + + /// Determine the total number density of the triangle that contains or is nearest to (x, y). + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 { + self.get_densities(x, y, z).iter().sum::() + } + + /// Find the concentrations of the triangle that contains or is nearest to (x, y). + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.concentrations[self.get_id(x, y, z)] + } + + /// Determines whether the point (x, y) is inside the mesh. + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2::new(x, y); + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2::new(x, y); + + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z); + + let dx = x_intersect - x; + let dy = y_intersect - y; + let dz = z_intersect - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + + if point_projection.is_inside { + (dx/mag, dy/mag, dz/mag) + } else { + (-dx/mag, -dy/mag, -dz/mag) + } + } + + fn new(geometry_input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> ParryMesh2D { + let length_unit: f64 = match geometry_input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "MM" => MM, + "ANGSTROM" => ANGSTROM, + "NM" => NM, + "M" => 1., + _ => geometry_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", + &geometry_input.length_unit.as_str() + ).as_str()), + }; + + let boundary_points_converted2: Vec> = geometry_input.points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + let number_boundary_points = boundary_points_converted2.len() as u32; + + for p in boundary_points_converted2.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let simulation_boundary_points_converted2: Vec> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + + for p in simulation_boundary_points_converted2.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let test_ccw = (0..number_boundary_points as usize) + .map(|i| (boundary_points_converted2[(i + 1) % number_boundary_points as usize].x - boundary_points_converted2[i].x)*(boundary_points_converted2[i].y + boundary_points_converted2[(i + 1) % number_boundary_points as usize].y)) + .sum::(); + + assert!(test_ccw <= 0.0, "Input error: signed area of boundary is less than or equal to zero"); + + let triangles = geometry_input.triangles.iter().map(|(i, j, k)| [*i as u32, *j as u32, *k as u32]).collect::>(); + let points_converted = geometry_input.points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect::>>(); + + let trimesh = TriMesh::new(points_converted, triangles); + + let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone(); + + let densities: Vec> = geometry_input.densities.iter().map(|density_list| density_list.iter().map(|density| density / (length_unit).powi(3)).collect::>()).collect::>>(); + + let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit; + + let concentrations: Vec> = densities.iter().map(|density_list| density_list.iter().map(|density| density / density_list.iter().sum::() ).collect::>()).collect::>>(); + + let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_boundary_points.push([number_boundary_points - 1, 0]); + let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_boundary_points)); + + let number_simulation_boundary_points = simulation_boundary_points_converted2.clone().len() as u32; + let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); + let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted2, Some(linked_simulation_boundary_points)); + + ParryMesh2D { + trimesh, + densities, + simulation_boundary: simulation_boundary2, + boundary: boundary2, + electronic_stopping_correction_factors, + energy_barrier_thickness, + concentrations + } + } +} + /// Triangular mesh for rustbca. #[derive(Clone)] pub struct Mesh2D { @@ -489,26 +841,6 @@ impl Geometry for Mesh2D { cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck)); } - /* - for ((coordinate_set, densities), ck) in triangles.iter().zip(densities).zip(electronic_stopping_correction_factors) { - let coordinate_set_converted = ( - coordinate_set.0*length_unit, - coordinate_set.1*length_unit, - coordinate_set.2*length_unit, - coordinate_set.3*length_unit, - coordinate_set.4*length_unit, - coordinate_set.5*length_unit, - ); - - let total_density: f64 = densities.iter().sum(); - let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); - - cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck)); - } - */ - - - let mut boundary_points_converted = Vec::with_capacity(material_boundary_point_indices.len()); for index in material_boundary_point_indices.iter() { boundary_points_converted.push((points[*index].0*length_unit, points[*index].1*length_unit)); @@ -547,8 +879,18 @@ 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)) { + if p.x() == x && p.y() == y { + dbg!(&self.boundary); + println!("single point"); + panic!("({} {} {}) ({} {} {})", p.x(), p.y(), z, x, y, z); + } (p.x(), p.y(), z) } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + if p.x() == x && p.y() == y { + dbg!(&self.boundary); + println!("intersection"); + panic!("({} {} {}) ({} {} {})", p.x(), p.y(), z, x, y, z); + } (p.x(), p.y(), z) } else { panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); @@ -615,6 +957,10 @@ impl Geometry for Mesh2D { fn inside(&self, x: f64, y: f64, z: f64) -> bool { self.boundary.contains(&Point::new(x, y)) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + panic!("Not implemented."); + } } /// A mesh cell that contains a triangle and the local number densities and concentrations. diff --git a/src/main.rs b/src/main.rs index 7da24c67..c371b7ca 100644 --- a/src/main.rs +++ b/src/main.rs @@ -58,7 +58,7 @@ pub use crate::consts::*; pub use crate::structs::*; 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, HomogeneousMesh2D}; +pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, HomogeneousMesh2D, ParryHomogeneousMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; pub use crate::physics::{physics_loop}; @@ -115,8 +115,8 @@ fn main() { physics_loop::(particle_input_array, material, options, output_units); } GeometryType::HOMOGENEOUS2D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); + let (particle_input_array, material, options, output_units) = input::input::(input_file); + physics_loop::(particle_input_array, material, options, output_units); } } } diff --git a/src/material.rs b/src/material.rs index e6931fa6..1367f564 100644 --- a/src/material.rs +++ b/src/material.rs @@ -300,7 +300,7 @@ impl Material { let stopping_power = match electronic_stopping_mode { //Biersack-Varelas Interpolation - ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low*ck), + ElectronicStoppingMode::INTERPOLATED => ck/(1./S_high + 1./S_low), //Oen-Robinson ElectronicStoppingMode::LOW_ENERGY_LOCAL => S_low*ck, //Lindhard-Scharff @@ -339,14 +339,28 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, let leaving = !inside_now & inside_old; let entering = inside_now & !inside_old; + assert!(!x.is_nan(), "Particle x is NaN before surface binding energy calculation."); + assert!(!y.is_nan(), "Particle y is NaN before surface binding energy calculation."); + assert!(!z.is_nan(), "Particle z is NaN before surface binding energy calculation."); + + assert!(!cosx.is_nan(), "Particle dirx is NaN before surface binding energy calculation."); + assert!(!cosy.is_nan(), "Particle diry is NaN before surface binding energy calculation."); + assert!(!cosz.is_nan(), "Particle dirz is NaN before surface binding energy calculation."); + if entering { if particle_1.backreflected { particle_1.backreflected = false; } else { - let (x2, y2, z2) = material.closest_point(x, y, z); - let dx = x2 - x; - let dy = y2 - y; - let dz = z2 - z; + let (x2, y2, z2) = material.closest_point(x_old, y_old, z_old); + assert!(!x2.is_nan(), "Closest x is NaN in normal calculation."); + assert!(!y2.is_nan(), "Closest y is NaN in normal calculation."); + assert!(!z2.is_nan(), "Closest z is NaN in normal calculation."); + let dx = x2 - x_old; + let dy = y2 - y_old; + let dz = z2 - z_old; + assert!(!dx.is_nan(), "dx is NaN in normal calculation."); + assert!(!dy.is_nan(), "dy is NaN in normal calculation."); + assert!(!dz.is_nan(), "dz is NaN in normal calculation."); let mag = (dx*dx + dy*dy + dz*dz).sqrt(); let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; @@ -360,10 +374,17 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, if leaving { let (x2, y2, z2) = material.closest_point(x, y, z); + assert!(!x2.is_nan(), "Closest x is NaN in normal calculation."); + assert!(!y2.is_nan(), "Closest y is NaN in normal calculation."); + assert!(!z2.is_nan(), "Closest z is NaN in normal calculation."); let dx = x2 - x; let dy = y2 - y; let dz = z2 - z; + assert!(!dx.is_nan(), "dx is NaN in normal calculation."); + assert!(!dy.is_nan(), "dy is NaN in normal calculation."); + assert!(!dz.is_nan(), "dz is NaN in normal calculation."); let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; let leaving_energy = match material.surface_binding_model { @@ -414,4 +435,6 @@ pub fn boundary_condition_planar(particle_1: &mut particle::Particl if !particle_1.stopped & !particle_1.left { surface_binding_energy(particle_1, material); } + + assert!(!particle_1.pos.x.is_nan(), "Particle position is NaN: ({}, {}, {}) old: ({}, {}, {}), left? {} stopped? {}", x, y, z, particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z, particle_1.left, particle_1.stopped); } diff --git a/src/parry.rs b/src/parry.rs index 86866836..2a337212 100644 --- a/src/parry.rs +++ b/src/parry.rs @@ -124,6 +124,15 @@ impl Geometry for ParryBall { let (x_, y_, z_) = (point_projection.point.x, point_projection.point.y, point_projection.point.z); (x_ as f64, y_ as f64, z_ as f64) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let r = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + let R = self.radius; + let ux = x/r; + let uy = y/r; + let uz = z/r; + (ux, uy, uz) + } } @@ -263,6 +272,22 @@ impl Geometry for ParryTriMesh { let (x_, y_, z_) = (point_projection.point.x, point_projection.point.y, point_projection.point.z); (x_ as f64, y_ as f64, z_ as f64) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point::new(x, y, z); + let point_projection = self.trimesh.project_local_point(&p, false); + + let dx = point_projection.point.x - x; + let dy = point_projection.point.y - y; + let dz = point_projection.point.z - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + + if point_projection.is_inside { + (dx/mag, dy/mag, dz/mag) + } else { + (-dx/mag, -dy/mag, -dz/mag) + } + } } fn inside_trimesh(trimesh: &TriMesh, x: f64, y: f64, z: f64) -> bool { diff --git a/src/particle.rs b/src/particle.rs index 9de3daa8..ba3bdeb3 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -65,7 +65,7 @@ pub struct ParticleInput { } /// Particle object. Particles in rustbca include incident ions and material atoms. -#[derive(Clone)] +#[derive(Clone, Debug)] pub struct Particle { pub m: f64, pub Z: f64, @@ -271,6 +271,9 @@ impl Particle { //In order to keep average denisty constant, must add back previous asymptotic deflection let distance_traveled = mfp + self.asymptotic_deflection - asymptotic_deflection; + assert!(!distance_traveled.is_nan(), "Travel distance is NaN, mfp: {}, asymp. defl: {}, old asymp. defl: {}", mfp, self.asymptotic_deflection, asymptotic_deflection); + assert!(!self.dir_old.x.is_nan(), "Particle direction is NaN, dir: ({}, {}, {})", self.dir_old.x, self.dir_old.y, self.dir_old.z); + //dir has been updated, so use previous direction to advance in space self.pos.x += self.dir_old.x*distance_traveled; self.pos.y += self.dir_old.y*distance_traveled; @@ -294,10 +297,20 @@ pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) { let u1x = (E/(E + Es)).sqrt()*particle.dir.x + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.x; let u1y = (E/(E + Es)).sqrt()*particle.dir.y + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.y; let u1z = (E/(E + Es)).sqrt()*particle.dir.z + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.z; + + particle.E += Es; + + if u1x.is_nan() { + println!("normal: ({}, {}, {})", normal.x, normal.y, normal.z); + dbg!(&particle); + + } + particle.dir.x = u1x; particle.dir.y = u1y; particle.dir.z = u1z; - particle.E += Es; + + assert!(!u1x.is_nan(), "Particle direction after surface refraction is NaN. dir: ({}, {}, {}), old dir: ({}, {}, {}), normal: ({}, {}, {}), costheta: {}", u1x, u1y, u1z, particle.dir_old.x, particle.dir_old.y, particle.dir_old.z, normal.x, normal.y, normal.z, costheta); } /// Calcualte the refraction angle based on the surface binding energy of the material. diff --git a/src/sphere.rs b/src/sphere.rs index 92553b61..6f4055ac 100644 --- a/src/sphere.rs +++ b/src/sphere.rs @@ -120,4 +120,13 @@ impl Geometry for Sphere { let uz = z/r; (ux*R, uy*R, uz*R) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let r = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + let R = self.radius; + let ux = x/r; + let uy = y/r; + let uz = z/r; + (ux, uy, uz) + } } diff --git a/src/structs.rs b/src/structs.rs index 47f20a0a..fa5993d5 100644 --- a/src/structs.rs +++ b/src/structs.rs @@ -1,5 +1,5 @@ /// 3D vector. -#[derive(Clone, Copy)] +#[derive(Clone, Copy, Debug)] pub struct Vector { pub x: f64, pub y: f64, @@ -45,7 +45,7 @@ impl Vector { } /// TrajectoryElement is a trajectory-tracking object that includes x, y, z, and the current energy. -#[derive(Clone)] +#[derive(Clone, Debug)] pub struct TrajectoryElement { pub E: f64, pub x: f64, @@ -65,7 +65,7 @@ impl TrajectoryElement { } /// Energy loss is an output tracker that tracks the separate nuclear and electronic energy losses. -#[derive(Clone)] +#[derive(Clone, Debug)] pub struct EnergyLoss { pub En: f64, pub Ee: f64, diff --git a/src/tests.rs b/src/tests.rs index fd263f00..d70e31c6 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -3,6 +3,121 @@ use super::*; #[cfg(test)] use float_cmp::*; +use geo::algorithm::contains::Contains; +use geo::{Polygon, LineString, Point, point, Closest}; +use geo::algorithm::closest_point::ClosestPoint; + +use parry2d_f64::shape::Polyline; +use parry2d_f64::query::{PointQuery, PointProjection, Ray, RayCast}; +use parry2d_f64::math::{Isometry}; +use parry2d_f64::math::Point as Point2d; +use parry2d_f64::bounding_volume::aabb; + +/* + let number_boundary_points = boundary_points_converted.clone().len() as u32; + let boundary = Polygon::new(LineString::from(boundary_points_converted), vec![]); + let mut linked_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); + let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_points)); +*/ + +#[test] +fn test_parry2d() { + let points: Vec> = vec![Point2d::new(-1.0, -1.0), Point2d::new(-1.0, 1.0), Point2d::new(1.0, 1.0), Point2d::new(1.0, -1.0)]; + //let mut indices: Vec<[u32; 2]> = vec![]; + //let mut indices: Vec<[u32; 2]> = vec![[0, 1], [1, 2], [2, 3], [3, 0]]; + let mut indices: Vec<[u32; 2]> = vec![[0, 3], [3, 2], [2, 1], [1, 0]]; + //indices.reverse(); + + //(0.0000006261797114005236, 0.000006009591179670447) + let query_point = Point2d::new(0.0, 0.0); + + let polyline = Polyline::new(points, Some(indices)); + + let (point_projection, (_, _)) = polyline.project_local_point_assuming_solid_interior_ccw( + query_point + ); + + assert!(point_projection.is_inside); + + assert!(polyline.aabb(&Isometry::identity()).contains_local_point(&query_point)); + //assert!(point_projection.is_inside); + + let test_geometry = vec![(0.25, 0.2), (0.75, 0.2), (0.8, 0.8), (0.5, 0.5), (0.2, 0.8)]; + let num_segments = 5; + + let geometry_input_homogeneous_2D = geometry::HomogeneousMesh2DInput { + length_unit: "M".to_string(), + points: test_geometry.clone(), + densities: vec![0.03, 0.03], + simulation_boundary_points: vec![(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)], + electronic_stopping_correction_factor: 1.0, + }; + + let material_parameters = material::MaterialParameters{ + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![0.0, 0.0], + Es: vec![2.0, 4.0], + Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], + Z: vec![29., 1.], + m: vec![63.54, 1.0008], + interaction_index: vec![0, 0], + surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET}, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); + + let query_points_outside = vec![[0.5, 0.1, 0.0], [0.9, 0.5, 0.0], [0.6, 0.7, 0.0], [0.4, 0.7, 0.0], [0.1, 0.5, 0.0]]; + let query_points_inside = vec![[0.5, 0.21, 0.0], [0.75, 0.21, 0.0], [0.51, 0.5, 0.0], [0.49, 0.5, 0.0], [0.25, 0.21, 0.0]]; + + for query_point in query_points_outside.clone() { + assert!(!material_homogeneous_2D.inside(query_point[0], query_point[1], query_point[2]), "Point ({}, {}, {}) failed outsidedness check.", query_point[0], query_point[1], query_point[2]); + } + + for query_point in query_points_inside.clone() { + assert!(material_homogeneous_2D.inside(query_point[0], query_point[1], query_point[2]), "Point ({}, {}, {}) failed outsidedness check.", query_point[0], query_point[1], query_point[2]); + } + + + let normals_outside = query_points_outside + .iter() + .map( |query_point| { + (material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2])) + } + ).collect::>(); + + let normals_inside = query_points_inside + .iter() + .map( |query_point| { + (material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2])) + } + ).collect::>(); + + let test_normals = (0..num_segments) + .zip(1..num_segments + 1) + .map(|(i, j)| { + let dx = test_geometry[i % num_segments].0 - test_geometry[j % num_segments].0; + let dy = test_geometry[i % num_segments].1 - test_geometry[j % num_segments].1; + let mag = (dx*dx + dy*dy).sqrt(); + (-dy/mag, dx/mag, 0.0) + }).collect::>(); + + for (test_normal, normal) in test_normals.iter().zip(normals_outside) { + assert!(approx_eq!(f64, normal.0, test_normal.0, epsilon=1E-9), "Normal vector check failed. Calculated Normal: ({}, {}, {}) Test: ({}, {}, {})", normal.0, normal.1, normal.2, test_normal.0, test_normal.1, test_normal.2); + assert!(approx_eq!(f64, normal.1, test_normal.1, epsilon=1E-9)); + assert!(approx_eq!(f64, normal.2, test_normal.2, epsilon=1E-9)); + } + + for (test_normal, normal) in test_normals.iter().zip(normals_inside) { + assert!(approx_eq!(f64, normal.0, test_normal.0, epsilon=1E-9), "Normal vector check failed. Calculated Normal: ({}, {}, {}) Test: ({}, {}, {})", normal.0, normal.1, normal.2, test_normal.0, test_normal.1, test_normal.2); + assert!(approx_eq!(f64, normal.1, test_normal.1, epsilon=1E-9)); + assert!(approx_eq!(f64, normal.2, test_normal.2, epsilon=1E-9)); + } + +} + #[test] #[cfg(feature = "cpr_rootfinder")] @@ -439,6 +554,68 @@ fn test_spherical_geometry() { } +#[test] +fn extended_test_2D_geometry() { + let mass = 1.; + let Z = 1.; + let E = 10.*EV; + let Ec = 1.*EV; + let Es = 5.76*EV; + let x = 1.5e-6; + let y = 10.0e-6; + let z = 0.; + let cosx = 0.0; + let cosy = -1.0; + let cosz = 0.; + + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); + + let material_parameters = material::MaterialParameters{ + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![0.0, 0.0], + Es: vec![2.0, 4.0], + Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], + Z: vec![29., 1.], + m: vec![63.54, 1.0008], + interaction_index: vec![0, 0], + surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET}, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let thickness: f64 = 1000.; + let depth: f64 = 1000.; + + let geometry_input_2D = geometry::Mesh2DInput { + length_unit: "MICRON".to_string(), + points: vec![(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (4.0, 1.0), (3.0, 1.0), (3.0, 10.0), (2.0, 10.0), (2.0, 1.0), (1.0, 1.0), (1.0, 10.0), (0.0, 10.0), (0.0, 1.0)], + triangles: vec![(13, 12, 11), (13, 11, 10), (0, 13, 10), (0, 10, 1), (1, 10, 9), (1, 9, 2), (9, 8, 7), (9, 7, 6), (2, 9, 6), (2, 6, 3), (3, 6, 5), (3, 5, 4)], + densities: vec![vec![3e10, 3e10]; 12], + boundary: vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 0], + simulation_boundary_points: vec![(-0.1, -0.1), (4.1, -0.1), (4.1, 10.1), (-0.1, 10.1), (-0.1, -0.1)], + energy_barrier_thickness: 4e-4, + electronic_stopping_correction_factors: vec![1.0; 12], + }; + + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + + if let Closest::Intersection(p) = material_2D.geometry.boundary.closest_point(&point!(x: 1.5e-6, y: 1.0001e-6)) { + println!("({} {})", p.x(), p.y()); + } else if let Closest::SinglePoint(p) = material_2D.geometry.boundary.closest_point(&point!(x: 1.5e-6, y: 1.0001e-6)) { + println!("({} {})", p.x(), p.y()); + } else { + panic!(); + } + + assert!(material_2D.inside(1.5e-6, 0.999e-6, 0.0)); + assert!(!material_2D.inside(1.5e-6, 1.0001e-6, 0.0)); + + material::surface_binding_energy(&mut particle_1, &material_2D); + + material::boundary_condition_planar(&mut particle_1, &material_2D); +} + #[test] fn test_geometry() { let mass = 1.; @@ -818,7 +995,7 @@ fn test_momentum_conservation() { for scattering_integral in vec![ScatteringIntegral::MENDENHALL_WELLER, ScatteringIntegral::GAUSS_MEHLER{n_points: 10}, ScatteringIntegral::GAUSS_LEGENDRE] { for root_finder in vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}] { - println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); + //println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); @@ -883,9 +1060,11 @@ fn test_momentum_conservation() { let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options); + /* println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal, binary_collision_geometries[0].impact_parameter/ANGSTROM, binary_collision_geometries[0].mfp/ANGSTROM); + */ let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1, &binary_collision_geometries[0], &options); @@ -898,12 +1077,13 @@ fn test_momentum_conservation() { let binary_collision_result = bca::calculate_binary_collision(&particle_1, &particle_2, &binary_collision_geometries[0], &options).unwrap(); + /* println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV, binary_collision_result.psi, binary_collision_result.psi_recoil); println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); - + */ //Energy transfer to recoil particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); @@ -924,11 +1104,13 @@ fn test_momentum_conservation() { let final_momentum = mom1_1.add(&mom2_1); + /* println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude()); println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude()); println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); println!(); + */ assert!(approx_eq!(f64, initial_momentum.x, final_momentum.x, epsilon = 1E-12)); assert!(approx_eq!(f64, initial_momentum.y, final_momentum.y, epsilon = 1E-12)); From b5c892ef94633d5c2d4e7fa8448e726b40436933 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 17 Jul 2024 10:32:02 -0700 Subject: [PATCH 02/13] Prelim draft of replacement of geo with parry2d. --- src/main.rs | 6 +++--- src/tests.rs | 45 +++++++++++++++++++++++++-------------------- 2 files changed, 28 insertions(+), 23 deletions(-) diff --git a/src/main.rs b/src/main.rs index c371b7ca..d29f6555 100644 --- a/src/main.rs +++ b/src/main.rs @@ -58,7 +58,7 @@ pub use crate::consts::*; pub use crate::structs::*; 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, HomogeneousMesh2D, ParryHomogeneousMesh2D}; +pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, ParryMesh2D, HomogeneousMesh2D, ParryHomogeneousMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; pub use crate::physics::{physics_loop}; @@ -97,8 +97,8 @@ fn main() { physics_loop::(particle_input_array, material, options, output_units); }, GeometryType::MESH2D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); + let (particle_input_array, material, options, output_units) = input::input::(input_file); + physics_loop::(particle_input_array, material, options, output_units); }, GeometryType::SPHERE => { let (particle_input_array, material, options, output_units) = input::input::(input_file); diff --git a/src/tests.rs b/src/tests.rs index d70e31c6..4712e75e 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -592,21 +592,14 @@ fn extended_test_2D_geometry() { points: vec![(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (4.0, 1.0), (3.0, 1.0), (3.0, 10.0), (2.0, 10.0), (2.0, 1.0), (1.0, 1.0), (1.0, 10.0), (0.0, 10.0), (0.0, 1.0)], triangles: vec![(13, 12, 11), (13, 11, 10), (0, 13, 10), (0, 10, 1), (1, 10, 9), (1, 9, 2), (9, 8, 7), (9, 7, 6), (2, 9, 6), (2, 6, 3), (3, 6, 5), (3, 5, 4)], densities: vec![vec![3e10, 3e10]; 12], - boundary: vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 0], - simulation_boundary_points: vec![(-0.1, -0.1), (4.1, -0.1), (4.1, 10.1), (-0.1, 10.1), (-0.1, -0.1)], + boundary: vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], + simulation_boundary_points: vec![(-0.1, -0.1), (4.1, -0.1), (4.1, 10.1), (-0.1, 10.1)], energy_barrier_thickness: 4e-4, electronic_stopping_correction_factors: vec![1.0; 12], }; - let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); - if let Closest::Intersection(p) = material_2D.geometry.boundary.closest_point(&point!(x: 1.5e-6, y: 1.0001e-6)) { - println!("({} {})", p.x(), p.y()); - } else if let Closest::SinglePoint(p) = material_2D.geometry.boundary.closest_point(&point!(x: 1.5e-6, y: 1.0001e-6)) { - println!("({} {})", p.x(), p.y()); - } else { - panic!(); - } assert!(material_2D.inside(1.5e-6, 0.999e-6, 0.0)); assert!(!material_2D.inside(1.5e-6, 1.0001e-6, 0.0)); @@ -654,7 +647,7 @@ fn test_geometry() { triangles: vec![(0, 1, 2), (0, 2, 3)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], boundary: vec![0, 1, 2, 3], - 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.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], energy_barrier_thickness: 10., electronic_stopping_correction_factors: vec![1.0, 1.0], }; @@ -663,7 +656,7 @@ fn test_geometry() { 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.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], electronic_stopping_correction_factor: 1.0, }; @@ -681,7 +674,7 @@ fn test_geometry() { electronic_stopping_correction_factor: 1.0 }; - let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); let material_1D: material::Material = material::Material::::new(&material_parameters, &geometry_input_1D); let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); @@ -706,7 +699,12 @@ fn test_geometry() { 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.)); + let (x2d, y2d, z2d) = material_2D.geometry.closest_point(0., 0., 0.); + let (x0d, y0d, z0d) = material_0D.geometry.closest_point(0., 0., 0.); + assert!(approx_eq!(f64, x2d, x0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, y2d, y0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, z2d, z0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + 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.)); @@ -714,7 +712,14 @@ fn test_geometry() { 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.)); + + let (x2d, y2d, z2d) = material_2D.geometry.closest_point(-10., 0., 5.); + let (x0d, y0d, z0d) = material_0D.geometry.closest_point(-10., 0., 5.); + assert!(approx_eq!(f64, x2d, x0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, y2d, y0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, z2d, z0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + + 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.)); } @@ -757,7 +762,7 @@ fn test_surface_binding_energy_barrier() { triangles: vec![(0, 1, 2), (0, 2, 3)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], boundary: vec![0, 1, 2, 3], - 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.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], energy_barrier_thickness: 10., electronic_stopping_correction_factors: vec![1.0, 1.0], }; @@ -768,7 +773,7 @@ fn test_surface_binding_energy_barrier() { electronic_stopping_correction_factor: 1.0 }; - let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; @@ -982,13 +987,13 @@ fn test_momentum_conservation() { triangles: vec![(0, 1, 2), (0, 2, 3)], points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], densities: vec![vec![0.06306], vec![0.06306]], - boundary: vec![0, 1, 2, 3], - 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.)], + boundary: vec![3, 2, 1, 0], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], electronic_stopping_correction_factors: vec![0.0, 0.0], energy_barrier_thickness: 0., }; - let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); + let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); for high_energy_free_flight_paths in vec![true, false] { for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { From 782fb04224b85366af2f318dc86676f0caaf2cb6 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 17 Jul 2024 10:46:26 -0700 Subject: [PATCH 03/13] Add updated Cargo.toml. --- Cargo.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index f4738f10..566c7236 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,13 +21,14 @@ toml = "0.7.4" anyhow = "1.0.71" itertools = "0.10.5" rayon = "1.7.0" -geo = {version = "0.25", optional = false} +geo = {version = "0.28", optional = false} indicatif = {version = "0.15.0", features=["rayon"]} serde = { version = "1.0.163", features = ["derive"] } hdf5 = {version = "0.7.1", optional = true} rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} ndarray = {version = "0.14.0", features = ["serde"], optional = true} parry3d-f64 = {optional = true, version="0.2.0"} +parry2d-f64 = {optional = false, version="*"} [dependencies.pyo3] version = "0.19.0" From d8e1109a349b3aeca07f3c8a4d9c17eee6623c24 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 17 Jul 2024 10:57:13 -0700 Subject: [PATCH 04/13] Fix to simulation boundary of layered geometry 2d example. --- examples/layered_geometry.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index 235cc5a6..0e48ddb6 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -54,7 +54,7 @@ points = [[0.0, -0.5], [0.01, -0.5], [0.04, -0.5], [0.5, -0.5], [0.5, 0.5], [0.0 triangles = [[0, 1, 6], [0, 6, 7], [1, 2, 5], [1, 5, 6], [2, 3, 4], [2, 4, 5]] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] boundary = [0, 3, 4, 7] -simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],] +simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],] electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] energy_barrier_thickness = 2.2E-4 From 9d0c1b371da8578cc9db4c12a20a7d6356bfff4f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 17 Jul 2024 10:57:56 -0700 Subject: [PATCH 05/13] Fix to multiple interaction potentials example. --- examples/multiple_interaction_potentials.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index f2194ab1..736b39d3 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -48,7 +48,7 @@ length_unit = "MICRON" triangles = [ [ 0.0, 0.01, 0.0, 0.5, 0.5, -0.5,], [ 0.0, 0.01, 0.01, -0.5, 0.5, -0.5,], [ 0.01, 0.01, 0.04, -0.5, 0.5, -0.5,], [ 0.01, 0.04, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.5, -0.5, 0.5, -0.5,],] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] material_boundary_points = [ [ 0.0, 0.5,], [ 0.5, 0.5,], [ 0.5, -0.5,], [ 0.0, -0.5,],] -simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],] +simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],] electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [material_parameters] From 3bd38c962930e792cb3341ed87b6b3d658520802 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 18 Jul 2024 11:27:00 -0700 Subject: [PATCH 06/13] Further fixes to geometry. --- examples/layered_geometry.toml | 5 +++-- src/geometry.rs | 6 ++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index 0e48ddb6..078c7993 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -1,7 +1,7 @@ #Use with 2D geometry option only [options] name = "2000.0eV_0.0001deg_He_TiO2_Al_Si" -track_trajectories = false +track_trajectories = true track_recoils = true track_recoil_trajectories = false write_buffer_size = 8000 @@ -54,7 +54,8 @@ points = [[0.0, -0.5], [0.01, -0.5], [0.04, -0.5], [0.5, -0.5], [0.5, 0.5], [0.0 triangles = [[0, 1, 6], [0, 6, 7], [1, 2, 5], [1, 5, 6], [2, 3, 4], [2, 4, 5]] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] boundary = [0, 3, 4, 7] -simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],] +#simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],] +simulation_boundary_points = [[0.6, 0.6], [-0.1, 0.6], [-0.1, -0.6], [0.6, -0.6]] electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] energy_barrier_thickness = 2.2E-4 diff --git a/src/geometry.rs b/src/geometry.rs index 8f41ebe3..b4c91cd1 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -719,10 +719,16 @@ impl Geometry for ParryMesh2D { let densities: Vec> = geometry_input.densities.iter().map(|density_list| density_list.iter().map(|density| density / (length_unit).powi(3)).collect::>()).collect::>>(); + dbg!(&densities); + let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit; let concentrations: Vec> = densities.iter().map(|density_list| density_list.iter().map(|density| density / density_list.iter().sum::() ).collect::>()).collect::>>(); + for concentration in &concentrations { + assert!((concentration.iter().sum::() - 1.0).abs() < 1e-6); + } + let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); linked_boundary_points.push([number_boundary_points - 1, 0]); let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_boundary_points)); From 5393183db7884891e1a740b7ad34ac9be2a0c0aa Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 19 Jul 2024 13:13:16 -0700 Subject: [PATCH 07/13] Updates workflow to include comments on cat-ed output files. --- .github/workflows/rustbca_compile_check.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 584476c3..382c9b8a 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -59,12 +59,12 @@ jobs: cargo run --release 0D examples/boron_nitride_0D.toml ./target/release/RustBCA 0D examples/titanium_dioxide_0D.toml ./target/release/RustBCA 1D examples/layered_geometry_1D.toml - cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output + echo "layered geometry 1D:"; cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output ./target/release/RustBCA examples/layered_geometry.toml - cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output + echo "layered geometry 2D:"; 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 ./target/release/RustBCA examples/boron_nitride_wire.toml - cat boron_nitride_summary.output + echo "boron nitride 2D:"; cat boron_nitride_summary.output ./target/release/RustBCA HOMOGENEOUS2D examples/boron_nitride_wire_homogeneous.toml - cat boron_nitride_h_summary.output + echo "boron nitride 2D Homogeneous:"; cat boron_nitride_h_summary.output From 5acfcd0439c86880c49f70d67261f12e1b7945f8 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 22 Jul 2024 11:17:54 -0700 Subject: [PATCH 08/13] Deprecation of geo calls. --- Cargo.toml | 5 +-- src/enums.rs | 4 +-- src/geometry.rs | 81 ++++++++++++++++++++++++++++++++----------------- src/lib.rs | 2 +- src/main.rs | 4 +-- src/tests.rs | 9 +++--- 6 files changed, 66 insertions(+), 39 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 566c7236..e05b0710 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,6 +5,8 @@ default-run = "RustBCA" authors = ["Jon Drobny ", "Jon Drobny "] edition = "2021" + + [[bin]] name = "RustBCA" path = "src/main.rs" @@ -21,7 +23,6 @@ toml = "0.7.4" anyhow = "1.0.71" itertools = "0.10.5" rayon = "1.7.0" -geo = {version = "0.28", optional = false} indicatif = {version = "0.15.0", features=["rayon"]} serde = { version = "1.0.163", features = ["derive"] } hdf5 = {version = "0.7.1", optional = true} @@ -42,7 +43,7 @@ float-cmp = "0.8.0" lto = "fat" codegen-units = 1 opt-level = 3 -debug = false +debug = 1 [features] hdf5_input = ["hdf5"] diff --git a/src/enums.rs b/src/enums.rs index 908c0409..d8582b6b 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -3,13 +3,13 @@ use super::*; pub enum MaterialType { MESH0D(material::Material), MESH1D(material::Material), - MESH2D(material::Material), + MESH2D(material::Material), SPHERE(material::Material), #[cfg(feature = "parry3d")] BALL(material::Material), #[cfg(feature = "parry3d")] TRIMESH(material::Material), - HOMOGENEOUS2D(material::Material) + HOMOGENEOUS2D(material::Material) } #[derive(Deserialize, PartialEq, Clone, Copy)] diff --git a/src/geometry.rs b/src/geometry.rs index b4c91cd1..55affcec 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -6,13 +6,9 @@ use parry2d_f64::shape::{Polyline, TriMesh, FeatureId::*}; use parry2d_f64::shape::SegmentPointLocation::{OnEdge, OnVertex}; use parry2d_f64::query::{PointQuery, Ray, RayCast}; use parry2d_f64::math::{Isometry, Vector}; -use parry2d_f64::math::Point as Point2; +use parry2d_f64::math::Point as Point2d; use parry2d_f64::bounding_volume::aabb; -use geo::algorithm::contains::Contains; -use geo::{Polygon, LineString, Point, point, Closest}; -use geo::algorithm::closest_point::ClosestPoint; - ///Trait for a Geometry object - all forms of geometry must implement these traits to be used pub trait Geometry { @@ -322,15 +318,23 @@ impl Geometry for ParryHomogeneousMesh2D { ).as_str()), }; - let boundary_points_converted2: Vec> = input.points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + let mut boundary_points_converted2: Vec> = input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); let number_boundary_points = boundary_points_converted2.len() as u32; for p in boundary_points_converted2.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") } - let simulation_boundary_points_converted2: Vec> = input.simulation_boundary_points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + let mut simulation_boundary_points_converted2: Vec> = input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_simulation_boundary_points = simulation_boundary_points_converted2.len() as u32; + + let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize) + .map(|i| (simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted2[i].x)*(simulation_boundary_points_converted2[i].y + simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].y)) + .sum::(); + if test_simulation_boundary_ccw > 0.0 { + simulation_boundary_points_converted2 = simulation_boundary_points_converted2.clone().into_iter().rev().collect::>>(); + } for p in simulation_boundary_points_converted2.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") @@ -340,7 +344,9 @@ impl Geometry for ParryHomogeneousMesh2D { .map(|i| (boundary_points_converted2[(i + 1) % number_boundary_points as usize].x - boundary_points_converted2[i].x)*(boundary_points_converted2[i].y + boundary_points_converted2[(i + 1) % number_boundary_points as usize].y)) .sum::(); - assert!(test_ccw <= 0.0, "Input error: signed area of boundary is less than or equal to zero"); + if test_ccw > 0.0 { + boundary_points_converted2 = boundary_points_converted2.clone().into_iter().rev().collect::>>(); + } let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; @@ -358,7 +364,6 @@ impl Geometry for ParryHomogeneousMesh2D { linked_boundary_points.push([number_boundary_points - 1, 0]); let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_boundary_points)); - let number_simulation_boundary_points = simulation_boundary_points_converted2.clone().len() as u32; let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted2, Some(linked_simulation_boundary_points)); @@ -390,13 +395,13 @@ impl Geometry for ParryHomogeneousMesh2D { } fn inside(&self, x: f64, y: f64, z: f64) -> bool { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); point_projection.is_inside } fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p); point_projection.is_inside @@ -406,20 +411,20 @@ impl Geometry for ParryHomogeneousMesh2D { if self.inside(x, y, z) { true } else { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); (self.boundary.distance_to_local_point(&p, true) as f64) < self.energy_barrier_thickness } } fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); let (x_, y_) = (point_projection.point.x, point_projection.point.y); (x_ as f64, y_ as f64, z) } fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z); @@ -446,6 +451,7 @@ pub struct HomogeneousMesh2DInput { pub electronic_stopping_correction_factor: f64 } +/* #[derive(Clone)] pub struct HomogeneousMesh2D { pub boundary: Polygon, @@ -554,6 +560,7 @@ impl Geometry for HomogeneousMesh2D { panic!("Not implemented.") } } +*/ /// Object that contains raw mesh input data. #[derive(Deserialize, Clone)] @@ -581,7 +588,7 @@ pub struct ParryMesh2D { impl ParryMesh2D { fn get_id(&self, x: f64, y: f64, z: f64) -> usize { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, feature_id) = self.trimesh.project_local_point_and_get_feature(&p); match feature_id { Vertex(vertex_id) => { @@ -609,7 +616,7 @@ impl Geometry for ParryMesh2D { if self.inside(x, y, z) { true } else { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); if self.trimesh.distance_to_local_point(&p, true) < self.energy_barrier_thickness { true } else { @@ -619,7 +626,7 @@ impl Geometry for ParryMesh2D { } fn inside_simulation_boundary (&self, x: f64, y: f64, z: f64) -> bool { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p); @@ -627,7 +634,7 @@ impl Geometry for ParryMesh2D { } fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); let (x_, y_) = (point_projection.point.x, point_projection.point.y); (x_ as f64, y_ as f64, z) @@ -653,13 +660,13 @@ impl Geometry for ParryMesh2D { /// Determines whether the point (x, y) is inside the mesh. fn inside(&self, x: f64, y: f64, z: f64) -> bool { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); point_projection.is_inside } fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { - let p = Point2::new(x, y); + let p = Point2d::new(x, y); let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z); @@ -691,14 +698,25 @@ impl Geometry for ParryMesh2D { ).as_str()), }; - let boundary_points_converted2: Vec> = geometry_input.points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + let mut boundary_points_converted2: Vec> = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); let number_boundary_points = boundary_points_converted2.len() as u32; - for p in boundary_points_converted2.iter().combinations(2) { - assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + let mut simulation_boundary_points_converted2: Vec> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_simulation_boundary_points = simulation_boundary_points_converted2.len() as u32; + + let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize) + .map(|i| (simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted2[i].x)*(simulation_boundary_points_converted2[i].y + simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].y)) + .sum::(); + + dbg!(test_simulation_boundary_ccw > 0.0); + + dbg!(&simulation_boundary_points_converted2); + + if test_simulation_boundary_ccw > 0.0 { + simulation_boundary_points_converted2 = simulation_boundary_points_converted2.clone().into_iter().rev().collect::>>(); } - let simulation_boundary_points_converted2: Vec> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect(); + dbg!(&simulation_boundary_points_converted2); for p in simulation_boundary_points_converted2.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") @@ -708,10 +726,16 @@ impl Geometry for ParryMesh2D { .map(|i| (boundary_points_converted2[(i + 1) % number_boundary_points as usize].x - boundary_points_converted2[i].x)*(boundary_points_converted2[i].y + boundary_points_converted2[(i + 1) % number_boundary_points as usize].y)) .sum::(); - assert!(test_ccw <= 0.0, "Input error: signed area of boundary is less than or equal to zero"); + if test_ccw > 0.0 { + boundary_points_converted2 = boundary_points_converted2.clone().into_iter().rev().collect::>>(); + } + + for p in boundary_points_converted2.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } let triangles = geometry_input.triangles.iter().map(|(i, j, k)| [*i as u32, *j as u32, *k as u32]).collect::>(); - let points_converted = geometry_input.points.iter().map(|(x, y)| Point2::new(x*length_unit, y*length_unit)).collect::>>(); + let points_converted = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect::>>(); let trimesh = TriMesh::new(points_converted, triangles); @@ -719,7 +743,6 @@ impl Geometry for ParryMesh2D { let densities: Vec> = geometry_input.densities.iter().map(|density_list| density_list.iter().map(|density| density / (length_unit).powi(3)).collect::>()).collect::>>(); - dbg!(&densities); let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit; @@ -750,6 +773,7 @@ impl Geometry for ParryMesh2D { } } +/* /// Triangular mesh for rustbca. #[derive(Clone)] pub struct Mesh2D { @@ -1012,6 +1036,7 @@ impl GeometryElement for Cell2D { self.electronic_stopping_correction_factor } } +*/ #[derive(Clone)] pub struct Layer1D { @@ -1061,6 +1086,7 @@ impl GeometryElement for Layer1D { } } +/* /// A triangle in 2D, with points (x1, y1), (x2, y2), (x3, y3), and the three line segments bewtween them. #[derive(Clone)] pub struct Triangle2D { @@ -1151,3 +1177,4 @@ impl Triangle2D { ((x - centroid.0)*(x - centroid.0) + (y - centroid.1)*(y - centroid.1)).sqrt() } } +*/ \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index e08c6793..2ba8f5d4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -75,7 +75,7 @@ pub use crate::consts::*; pub use crate::structs::*; 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, ParryMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; #[cfg(feature = "parry3d")] diff --git a/src/main.rs b/src/main.rs index d29f6555..a66f4d7d 100644 --- a/src/main.rs +++ b/src/main.rs @@ -58,7 +58,7 @@ pub use crate::consts::*; pub use crate::structs::*; 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, ParryMesh2D, HomogeneousMesh2D, ParryHomogeneousMesh2D}; +pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, ParryMesh2D, ParryHomogeneousMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; pub use crate::physics::{physics_loop}; @@ -84,7 +84,7 @@ fn main() { "HOMOGENEOUS2D" => GeometryType::HOMOGENEOUS2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), - _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), + _ => panic!("Too many command line arguments. RustBCA accepts 0 (defaulting to 'input.toml' and 2D mesh), 1 ( defaulting to 2D mesh), or 2 ( )"), }; match geometry_type { diff --git a/src/tests.rs b/src/tests.rs index 4712e75e..4e26124b 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -3,10 +3,6 @@ use super::*; #[cfg(test)] use float_cmp::*; -use geo::algorithm::contains::Contains; -use geo::{Polygon, LineString, Point, point, Closest}; -use geo::algorithm::closest_point::ClosestPoint; - use parry2d_f64::shape::Polyline; use parry2d_f64::query::{PointQuery, PointProjection, Ray, RayCast}; use parry2d_f64::math::{Isometry}; @@ -675,7 +671,7 @@ fn test_geometry() { }; let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); - let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); + let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); let material_1D: material::Material = material::Material::::new(&material_parameters, &geometry_input_1D); let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; @@ -844,6 +840,8 @@ fn test_surface_binding_energy_barrier() { assert!(material_0D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); } +//Deprecated with geo +/* #[test] fn test_triangle_contains() { let triangle_1 = geometry::Triangle2D::new((0., 2., 0., 2., 0., 0.)); @@ -867,6 +865,7 @@ fn test_triangle_distance_to() { assert!(approx_eq!(f64, triangle_1.distance_to(2., 0.), 0., epsilon=1E-12), "{}", triangle_1.distance_to(2., 0.)); assert!(approx_eq!(f64, triangle_1.distance_to(0., 2.), 0., epsilon=1E-12), "{}", triangle_1.distance_to(0., 2.)); } +*/ #[test] fn test_surface_refraction() { From e5e8aa4a1eea42b77d2e3a18c433356213da66fb Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 22 Jul 2024 13:30:33 -0700 Subject: [PATCH 09/13] Cleanup of new 2D geometry options. --- src/geometry.rs | 78 ++++++++++++++++++++++++++----------------------- src/input.rs | 5 ---- src/lib.rs | 1 + src/tests.rs | 10 +++---- 4 files changed, 46 insertions(+), 48 deletions(-) diff --git a/src/geometry.rs b/src/geometry.rs index 55affcec..81f9f33e 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -3,11 +3,9 @@ use super::*; use itertools::Itertools; use parry2d_f64::shape::{Polyline, TriMesh, FeatureId::*}; -use parry2d_f64::shape::SegmentPointLocation::{OnEdge, OnVertex}; -use parry2d_f64::query::{PointQuery, Ray, RayCast}; -use parry2d_f64::math::{Isometry, Vector}; +use parry2d_f64::query::PointQuery; use parry2d_f64::math::Point as Point2d; -use parry2d_f64::bounding_volume::aabb; +use parry2d_f64::math::Isometry; ///Trait for a Geometry object - all forms of geometry must implement these traits to be used pub trait Geometry { @@ -318,34 +316,34 @@ impl Geometry for ParryHomogeneousMesh2D { ).as_str()), }; - let mut boundary_points_converted2: Vec> = input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); - let number_boundary_points = boundary_points_converted2.len() as u32; + let mut boundary_points_converted: Vec> = input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_boundary_points = boundary_points_converted.len() as u32; - for p in boundary_points_converted2.iter().combinations(2) { + for p in boundary_points_converted.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") } - let mut simulation_boundary_points_converted2: Vec> = input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); - let number_simulation_boundary_points = simulation_boundary_points_converted2.len() as u32; + let mut simulation_boundary_points_converted: Vec> = input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_simulation_boundary_points = simulation_boundary_points_converted.len() as u32; let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize) - .map(|i| (simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted2[i].x)*(simulation_boundary_points_converted2[i].y + simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].y)) + .map(|i| (simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted[i].x)*(simulation_boundary_points_converted[i].y + simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].y)) .sum::(); if test_simulation_boundary_ccw > 0.0 { - simulation_boundary_points_converted2 = simulation_boundary_points_converted2.clone().into_iter().rev().collect::>>(); + simulation_boundary_points_converted = simulation_boundary_points_converted.clone().into_iter().rev().collect::>>(); } - for p in simulation_boundary_points_converted2.iter().combinations(2) { + for p in simulation_boundary_points_converted.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") } let test_ccw = (0..number_boundary_points as usize) - .map(|i| (boundary_points_converted2[(i + 1) % number_boundary_points as usize].x - boundary_points_converted2[i].x)*(boundary_points_converted2[i].y + boundary_points_converted2[(i + 1) % number_boundary_points as usize].y)) + .map(|i| (boundary_points_converted[(i + 1) % number_boundary_points as usize].x - boundary_points_converted[i].x)*(boundary_points_converted[i].y + boundary_points_converted[(i + 1) % number_boundary_points as usize].y)) .sum::(); if test_ccw > 0.0 { - boundary_points_converted2 = boundary_points_converted2.clone().into_iter().rev().collect::>>(); + boundary_points_converted = boundary_points_converted.clone().into_iter().rev().collect::>>(); } let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; @@ -358,15 +356,13 @@ impl Geometry for ParryHomogeneousMesh2D { let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); - - let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); linked_boundary_points.push([number_boundary_points - 1, 0]); - let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_boundary_points)); + let boundary2 = Polyline::new(boundary_points_converted, Some(linked_boundary_points)); let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); - let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted2, Some(linked_simulation_boundary_points)); + let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points)); ParryHomogeneousMesh2D { densities, @@ -396,8 +392,12 @@ impl Geometry for ParryHomogeneousMesh2D { fn inside(&self, x: f64, y: f64, z: f64) -> bool { let p = Point2d::new(x, y); - let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); - point_projection.is_inside + if self.boundary.aabb(&Isometry::identity()).contains_local_point(&p) { + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + } else { + false + } } fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { @@ -661,8 +661,12 @@ impl Geometry for ParryMesh2D { /// Determines whether the point (x, y) is inside the mesh. fn inside(&self, x: f64, y: f64, z: f64) -> bool { let p = Point2d::new(x, y); - let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); - point_projection.is_inside + if self.boundary.aabb(&Isometry::identity()).contains_local_point(&p) { + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + } else { + false + } } fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { @@ -698,39 +702,39 @@ impl Geometry for ParryMesh2D { ).as_str()), }; - let mut boundary_points_converted2: Vec> = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); - let number_boundary_points = boundary_points_converted2.len() as u32; + let mut boundary_points_converted: Vec> = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_boundary_points = boundary_points_converted.len() as u32; - let mut simulation_boundary_points_converted2: Vec> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); - let number_simulation_boundary_points = simulation_boundary_points_converted2.len() as u32; + let mut simulation_boundary_points_converted: Vec> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_simulation_boundary_points = simulation_boundary_points_converted.len() as u32; let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize) - .map(|i| (simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted2[i].x)*(simulation_boundary_points_converted2[i].y + simulation_boundary_points_converted2[(i + 1) % number_simulation_boundary_points as usize].y)) + .map(|i| (simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted[i].x)*(simulation_boundary_points_converted[i].y + simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].y)) .sum::(); dbg!(test_simulation_boundary_ccw > 0.0); - dbg!(&simulation_boundary_points_converted2); + dbg!(&simulation_boundary_points_converted); if test_simulation_boundary_ccw > 0.0 { - simulation_boundary_points_converted2 = simulation_boundary_points_converted2.clone().into_iter().rev().collect::>>(); + simulation_boundary_points_converted = simulation_boundary_points_converted.clone().into_iter().rev().collect::>>(); } - dbg!(&simulation_boundary_points_converted2); + dbg!(&simulation_boundary_points_converted); - for p in simulation_boundary_points_converted2.iter().combinations(2) { + for p in simulation_boundary_points_converted.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") } let test_ccw = (0..number_boundary_points as usize) - .map(|i| (boundary_points_converted2[(i + 1) % number_boundary_points as usize].x - boundary_points_converted2[i].x)*(boundary_points_converted2[i].y + boundary_points_converted2[(i + 1) % number_boundary_points as usize].y)) + .map(|i| (boundary_points_converted[(i + 1) % number_boundary_points as usize].x - boundary_points_converted[i].x)*(boundary_points_converted[i].y + boundary_points_converted[(i + 1) % number_boundary_points as usize].y)) .sum::(); if test_ccw > 0.0 { - boundary_points_converted2 = boundary_points_converted2.clone().into_iter().rev().collect::>>(); + boundary_points_converted = boundary_points_converted.clone().into_iter().rev().collect::>>(); } - for p in boundary_points_converted2.iter().combinations(2) { + for p in boundary_points_converted.iter().combinations(2) { assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") } @@ -754,12 +758,12 @@ impl Geometry for ParryMesh2D { let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); linked_boundary_points.push([number_boundary_points - 1, 0]); - let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_boundary_points)); + let boundary2 = Polyline::new(boundary_points_converted, Some(linked_boundary_points)); - let number_simulation_boundary_points = simulation_boundary_points_converted2.clone().len() as u32; + let number_simulation_boundary_points = simulation_boundary_points_converted.clone().len() as u32; let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); - let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted2, Some(linked_simulation_boundary_points)); + let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points)); ParryMesh2D { trimesh, diff --git a/src/input.rs b/src/input.rs index 8e2b81dd..af04e1f9 100644 --- a/src/input.rs +++ b/src/input.rs @@ -172,11 +172,6 @@ fn one_u64() -> u64 { 1 } -///This helper function is a workaround to issue #368 in serde -fn three() -> usize { - 3 -} - fn zero_usize() -> usize{ 0 } diff --git a/src/lib.rs b/src/lib.rs index 2ba8f5d4..dd255697 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -56,6 +56,7 @@ use pyo3::types::*; //Load internal modules pub mod material; pub mod particle; +#[cfg(test)] pub mod tests; pub mod interactions; pub mod bca; diff --git a/src/tests.rs b/src/tests.rs index 4e26124b..5bd32003 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -4,10 +4,8 @@ use super::*; use float_cmp::*; use parry2d_f64::shape::Polyline; -use parry2d_f64::query::{PointQuery, PointProjection, Ray, RayCast}; -use parry2d_f64::math::{Isometry}; +use parry2d_f64::math::Isometry; use parry2d_f64::math::Point as Point2d; -use parry2d_f64::bounding_volume::aabb; /* let number_boundary_points = boundary_points_converted.clone().len() as u32; @@ -21,7 +19,7 @@ fn test_parry2d() { let points: Vec> = vec![Point2d::new(-1.0, -1.0), Point2d::new(-1.0, 1.0), Point2d::new(1.0, 1.0), Point2d::new(1.0, -1.0)]; //let mut indices: Vec<[u32; 2]> = vec![]; //let mut indices: Vec<[u32; 2]> = vec![[0, 1], [1, 2], [2, 3], [3, 0]]; - let mut indices: Vec<[u32; 2]> = vec![[0, 3], [3, 2], [2, 1], [1, 0]]; + let indices: Vec<[u32; 2]> = vec![[0, 3], [3, 2], [2, 1], [1, 0]]; //indices.reverse(); //(0.0000006261797114005236, 0.000006009591179670447) @@ -80,14 +78,14 @@ fn test_parry2d() { let normals_outside = query_points_outside .iter() .map( |query_point| { - (material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2])) + material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2]) } ).collect::>(); let normals_inside = query_points_inside .iter() .map( |query_point| { - (material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2])) + material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2]) } ).collect::>(); From d8ef8f9cacfc0639b262ebcdb63ca54155551b1e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 7 Feb 2025 15:29:12 -0800 Subject: [PATCH 10/13] Fix expect-unwrap error in trimesh creation. --- src/geometry.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geometry.rs b/src/geometry.rs index 814f7e52..dae8e687 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -768,7 +768,7 @@ impl Geometry for ParryMesh2D { let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points)); ParryMesh2D { - trimesh, + trimesh.expect(""), densities, simulation_boundary: simulation_boundary2, boundary: boundary2, From 54e7286e04cbfcf842307be3680636191c1018d2 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 7 Feb 2025 15:35:28 -0800 Subject: [PATCH 11/13] Reduce version number of parry2d to avoid experimantal compiler features. --- Cargo.toml | 2 +- src/geometry.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index e05b0710..3deaad9a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -29,7 +29,7 @@ hdf5 = {version = "0.7.1", optional = true} rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} ndarray = {version = "0.14.0", features = ["serde"], optional = true} parry3d-f64 = {optional = true, version="0.2.0"} -parry2d-f64 = {optional = false, version="*"} +parry2d-f64 = {optional = false, version="0.17.0"} [dependencies.pyo3] version = "0.19.0" diff --git a/src/geometry.rs b/src/geometry.rs index dae8e687..814f7e52 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -768,7 +768,7 @@ impl Geometry for ParryMesh2D { let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points)); ParryMesh2D { - trimesh.expect(""), + trimesh, densities, simulation_boundary: simulation_boundary2, boundary: boundary2, From d240c4238db279a38f023b8460092ee2b4200480 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 8 Apr 2025 16:12:16 -0700 Subject: [PATCH 12/13] Fix comment at top of boron_nitride_wire_homogeneous.toml --- examples/boron_nitride_wire_homogeneous.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/boron_nitride_wire_homogeneous.toml b/examples/boron_nitride_wire_homogeneous.toml index 57156b74..9afd461a 100644 --- a/examples/boron_nitride_wire_homogeneous.toml +++ b/examples/boron_nitride_wire_homogeneous.toml @@ -1,4 +1,4 @@ -#Use with 0D geometry option only +#Use with HOMOGENEOUS geometry option only [options] name = "boron_nitride_h_" weak_collision_order = 0 From 21cfff4ef80cb506aa9d343acad532c5e1d446bc Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 8 Apr 2025 16:18:02 -0700 Subject: [PATCH 13/13] Update readme --- README.md | 32 ++++++++++++-------------------- 1 file changed, 12 insertions(+), 20 deletions(-) diff --git a/README.md b/README.md index 6e4223ac..d5b5ee9c 100644 --- a/README.md +++ b/README.md @@ -12,9 +12,9 @@ By discretizing the collision cascade into a sequence of binary collisions, between an energetic ion and a target material. This includes reflection, implantation, and transmission of the incident ion, as well as sputtering and displacement damage of the target. Generally, [BCA] codes can be -valid for incident ion energies between approximately ~1 eV/nucleon +valid for incident ion energies between approximately ~1 eV/nucleon to <1 GeV/nucleon. Improvements to RustBCA have expanded the regime -of validity for some quantities, such as reflection coefficients, below +of validity for some quantities, such as reflection coefficients, below 1 eV/nucleon. Check out the `RustBCA` [Wiki] for detailed information, installation @@ -160,7 +160,7 @@ Windows, MacOS, and Linux systems. * [rcpr], a CPR and polynomial rootfinder, required for using attractive-repulsive interaction potentials such as Lennard-Jones or Morse. * For manipulating input files and running associated scripts, the following are required: * [Python] 3.6+ - * The [Python] libraries: `numpy`, `matplotlib`, `toml` (must build from source), `shapely`, and `scipy`. + * The [Python] libraries: `numpy`, `matplotlib`, `toml`, `shapely`, and `scipy`. ### Detailed instructions for Ubuntu 18.04 LTS @@ -187,36 +187,28 @@ git clone https://github.com/uiri/toml.git cd toml python3 setup.py install ``` -7. (Optional) Install software for [rcpr]: -```bash -sudo apt-get install gcc gfortran build-essential cmake liblapack-dev libblas-dev liblapacke-dev -``` -8. (Optional - should come with rustup) Install `cargo`: +7. (Optional - should come with rustup) Install `cargo`: ```bash sudo apt-get install cargo ``` -9. Build `RustBCA`: +8. Build `RustBCA`: ```bash git clone https://github.com/lcpp-org/rustBCA cd RustBCA cargo build --release ``` -10. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr` (with your choice of backend: `openblas`, `netlib`, or `intel-mkl`): +9. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr`: ```bash -cargo build --release --features cpr_rootfinder_netlib,hdf5_input -cargo build --release --features cpr_rootfinder_openblas,hdf5_input -cargo build --release --features cpr_rootfinder_intel_mkl,hdf5_input +cargo build --release --features cpr_rootfinder,hdf5_input ``` -11. `input.toml` is the input file - see [Usage](https://github.com/lcpp-org/RustBCA/wiki/Usage,-Input-File,-and-Output-Files) for more information -12. Run the required tests using: +10. `input.toml` is the input file - see [Usage](https://github.com/lcpp-org/RustBCA/wiki/Usage,-Input-File,-and-Output-Files) for more information +11. Run the required tests using: ```bash cargo test ``` -13. (Optional) Run the required and optional tests for the desired backend(s): +12. (Optional) Run the required and optional tests for the desired backend(s): ```bash -cargo test --features cpr_rootfinder_netlib -cargo test --features cpr_rootfinder_openblas -cargo test --features cpr_rootfinder_intel_mkl +cargo test --features cpr_rootfinder ``` ### Detailed instructions for Fedora 33 @@ -275,7 +267,7 @@ command line argument. ./RustBCA /path/to/input.toml ``` -Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `TRIMESH`, `SPHERE` - see the wiki for more details): +Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `TRIMESH`, `SPHERE`, `HOMOGENEOUS2D` - see the wiki for more details): ```bash ./RustBCA 0D /path/to/input.toml ```