From 3a409c686316be64682e01f623cb3b70997fbf54 Mon Sep 17 00:00:00 2001 From: jellepoland Date: Tue, 13 Aug 2024 16:05:32 +0200 Subject: [PATCH] Resolved a bug, which was due to newly implemented spring damper force calculation --- Simulations/Gao_et_al.py | 14 +- Simulations/Gao_et_al_sweep_dirty.py | 14 +- Simulations/Simple_config_round.py | 14 +- Tutourial/Poisson's Ratio.ipynb | 2 +- Tutourial/Tutorial 1.ipynb | 20 +- Tutourial/Tutourial2.ipynb | 16 +- .../airbag_problem/airbag_square_cross.py | 6 +- .../hencky_problem/hencky_problem.py | 2 +- .../inflated_membrane.py | 6 +- code_Validation/saddle_form/saddle_form.py | 2 +- .../tether_deflection_gravity.py | 2 +- .../tether_deflection_windFlow.py | 2 +- .../tether_longitudal_oscillations.py | 4 +- .../damping_force/damping_force.py | 4 +- .../gravitational_force.py | 2 +- .../spring_force/spring_force.py | 2 +- src/PSS/ExternalForces/LaserBeam.py | 2 +- .../ExternalForces/OpticalForceCalculator.py | 318 ++++++------ .../optical_interpolators/interpolators.py | 4 +- src/PSS/Mesh/Mesher.py | 6 +- src/PSS/Mesh/mesh_functions.py | 4 +- src/PSS/Mesh/mesh_tester.py | 6 +- src/PSS/Sim/simulations.py | 4 +- src/PSS/particleSystem/Force.py | 2 +- src/PSS/particleSystem/ImplicitForce.py | 4 +- src/PSS/particleSystem/ParticleSystem.py | 4 +- .../TestOpticalForceCalculator.py | 470 ++++++++++-------- test/src/particleSystem/TestParticleSystem.py | 6 +- .../src/particleSystem/test_ParticleSystem.py | 6 +- 29 files changed, 499 insertions(+), 449 deletions(-) diff --git a/Simulations/Gao_et_al.py b/Simulations/Gao_et_al.py index 4da0f0b..97765e5 100644 --- a/Simulations/Gao_et_al.py +++ b/Simulations/Gao_et_al.py @@ -8,15 +8,15 @@ import matplotlib.pyplot as plt from scipy.constants import c -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import Simulate_Lightsail -import Particle_System_Simulator.Mesh.mesh_functions as MF -import Particle_System_Simulator.ExternalForces.optical_interpolators.interpolators as interp -from Particle_System_Simulator.ExternalForces.LaserBeam import LaserBeam -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import Simulate_Lightsail +import PSS.Mesh.mesh_functions as MF +import PSS.ExternalForces.optical_interpolators.interpolators as interp +from PSS.ExternalForces.LaserBeam import LaserBeam +from PSS.ExternalForces.OpticalForceCalculator import ( OpticalForceCalculator, ) -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.ExternalForces.OpticalForceCalculator import ( ParticleOpticalPropertyType, ) diff --git a/Simulations/Gao_et_al_sweep_dirty.py b/Simulations/Gao_et_al_sweep_dirty.py index aabd4ae..44c68e8 100644 --- a/Simulations/Gao_et_al_sweep_dirty.py +++ b/Simulations/Gao_et_al_sweep_dirty.py @@ -13,15 +13,15 @@ import matplotlib.pyplot as plt from scipy.constants import c -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import Simulate_Lightsail -import Particle_System_Simulator.Mesh.mesh_functions as MF -import Particle_System_Simulator.ExternalForces.optical_interpolators.interpolators as interp -from Particle_System_Simulator.ExternalForces.LaserBeam import LaserBeam -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import Simulate_Lightsail +import PSS.Mesh.mesh_functions as MF +import PSS.ExternalForces.optical_interpolators.interpolators as interp +from PSS.ExternalForces.LaserBeam import LaserBeam +from PSS.ExternalForces.OpticalForceCalculator import ( OpticalForceCalculator, ) -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.ExternalForces.OpticalForceCalculator import ( ParticleOpticalPropertyType, ) diff --git a/Simulations/Simple_config_round.py b/Simulations/Simple_config_round.py index 47a3db0..a526ef6 100644 --- a/Simulations/Simple_config_round.py +++ b/Simulations/Simple_config_round.py @@ -19,15 +19,15 @@ from scipy.interpolate import griddata from scipy.interpolate import interp1d -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import Simulate_Lightsail -import Particle_System_Simulator.Mesh.mesh_functions as MF -import Particle_System_Simulator.ExternalForces.optical_interpolators.interpolators as interp -from Particle_System_Simulator.ExternalForces.LaserBeam import LaserBeam -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import Simulate_Lightsail +import PSS.Mesh.mesh_functions as MF +import PSS.ExternalForces.optical_interpolators.interpolators as interp +from PSS.ExternalForces.LaserBeam import LaserBeam +from PSS.ExternalForces.OpticalForceCalculator import ( OpticalForceCalculator, ) -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.ExternalForces.OpticalForceCalculator import ( ParticleOpticalPropertyType, ) diff --git a/Tutourial/Poisson's Ratio.ipynb b/Tutourial/Poisson's Ratio.ipynb index 5bcfa9a..6d92e3d 100644 --- a/Tutourial/Poisson's Ratio.ipynb +++ b/Tutourial/Poisson's Ratio.ipynb @@ -28,7 +28,7 @@ "\n", "sys.path.append(os.path.abspath(\"../.\"))\n", "sys.path.append(os.path.abspath(\"../..\"))\n", - "from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem\n", + "from PSS.particleSystem.ParticleSystem import ParticleSystem\n", "\n", "matplotlib.rcParams[\"figure.figsize\"] = [10, 5]\n", "# %matplotlib qt\n", diff --git a/Tutourial/Tutorial 1.ipynb b/Tutourial/Tutorial 1.ipynb index 74fc55f..7f3d12d 100644 --- a/Tutourial/Tutorial 1.ipynb +++ b/Tutourial/Tutorial 1.ipynb @@ -12,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 1, "id": "db865830", "metadata": {}, "outputs": [], @@ -23,7 +23,7 @@ "import sys, os\n", "\n", "sys.path.append(os.path.abspath(\"..\"))\n", - "from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem\n", + "from PSS.particleSystem.ParticleSystem import ParticleSystem\n", "\n", "# dictionary of required parameters\n", "params = {\n", @@ -77,7 +77,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 2, "id": "ff126088", "metadata": {}, "outputs": [ @@ -161,7 +161,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 3, "id": "d037c6dc", "metadata": {}, "outputs": [ @@ -171,7 +171,7 @@ "Text(0.5, 0.92, 'Initial state')" ] }, - "execution_count": 24, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, @@ -229,7 +229,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 4, "id": "5a3e4150", "metadata": {}, "outputs": [ @@ -287,7 +287,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 5, "id": "f0ccc10f", "metadata": {}, "outputs": [ @@ -297,7 +297,7 @@ "(0.1, 54.7)" ] }, - "execution_count": 13, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, @@ -356,7 +356,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 6, "id": "9f376722", "metadata": {}, "outputs": [ @@ -366,7 +366,7 @@ "Text(0.5, 0.92, 'Final state')" ] }, - "execution_count": 14, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, diff --git a/Tutourial/Tutourial2.ipynb b/Tutourial/Tutourial2.ipynb index 4b61604..2c42b2d 100644 --- a/Tutourial/Tutourial2.ipynb +++ b/Tutourial/Tutourial2.ipynb @@ -51,15 +51,15 @@ "sys.path.append(os.path.abspath(\"..\"))\n", "\n", "# LLS specific imports\n", - "from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem\n", - "from Particle_System_Simulator.Sim.simulations import Simulate_Lightsail\n", - "import Particle_System_Simulator.Mesh.mesh_functions as MF\n", - "import Particle_System_Simulator.ExternalForces.optical_interpolators.interpolators as interp\n", - "from Particle_System_Simulator.ExternalForces.LaserBeam import LaserBeam\n", - "from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import (\n", + "from PSS.particleSystem.ParticleSystem import ParticleSystem\n", + "from PSS.Sim.simulations import Simulate_Lightsail\n", + "import PSS.Mesh.mesh_functions as MF\n", + "import PSS.ExternalForces.optical_interpolators.interpolators as interp\n", + "from PSS.ExternalForces.LaserBeam import LaserBeam\n", + "from PSS.ExternalForces.OpticalForceCalculator import (\n", " OpticalForceCalculator,\n", ")\n", - "from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import (\n", + "from PSS.ExternalForces.OpticalForceCalculator import (\n", " ParticleOpticalPropertyType,\n", ")" ] @@ -681,7 +681,7 @@ "outputs": [], "source": [ "# Importing necessary modules for optical interpolators\n", - "from Particle_System_Simulator.ExternalForces.optical_interpolators import (\n", + "from PSS.ExternalForces.optical_interpolators import (\n", " interpolators as interp,\n", ")\n", "\n", diff --git a/code_Validation/airbag_problem/airbag_square_cross.py b/code_Validation/airbag_problem/airbag_square_cross.py index a81d0cb..4804e3d 100644 --- a/code_Validation/airbag_problem/airbag_square_cross.py +++ b/code_Validation/airbag_problem/airbag_square_cross.py @@ -15,9 +15,9 @@ import numpy as np import matplotlib.pyplot as plt -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import Simulate_airbag -import Particle_System_Simulator.Mesh.mesh_functions as MF +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import Simulate_airbag +import PSS.Mesh.mesh_functions as MF params = { # model parameters diff --git a/code_Validation/hencky_problem/hencky_problem.py b/code_Validation/hencky_problem/hencky_problem.py index 7a9e8a3..315567a 100644 --- a/code_Validation/hencky_problem/hencky_problem.py +++ b/code_Validation/hencky_problem/hencky_problem.py @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt import pandas as pd import time -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem +from PSS.particleSystem.ParticleSystem import ParticleSystem def instantiate_ps(): diff --git a/code_Validation/inflated_SiN_membrane/inflated_membrane.py b/code_Validation/inflated_SiN_membrane/inflated_membrane.py index 732bcf8..8ab524a 100644 --- a/code_Validation/inflated_SiN_membrane/inflated_membrane.py +++ b/code_Validation/inflated_SiN_membrane/inflated_membrane.py @@ -55,9 +55,9 @@ import numpy as np import matplotlib.pyplot as plt -from Particle_System_Simulator.particleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import Simulate_airbag -import Particle_System_Simulator.Mesh.mesh_functions as MF +from PSS.particleSystem import ParticleSystem +from PSS.Sim.simulations import Simulate_airbag +import PSS.Mesh.mesh_functions as MF curve_fitting_parameters = { "125": {"A1/s": 1.094e-02, "A2/s": 1.630e-04, "A1/r": 6.462e-03, "A2/r": 9.490e-05}, diff --git a/code_Validation/saddle_form/saddle_form.py b/code_Validation/saddle_form/saddle_form.py index d3e73bb..7fb579e 100644 --- a/code_Validation/saddle_form/saddle_form.py +++ b/code_Validation/saddle_form/saddle_form.py @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt import pandas as pd import time -from Particle_System_Simulator.particleSystem import ParticleSystem +from PSS.particleSystem import ParticleSystem def instantiate_ps(): diff --git a/code_Validation/tether_deflection_gravity/tether_deflection_gravity.py b/code_Validation/tether_deflection_gravity/tether_deflection_gravity.py index ef52a34..6cf290b 100644 --- a/code_Validation/tether_deflection_gravity/tether_deflection_gravity.py +++ b/code_Validation/tether_deflection_gravity/tether_deflection_gravity.py @@ -10,7 +10,7 @@ import pandas as pd import sys import time -from Particle_System_Simulator.particleSystem import ParticleSystem +from PSS.particleSystem import ParticleSystem from sympy import * diff --git a/code_Validation/tether_deflection_windFlow/tether_deflection_windFlow.py b/code_Validation/tether_deflection_windFlow/tether_deflection_windFlow.py index bbba31e..2e1cad7 100644 --- a/code_Validation/tether_deflection_windFlow/tether_deflection_windFlow.py +++ b/code_Validation/tether_deflection_windFlow/tether_deflection_windFlow.py @@ -10,7 +10,7 @@ import pandas as pd import sys import time -from Particle_System_Simulator.particleSystem import ParticleSystem +from PSS.particleSystem import ParticleSystem def instantiate_ps(): diff --git a/code_Validation/tether_longitudal_oscillations/tether_longitudal_oscillations.py b/code_Validation/tether_longitudal_oscillations/tether_longitudal_oscillations.py index 7b75e8c..58c4ff8 100644 --- a/code_Validation/tether_longitudal_oscillations/tether_longitudal_oscillations.py +++ b/code_Validation/tether_longitudal_oscillations/tether_longitudal_oscillations.py @@ -10,8 +10,8 @@ import pandas as pd import sys import time -from Particle_System_Simulator.particleSystem import ParticleSystem -from Particle_System_Simulator.AnalysisModules import system_energy +from PSS.particleSystem import ParticleSystem +from PSS.AnalysisModules import system_energy def instantiate_ps(): diff --git a/code_Verification/damping_force/damping_force.py b/code_Verification/damping_force/damping_force.py index ca05d0e..df21fc7 100644 --- a/code_Verification/damping_force/damping_force.py +++ b/code_Verification/damping_force/damping_force.py @@ -10,8 +10,8 @@ import pandas as pd from scipy.integrate import solve_ivp import sys -from Particle_System_Simulator.particleSystem import ParticleSystem -from Particle_System_Simulator.AnalysisModules import system_energy +from PSS.particleSystem import ParticleSystem +from PSS.AnalysisModules import system_energy def instantiate_ps(): diff --git a/code_Verification/gravitational_force/gravitational_force.py b/code_Verification/gravitational_force/gravitational_force.py index 7174342..fc2657c 100644 --- a/code_Verification/gravitational_force/gravitational_force.py +++ b/code_Verification/gravitational_force/gravitational_force.py @@ -8,7 +8,7 @@ import matplotlib.pyplot as plt import pandas as pd import sys -from Particle_System_Simulator.particleSystem import ParticleSystem +from PSS.particleSystem import ParticleSystem def instantiate_ps(): diff --git a/code_Verification/spring_force/spring_force.py b/code_Verification/spring_force/spring_force.py index 59de011..51455d3 100644 --- a/code_Verification/spring_force/spring_force.py +++ b/code_Verification/spring_force/spring_force.py @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt import pandas as pd import sys -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem +from PSS.particleSystem.ParticleSystem import ParticleSystem def connectivity_matrix(): diff --git a/src/PSS/ExternalForces/LaserBeam.py b/src/PSS/ExternalForces/LaserBeam.py index b623d82..fc878ee 100644 --- a/src/PSS/ExternalForces/LaserBeam.py +++ b/src/PSS/ExternalForces/LaserBeam.py @@ -6,7 +6,7 @@ import numpy as np import matplotlib.pyplot as plt -from Particle_System_Simulator.particleSystem.SystemObject import SystemObject +from PSS.particleSystem.SystemObject import SystemObject class LaserBeam(SystemObject): diff --git a/src/PSS/ExternalForces/OpticalForceCalculator.py b/src/PSS/ExternalForces/OpticalForceCalculator.py index 40131cc..090d008 100644 --- a/src/PSS/ExternalForces/OpticalForceCalculator.py +++ b/src/PSS/ExternalForces/OpticalForceCalculator.py @@ -15,25 +15,28 @@ import scipy as sp from scipy.constants import c from scipy.spatial.transform import Rotation -from Particle_System_Simulator.particleSystem.Force import Force +from PSS.particleSystem.Force import Force import logging + class OpticalForceCalculator(Force): """ Handles the calculation of forces arising from optical pressure """ + def __init__(self, ParticleSystem, LaserBeam): self.ParticleSystem = ParticleSystem - self.PS = self.ParticleSystem #alias for convenience + self.PS = self.ParticleSystem # alias for convenience self.LaserBeam = LaserBeam - if not hasattr(self.ParticleSystem.particles[0],'optical_type'): - raise AttributeError("ParticleSystem does not have any optical properties set!") + if not hasattr(self.ParticleSystem.particles[0], "optical_type"): + raise AttributeError( + "ParticleSystem does not have any optical properties set!" + ) super().__init__() return - def __str__(self): print("OpticalForceCalculator object instantiated with attributes:") print(f"ParticleSystem: \n {self.ParticleSystem}") @@ -57,38 +60,42 @@ def force_value(self): locations, _ = PS.x_v_current_3D forces = np.zeros(locations.shape) - if not hasattr(self, 'optical_type_mask'): + if not hasattr(self, "optical_type_mask"): self.create_optical_type_mask() # ! Note ! This bakes in implicitly that the orientation of the light # vector is in z+ direction - intensity_vectors = np.array([[0,0,LB.intensity_profile(x,y)] for x,y,z in locations]) - polarisation_vectors = LB.polarization_map(locations[:,0],locations[:,1]) + intensity_vectors = np.array( + [[0, 0, LB.intensity_profile(x, y)] for x, y, z in locations] + ) + polarisation_vectors = LB.polarization_map(locations[:, 0], locations[:, 1]) for optical_type in self.optical_type_mask.keys(): if optical_type == ParticleOpticalPropertyType.SPECULAR: mask = self.optical_type_mask[optical_type] - forces[mask] = self.calculate_specular_force(area_vectors[mask], - intensity_vectors[mask]) + forces[mask] = self.calculate_specular_force( + area_vectors[mask], intensity_vectors[mask] + ) elif optical_type == ParticleOpticalPropertyType.AXICONGRATING: mask = self.optical_type_mask[optical_type] filtered_particles = compress(PS.particles, mask) axicon_angle = [p.axicon_angle for p in filtered_particles] - forces[mask] = self.calculate_axicongrating_force(area_vectors[mask], - intensity_vectors[mask], - axicon_angle) + forces[mask] = self.calculate_axicongrating_force( + area_vectors[mask], intensity_vectors[mask], axicon_angle + ) elif optical_type == ParticleOpticalPropertyType.ARBITRARY_PHC: mask = self.optical_type_mask[optical_type] - if not hasattr(self, 'optical_interpolators'): + if not hasattr(self, "optical_interpolators"): self.create_phc_map(mask) - - forces[mask] = self.calculate_arbitrary_phc_force(area_vectors[mask], - intensity_vectors[mask], - polarisation_vectors[mask], - self.optical_interpolators) + forces[mask] = self.calculate_arbitrary_phc_force( + area_vectors[mask], + intensity_vectors[mask], + polarisation_vectors[mask], + self.optical_interpolators, + ) return forces def create_phc_map(self, mask): @@ -101,7 +108,7 @@ def create_phc_map(self, mask): self.phc_dict = defaultdict(list) for i, p in enumerate(filtered_particles): self.optical_interpolators.append(p.optical_interpolator) - #build dict of mask entries and which PhC they belong to. + # build dict of mask entries and which PhC they belong to. self.phc_dict[p.optical_interpolator].append(original_indices[i]) for key in self.phc_dict.keys(): @@ -114,7 +121,6 @@ def create_submask(self, indice_list): submask[i] = True return submask - def calculate_specular_force(self, area_vectors, intensity_vectors): """ Calculates forces for particles of optical type 'specular' @@ -134,8 +140,8 @@ def calculate_specular_force(self, area_vectors, intensity_vectors): flattened array of external forces of length 3 * n_particles. """ # First we compute the incident power on the particle areas - abs_area_vectors = area_vectors[:,2] # assumes z+ poynting vector - abs_intensity_vectors = intensity_vectors[:,2] # assumes z+ poynting vector + abs_area_vectors = area_vectors[:, 2] # assumes z+ poynting vector + abs_intensity_vectors = intensity_vectors[:, 2] # assumes z+ poynting vector incident_power = abs_area_vectors * abs_intensity_vectors # To get the direction of the forces we need to normalise the area @@ -144,14 +150,13 @@ def calculate_specular_force(self, area_vectors, intensity_vectors): norms = np.linalg.norm(area_vectors, axis=1) forces = area_vectors.copy() for i in range(3): - forces[:,i] *= 2*incident_power / (c*norms) + forces[:, i] *= 2 * incident_power / (c * norms) return forces - def calculate_axicongrating_force(self, - area_vectors, - intensity_vectors, - axicon_angle): + def calculate_axicongrating_force( + self, area_vectors, intensity_vectors, axicon_angle + ): """ Calculates forces for particles of optical type 'axicon grating' @@ -175,26 +180,28 @@ def calculate_axicongrating_force(self, forces = self.calculate_specular_force(area_vectors, intensity_vectors) forces = rotation_super_matrix.dot(np.hstack(forces).T) - forces = np.reshape(forces, [int(forces.shape[0]/3),3]) + forces = np.reshape(forces, [int(forces.shape[0] / 3), 3]) # The forces need to be scaled to account for the fact that # |[1,1]| != |[1]|+|[1]| # We don't have acces to the angle, but we can make use of the cosine # rule: cos(alpha) = A.dot(B) / (|A| |B|) to get the angle between # z+ and the line of action of the force. - unit_z = np.array([0,0,1]) + unit_z = np.array([0, 0, 1]) scaling_factor = np.matmul(axicon_angle, unit_z).dot(unit_z) for i in range(3): - forces[:,i]*=scaling_factor + forces[:, i] *= scaling_factor return forces - def calculate_arbitrary_phc_force(self, - area_vectors, - intensity_vectors, - polarisation_vectors, - optical_interpolators): + def calculate_arbitrary_phc_force( + self, + area_vectors, + intensity_vectors, + polarisation_vectors, + optical_interpolators, + ): """ Calculates forces for particles of optical type 'axicon grating' @@ -223,16 +230,13 @@ def calculate_arbitrary_phc_force(self, # azimuth_angles += np.pi # azimuth_angles %= 2*np.pi - # Condition them for the interpolator: - wrapped_coordinates = wrap_spherical_coordinates(polar_angles, - azimuth_angles, - polarisation_angles) + wrapped_coordinates = wrap_spherical_coordinates( + polar_angles, azimuth_angles, polarisation_angles + ) polar_angles, azimuth_angles, polarisation_angles = wrapped_coordinates - incoming_ray = np.vstack((polar_angles, - azimuth_angles, - polarisation_angles)).T + incoming_ray = np.vstack((polar_angles, azimuth_angles, polarisation_angles)).T # Find directions of outgoing rays # Interpolator([polar_in, azimuth_in, polarization_in])->[polar_out, azimuth_out, magnitude] @@ -242,7 +246,6 @@ def calculate_arbitrary_phc_force(self, submask = self.phc_dict[phc] reflected_ray[submask] = phc(incoming_ray[submask]) - # reflected_ray = [interp(incoming_ray[i]) # reflected_ray = [interp(tuple(incoming_ray[i])) # for i, interp @@ -252,19 +255,24 @@ def calculate_arbitrary_phc_force(self, # Switch reference frame again; made easier becasue we are going to [0,0,1] polar_angles_out += polar_angles - reflected_vectors = spherical_to_cartesian(polar_angles_out, - azimuth_angles_out, - magnitudes) + reflected_vectors = spherical_to_cartesian( + polar_angles_out, azimuth_angles_out, magnitudes + ) # Compute the incident power on the particle areas - abs_area_vectors = area_vectors[:,2] # assumes z+ poynting vector - abs_intensity_vectors = intensity_vectors[:,2] # assumes z+ poynting vector + abs_area_vectors = area_vectors[:, 2] # assumes z+ poynting vector + abs_intensity_vectors = intensity_vectors[:, 2] # assumes z+ poynting vector incident_power = abs_area_vectors * abs_intensity_vectors - scattered_power = reflected_vectors*incident_power[:,np.newaxis] + scattered_power = reflected_vectors * incident_power[:, np.newaxis] - net_power = np.hstack((np.zeros((incident_power.shape[0],2)),incident_power[:,np.newaxis])) + scattered_power - forces = net_power/c + net_power = ( + np.hstack( + (np.zeros((incident_power.shape[0], 2)), incident_power[:, np.newaxis]) + ) + + scattered_power + ) + forces = net_power / c return forces @@ -284,26 +292,27 @@ def create_optical_type_mask(self): optical_type_list = [] error_index_list = [] for i, particle in enumerate(self.ParticleSystem.particles): - if hasattr(particle, 'optical_type'): + if hasattr(particle, "optical_type"): optical_type_list.append(particle.optical_type) else: error_index_list.append(i) - if len(error_index_list)>0: - raise AttributeError("All particles should have an optical type" - " set prior to calculation of optical forces." - " Currently the particles with indices" - f" {error_index_list} have no property set") + if len(error_index_list) > 0: + raise AttributeError( + "All particles should have an optical type" + " set prior to calculation of optical forces." + " Currently the particles with indices" + f" {error_index_list} have no property set" + ) optical_type_list = np.array(optical_type_list) self.optical_type_mask = {} for optical_type in ParticleOpticalPropertyType: mask = optical_type_list == optical_type - if sum(mask)>0: + if sum(mask) > 0: self.optical_type_mask[optical_type] = mask - - def calculate_stability_coefficients(self, displacement_range = [0.1, 5]): + def calculate_stability_coefficients(self, displacement_range=[0.1, 5]): """ Calculates the stability coefficients for the particle system @@ -324,21 +333,24 @@ def calculate_stability_coefficients(self, displacement_range = [0.1, 5]): """ q, alpha = displacement_range - displacement_vectors = np.array([[q,0,0,0,0,0], - [0,q,0,0,0,0], - [0,0,q,0,0,0], - [0,0,0,alpha,0,0], - [0,0,0,0,alpha,0], - [0,0,0,0,0,alpha]]) - - jacobian = np.zeros((6,6)) + displacement_vectors = np.array( + [ + [q, 0, 0, 0, 0, 0], + [0, q, 0, 0, 0, 0], + [0, 0, q, 0, 0, 0], + [0, 0, 0, alpha, 0, 0], + [0, 0, 0, 0, alpha, 0], + [0, 0, 0, 0, 0, alpha], + ] + ) + + jacobian = np.zeros((6, 6)) for i, vector in enumerate(displacement_vectors): - jacobian[:,i] =np.hstack(self.calculate_force_gradient(vector)) + jacobian[:, i] = np.hstack(self.calculate_force_gradient(vector)) return jacobian - - def calculate_force_gradient(self, displacement_vector : npt.ArrayLike): + def calculate_force_gradient(self, displacement_vector: npt.ArrayLike): """ Calculates force and moment coefficients of ParticleSystems based on a 1 DOF displacement @@ -359,22 +371,23 @@ def calculate_force_gradient(self, displacement_vector : npt.ArrayLike): k_rot : TYPE lenght 3 list of translational reaction coefficients [dM_x/dx__i, dM_y/dx__i, dM_z/dx__i] """ - displacement = displacement_vector[displacement_vector !=0] - if len(displacement)>1: - raise AttributeError("Expected vector with only one nonzero value," - f"instead got {displacement_vector}") + displacement = displacement_vector[displacement_vector != 0] + if len(displacement) > 1: + raise AttributeError( + "Expected vector with only one nonzero value," + f"instead got {displacement_vector}" + ) original = self.calculate_restoring_forces() self.PS.displace(displacement_vector) reaction = self.calculate_restoring_forces() self.PS.un_displace() - k_trans = (reaction[0] - original[0])/displacement - k_rot = (reaction[1] - original[1])/displacement + k_trans = (reaction[0] - original[0]) / displacement + k_rot = (reaction[1] - original[1]) / displacement return k_trans, k_rot - - def calculate_restoring_forces(self, forces : npt.ArrayLike= None): + def calculate_restoring_forces(self, forces: npt.ArrayLike = None): """ calculates net forces and moments around the center of mass @@ -395,16 +408,19 @@ def calculate_restoring_forces(self, forces : npt.ArrayLike= None): PS = self.ParticleSystem if type(forces) == type(None): forces = self.force_value() - net_force = np.sum(forces,axis=0) + net_force = np.sum(forces, axis=0) - COM =PS.calculate_center_of_mass() + COM = PS.calculate_center_of_mass() locations, _ = PS.x_v_current_3D - moment_arms = PS.translate_mesh(locations, -COM) # note: this doesn't displace the PS, just applies a transformation on the 'locations' variable + moment_arms = PS.translate_mesh( + locations, -COM + ) # note: this doesn't displace the PS, just applies a transformation on the 'locations' variable moments = np.cross(moment_arms, forces) - net_moments = np.sum(moments,axis=0) + net_moments = np.sum(moments, axis=0) return net_force, net_moments + class ParticleOpticalPropertyType(Enum): """ Enumeration representing the various types of optical properties for the Particles @@ -429,11 +445,13 @@ class ParticleOpticalPropertyType(Enum): AXICONGRATING = "axicongrating" ARBITRARY_PHC = "ARBITRARY_PHC" -vectorized_optical_type_retriever = np.vectorize(lambda p: p.optical_type) -def compute_spherical_coordinates(area_vectors: npt.NDArray, - polarisation_vectors: npt.NDArray) -> ( - npt.NDArray, npt.NDArray, npt.NDArray): +vectorized_optical_type_retriever = np.vectorize(lambda p: p.optical_type) + + +def compute_spherical_coordinates( + area_vectors: npt.NDArray, polarisation_vectors: npt.NDArray +) -> (npt.NDArray, npt.NDArray, npt.NDArray): """ Computes the polar angles, azimuth angles, and polarisation angles of the incoming ray and its polarisation relative to the orientation of area elements represented by area vectors. @@ -457,24 +475,25 @@ def compute_spherical_coordinates(area_vectors: npt.NDArray, orthogonal to the area vectors [rad]. """ # Normalize the area vectors - norm_area_vectors = area_vectors / np.linalg.norm(area_vectors, axis=1)[:, np.newaxis] + norm_area_vectors = ( + area_vectors / np.linalg.norm(area_vectors, axis=1)[:, np.newaxis] + ) # Compute polar angles using the dot product between area vectors and the z-axis polar_angles = np.arccos(norm_area_vectors[:, 2]) # Compute azimuth angles - azimuth_angles = np.arctan2(norm_area_vectors[:,1], norm_area_vectors[:,0]) + azimuth_angles = np.arctan2(norm_area_vectors[:, 1], norm_area_vectors[:, 0]) # Compute the polarisation angle in cartesian space: - polarisation_angles = np.arccos(polarisation_vectors[:,0]) - + polarisation_angles = np.arccos(polarisation_vectors[:, 0]) return polar_angles, azimuth_angles, polarisation_angles -def spherical_to_cartesian(polar_angles: npt.NDArray, - azimuth_angles: npt.NDArray, - magnitudes: npt.NDArray) -> npt.NDArray: +def spherical_to_cartesian( + polar_angles: npt.NDArray, azimuth_angles: npt.NDArray, magnitudes: npt.NDArray +) -> npt.NDArray: """ Converts spherical coordinates back to Cartesian coordinates in the global frame, using the area vectors to define the local reference frames. Scales the resulting vectors by the magnitudes. @@ -501,11 +520,12 @@ def spherical_to_cartesian(polar_angles: npt.NDArray, y = magnitudes * np.sin(polar_angles) * np.sin(azimuth_angles) z = magnitudes * np.cos(polar_angles) - cartesian_vectors= np.vstack((x,y,z)).T + cartesian_vectors = np.vstack((x, y, z)).T return cartesian_vectors -def cartesian_to_sphereical(vectors:npt.NDArray) -> npt.NDArray: + +def cartesian_to_sphereical(vectors: npt.NDArray) -> npt.NDArray: """ Converts a set of vectors from cartesian to spherical coordinates @@ -523,16 +543,16 @@ def cartesian_to_sphereical(vectors:npt.NDArray) -> npt.NDArray: magnitudes : npt.NDArray An array of magnitudes of the vectors. """ - x,y,z = vectors[:,0], vectors[:,1], vectors[:,2] + x, y, z = vectors[:, 0], vectors[:, 1], vectors[:, 2] magnitudes = np.linalg.norm(vectors, axis=1) - polar_angles = np.arccos(z/magnitudes) - azimuth_angles = np.arctan2(y,x) + polar_angles = np.arccos(z / magnitudes) + azimuth_angles = np.arctan2(y, x) return polar_angles, azimuth_angles, magnitudes -def wrap_spherical_coordinates(theta: npt.NDArray, - phi: npt.NDArray, - pol: npt.NDArray=None): +def wrap_spherical_coordinates( + theta: npt.NDArray, phi: npt.NDArray, pol: npt.NDArray = None +): """ wraps points in spherical coordinates to always stay within the interpolators defined range @@ -555,19 +575,19 @@ def wrap_spherical_coordinates(theta: npt.NDArray, pol : npt.NDArray polarisation angle """ - phi[theta>np.pi] += np.pi - theta[theta>np.pi]= np.pi - theta[theta>np.pi]%np.pi + phi[theta > np.pi] += np.pi + theta[theta > np.pi] = np.pi - theta[theta > np.pi] % np.pi - phi[theta<0] += np.pi - theta[theta<0] *=-1 - phi %= 2*np.pi + phi[theta < 0] += np.pi + theta[theta < 0] *= -1 + phi %= 2 * np.pi - phi[abs(phi-2*np.pi)<1e-5]=0 # wraps values that are _almost_ 2*np.pi + phi[abs(phi - 2 * np.pi) < 1e-5] = 0 # wraps values that are _almost_ 2*np.pi if np.any(pol != None): x = abs(np.cos(pol)) y = abs(np.sin(pol)) - pol = np.arctan(y/x) + pol = np.arctan(y / x) return theta, phi, pol else: @@ -575,36 +595,42 @@ def wrap_spherical_coordinates(theta: npt.NDArray, # quick test # !!! todo move to testing file - theta = np.random.random(100)*3*np.pi - np.pi - phi = np.random.random(100)*3*np.pi - np.pi - pol = np.random.random(100)*3*np.pi - np.pi + theta = np.random.random(100) * 3 * np.pi - np.pi + phi = np.random.random(100) * 3 * np.pi - np.pi + pol = np.random.random(100) * 3 * np.pi - np.pi + def test_wrap_spherical_coordinates(dat): mags = np.ones(dat[0].shape) t1 = spherical_to_cartesian(*wrap_spherical_coordinates(*dat)[:2], mags) t2 = spherical_to_cartesian(*dat[:2], mags) - return (t1==t2).all() - test_wrap_spherical_coordinates((theta,phi,pol)) + return (t1 == t2).all() + + test_wrap_spherical_coordinates((theta, phi, pol)) + if __name__ == "__main__": from code_Validation.saddle_form import saddle_form - from Particle_System_Simulator.ExternalForces.LaserBeam import LaserBeam + from PSS.ExternalForces.LaserBeam import LaserBeam import matplotlib.pyplot as plt - PS = saddle_form.instantiate_ps() - #PS.stress_self() - #for i in range(10): PS.simulate() + # PS.stress_self() + # for i in range(10): PS.simulate() for particle in PS.particles: - particle.x[2]= 0 + particle.x[2] = 0 - I_0 = 100e9 /(10*10) + I_0 = 100e9 / (10 * 10) mu_x = 5 mu_y = 5 sigma = 5 - LB = LaserBeam(lambda x, y: I_0 * np.exp(-1/2 *((x-mu_x)/sigma)**2 - -1/2 *((y-mu_y)/sigma)**2), - lambda x,y: [0,1]) - LB = LaserBeam(lambda x, y: np.ones(x.shape)*I_0, lambda x,y: [0,1]) + LB = LaserBeam( + lambda x, y: I_0 + * np.exp( + -1 / 2 * ((x - mu_x) / sigma) ** 2 - 1 / 2 * ((y - mu_y) / sigma) ** 2 + ), + lambda x, y: [0, 1], + ) + LB = LaserBeam(lambda x, y: np.ones(x.shape) * I_0, lambda x, y: [0, 1]) # One half of example will be 45 deg axicon angle directed towards (5,5) # other half will be specular reflection @@ -613,40 +639,32 @@ def test_wrap_spherical_coordinates(dat): for particle in PS.particles: particle.optical_type = ParticleOpticalPropertyType.SPECULAR - if (particle.x[0]-5)**2 + (particle.x[1]-5)**2>= 3**2: + if (particle.x[0] - 5) ** 2 + (particle.x[1] - 5) ** 2 >= 3**2: roty = 45 - rotz = np.rad2deg(np.arctan2((particle.x[1]-5), (particle.x[0]-4.999))) + rotz = np.rad2deg(np.arctan2((particle.x[1] - 5), (particle.x[0] - 4.999))) particle.optical_type = ParticleOpticalPropertyType.AXICONGRATING - #particle.axicon_angle = Rotation.from_euler('yz', [roty, rotz], degrees=True).as_matrix() - particle.axicon_angle = Rotation.from_euler('yz', [roty, rotz], degrees=True).as_matrix() - rots.append((roty,rotz%360)) + # particle.axicon_angle = Rotation.from_euler('yz', [roty, rotz], degrees=True).as_matrix() + particle.axicon_angle = Rotation.from_euler( + "yz", [roty, rotz], degrees=True + ).as_matrix() + rots.append((roty, rotz % 360)) OFC = OpticalForceCalculator(PS, LB) forces = OFC.force_value() - ax = PS.plot() points, _ = PS.x_v_current_3D - x,y,z = points[:,0], points[:,1], points[:,2] - a_u = forces[:,0] - a_v = forces[:,1] - a_w = forces[:,2] - ax.scatter(x,y,z) - ax.quiver(x,y,z,a_u,a_v,a_w, length = 0.1) - ax.set_box_aspect([1,1,1]) + x, y, z = points[:, 0], points[:, 1], points[:, 2] + a_u = forces[:, 0] + a_v = forces[:, 1] + a_w = forces[:, 2] + ax.scatter(x, y, z) + ax.quiver(x, y, z, a_u, a_v, a_w, length=0.1) + ax.set_box_aspect([1, 1, 1]) ax.set_zlim(-5, 5) - #ax.set_zscale('symlog') - - #ax2 = fig.add_subplot(projection='3d') - #LB.plot(ax2, x_range = [0,10], y_range=[0,10]) - - - - - - - - + # ax.set_zscale('symlog') + # ax2 = fig.add_subplot(projection='3d') + # LB.plot(ax2, x_range = [0,10], y_range=[0,10]) diff --git a/src/PSS/ExternalForces/optical_interpolators/interpolators.py b/src/PSS/ExternalForces/optical_interpolators/interpolators.py index 891136c..f33669d 100644 --- a/src/PSS/ExternalForces/optical_interpolators/interpolators.py +++ b/src/PSS/ExternalForces/optical_interpolators/interpolators.py @@ -16,8 +16,8 @@ from scipy.spatial import KDTree import matplotlib.pyplot as plt -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import ( +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.ExternalForces.OpticalForceCalculator import ( wrap_spherical_coordinates, ) diff --git a/src/PSS/Mesh/Mesher.py b/src/PSS/Mesh/Mesher.py index 565db19..e31158f 100644 --- a/src/PSS/Mesh/Mesher.py +++ b/src/PSS/Mesh/Mesher.py @@ -91,9 +91,9 @@ def mesh_square_cross(self): if __name__ == "__main__": - from Particle_System_Simulator.particleSystem.Particle import Particle - from Particle_System_Simulator.particleSystem.SpringDamper import SpringDamper - from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem + from PSS.particleSystem.Particle import Particle + from PSS.particleSystem.SpringDamper import SpringDamper + from PSS.particleSystem.ParticleSystem import ParticleSystem mesh = Mesher() mesh.add_rectangle(0, 0, "square") diff --git a/src/PSS/Mesh/mesh_functions.py b/src/PSS/Mesh/mesh_functions.py index 7c4b727..7d99424 100644 --- a/src/PSS/Mesh/mesh_functions.py +++ b/src/PSS/Mesh/mesh_functions.py @@ -7,7 +7,7 @@ import numpy as np from scipy.spatial.transform import Rotation as R -from Particle_System_Simulator.particleSystem.SpringDamper import SpringDamperType +from PSS.particleSystem.SpringDamper import SpringDamperType params = { @@ -618,7 +618,7 @@ def ps_find_strip_dimentions(ParticleSystem, midstrip): if __name__ == "__main__": - from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem + from PSS.particleSystem.ParticleSystem import ParticleSystem import matplotlib.pyplot as plt params = { diff --git a/src/PSS/Mesh/mesh_tester.py b/src/PSS/Mesh/mesh_tester.py index e7838db..f6d717c 100644 --- a/src/PSS/Mesh/mesh_tester.py +++ b/src/PSS/Mesh/mesh_tester.py @@ -11,9 +11,9 @@ import numpy as np import matplotlib.pyplot as plt -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import Simulate_1d_Stretch -import Particle_System_Simulator.Mesh.mesh_functions as MF +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import Simulate_1d_Stretch +import PSS.Mesh.mesh_functions as MF class AnalyseMaterial: diff --git a/src/PSS/Sim/simulations.py b/src/PSS/Sim/simulations.py index 14e11a2..b8b5ec8 100644 --- a/src/PSS/Sim/simulations.py +++ b/src/PSS/Sim/simulations.py @@ -15,8 +15,8 @@ import numpy as np import matplotlib.pyplot as plt -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Mesh import mesh_functions as MF +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Mesh import mesh_functions as MF class Simulate: diff --git a/src/PSS/particleSystem/Force.py b/src/PSS/particleSystem/Force.py index 894bcc9..06676dd 100644 --- a/src/PSS/particleSystem/Force.py +++ b/src/PSS/particleSystem/Force.py @@ -2,7 +2,7 @@ Child Abstract Base Class 'Force', for force objects to be instantiated in ParticleSystem """ -from Particle_System_Simulator.particleSystem.SystemObject import SystemObject +from PSS.particleSystem.SystemObject import SystemObject from abc import abstractmethod diff --git a/src/PSS/particleSystem/ImplicitForce.py b/src/PSS/particleSystem/ImplicitForce.py index 615ada1..e825d23 100644 --- a/src/PSS/particleSystem/ImplicitForce.py +++ b/src/PSS/particleSystem/ImplicitForce.py @@ -2,8 +2,8 @@ Child Abstract Base Class 'ImplicitForce', for implicit force objects to be instantiated in ParticleSystem """ -from Particle_System_Simulator.particleSystem.Force import Force -from Particle_System_Simulator.particleSystem.Particle import Particle +from PSS.particleSystem.Force import Force +from PSS.particleSystem.Particle import Particle from abc import abstractmethod from abc import abstractproperty diff --git a/src/PSS/particleSystem/ParticleSystem.py b/src/PSS/particleSystem/ParticleSystem.py index c23ebee..15d9534 100644 --- a/src/PSS/particleSystem/ParticleSystem.py +++ b/src/PSS/particleSystem/ParticleSystem.py @@ -419,14 +419,14 @@ def __one_d_force_vector(self): delta_length_pulley_other_line = norm_p3p4 - rest_length_p3p4 f_int = SD.force_value(delta_length_pulley_other_line) - i, j = self.__connectivity_matrix[idx] + i, j, *_ = self.__connectivity_matrix[idx] self.__f[i * 3 : i * 3 + 3] += f_int self.__f[j * 3 : j * 3 + 3] -= f_int # if regular spring-damper else: f_int = self.__springdampers[idx].force_value() - i, j = self.__connectivity_matrix[idx] + i, j, *_ = self.__connectivity_matrix[idx] self.__f[i * 3 : i * 3 + 3] += f_int self.__f[j * 3 : j * 3 + 3] -= f_int diff --git a/test/src/ExternalForces/TestOpticalForceCalculator.py b/test/src/ExternalForces/TestOpticalForceCalculator.py index dddcd59..469ea79 100644 --- a/test/src/ExternalForces/TestOpticalForceCalculator.py +++ b/test/src/ExternalForces/TestOpticalForceCalculator.py @@ -9,176 +9,192 @@ import numpy as np -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.ExternalForces.OpticalForceCalculator import OpticalForceCalculator, ParticleOpticalPropertyType -from Particle_System_Simulator.ExternalForces.LaserBeam import LaserBeam +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.ExternalForces.OpticalForceCalculator import ( + OpticalForceCalculator, + ParticleOpticalPropertyType, +) +from PSS.ExternalForces.LaserBeam import LaserBeam from scipy.constants import c from scipy.spatial.transform import Rotation -import Particle_System_Simulator.Mesh.mesh_functions as MF +import PSS.Mesh.mesh_functions as MF + class TestOpticalForceCalculator(unittest.TestCase): def setUp(self): - + # Set up ParticleSystem self.params = { # model parameters "k": 1, # [N/m] spring stiffness "k_d": 1, # [N/m] spring stiffness for diagonal elements "c": 10, # [N s/m] damping coefficient - "m_segment": 1, # [kg] mass of each node - + "m_segment": 1, # [kg] mass of each node # simulation settings "dt": 0.1, # [s] simulation timestep "t_steps": 1000, # [-] number of simulated time steps "abs_tol": 1e-50, # [m/s] absolute error tolerance iterative solver "rel_tol": 1e-5, # [-] relative error tolerance iterative solver "max_iter": 1e5, # [-] maximum number of iterations] - "convergence_threshold": 1e-6, # [-] - "min_iterations": 10, # [-] - } - self.initial_conditions, self.connectivity_matrix = MF.mesh_square(1, 1, 0.02, self.params) - self.PS = ParticleSystem(self.connectivity_matrix, - self.initial_conditions, - self.params, - clean_particles = False) - - self.PS.initialize_find_surface() - self.PS.calculate_correct_masses(1e-3, # [m] thickness - 1000) # [kg/m3] density - self.mass = 1*1*1e-3 * 1000 - + "convergence_threshold": 1e-6, # [-] + "min_iterations": 10, # [-] + } + self.initial_conditions, self.connectivity_matrix = MF.mesh_square( + 1, 1, 0.02, self.params + ) + self.PS = ParticleSystem( + self.connectivity_matrix, + self.initial_conditions, + self.params, + clean_particles=False, + ) + + self.PS.initialize_find_surface() + self.PS.calculate_correct_masses(1e-3, 1000) # [m] thickness # [kg/m3] density + self.mass = 1 * 1 * 1e-3 * 1000 + # Set up some sample LaserBeam instances - I_0 = 100e9 /(10*10) # 100 GW divided over 100 square meters - self.I_0 = I_0 #needed for calculations later + I_0 = 100e9 / (10 * 10) # 100 GW divided over 100 square meters + self.I_0 = I_0 # needed for calculations later mu_x = 0.5 mu_y = 0.5 sigma = 0.25 - self.LB_gaussian = LaserBeam(lambda x, y: I_0 * np.exp(-1/2 *((x-mu_x)/sigma)**2 - -1/2 *((y-mu_y)/sigma)**2), - lambda x,y: [0,1]) - self.LB_flat = LaserBeam(lambda x, y: np.ones(x.shape)*I_0, lambda x,y: [0,1]) - self.LB_flat_bounded = LaserBeam(laser_intensity_bounded, lambda x,y: [0,1]) + self.LB_gaussian = LaserBeam( + lambda x, y: I_0 + * np.exp( + -1 / 2 * ((x - mu_x) / sigma) ** 2 - 1 / 2 * ((y - mu_y) / sigma) ** 2 + ), + lambda x, y: [0, 1], + ) + self.LB_flat = LaserBeam( + lambda x, y: np.ones(x.shape) * I_0, lambda x, y: [0, 1] + ) + self.LB_flat_bounded = LaserBeam(laser_intensity_bounded, lambda x, y: [0, 1]) # Set up general purpose OpticalForceCalculator for particle in self.PS.particles: particle.optical_type = ParticleOpticalPropertyType.SPECULAR - + self.OpticalForces = OpticalForceCalculator(self.PS, self.LB_flat) - - - + def test_specular_flat(self): - expected_force = 2*self.I_0 / c + expected_force = 2 * self.I_0 / c OpticalForces = self.OpticalForces forces = OpticalForces.force_value() net_force = sum(np.linalg.norm(forces, axis=1)) self.assertAlmostEqual(net_force, expected_force) - + def test_specular_45deg_x(self): - expected_force =2*self.I_0 / c + expected_force = 2 * self.I_0 / c OpticalForces = self.OpticalForces - - # Rotate the PS, calculate forces, put it back - OpticalForces.displace_particle_system([0,0,0,45,0,0]) + + # Rotate the PS, calculate forces, put it back + OpticalForces.displace_particle_system([0, 0, 0, 45, 0, 0]) forces = OpticalForces.force_value() net_force = sum(np.linalg.norm(forces, axis=1)) OpticalForces.un_displace_particle_system() - + # Calculate summary numbers - vertical_force = sum(forces[:,2]) - horizontal_force = sum(forces[:,1]) - null_force = sum(forces[:,0]) - + vertical_force = sum(forces[:, 2]) + horizontal_force = sum(forces[:, 1]) + null_force = sum(forces[:, 0]) + # Perform checks # Scale by 1/sqrt(2) to account for cosine factor due to rotation with self.subTest(i=0): - self.assertAlmostEqual(net_force, expected_force/ np.sqrt(2)) - + self.assertAlmostEqual(net_force, expected_force / np.sqrt(2)) + # Scale twice by 1/sqrt(2), once for cosine factor, once for vertical component - with self.subTest(i=1): - self.assertAlmostEqual(vertical_force,expected_force/2) - + with self.subTest(i=1): + self.assertAlmostEqual(vertical_force, expected_force / 2) + # Misc. checks - with self.subTest(i=2): + with self.subTest(i=2): self.assertAlmostEqual(abs(vertical_force), abs(horizontal_force)) - with self.subTest(i=3): - self.assertEqual(null_force,0) - + with self.subTest(i=3): + self.assertEqual(null_force, 0) + def test_specular_45deg_y(self): - expected_force = 2*self.I_0 / c + expected_force = 2 * self.I_0 / c OpticalForces = self.OpticalForces - - # Rotate the PS, calculate forces, put it back - OpticalForces.displace_particle_system([0,0,0,0,45,0]) + + # Rotate the PS, calculate forces, put it back + OpticalForces.displace_particle_system([0, 0, 0, 0, 45, 0]) forces = OpticalForces.force_value() net_force = sum(np.linalg.norm(forces, axis=1)) OpticalForces.un_displace_particle_system() - + # Calculate summary numbers - vertical_force = sum(forces[:,2]) - horizontal_force = sum(forces[:,0]) - null_force = sum(forces[:,1]) - + vertical_force = sum(forces[:, 2]) + horizontal_force = sum(forces[:, 0]) + null_force = sum(forces[:, 1]) + # Perform checks # Scale by 1/sqrt(2) to account for cosine factor due to rotation - with self.subTest(i=0): - self.assertAlmostEqual(net_force, expected_force/ np.sqrt(2)) - + with self.subTest(i=0): + self.assertAlmostEqual(net_force, expected_force / np.sqrt(2)) + # Scale twice by 1/sqrt(2), once for cosine factor, once for vertical component - with self.subTest(i=1): - self.assertAlmostEqual(vertical_force,expected_force/2) - + with self.subTest(i=1): + self.assertAlmostEqual(vertical_force, expected_force / 2) + # Misc. checks - with self.subTest(i=2): + with self.subTest(i=2): self.assertAlmostEqual(abs(vertical_force), abs(horizontal_force)) - with self.subTest(i=3): - self.assertEqual(null_force,0) + with self.subTest(i=3): + self.assertEqual(null_force, 0) - def test_specular_pyramid(self): - expected_force = 2*self.I_0 / c + expected_force = 2 * self.I_0 / c OpticalForces = self.OpticalForces # Stretch PS into pyramid shape, which incrases surface area by sqrt(2) height = lambda x, y: min(0.5 - abs(x - 0.5), 0.5 - abs(y - 0.5)) for particle in self.PS.particles: - x,y,_ = particle.x - z = height(x,y) - particle.x[2]= z - + x, y, _ = particle.x + z = height(x, y) + particle.x[2] = z + forces = OpticalForces.force_value() - - # Net force should be increased due to increase in surface area, but + + # Net force should be increased due to increase in surface area, but # decreased due to angle of incidence, effects should cancel net_force = sum(np.linalg.norm(forces, axis=1)) - + # There is a know discrepancy due to an error in the surface calculation # Therefore, these are tested to be less than two percent different - with self.subTest(i=0): - #self.assertAlmostEqual(net_force, expected_force, places = 2) - self.assertGreater(0.02, abs(net_force/expected_force-1)) + with self.subTest(i=0): + # self.assertAlmostEqual(net_force, expected_force, places = 2) + self.assertGreater(0.02, abs(net_force / expected_force - 1)) - vertical_force = sum(forces[:,2]) + vertical_force = sum(forces[:, 2]) with self.subTest(i=1): - expected_vertical_force = expected_force/ np.sqrt(2) - self.assertGreater(0.02, abs(vertical_force/expected_vertical_force-1)) + expected_vertical_force = expected_force / np.sqrt(2) + self.assertGreater(0.02, abs(vertical_force / expected_vertical_force - 1)) - def test_axicon_flat(self): for i in range(4): with self.subTest(i=i): - directing_angle = i*15 + directing_angle = i * 15 directing_angle_rad = np.deg2rad(directing_angle) - + force_absorption = np.array([0, self.I_0 / c]) - force_emission = np.array([np.sin(directing_angle_rad*2), - np.cos(directing_angle_rad*2)] - )*self.I_0 / c - expected_force = np.linalg.norm(force_absorption+force_emission) - + force_emission = ( + np.array( + [ + np.sin(directing_angle_rad * 2), + np.cos(directing_angle_rad * 2), + ] + ) + * self.I_0 + / c + ) + expected_force = np.linalg.norm(force_absorption + force_emission) + for particle in self.PS.particles: particle.optical_type = ParticleOpticalPropertyType.AXICONGRATING - particle.axicon_angle = Rotation.from_euler('y', directing_angle, degrees=True).as_matrix() - - + particle.axicon_angle = Rotation.from_euler( + "y", directing_angle, degrees=True + ).as_matrix() + OpticalForces = OpticalForceCalculator(self.PS, self.LB_flat) forces = OpticalForces.force_value() net_force = sum(np.linalg.norm(forces, axis=1)) @@ -187,243 +203,259 @@ def test_axicon_flat(self): def test_find_center_of_mass(self): expected_COM = [0.5, 0.5, 0] OpticalForces = self.OpticalForces - + COM = OpticalForces.find_center_of_mass() for i in range(3): with self.subTest(i=i): self.assertAlmostEqual(COM[i], expected_COM[i]) - + # Now we inflate and re-test - expected_COM = [0.5, 0.5, 0.5/3] + expected_COM = [0.5, 0.5, 0.5 / 3] height = lambda x, y: min(0.5 - abs(x - 0.5), 0.5 - abs(y - 0.5)) for particle in self.PS.particles: - x,y,_ = particle.x - z = height(x,y) - particle.x[2]= z - + x, y, _ = particle.x + z = height(x, y) + particle.x[2] = z + COM = OpticalForces.find_center_of_mass() for i in range(3): - with self.subTest(i=i+3): - self.assertAlmostEqual(COM[i], expected_COM[i], places = 3) - + with self.subTest(i=i + 3): + self.assertAlmostEqual(COM[i], expected_COM[i], places=3) + def test_translate_mesh(self): OpticalForces = OpticalForceCalculator(self.PS, self.LB_flat) - mesh = OpticalForces.translate_mesh(np.zeros([10,3]), np.ones([3])) - - self.assertEqual(np.all(mesh == np.ones([10,3])), True) - + mesh = OpticalForces.translate_mesh(np.zeros([10, 3]), np.ones([3])) + + self.assertEqual(np.all(mesh == np.ones([10, 3])), True) + def test_rotate_mesh(self): OpticalForces = self.OpticalForces - mesh = np.meshgrid(np.linspace(0, 1, 6), - np.linspace(0, 1, 6)) - mesh = np.column_stack(list(zip(mesh[0],mesh[1]))).T - mesh = np.column_stack((mesh,np.zeros(len(mesh)).T)) - + mesh = np.meshgrid(np.linspace(0, 1, 6), np.linspace(0, 1, 6)) + mesh = np.column_stack(list(zip(mesh[0], mesh[1]))).T + mesh = np.column_stack((mesh, np.zeros(len(mesh)).T)) + original = mesh.copy() - + for i in range(3): - mesh = OpticalForces.rotate_mesh(mesh, [0,0,120]) - - diff = np.sum(mesh-original) + mesh = OpticalForces.rotate_mesh(mesh, [0, 0, 120]) + + diff = np.sum(mesh - original) self.assertAlmostEqual(diff, 0) - + def test_displace_particle_system(self): OpticalForces = self.OpticalForces - + with self.subTest(i=0): expected_COM = np.array([0.5, 0.5, 0]) + 1 - OpticalForces.displace_particle_system([1,1,1,0,0,0]) + OpticalForces.displace_particle_system([1, 1, 1, 0, 0, 0]) COM = OpticalForces.find_center_of_mass() OpticalForces.un_displace_particle_system() - self.assertAlmostEqual(np.sum(expected_COM-COM), 0) - + self.assertAlmostEqual(np.sum(expected_COM - COM), 0) + with self.subTest(i=1): - x_0,_ = self.PS.x_v_current_3D + x_0, _ = self.PS.x_v_current_3D logging.warning("Two displacement warnings expected as part of test") for i in range(3): # this function will issue a warning when called the second and third times - OpticalForces.displace_particle_system([0,0,0,0,0,120]) - self.PS.current_displacement = [0,0,0,0,0,0] - x_1,_ = self.PS.x_v_current_3D + OpticalForces.displace_particle_system([0, 0, 0, 0, 0, 120]) + self.PS.current_displacement = [0, 0, 0, 0, 0, 0] + x_1, _ = self.PS.x_v_current_3D error = x_1 - x_0 self.assertAlmostEqual(np.sum(error), 0) - + def test_un_displace_particle_system(self): OpticalForces = self.OpticalForces expected_COM = np.array([0.5, 0.5, 0]) - OpticalForces.displace_particle_system([0,0,0,45,45,45]) + OpticalForces.displace_particle_system([0, 0, 0, 45, 45, 45]) COM = OpticalForces.find_center_of_mass() OpticalForces.un_displace_particle_system() - self.assertAlmostEqual(np.sum(expected_COM-COM), 0) - - + self.assertAlmostEqual(np.sum(expected_COM - COM), 0) + def test_stability_unit_mesh_uniform_beam(self): # For context, J maps changes in x,y,z,rx,ry,rz to changes in forces in those DOF's - expected_force = 2*self.I_0 / c - + expected_force = 2 * self.I_0 / c + J = self.OpticalForces.calculate_stability_coefficients() with self.subTest(i=0): # Check that translations do not affect other translations nor rotations - self.assertTrue(np.allclose(J[0:6,0:3], np.zeros((6,3)))) - with self.subTest(i=1): + self.assertTrue(np.allclose(J[0:6, 0:3], np.zeros((6, 3)))) + with self.subTest(i=1): # Check that rotating about x does that same as about y - self.assertAlmostEqual(J[1,3],-J[0,4]) - with self.subTest(i=2): + self.assertAlmostEqual(J[1, 3], -J[0, 4]) + with self.subTest(i=2): # Check that rotating about x and y gives right force gradient in x or y - # If the plane is rotated 5 deg then the horizontal component of the foce will be - force_horizontal = expected_force*np.sin(np.deg2rad(5))*np.cos(np.deg2rad(5)) - d_force_horizontal_d_deg = force_horizontal /5 - self.assertAlmostEqual(J[1,3],-d_force_horizontal_d_deg) - with self.subTest(i=3): + # If the plane is rotated 5 deg then the horizontal component of the foce will be + force_horizontal = ( + expected_force * np.sin(np.deg2rad(5)) * np.cos(np.deg2rad(5)) + ) + d_force_horizontal_d_deg = force_horizontal / 5 + self.assertAlmostEqual(J[1, 3], -d_force_horizontal_d_deg) + with self.subTest(i=3): # Check that rotating about x and y gives right force gradient in z # Cos is squared, once for cosine factor, once for angle change of force - d_force_vertical_d_deg = -expected_force*(1-np.cos(np.deg2rad(5))**2)/5 - self.assertAlmostEqual(J[2,3],d_force_vertical_d_deg) - with self.subTest(i=4): + d_force_vertical_d_deg = ( + -expected_force * (1 - np.cos(np.deg2rad(5)) ** 2) / 5 + ) + self.assertAlmostEqual(J[2, 3], d_force_vertical_d_deg) + with self.subTest(i=4): # Check that rotating about x or y has same effect on f_z - self.assertAlmostEqual(J[2,3],J[2,4]) - with self.subTest(i=5): + self.assertAlmostEqual(J[2, 3], J[2, 4]) + with self.subTest(i=5): # Check that rotation about z affects nothing - self.assertTrue(np.allclose(J[:,5], np.zeros(6))) - with self.subTest(i=6): + self.assertTrue(np.allclose(J[:, 5], np.zeros(6))) + with self.subTest(i=6): # Check that rotations about x or y don't affect any other rotations - self.assertTrue(np.allclose(J[3:6,3:5], np.zeros((3,2)))) - + self.assertTrue(np.allclose(J[3:6, 3:5], np.zeros((3, 2)))) + def test_stability_unit_mesh_bounded_uniform_beam(self): # For context, J maps changes in x,y,z,rx,ry,rz to changes in forces in those DOF's - expected_force = 2*self.I_0 / c - + expected_force = 2 * self.I_0 / c + OpticalForces = OpticalForceCalculator(self.PS, self.LB_flat_bounded) J = OpticalForces.calculate_stability_coefficients() with self.subTest(i=0): # Check that translations in x,y,z do not affect forces in x,y - self.assertTrue(np.allclose(J[0:2,0:3], np.zeros((2,3)))) - with self.subTest(i=1): + self.assertTrue(np.allclose(J[0:2, 0:3], np.zeros((2, 3)))) + with self.subTest(i=1): # Check that rotating about x does that same as about y - self.assertAlmostEqual(J[1,3],-J[0,4]) - with self.subTest(i=2): + self.assertAlmostEqual(J[1, 3], -J[0, 4]) + with self.subTest(i=2): # Check that rotating about x or y gives right force gradient in x or y - # If the plane is rotated 5 deg then the horizontal component of the foce will be - force_horizontal = expected_force*np.sin(np.deg2rad(5))*np.cos(np.deg2rad(5)) - d_force_horizontal_d_deg = force_horizontal /5 - self.assertAlmostEqual(J[1,3],-d_force_horizontal_d_deg) - with self.subTest(i=3): + # If the plane is rotated 5 deg then the horizontal component of the foce will be + force_horizontal = ( + expected_force * np.sin(np.deg2rad(5)) * np.cos(np.deg2rad(5)) + ) + d_force_horizontal_d_deg = force_horizontal / 5 + self.assertAlmostEqual(J[1, 3], -d_force_horizontal_d_deg) + with self.subTest(i=3): # Check that rotating about x and y gives right force gradient in z # Cos is squared, once for cosine factor, once for angle change of force - d_force_vertical_d_deg = -expected_force*(1-np.cos(np.deg2rad(5))**2)/5 - self.assertAlmostEqual(J[2,3],d_force_vertical_d_deg) - with self.subTest(i=4): + d_force_vertical_d_deg = ( + -expected_force * (1 - np.cos(np.deg2rad(5)) ** 2) / 5 + ) + self.assertAlmostEqual(J[2, 3], d_force_vertical_d_deg) + with self.subTest(i=4): # Check that rotating about x or y has same effect on f_z - self.assertAlmostEqual(J[2,3],J[2,4]) - with self.subTest(i=5): + self.assertAlmostEqual(J[2, 3], J[2, 4]) + with self.subTest(i=5): # Check that rotation about z affects nothing except force in z mask = np.zeros(6) - mask[2] = 1 - mask = mask==0 - self.assertTrue(np.allclose(J[:,5][mask], np.zeros(5))) - with self.subTest(i=6): + mask[2] = 1 + mask = mask == 0 + self.assertTrue(np.allclose(J[:, 5][mask], np.zeros(5))) + with self.subTest(i=6): # Check that rotations about x or y don't affect any other rotations - self.assertTrue(np.allclose(J[3:6,3:5], np.zeros((3,2)))) + self.assertTrue(np.allclose(J[3:6, 3:5], np.zeros((3, 2)))) with self.subTest(i=7): # Check that translations in x or y affect the force in z correctly - force_vertical = expected_force*(1-0.1) # displacement 0.1 out of width 1 - self.assertAlmostEqual(J[2,0],J[2,1]) + force_vertical = expected_force * ( + 1 - 0.1 + ) # displacement 0.1 out of width 1 + self.assertAlmostEqual(J[2, 0], J[2, 1]) with self.subTest(i=8): # and that they're equal - self.assertAlmostEqual(J[2,0],J[2,1]) + self.assertAlmostEqual(J[2, 0], J[2, 1]) with self.subTest(i=9): # Check that translation in x or y results in correct moment about y or x # inbalance of 10% of the force on the outer 10% of the sail - moment_arm = 0.1/2 + 0.4 # distance to center of strip + distance strip+COM - force = expected_force*0.1 - moment = force*moment_arm - d_moment_d_translation = moment/0.1 - self.assertAlmostEqual(J[4,0],-d_moment_d_translation) + moment_arm = ( + 0.1 / 2 + 0.4 + ) # distance to center of strip + distance strip+COM + force = expected_force * 0.1 + moment = force * moment_arm + d_moment_d_translation = moment / 0.1 + self.assertAlmostEqual(J[4, 0], -d_moment_d_translation) with self.subTest(i=10): # and that they're equal - self.assertAlmostEqual(J[4,0],-J[3,1]) - + self.assertAlmostEqual(J[4, 0], -J[3, 1]) + def test_stability_unit_mesh_gaussian_beam(self): # Goal is to compare this against gaussian beam analytical results # Computations made with wolfram alpha - + # Set up the gaussian beam system OpticalForces = OpticalForceCalculator(self.PS, self.LB_gaussian) J = OpticalForces.calculate_stability_coefficients() with self.subTest(i=0): # To start out check if the forces are correct - # integral of gaussian beam with constants supplied above is + # integral of gaussian beam with constants supplied above is P_gauss = self.I_0 * 0.35776 - f_expected = 2* P_gauss/c - + f_expected = 2 * P_gauss / c + forces = OpticalForces.force_value() - f_abs = np.linalg.norm(forces,axis=1) + f_abs = np.linalg.norm(forces, axis=1) f_total = np.sum(f_abs) - self.assertAlmostEqual(f_total, f_expected,places=3) - + self.assertAlmostEqual(f_total, f_expected, places=3) + # Check effect of displacement # integral of gaussian over shifted square: f_unshifted = f_total P_gauss = self.I_0 * 0.351218 - f_expected = 2* P_gauss/c - + f_expected = 2 * P_gauss / c + forces = OpticalForces.force_value() - f_abs = np.linalg.norm(forces,axis=1) + f_abs = np.linalg.norm(forces, axis=1) with self.subTest(i=1): # change in z force - d_f_d_x = (f_expected-f_unshifted)/0.1 - self.assertAlmostEqual(J[2,0], d_f_d_x,places = 3) + d_f_d_x = (f_expected - f_unshifted) / 0.1 + self.assertAlmostEqual(J[2, 0], d_f_d_x, places=3) def test_calculate_restoring_forces(self): # testing if forces for positive and negative displacements are the same # First translation, then rotation OpticalForces = OpticalForceCalculator(self.PS, self.LB_gaussian) - - displacement = np.array([0.1,0,0,0,0,0]) - #Positive translation + + displacement = np.array([0.1, 0, 0, 0, 0, 0]) + # Positive translation OpticalForces.displace_particle_system(displacement) net_force_pos, net_moments_pos = OpticalForces.calculate_restoring_forces() OpticalForces.un_displace_particle_system() - - #Negative translation + + # Negative translation OpticalForces.displace_particle_system(-displacement) net_force_neg, net_moments_neg = OpticalForces.calculate_restoring_forces() OpticalForces.un_displace_particle_system() - + with self.subTest(i=0): - self.assertTrue(np.allclose(net_force_pos,net_force_neg)) + self.assertTrue(np.allclose(net_force_pos, net_force_neg)) with self.subTest(i=1): - self.assertTrue(np.allclose(np.abs(net_moments_pos),np.abs(net_moments_neg))) - - displacement = np.array([0,0,0,5,0,0]) - #Positive rotation + self.assertTrue( + np.allclose(np.abs(net_moments_pos), np.abs(net_moments_neg)) + ) + + displacement = np.array([0, 0, 0, 5, 0, 0]) + # Positive rotation OpticalForces.displace_particle_system(displacement) net_force_pos, net_moments_pos = OpticalForces.calculate_restoring_forces() OpticalForces.un_displace_particle_system() - - #Negative rotation + + # Negative rotation OpticalForces.displace_particle_system(-displacement) net_force_neg, net_moments_neg = OpticalForces.calculate_restoring_forces() OpticalForces.un_displace_particle_system() - + with self.subTest(i=2): - self.assertTrue(np.allclose(np.abs(net_force_pos),np.abs(net_force_neg))) + self.assertTrue(np.allclose(np.abs(net_force_pos), np.abs(net_force_neg))) with self.subTest(i=3): - self.assertTrue(np.allclose(np.abs(net_moments_pos),np.abs(net_moments_neg))) + self.assertTrue( + np.allclose(np.abs(net_moments_pos), np.abs(net_moments_neg)) + ) + -def laser_intensity_bounded(x,y): - I_0 = 100e9 /(10*10) +def laser_intensity_bounded(x, y): + I_0 = 100e9 / (10 * 10) intensity = np.zeros(x.shape) - intensity[(x>-0.001)*(x<1.001)*(y>-0.001)*(y<1.001)] = I_0 - return intensity + intensity[(x > -0.001) * (x < 1.001) * (y > -0.001) * (y < 1.001)] = I_0 + return intensity -if __name__ == '__main__': + +if __name__ == "__main__": unittest.main() - + debug = False if debug: Logger = logging.getLogger() - Logger.setLevel(10) # Use Logger.setLevel(30) to set it back to default + Logger.setLevel(10) # Use Logger.setLevel(30) to set it back to default T = TestOpticalForceCalculator() T.setUp() diff --git a/test/src/particleSystem/TestParticleSystem.py b/test/src/particleSystem/TestParticleSystem.py index 31f1eaf..6d496ee 100644 --- a/test/src/particleSystem/TestParticleSystem.py +++ b/test/src/particleSystem/TestParticleSystem.py @@ -9,10 +9,10 @@ import numpy as np -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import SimulateTripleChainWithMass +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import SimulateTripleChainWithMass from scipy.spatial.transform import Rotation -import Particle_System_Simulator.Mesh.mesh_functions as MF +import PSS.Mesh.mesh_functions as MF class TestParticleSystem(unittest.TestCase): diff --git a/test/src/particleSystem/test_ParticleSystem.py b/test/src/particleSystem/test_ParticleSystem.py index 6c6fa08..61f84d5 100644 --- a/test/src/particleSystem/test_ParticleSystem.py +++ b/test/src/particleSystem/test_ParticleSystem.py @@ -2,9 +2,9 @@ from scipy.spatial.transform import Rotation import pytest -from Particle_System_Simulator.particleSystem.ParticleSystem import ParticleSystem -from Particle_System_Simulator.Sim.simulations import SimulateTripleChainWithMass -import Particle_System_Simulator.Mesh.mesh_functions as MF +from PSS.particleSystem.ParticleSystem import ParticleSystem +from PSS.Sim.simulations import SimulateTripleChainWithMass +import PSS.Mesh.mesh_functions as MF @pytest.fixture