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 + + +