From 0e1f64c6de34ac361533b22bd5bcb907782621e0 Mon Sep 17 00:00:00 2001 From: LyceanEM Date: Wed, 4 Dec 2024 12:16:23 +0000 Subject: [PATCH] Reformatting with Black --- lyceanem/__init__.py | 1 + lyceanem/base_classes.py | 33 +++--- lyceanem/electromagnetics/beamforming.py | 6 +- lyceanem/electromagnetics/emfunctions.py | 131 +++++++++++---------- lyceanem/electromagnetics/empropagation.py | 21 +++- lyceanem/geometry/geometryfunctions.py | 30 ++--- lyceanem/geometry/targets.py | 69 +++++++---- lyceanem/models/frequency_domain.py | 62 ++++++---- lyceanem/raycasting/rayfunctions.py | 24 ++-- lyceanem/tests/reflectordata.py | 14 +-- lyceanem/utility/math_functions.py | 12 +- lyceanem/utility/mesh_functions.py | 48 ++++---- 12 files changed, 256 insertions(+), 195 deletions(-) diff --git a/lyceanem/__init__.py b/lyceanem/__init__.py index 380b78e..2611834 100644 --- a/lyceanem/__init__.py +++ b/lyceanem/__init__.py @@ -1,4 +1,5 @@ """LyceanEM : Electromagnetics Modelling for Antenna and Antenna Array Development on Complex Platforms""" + try: from importlib import metadata except ImportError: # for Python<3.8 diff --git a/lyceanem/base_classes.py b/lyceanem/base_classes.py index ea2e56d..9847bc9 100644 --- a/lyceanem/base_classes.py +++ b/lyceanem/base_classes.py @@ -389,7 +389,7 @@ def triangles_base_raycaster(self): """ triangles = np.empty((0), dtype=base_types.triangle_t) for item in range(len(self.solids)): - #temp_object = copy.deepcopy(self.solids[item]) + # temp_object = copy.deepcopy(self.solids[item]) temp_object = GF.mesh_transform(self.solids[item], self.pose, False) triangles = np.append(triangles, RF.convertTriangles(temp_object)) @@ -429,15 +429,16 @@ def export_all_points(self, point_index=None): return point_cloud - def excitation_function(self, - desired_e_vector, - point_index=None, - phase_shift="none", - wavelength=1.0, - steering_vector=np.zeros((1, 3)), - transmit_power=1.0, - local_projection=True - ): + def excitation_function( + self, + desired_e_vector, + point_index=None, + phase_shift="none", + wavelength=1.0, + steering_vector=np.zeros((1, 3)), + transmit_power=1.0, + local_projection=True, + ): # generate the local excitation function and then convert into the global coordinate frame. if point_index == None: aperture_points = self.export_all_points() @@ -446,23 +447,27 @@ def excitation_function(self, from .electromagnetics.emfunctions import excitation_function - excitation_weights=excitation_function( + excitation_weights = excitation_function( aperture_points, desired_e_vector, phase_shift=phase_shift, wavelength=wavelength, steering_vector=steering_vector, transmit_power=transmit_power, - local_projection=local_projection) + local_projection=local_projection, + ) return excitation_weights + def export_all_structures(self): - #objects = copy.deepcopy(self.structures.solids) + # objects = copy.deepcopy(self.structures.solids) for item in range(len(self.structures.solids)): if self.structures.solids[item] is None: print("Structure does not exist") else: - self.structures.solids[item] = GF.mesh_transform(self.structures.solids[item], self.pose, False) + self.structures.solids[item] = GF.mesh_transform( + self.structures.solids[item], self.pose, False + ) return self.structures.solids diff --git a/lyceanem/electromagnetics/beamforming.py b/lyceanem/electromagnetics/beamforming.py index b426396..2063b36 100644 --- a/lyceanem/electromagnetics/beamforming.py +++ b/lyceanem/electromagnetics/beamforming.py @@ -443,7 +443,11 @@ def MaximumDirectivityMap( ) WS_weights = WavefrontWeights(source_points, steering_vector, wavelength) EGC_weights = EGCWeights( - Etheta.astype(np.complex64), Ephi.astype(np.complex64), command_angles, az_range=az_range, elev_range=elev_range + Etheta.astype(np.complex64), + Ephi.astype(np.complex64), + command_angles, + az_range=az_range, + elev_range=elev_range, ) EGC_weights2 = EGCWeights( Etheta.astype(np.complex64), diff --git a/lyceanem/electromagnetics/emfunctions.py b/lyceanem/electromagnetics/emfunctions.py index 1b5c2fe..ae854c3 100644 --- a/lyceanem/electromagnetics/emfunctions.py +++ b/lyceanem/electromagnetics/emfunctions.py @@ -3,37 +3,39 @@ def excitation_function( - aperture_points, - desired_e_vector, - phase_shift="none", - wavelength=1.0, - steering_vector=np.zeros((1, 3)), - transmit_power=1.0, - local_projection=True + aperture_points, + desired_e_vector, + phase_shift="none", + wavelength=1.0, + steering_vector=np.zeros((1, 3)), + transmit_power=1.0, + local_projection=True, ): """ Calculate the excitation function for the given aperture points, desired electric field vector, phase shift, wavelength, steering vector, transmit power, and local projection. This will generate the normalised field intensities for the desired total transmit power, and beamforming algorithm. The aperture mesh should have Normals and Area associated with each point for full functionality. """ - if local_projection: # as export all points imposes the transformation from local to global frame on the points and associated normal vectors, no rotation is required within calculate_conformalVectors from .empropagation import calculate_conformalVectors + aperture_weights = calculate_conformalVectors( desired_e_vector, aperture_points.point_data["Normals"], np.eye(3) ) else: - aperture_weights = np.repeat(desired_e_vector, aperture_points.points.shape[0], axis=0) + aperture_weights = np.repeat( + desired_e_vector, aperture_points.points.shape[0], axis=0 + ) if phase_shift == "wavefront": from .beamforming import WavefrontWeights + source_points = np.asarray(aperture_points.points) - phase_weights = WavefrontWeights( - source_points, steering_vector, wavelength - ) + phase_weights = WavefrontWeights(source_points, steering_vector, wavelength) aperture_weights = aperture_weights * phase_weights.reshape(-1, 1) if phase_shift == "coherence": from .beamforming import ArbitaryCoherenceWeights + source_points = np.asarray(aperture_points.points) phase_weights = ArbitaryCoherenceWeights( source_points, steering_vector, wavelength @@ -52,6 +54,8 @@ def excitation_function( aperture_weights, areas, transmit_power ) return calibrated_amplitude_weights + + def fresnel_zone(pointA, pointB, wavelength, zone=1): # based on the provided points, wavelength, and zone number, calculate the fresnel zone. This is defined as an ellipsoid for which the difference in distance between the line AB (line of sight), and AP+PB (reflected wave) is a constant multiple of ($n\dfrac{\lambda}{2}$). foci = np.stack([pointA, pointB]) @@ -90,15 +94,9 @@ def field_magnitude_phase(field_data): ) ): # Exyz exsists in the dataset - Ex = ( - field_data.point_data["Ex-Real"] + 1j * field_data.point_data["Ex-Imag"] - ) - Ey = ( - field_data.point_data["Ey-Real"] + 1j * field_data.point_data["Ey-Imag"] - ) - Ez = ( - field_data.point_data["Ez-Real"] + 1j * field_data.point_data["Ez-Imag"] - ) + Ex = field_data.point_data["Ex-Real"] + 1j * field_data.point_data["Ex-Imag"] + Ey = field_data.point_data["Ey-Real"] + 1j * field_data.point_data["Ey-Imag"] + Ez = field_data.point_data["Ez-Real"] + 1j * field_data.point_data["Ez-Imag"] field_data.point_data["Ex-Magnitude"] = np.abs(Ex) field_data.point_data["Ex-Phase"] = np.angle(Ex) field_data.point_data["Ey-Magnitude"] = np.abs(Ey) @@ -217,7 +215,7 @@ def Exyz_to_EthetaEphi(field_data): from lyceanem.geometry.geometryfunctions import theta_phi_r field_data = theta_phi_r(field_data) - + electric_fields = extract_electric_fields(field_data) theta = field_data.point_data["theta_(Radians)"] phi = field_data.point_data["phi_(Radians)"] @@ -241,12 +239,9 @@ def Exyz_to_EthetaEphi(field_data): def field_vectors(field_data): fields = np.array( [ - field_data.point_data["Ex-Real"] - + 1j * field_data.point_data["Ex-Imag"], - field_data.point_data["Ey-Real"] - + 1j * field_data.point_data["Ey-Imag"], - field_data.point_data["Ez-Real"] - + 1j * field_data.point_data["Ez-Imag"], + field_data.point_data["Ex-Real"] + 1j * field_data.point_data["Ex-Imag"], + field_data.point_data["Ey-Real"] + 1j * field_data.point_data["Ey-Imag"], + field_data.point_data["Ez-Real"] + 1j * field_data.point_data["Ez-Imag"], ] ).transpose() directions = np.abs(fields) @@ -267,20 +262,23 @@ def transform_em(field_data, r): """ if all( - k in field_data.point_data.keys() - for k in ("Ex-Real", "Ey-Real", "Ez-Real") + k in field_data.point_data.keys() for k in ("Ex-Real", "Ey-Real", "Ez-Real") ): # print("Rotating Ex,Ey,Ez") - fields = np.array( - [ - field_data.point_data["Ex-Real"] - + 1j * field_data.point_data["Ex-Imag"], - field_data.point_data["Ey-Real"] - + 1j * field_data.point_data["Ey-Imag"], - field_data.point_data["Ez-Real"] - + 1j * field_data.point_data["Ez-Imag"], - ] - ).reshape(3,-1).transpose() + fields = ( + np.array( + [ + field_data.point_data["Ex-Real"] + + 1j * field_data.point_data["Ex-Imag"], + field_data.point_data["Ey-Real"] + + 1j * field_data.point_data["Ey-Imag"], + field_data.point_data["Ez-Real"] + + 1j * field_data.point_data["Ez-Imag"], + ] + ) + .reshape(3, -1) + .transpose() + ) rot_fields = r.apply(fields) field_data.point_data["Ex-Real"] = np.real( rot_fields[:, 0] @@ -299,16 +297,20 @@ def transform_em(field_data, r): field_data = theta_phi_r(field_data) field_data = EthetaEphi_to_Exyz(field_data) - fields = np.array( - [ - field_data.point_data["Ex-Real"] - + 1j * field_data.point_data["Ex-Imag"], - field_data.point_data["Ey-Real"] - + 1j * field_data.point_data["Ey-Imag"], - field_data.point_data["Ez-Real"] - + 1j * field_data.point_data["Ez-Imag"], - ] - ).reshape(3,-1).transpose() + fields = ( + np.array( + [ + field_data.point_data["Ex-Real"] + + 1j * field_data.point_data["Ex-Imag"], + field_data.point_data["Ey-Real"] + + 1j * field_data.point_data["Ey-Imag"], + field_data.point_data["Ez-Real"] + + 1j * field_data.point_data["Ez-Imag"], + ] + ) + .reshape(3, -1) + .transpose() + ) rot_fields = r.apply(fields) field_data.point_data["Ex-Real"] = np.real(rot_fields[:, 0]) field_data.point_data["Ey-Real"] = np.real(rot_fields[:, 1]) @@ -335,7 +337,7 @@ def update_electric_fields(field_data, ex, ey, ez): return field_data -def PoyntingVector(field_data,measurement=False,aperture=None): +def PoyntingVector(field_data, measurement=False, aperture=None): """ Calculate the poynting vector for the given field data. If the magnetic field data is present, then the poynting vector is calculated using the cross product of the electric and magnetic field vectors. If the magnetic field data is not present, then the poynting vector is calculated using the electric field vectors and the impedance of the material. If material parameters are not included in the field data, then the impedance is assumed to be that of free space. @@ -384,8 +386,7 @@ def PoyntingVector(field_data,measurement=False,aperture=None): ] ).transpose() if all( - k in field_data.point_data.keys() - for k in ("Hx-Real", "Hy-Real", "Hz-Real") + k in field_data.point_data.keys() for k in ("Hx-Real", "Hy-Real", "Hz-Real") ): # magnetic field data present, so use magnetic_field_vectors = np.array( @@ -410,7 +411,7 @@ def PoyntingVector(field_data,measurement=False,aperture=None): ).reshape(-1, 1) if measurement: - poynting_vector_complex = poynting_vector_complex/aperture + poynting_vector_complex = poynting_vector_complex / aperture field_data.point_data["Poynting_Vector_(Magnitude)"] = np.zeros( (field_data.points.shape[0], 1) @@ -444,10 +445,9 @@ def Directivity(field_data): # E(theta) and E(phi) are not present in the data field_data = Exyz_to_EthetaEphi(field_data) - if not all( - k in field_data.point_data.keys() for k in ("Area") - ): + if not all(k in field_data.point_data.keys() for k in ("Area")): from lyceanem.geometry.geometryfunctions import compute_areas + # E(theta) and E(phi) are not present in the data field_data = compute_areas(field_data) @@ -465,17 +465,20 @@ def Directivity(field_data): ) ** 2 ) - #Calculate Solid Angle - solid_angle=field_data.point_data["Area"]/field_data.point_data["Radial_Distance_(m)"]**2 + # Calculate Solid Angle + solid_angle = ( + field_data.point_data["Area"] + / field_data.point_data["Radial_Distance_(m)"] ** 2 + ) Utotal = Utheta + Uphi - Uav=(np.nansum(Utotal.ravel()*solid_angle)/(4*np.pi)) - #sin_factor = np.abs( + Uav = np.nansum(Utotal.ravel() * solid_angle) / (4 * np.pi) + # sin_factor = np.abs( # np.sin(field_data.point_data["theta (Radians)"]) - #).reshape(-1,1) # only want positive factor - #power_sum = np.sum(np.abs(Utheta * sin_factor)) + np.sum(np.abs(Uphi * sin_factor)) + # ).reshape(-1,1) # only want positive factor + # power_sum = np.sum(np.abs(Utheta * sin_factor)) + np.sum(np.abs(Uphi * sin_factor)) # need to dynamically account for the average contribution of each point, this is only true for a theta step of 1 degree, and phi step of 10 degrees - #Uav = (power_sum * (np.radians(1.0) * np.radians(10.0))) / (4 * np.pi) + # Uav = (power_sum * (np.radians(1.0) * np.radians(10.0))) / (4 * np.pi) Dtheta = Utheta / Uav Dphi = Uphi / Uav Dtot = Utotal / Uav diff --git a/lyceanem/electromagnetics/empropagation.py b/lyceanem/electromagnetics/empropagation.py index 6f4bc2b..7c5bc1e 100644 --- a/lyceanem/electromagnetics/empropagation.py +++ b/lyceanem/electromagnetics/empropagation.py @@ -10,7 +10,16 @@ import cupy as cp import numpy as np import scipy.stats -from numba import cuda, float32, float64, complex64, complex128, njit, guvectorize, complex128 +from numba import ( + cuda, + float32, + float64, + complex64, + complex128, + njit, + guvectorize, + complex128, +) from numpy.linalg import norm import lyceanem.base_types as base_types @@ -784,8 +793,8 @@ def lossy_propagation(point1, point2, alpha, beta): normal[2] = point2["nz"] projection_dot = dot_vec(outgoing_dir, normal) front = -(1 / (2 * cmath.pi)) - s=2.5 - distance_loss=1.0/((1+length[0]**s)**(1/s)) + s = 2.5 + distance_loss = 1.0 / ((1 + length[0] ** s) ** (1 / s)) G = (cmath.exp(-(alpha[0] + 1j * beta[0]) * length[0])) * distance_loss dG = (-(alpha[0] + 1j * beta[0]) - complex64((distance_loss))) * G @@ -2208,7 +2217,7 @@ def EMGPUFreqDomain( """ # ctx = cuda.current_context() # ctx.reset() - #clear GPU memory + # clear GPU memory cuda.current_context().memory_manager.deallocations.clear() free_mem, total_mem = cuda.current_context().get_memory_info() max_mem = np.ceil(free_mem).astype(np.int64) @@ -2225,9 +2234,9 @@ def EMGPUFreqDomain( # print("Number of Chunks",np.ceil(memory_requirements/max_mem).astype(int)+1) # create chunks based upon number of chunks required num_chunks = np.ceil(memory_requirements / max_mem).astype(int) + 1 - if num_chunks<0: + if num_chunks < 0: print(num_chunks) - + source_chunking = np.linspace(0, source_num, num_chunks + 1).astype(np.int32) scattering_network = np.zeros( (source_num, sink_num, 3, 2), diff --git a/lyceanem/geometry/geometryfunctions.py b/lyceanem/geometry/geometryfunctions.py index a6380b6..77e48d2 100644 --- a/lyceanem/geometry/geometryfunctions.py +++ b/lyceanem/geometry/geometryfunctions.py @@ -11,7 +11,6 @@ from ..electromagnetics.emfunctions import transform_em - def cell_centroids(field_data): """ In order to calculate the centroid of the triangle, take vertices from meshio triangle mesh, then put each triangle @@ -92,9 +91,7 @@ def mesh_rotate(mesh, rotation, rotation_centre=np.zeros((1, 3), dtype=np.float3 mesh_return.cell_data = cell_data # if field data is present, rotate fields - if all( - k in mesh.point_data.keys() for k in ("Ex-Real", "Ey-Real", "Ez-Real") - ): + if all(k in mesh.point_data.keys() for k in ("Ex-Real", "Ey-Real", "Ez-Real")): # rotate field data fields = transform_em(copy.deepcopy(mesh), r) for key in ( @@ -143,28 +140,31 @@ def mesh_transform(mesh, transform_matrix, rotate_only): )[:3] if "Normals" in mesh.cell_data: if isinstance(mesh.cell_data["Normals"], list): - #multiple cell types + # multiple cell types for i in range(len(mesh.cell_data["Normals"])): for j in range(mesh.cell_data["Normals"][i].shape[0]): return_mesh.cell_data["Normals"][i][j] = np.dot( - transform_matrix, np.append(mesh.cell_data["Normals"][i][j], 0) + transform_matrix, + np.append(mesh.cell_data["Normals"][i][j], 0), )[:3] else: - #single cell type + # single cell type for i in range(len(mesh.cell_data["Normals"])): - return_mesh.cell_data["Normals"][i] = np.dot( - transform_matrix, np.append(mesh.cell_data["Normals"][i], 0) - )[:3] + return_mesh.cell_data["Normals"][i] = np.dot( + transform_matrix, np.append(mesh.cell_data["Normals"][i], 0) + )[:3] return return_mesh -def locate_cell_index(field_data,cell_type='triangle'): + +def locate_cell_index(field_data, cell_type="triangle"): for inc, cell in enumerate(field_data.cells): - if cell.type==cell_type: - desired_index=inc - + if cell.type == cell_type: + desired_index = inc + return desired_index - + + def compute_areas(field_data): cell_areas = [] for inc, cell in enumerate(field_data.cells): diff --git a/lyceanem/geometry/targets.py b/lyceanem/geometry/targets.py index 14251c2..544915b 100755 --- a/lyceanem/geometry/targets.py +++ b/lyceanem/geometry/targets.py @@ -117,14 +117,16 @@ def rectReflector(majorsize, minorsize, thickness): pv_mesh = pv_mesh.triangulate() pv_mesh.compute_normals(inplace=True, consistent_normals=False) from ..utility.mesh_functions import pyvista_to_meshio + triangles = np.reshape(np.array(pv_mesh.faces), (12, 4)) triangles = triangles[:, 1:] - #mesh = meshio.Mesh(pv_mesh.points, {"triangle": triangles}) - mesh=pyvista_to_meshio(pv_mesh) + # mesh = meshio.Mesh(pv_mesh.points, {"triangle": triangles}) + mesh = pyvista_to_meshio(pv_mesh) from .geometryfunctions import compute_normals - mesh=compute_normals(mesh) + + mesh = compute_normals(mesh) red = np.zeros((8, 1), dtype=np.float32) green = np.ones((8, 1), dtype=np.float32) * 0.259 blue = np.ones((8, 1), dtype=np.float32) * 0.145 @@ -201,7 +203,8 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle): mesh.point_data["Normals"] = np.asarray(pv_mesh.point_normals) mesh.cell_data["Normals"] = [np.asarray(pv_mesh.cell_normals)] from .geometryfunctions import compute_areas - mesh=compute_areas(mesh) + + mesh = compute_areas(mesh) red = np.zeros((8, 1), dtype=np.float32) green = np.ones((8, 1), dtype=np.float32) * 0.259 blue = np.ones((8, 1), dtype=np.float32) * 0.145 @@ -209,7 +212,7 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle): mesh.point_data["red"] = red mesh.point_data["green"] = green mesh.point_data["blue"] = blue - + return mesh @@ -324,7 +327,9 @@ def parabola(x): point_num = 15 x_pos = np.linspace(-0.5 * diameter, 0.5 * diameter, point_num) z_pos = parabola(x_pos) - coords = np.array([x_pos.ravel(), np.zeros((point_num)), z_pos.ravel()]).transpose() + coords = np.array( + [x_pos.ravel(), np.zeros((point_num)), z_pos.ravel()] + ).transpose() points_list = [] for inc in range(point_num): points_list.append(geom.add_point(coords[inc, :].tolist())) @@ -457,7 +462,8 @@ def parabola(x): return mesh, aperture_points -def spherical_field(az_range,elev_range,outward_normals=False,field_radius=1.0): + +def spherical_field(az_range, elev_range, outward_normals=False, field_radius=1.0): """ Create a spherical field of points, with normals, a triangle mesh, and areas calculated for each triangle and adjusted for each field point. @@ -477,25 +483,36 @@ def spherical_field(az_range,elev_range,outward_normals=False,field_radius=1.0): mesh : meshio object spherical field of points at specified azimuth and elevation angles, with meshed triangles """ - vista_pattern = pv.grid_from_sph_coords(az_range, (90-elev_range), field_radius).extract_surface() + vista_pattern = pv.grid_from_sph_coords( + az_range, (90 - elev_range), field_radius + ).extract_surface() if outward_normals: - vista_pattern.point_data['Normals']=vista_pattern.points/(np.linalg.norm(vista_pattern.points,axis=1).reshape(-1,1)) + vista_pattern.point_data["Normals"] = vista_pattern.points / ( + np.linalg.norm(vista_pattern.points, axis=1).reshape(-1, 1) + ) else: - vista_pattern.point_data['Normals'] = (vista_pattern.points / ( - np.linalg.norm(vista_pattern.points, axis=1).reshape(-1, 1)))*-1.0 + vista_pattern.point_data["Normals"] = ( + vista_pattern.points + / (np.linalg.norm(vista_pattern.points, axis=1).reshape(-1, 1)) + ) * -1.0 from ..utility.mesh_functions import pyvista_to_meshio - mesh=pyvista_to_meshio(vista_pattern.triangulate()) - from ..geometry.geometryfunctions import compute_areas,theta_phi_r - mesh=theta_phi_r(compute_areas(mesh)) + + mesh = pyvista_to_meshio(vista_pattern.triangulate()) + from ..geometry.geometryfunctions import compute_areas, theta_phi_r + + mesh = theta_phi_r(compute_areas(mesh)) return mesh -def linear_parabolic_aperture(diameter, focal_length, height, thickness, mesh_size, lip=False, lip_width=1e-3): + +def linear_parabolic_aperture( + diameter, focal_length, height, thickness, mesh_size, lip=False, lip_width=1e-3 +): # Define function for parabola equation (y^2 = 4*focal_length*x) def parabola(x): - return (1 / (4 * focal_length)) * x ** 2 + return (1 / (4 * focal_length)) * x**2 with pygmsh.occ.Geometry() as geom: geom.characteristic_length_max = mesh_size * 0.8 @@ -503,7 +520,9 @@ def parabola(x): point_num = 15 x_pos = np.linspace(-0.5 * diameter, 0.5 * diameter, point_num) z_pos = parabola(x_pos) - coords = np.array([x_pos.ravel(), np.zeros((point_num)), z_pos.ravel()]).transpose() + coords = np.array( + [x_pos.ravel(), np.zeros((point_num)), z_pos.ravel()] + ).transpose() points_list = [] for inc in range(point_num): points_list.append(geom.add_point(coords[inc, :].tolist())) @@ -525,9 +544,13 @@ def parabola(x): # volume_list.append(b2) if lip: - start_point = np.array([0.5 * diameter, 0.0, parabola(diameter * 0.5) - thickness]) + start_point = np.array( + [0.5 * diameter, 0.0, parabola(diameter * 0.5) - thickness] + ) r1 = geom.add_box(start_point, [lip_width, height, thickness]) - start_point2 = np.array([-0.5 * diameter - lip_width, 0.0, parabola(diameter * 0.5) - thickness]) + start_point2 = np.array( + [-0.5 * diameter - lip_width, 0.0, parabola(diameter * 0.5) - thickness] + ) r2 = geom.add_box(start_point2, [lip_width, height, thickness]) volume_list.append(r1) @@ -538,10 +561,11 @@ def parabola(x): mesh_temp = geom.generate_mesh(dim=2) for inc, cell in enumerate(mesh_temp.cells): - if cell.type == 'triangle': + if cell.type == "triangle": triangle_index = inc import meshio + triangle_cells = [("triangle", mesh_temp.cells[triangle_index].data)] mesh = meshio.Mesh(mesh_temp.points, triangle_cells) mesh = GF.compute_normals(mesh) @@ -549,6 +573,7 @@ def parabola(x): aperture_points = GF.cell_centroids(mesh) return mesh, aperture_points + def gridedReflectorPoints( majorsize, minorsize, thickness, grid_resolution, sides="all" ): @@ -713,6 +738,8 @@ def gridedReflectorPoints( ], point_data={"Normals": mesh_normals}, ) - mesh_points.point_data['Area']=np.ones((mesh_points.points.shape[0]))*(grid_resolution**2) + mesh_points.point_data["Area"] = np.ones((mesh_points.points.shape[0])) * ( + grid_resolution**2 + ) return mesh_points diff --git a/lyceanem/models/frequency_domain.py b/lyceanem/models/frequency_domain.py index 5ec9fe4..0bd6683 100644 --- a/lyceanem/models/frequency_domain.py +++ b/lyceanem/models/frequency_domain.py @@ -54,22 +54,23 @@ def aperture_projection( ) triangle_centroids = GF.cell_centroids(aperture) aperture = GF.compute_areas(aperture) - aperture=GF.compute_normals(aperture) - triangle_cell_index=GF.locate_cell_index(aperture) + aperture = GF.compute_normals(aperture) + triangle_cell_index = GF.locate_cell_index(aperture) triangle_normals = aperture.cell_data["Normals"][triangle_cell_index] # Ensure no clash with triangles in raycaster - triangle_centroids.points+=1e-6*triangle_normals + triangle_centroids.points += 1e-6 * triangle_normals import pyvista as pv - pl=pv.Plotter() - pl.add_mesh(pv.from_meshio(aperture),scalars="Area") - pl.add_mesh(pv.from_meshio(environment.solids[0]),color="red") - pl.add_mesh(pv.from_meshio(triangle_centroids),color="green") + + pl = pv.Plotter() + pl.add_mesh(pv.from_meshio(aperture), scalars="Area") + pl.add_mesh(pv.from_meshio(environment.solids[0]), color="red") + pl.add_mesh(pv.from_meshio(triangle_centroids), color="green") pl.show() visible_patterns, pcd = RF.visiblespace( triangle_centroids, triangle_normals, blocking_triangles, - vertex_area=aperture.point_data['Area'], + vertex_area=aperture.point_data["Area"], az_range=az_range, elev_range=elev_range, shell_range=farfield_distance, @@ -143,30 +144,39 @@ def calculate_farfield( # create sink points for the model from ..geometry.targets import spherical_field + sink_coords = spherical_field(az_range, el_range, farfield_distance) - Ex, Ey, Ez=calculate_scattering(aperture_coords, - sink_coords, - antenna_solid, - desired_E_axis, - scatter_points=scatter_points, - wavelength=wavelength, - scattering=scattering, - elements=elements, - los=los, - mesh_resolution=mesh_resolution, - project_vectors=project_vectors, - antenna_axes=antenna_axes, - multiE=False, - alpha=alpha, - beta=beta) + Ex, Ey, Ez = calculate_scattering( + aperture_coords, + sink_coords, + antenna_solid, + desired_E_axis, + scatter_points=scatter_points, + wavelength=wavelength, + scattering=scattering, + elements=elements, + los=los, + mesh_resolution=mesh_resolution, + project_vectors=project_vectors, + antenna_axes=antenna_axes, + multiE=False, + alpha=alpha, + beta=beta, + ) # convert to etheta,ephi etheta = ( - Ex * np.cos(sink_coords.point_data["phi_(Radians)"]) * np.cos(sink_coords.point_data["theta_(Radians)"]) - + Ey * np.sin(sink_coords.point_data["phi_(Radians)"]) * np.cos(sink_coords.point_data["theta_(Radians)"]) + Ex + * np.cos(sink_coords.point_data["phi_(Radians)"]) + * np.cos(sink_coords.point_data["theta_(Radians)"]) + + Ey + * np.sin(sink_coords.point_data["phi_(Radians)"]) + * np.cos(sink_coords.point_data["theta_(Radians)"]) - Ez * np.sin(sink_coords.point_data["theta_(Radians)"]) ) - ephi = -Ex * np.sin(sink_coords.point_data["phi_(Radians)"]) + Ey * np.cos(sink_coords.point_data["phi_(Radians)"]) + ephi = -Ex * np.sin(sink_coords.point_data["phi_(Radians)"]) + Ey * np.cos( + sink_coords.point_data["phi_(Radians)"] + ) return etheta, ephi diff --git a/lyceanem/raycasting/rayfunctions.py b/lyceanem/raycasting/rayfunctions.py index 4096948..f57cb8d 100644 --- a/lyceanem/raycasting/rayfunctions.py +++ b/lyceanem/raycasting/rayfunctions.py @@ -61,7 +61,7 @@ def hit(ray, triangle): # if A is less than zero, then the ray is coming from behind the triangle # don't cull backface triangles - #if A < 0: + # if A < 0: # return False, math.inf # calculate distance from vertice 0 to ray origin tvecx = ray.ox - triangle.v0x # s @@ -386,7 +386,7 @@ def visiblespace( ) # angles[np.isnan(angles)]=0 - vertex_area[np.isnan(vertex_area)]=0 + vertex_area[np.isnan(vertex_area)] = 0 # visible_patterns=quickpatterncreator(az_range,elev_range,source_coords,angles,vertex_area,hit_index) if len(vertex_area) == 1: if vertex_area == 0: @@ -1842,9 +1842,9 @@ def charge_rays_environment1Dv2(sources, sinks, environment_points, point_indexi ray payload to be sent to GPU """ - #print("sources shape", sources.shape) - #print("sinks shape", sinks.shape) - #print("environment_points shape", environment_points.shape) + # print("sources shape", sources.shape) + # print("sinks shape", sinks.shape) + # print("environment_points shape", environment_points.shape) unified_model = np.append( np.append(sources, sinks, axis=0), environment_points, axis=0 @@ -2140,10 +2140,12 @@ def chunkingRaycaster1Dv3( ): start = timer() cuda.current_context().memory_manager.deallocations.clear() - + free_mem, total_mem = cuda.current_context().get_memory_info() max_mem = np.ceil(free_mem * 0.5).astype(np.int64) - ray_limit = np.floor(((max_mem - environment_local.nbytes) / base_types.ray_t.size)).astype(np.int64) + ray_limit = np.floor( + ((max_mem - environment_local.nbytes) / base_types.ray_t.size) + ).astype(np.int64) sink_index = np.arange( sources.shape[0] + 1, sources.shape[0] + 1 + sinks.shape[0] ).reshape(sinks.shape[0], 1) @@ -2770,7 +2772,7 @@ def workchunkingv2( ) if io_indexing.shape[0] >= ray_limit: # need to split the array and process seperatly - + sub_io = np.array_split( io_indexing, np.ceil(io_indexing.shape[0] / ray_limit).astype(int) ) @@ -2926,8 +2928,10 @@ def workchunkingv2( sources, sinks, scattering_points, filtered_index, environment, False ) if filtered_index2.shape[0] * sinks.shape[0] > 2e8: - temp_chunks = np.ceil((filtered_index2.shape[0] * sinks.shape[0]) / 2e8).astype(int) - + temp_chunks = np.ceil( + (filtered_index2.shape[0] * sinks.shape[0]) / 2e8 + ).astype(int) + temp_chunking = np.linspace( 0, filtered_index2.shape[0], temp_chunks, dtype=np.int32 ) diff --git a/lyceanem/tests/reflectordata.py b/lyceanem/tests/reflectordata.py index 7ebe156..98bd6f6 100644 --- a/lyceanem/tests/reflectordata.py +++ b/lyceanem/tests/reflectordata.py @@ -54,18 +54,18 @@ def exampleUAV(frequency): array, np.array([0.25, 0, 0]) + np.array([-0.18, 0, 0.0125]) ) - #def structure_cells(array): - ## add collumn of 3s to beggining of each row + # def structure_cells(array): + ## add collumn of 3s to beggining of each row # array = np.append( # np.ones((array.shape[0], 1), dtype=np.int32) * 3, array, axis=1 # ) # return array - #pyvista_array = pv.PolyData(array.points, structure_cells(array.cells[0].data)) - #pyvista_body = pv.PolyData(body.points, structure_cells(body.cells[0].data)) - #pyvista_array.compute_normals(inplace=True) - #pyvista_body.compute_normals(inplace=True) - + # pyvista_array = pv.PolyData(array.points, structure_cells(array.cells[0].data)) + # pyvista_body = pv.PolyData(body.points, structure_cells(body.cells[0].data)) + # pyvista_array.compute_normals(inplace=True) + # pyvista_body.compute_normals(inplace=True) + array = GF.compute_normals(array) body = GF.compute_normals(body) diff --git a/lyceanem/utility/math_functions.py b/lyceanem/utility/math_functions.py index 0f36e6f..4b9ad9d 100644 --- a/lyceanem/utility/math_functions.py +++ b/lyceanem/utility/math_functions.py @@ -105,11 +105,13 @@ def discrete_transmit_power( power_at_point * power_normalisation ) / element_area.reshape(-1, 1) # calculate amplitude (V/m) - transmit_amplitude = ((transmit_power_density * impedance) ** 0.5)*element_area.reshape(-1,1) - - polarized_amplitudes=(weights/np.linalg.norm(weights, axis=1).reshape( - -1, 1 - ))*transmit_amplitude + transmit_amplitude = ( + (transmit_power_density * impedance) ** 0.5 + ) * element_area.reshape(-1, 1) + + polarized_amplitudes = ( + weights / np.linalg.norm(weights, axis=1).reshape(-1, 1) + ) * transmit_amplitude # transmit_excitation=transmit_amplitude_density.reshape(-1,1)*element_area.reshape(-1,1) return polarized_amplitudes diff --git a/lyceanem/utility/mesh_functions.py b/lyceanem/utility/mesh_functions.py index 13a1cee..7c7cae5 100644 --- a/lyceanem/utility/mesh_functions.py +++ b/lyceanem/utility/mesh_functions.py @@ -1,15 +1,16 @@ import pyvista as pv import meshio + def pyvista_to_meshio(polydata_object): """ Convert a pyvista object to a meshio object """ - #extract only the triangles - if type(polydata_object )==pv.core.pointset.UnstructuredGrid: - cells =id_cells(polydata_object.cells) + # extract only the triangles + if type(polydata_object) == pv.core.pointset.UnstructuredGrid: + cells = id_cells(polydata_object.cells) else: - cells =id_cells(polydata_object.faces) + cells = id_cells(polydata_object.faces) meshio_object = meshio.Mesh( points=polydata_object.points, cells=cells, @@ -17,28 +18,23 @@ def pyvista_to_meshio(polydata_object): ) return meshio_object + def id_cells(faces): import copy - #temp_faces=copy.deepcopy(faces) + + # temp_faces=copy.deepcopy(faces) import numpy as np - cell_types ={1 :"vertex", - 2 :"line", - 3 :"triangle", - 4 :"quad" - } - cells ={"vertex" :[], - "line" :[], - "triangle" :[], - "quad" :[] - } - - moving_index=0 - while moving_index<=faces.shape[0]-1: - face_num=faces[moving_index] - temp_array=faces[moving_index+1:moving_index+face_num+1] + + cell_types = {1: "vertex", 2: "line", 3: "triangle", 4: "quad"} + cells = {"vertex": [], "line": [], "triangle": [], "quad": []} + + moving_index = 0 + while moving_index <= faces.shape[0] - 1: + face_num = faces[moving_index] + temp_array = faces[moving_index + 1 : moving_index + face_num + 1] cells[cell_types[face_num]].append(temp_array.tolist()) - moving_index+=face_num+1 - + moving_index += face_num + 1 + # while (temp_faces.shape[0]>=1): # trim_num =temp_faces[0] # if trim_num>3: @@ -47,9 +43,9 @@ def id_cells(faces): # cells[cell_types[trim_num]].append(temp_array.tolist()) # temp_faces =np.delete(temp_faces ,np.arange(0 ,trim_num +1)) - meshio_cells =[] + meshio_cells = [] for key in cells: - if len(cells[key] ) >0: - meshio_cells.append((key ,cells[key])) + if len(cells[key]) > 0: + meshio_cells.append((key, cells[key])) - return meshio_cells \ No newline at end of file + return meshio_cells