Module traceon.mesher

Classes

class GeometricObject

The Mesh class (and the classes defined in traceon.geometry) are subclasses of GeometricObject. This means that they all can be moved, rotated, mirrored.

Expand source code

class GeometricObject(ABC):
    """The Mesh class (and the classes defined in `traceon.geometry`) are subclasses
    of GeometricObject. This means that they all can be moved, rotated, mirrored."""
    
    @abstractmethod
    def map_points(self, fun: Callable[[np.ndarray], np.ndarray]) -> Any:
        """Create a new geometric object, by mapping each point by a function.
        
        Parameters
        -------------------------
        fun: (3,) float -> (3,) float
            Function taking a three dimensional point and returning a 
            three dimensional point.

        Returns
        ------------------------
        GeometricObject

        This function returns the same type as the object on which this method was called."""
        ...
    
    def move(self, dx=0., dy=0., dz=0.):
        """Move along x, y or z axis.

        Parameters
        ---------------------------
        dx: float
            Amount to move along the x-axis.
        dy: float
            Amount to move along the y-axis.
        dz: float
            Amount to move along the z-axis.

        Returns
        ---------------------------
        GeometricObject
        
        This function returns the same type as the object on which this method was called."""
    
        assert all([isinstance(d, float) or isinstance(d, int) for d in [dx, dy, dz]])
        return self.map_points(lambda p: p + np.array([dx, dy, dz]))
     
    def rotate(self, Rx=0., Ry=0., Rz=0., origin=[0., 0., 0.]):
        """Rotate counter clockwise around the x, y or z axis. Only one axis supported at the same time
        (rotations do not commute).

        Parameters
        ------------------------------------
        Rx: float
            Amount to rotate around the x-axis (radians).
        Ry: float
            Amount to rotate around the y-axis (radians).
        Rz: float
            Amount to rotate around the z-axis (radians).
        origin: (3,) float
            Point around which to rotate, which is the origin by default.

        Returns
        --------------------------------------
        GeometricObject
        
        This function returns the same type as the object on which this method was called."""
        
        assert sum([Rx==0., Ry==0., Rz==0.]) >= 2, "Only supply one axis of rotation"
        origin = np.array(origin)
        assert origin.shape == (3,), "Please supply a 3D point for origin"
         
        if Rx != 0.:
            matrix = np.array([[1, 0, 0],
                [0, np.cos(Rx), -np.sin(Rx)],
                [0, np.sin(Rx), np.cos(Rx)]])
        elif Ry != 0.:
            matrix = np.array([[np.cos(Ry), 0, np.sin(Ry)],
                [0, 1, 0],
                [-np.sin(Ry), 0, np.cos(Ry)]])
        elif Rz != 0.:
            matrix = np.array([[np.cos(Rz), -np.sin(Rz), 0],
                [np.sin(Rz), np.cos(Rz), 0],
                [0, 0, 1]])

        return self.map_points(lambda p: origin + matrix @ (p - origin))

    def mirror_xz(self):
        """Mirror object in the XZ plane.

        Returns
        --------------------------------------
        GeometricObject
        
        This function returns the same type as the object on which this method was called."""
        return self.map_points(lambda p: np.array([p[0], -p[1], p[2]]))
     
    def mirror_yz(self):
        """Mirror object in the YZ plane.

        Returns
        --------------------------------------
        GeometricObject
        This function returns the same type as the object on which this method was called."""
        return self.map_points(lambda p: np.array([-p[0], p[1], p[2]]))
    
    def mirror_xy(self):
        """Mirror object in the XY plane.

        Returns
        --------------------------------------
        GeometricObject
        
        This function returns the same type as the object on which this method was called."""
        return self.map_points(lambda p: np.array([p[0], p[1], -p[2]]))

Ancestors

• abc.ABC

Subclasses

• Path
• PathCollection
• Surface
• SurfaceCollection
• Mesh

Methods

def map_points(self, fun: Callable[[numpy.ndarray], numpy.ndarray]) ‑> Any

Create a new geometric object, by mapping each point by a function.

Parameters

fun : (3,) float -> (3,) float

Function taking a three dimensional point and returning a three dimensional point.

Returns

GeometricObject

 

