From 2d9f1bbe0146faf29575a5510e75b3e61052dbf7 Mon Sep 17 00:00:00 2001 From: LyceanEM <60020395+LyceanEM@users.noreply.github.com> Date: Mon, 20 May 2024 12:54:48 +0100 Subject: [PATCH 1/4] Shifting the different propagation kernals to the new algorithm. --- .../03_frequency_domain_channel_modelling.py | 10 ++-- .../04_time_domain_channel_modelling.py | 2 +- lyceanem/electromagnetics/empropagation.py | 46 ++++++++++++++++--- 3 files changed, 46 insertions(+), 12 deletions(-) diff --git a/docs/examples/03_frequency_domain_channel_modelling.py b/docs/examples/03_frequency_domain_channel_modelling.py index a8fb6f5..3c917d3 100644 --- a/docs/examples/03_frequency_domain_channel_modelling.py +++ b/docs/examples/03_frequency_domain_channel_modelling.py @@ -21,7 +21,7 @@ # Frequency and Mesh Resolution # ------------------------------ # -freq = np.asarray(24.0e9) +freq = np.asarray(16.0e9) wavelength = 3e8 / freq mesh_resolution = 0.5 * wavelength @@ -179,7 +179,7 @@ def structure_cells(array): -angle_values = np.linspace(0, 90, 91) +angle_values = np.linspace(0, 90, 181) angle_increment = np.diff(angle_values)[0] responsex = np.zeros((len(angle_values)), dtype="complex") responsey = np.zeros((len(angle_values)), dtype="complex") @@ -223,9 +223,9 @@ def structure_cells(array): scattering=1, project_vectors=False ) - responsex[angle_inc] = Ex - responsey[angle_inc] = Ey - responsez[angle_inc] = Ez + responsex[angle_inc] = np.sum(Ex) + responsey[angle_inc] = np.sum(Ey) + responsez[angle_inc] = np.sum(Ez) diff --git a/docs/examples/04_time_domain_channel_modelling.py b/docs/examples/04_time_domain_channel_modelling.py index 6ebde72..82ee9ba 100644 --- a/docs/examples/04_time_domain_channel_modelling.py +++ b/docs/examples/04_time_domain_channel_modelling.py @@ -30,7 +30,7 @@ noise_power = -80 # dbw mean_noise = 0 -model_freq = 16e9 +model_freq = 24e9 wavelength = 3e8 / model_freq # %% diff --git a/lyceanem/electromagnetics/empropagation.py b/lyceanem/electromagnetics/empropagation.py index 70118fe..d37acdb 100644 --- a/lyceanem/electromagnetics/empropagation.py +++ b/lyceanem/electromagnetics/empropagation.py @@ -1266,10 +1266,27 @@ def freqdomainisokernal( # print(cu_ray_num,source_sink_index[cu_ray_num,0],source_sink_index[cu_ray_num,1]) wave_vector = (2.0 * cmath.pi) / wavelength # scatter_coefficient=(1/(4*cmath.pi))**(complex(scatter_index)) - loss = cmath.exp(-lengths * wave_vector * 1j) * ( - wavelength / (4 * cmath.pi * lengths) + alpha = 0.0 + beta = (2.0 * cmath.pi) / wavelength[0] + loss = lossy_propagation( + calc_sep( + point_information[network_index[cu_ray_num, 0] - 1], + point_information[network_index[cu_ray_num, 1] - 1], + float(0) + ), + alpha, + beta ) - ray_component[0] = loss + for i in range(1,network_index.shape[1] - 1): + if network_index[cu_ray_num, i + 1] != 0: + + loss *= lossy_propagation(calc_sep(point_information[network_index[cu_ray_num, i] - 1], + point_information[network_index[cu_ray_num, i + 1] - 1], + float(0)), + alpha, + beta + ) + ray_component[0] *= loss # ray_component[1]=loss # ray_component[2]=loss # print(ray_component[0].real,ray_component[1].real,ray_component[2].real) @@ -1419,9 +1436,26 @@ def timedomainkernal( i = i + 1 # print(cu_ray_num,source_sink_index[cu_ray_num,0],source_sink_index[cu_ray_num,1]) - wave_vector = (2.0 * cmath.pi) / wavelength[0] - # scatter_coefficient=(1/(4*cmath.pi))**(complex(scatter_index)) - loss = wavelength[0] / (4 * cmath.pi * lengths) + alpha = 0.0 + beta = (2.0 * cmath.pi) / wavelength[0] + loss = lossy_propagation( + calc_sep( + point_information[full_index[cu_ray_num, 0] - 1], + point_information[full_index[cu_ray_num, 1] - 1], + float(0) + ), + alpha, + beta + ) + for i in range(1,full_index.shape[1] - 1): + if full_index[cu_ray_num, i + 1] != 0: + + loss *= lossy_propagation(calc_sep(point_information[full_index[cu_ray_num, i] - 1], + point_information[full_index[cu_ray_num, i + 1] - 1], + float(0)), + alpha, + beta + ) ray_component[0] *= loss ray_component[1] *= loss From 6db27ad8b091e80fc4787391f242c5290bf44a7f Mon Sep 17 00:00:00 2001 From: LyceanEM <60020395+LyceanEM@users.noreply.github.com> Date: Mon, 20 May 2024 22:27:15 +0100 Subject: [PATCH 2/4] Rolling parabolic reflector function using pygmsh and meshio into targets.py --- lyceanem/geometry/geometryfunctions.py | 38 +++++++++ lyceanem/geometry/targets.py | 111 ++++++++++++++++++++++++- 2 files changed, 147 insertions(+), 2 deletions(-) diff --git a/lyceanem/geometry/geometryfunctions.py b/lyceanem/geometry/geometryfunctions.py index 31597c7..8f9213f 100644 --- a/lyceanem/geometry/geometryfunctions.py +++ b/lyceanem/geometry/geometryfunctions.py @@ -92,6 +92,44 @@ def mesh_transform(mesh, transform_matrix, rotate_only): return return_mesh +def compute_normals(mesh): + """ + Computes the Cell Normals for meshio mesh objects, this does not currently calculate the point normals, but this will be done soon. + + Parameters + ---------- + mesh + + Returns + ------- + + """ + cell_normal_list = [] + for inc, cell in enumerate(mesh.cells): + print(cell.type, cell.data.shape[0]) + if cell.type == 'vertex': + vertex_normals = np.zeros((cell.data.shape[0], 3)) + cell_normal_list.append(vertex_normals) + if cell.type == 'line': + line_normals = np.zeros((cell.data.shape[0], 3)) + cell_normal_list.append(line_normals) + if cell.type == 'triangle': + # print(inc) + edge1 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 1], :] + edge2 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 2], :] + tri_cell_normals = np.cross(edge1, edge2) + tri_cell_normals *= (1 / np.linalg.norm(tri_cell_normals, axis=1)).reshape(-1, 1) + cell_normal_list.append(tri_cell_normals) + if cell.type == 'tetra': + # print(inc) + edge1 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 1], :] + edge2 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 2], :] + tetra_cell_normals = np.cross(edge1, edge2) + tetra_cell_normals *= (1 / np.linalg.norm(tetra_cell_normals, axis=1)).reshape(-1, 1) + cell_normal_list.append(tetra_cell_normals) + + mesh.cell_data['Normals'] = cell_normal_list + return mesh def mesh_conversion(conversion_object): """ diff --git a/lyceanem/geometry/targets.py b/lyceanem/geometry/targets.py index 20345c0..ca9e5ad 100755 --- a/lyceanem/geometry/targets.py +++ b/lyceanem/geometry/targets.py @@ -5,8 +5,7 @@ import numpy as np import meshio import pyvista as pv -import scipy.stats -import solid as sd +import pygmsh from importlib_resources import files from scipy.spatial.transform import Rotation as R @@ -305,6 +304,114 @@ def meshedHorn( return structure, mesh_points +def parabolic_aperture(diameter, focal_length, thickness, mesh_size, sides='front'): + # Define function for parabola equation (y^2 = 4*focal_length*x) + def parabola(x): + return (1 / (4 * focal_length)) * x ** 2 + + with pygmsh.occ.Geometry() as geom: + geom.characteristic_length_max = mesh_size + # Define points + cp1 = geom.add_point([0, 0, 0]) # Center point + cp2 = geom.add_point([diameter * 0.5 * (1 / 6), 0, parabola(diameter * 0.5 * (1 / 6))]) + cp3 = geom.add_point([diameter * 0.5 * (2 / 6), 0, parabola(diameter * 0.5 * (2 / 6))]) + cp4 = geom.add_point([diameter * 0.5 * (3 / 6), 0, parabola(diameter * 0.5 * (3 / 6))]) + cp5 = geom.add_point([diameter * 0.5 * (4 / 6), 0, parabola(diameter * 0.5 * (4 / 6))]) + cp6 = geom.add_point([diameter * 0.5 * (5 / 6), 0, parabola(diameter * 0.5 * (5 / 6))]) + cp7 = geom.add_point([diameter * 0.5 * (6 / 6), 0, parabola(diameter * 0.5 * (6 / 6))]) + + # Define top line based on points + line = geom.add_bspline([cp1, cp2, cp3, cp4, cp5, cp6, cp7]) + + _, surface, _ = geom.extrude(line, translation_axis=[0.0, 0.0, -thickness]) + + # Revolve line to create revolution surface + volume_list = [] + _, b, _ = geom.revolve(surface, rotation_axis=[0.0, 0.0, 1.0], point_on_axis=[0.0, 0.0, 0.0], + angle=0.25 * np.pi) + volume_list.append(b) + for inc in range(7): + + geom.rotate(surface, point=[0.0, 0.0, 0.0], angle=(1 / 4) * np.pi, axis=[0.0, 0.0, 1.0]) + _, b2, _ = geom.revolve(surface, rotation_axis=[0.0, 0.0, 1.0], point_on_axis=[0.0, 0.0, 0.0], + angle=0.25 * np.pi) + volume_list.append(b2) + + full_reflector = geom.boolean_union(volume_list) + + mesh_temp = geom.generate_mesh(dim=2) + for inc, cell in enumerate(mesh_temp.cells): + 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) + + x_space = np.linspace( + mesh_size, + (diameter / 2), + int(np.max(np.asarray([2, np.ceil((diameter * 0.5) / (mesh_size))]))), + ) + z_space = (1 / (4 * focal_length)) * x_space ** 2 + c_space = np.ceil((2 * np.pi * x_space) / mesh_size).astype(int) + normal_gradiant_vector = np.array( + [ + np.ones((len(x_space))), + np.zeros((len(x_space))), + -1 / (1 / (2 * focal_length) * x_space), + ] + ) + source = np.array([x_space, np.zeros((len(x_space))), z_space]).transpose() + target = (normal_gradiant_vector * -x_space).transpose() + base_directions = np.zeros((x_space.shape[0], 3), dtype=np.float32) + norm_length = np.zeros((x_space.shape[0], 1), dtype=np.float32) + base_directions, norm_length = math_functions.calc_dv_norm( + source, target, base_directions, norm_length + ) + source_coords = np.empty(((1, 3)), dtype=np.float32) + source_coords[0, :] = 0 + source_normals = np.empty(((1, 3)), dtype=np.float32) + source_normals[0, :] = 0 + source_normals[0, 2] = 1 + for r_index in range(x_space.shape[0]): + source_coords = np.append( + source_coords, + np.array( + [ + x_space[r_index] + * np.cos(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]), + x_space[r_index] + * np.sin(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]), + z_space[r_index] * np.ones(c_space[r_index] - 1), + ] + ).transpose(), + axis=0, + ) + source_normals = np.append( + source_normals, + np.array( + [ + base_directions[r_index, 0] + * np.cos(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]), + base_directions[r_index, 0] + * np.sin(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]), + base_directions[r_index, 2] * np.ones(c_space[r_index] - 1), + ] + ).transpose(), + axis=0, + ) + + mesh_vertices = source_coords + np.array([0, 0, 1e-6]) + mesh_normals = source_normals + aperture_points = meshio.Mesh(points=mesh_vertices, + cells=[("vertex", np.array([[i, ] for i in range(len(mesh_vertices))]))], + point_data={'Normals': mesh_normals}) + + return mesh, aperture_points + + + From 91d1b8e65300fbe67c1b3ffd5724bf281764b698 Mon Sep 17 00:00:00 2001 From: LyceanEM Date: Thu, 23 May 2024 16:29:28 +0100 Subject: [PATCH 3/4] Changing "normals" to "Normals" as this is vtk convention. Also adapted convert triangles to find the triangle cells in a multicell mesh. --- docs/examples/01_aperture_projection.py | 2 +- .../examples/02_coherently_polarised_array.py | 2 +- .../07_aperture_farfield_polarisation.py | 39 +------------------ lyceanem/tests/reflectordata.py | 10 ++--- 4 files changed, 9 insertions(+), 44 deletions(-) diff --git a/docs/examples/01_aperture_projection.py b/docs/examples/01_aperture_projection.py index 9799548..1a2504b 100644 --- a/docs/examples/01_aperture_projection.py +++ b/docs/examples/01_aperture_projection.py @@ -63,7 +63,7 @@ surface_array = copy.deepcopy(array) surface_array.cells[0].data = np.asarray(array.cells[0].data)[: (array.cells[0].data).shape[0] // 2, :] -surface_array.cell_data["normals"] = np.array(array.cell_data["normals"])[: (array.cells[0].data).shape[0] // 2] +surface_array.cell_data["Normals"] = np.array(array.cell_data["Normals"])[: (array.cells[0].data).shape[0] // 2] # %% # Structures diff --git a/docs/examples/02_coherently_polarised_array.py b/docs/examples/02_coherently_polarised_array.py index 19636d7..148b933 100644 --- a/docs/examples/02_coherently_polarised_array.py +++ b/docs/examples/02_coherently_polarised_array.py @@ -48,7 +48,7 @@ surface_array = copy.deepcopy(array) surface_array.cells[0].data = np.asarray(array.cells[0].data)[: (array.cells[0].data).shape[0] // 2, :] -surface_array.cell_data["normals"] = np.array(array.cell_data["normals"])[: (array.cells[0].data).shape[0] // 2] +surface_array.cell_data["Normals"] = np.array(array.cell_data["Normals"])[: (array.cells[0].data).shape[0] // 2] from lyceanem.base_classes import structures diff --git a/docs/examples/07_aperture_farfield_polarisation.py b/docs/examples/07_aperture_farfield_polarisation.py index 75ed619..7fdbaca 100644 --- a/docs/examples/07_aperture_farfield_polarisation.py +++ b/docs/examples/07_aperture_farfield_polarisation.py @@ -37,43 +37,8 @@ from lyceanem.geometry.targets import meshedHorn structure,array_points=meshedHorn(3*wavelength, 1*wavelength, 4*wavelength, 1*wavelength,np.radians(10),wavelength*0.5) -def planar_aperture(u_size,v_size,wall_thickness=2e-3,depth=2e-3,mesh_size=1.0): - with pygmsh.occ.Geometry() as geom: - geom.characteristic_length_max=mesh_size - u_points=np.linspace(-u_size/2,u_size/2,np.ceil(u_size/mesh_size).astype(int)) - v_points=np.linspace(-v_size/2,v_size/2,np.ceil(v_size/mesh_size).astype(int)) - n_points=np.linspace(0,0,1) - u_mesh,v_mesh,n_mesh=np.meshgrid(u_points,v_points,n_points,indexing='ij') - aperture_points=np.array([u_mesh.ravel(),v_mesh.ravel(),n_mesh.ravel()]).transpose() - source_normals = np.empty(((aperture_points.shape[0], 3)), dtype=np.float32) - source_normals[:] = 0.0 - source_normals[:, 2] = 1.0 - ground_plane=geom.add_box((-(u_size+wall_thickness*2)/2.0,-(v_size+wall_thickness*2)/2.0,-depth), extents=(u_size+wall_thickness*2,v_size+wall_thickness*2,depth),mesh_size=mesh_size) - point_list=[] - for inc in range(aperture_points.shape[0]): - point_list.append(geom.add_point(aperture_points[inc,:])) - geom.in_surface(point_list[-1],ground_plane) - - #aperture=geom.add_rectangle((-u_size/2.0,-v_size/2.0,0.0), u_size, v_size,mesh_size=mesh_size) - #ground_plane=geom.add_box((-(u_size+wall_thickness*2)/2.0,-(v_size+wall_thickness*2)/2.0,-depth), extents=(u_size+wall_thickness*2,v_size+wall_thickness*2,depth),mesh_size=mesh_size) - - - mesh = geom.generate_mesh() - end_point=mesh.points.shape[0] - mesh.point_data['Normals']=np.zeros((mesh.points.shape[0],3)) - mesh.point_data['Normals'][mesh.points[:,2]>-depth*0.5,2]=1.0 - mesh.point_data['Normals'][mesh.points[:,2]>-depth*0.5,2]=-1.0 - mesh.points=np.append(mesh.points,aperture_points,axis=0) - mesh.point_data['Normals']=np.append(mesh.point_data['Normals'],source_normals,axis=0) - mesh.cell_sets['Aperture']=[None,None,None,np.linspace(end_point,aperture_points.shape[0]+end_point-1,aperture_points.shape[0]).astype(int)] - mesh.cells[3].data=np.append(mesh.cells[3].data,np.linspace(end_point,aperture_points.shape[0]+end_point-1,aperture_points.shape[0]).astype(int).reshape(-1,1),axis=0) - return mesh - -mesh2=planar_aperture(3*wavelength,1*wavelength,mesh_size=wavelength*0.5) - - - -horn_antenna=antenna_structures(structures(solids=[mesh2]), points(points=[array_points])) + +horn_antenna=antenna_structures(structures(solids=[structure]), points(points=[array_points])) from lyceanem.models.frequency_domain import calculate_farfield diff --git a/lyceanem/tests/reflectordata.py b/lyceanem/tests/reflectordata.py index 628c377..fd269de 100644 --- a/lyceanem/tests/reflectordata.py +++ b/lyceanem/tests/reflectordata.py @@ -61,10 +61,10 @@ def structure_cells(array): pyvista_array.compute_normals(inplace=True) pyvista_body.compute_normals(inplace=True) - array.point_data['normals'] = pyvista_array.point_normals - body.point_data['normals'] = pyvista_body.point_normals - array.cell_data['normals'] = pyvista_array.cell_normals - body.cell_data['normals'] = pyvista_body.cell_normals + array.point_data['Normals'] = pyvista_array.point_normals + body.point_data['Normals'] = pyvista_body.point_normals + array.cell_data['Normals'] = pyvista_array.cell_normals + body.cell_data['Normals'] = pyvista_body.cell_normals wavelength = 3e8 / frequency @@ -195,7 +195,7 @@ def structure_cells(array): axis=0, ) + np.array([0.1025, 0, -0.025]) source_pcd = RF.points2pointcloud(total_array) - source_pcd.point_data['normals'] = total_array_normals + source_pcd.point_data['Normals'] = total_array_normals source_pcd = GF.translate_mesh(source_pcd, np.array([-0.18, 0, 0.0125])) From e7bed534abaa29ddb06c7431154827e4ab61a4f1 Mon Sep 17 00:00:00 2001 From: LyceanEM <60020395+LyceanEM@users.noreply.github.com> Date: Fri, 24 May 2024 10:39:37 +0100 Subject: [PATCH 4/4] Corrected file addresses for test assets. --- lyceanem/tests/test_base_classes.py | 5 +++-- lyceanem/tests/test_geometryfunctions.py | 22 ++++++++++++---------- lyceanem/tests/test_ray_functions.py | 7 ++++--- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/lyceanem/tests/test_base_classes.py b/lyceanem/tests/test_base_classes.py index 104932a..c8c104f 100644 --- a/lyceanem/tests/test_base_classes.py +++ b/lyceanem/tests/test_base_classes.py @@ -33,14 +33,15 @@ def cube(): [1, 4, 5]] ## put into meshio mesh cube = meshio.Mesh(points=cube_points, cells=[("triangle", cube_cells)]) - + from ..geometry.geometryfunctions import compute_normals + cube=compute_normals(cube) return cube def point(): # a single point on the +x center of the cube with consistent normal vector - pc_mesh = meshio.Mesh(points=np.array([[0.5, 0, 0]]), cells=[], point_data={"normals": np.array([[1.0, 0, 0]])}) + pc_mesh = meshio.Mesh(points=np.array([[0.5, 0, 0]]), cells=[], point_data={"Normals": np.array([[1.0, 0, 0]])}) return pc_mesh diff --git a/lyceanem/tests/test_geometryfunctions.py b/lyceanem/tests/test_geometryfunctions.py index bb2fa01..596f935 100644 --- a/lyceanem/tests/test_geometryfunctions.py +++ b/lyceanem/tests/test_geometryfunctions.py @@ -1,17 +1,19 @@ import numpy as np - +import lyceanem import pytest from numpy.testing import assert_allclose from scipy.spatial.transform import Rotation as R import meshio - +from importlib_resources import files from ..geometry import geometryfunctions as GF +import lyceanem.tests.data @pytest.fixture def standard_cube(): - cube = meshio.read('data/cube.ply') + cube_location= files(lyceanem.tests.data).joinpath("cube.ply") + cube = meshio.read(cube_location) return cube @@ -77,21 +79,21 @@ def are_meshes_equal(mesh1, mesh2, rtol=1e-5): def test_tri_areas_2(): - refrence = np.load('data/areas.npy') - input = meshio.read('data/receive_horn.ply') + refrence = np.load(files(lyceanem.tests.data).joinpath("areas.npy")) + input = meshio.read(files(lyceanem.tests.data).joinpath("receive_horn.ply")) result = GF.tri_areas(input) assert np.allclose(result, refrence) def test_tri_centroids(): - refrence = np.load('data/centroid_numpy.npy') - refrenece_mesh = meshio.read('data/centroid_mesh.ply') - input = meshio.read('data/receive_horn.ply') + refrence = np.load(files(lyceanem.tests.data).joinpath('centroid_numpy.npy')) + refrenece_mesh = meshio.read(files(lyceanem.tests.data).joinpath('centroid_mesh.ply')) + input = meshio.read(files(lyceanem.tests.data).joinpath('receive_horn.ply')) result_numpy, result_mesh = GF.tri_centroids(input) assert np.allclose(result_numpy, refrence) assert are_meshes_equal(result_mesh, refrenece_mesh) def test_mesh_rotate(): - refrence = meshio.read('data/rotated_recieve.ply') - input = meshio.read('data/receive_horn.ply') + refrence = meshio.read(files(lyceanem.tests.data).joinpath('rotated_recieve.ply')) + input = meshio.read(files(lyceanem.tests.data).joinpath('receive_horn.ply')) rotation_vector1 = np.radians(np.asarray([90.0, 0.0, 0.0])) centre = np.array([1.0, 1.0, 1.0]) diff --git a/lyceanem/tests/test_ray_functions.py b/lyceanem/tests/test_ray_functions.py index 8cb4322..fc67be5 100644 --- a/lyceanem/tests/test_ray_functions.py +++ b/lyceanem/tests/test_ray_functions.py @@ -2,12 +2,13 @@ import numpy as np import pytest import lyceanem.base_types as base_types - +from importlib_resources import files +import lyceanem.tests.data import meshio def test_convertTriangles(): - reference = np.load("data/triangle_type_array_ref.npy") - mesh = meshio.read("data/receive_horn.ply") + reference = np.load(files(lyceanem.tests.data).joinpath("triangle_type_array_ref.npy")) + mesh = meshio.read(files(lyceanem.tests.data).joinpath("receive_horn.ply")) triangles = RF.convertTriangles(mesh) for i in range(triangles.size): assert np.allclose(triangles[i]["v0x"], reference[i]["v0x"])