Module traceon.mesher

Classes

class GeometricObject

Expand source code

class GeometricObject(ABC):
    """The Mesh class (and the classes defined in `traceon.geometry`) are subclasses
    of `traceon.mesher.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 counterclockwise 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"

        if Rx !=0.:
            return self.rotate_around_axis([1,0,0], angle=Rx, origin=origin)
        if Ry !=0.:
            return self.rotate_around_axis([0,1,0], angle=Ry, origin=origin)
        if Rz !=0.:
            return self.rotate_around_axis([0,0,1], angle=Rz, origin=origin)
        

    
    def rotate_around_axis(self, axis=[0., 0, 1.], angle=0., origin=[0., 0., 0.]):
        """
        Rotate  counterclockwise around a general axis defined by a vector.
        
        Parameters
        ------------------------------------
        axis: (3,) float
            Vector defining the axis of rotation. Must be non-zero.
        angle: float
            Amount to rotate around the 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.
        """
        origin = np.array(origin, dtype=float)
        assert origin.shape == (3,), "Please supply a 3D point for origin"
        axis = np.array(axis, dtype=float)
        assert axis.shape == (3,), "Please supply a 3D vector for axis"
        axis = axis / np.linalg.norm(axis)

        if angle == 0:
            return self.map_points(lambda x: x)
        
        # Rotation is implemented using Rodrigues' rotation formula
        # See: https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
        K = np.array([
            [0, -axis[2], axis[1]],
            [axis[2], 0, -axis[0]],
            [-axis[1], axis[0], 0]])
        
        K2 = K @ K # Transformation matrix that maps a vector to its cross product with the rotation axis
        I = np.eye(3)
        R = I + np.sin(angle) * K + (1 - np.cos(angle)) * K2

        return self.map_points(lambda p: origin + R @ (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]]))

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

Ancestors

• abc.ABC

Subclasses

• Path
• PathCollection
• Surface
• SurfaceCollection
• Mesh

Methods

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

Expand source code

@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."""
    ...

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)

Expand source code

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]]))

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)

Expand source code

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]]))

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)

Expand source code

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]]))

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)

Expand source code

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]))

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])

Expand source code

def rotate(self, Rx=0., Ry=0., Rz=0., origin=[0., 0., 0.]):
    """Rotate counterclockwise 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"

    if Rx !=0.:
        return self.rotate_around_axis([1,0,0], angle=Rx, origin=origin)
    if Ry !=0.:
        return self.rotate_around_axis([0,1,0], angle=Ry, origin=origin)
    if Rz !=0.:
        return self.rotate_around_axis([0,0,1], angle=Rz, origin=origin)

Rotate counterclockwise 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.

def rotate_around_axis(self, axis=[0.0, 0, 1.0], angle=0.0, origin=[0.0, 0.0, 0.0])

Expand source code

def rotate_around_axis(self, axis=[0., 0, 1.], angle=0., origin=[0., 0., 0.]):
    """
    Rotate  counterclockwise around a general axis defined by a vector.
    
    Parameters
    ------------------------------------
    axis: (3,) float
        Vector defining the axis of rotation. Must be non-zero.
    angle: float
        Amount to rotate around the 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.
    """
    origin = np.array(origin, dtype=float)
    assert origin.shape == (3,), "Please supply a 3D point for origin"
    axis = np.array(axis, dtype=float)
    assert axis.shape == (3,), "Please supply a 3D vector for axis"
    axis = axis / np.linalg.norm(axis)

    if angle == 0:
        return self.map_points(lambda x: x)
    
    # Rotation is implemented using Rodrigues' rotation formula
    # See: https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
    K = np.array([
        [0, -axis[2], axis[1]],
        [axis[2], 0, -axis[0]],
        [-axis[1], axis[0], 0]])
    
    K2 = K @ K # Transformation matrix that maps a vector to its cross product with the rotation axis
    I = np.eye(3)
    R = I + np.sin(angle) * K + (1 - np.cos(angle)) * K2

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

Rotate counterclockwise around a general axis defined by a vector.

Parameters

axis : (3,) float

Vector defining the axis of rotation. Must be non-zero.

angle : float

Amount to rotate around the 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)

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)

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.

Ancestors

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

Static methods

def from_meshio(mesh: meshio._mesh.Mesh)

Expand source code

@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)

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)

Expand source code

@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

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)

Expand source code

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)

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)

Expand source code

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)

def ensure_outward_normals(self, electrode)

Expand source code

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 extract_physical_group(self, name)

Expand source code

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)

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)

Expand source code

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)

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

Returns

Mesh

 

def get_electrodes(self)

Expand source code

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())

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

Returns

str iterable

 

Names

 

def is_2d(self)

Expand source code

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.)

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)

Expand source code

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.)

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)

Expand source code

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

Whether the mesh contains higher order elements.

Returns

bool

 

def map_points(self, fun)

Expand source code

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)

See GeometricObject

def remove_lines(self)

Expand source code

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)

Remove all the lines from the mesh.

Returns

Mesh

 

def remove_triangles(self)

Expand source code

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)

Remove all triangles from the mesh.

Returns

Mesh

 

def to_meshio(self)

Expand source code

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)

Convert the Mesh to a meshio object.

Returns

meshio.Mesh

 

def write(self, filename)

Expand source code

def write(self, filename):
    self.write_file(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)

Expand source code

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)

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
  • rotate_around_axis