This function returns the same type as the object on which this method was called.

def mirror_xy(self)

Mirror object in the XY plane.

Returns

GeometricObject

 

This function returns the same type as the object on which this method was called.

def mirror_xz(self)

Mirror object in the XZ plane.

Returns

GeometricObject

 

This function returns the same type as the object on which this method was called.

def mirror_yz(self)

Mirror object in the YZ plane.

Returns

GeometricObject

 

This function returns the same type as the object on which this method was called.

def move(self, dx=0.0, dy=0.0, dz=0.0)

Move along x, y or z axis.

Parameters

dx : float

Amount to move along the x-axis.

dy : float

Amount to move along the y-axis.

dz : float

Amount to move along the z-axis.

Returns

GeometricObject

 

This function returns the same type as the object on which this method was called.

def rotate(self, Rx=0.0, Ry=0.0, Rz=0.0, origin=[0.0, 0.0, 0.0])

Rotate counter clockwise around the x, y or z axis. Only one axis supported at the same time (rotations do not commute).

Parameters

Rx : float

Amount to rotate around the x-axis (radians).

Ry : float

Amount to rotate around the y-axis (radians).

Rz : float

Amount to rotate around the z-axis (radians).

origin : (3,) float

Point around which to rotate, which is the origin by default.

Returns

GeometricObject

 

This function returns the same type as the object on which this method was called.

class Mesh (points=[], lines=[], triangles=[], physical_to_lines={}, physical_to_triangles={}, ensure_outward_normals=True)

Mesh containing lines and triangles. Groups of lines or triangles can be named. These names are later used to apply the correct excitation. Line elements can be curved (or 'higher order'), in which case they are represented by four points per element. Note that Mesh is a subclass of GeometricObject, and therefore can be easily moved and rotated.

Expand source code

