diff --git a/docs/examples/07_aperture_farfield_polarisation.py b/docs/examples/07_aperture_farfield_polarisation.py index cfcfb6e..75ed619 100644 --- a/docs/examples/07_aperture_farfield_polarisation.py +++ b/docs/examples/07_aperture_farfield_polarisation.py @@ -11,6 +11,7 @@ """ import numpy as np +import pygmsh # %% # Setting Farfield Resolution and Wavelength @@ -34,9 +35,45 @@ from lyceanem.base_classes import points,structures,antenna_structures from lyceanem.geometry.targets import meshedHorn -structure,array_points=meshedHorn(3*wavelength, 1*wavelength, 4*wavelength, 0.05*wavelength,np.radians(10),wavelength*0.5) - -horn_antenna=antenna_structures(structures(solids=[structure]), points(points=[array_points])) +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])) from lyceanem.models.frequency_domain import calculate_farfield diff --git a/lyceanem/base_classes.py b/lyceanem/base_classes.py index 9fec916..dcb8b8c 100644 --- a/lyceanem/base_classes.py +++ b/lyceanem/base_classes.py @@ -390,7 +390,7 @@ def excitation_function( else: aperture_points = self.export_all_points(point_index=point_index) aperture_weights = EM.calculate_conformalVectors( - desired_e_vector, aperture_points.point_data['normals'], self.pose[:3, :3] + desired_e_vector, aperture_points.point_data['Normals'], self.pose[:3, :3] ) if phase_shift == "wavefront": source_points = np.asarray(aperture_points.points) diff --git a/lyceanem/geometry/geometryfunctions.py b/lyceanem/geometry/geometryfunctions.py index fbad032..31597c7 100644 --- a/lyceanem/geometry/geometryfunctions.py +++ b/lyceanem/geometry/geometryfunctions.py @@ -82,12 +82,12 @@ def mesh_transform(mesh, transform_matrix, rotate_only): if rotate_only: for i in range(mesh.points.shape[0]): return_mesh.points[i] = np.dot(transform_matrix, np.append(mesh.points[i], 0))[:3] - return_mesh.point_data['normals'][i] = np.dot(transform_matrix, np.append(mesh.point_data['normals'][i], 0))[:3] + return_mesh.point_data['Normals'][i] = np.dot(transform_matrix, np.append(mesh.point_data['Normals'][i], 0))[:3] else: for i in range(mesh.points.shape[0]): return_mesh.points[i] = np.dot(transform_matrix, np.append(mesh.points[i], 1))[:3] - return_mesh.point_data['normals'][i]= np.dot(transform_matrix, np.append(mesh.point_data['normals'][i], 0))[:3] + return_mesh.point_data['Normals'][i]= np.dot(transform_matrix, np.append(mesh.point_data['Normals'][i], 0))[:3] return return_mesh diff --git a/lyceanem/geometry/targets.py b/lyceanem/geometry/targets.py index f5279e3..435ca97 100755 --- a/lyceanem/geometry/targets.py +++ b/lyceanem/geometry/targets.py @@ -125,8 +125,8 @@ def rectReflector(majorsize, minorsize, thickness): - mesh.point_data["normals"] = pv_mesh.point_normals - mesh.cell_data["normals"] = pv_mesh.cell_normals + mesh.point_data["Normals"] = pv_mesh.point_normals + mesh.cell_data["Normals"] = pv_mesh.cell_normals 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 @@ -199,8 +199,8 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle): pv_mesh = pv.PolyData( mesh_vertices, faces = triangle_list) pv_mesh.compute_normals(inplace=True,consistent_normals=False) - mesh.point_data["normals"] = np.asarray(pv_mesh.point_normals) - mesh.cell_data["normals"] = np.asarray(pv_mesh.cell_normals) + mesh.point_data["Normals"] = np.asarray(pv_mesh.point_normals) + mesh.cell_data["Normals"] = np.asarray(pv_mesh.cell_normals) red = np.zeros((8, 1), dtype=np.float32) green = np.ones((8, 1), dtype=np.float32) * 0.259 @@ -476,6 +476,6 @@ def gridedReflectorPoints( mesh_vertices = source_coords mesh_normals = source_normals - mesh_points = meshio.Mesh(points=mesh_vertices, cells=[], point_data={"normals": mesh_normals}) + mesh_points = meshio.Mesh(points=mesh_vertices, cells=[], point_data={"Normals": mesh_normals}) return mesh_points diff --git a/lyceanem/models/frequency_domain.py b/lyceanem/models/frequency_domain.py index a17713e..e50d5c0 100644 --- a/lyceanem/models/frequency_domain.py +++ b/lyceanem/models/frequency_domain.py @@ -54,7 +54,7 @@ def aperture_projection( ) triangle_centroids, _ = GF.tri_centroids(aperture) triangle_areas = GF.tri_areas(aperture) - triangle_normals = np.asarray(aperture.cell_data["normals"]) + triangle_normals = np.asarray(aperture.cell_data["Normals"]) visible_patterns, pcd = RF.visiblespace( triangle_centroids, triangle_normals, @@ -143,17 +143,17 @@ def calculate_farfield( sinks, np.zeros((len(sinks), 3), dtype=np.float32), sink_normals, lengths ) sink_cloud = RF.points2pointcloud(sinks) - sink_cloud.point_data["normals"] = sink_normals + sink_cloud.point_data["Normals"] = sink_normals num_sources = len(np.asarray(aperture_coords.points)) num_sinks = len(np.asarray(sink_cloud.points)) environment_triangles=GF.mesh_conversion(antenna_solid) if project_vectors: conformal_E_vectors = EM.calculate_conformalVectors( - desired_E_axis, np.asarray(aperture_coords.point_data["normals"]), antenna_axes + desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes ) else: - if desired_E_axis.shape[0]==np.asarray(aperture_coords.point_data["normals"]).shape[0]: + if desired_E_axis.shape[0]==np.asarray(aperture_coords.point_data["Normals"]).shape[0]: conformal_E_vectors=copy.deepcopy(desired_E_axis) else: conformal_E_vectors = np.repeat( @@ -169,8 +169,8 @@ def calculate_farfield( axis=0, ) unified_normals = np.append( - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), - np.asarray(sink_cloud.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), + np.asarray(sink_cloud.point_data["Normals"]).astype(np.float32), axis=0, ) unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) @@ -199,7 +199,7 @@ def calculate_farfield( point_informationv2[0:num_sources]["pz"] = np.asarray( aperture_coords.points ).astype(np.float32)[:, 2] - normals = np.asarray(aperture_coords.point_data["normals"]) + normals = np.asarray(aperture_coords.point_data["Normals"]) point_informationv2[0:num_sources]["nx"] = normals[:, 0] point_informationv2[0:num_sources]["ny"] = normals[:, 1] point_informationv2[0:num_sources]["nz"] = normals[:, 2] @@ -242,11 +242,11 @@ def calculate_farfield( ) unified_normals = np.append( np.append( - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), - np.asarray(sink_cloud.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), + np.asarray(sink_cloud.point_data["Normals"]).astype(np.float32), axis=0, ), - np.asarray(scatter_points.point_data["normals"]).astype(np.float32), + np.asarray(scatter_points.point_data["Normals"]).astype(np.float32), axis=0, ) unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) @@ -278,13 +278,13 @@ def calculate_farfield( aperture_coords.points ).astype(np.float32)[:, 2] point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 0] point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 1] point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 2] # point_informationv2[0:num_sources]['ex']=unified_weights[0:num_sources,0] # point_informationv2[0:num_sources]['ey']=unified_weights[0:num_sources,1] @@ -314,7 +314,7 @@ def calculate_farfield( point_informationv2[(num_sources + num_sinks) :]["pz"] = np.asarray( scatter_points.points ).astype(np.float32)[:, 2] - scatter_normals = np.asarray(scatter_points.point_data["normals"]) + scatter_normals = np.asarray(scatter_points.point_data["Normals"]) point_informationv2[(num_sources + num_sinks) :]["nx"] = scatter_normals[:, 0] point_informationv2[(num_sources + num_sinks) :]["ny"] = scatter_normals[:, 1] point_informationv2[(num_sources + num_sinks) :]["nz"] = scatter_normals[:, 2] @@ -486,11 +486,11 @@ def calculate_scattering( if not multiE: if project_vectors: conformal_E_vectors = EM.calculate_conformalVectors( - desired_E_axis, np.asarray(aperture_coords.point_data["normals"]), antenna_axes + desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes ) else: print("hi from here", aperture_coords.cell_data) - if desired_E_axis.shape[0] == np.asarray(aperture_coords.point_data["normals"]).shape[0]: + if desired_E_axis.shape[0] == np.asarray(aperture_coords.point_data["Normals"]).shape[0]: conformal_E_vectors = copy.deepcopy(desired_E_axis) else: conformal_E_vectors = np.repeat( @@ -499,10 +499,10 @@ def calculate_scattering( else: if project_vectors: conformal_E_vectors = EM.calculate_conformalVectors( - desired_E_axis, np.asarray(aperture_coords.point_data["normals"]), antenna_axes + desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes ) else: - if desired_E_axis.shape[0] == np.asarray(aperture_coords.point_data["normals"]).shape[0]: + if desired_E_axis.shape[0] == np.asarray(aperture_coords.point_data["Normals"]).shape[0]: conformal_E_vectors = copy.deepcopy(desired_E_axis) else: conformal_E_vectors = np.repeat( @@ -518,8 +518,8 @@ def calculate_scattering( axis=0, ) unified_normals = np.append( - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), - np.asarray(sink_coords.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), + np.asarray(sink_coords.point_data["Normals"]).astype(np.float32), axis=0, ) unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) @@ -545,13 +545,13 @@ def calculate_scattering( aperture_coords.points ).astype(np.float32)[:, 2] point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 0] point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 1] point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 2] # set position and velocity of sinks point_informationv2[num_sources : (num_sources + num_sinks)]["px"] = np.asarray( @@ -564,13 +564,13 @@ def calculate_scattering( sink_coords.points ).astype(np.float32)[:, 2] point_informationv2[num_sources : (num_sources + num_sinks)]["nx"] = np.asarray( - sink_coords.point_data["normals"] + sink_coords.point_data["Normals"] ).astype(np.float32)[:, 0] point_informationv2[num_sources : (num_sources + num_sinks)]["ny"] = np.asarray( - sink_coords.point_data["normals"] + sink_coords.point_data["Normals"] ).astype(np.float32)[:, 1] point_informationv2[num_sources : (num_sources + num_sinks)]["nz"] = np.asarray( - sink_coords.point_data["normals"] + sink_coords.point_data["Normals"] ).astype(np.float32)[:, 2] point_informationv2[:]["ex"] = unified_weights[:, 0] @@ -591,7 +591,7 @@ def calculate_scattering( if project_vectors: conformal_E_vectors = EM.calculate_conformalVectors( desired_E_axis[0, :].reshape(1, 3), - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), ) else: conformal_E_vectors = np.repeat( @@ -603,7 +603,7 @@ def calculate_scattering( if project_vectors: conformal_E_vectors = EM.calculate_conformalVectors( desired_E_axis[0, :].reshape(1, 3), - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), ) else: if desired_E_axis.size == 3: @@ -626,11 +626,11 @@ def calculate_scattering( ) unified_normals = np.append( np.append( - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), - np.asarray(sink_coords.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), + np.asarray(sink_coords.point_data["Normals"]).astype(np.float32), axis=0, ), - np.asarray(scatter_points.point_data["normals"]).astype(np.float32), + np.asarray(scatter_points.point_data["Normals"]).astype(np.float32), axis=0, ) unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) @@ -659,13 +659,13 @@ def calculate_scattering( aperture_coords.points ).astype(np.float32)[:, 2] point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 0] point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 1] point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.point_data["normals"] + aperture_coords.point_data["Normals"] ).astype(np.float32)[:, 2] # point_informationv2[0:num_sources]['ex']=unified_weights[0:num_sources,0] # point_informationv2[0:num_sources]['ey']=unified_weights[0:num_sources,1] @@ -684,13 +684,13 @@ def calculate_scattering( # point_informationv2[num_sources:(num_sources+num_sinks)]['vy']=0.0 # point_informationv2[num_sources:(num_sources+num_sinks)]['vz']=0.0 point_informationv2[num_sources : (num_sources + num_sinks)]["nx"] = np.asarray( - sink_coords.point_data["normals"] + sink_coords.point_data["Normals"] ).astype(np.float32)[:, 0] point_informationv2[num_sources : (num_sources + num_sinks)]["ny"] = np.asarray( - sink_coords.point_data["normals"] + sink_coords.point_data["Normals"] ).astype(np.float32)[:, 1] point_informationv2[num_sources : (num_sources + num_sinks)]["nz"] = np.asarray( - sink_coords.point_data["normals"] + sink_coords.point_data["Normals"] ).astype(np.float32)[:, 2] point_informationv2[(num_sources + num_sinks) :]["px"] = np.asarray( scatter_points.points @@ -702,13 +702,13 @@ def calculate_scattering( scatter_points.points ).astype(np.float32)[:, 2] point_informationv2[(num_sources + num_sinks) :]["nx"] = np.asarray( - scatter_points.point_data["normals"] + scatter_points.point_data["Normals"] ).astype(np.float32)[:, 0] point_informationv2[(num_sources + num_sinks) :]["ny"] = np.asarray( - scatter_points.point_data["normals"] + scatter_points.point_data["Normals"] ).astype(np.float32)[:, 1] point_informationv2[(num_sources + num_sinks) :]["nz"] = np.asarray( - scatter_points.point_data["normals"] + scatter_points.point_data["Normals"] ).astype(np.float32)[:, 2] point_informationv2[:]["ex"] = unified_weights[:, 0] point_informationv2[:]["ey"] = unified_weights[:, 1] @@ -741,7 +741,7 @@ def calculate_scattering( for e_inc in range(desired_E_axis.shape[0]): conformal_E_vectors = EM.calculate_conformalVectors( desired_E_axis[e_inc, :], - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), ) unified_weights[0:num_sources, :] = conformal_E_vectors# / num_sources point_informationv2[:]["ex"] = unified_weights[:, 0] @@ -772,7 +772,7 @@ def calculate_scattering( for e_inc in range(desired_E_axis.shape[1]): conformal_E_vectors = EM.calculate_conformalVectors( desired_E_axis[e_inc, :], - np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), + np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), ) for element in range(num_sources): point_informationv2[0:num_sources]["ex"] = 0.0 diff --git a/lyceanem/raycasting/rayfunctions.py b/lyceanem/raycasting/rayfunctions.py index 3a9596a..4a594b2 100644 --- a/lyceanem/raycasting/rayfunctions.py +++ b/lyceanem/raycasting/rayfunctions.py @@ -528,10 +528,13 @@ def convertTriangles(triangle_object): if triangle_object == None: triangles = np.empty(0, dtype=base_types.triangle_t) else: - assert triangle_object.cells[0].type == "triangle", "Not a triangle mesh" - + #assert triangle_object.cells[1].type == "triangle", "Not a triangle mesh" + for idx in range(len(triangle_object.cells)): + if (triangle_object.cells[idx].type=="triangle"): + triangle_cell_index=idx + vertices = np.asarray(triangle_object.points) - tri_index = np.asarray(triangle_object.cells[0].data) + tri_index = np.asarray(triangle_object.cells[triangle_cell_index].data) triangles = np.empty(len(tri_index), dtype=base_types.triangle_t) for idx in range(len(tri_index)): triangles[idx]["v0x"] = np.single(vertices[tri_index[idx, 0], 0])