Merge branch 'fix-examples'

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)
 

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]
 

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,14 +724,19 @@ 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 = {}
 
         return Mesh(points=points, lines=lines, physical_to_lines=physical_to_lines)
 
+    def __str__(self):
+        return f"<Path name:{self.name}, length:{self.path_length:.1e}, number of breakpoints:{len(self.breakpoints)}>"
+
 
 class PathCollection(GeometricObject):
     """A PathCollection is a collection of `Path`. It can be created using the + operator (for example path1+path2).
@@ -740,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
 
@@ -790,7 +806,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"<PathCollection with {len(self.paths)} surfaces, name: {self.name}>"
 
 
 class Surface(GeometricObject):
@@ -861,6 +879,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
@@ -1022,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)
@@ -1031,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"<Surface with name: {self.name}>"
 
 
 class SurfaceCollection(GeometricObject):
@@ -1043,19 +1090,29 @@ 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])
 
     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
 
@@ -1075,6 +1132,9 @@ def __iadd__(self, other):
             self.surfaces.append(other)
         else:
             self.surfaces.extend(other.surfaces)
+
+    def __str__(self):
+        return f"<SurfaceCollection with {len(self.surfaces)} surfaces, name: {self.name}>"