+
+ Expand source code
+
+ class Path(GeometricObject):
+ """A path is a mapping from a number in the range [0, path_length] to a three dimensional point. Note that `Path` is a
+ subclass of `traceon.mesher.GeometricObject`, and therefore can be easily moved and rotated."""
+
+ def __init__(self, fun, path_length, breakpoints=[], name=None):
+ # Assumption: fun takes in p, the path length
+ # and returns the point on the path
+ self.fun = fun
+ self.path_length = path_length
+ self.breakpoints = breakpoints
+ self.name = name
+
+ @staticmethod
+ def from_irregular_function(to_point, N=100, breakpoints=[]):
+ """Construct a path from a function that is of the form u -> point, where 0 <= u <= 1.
+ The length of the path is determined by integration.
+
+ Parameters
+ ---------------------------------
+ to_point: callable
+ A function accepting a number in the range [0, 1] and returns a the dimensional point.
+ N: int
+ Number of samples to use in the cubic spline interpolation.
+ breakpoints: float iterable
+ Points (0 <= u <= 1) on the path where the function is non-differentiable. These points
+ are always included in the resulting mesh.
+
+ Returns
+ ---------------------------------
+ Path"""
+
+ # path length = integrate |f'(x)|
+ fun = lambda u: np.array(to_point(u))
+
+ u = np.linspace(0, 1, N)
+ samples = CubicSpline(u, [fun(u_) for u_ in u])
+ derivatives = samples.derivative()(u)
+ norm_derivatives = np.linalg.norm(derivatives, axis=1)
+ path_lengths = CubicSpline(u, norm_derivatives).antiderivative()(u)
+ interpolation = CubicSpline(path_lengths, u) # Path length to [0,1]
+
+ path_length = path_lengths[-1]
+
+ return Path(lambda pl: fun(interpolation(pl)), path_length, breakpoints=[b*path_length for b in breakpoints])
+
+ @staticmethod
+ def spline_through_points(points, N=100):
+ """Construct a path by fitting a cubic spline through the given points.
+
+ Parameters
+ -------------------------
+ points: (N, 3) ndarray of float
+ Three dimensional points through which the spline is fitted.
+
+ Returns
+ -------------------------
+ Path"""
+
+ x = np.linspace(0, 1, len(points))
+ interp = CubicSpline(x, points)
+ return Path.from_irregular_function(interp, N=N)
+
+ def average(self, fun):
+ """Average a function along the path, by integrating 1/l * fun(path(l)) with 0 <= l <= path length.
+
+ Parameters
+ --------------------------
+ fun: callable (3,) -> float
+ A function taking a three dimensional point and returning a float.
+
+ Returns
+ -------------------------
+ float
+
+ The average value of the function along the point."""
+ return quad(lambda s: fun(self(s)), 0, self.path_length, points=self.breakpoints)[0]/self.path_length
+
+ def map_points(self, fun):
+ """Return a new function by mapping a function over points along the path (see `traceon.mesher.GeometricObject`).
+ The path length is assumed to stay the same after this operation.
+
+ Parameters
+ ----------------------------
+ fun: callable (3,) -> (3,)
+ Function taking three dimensional points and returning three dimensional points.
+
+ Returns
+ ---------------------------
+
+ Path"""
+ return Path(lambda u: fun(self(u)), self.path_length, self.breakpoints, name=self.name)
+
+ def __call__(self, t):
+ """Evaluate a point along the path.
+
+ Parameters
+ ------------------------
+ t: float
+ The length along the path.
+
+ Returns
+ ------------------------
+ (3,) float
+
+ Three dimensional point."""
+ return self.fun(t)
+
+ def is_closed(self):
+ """Determine whether the path is closed, by comparing the starting and endpoint.
+
+ Returns
+ ----------------------
+ bool: True if the path is closed, False otherwise."""
+ return _points_close(self.starting_point(), self.endpoint())
+
+ def add_phase(self, l):
+ """Add a phase to a closed path. A path is closed when the starting point is equal to the
+ end point. A phase of length l means that the path starts 'further down' the closed path.
+
+ Parameters
+ --------------------
+ l: float
+ The phase (expressed as a path length). The resulting path starts l distance along the
+ original path.
+
+ Returns
+ --------------------
+ Path"""
+ assert self.is_closed()
+
+ def fun(u):
+ return self( (l + u) % self.path_length )
+
+ return Path(fun, self.path_length, sorted([(b-l)%self.path_length for b in self.breakpoints + [0.]]), name=self.name)
+
+ def __rshift__(self, other):
+ """Combine two paths to create a single path. The endpoint of the first path needs
+ to match the starting point of the second path. This common point is marked as a breakpoint and
+ always included in the mesh. To use this function use the right shift operator (p1 >> p2).
+
+ Parameters
+ -----------------------
+ other: Path
+ The second path, to extend the current path.
+
+ Returns
+ -----------------------
+ Path"""
+
+ assert isinstance(other, Path), "Exteding path with object that is not actually a Path"
+
+ assert _points_close(self.endpoint(), other.starting_point())
+
+ total = self.path_length + other.path_length
+
+ def f(t):
+ assert 0 <= t <= total
+
+ if t <= self.path_length:
+ return self(t)
+ else:
+ return other(t - self.path_length)
+
+ return Path(f, total, self.breakpoints + [self.path_length] + other.breakpoints, name=self.name)
+
+ def starting_point(self):
+ """Returns the starting point of the path.
+
+ Returns
+ ---------------------
+ (3,) float
+
+ The starting point of the path."""
+ return self(0.)
+
+ def middle_point(self):
+ """Returns the midpoint of the path (in terms of length along the path.)
+
+ Returns
+ ----------------------
+ (3,) float
+
+ The point at the middle of the path."""
+ return self(self.path_length/2)
+
+ def endpoint(self):
+ """Returns the endpoint of the path.
+
+ Returns
+ ------------------------
+ (3,) float
+
+ The endpoint of the path."""
+ return self(self.path_length)
+
+ def line_to(self, point):
+ """Extend the current path by a line from the current endpoint to the given point.
+ The given point is marked a breakpoint.
+
+ Parameters
+ ----------------------
+ point: (3,) float
+ The new endpoint.
+
+ Returns
+ ---------------------
+ Path"""
+ warnings.warn("line_to() is deprecated and will be removed in version 0.8.0."
+ "Use extend_with_line() instead.",
+ DeprecationWarning,
+ stacklevel=2)
+
+ point = np.array(point)
+ assert point.shape == (3,), "Please supply a three dimensional point to .line_to(...)"
+ l = Path.line(self.endpoint(), point)
+ return self >> l
+
+ def extend_with_line(self, point):
+ """Extend the current path by a line from the current endpoint to the given point.
+ The given point is marked a breakpoint.
+
+ Parameters
+ ----------------------
+ point: (3,) float
+ The new endpoint.
+
+ Returns
+ ---------------------
+ Path"""
+ point = np.array(point)
+ assert point.shape == (3,), "Please supply a three dimensional point to .extend_with_line(...)"
+ l = Path.line(self.endpoint(), point)
+ return self >> l
+
+ @staticmethod
+ def circle_xz(x0, z0, radius, angle=2*pi):
+ """Returns (part of) a circle in the XZ plane around the x-axis. Starting on the positive x-axis.
+
+ Parameters
+ --------------------------------
+ x0: float
+ x-coordinate of the center of the circle
+ z0: float
+ z-coordinate of the center of the circle
+ radius: float
+ radius of the circle
+ angle: float
+ The circumference of the circle in radians. The default of 2*pi gives a full circle.
+
+ Returns
+ ---------------------------------
+ Path"""
+ def f(u):
+ theta = u / radius
+ return np.array([radius*cos(theta), 0., radius*sin(theta)])
+ return Path(f, angle*radius).move(dx=x0, dz=z0)
+
+ @staticmethod
+ def circle_yz(y0, z0, radius, angle=2*pi):
+ """Returns (part of) a circle in the YZ plane around the x-axis. Starting on the positive y-axis.
+
+ Parameters
+ --------------------------------
+ y0: float
+ x-coordinate of the center of the circle
+ z0: float
+ z-coordinate of the center of the circle
+ radius: float
+ radius of the circle
+ angle: float
+ The circumference of the circle in radians. The default of 2*pi gives a full circle.
+
+ Returns
+ ---------------------------------
+ Path"""
+ def f(u):
+ theta = u / radius
+ return np.array([0., radius*cos(theta), radius*sin(theta)])
+ return Path(f, angle*radius).move(dy=y0, dz=z0)
+
+ @staticmethod
+ def circle_xy(x0, y0, radius, angle=2*pi):
+ """Returns (part of) a circle in the XY plane around the z-axis. Starting on the positive X-axis.
+
+ Parameters
+ --------------------------------
+ x0: float
+ x-coordinate of the center of the circle
+ y0: float
+ y-coordinate of the center of the circle
+ radius: float
+ radius of the circle
+ angle: float
+ The circumference of the circle in radians. The default of 2*pi gives a full circle.
+
+ Returns
+ ---------------------------------
+ Path"""
+ def f(u):
+ theta = u / radius
+ return np.array([radius*cos(theta), radius*sin(theta), 0.])
+ return Path(f, angle*radius).move(dx=x0, dy=y0)
+
+ def arc_to(self, center, end, reverse=False):
+ """Extend the current path using an arc.
+
+ Parameters
+ ----------------------------
+ center: (3,) float
+ The center point of the arc.
+ end: (3,) float
+ The endpoint of the arc, shoud lie on a circle determined
+ by the given centerpoint and the current endpoint.
+
+ Returns
+ -----------------------------
+ Path"""
+ warnings.warn("arc_to() is deprecated and will be removed in version 0.8.0."
+ "Use extend_with_arc() instead.",
+ DeprecationWarning,
+ stacklevel=2)
+
+ start = self.endpoint()
+ return self >> Path.arc(center, start, end, reverse=reverse)
+
+ def extend_with_arc(self, center, end, reverse=False):
+ """Extend the current path using an arc.
+
+ Parameters
+ ----------------------------
+ center: (3,) float
+ The center point of the arc.
+ end: (3,) float
+ The endpoint of the arc, shoud lie on a circle determined
+ by the given centerpoint and the current endpoint.
+
+ Returns
+ -----------------------------
+ Path"""
+ start = self.endpoint()
+ return self >> Path.arc(center, start, end, reverse=reverse)
+
+ @staticmethod
+ def arc(center, start, end, reverse=False):
+ """Return an arc by specifying the center, start and end point.
+
+ Parameters
+ ----------------------------
+ center: (3,) float
+ The center point of the arc.
+ start: (3,) float
+ The start point of the arc.
+ end: (3,) float
+ The endpoint of the arc.
+
+ Returns
+ ----------------------------
+ Path"""
+ start_arr, center_arr, end_arr = np.array(start), np.array(center), np.array(end)
+
+ x_unit = start_arr - center_arr
+ x_unit /= np.linalg.norm(x_unit)
+
+ vector = end_arr - center_arr
+
+ y_unit = vector - np.dot(vector, x_unit) * x_unit
+ y_unit /= np.linalg.norm(y_unit)
+
+ radius = np.linalg.norm(start_arr - center_arr)
+ theta_max = atan2(np.dot(vector, y_unit), np.dot(vector, x_unit))
+
+ if reverse:
+ theta_max = theta_max - 2*pi
+
+ path_length = abs(theta_max * radius)
+
+ def f(l):
+ theta = l/path_length * theta_max
+ return center + radius*cos(theta)*x_unit + radius*sin(theta)*y_unit
+
+ return Path(f, path_length)
+
+ def revolve_x(self, angle=2*pi):
+ """Create a surface by revolving the path anti-clockwise around the x-axis.
+
+ Parameters
+ -----------------------
+ angle: float
+ The angle by which to revolve. THe default 2*pi gives a full revolution.
+
+ Returns
+ -----------------------
+ Surface"""
+
+ r_avg = self.average(lambda p: sqrt(p[1]**2 + p[2]**2))
+ length2 = 2*pi*r_avg
+
+ def f(u, v):
+ p = self(u)
+ theta = atan2(p[2], p[1])
+ r = sqrt(p[1]**2 + p[2]**2)
+ return np.array([p[0], r*cos(theta + v/length2*angle), r*sin(theta + v/length2*angle)])
+
+ return Surface(f, self.path_length, length2, self.breakpoints, name=self.name)
+
+ def revolve_y(self, angle=2*pi):
+ """Create a surface by revolving the path anti-clockwise around the y-axis.
+
+ Parameters
+ -----------------------
+ angle: float
+ The angle by which to revolve. THe default 2*pi gives a full revolution.
+
+ Returns
+ -----------------------
+ Surface"""
+
+ r_avg = self.average(lambda p: sqrt(p[0]**2 + p[2]**2))
+ length2 = 2*pi*r_avg
+
+ def f(u, v):
+ p = self(u)
+ theta = atan2(p[2], p[0])
+ r = sqrt(p[0]*p[0] + p[2]*p[2])
+ return np.array([r*cos(theta + v/length2*angle), p[1], r*sin(theta + v/length2*angle)])
+
+ return Surface(f, self.path_length, length2, self.breakpoints, name=self.name)
+
+ def revolve_z(self, angle=2*pi):
+ """Create a surface by revolving the path anti-clockwise around the z-axis.
+
+ Parameters
+ -----------------------
+ angle: float
+ The angle by which to revolve. THe default 2*pi gives a full revolution.
+
+ Returns
+ -----------------------
+ Surface"""
+
+ r_avg = self.average(lambda p: sqrt(p[0]**2 + p[1]**2))
+ length2 = 2*pi*r_avg
+
+ def f(u, v):
+ p = self(u)
+ theta = atan2(p[1], p[0])
+ r = sqrt(p[0]*p[0] + p[1]*p[1])
+ return np.array([r*cos(theta + v/length2*angle), r*sin(theta + v/length2*angle), p[2]])
+
+ return Surface(f, self.path_length, length2, self.breakpoints, name=self.name)
+
+ def extrude(self, vector):
+ """Create a surface by extruding the path along a vector. The vector gives both
+ the length and the direction of the extrusion.
+
+ Parameters
+ -------------------------
+ vector: (3,) float
+ The direction and length (norm of the vector) to extrude by.
+
+ Returns
+ -------------------------
+ Surface"""
+ vector = np.array(vector)
+ length = np.linalg.norm(vector)
+
+ def f(u, v):
+ return self(u) + v/length*vector
+
+ return Surface(f, self.path_length, length, self.breakpoints, name=self.name)
+
+ def extrude_by_path(self, p2):
+ """Create a surface by extruding the path along a second path. The second
+ path does not need to start along the first path. Imagine the surface created
+ by moving the first path along the second path.
+
+ Parameters
+ -------------------------
+ p2: Path
+ The (second) path defining the extrusion.
+
+ Returns
+ ------------------------
+ Surface"""
+ p0 = p2.starting_point()
+
+ def f(u, v):
+ return self(u) + p2(v) - p0
+
+ return Surface(f, self.path_length, p2.path_length, self.breakpoints, p2.breakpoints, name=self.name)
+
+ def close(self):
+ """Close the path, by making a straight line to the starting point.
+
+ Returns
+ -------------------
+ Path"""
+ return self.extend_with_line(self.starting_point())
+
+ @staticmethod
+ def ellipse(major, minor):
+ """Create a path along the outline of an ellipse. The ellipse lies
+ in the XY plane, and the path starts on the positive x-axis.
+
+ Parameters
+ ---------------------------
+ major: float
+ The major axis of the ellipse (lies along the x-axis).
+ minor: float
+ The minor axis of the ellipse (lies along the y-axis).
+
+ Returns
+ ---------------------------
+ Path"""
+ # Crazy enough there is no closed formula
+ # to go from path length to a point on the ellipse.
+ # So we have to use `from_irregular_function`
+ def f(u):
+ return np.array([major*cos(2*pi*u), minor*sin(2*pi*u), 0.])
+ return Path.from_irregular_function(f)
+
+ @staticmethod
+ def line(from_, to):
+ """Create a straight line between two points.
+
+ Parameters
+ ------------------------------
+ from_: (3,) float
+ The starting point of the path.
+ to: (3,) float
+ The endpoint of the path.
+
+ Returns
+ ---------------------------
+ Path"""
+ from_, to = np.array(from_), np.array(to)
+ length = np.linalg.norm(from_ - to)
+ return Path(lambda pl: (1-pl/length)*from_ + pl/length*to, length)
+
+ def cut(self, length):
+ """Cut the path in two at a specific length along the path.
+
+ Parameters
+ --------------------------------------
+ length: float
+ The length along the path at which to cut.
+
+ Returns
+ -------------------------------------
+ (Path, Path)
+
+ A tuple containing two paths. The first path contains the path upto length, while the second path contains the rest."""
+ return (Path(self.fun, length, [b for b in self.breakpoints if b <= length], name=self.name),
+ Path(lambda l: self.fun(l + length), self.path_length - length, [b - length for b in self.breakpoints if b >= length], name=self.name))
+
+ @staticmethod
+ def rectangle_xz(xmin, xmax, zmin, zmax):
+ """Create a rectangle in the XZ plane. The path starts at (xmin, 0, zmin), and is
+ counter clockwise around the y-axis.
+
+ Parameters
+ ------------------------
+ xmin: float
+ Minimum x-coordinate of the corner points.
+ xmax: float
+ Maximum x-coordinate of the corner points.
+ zmin: float
+ Minimum z-coordinate of the corner points.
+ zmax: float
+ Maximum z-coordinate of the corner points.
+
+ Returns
+ -----------------------
+ Path"""
+ return Path.line([xmin, 0., zmin], [xmax, 0, zmin]) \
+ .extend_with_line([xmax, 0, zmax]).extend_with_line([xmin, 0., zmax]).close()
+
+ @staticmethod
+ def rectangle_yz(ymin, ymax, zmin, zmax):
+ """Create a rectangle in the YZ plane. The path starts at (0, ymin, zmin), and is
+ counter clockwise around the x-axis.
+
+ Parameters
+ ------------------------
+ ymin: float
+ Minimum y-coordinate of the corner points.
+ ymax: float
+ Maximum y-coordinate of the corner points.
+ zmin: float
+ Minimum z-coordinate of the corner points.
+ zmax: float
+ Maximum z-coordinate of the corner points.
+
+ Returns
+ -----------------------
+ Path"""
+
+ return Path.line([0., ymin, zmin], [0, ymin, zmax]) \
+ .extend_with_line([0., ymax, zmax]).extend_with_line([0., ymax, zmin]).close()
+
+ @staticmethod
+ def rectangle_xy(xmin, xmax, ymin, ymax):
+ """Create a rectangle in the XY plane. The path starts at (xmin, ymin, 0), and is
+ counter clockwise around the z-axis.
+
+ Parameters
+ ------------------------
+ xmin: float
+ Minimum x-coordinate of the corner points.
+ xmax: float
+ Maximum x-coordinate of the corner points.
+ ymin: float
+ Minimum y-coordinate of the corner points.
+ ymax: float
+ Maximum y-coordinate of the corner points.
+
+ Returns
+ -----------------------
+ Path"""
+ return Path.line([xmin, ymin, 0.], [xmin, ymax, 0.]) \
+ .extend_with_line([xmax, ymax, 0.]).extend_with_line([xmax, ymin, 0.]).close()
+
+ @staticmethod
+ def aperture(height, radius, extent, z=0.):
+ """Create an 'aperture'. Note that in a radially symmetric geometry
+ an aperture is basically a rectangle with the right side 'open'. Revolving
+ this path around the z-axis would generate a cylindircal hole in the center.
+ This is the most basic model of an aperture.
+
+ Parameters
+ ------------------------
+ height: float
+ The height of the aperture
+ radius: float
+ The radius of the aperture hole (distance to the z-axis)
+ extent: float
+ The maximum x value
+ z: float
+ The z-coordinate of the center of the aperture
+
+ Returns
+ ------------------------
+ Path"""
+ return Path.line([extent, 0., -height/2], [radius, 0., -height/2])\
+ .extend_with_line([radius, 0., height/2]).extend_with_line([extent, 0., height/2]).move(dz=z)
+
+ @staticmethod
+ def polar_arc(radius, angle, start, direction, plane_normal=[0,1,0]):
+ """Return an arc specified by polar coordinates. The arc lies in a plane defined by the
+ provided normal vector and curves from the start point in the specified direction
+ counterclockwise around the normal.
+
+ Parameters
+ ---------------------------
+ radius : float
+ The radius of the arc.
+ angle : float
+ The angle subtended by the arc (in radians)
+ start: (3,) float
+ The start point of the arc
+ plane_normal : (3,) float
+ The normal vector of the plane containing the arc
+ direction : (3,) float
+ A tangent of the arc at the starting point.
+ Must lie in the specified plane. Does not need to be normalized.
+ Returns
+ ----------------------------
+ Path"""
+ start = np.array(start, dtype=float)
+ plane_normal = np.array(plane_normal, dtype=float)
+ direction = np.array(direction, dtype=float)
+
+ direction_unit = direction / np.linalg.norm(direction)
+ plane_normal_unit = plane_normal / np.linalg.norm(plane_normal)
+
+ if not np.isclose(np.dot(direction_unit, plane_normal_unit), 0., atol=1e-7):
+ corrected_direction = direction - np.dot(direction, plane_normal_unit) * plane_normal_unit
+ raise AssertionError(
+ f"The provided direction {direction} does not lie in the specified plane. \n"
+ f"The closed valid direction is {np.round(corrected_direction, 10)}.")
+
+ if angle < 0:
+ direction, angle = -direction, -angle
+
+ center = start - radius * np.cross(direction, plane_normal)
+ center_to_start = start - center
+
+ def f(l):
+ theta = l/radius
+ return center + np.cos(theta) * center_to_start + np.sin(theta)*np.cross(plane_normal, center_to_start)
+
+ return Path(f, radius*angle)
+
+ def extend_with_polar_arc(self, radius, angle, plane_normal=[0, 1, 0]):
+ """Extend the current path by a smooth arc using polar coordinates.
+ The arc is defined by a specified radius and angle and rotates counterclockwise
+ around around the normal that defines the arcing plane.
+
+ Parameters
+ ---------------------------
+ radius : float
+ The radius of the arc
+ angle : float
+ The angle subtended by the arc (in radians)
+ plane_normal : (3,) float
+ The normal vector of the plane containing the arc
+
+ Returns
+ ----------------------------
+ Path"""
+ plane_normal = np.array(plane_normal, dtype=float)
+ start_point = self.endpoint()
+ direction = self.velocity_vector(self.path_length)
+
+ plane_normal_unit = plane_normal / np.linalg.norm(plane_normal)
+ direction_unit = direction / np.linalg.norm(direction)
+
+ if not np.isclose(np.dot(plane_normal_unit, direction_unit), 0,atol=1e-7):
+ corrected_normal = plane_normal - np.dot(direction_unit, plane_normal) * direction_unit
+ raise AssertionError(
+ f"The provided plane normal {plane_normal} is not orthogonal to the direction {direction} \n"
+ f"of the path at the endpoint so no smooth arc can be made. The closest valid normal is "
+ f"{np.round(corrected_normal, 10)}.")
+
+ return self >> Path.polar_arc(radius, angle, start_point, direction, plane_normal)
+
+ def reverse(self):
+ """Generate a reversed version of the current path.
+ The reversed path is created by inverting the traversal direction,
+ such that the start becomes the end and vice versa.
+
+ Returns
+ ----------------------------
+ Path"""
+ return Path(lambda t: self(self.path_length-t), self.path_length,
+ [self.path_length - b for b in self.breakpoints], self.name)
+
+ def velocity_vector(self, t):
+ """Calculate the velocity (tangent) vector at a specific point on the path
+ using cubic spline interpolation.
+
+ Parameters
+ ----------------------------
+ t : float
+ The point on the path at which to calculate the velocity
+ num_splines : int
+ The number of samples used for cubic spline interpolation
+
+ Returns
+ ----------------------------
+ (3,) np.ndarray of float"""
+
+ samples = np.linspace(t - self.path_length*1e-3, t + self.path_length*1e-3, 7) # Odd number to include t
+ samples_on_path = [s for s in samples if 0 <= s <= self.path_length]
+ assert len(samples_on_path), "Please supply a point that lies on the path"
+ return CubicSpline(samples_on_path, [self(s) for s in samples_on_path])(t, nu=1)
+
+
+ def __add__(self, other):
+ """Add two paths to create a PathCollection. Note that a PathCollection supports
+ a subset of the methods of Path (for example, movement, rotation and meshing). Use
+ the + operator to combine paths into a path collection: path1 + path2 + path3.
+
+ Returns
+ -------------------------
+ PathCollection"""
+
+ if isinstance(other, Path):
+ return PathCollection([self, other])
+
+ if isinstance(other, PathCollection):
+ return PathCollection([self] + [other.paths])
+
+ return NotImplemented
+
+ def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False, name=None, ensure_outward_normals=True):
+ """Mesh the path, so it can be used in the BEM solver. The result of meshing a path
+ are (possibly curved) line elements.
+
+ Parameters
+ --------------------------
+ mesh_size: float
+ Determines amount of elements in the mesh. A smaller
+ mesh size leads to more elements.
+ mesh_size_factor: float
+ Alternative way to specify the mesh size, which scales
+ with the dimensions of the geometry, and therefore more
+ easily translates between different geometries.
+ higher_order: bool
+ Whether to generate a higher order mesh. A higher order
+ produces curved line elements (determined by 4 points on
+ each curved element). The BEM solver supports higher order
+ elements in radial symmetric geometries only.
+ name: str
+ Assign this name to the mesh, instead of the name value assinged to Surface.name
+
+ Returns
+ ----------------------------
+ `traceon.mesher.Mesh`"""
+ u = discretize_path(self.path_length, self.breakpoints, mesh_size, mesh_size_factor, N_factor=3 if higher_order else 1)
+
+ N = len(u)
+ points = np.zeros( (N, 3) )
+
+ for i in range(N):
+ points[i] = self(u[i])
+
+ if not higher_order:
+ lines = np.array([np.arange(N-1), np.arange(1, N)]).T
+ else:
+ assert N % 3 == 1
+ r = np.arange(N)
+ p0 = r[0:-1:3]
+ p1 = r[3::3]
+ p2 = r[1::3]
+ p3 = r[2::3]
+ lines = np.array([p0, p1, p2, p3]).T
+
+ assert lines.dtype == np.int64 or lines.dtype == np.int32
+
+ name = self.name if name is None else name
+
+ if name is not None:
+ physical_to_lines = {name:np.arange(len(lines))}
+ else:
+ physical_to_lines = {}
+
+ return Mesh(points=points, lines=lines, physical_to_lines=physical_to_lines, ensure_outward_normals=ensure_outward_normals)
+
+ def __str__(self):
+ return f"<Path name:{self.name}, length:{self.path_length:.1e}, number of breakpoints:{len(self.breakpoints)}>"
+
+
+