class Mesh(Saveable, GeometricObject):
    """Mesh containing lines and triangles. Groups of lines or triangles can be named. These
    names are later used to apply the correct excitation. Line elements can be curved (or 'higher order'), 
    in which case they are represented by four points per element.  Note that `Mesh` is a subclass of
    `traceon.mesher.GeometricObject`, and therefore can be easily moved and rotated."""
     
    def __init__(self,
            points=[],
            lines=[],
            triangles=[],
            physical_to_lines={},
            physical_to_triangles={},
            ensure_outward_normals=True):
        
        # Ensure the correct shape even if empty arrays
        if len(points):
            self.points = np.array(points, dtype=np.float64)
        else:
            self.points = np.empty((0,3), dtype=np.float64)
         
        if len(lines) or (isinstance(lines, np.ndarray) and len(lines.shape)==2):
            self.lines = np.array(lines, dtype=np.uint64)
        else:
            self.lines = np.empty((0,2), dtype=np.uint64)
    
        if len(triangles):
            self.triangles = np.array(triangles, dtype=np.uint64)
        else:
            self.triangles = np.empty((0, 3), dtype=np.uint64)
         
        self.physical_to_lines = physical_to_lines.copy()
        self.physical_to_triangles = physical_to_triangles.copy()

        self._remove_degenerate_triangles()
        self._deduplicate_points()

        if ensure_outward_normals:
            for el in self.get_electrodes():
                self.ensure_outward_normals(el)
         
        assert np.all( (0 <= self.lines) & (self.lines < len(self.points)) ), "Lines reference points outside points array"
        assert np.all( (0 <= self.triangles) & (self.triangles < len(self.points)) ), "Triangles reference points outside points array"
        assert np.all([np.all( (0 <= group) & (group < len(self.lines)) ) for group in self.physical_to_lines.values()])
        assert np.all([np.all( (0 <= group) & (group < len(self.triangles)) ) for group in self.physical_to_triangles.values()])
        assert not len(self.lines) or self.lines.shape[1] in [2,4], "Lines should contain either 2 or 4 points."
        assert not len(self.triangles) or self.triangles.shape[1] in [3,6], "Triangles should contain either 3 or 6 points"
    
    def is_higher_order(self):
        """Whether the mesh contains higher order elements.

        Returns
        ----------------------------
        bool"""
        return isinstance(self.lines, np.ndarray) and len(self.lines.shape) == 2 and self.lines.shape[1] == 4
    
    def map_points(self, fun):
        """See `GeometricObject`

        """
        new_points = np.empty_like(self.points)
        for i in range(len(self.points)):
            new_points[i] = fun(self.points[i])
        assert new_points.shape == self.points.shape and new_points.dtype == self.points.dtype
        
        return Mesh(new_points, self.lines, self.triangles, self.physical_to_lines, self.physical_to_triangles)
    
    def _remove_degenerate_triangles(self):
        areas = triangle_areas(self.points[self.triangles[:,:3]])
        degenerate = areas < 1e-12
        map_index = np.arange(len(self.triangles)) - np.cumsum(degenerate)
         
        self.triangles = self.triangles[~degenerate]
        
        for k in self.physical_to_triangles.keys():
            v = self.physical_to_triangles[k]
            self.physical_to_triangles[k] = map_index[v[~degenerate[v]]]
         
        if np.any(degenerate):
            log_debug(f'Removed {sum(degenerate)} degenerate triangles')

    def _deduplicate_points(self):
        if not len(self.points):
            return
         
        # Step 1: Make a copy of the points array using np.array
        points_copy = np.array(self.points, dtype=np.float64)

        # Step 2: Zero the low 16 bits of the mantissa of the X, Y, Z coordinates
        points_copy.view(np.uint64)[:] &= np.uint64(0xFFFFFFFFFFFF0000)

        # Step 3: Use Numpy lexsort directly on points_copy
        sorted_indices = np.lexsort(points_copy.T)
        points_sorted = points_copy[sorted_indices]

        # Step 4: Create a mask to identify unique points
        equal_to_previous = np.all(points_sorted[1:] == points_sorted[:-1], axis=1)
        keep_mask = np.concatenate(([True], ~equal_to_previous))

        # Step 5: Compute new indices for the unique points
        new_indices_in_sorted_order = np.cumsum(keep_mask) - 1

        # Map old indices to new indices
        old_to_new_indices = np.empty(len(points_copy), dtype=np.uint64)
        old_to_new_indices[sorted_indices] = new_indices_in_sorted_order
        
        # Step 6: Update the points array with unique points
        self.points = points_sorted[keep_mask]

        # Step 7: Update all indices
        if len(self.triangles):
            self.triangles = old_to_new_indices[self.triangles]
        if len(self.lines):
            self.lines = old_to_new_indices[self.lines]
    
    @staticmethod
    def _merge_dicts(dict1, dict2):
        dict_: dict[str, np.ndarray] = {}
        
        for (k, v) in chain(dict1.items(), dict2.items()):
            if k in dict_:
                dict_[k] = np.concatenate( (dict_[k], v), axis=0)
            else:
                dict_[k] = v

        return dict_
     
    def __add__(self, other):
        """Add meshes together, using the + operator (mesh1 + mesh2).
        
        Returns
        ------------------------------
        Mesh

        A new mesh consisting of the elements of the added meshes"""
        if not isinstance(other, Mesh):
            return NotImplemented
         
        N_points = len(self.points)
        N_lines = len(self.lines)
        N_triangles = len(self.triangles)
         
        points = _concat_arrays(self.points, other.points)
        lines = _concat_arrays(self.lines, other.lines+N_points)
        triangles = _concat_arrays(self.triangles, other.triangles+N_points)

        other_physical_to_lines = {k:(v+N_lines) for k, v in other.physical_to_lines.items()}
        other_physical_to_triangles = {k:(v+N_triangles) for k, v in other.physical_to_triangles.items()}
         
        physical_lines = Mesh._merge_dicts(self.physical_to_lines, other_physical_to_lines)
        physical_triangles = Mesh._merge_dicts(self.physical_to_triangles, other_physical_to_triangles)
         
        return Mesh(points=points,
                    lines=lines,
                    triangles=triangles,
                    physical_to_lines=physical_lines,
                    physical_to_triangles=physical_triangles)
     
    def extract_physical_group(self, name):
        """Extract a named group from the mesh.

        Parameters
        ---------------------------
        name: str
            Name of the group of elements

        Returns
        --------------------------
        Mesh

        Subset of the mesh consisting only of the elements with the given name."""
        assert name in self.physical_to_lines or name in self.physical_to_triangles, "Physical group not in mesh, so cannot extract"

        if name in self.physical_to_lines:
            elements = self.lines
            physical = self.physical_to_lines
        elif name in self.physical_to_triangles:
            elements = self.triangles
            physical = self.physical_to_triangles
         
        elements_indices = np.unique(physical[name])
        elements = elements[elements_indices]
          
        points_mask = np.full(len(self.points), False)
        points_mask[elements] = True
        points = self.points[points_mask]
          
        new_index = np.cumsum(points_mask) - 1
        elements = new_index[elements]
        physical_to_elements = {name:np.arange(len(elements))}
         
        if name in self.physical_to_lines:
            return Mesh(points=points, lines=elements, physical_to_lines=physical_to_elements)
        elif name in self.physical_to_triangles:
            return Mesh(points=points, triangles=elements, physical_to_triangles=physical_to_elements)
     
    @staticmethod
    def read_file(filename,  name=None):
        """Create a mesh from a given file. All formats supported by meshio are accepted.

        Parameters
        ------------------------------
        filename: str
            Path of the file to convert to Mesh
        name: str
            (optional) name to assign to all elements readed

        Returns
        -------------------------------
        Mesh"""
        meshio_obj = meshio.read(filename)
        mesh = Mesh.from_meshio(meshio_obj)
         
        if name is not None:
            mesh.physical_to_lines[name] = np.arange(len(mesh.lines))
            mesh.physical_to_triangles[name] = np.arange(len(mesh.triangles))
         
        return mesh
     
    def write_file(self, filename):
        """Write a mesh to a given file. The format is determined from the file extension.
        All formats supported by meshio are supported.

        Parameters
        ----------------------------------
        filename: str
            The name of the file to write the mesh to."""
        meshio_obj = self.to_meshio()
        meshio_obj.write(filename)
    
    def write(self, filename):
        self.write_file(filename)
        
     
    def to_meshio(self):
        """Convert the Mesh to a meshio object.

        Returns
        ------------------------------------
        meshio.Mesh"""
        to_write = []
        
        if len(self.lines):
            line_type = 'line' if self.lines.shape[1] == 2 else 'line4'
            to_write.append( (line_type, self.lines) )
        
        if len(self.triangles):
            triangle_type = 'triangle' if self.triangles.shape[1] == 3 else 'triangle6'
            to_write.append( (triangle_type, self.triangles) )
        
        return meshio.Mesh(self.points, to_write)
     
    @staticmethod
    def from_meshio(mesh: meshio.Mesh):
        """Create a Traceon mesh from a meshio.Mesh object.

        Parameters
        --------------------------
        mesh: meshio.Mesh
            The mesh to convert to a Traceon mesh

        Returns
        -------------------------
        Mesh"""
        def extract(type_):
            elements = mesh.cells_dict[type_]
            physical = {k:v[type_] for k,v in mesh.cell_sets_dict.items() if type_ in v}
            return elements, physical
        
        lines, physical_lines = [], {}
        triangles, physical_triangles = [], {}
        
        if 'line' in mesh.cells_dict:
            lines, physical_lines = extract('line')
        elif 'line4' in mesh.cells_dict:
            lines, physical_lines = extract('line4')
        
        if 'triangle' in mesh.cells_dict:
            triangles, physical_triangles = extract('triangle')
        elif 'triangle6' in mesh.cells_dict:
            triangles, physical_triangles = extract('triangle6')
        
        return Mesh(points=mesh.points,
            lines=lines, physical_to_lines=physical_lines,
            triangles=triangles, physical_to_triangles=physical_triangles)
     
    def is_3d(self):
        """Check if the mesh is three dimensional by checking whether any z coordinate is non-zero.

        Returns
        ----------------
        bool

        Whether the mesh is three dimensional"""
        return np.any(self.points[:, 1] != 0.)
    
    def is_2d(self):
        """Check if the mesh is two dimensional, by checking that all z coordinates are zero.
        
        Returns
        ----------------
        bool

        Whether the mesh is two dimensional"""
        return np.all(self.points[:, 1] == 0.)
    
    def flip_normals(self):
        """Flip the normals in the mesh by inverting the 'orientation' of the elements.

        Returns
        ----------------------------
        Mesh"""
        lines = self.lines
        triangles = self.triangles
        
        # Flip the orientation of the lines
        if lines.shape[1] == 4:
            p0, p1, p2, p3 = lines.T
            lines = np.array([p1, p0, p3, p2]).T
        else:
            p0, p1 = lines.T
            lines = np.array([p1, p0]).T
          
        # Flip the orientation of the triangles
        if triangles.shape[1] == 6:
            p0, p1, p2, p3, p4, p5 = triangles.T
            triangles = np.array([p0, p2, p1, p5, p4, p3]).T
        else:
            p0, p1, p2 = triangles.T
            triangles = np.array([p0, p2, p1]).T
        
        return Mesh(self.points, lines, triangles,
            physical_to_lines=self.physical_to_lines,
            physical_to_triangles=self.physical_to_triangles)
     
    def remove_lines(self):
        """Remove all the lines from the mesh.

        Returns
        -----------------------------
        Mesh"""
        return Mesh(self.points, triangles=self.triangles, physical_to_triangles=self.physical_to_triangles)
    
    def remove_triangles(self):
        """Remove all triangles from the mesh.

        Returns
        -------------------------------------
        Mesh"""
        return Mesh(self.points, lines=self.lines, physical_to_lines=self.physical_to_lines)
     
    def get_electrodes(self):
        """Get the names of all the named groups (i.e. electrodes) in the mesh
         
        Returns
        ---------
        str iterable

        Names
        """
        return list(self.physical_to_lines.keys()) + list(self.physical_to_triangles.keys())
     
    @staticmethod
    def _lines_to_higher_order(points, elements):
        N_elements = len(elements)
        N_points = len(points)
         
        v0, v1 = elements.T
        p2 = points[v0] + (points[v1] - points[v0]) * 1/3
        p3 = points[v0] + (points[v1] - points[v0]) * 2/3
         
        assert all(p.shape == (N_elements, points.shape[1]) for p in [p2, p3])
         
        points = np.concatenate( (points, p2, p3), axis=0)
          
        elements = np.array([
            elements[:, 0], elements[:, 1], 
            np.arange(N_points, N_points + N_elements, dtype=np.uint64),
            np.arange(N_points + N_elements, N_points + 2*N_elements, dtype=np.uint64)]).T
         
        assert np.allclose(p2, points[elements[:, 2]]) and np.allclose(p3, points[elements[:, 3]])
        return points, elements


    def _to_higher_order_mesh(self):
        # The matrix solver currently only works with higher order meshes.
        # We can however convert a simple mesh easily to a higher order mesh, and solve that.
        
        points, lines, triangles = self.points, self.lines, self.triangles

        if not len(lines):
            lines = np.empty( (0, 4), dtype=np.float64)
        elif len(lines) and lines.shape[1] == 2:
            points, lines = Mesh._lines_to_higher_order(points, lines)
        
        assert lines.shape == (len(lines), 4)

        return Mesh(points=points,
            lines=lines, physical_to_lines=self.physical_to_lines,
            triangles=triangles, physical_to_triangles=self.physical_to_triangles)
     
    def __str__(self):
        physical_lines = ', '.join(self.physical_to_lines.keys())
        physical_lines_nums = ', '.join([str(len(self.physical_to_lines[n])) for n in self.physical_to_lines.keys()])
        physical_triangles = ', '.join(self.physical_to_triangles.keys())
        physical_triangles_nums = ', '.join([str(len(self.physical_to_triangles[n])) for n in self.physical_to_triangles.keys()])
        
        return f'<Traceon Mesh,\n' \
            f'\tNumber of points: {len(self.points)}\n' \
            f'\tNumber of lines: {len(self.lines)}\n' \
            f'\tNumber of triangles: {len(self.triangles)}\n' \
            f'\tPhysical lines: {physical_lines}\n' \
            f'\tElements in physical line groups: {physical_lines_nums}\n' \
            f'\tPhysical triangles: {physical_triangles}\n' \
            f'\tElements in physical triangle groups: {physical_triangles_nums}>'

    def _ensure_normal_orientation_triangles(self, electrode, outwards):
        assert electrode in self.physical_to_triangles, "electrode should be part of mesh"
        
        triangle_indices = self.physical_to_triangles[electrode]
        electrode_triangles = self.triangles[triangle_indices]
          
        if not len(electrode_triangles):
            return
        
        connected_indices = _get_connected_elements(electrode_triangles)
        
        for indices in connected_indices:
            connected_triangles = electrode_triangles[indices]
            _ensure_triangle_orientation(connected_triangles, self.points, outwards)
            electrode_triangles[indices] = connected_triangles

        self.triangles[triangle_indices] = electrode_triangles
     
    def _ensure_normal_orientation_lines(self, electrode, outwards):
        assert electrode in self.physical_to_lines, "electrode should be part of mesh"
        
        line_indices = self.physical_to_lines[electrode]
        electrode_lines = self.lines[line_indices]
          
        if not len(electrode_lines):
            return
        
        connected_indices = _get_connected_elements(electrode_lines)
        
        for indices in connected_indices:
            connected_lines = electrode_lines[indices]
            _ensure_line_orientation(connected_lines, self.points, outwards)
            electrode_lines[indices] = connected_lines

        self.lines[line_indices] = electrode_lines
     
    def ensure_outward_normals(self, electrode):
        if electrode in self.physical_to_triangles:
            self._ensure_normal_orientation_triangles(electrode, True)
        
        if electrode in self.physical_to_lines:
            self._ensure_normal_orientation_lines(electrode, True)
     
    def ensure_inward_normals(self, electrode):
        if electrode in self.physical_to_triangles:
            self._ensure_normal_orientation_triangles(electrode, False)
         
        if electrode in self.physical_to_lines:
            self._ensure_normal_orientation_lines(electrode, False)

