From 7a3fba8fae5bd00aa696844abe13f891840ca319 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9on=20van=20Velzen?= Date: Thu, 31 Oct 2024 16:22:09 +0100 Subject: [PATCH 1/5] Add Surface.box method --- traceon/geometry.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/traceon/geometry.py b/traceon/geometry.py index dc3e93b3..90c0fb99 100644 --- a/traceon/geometry.py +++ b/traceon/geometry.py @@ -861,6 +861,30 @@ def f(u, v): return Surface(f, length1, length2) + @staticmethod + def box(p0, p1): + x0, y0, z0 = p0 + x1, y1, z1 = p1 + + xmin, ymin, zmin = min(x0, x1), min(y0, y1), min(z0, z1) + xmax, ymax, zmax = max(x0, x1), max(y0, y1), max(z0, z1) + + path1 = Path.line([xmin, ymin, zmax], [xmax, ymin, zmax]) + path2 = Path.line([xmin, ymin, zmin], [xmax, ymin, zmin]) + path3 = Path.line([xmin, ymax, zmax], [xmax, ymax, zmax]) + path4 = Path.line([xmin, ymax, zmin], [xmax, ymax, zmin]) + + side_path = Path.line([xmin, ymin, zmin], [xmax, ymin, zmin])\ + .line_to([xmax, ymin, zmax])\ + .line_to([xmin, ymin, zmax])\ + .close() + + side_surface = side_path.extrude([0.0, ymax-ymin, 0.0]) + top = Surface.spanned_by_paths(path1, path2) + bottom = Surface.spanned_by_paths(path4, path3) + + return (top + bottom + side_surface) + @staticmethod def from_boundary_paths(p1, p2, p3, p4): path_length_p1_and_p3 = (p1.path_length + p3.path_length)/2 From 04e93d0b1ef24e9ec58b3ae234ec96a7769628fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9on=20van=20Velzen?= Date: Thu, 31 Oct 2024 16:26:12 +0100 Subject: [PATCH 2/5] Add more __str__ methods in geometry.py --- traceon/geometry.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/traceon/geometry.py b/traceon/geometry.py index 90c0fb99..c58efc75 100644 --- a/traceon/geometry.py +++ b/traceon/geometry.py @@ -732,6 +732,9 @@ def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False): return Mesh(points=points, lines=lines, physical_to_lines=physical_to_lines) + def __str__(self): + return f"" + class PathCollection(GeometricObject): """A PathCollection is a collection of `Path`. It can be created using the + operator (for example path1+path2). @@ -790,7 +793,9 @@ def extrude(self, vector): return self._map_to_surfaces(Path.extrude, vector) def extrude_by_path(self, p2): return self._map_to_surfaces(Path.extrude_by_path, p2) - + + def __str__(self): + return f"" class Surface(GeometricObject): @@ -1099,6 +1104,9 @@ def __iadd__(self, other): self.surfaces.append(other) else: self.surfaces.extend(other.surfaces) + + def __str__(self): + return f"" From e395aad809cee79da287ea74d8426a5155dc693f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9on=20van=20Velzen?= Date: Thu, 31 Oct 2024 16:26:35 +0100 Subject: [PATCH 3/5] Fix the way name is handled in geometry.py using properties --- traceon/geometry.py | 60 +++++++++++++++++++++++++++++++++------------ 1 file changed, 44 insertions(+), 16 deletions(-) diff --git a/traceon/geometry.py b/traceon/geometry.py index c58efc75..2073060a 100644 --- a/traceon/geometry.py +++ b/traceon/geometry.py @@ -683,7 +683,7 @@ def __add__(self, other): elif isinstance(other, PathCollection): return PathCollection([self] + [other.paths]) - def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False): + def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False, name=None): """Mesh the path, so it can be used in the BEM solver. Parameters @@ -724,9 +724,11 @@ def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False): 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 self.name is not None: - physical_to_lines = {self.name:np.arange(len(lines))} + if name is not None: + physical_to_lines = {name:np.arange(len(lines))} else: physical_to_lines = {} @@ -743,18 +745,29 @@ class PathCollection(GeometricObject): def __init__(self, paths): assert all([isinstance(p, Path) for p in paths]) self.paths = paths - self.name = None + self._name = None + @property + def name(self): + return self._name + + @name.setter + def name(self, name): + self._name = name + + for path in self.paths: + path.name = name + def map_points(self, fun): return PathCollection([p.map_points(fun) for p in self.paths]) - def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False): + def mesh(self, mesh_size=None, mesh_size_factor=None, higher_order=False, name=None): mesh = Mesh() + name = self.name if name is None else name + for p in self.paths: - if self.name is not None: - p.name = self.name - mesh = mesh + p.mesh(mesh_size=mesh_size, mesh_size_factor=mesh_size_factor, higher_order=higher_order) + mesh = mesh + p.mesh(mesh_size=mesh_size, mesh_size_factor=mesh_size_factor, higher_order=higher_order, name=name) return mesh @@ -1051,7 +1064,7 @@ def __add__(self, other): else: return SurfaceCollection([self] + other.surfaces) - def mesh(self, mesh_size=None, mesh_size_factor=None): + def mesh(self, mesh_size=None, mesh_size_factor=None, name=None): if mesh_size is None: path_length = min(self.path_length1, self.path_length2) @@ -1060,9 +1073,14 @@ def mesh(self, mesh_size=None, mesh_size_factor=None): if mesh_size_factor is not None: mesh_size /= sqrt(mesh_size_factor) - - return _mesh(self, mesh_size, name=self.name) + name = self.name if name is None else name + print('meshing surface with name ', name) + + return _mesh(self, mesh_size, name=name) + + def __str__(self): + return f"" class SurfaceCollection(GeometricObject): @@ -1072,7 +1090,18 @@ class SurfaceCollection(GeometricObject): def __init__(self, surfaces): assert all([isinstance(s, Surface) for s in surfaces]) self.surfaces = surfaces - self.name = None + self._name = None + + @property + def name(self): + return self._name + + @name.setter + def name(self, name): + self._name = name + + for surf in self.surfaces: + surf.name = name def map_points(self, fun): return SurfaceCollection([s.map_points(fun) for s in self.surfaces]) @@ -1080,11 +1109,10 @@ def map_points(self, fun): def mesh(self, mesh_size=None, mesh_size_factor=None, name=None): mesh = Mesh() + name = self.name if name is None else name + for s in self.surfaces: - if self.name is not None: - s.name = self.name - - mesh = mesh + s.mesh(mesh_size=mesh_size, mesh_size_factor=mesh_size_factor) + mesh = mesh + s.mesh(mesh_size=mesh_size, mesh_size_factor=mesh_size_factor, name=name) return mesh From 3ccbbba3299e5392e5e779f1f168dda5e58f65d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9on=20van=20Velzen?= Date: Thu, 31 Oct 2024 16:26:51 +0100 Subject: [PATCH 4/5] Update examples/deflector-3d.py to new API --- examples/deflector-3d.py | 71 +++++++++++++++------------------------- 1 file changed, 27 insertions(+), 44 deletions(-) diff --git a/examples/deflector-3d.py b/examples/deflector-3d.py index f6195632..9f095e03 100644 --- a/examples/deflector-3d.py +++ b/examples/deflector-3d.py @@ -21,55 +21,44 @@ # This is a 'disk' or round electrode with a hole in the middle. # The top and bottom electrodes are round grounded electrodes # which shield the microscope from the deflector fields. -def round_electrode(geom, z0, name): +def round_electrode(z0, name): - points = [ [RADIUS, 0.0, z0], [RADIUS, 0.0, z0+THICKNESS], - [RADIUS+ELECTRODE_WIDTH, 0.0, z0+THICKNESS], [RADIUS+ELECTRODE_WIDTH, 0.0, z0] ] + path = G.Path.line([RADIUS, 0.0, z0], [RADIUS, 0.0, z0+THICKNESS])\ + .line_to([RADIUS+ELECTRODE_WIDTH, 0.0, z0+THICKNESS])\ + .line_to([RADIUS+ELECTRODE_WIDTH, 0.0, z0])\ + .close() - points = [geom.add_point(p) for p in points] + surf = path.revolve_z() + surf.name = name - l1 = geom.add_line(points[0], points[1]) - l2 = geom.add_line(points[1], points[2]) - l3 = geom.add_line(points[2], points[3]) - l4 = geom.add_line(points[3], points[0]) - - revolved = G.revolve_around_optical_axis(geom, [l1, l2, l3, l4]) - geom.add_physical(revolved, 'ground') + return surf # Create a simple electrode consisting of a extruded rectangle. # A deflector consists of two rectangle electrodes 'facing' each other. -def rectangle_electrode(geom, x, z, name): +def rectangle_electrode(x, z, name): + + p0 = [x, -LENGTH_DEFLECTORS/2, z] + p1 = [x+THICKNESS, +LENGTH_DEFLECTORS/2, z+THICKNESS] - points = [ [x, -LENGTH_DEFLECTORS/2, z], [x, -LENGTH_DEFLECTORS/2, z+THICKNESS], - [x+THICKNESS, -LENGTH_DEFLECTORS/2, z+THICKNESS], [x+THICKNESS, -LENGTH_DEFLECTORS/2, z] ] + box = G.Surface.box(p0, p1) + box.name = name - poly = geom.add_polygon(points) + return box - top, extruded, lateral = geom.extrude(poly, [0.0, LENGTH_DEFLECTORS, 0.0]) - geom.add_physical(lateral, name) - geom.add_physical(top, name) - geom.add_physical(poly, name) +electrode_top = round_electrode(z0, 'ground') +defl_pos = rectangle_electrode(RADIUS, z0+THICKNESS+SPACING, 'deflector_positive') +defl_neg = rectangle_electrode(-RADIUS-THICKNESS, z0+THICKNESS+SPACING, 'deflector_negative') +electrode_bottom = round_electrode(z0+THICKNESS+SPACING+THICKNESS+SPACING, 'ground') +electrode_mesh = (electrode_top + electrode_bottom).mesh(mesh_size_factor=24) +defl_mesh = (defl_pos + defl_neg).mesh(mesh_size_factor=3) -# Create the actual geometry using the utility functions above. -with G.Geometry(G.Symmetry.THREE_D, size_from_distance=True) as geom: - round_electrode(geom, z0, 'ground') - - rectangle_electrode(geom, RADIUS, z0+THICKNESS+SPACING, 'deflector_positive') - rectangle_electrode(geom, -RADIUS-THICKNESS, z0+THICKNESS+SPACING, 'deflector_negative') - - round_electrode(geom, z0+THICKNESS+SPACING+THICKNESS+SPACING, 'ground') - - # The higher the mesh factor, the more triangles are used. This improves - # accuracy at the expense of computation time. - geom.set_mesh_size_factor(250) - - mesh = geom.generate_triangle_mesh() +mesh = electrode_mesh + defl_mesh # Show the generated triangle mesh. -P.plot_mesh(mesh, ground='green', deflector_positive='red', deflector_negative='blue') +P.plot_mesh(mesh, ground='green', deflector_positive='red', deflector_negative='blue', show_normals=True) -excitation = E.Excitation(mesh) +excitation = E.Excitation(mesh, E.Symmetry.THREE_D) # Apply the correct voltages. Here we set one deflector electrode to 5V and # the other electrode to -5V. @@ -79,16 +68,10 @@ def rectangle_electrode(geom, x, z, name): # the surface charges gives rise to a electrostatic field. field = S.solve_direct(excitation) -# But using an integration over the surface charges to calculate the electron -# trajectories is inherently slow. Instead, use an interpolation technique -# in which we use the multipole coefficients of the potential along the potential axis. -# The complicated mathematics are all abstracted away from the user. -field_axial = field.axial_derivative_interpolation(-2, 2, 200) - # An instance of the tracer class allows us to easily find the trajectories of -# electrons. Here we specify that the interpolated field should be used, and that -# the tracing should stop if the x,y value goes outside ±RADIUS/2 or the z value outside ±7 mm. -tracer = field_axial.get_tracer( [(-RADIUS/2, RADIUS/2), (-RADIUS/2, RADIUS/2), (-7, 7)] ) +# electrons. Here we specify that the tracing should stop if the x,y values +# go outside ±RADIUS/2 or the z value outside ±7 mm. +tracer = field.get_tracer( [(-RADIUS/2, RADIUS/2), (-RADIUS/2, RADIUS/2), (-7, 7)] ) r_start = np.linspace(-RADIUS/8, RADIUS/8, 5) From 2271903425e7a446d27eb047f7d3fa1adceccefe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9on=20van=20Velzen?= Date: Sun, 3 Nov 2024 16:54:13 +0100 Subject: [PATCH 5/5] Fix examples/einzel-lens.py --- examples/einzel-lens.py | 37 ++++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/examples/einzel-lens.py b/examples/einzel-lens.py index 91f4c374..92e28c1c 100644 --- a/examples/einzel-lens.py +++ b/examples/einzel-lens.py @@ -6,6 +6,7 @@ import traceon.excitation as E import traceon.plotting as P import traceon.tracing as T +from traceon.interpolation import FieldRadialAxial # Dimensions of the einzel lens. THICKNESS = 0.5 @@ -16,26 +17,28 @@ # lens is at z = 0mm. z0 = -THICKNESS - SPACING - THICKNESS/2 -with G.MEMSStack(z0=z0, zmin=-1, zmax=2, size_from_distance=True) as geom: - - # Einzel lens consists of three electrodes: - # a lens electrode sandwiched between two ground electrodes. - geom.add_electrode(RADIUS, THICKNESS, 'ground') - geom.add_spacer(THICKNESS) - geom.add_electrode(RADIUS, THICKNESS, 'lens') - geom.add_spacer(THICKNESS) - geom.add_electrode(RADIUS, THICKNESS, 'ground') - - geom.set_mesh_size_factor(80) +boundary = G.Path.line([0., 0., 1.75], [2.0, 0., 1.75])\ + .line_to([2.0, 0., -1.75]).line_to([0., 0., -1.75]) + + +margin_right = 0.1 +extent = 2.0 - margin_right + +bottom = G.Path.aperture(THICKNESS, RADIUS, extent, -THICKNESS - SPACING) +middle = G.Path.aperture(THICKNESS, RADIUS, extent) +top = G.Path.aperture(THICKNESS, RADIUS, extent, THICKNESS + SPACING) - # Actually generate the mesh, which takes the boundaries in the - # geometry and produces many line elements. - mesh = geom.generate_line_mesh(False) +boundary.name = 'boundary' +bottom.name = 'ground' +middle.name = 'lens' +top.name = 'ground' +mesh = (boundary + bottom + middle + top).mesh(mesh_size_factor=45) + # Show the generated mesh, with the given electrode colors. P.plot_mesh(mesh, ground='green', lens='blue', show_normals=True) -excitation = E.Excitation(mesh) +excitation = E.Excitation(mesh, E.Symmetry.RADIAL) # Excite the geometry, put ground at 0V and the lens electrode at 1000V. excitation.add_voltage(ground=0.0, lens=1000) @@ -49,12 +52,12 @@ # trajectories is inherently slow. Instead, use an interpolation technique # in which we use the derivatives of the potential along the potential axis. # The complicated mathematics are all abstracted away from the user. -field_axial = field.axial_derivative_interpolation(-1, 2, 150) +field_axial = FieldRadialAxial(field, -1.5, 1.5, 150) # Plot the potential along the optical axis to show that the interpolated # potential is very close to the potential found by an integration over the # surface charge. -z = np.linspace(-1, 2, 150) +z = np.linspace(-1.5, 1.5, 150) pot = [field.potential_at_point(np.array([0.0, z_])) for z_ in z] pot_axial = [field_axial.potential_at_point(np.array([0.0, z_])) for z_ in z]