Ancestors

• traceon.util.Saveable
• GeometricObject
• abc.ABC

Static methods

def from_meshio(mesh: meshio._mesh.Mesh)

Create a Traceon mesh from a meshio.Mesh object.

Parameters

mesh : meshio.Mesh

The mesh to convert to a Traceon mesh

Returns

Mesh

 

def read_file(filename, name=None)

Create a mesh from a given file. All formats supported by meshio are accepted.

Parameters

filename : str

Path of the file to convert to Mesh

name : str

(optional) name to assign to all elements readed

Returns

Mesh

 

Methods

def __add__(self, other)

Add meshes together, using the + operator (mesh1 + mesh2).

Returns

Mesh

 

A new mesh consisting of the elements of the added meshes

 

def ensure_inward_normals(self, electrode)

def ensure_outward_normals(self, electrode)

def extract_physical_group(self, name)

Extract a named group from the mesh.

Parameters

name : str

Name of the group of elements

Returns

Mesh

 

Subset of the mesh consisting only of the elements with the given name.

def flip_normals(self)

Flip the normals in the mesh by inverting the 'orientation' of the elements.

Returns

Mesh

 

def get_electrodes(self)

Get the names of all the named groups (i.e. electrodes) in the mesh

Returns

str iterable

 

Names

 

def is_2d(self)

Check if the mesh is two dimensional, by checking that all z coordinates are zero.

Returns

bool

 

Whether the mesh is two dimensional

 

def is_3d(self)

Check if the mesh is three dimensional by checking whether any z coordinate is non-zero.

Returns

bool

 

Whether the mesh is three dimensional

 

def is_higher_order(self)

Whether the mesh contains higher order elements.

Returns

bool

 

def map_points(self, fun)

See GeometricObject

def remove_lines(self)

Remove all the lines from the mesh.

Returns

Mesh

 

def remove_triangles(self)

Remove all triangles from the mesh.

Returns

Mesh

 

def to_meshio(self)

Convert the Mesh to a meshio object.

Returns

meshio.Mesh

 

def write(self, filename)

Write a mesh to a file. The pickle module will be used to save the Geometry object.

Args

filename

name of the file

def write_file(self, filename)

Write a mesh to a given file. The format is determined from the file extension. All formats supported by meshio are supported.

Parameters

filename : str

The name of the file to write the mesh to.

Inherited members

• GeometricObject:
  • mirror_xy
  • mirror_xz
  • mirror_yz
  • move
  • rotate