geometry.py

"""Classes describing geometry in sectionproperties."""

from __future__ import annotations

import copy
import math
import pathlib
from typing import TYPE_CHECKING, Any, cast

import more_itertools
import numpy as np
import numpy.typing as npt
from shapely import (
    GeometryCollection,
    LinearRing,
    LineString,
    MultiLineString,
    MultiPolygon,
    Point,
    Polygon,
    affinity,
    box,  # pyright: ignore [reportUnknownVariableType]
    line_merge,  # pyright: ignore [reportUnknownVariableType]
    shared_paths,  # pyright: ignore [reportUnknownVariableType]
)
from shapely.ops import split, unary_union

import sectionproperties.post.post as post
import sectionproperties.pre.bisect_section as bisect
import sectionproperties.pre.pre as pre

if TYPE_CHECKING:
    import matplotlib.axes


SCALE_CONSTANT = 1e-9


class Geometry:
    """Class for defining the geometry of a contiguous section of a single material.

    Provides an interface for the user to specify the geometry defining a section. A
    method is provided for generating a triangular mesh, transforming the section (e.g.
    translation, rotation, perimeter offset, mirroring), aligning the geometry to
    another geometry, and designating stress recovery points.
    """

    def __init__(
        self,
        geom: Polygon,
        material: pre.Material | None = pre.DEFAULT_MATERIAL,
        control_points: Point | tuple[float, float] | None = None,
        tol: int = 12,
    ) -> None:
        """Inits the Geometry class.

        Args:
            geom: A Shapely Polygon object that defines the geometry
            material: A material to associate with this geometry. Defaults to
                ``pre.DEFAULT_MATERIAL``.
            control_points: An ``(x, y)`` coordinate within the geometry that represents
                a pre-assigned control point (aka, a region identification point) to be
                used instead of the automatically assigned control point generated with
                :meth:`shapely.Polygon.representative_point`. Defaults to ``None``.
            tol: Number of decimal places to round the geometry vertices to. A lower
                value may reduce accuracy of geometry but increases precision when
                aligning geometries to each other. Defaults to ``12``.

        Raises:
            ValueError: If ``geom`` is not valid, i.e. not a shapely object, or a
                MultiPolygon object
        """
        if isinstance(geom, MultiPolygon):
            msg = "Use CompoundGeometry(...) for a MultiPolygon object."
            raise ValueError(msg)

        if not isinstance(geom, Polygon):  # pyright: ignore [reportUnnecessaryIsInstance]
            msg = f"Argument is not a valid shapely.Polygon object: {geom!r}"
            raise ValueError(msg)

        self.assigned_control_point = None

        if control_points is not None and (
            isinstance(control_points, Point) or len(control_points) == 2
        ):
            self.assigned_control_point = Point(control_points)

        self.tol = tol  # num of decimal places of precision for point locations
        self.geom = round_polygon_vertices(geom, self.tol)
        self.material = pre.DEFAULT_MATERIAL if material is None else material
        self.control_points: list[tuple[float, float]] = []
        self.points: list[tuple[float, float]] = []
        self.facets: list[tuple[int, int]] = []
        self.holes: list[tuple[float, float]] = []
        self._recovery_points: list[tuple[float, float]] | list[Point] = []
        self.mesh: dict[str, Any] | None = None

        self.compile_geometry()

    def _repr_svg_(self) -> str:
        """Generates an svg of the geometry.

        Returns:
            Representative svg
        """
        print("sectionproperties.pre.geometry.Geometry")
        print(f"object at: {hex(id(self))}")
        print(f"Material: {self.material.name}")
        return str(self.geom._repr_svg_())

    def assign_control_point(
        self,
        control_point: tuple[float, float],
    ) -> Geometry:
        """Assigns a control point to the geometry.

        Returns a new Geometry object with ``control_point`` assigned as the control
        point for the new Geometry. The assignment of a control point is intended to
        replace the control point automatically generated by
        :meth:`shapely.Polygon.representative_point()`.

        An assigned control point is carried through and transformed with the Geometry
        whenever it is shifted, aligned, mirrored, unioned, and/or rotated. If a
        perimeter_offset operation is applied, a check is performed to see if the
        assigned control point is still valid (within the new region) and, if so, it is
        kept. If not, a new control point is auto-generated.

        The same check is performed when the geometry undergoes a difference operation
        (with the ``-`` operator) or a shift_points operation. If the assigned control
        point is valid, it is kept. If not, a new one is auto-generated.

        For all other operations (e.g. symmetric difference, intersection, split), the
        assigned control point is discarded and a new one auto-generated.

        Args:
            control_point: An (``x``, ``y``) coordinate that describes the distinct,
                contiguous, region of a single material within the geometry.

        Returns:
            New Geometry object with ``control_point`` assigned as the control point
        """
        return Geometry(
            geom=self.geom,
            material=self.material,
            control_points=control_point,
            tol=self.tol,
        )

    @staticmethod
    def from_points(
        points: list[tuple[float, float]],
        facets: list[tuple[int, int]],
        control_points: list[tuple[float, float]],
        holes: list[tuple[float, float]] | None = None,
        material: pre.Material = pre.DEFAULT_MATERIAL,
    ) -> Geometry:
        """Creates Geometry from points, facets, a control point and holes.

        Args:
            points: List of points (``x``, ``y``) defining the vertices of the section
                geometry. If the geometry simply contains a continuous list of exterior
                points, consider creating a :class:`shapely.Polygon` object (only
                requiring points), and create a ``Geometry`` object using the
                constructor.
            facets: A list of (``start``, ``end``) indices of vertices defining the
                edges of the section geoemtry. Can be used to define both external and
                internal perimeters of holes. Facets are assumed to be described in the
                order of exterior perimeter, interior perimeter 1, interior perimeter 2,
                etc.
            control_points: An (``x``, ``y``) coordinate that describes the distinct,
                contiguous, region of a single material within the geometry. Must be
                entered as a list of coordinates, e.g. [(0.5, 3.2)]. Exactly one point
                is required for each geometry with a distinct material. If there are
                multiple distinct regions, then use
                :meth:`sectionproperties.pre.geometry.CompoundGeometry.from_points()`
            holes: A list of points (``x``, ``y``) that define interior regions as
                being holes or voids. The point can be located anywhere within the hole
                region. Only one point is required per hole region. Defaults to
                ``None``.
            material: A :class:`~sectionproperties.pre.pre.Material` object that is to
                be assigned. Defaults to ``pre.DEFAULT_MATERIAL``.

        Raises:
            ValueError: If there is not exactly one control point specified
            RuntimeError: If geometry generation failed

        Returns:
            Geometry object

        Example:
            .. plot::
                :include-source: True
                :caption: Geometry created from points, facets and holes.

                from sectionproperties.pre import Geometry

                points = [(0, 0), (10, 5), (15, 15), (5, 10), (6, 6), (9, 7), (7, 9)]
                facets = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 4)]
                control_points = [(4, 4)]
                holes = [(7, 7)]

                Geometry.from_points(
                    points=points,
                    facets=facets,
                    control_points=control_points,
                    holes=holes,
                ).plot_geometry()
        """
        if len(control_points) != 1:
            msg = "Control points for Geometry instances must have exactly "
            msg += "one x, y coordinate and entered as a list of 2D float tuples, "
            msg += "e.g. [(0.1, 3.4)]. CompoundGeometry.from_points() can accept "
            msg += "multiple control points\n"
            msg += f"Control points received: {control_points}"
            raise ValueError(msg)

        if holes is None:
            holes = []

        prev_facet: tuple[int, int] | None = None

        # Initialize the total number of accumulators needed
        # Always an exterior, plus, a separate accumulator for each interior region
        exterior: list[tuple[float, float]] = []
        interiors: list[list[tuple[float, float]]] = [[] for _ in holes]
        interior_counter = 0  # To keep track of interior regions
        active_list = exterior  # The active_list is the list being accumulated on

        for facet in facets:  # Loop through facets for graph connectivity
            i_idx, _ = facet

            if not prev_facet:  # Add the first facet vertex to exterior and move on
                vertex = points[i_idx]
                active_list.append(vertex)
                prev_facet = facet
                continue

            # Look at the last j_idx to test for a break in the chain of edges
            prev_j_idx = prev_facet[1]

            # If there is a break in the chain of edges...
            if i_idx != prev_j_idx and holes:
                # ...and we are still accumulating on the exterior...
                if active_list == exterior:
                    active_list = interiors[interior_counter]
                    # ... then move to the interior accumulator
                # ...or if we are already in the interior accumulator...
                else:
                    # ...then start the next interior accumulator for a new hole.
                    interior_counter += 1
                    active_list = interiors[interior_counter]
                active_list.append(points[i_idx])
            else:
                # Only need i_idx b/c shapely auto-closes polygons
                active_list.append(points[i_idx])

            prev_facet = facet

        exterior_geometry = Polygon(exterior)
        interior_polys = [Polygon(interior) for interior in interiors]
        interior_geometry = MultiPolygon(interior_polys)
        sub_poly = exterior_geometry - interior_geometry

        if isinstance(sub_poly, Polygon):
            return Geometry(
                geom=sub_poly, material=material, control_points=control_points[0]
            )
        else:
            msg = "Geometry generation failed."
            raise RuntimeError(msg)

    @staticmethod
    def from_dxf(
        dxf_filepath: str | pathlib.Path,
        spline_delta: float = 0.1,
        degrees_per_segment: float = 1.0,
    ) -> Geometry | CompoundGeometry:
        """An interface for the creation of Geometry objects from CAD .dxf files.

        Args:
            dxf_filepath: A path-like object for the dxf file
            spline_delta: Splines are not supported in ``shapely``, so they are
                approximated as polylines, this argument affects the spline sampling
                rate. Defaults to ``0.1``.
            degrees_per_segment: The number of degrees discretised as a single line
                segment. Defaults to ``1.0``.

        Returns:
            Geometry or CompoundGeometry object
        """
        return load_dxf(
            dxf_filepath=dxf_filepath,
            spline_delta=spline_delta,
            degrees_per_segment=degrees_per_segment,
        )

    @classmethod
    def from_3dm(
        cls,
        filepath: str | pathlib.Path,
        **kwargs: Any,
    ) -> Geometry:
        """Creates a Geometry object from a Rhino ``.3dm`` file.

        Args:
            filepath: File path to the rhino ``.3dm`` file.
            kwargs: See below.

        Keyword Args:
            refine_num (Optional[int]):  Bézier curve interpolation number. In Rhino a
                surface's edges are nurb based curves. Shapely does not support nurbs,
                so the individual Bézier curves are interpolated using straight lines.
                This parameter sets the number of straight lines used in the
                interpolation. Default is 1.
            vec1 (Optional[numpy.ndarray]): A 3d vector in the Shapely plane. Rhino is a
                3D geometry environment. Shapely is a 2D geometric library. Thus a 2D
                plane needs to be defined in Rhino that represents the Shapely
                coordinate system. ``vec1`` represents the 1st vector of this plane. It
                will be used as Shapely's x direction. Default is [1,0,0].
            vec2 (Optional[numpy.ndarray]): Continuing from ``vec1``, ``vec2`` is
                another vector to define the Shapely plane. It must not be [0,0,0] and
                it's only requirement is that it is any vector in the Shapely plane (but
                not equal to ``vec1``). Default is [0,1,0].
            plane_distance (Optional[float]): The distance to the Shapely plane. Default
                is 0.
            project (Optional[bool]): Controls if the breps are projected onto the plane
                in the direction of the Shapley plane's normal. Default is True.
            parallel (Optional[bool]): Controls if only the rhino surfaces that have the
                same normal as the Shapely plane are yielded. If true, all non parallel
                surfaces are filtered out. Default is False.

        Raises:
            RuntimeError: A RuntimeError is raised if two or more polygons are found.
                This is dependent on the keyword arguments. Try adjusting the keyword
                arguments if this error is raised.
            ImportError: If ``rhino3dm`` is not installed. To enable rhino features use
                ``pip install sectionproperties[rhino]``.

        Returns:
            A Geometry object.
        """
        try:
            import sectionproperties.pre.rhino as rhino_importer
        except ImportError as e:
            print(e)
            msg = "There is something wrong with your rhino library installation. "
            msg += "Please report this error at "
            msg += "https://github.com/robbievanleeuwen/section-properties/issues"
            raise ImportError(msg) from e

        geom = None
        list_poly = rhino_importer.load_3dm(filepath, **kwargs)

        if len(list_poly) == 1:
            geom = cls(geom=list_poly[0])
        else:
            msg = "Multiple surfaces extracted from the file. Either use "
            msg += "CompoundGeometry or extract individual surfaces manually via "
            msg += "pre.rhino."
            raise RuntimeError(msg)

        return geom

    @classmethod
    def from_rhino_encoding(
        cls,
        r3dm_brep: str,
        **kwargs: Any,
    ) -> Geometry:
        """Load an encoded single surface planer brep.

        Args:
            r3dm_brep: A Rhino3dm.Brep encoded as a string.
            kwargs: See below.

        Keyword Args:
            refine_num (Optional[int]):  Bézier curve interpolation number. In Rhino a
                surface's edges are nurb based curves. Shapely does not support nurbs,
                so the individual Bézier curves are interpolated using straight lines.
                This parameter sets the number of straight lines used in the
                interpolation. Default is 1.
            vec1 (Optional[numpy.ndarray]): A 3d vector in the Shapely plane. Rhino is a
                3D geometry environment. Shapely is a 2D geometric library. Thus a 2D
                plane needs to be defined in Rhino that represents the Shapely
                coordinate system. ``vec1`` represents the 1st vector of this plane. It
                will be used as Shapely's x direction. Default is [1,0,0].
            vec2 (Optional[numpy.ndarray]): Continuing from ``vec1``, ``vec2`` is
                another vector to define the Shapely plane. It must not be [0,0,0] and
                it's only requirement is that it is any vector in the Shapely plane (but
                not equal to ``vec1``). Default is [0,1,0].
            plane_distance (Optional[float]): The distance to the Shapely plane. Default
                is 0.
            project (Optional[bool]): Controls if the breps are projected onto the plane
                in the direction of the Shapley plane's normal. Default is True.
            parallel (Optional[bool]): Controls if only the rhino surfaces that have the
                same normal as the Shapely plane are yielded. If true, all non parallel
                surfaces are filtered out. Default is False.

        Raises:
            ImportError: If ``rhino3dm`` is not installed. To enable rhino features use
                ``pip install sectionproperties[rhino]``.

        Returns:
            A Geometry object found in the encoded string.
        """
        try:
            import sectionproperties.pre.rhino as rhino_importer
        except ImportError as e:
            msg = "There is something wrong with your rhino library installation. "
            msg += "Please report this error at "
            msg += "https://github.com/robbievanleeuwen/section-properties/issues"
            raise ImportError(msg) from e

        return cls(geom=rhino_importer.load_brep_encoding(r3dm_brep, **kwargs)[0])

    def create_facets_and_control_points(self) -> None:
        """Creates geometry data from shapely data.

        Generates points, facets, control points, holes and perimeter from the
        shapely geometry.
        """
        self.holes = []
        self.points = []
        self.facets = []
        self.points, self.facets = create_points_and_facets(
            shape=self.geom,
            tol=self.tol,
        )

        if not self.assigned_control_point:
            self.control_points = list(self.geom.representative_point().coords)  # pyright: ignore [reportAttributeAccessIssue]
        else:
            self.control_points = list(self.assigned_control_point.coords)  # pyright: ignore [reportAttributeAccessIssue]

        for hole in self.geom.interiors:
            hole_polygon = Polygon(hole)
            self.holes += tuple(hole_polygon.representative_point().coords)  # pyright: ignore [reportOperatorIssue]

    def compile_geometry(self) -> None:
        """Alias for ``create_facets_and_control_points()``.

        Alters attributes .points, .facets, .holes, .control_points to represent
        the data in the shapely geometry.
        """
        self.create_facets_and_control_points()

    def create_mesh(
        self,
        mesh_sizes: float | list[float],
        min_angle: float = 30.0,
        coarse: bool = False,
    ) -> Geometry:
        r"""Creates a quadratic triangular mesh from the Geometry object.

        Args:
            mesh_sizes: A float describing the maximum mesh element area to be used
                within the Geometry-object finite-element mesh (may also be a list of
                length 1)
            min_angle: The meshing algorithm adds vertices to the mesh to ensure that no
                angle smaller than the minimum angle (in degrees, rounded to 1 decimal
                place). Note that small angles between input segments cannot be
                eliminated. If the minimum angle is 20.7 deg or smaller, the
                triangulation algorithm is theoretically guaranteed to terminate (given
                sufficient precision). The algorithm often doesn't terminate for angles
                greater than 33 deg. Some meshes may require angles well below 20 deg to
                avoid problems associated with insufficient floating-point precision.
                Defaults to ``30.0``.
            coarse: If set to True, will create a coarse mesh (no area or quality
                constraints). Defaults to ``False``.

        Raises:
            ValueError: ``mesh_sizes`` is not valid

        Returns:
            ``Geometry`` object with mesh data stored in ``.mesh`` attribute. Returned
            ``Geometry`` object is self, not a new instance.

        Example:
            The following example creates a circular cross-section with a diameter of
            50 mm with 64 points, and generates a mesh with a maximum triangular area of
            2.5 mm\ :sup:`2`.

            .. plot::
                :include-source: True
                :caption: Mesh for a 50 mm diameter circle

                from sectionproperties.pre.library import circular_section
                from sectionproperties.analysis import Section

                geom = circular_section(d=50, n=64)
                geom = geom.create_mesh(mesh_sizes=2.5)
                Section(geom).plot_mesh(materials=False)
        """
        if isinstance(mesh_sizes, list) and len(mesh_sizes) == 1:
            mesh_size = mesh_sizes[0]
        elif isinstance(mesh_sizes, float | int):
            mesh_size = mesh_sizes
        else:
            msg = "Argument 'mesh_sizes' for a Geometry must be either a float, or a "
            msg += f"list of float with length of 1, not {mesh_sizes}."
            raise ValueError(msg)

        self.mesh = pre.create_mesh(
            points=self.points,
            facets=self.facets,
            holes=self.holes,
            control_points=self.control_points,
            mesh_sizes=mesh_size,
            min_angle=min_angle,
            coarse=coarse,
        )

        return self

    def align_to(
        self,
        other: Geometry | CompoundGeometry | tuple[float, float],
        on: str,
        inner: bool = False,
    ) -> Geometry:
        """Aligns the geometry to another Geometry object.

        Returns a new Geometry object, representing ``self`` translated so that is
        aligned ``on`` one of the outer bounding box edges of ``other``.

        If ``other`` is a tuple representing an (``x``, ``y``) coordinate, then the new
        Geometry object will represent 'self' translated so that it is aligned ``on``
        that side of the point.

        Args:
            other: Either another Geometry or a tuple representing an (``x``, ``y``)
                coordinate point that ``self`` should align to.
            on: A str of either "left", "right", "bottom", or "top" indicating which
                side of ``other`` that ``self`` should be aligned to.
            inner: If True, align ``self`` to ``other`` in such a way that ``self`` is
                aligned to the "inside" of ``other``. In other words, align ``self`` to
                ``other`` on the specified edge so they overlap. Defaults to ``False``.

        Returns:
             Geometry object translated to alignment location

        Example:
            .. plot::
                :include-source: True
                :caption: A triangle aligned to the top of a rectangle.

                from sectionproperties.pre.library import rectangular_section
                from sectionproperties.pre.library import triangular_section

                rect = rectangular_section(b=100, d=50)
                tri = triangular_section(b=50, h=50)
                tri = tri.align_to(other=rect, on="top")
                (rect + tri).plot_geometry()
        """
        # Mappings are for indexes in the list of bbox extents of both
        # 'self' and 'align_to'. i.e. a mapping of which "word" corresponds
        # to which bounding box coordinate
        align_self_map = {
            "left": 1,
            "right": 0,
            "bottom": 3,
            "top": 2,
        }

        other_as_geom_map = {
            "left": 0,
            "right": 1,
            "bottom": 2,
            "top": 3,
        }

        other_as_point_map = {
            "left": 0,
            "right": 0,
            "bottom": 1,
            "top": 1,
        }

        self_align_idx = align_self_map[on]

        if isinstance(other, Geometry | CompoundGeometry):
            align_to_idx = other_as_geom_map[on]
            align_to_coord = other.calculate_extents()[align_to_idx]
        else:
            align_to_idx = other_as_point_map[on]
            align_to_coord = other[align_to_idx]

        self_extents = self.calculate_extents()  # min x, max x, min y, max y
        self_align_coord = self_extents[self_align_idx]

        if inner:
            self_align_coord = self_extents[align_to_idx]

        offset = align_to_coord - self_align_coord
        arg = "x_offset"

        if on in ["top", "bottom"]:
            arg = "y_offset"

        kwargs = {arg: offset}

        return self.shift_section(**kwargs)

    def align_center(
        self,
        align_to: Geometry | tuple[float, float] | None = None,
    ) -> Geometry:
        """Aligns the geometry to a center point.

        Returns a new Geometry object, translated in both ``x`` and ``y``, so that the
        the new object's centroid will be aligned with the centroid of the object in
        ``align_to``. If ``align_to`` is an (``x``, ``y``) coordinate, then the centroid
        will be aligned to the coordinate. If ``align_to`` is ``None`` then the new
        object will be aligned with its centroid at the origin.

        Args:
            align_to: Another Geometry to align to, an (``x``, ``y``) coordinate or
                ``None``. Defaults to ``None``.

        Raises:
            ValueError: ``align_to`` is not valid

        Returns:
            Geometry object translated to new alignment

        Example:
            .. plot::
                :include-source: True
                :caption: A triangular hole aligned to the centre of a square.

                from sectionproperties.pre.library import rectangular_section
                from sectionproperties.pre.library import triangular_section

                rect = rectangular_section(b=200, d=200)
                tri = triangular_section(b=50, h=50)
                tri = tri.align_center(align_to=rect)
                geom = rect + tri
                geom.holes = [(100, 100)]
                geom.control_points = [(25, 25)]
                geom.plot_geometry()
        """
        cx, cy = next(iter(self.geom.centroid.coords))

        # Suggested by @Agent6-6-6: Hard-rounding of cx and cy allows
        # for greater precision in placing geometry with its centroid
        # near [0, 0]. True [0, 0] placement will not be possible due
        # to floating point errors.
        if align_to is None:
            shift_x, shift_y = round(-cx, self.tol), round(-cy, self.tol)
        elif isinstance(align_to, Geometry):
            align_cx, align_cy = next(iter(align_to.geom.centroid.coords))  # RUF015
            shift_x = round(align_cx - cx, self.tol)
            shift_y = round(align_cy - cy, self.tol)
        else:
            try:
                point_x, point_y = align_to
                shift_x = round(point_x - cx, self.tol)
                shift_y = round(point_y - cy, self.tol)
            except (
                TypeError,
                ValueError,
            ) as e:  # align_to not subscriptable, incorrect length, etc.
                msg = "align_to must be either a Geometry object, an x, y coordinate, "
                msg += f"or None, not {align_to}."
                raise ValueError(msg) from e

        return self.shift_section(x_offset=shift_x, y_offset=shift_y)

    def shift_section(
        self,
        x_offset: float = 0.0,
        y_offset: float = 0.0,
    ) -> Geometry:
        """Returns a new Geometry object translated by (``x_offset``, ``y_offset``).

        Args:
            x_offset: Distance in x-direction by which to shift the geometry. Defaults
                to ``0.0``.
            y_offset: Distance in y-direction by which to shift the geometry. Defaults
                to ``0.0``.

        Returns:
            New Geometry object shifted by ``x_offset`` and ``y_offset``

        Example:
            .. plot::
                :include-source: True
                :caption: Using ``shift_section`` to shift geometry.

                from sectionproperties.pre.library import rectangular_section

                rect1 = rectangular_section(b=200, d=100)
                rect2 = rect1.shift_section(x_offset=100, y_offset=100)
                (rect1 + rect2).plot_geometry()
        """
        # Move assigned control point
        new_ctrl_point: tuple[float, float] | None = None

        if self.assigned_control_point:
            new_ctrl_point_geom = affinity.translate(
                self.assigned_control_point, x_offset, y_offset
            )
            new_ctrl_point = new_ctrl_point_geom.x, new_ctrl_point_geom.y

        return Geometry(
            geom=affinity.translate(self.geom, x_offset, y_offset),
            material=self.material,
            control_points=new_ctrl_point,
            tol=self.tol,
        )

    def rotate_section(
        self,
        angle: float,
        rot_point: tuple[float, float] | str = "center",
        use_radians: bool = False,
    ) -> Geometry:
        """Rotate the Geometry object.

        Rotates the geometry and specified angle about a point. If the rotation point is
        not provided, rotates the section about the center of the geometry's bounding
        box.

        Args:
            angle: Angle (degrees by default) by which to rotate the section. A positive
                angle leads to a counter-clockwise rotation.
            rot_point: Point (``x``, ``y``) about which to rotate the section. If not
                provided, will rotate about the "center" of the geometry's bounding box.
                Defaults to ``"center"``.
            use_radians: Boolean to indicate whether ``angle`` is in degrees or radians.
                If True, ``angle`` is interpreted as radians. Defaults to ``False``.

        Returns:
            New Geometry object rotated by ``angle`` about ``rot_point``

        Example:
            .. plot::
                :include-source: True
                :caption: A 200UB25 section rotated clockwise by 30 degrees

                from sectionproperties.pre.library import i_section

                geom = i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
                geom.rotate_section(angle=-30).plot_geometry()
        """
        new_ctrl_point = None

        if self.assigned_control_point:
            if rot_point == "center":
                centre = box(*self.geom.bounds).centroid
                rotate_point = centre.x, centre.y
            else:
                rotate_point = rot_point

            new_ctrl_point_geom = affinity.rotate(
                self.assigned_control_point,
                angle,
                rotate_point,  # pyright: ignore [reportArgumentType]
                use_radians,
            )
            new_ctrl_point = new_ctrl_point_geom.x, new_ctrl_point_geom.y

        return Geometry(
            geom=affinity.rotate(self.geom, angle, rot_point, use_radians),  # pyright: ignore [reportArgumentType]
            material=self.material,
            control_points=new_ctrl_point,
            tol=self.tol,
        )

    def mirror_section(
        self,
        axis: str = "x",
        mirror_point: tuple[float, float] | str = "center",
    ) -> Geometry:
        """Mirrors the geometry about a point on either the x or y-axis.

        Args:
            axis: Axis about which to mirror the geometry, ``"x"`` or ``"y"``. Defaults
                to ``"x"``.
            mirror_point: Point about which to mirror the geometry (``x``, ``y``). If no
                point is provided, mirrors the geometry about the center of the shape's
                bounding box. Defaults to ``"center"``.

        Returns:
            New Geometry object mirrored on ``axis`` about ``mirror_point``

        Example:
            The following example mirrors a 200PFC section about the y-axis:

            .. plot::
                :include-source: True
                :caption: A 200PFC section mirrored about the y-axis and the point
                    ``(0, 0)``

                from sectionproperties.pre.library import channel_section

                geom = channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
                geom.mirror_section(axis="y", mirror_point=(0, 0)).plot_geometry()
        """
        x_mirror = 1.0
        y_mirror = 1.0
        m_pt: tuple[float, float] | str

        if mirror_point != "center" and not isinstance(mirror_point, str):
            x, y = mirror_point
            m_pt = (x, y)
        else:
            m_pt = mirror_point

        if axis == "x":
            x_mirror = -x_mirror
        elif axis == "y":
            y_mirror = -y_mirror

        mirrored_geom = affinity.scale(
            self.geom,
            xfact=y_mirror,
            yfact=x_mirror,
            zfact=1.0,
            origin=m_pt,  # pyright: ignore [reportArgumentType]
        )

        new_ctrl_point: tuple[float, float] | None = None

        if self.assigned_control_point:
            new_ctrl_point_geom = affinity.scale(
                self.assigned_control_point,
                xfact=y_mirror,
                yfact=x_mirror,
                zfact=1.0,
                origin=mirror_point,  # pyright: ignore [reportArgumentType]
            )
            new_ctrl_point = new_ctrl_point_geom.x, new_ctrl_point_geom.y

        return Geometry(
            geom=mirrored_geom,
            material=self.material,
            control_points=new_ctrl_point,
            tol=self.tol,
        )

    def split_section(
        self,
        point_i: tuple[float, float],
        point_j: tuple[float, float] | None = None,
        vector: tuple[float, float] | npt.NDArray[np.float64] | None = None,
    ) -> tuple[list[Geometry], list[Geometry]]:
        """Splits geometry about a line.

        Splits, or bisects, the geometry about a line, as defined by two points on the
        line or by one point on the line and a vector. Either ``point_j`` or ``vector``
        must be given. If ``point_j`` is given, ``vector`` is ignored.

        Returns a tuple of two lists each containing new Geometry instances representing
        the ``"top"`` and ``"bottom"`` portions, respectively, of the bisected geometry.

        If the line is a vertical line then the ``"right"`` and ``"left"`` portions,
        respectively, are returned.

        Args:
            point_i: A tuple of (``x``, ``y``) coordinates to define a first point on
                the line
            point_j: A tuple of (``x``, ``y``) coordinates to define a second point on
                the line. Defaults to ``None``.
            vector: A tuple or numpy array of (``x``, ``y``) components to define the
                line direction. Defaults to ``None``.

        Raises:
            ValueError: Line definition is invalid

        Returns:
            A tuple of lists containing Geometry objects that are bisected about the
            line defined by the two given points. The first item in the tuple represents
            the geometries on the ``"top"`` of the line (or to the ``"right"`` of the
            line, if vertical) and the second item represents the geometries to the
            ``"bottom"`` of the line (or to the ``"left"`` of the line, if vertical).

        Example:
            The following example splits a 200PFC section about the x-axis at ``y=100``:

            .. plot::
                :include-source: True
                :caption: A 200PFC section split about the y-axis.

                from sectionproperties.pre.library import channel_section

                geom = channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
                top_geoms, bot_geoms = geom.split_section(
                    point_i=(0, 100),
                    point_j=(1, 100),
                )
                (top_geoms[0] + bot_geoms[0]).plot_geometry()
        """
        if point_j:
            vector = np.array(
                [
                    point_j[0] - point_i[0],
                    point_j[1] - point_i[1],
                ]
            )
        elif vector is not None:
            vector = np.array(vector)
        else:
            msg = "Either a second point or a vector must be given to define the line."
            raise ValueError(msg)

        bounds = self.calculate_extents()
        line_segment = bisect.create_line_segment(
            point_on_line=point_i,
            vector=vector,
            bounds=bounds,
        )
        top_right_polys, bottom_left_polys = bisect.group_top_and_bottom_polys(
            polys=split(self.geom, line_segment),
            line=line_segment,
        )

        # Create new Geometrys from polys, preserve original material assignments
        top_right_geoms = [Geometry(poly, self.material) for poly in top_right_polys]
        bottom_left_geoms = [
            Geometry(poly, self.material) for poly in bottom_left_polys
        ]

        return top_right_geoms, bottom_left_geoms

    def offset_perimeter(
        self,
        amount: float = 0.0,
        where: str = "exterior",
        resolution: int = 12,
    ) -> Geometry | CompoundGeometry:
        """Dilates or erodes the section perimeter by a discrete amount.

        Args:
            amount: Distance to offset the section by. A negative value "erodes" the
                section. A positive value "dilates" the section. Defaults to ``0.0``.
            where: One of either ``"exterior"``, ``"interior"``, or ``"all"`` to specify
                which edges of the geometry to offset. If geometry has no interiors,
                then this parameter has no effect. Defaults to ``"exterior"``.
            resolution: Number of segments used to approximate a quarter circle around a
                point. Defaults to ``12``.

        Raises:
            ValueError: ``where`` is invalid
            ValueError: Attempted to offset internally where there are no holes
            RuntimeError: If geometry generation fails

        Returns:
            Geometry object translated to new alignment

        Example:
            The following example erodes a 200PFC section by 2 mm:

            .. plot::
                :include-source: True
                :caption: A 200PFC eroded by 2 mm.

                from sectionproperties.pre.library import channel_section

                geom = channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
                geom.offset_perimeter(amount=-2.0).plot_geometry()
        """
        if self.geom.interiors and where == "interior":
            exterior_polygon = Polygon(self.geom.exterior)

            for interior in self.geom.interiors:
                buffered_interior = buffer_polygon(
                    Polygon(interior), amount, resolution
                )
                exterior_polygon = exterior_polygon - buffered_interior

            if isinstance(exterior_polygon, MultiPolygon):
                return CompoundGeometry(
                    geoms=[
                        Geometry(poly, self.material) for poly in exterior_polygon.geoms
                    ]
                )
            elif not isinstance(exterior_polygon, Polygon):
                msg = "Geometry generation failed."
                raise RuntimeError(msg)

            # Check to see if assigned_control_point is still valid
            if self.assigned_control_point and exterior_polygon.contains(
                self.assigned_control_point
            ):
                return Geometry(
                    geom=exterior_polygon,
                    material=self.material,
                    control_points=self.control_points[0],
                    tol=self.tol,
                )

            return Geometry(
                geom=exterior_polygon,
                material=self.material,
                tol=self.tol,
            )
        elif not self.geom.interiors and where == "interior":
            msg = "Cannot buffer interior of Geometry object if it has no holes."
            raise ValueError(msg)
        elif where == "exterior":
            exterior_polygon = Polygon(self.geom.exterior)
            buffered_exterior = buffer_polygon(exterior_polygon, amount, resolution)

            for interior in self.geom.interiors:
                interior_poly = Polygon(interior)
                buffered_exterior = buffered_exterior - interior_poly

            if isinstance(buffered_exterior, MultiPolygon):
                return CompoundGeometry(
                    geoms=[
                        Geometry(poly, self.material)
                        for poly in buffered_exterior.geoms
                    ]
                )
            elif not isinstance(buffered_exterior, Polygon):
                msg = "Geometry generation failed."
                raise RuntimeError(msg)

            # Check to see if assigned_control_point is still valid
            if self.assigned_control_point and buffered_exterior.contains(
                self.assigned_control_point
            ):
                return Geometry(
                    geom=buffered_exterior,
                    material=self.material,
                    control_points=self.control_points[0],
                    tol=self.tol,
                )

            return Geometry(
                geom=buffered_exterior,
                material=self.material,
                tol=self.tol,
            )
        elif where == "all":
            buffered_geom = buffer_polygon(self.geom, amount, resolution)

            if isinstance(buffered_geom, MultiPolygon):
                compound_geom = CompoundGeometry(
                    geoms=[
                        Geometry(poly, self.material) for poly in buffered_geom.geoms
                    ]
                )
                return compound_geom

            # Check to see if assigned_control_point is still valid
            if self.assigned_control_point and buffered_geom.contains(
                self.assigned_control_point
            ):
                return Geometry(
                    geom=buffered_geom,
                    material=self.material,
                    control_points=self.control_points[0],
                    tol=self.tol,
                )
            return Geometry(
                geom=buffered_geom,
                material=self.material,
                tol=self.tol,
            )
        else:
            msg = "'where' must be 'exterior', 'interior' or 'all'."
            raise ValueError(msg)

    def shift_points(
        self,
        point_idxs: int | list[int],
        dx: float = 0.0,
        dy: float = 0.0,
        abs_x: float | None = None,
        abs_y: float | None = None,
    ) -> Geometry:
        """Shifts the points in the geometry.

        Translates one (or many points) in the geometry by either a relative amount or
        to a new absolute location. Returns a new Geometry representing the original
        with the selected point(s) shifted to the new location.

        Points are identified by their index, their relative location within the points
        list found in ``self.points``. You can call
        ``self.plot_geometry(labels="points")`` to see a plot with the points labeled to
        find the appropriate point indexes.

        Args:
            point_idxs: An integer representing an index location or a list of integer
                index locations.
            dx: The number of units in the x-direction to shift the point(s) by.
                Defaults to ``0.0``.
            dy: The number of units in the y-direction to shift the point(s) by.
                Defaults to ``0.0``.
            abs_x: Absolute x-coordinate in coordinate system to shift the point(s) to.
                If ``abs_x`` is provided, ``dx`` is ignored. If providing a list to
                ``point_idxs``, all points will be moved to this absolute location.
                Defaults to ``None``.
            abs_y: Absolute y-coordinate in coordinate system to shift the point(s) to.
                If ``abs_y`` is provided, ``dy`` is ignored. If providing a list to
                ``point_idxs``, all points will be moved to this absolute location.
                Defaults to ``None``.

        Returns:
            Geometry object with selected points translated to the new location

        Example:
            The following example expands the sides of a rectangle, one point at a time,
            to make it a square:

            .. plot::
                :include-source: True
                :caption: Rectangle transformed to a square with ``shift_points``

                from sectionproperties.pre.library import rectangular_section

                geom = rectangular_section(d=200, b=150)

                # using relative shifting
                geom_step_1 = geom.shift_points(point_idxs=1, dx=50)

                # using absolute relocation
                geom_step_1.shift_points(point_idxs=2, abs_x=200).plot_geometry()
        """
        current_points = copy.copy(self.points)
        current_facets = copy.copy(self.facets)
        current_holes = copy.copy(self.holes)
        current_material = copy.copy(self.material)
        current_control_points = copy.copy(self.control_points)

        if isinstance(point_idxs, int):
            point_idxs = [point_idxs]

        new_x, new_y = None, None

        for point_idx in point_idxs:
            current_x, current_y = current_points[point_idx]
            new_x = abs_x if abs_x else current_x + dx
            new_y = abs_y if abs_y else current_y + dy
            current_points[point_idx] = (new_x, new_y)

        new_geom = Geometry.from_points(
            points=current_points,
            facets=current_facets,
            control_points=current_control_points,
            holes=current_holes,
            material=current_material,
        )
        if self.assigned_control_point and new_geom.geom.contains(
            self.assigned_control_point
        ):
            new_geom = Geometry.from_points(
                points=current_points,
                facets=current_facets,
                control_points=current_control_points,
                holes=current_holes,
                material=current_material,
            )
        return new_geom

    def plot_geometry(
        self,
        labels: tuple[str] = ("control_points",),
        title: str = "Cross-Section Geometry",
        cp: bool = True,
        legend: bool = True,
        **kwargs: Any,
    ) -> matplotlib.axes.Axes:
        """Plots the geometry defined by the input section.

        Args:
            labels: A tuple of str which indicate which labels to plot. Can be one or a
                combination of ``"points"``, ``"facets"``, ``"control_points"``, or an
                empty tuple to indicate no labels. Defaults to ``("control_points")``.
            title: Plot title. Defaults to ``"Cross-Section Geometry"``.
            cp: If set to True, plots the control points. Defaults to ``True``.
            legend: If set to True, plots the legend. Defaults to ``True``.
            kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context()`

        Returns:
            Matplotlib axes object
        """
        # create plot and setup the plot
        with post.plotting_context(title=title, **kwargs) as (_, ax):
            if not ax:
                msg = "Matplotlib axes not created."
                raise RuntimeError(msg)

            # plot the points and facets
            for i, f in enumerate(self.facets):
                label = "Points & Facets" if i == 0 else None
                ax.plot(  # pyright: ignore
                    [self.points[f[0]][0], self.points[f[1]][0]],
                    [self.points[f[0]][1], self.points[f[1]][1]],
                    "ko-",
                    markersize=2,
                    linewidth=1.5,
                    label=label,
                )

            # plot the holes
            for i, h in enumerate(self.holes):
                label = "Holes" if i == 0 else None
                ax.plot(h[0], h[1], "rx", markersize=5, markeredgewidth=1, label=label)  # pyright: ignore

            if cp:
                # plot the control points
                for i, cp_i in enumerate(self.control_points):
                    label = "Control Points" if i == 0 else None
                    ax.plot(cp_i[0], cp_i[1], "bo", markersize=5, label=label)  # pyright: ignore

            # display the labels
            for label in labels:
                # plot control_point labels
                if label == "control_points":
                    for i, pt in enumerate(self.control_points):
                        ax.annotate(str(i), xy=pt, color="b")  # pyright: ignore

                # plot point labels
                if label == "points":
                    for i, pt in enumerate(self.points):
                        ax.annotate(str(i), xy=pt, color="r")  # pyright: ignore

                # plot facet labels
                if label == "facets":
                    for i, fct in enumerate(self.facets):
                        pt1 = self.points[fct[0]]
                        pt2 = self.points[fct[1]]
                        xy = ((pt1[0] + pt2[0]) / 2, (pt1[1] + pt2[1]) / 2)

                        ax.annotate(str(i), xy=xy, color="b")  # pyright: ignore

                # plot hole labels
                if label == "holes":
                    for i, pt in enumerate(self.holes):
                        ax.annotate(str(i), xy=pt, color="r")  # pyright: ignore

            # display the legend
            if legend:
                ax.legend(loc="center left", bbox_to_anchor=(1, 0.5))  # pyright: ignore

        return ax

    def calculate_extents(self) -> tuple[float, float, float, float]:
        """Calculate geometry extents.

        Calculates the minimum and maximum x and y-values amongst the list of points;
        the points that describe the bounding box of the Geometry instance.

        Returns:
            Minimum and maximum x and y-values (``x_min``, ``x_max``, ``y_min``,
            ``y_max``)
        """
        min_x, min_y, max_x, max_y = self.geom.bounds
        return min_x, max_x, min_y, max_y

    def calculate_perimeter(self) -> float:
        """Calculates the exterior perimeter of the geometry.

        Returns:
            Geometry perimeter
        """
        return float(self.geom.exterior.length)

    def calculate_area(self) -> float:
        """Calculates the area of the geometry.

        Returns:
            Geometry area
        """
        return float(self.geom.area)

    def calculate_centroid(self) -> tuple[float, float]:
        """Calculates the centroid of the geometry.

        Returns:
            Geometry centroid
        """
        return self.geom.centroid.x, self.geom.centroid.y

    @property
    def recovery_points(self) -> list[tuple[float, float]] | list[Point]:
        """Recovery points.

        Returns four stress recovery points for the section geometry. If the Geometry
        instance was created by a NASTRAN geometry function,
        e.g. :func:`~sectionproperties.pre.nastran_sections.nastran_bar()`, then the
        recovery points will be pre-set on the Geometry instance.

        Returns:
            Stress recovery poiints.
        """
        return self._recovery_points

    @recovery_points.setter
    def recovery_points(
        self,
        new_points: list[tuple[float, float]] | list[Point],
    ) -> None:
        """Sets the recovery points.

        Args:
            new_points: Recovering points to set
        """
        self._recovery_points = new_points

    @recovery_points.getter
    def recovery_points(self) -> list[tuple[float, float]] | list[Point]:
        """Returns the recovery points.

        Returns:
            Recovery points
        """
        return self._recovery_points

    def __or__(
        self,
        other: Geometry | CompoundGeometry,
    ) -> Geometry | CompoundGeometry:
        """Performs a union on Geometry objects with the ``|`` operator.

        Args:
            other: Geometry object to perform the union with

        Raises:
            ValueError: Unable to perform union

        Returns:
            New Geometry object

        Example:
            The following example combines two rectangles using the ``|`` operator:

            .. plot::
                :include-source: True
                :caption: Union of two rectangles

                from sectionproperties.pre.library import rectangular_section

                rect1 = rectangular_section(d=100, b=200)
                rect2 = rectangular_section(d=100, b=200).shift_section(
                    x_offset=150,
                    y_offset=70,
                )
                (rect1 | rect2).plot_geometry()
        """
        material = self.material or other.material

        try:
            or_poly = self.geom | other.geom

            if isinstance(
                or_poly,
                GeometryCollection | LineString | Point | Polygon | MultiPolygon,
            ):
                new_polygon = filter_non_polygons(input_geom=or_poly)  # pyright: ignore [reportUnknownArgumentType]
            else:
                msg = "Geometry generation failed."
                raise RuntimeError(msg)

            if isinstance(new_polygon, MultiPolygon):
                return CompoundGeometry(
                    geoms=[Geometry(polygon, material) for polygon in new_polygon.geoms]
                )
            return Geometry(
                geom=new_polygon,
                material=material,
                control_points=self.control_points[0],
                tol=self.tol,
            )
        except Exception as e:
            msg = f"Cannot perform 'union' on these two objects: {self} | {other}"
            raise ValueError(msg) from e

    def __xor__(
        self,
        other: Geometry | CompoundGeometry,
    ) -> Geometry | CompoundGeometry:
        """Performs a symmetric difference on Geometry objects with the ``^`` operator.

        Returns the regions of geometry that are not overlapping.

        Args:
            other: Geometry object to perform the symmetric difference with

        Raises:
            ValueError: Unable to perform symmetric difference

        Returns:
            New Geometry object

        Example:
            The following example performs a symmetric difference on two circles with
            the ``^`` operator:

            .. plot::
                :include-source: True
                :caption: Symmetric difference of two circles

                from sectionproperties.pre.library import circular_section
                from sectionproperties.analysis import Section

                circ1 = circular_section(d=100, n=64)
                circ2 = circular_section(d=100, n=64).shift_section(x_offset=35)
                (circ1 ^ circ2).plot_geometry()
        """
        material = self.material or other.material

        try:
            xor_poly = self.geom ^ other.geom

            if isinstance(
                xor_poly,
                GeometryCollection | LineString | Point | Polygon | MultiPolygon,
            ):
                new_polygon = filter_non_polygons(xor_poly)  # pyright: ignore [reportUnknownArgumentType]
            else:
                msg = "Geometry generation failed."
                raise RuntimeError(msg)

            if isinstance(new_polygon, MultiPolygon):
                return CompoundGeometry(
                    geoms=[Geometry(polygon, material) for polygon in new_polygon.geoms]
                )

            if isinstance(new_polygon, GeometryCollection):
                msg = "Cannot perform 'symmetric difference' on these two objects: "
                msg += f"{self} ^ {other}"
                raise ValueError(msg)

            return Geometry(
                geom=new_polygon,
                material=material,
                tol=self.tol,
            )
        except Exception as e:
            msg = "Cannot perform 'symmetric difference' on these two objects: "
            msg += f"{self} ^ {other}"
            raise ValueError(msg) from e

    def __sub__(
        self,
        other: Geometry | CompoundGeometry,
    ) -> Geometry | CompoundGeometry:
        """Performs a difference operation on Geometry objects with the ``-`` operator.

        Subtracts the second geometry from the first geometry.

        Args:
            other: Geometry object to perform the difference operation with

        Raises:
            ValueError: Unable to perform difference

        Returns:
            New Geometry object

        Example:
            The following example creates a hollow box using the ``-`` operator:

            .. plot::
                :include-source: True
                :caption: Hollow box

                from sectionproperties.pre.library import rectangular_section

                rect1 = rectangular_section(d=400, b=200)
                rect2 = rectangular_section(d=300, b=100).align_center(align_to=rect1)
                (rect1 - rect2).plot_geometry()
        """
        if isinstance(self, CompoundGeometry):
            subs_geom_acc: list[Geometry | CompoundGeometry] = []

            for geom in self.geoms:
                new_geom = geom - other
                subs_geom_acc.append(new_geom)

            return CompoundGeometry(geoms=subs_geom_acc)
        else:
            material = self.material or pre.DEFAULT_MATERIAL

            try:
                sub_poly = self.geom - other.geom

                if isinstance(
                    sub_poly,
                    GeometryCollection | LineString | Point | Polygon | MultiPolygon,
                ):
                    new_polygon = filter_non_polygons(sub_poly)  # pyright: ignore [reportUnknownArgumentType]
                else:
                    msg = "Geometry generation failed."
                    raise RuntimeError(msg)

                if isinstance(new_polygon, GeometryCollection):  # Non-polygon results
                    msg = "Cannot perform 'difference' on these two objects: "
                    msg += f"{self} - {other}"
                    raise ValueError(msg)
                elif isinstance(new_polygon, MultiPolygon):
                    return CompoundGeometry(
                        geoms=[
                            Geometry(polygon, material) for polygon in new_polygon.geoms
                        ]
                    )
                elif isinstance(new_polygon, Polygon):  # pyright: ignore [reportUnnecessaryIsInstance]
                    if self.assigned_control_point and new_polygon.contains(
                        self.assigned_control_point
                    ):
                        return Geometry(
                            geom=new_polygon,
                            material=material,
                            control_points=self.assigned_control_point,
                            tol=self.tol,
                        )
                    else:
                        return Geometry(
                            geom=new_polygon,
                            material=material,
                            tol=self.tol,
                        )
                else:
                    msg = "Cannot perform 'difference' on these two objects: "
                    msg += f"{self} - {other}"
                    raise ValueError(msg)
            except Exception as e:
                msg = "Cannot perform 'difference' on these two objects: "
                msg += f"{self} - {other}"
                raise ValueError(msg) from e

    def __add__(
        self,
        other: Geometry | CompoundGeometry,
    ) -> CompoundGeometry:
        """Combine Geometry objects into a CompoundGeometry using the ``+`` operator.

        Args:
            other: Geometry object to perform the combination with

        Raises:
            ValueError: Unable to perform combination operation

        Returns:
            New Geometry object

        Example:
            The following example creates a tee section the ``+`` operator:

            .. plot::
                :include-source: True
                :caption: Tee section

                from sectionproperties.pre.library import rectangular_section

                flange = rectangular_section(d=16, b=200)
                web = (
                    rectangular_section(d=284, b=16)
                    .align_center(align_to=flange)
                    .align_to(other=flange, on="bottom")
                )
                (flange + web).plot_geometry()
        """
        try:
            return CompoundGeometry(geoms=[self, other])
        except Exception as e:
            msg = "Cannot create new CompoundGeometry with these objects: "
            msg += f"{self} + {other}"
            raise ValueError(msg) from e

    def __and__(
        self,
        other: Geometry | CompoundGeometry,
    ) -> Geometry | CompoundGeometry:
        """Performs an intersection on Geometry objects with the ``&`` operator.

        Returns the regions of geometry common to both geometries.

        Args:
            other: Geometry object to perform the intersection with

        Raises:
            ValueError: Unable to perform intersection

        Returns:
            New Geometry object

        Example:
            The following example performs an intersection of a square and circle using
            the ``&`` operator:

            .. plot::
                :include-source: True
                :caption: Intersection of a square and circle

                from sectionproperties.pre.library import rectangular_section
                from sectionproperties.pre.library import circular_section

                rect = rectangular_section(d=200, b=200).align_center(align_to=(0, 0))
                circle = circular_section(d=250, n=64)
                (rect & circle).plot_geometry()
        """
        material = self.material or other.material

        try:
            and_poly = self.geom & other.geom

            if isinstance(
                and_poly,
                GeometryCollection | LineString | Point | Polygon | MultiPolygon,
            ):
                new_polygon = filter_non_polygons(and_poly)  # pyright: ignore [reportUnknownArgumentType]
            else:
                msg = "Geometry generation failed."
                raise RuntimeError(msg)

            if isinstance(new_polygon, MultiPolygon):
                return CompoundGeometry(
                    geoms=[Geometry(polygon, material) for polygon in new_polygon.geoms]
                )

            return Geometry(
                geom=new_polygon,
                material=material,
                tol=self.tol,
            )
        except Exception as e:
            msg = "Cannot perform 'intersection' on these two Geometry instances: "
            msg += f"{self} & {other}"
            raise ValueError(msg) from e


class CompoundGeometry(Geometry):
    """Class for defining a geometry of multiple distinct regions.

    CompoundGeometry instances are composed of multiple Geometry objects. As with
    Geometry objects, CompoundGeometry objects have methods for generating a triangular
    mesh over all geometries, transforming the collection of geometries as though they
    were one (e.g. translation, rotation, and mirroring), and aligning the
    CompoundGeometry to another Geometry (or to another CompoundGeometry).

    Each Geometry object may have different material properties.

    CompoundGeometry objects can be created directly between two or more Geometry
    objects by using the ``+`` operator.
    """

    def __init__(
        self,
        geoms: MultiPolygon | list[Geometry],
    ) -> None:
        """Inits the CompoundGeometry class.

        Args:
            geoms: Either a list of Geometry objects or a :class:`shapely.MultiPolygon`
                instance.
        """
        if isinstance(geoms, MultiPolygon):
            self.geoms = [
                Geometry(poly, material=pre.DEFAULT_MATERIAL) for poly in geoms.geoms
            ]
            self.geom = geoms
        elif isinstance(geoms, list):  # pyright: ignore [reportUnnecessaryIsInstance]
            processed_geoms: list[Geometry] = []

            for item in geoms:
                if isinstance(item, CompoundGeometry):
                    # "Flatten" any CompoundGeometry objects in the 'geoms' list
                    processed_geoms += item.geoms
                elif isinstance(item, Geometry):  # pyright: ignore [reportUnnecessaryIsInstance]
                    processed_geoms.append(item)

            self.geoms = processed_geoms
            self.geom = MultiPolygon([geom.geom for geom in processed_geoms])

        self.control_points = []
        self.points = []
        self.facets = []
        self.holes = []
        self.material = None
        self.compile_geometry()
        self.tol = 12
        self.mesh: dict[str, Any] | None = None

    def _repr_svg_(self) -> str:
        """Returns an svg representation of the CompoundGeometry.

        Wraps :func:`shapely.MultiPolygon._repr_svg_()` by returning
        ``self.geom._repr_svg_()``

        Returns:
            Representative svg
        """
        materials_list = [geom.material.name for geom in self.geoms]
        print("sectionproperties.pre.geometry.CompoundGeometry")
        print(f"object at: {hex(id(self))}")
        print(f"Materials incl.: {list(set(materials_list))}")

        return str(self.geom._repr_svg_())

    @staticmethod
    def from_points(  # pyright: ignore [reportIncompatibleMethodOverride]
        points: list[tuple[float, float]],
        facets: list[tuple[int, int]],
        control_points: list[tuple[float, float]],
        holes: list[tuple[float, float]] | None = None,
        materials: list[pre.Material] | None = None,
    ) -> CompoundGeometry:
        """Creates CompoundGeometry from points, facets, control points and holes.

        An interface for the creation of CompoundGeometry objects through the definition
        of points, facets, control points and holes. Geometries created through this
        method are expected to be non-ambiguous meaning that no "overlapping" geometries
        exists and that nodal connectivity is maintained (e.g. there are no nodes
        "overlapping" with facets without nodal connectivity).

        Args:
            points: List of points (``x``, ``y``) defining the vertices of the
                section geometry.
            facets: A list of (``start``, ``end``) indices of vertices defining the
                edges of the section geoemtry. Can be used to define both external and
                internal perimeters of holes. Facets are assumed to be described in the
                order of exterior perimeter, interior perimeter 1, interior perimeter 2,
                etc.
            control_points: A list of points (``x``, ``y``) that define regions as being
                distinct, contiguous, and having one material. The point can be located
                anywhere within region. Only one point is permitted per region. The
                order of ``control_points`` must be given in the same order as the order
                that polygons are created by ``facets``.
            holes: A list of points (``x``, ``y``) that define interior regions as
                being holes or voids. The point can be located anywhere within the hole
                region. Only one point is required per hole region. Defaults to
                ``None``.
            materials: A list of :class:`~sectionproperties.pre.pre.Material` objects
                that are to be assigned, in order, to the regions defined by the given
                ``control_points``. If not given, then the
                :class:`~sectionproperties.pre.DEFAULT_MATERIAL` will be used for
                each region. Defaults to ``None``.

        Raises:
            ValueError: If there are materials provided without control points
            ValueError: If the number of materials does not equal the number of control
                points
            ValueError: If the number of exterior regions doesn't match the number of
                control points
            ValueError: If control points are not contained within geometries with holes

        Returns:
            CompoundGeometry object from points

        Example:
            This example creates two regions with different material properties:

            .. plot::
                :include-source: True
                :caption: Composite CompoundGeometry created from points and facets.

                from sectionproperties.pre import Material, CompoundGeometry
                from sectionproperties.analysis import Section

                mat1 = Material(
                    name="mat1",
                    elastic_modulus=1.0,
                    poissons_ratio=0.0,
                    density=1.0,
                    yield_strength=1.0,
                    color="tab:olive",
                )

                mat2 = Material(
                    name="mat2",
                    elastic_modulus=2.0,
                    poissons_ratio=0.0,
                    density=2.0,
                    yield_strength=2.0,
                    color="tab:purple",
                )

                points = [(0, 0), (10, 0), (15, 10), (-12, -5)]
                facets = [(0, 1), (1, 2), (2, 0), (1, 3), (3, 0), (0, 1)]
                control_points = [(1, 1), (1, -1)]
                materials = [mat1, mat2]

                geom = CompoundGeometry.from_points(
                    points=points,
                    facets=facets,
                    control_points=control_points,
                    materials=materials,
                )
                geom.create_mesh(mesh_sizes=[0])
                Section(geometry=geom).plot_mesh()
        """
        if materials is None:
            # assign default material to each region
            materials = [pre.DEFAULT_MATERIAL] * len(control_points)

        if materials and not control_points:
            msg = "Materials cannot be assigned without control_points. Please provide "
            msg += "corresponding control_points for each material."
            raise ValueError(msg)

        if holes is None:
            holes = []

        if len(materials) != len(control_points):
            msg = "If materials are provided, the number of materials in the list "
            msg += "must match the number of control_points provided."
            msg += f"len(materials)=={len(materials)}, "
            msg += f"len(control_points)=={len(control_points)}."
            raise ValueError(msg)

        # First, generate all invidual polygons from points and facets
        current_polygon_points = []
        all_polygons: list[Polygon] = []
        prev_facet = None

        # initialise indices
        i_idx, j_idx = 0, 0

        for facet in facets:
            i_idx, j_idx = facet

            if not prev_facet:  # Add the first facet vertex to exterior and move on
                current_polygon_points.append(points[i_idx])
                prev_facet = facet
                continue

            prev_j_idx = prev_facet[1]

            if i_idx != prev_j_idx:  # If there is a break in the chain of edges...
                # ... then add the last point, close off the polygon,
                # and add the polygon to the all_polygons accumulator....
                current_polygon_points.append(points[prev_j_idx])
                all_polygons.append(Polygon(current_polygon_points))
                # Then start collecting the points of the new polygon
                current_polygon_points = [points[i_idx]]
            else:
                current_polygon_points.append(
                    points[i_idx]
                )  # Only need i_idx b/c shapely auto-closes polygons
            prev_facet = facet
        # Use the for...else clause to add the last point and close the last polygon.
        else:
            current_polygon_points.append(points[j_idx])
            all_polygons.append(Polygon(current_polygon_points))

        # Then classify all of the collected polygons as either "exterior" or "interior"
        exteriors: list[Polygon] = []
        interiors: list[Polygon] = []

        for polygon in all_polygons:
            hole_coord_in_polygon = [
                hole_coord
                for hole_coord in holes
                if polygon.contains(Point(hole_coord))
            ]
            ctrl_coord_in_polygon = [
                ctrl_coord
                for ctrl_coord in control_points
                if polygon.contains(Point(ctrl_coord))
            ]

            if any(hole_coord_in_polygon) and not any(ctrl_coord_in_polygon):
                interiors.append(polygon)
            else:
                exteriors.append(polygon)

        # Create the holes by subtracting interior regions from exterior regions
        if len(exteriors) != len(control_points):
            msg = f"The number of exterior regions ({len(exteriors)}) "
            msg += "does not match the number of control_points given "
            msg += f"({len(control_points)})."
            raise ValueError(msg)

        if not interiors:
            return CompoundGeometry(
                geoms=[
                    Geometry(
                        geom=exterior,
                        control_points=control_points[idx],
                        material=materials[idx],
                    )
                    for idx, exterior in enumerate(exteriors)
                ]
            )
        else:
            # "Punch" all holes through each exterior geometry
            punched_exterior_geometries: list[Geometry] = []

            for idx, exterior in enumerate(exteriors):
                punched_exterior = exterior
                exterior_control_point = 0.0, 0.0  # initialise

                for interior in interiors:
                    punched_exterior = punched_exterior - interior
                    try:
                        exterior_control_point = next(
                            control_point
                            for control_point in control_points
                            if punched_exterior.contains(Point(control_point))
                        )
                    except StopIteration as e:
                        msg = "Control points given are not contained within the "
                        msg += f"geometry once holes are subtracted: {control_points}"
                        raise ValueError(msg) from e

                if not isinstance(punched_exterior, Polygon):
                    msg = "Geometry generation failed."
                    raise RuntimeError(msg)

                exterior_geometry = Geometry(
                    geom=punched_exterior,
                    control_points=exterior_control_point,
                    material=materials[idx],
                )
                punched_exterior_geometries.append(exterior_geometry)

            return CompoundGeometry(punched_exterior_geometries)

    @classmethod
    def from_3dm(
        cls,
        filepath: str | pathlib.Path,
        **kwargs: Any,
    ) -> CompoundGeometry:
        """Creates a CompoundGeometry object from the objects in a Rhino ``3dm`` file.

        Args:
            filepath: File path to the rhino ``.3dm`` file.
            kwargs: See below.

        Keyword Args:
            refine_num (Optional[int]):  Bézier curve interpolation number. In Rhino a
                surface's edges are nurb based curves. Shapely does not support nurbs,
                so the individual Bézier curves are interpolated using straight lines.
                This parameter sets the number of straight lines used in the
                interpolation. Default is 1.
            vec1 (Optional[numpy.ndarray]): A 3d vector in the Shapely plane. Rhino is a
                3D geometry environment. Shapely is a 2D geometric library. Thus a 2D
                plane needs to be defined in Rhino that represents the Shapely
                coordinate system. ``vec1`` represents the 1st vector of this plane. It
                will be used as Shapely's x direction. Default is [1,0,0].
            vec2 (Optional[numpy.ndarray]): Continuing from ``vec1``, ``vec2`` is
                another vector to define the Shapely plane. It must not be [0,0,0] and
                it's only requirement is that it is any vector in the Shapely plane (but
                not equal to ``vec1``). Default is [0,1,0].
            plane_distance (Optional[float]): The distance to the Shapely plane. Default
                is 0.
            project (Optional[bool]): Controls if the breps are projected onto the plane
                in the direction of the Shapley plane's normal. Default is True.
            parallel (Optional[bool]): Controls if only the rhino surfaces that have the
                same normal as the Shapely plane are yielded. If true, all non parallel
                surfaces are filtered out. Default is False.

        Raises:
            ImportError: If ``rhino3dm`` is not installed. To enable rhino features use
                ``pip install sectionproperties[rhino]``.

        Returns:
            CompoundGeometry object
        """
        try:
            import sectionproperties.pre.rhino as rhino_importer
        except ImportError as e:
            msg = "There is something wrong with your rhino library installation. "
            msg += "Please report this error at "
            msg += "https://github.com/robbievanleeuwen/section-properties/issues"
            raise ImportError(msg) from e

        list_poly = rhino_importer.load_3dm(filepath, **kwargs)

        return cls(geoms=MultiPolygon(list_poly))

    def create_mesh(
        self,
        mesh_sizes: float | list[float],
        min_angle: float = 30.0,
        coarse: bool = False,
    ) -> CompoundGeometry:
        """Creates a quadratic triangular mesh from the CompoundGeometry object.

        Args:
            mesh_sizes: A float describing the maximum mesh element area to be used in
                the finite-element mesh for each Geometry object within the
                CompoundGeometry object. If a list of length 1 or a float is passed,
                then the one size will be applied to all constituent Geometry meshes.
            min_angle: The meshing algorithm adds vertices to the mesh to ensure that no
                angle smaller than the minimum angle (in degrees, rounded to 1 decimal
                place). Note that small angles between input segments cannot be
                eliminated. If the minimum angle is 20.7 deg or smaller, the
                triangulation algorithm is theoretically guaranteed to terminate (given
                sufficient precision). The algorithm often doesn't terminate for angles
                greater than 33 deg. Some meshes may require angles well below 20 deg to
                avoid problems associated with insufficient floating-point precision.
                Defaults to ``30.0``.
            coarse: If set to True, will create a coarse mesh (no area or quality
                constraints). Defaults to ``False``.

        Returns:
            ``CompoundGeometry`` object with mesh data stored in ``.mesh`` attribute.
            Returned ``CompoundGeometry`` object is self, not a new instance.

        Example:
            The following example creates a mesh for a plate with a hole, with a refined
            mesh around the hole region:

            .. plot::
                :include-source: True
                :caption: Mesh refinement around hole

                from sectionproperties.pre.library import rectangular_section
                from sectionproperties.pre.library import circular_section
                from sectionproperties.analysis import Section

                outer_rect = rectangular_section(d=100, b=200)
                inner_rect = rectangular_section(d=50, b=50).align_center(
                    align_to=outer_rect
                )
                circle = circular_section(d=20, n=32).align_center(align_to=outer_rect)

                geom = outer_rect - inner_rect + inner_rect - circle
                geom.create_mesh(mesh_sizes=[100, 5])
                Section(geometry=geom).plot_mesh(materials=False)
        """
        if isinstance(mesh_sizes, float | int):
            mesh_sizes = [mesh_sizes]

        if len(mesh_sizes) == 1:
            mesh_sizes = mesh_sizes * len(self.control_points)

        self.mesh = pre.create_mesh(
            points=self.points,
            facets=self.facets,
            holes=self.holes,
            control_points=self.control_points,
            mesh_sizes=mesh_sizes,
            min_angle=min_angle,
            coarse=coarse,
        )

        return self

    def shift_section(
        self,
        x_offset: float = 0.0,
        y_offset: float = 0.0,
    ) -> CompoundGeometry:
        """Returns a CompoundGeometry object translated by (``x_offset``, ``y_offset``).

        Args:
            x_offset: Distance in x-direction by which to shift the geometry. Defaults
                to ``0.0``.
            y_offset: Distance in y-direction by which to shift the geometry. Defaults
                to ``0.0``.

        Returns:
            New Geometry object shifted by ``x_offset`` and ``y_offset``
        """
        geoms_acc = [
            geom.shift_section(x_offset=x_offset, y_offset=y_offset)
            for geom in self.geoms
        ]

        return CompoundGeometry(geoms=geoms_acc)

    def rotate_section(
        self,
        angle: float,
        rot_point: tuple[float, float] | str = "center",
        use_radians: bool = False,
    ) -> CompoundGeometry:
        """Rotates a CompoundGeometry object.

        Rotates the CompoundGeometry and specified angle about a point. If the rotation
        point is not provided, rotates the section about the center of the
        CompoundGeometry's bounding box.

        Args:
            angle: Angle (degrees by default) by which to rotate the section. A
                positive angle leads to a counter-clockwise rotation.
            rot_point: Point (``x``, ``y``) about which to rotate the section. If not
                provided, will rotate about the center of the compound geometry's
                bounding box. Defaults to ``"center"``.
            use_radians: Boolean to indicate whether ``angle`` is in degrees or radians.
                If True, ``angle`` is interpreted as radians. Defaults to ``False``.

        Returns:
            CompoundGeometry object rotated by ``angle`` about ``rot_point``

        Example:
            The following example rotates a 200UB25 section with a plate clockwise by 30
            degrees:

            .. plot::
                :include-source: True
                :caption: Reinforced 200UB25 rotated.

                from sectionproperties.pre.library import rectangular_section, i_section

                geom1 = i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
                geom2 = rectangular_section(d=20, b=133)
                compound = geom1 + geom2.align_center(align_to=geom1).align_to(
                    other=geom1, on="top"
                )
                compound.rotate_section(angle=-30).plot_geometry()
        """
        if rot_point == "center":
            centre = box(*MultiPolygon(self.geom).bounds).centroid
            rp = centre.x, centre.y
        else:
            rp = rot_point

        geoms_acc = [
            geom.rotate_section(angle=angle, rot_point=rp, use_radians=use_radians)
            for geom in self.geoms
        ]

        return CompoundGeometry(geoms=geoms_acc)

    def mirror_section(
        self,
        axis: str = "x",
        mirror_point: tuple[float, float] | str = "center",
    ) -> CompoundGeometry:
        """Mirrors the CompoundGeometry object about a point on either the x or y-axis.

        Args:
            axis: Axis about which to mirror the geometry, ``"x"`` or ``"y"``. Defaults
                to ``"x"``.
            mirror_point: Point about which to mirror the geometry (``x``, ``y``). If no
                point is provided, mirrors the geometry about the center of the shape's
                bounding box. Defaults to ``"center"``.

        Returns:
            CompoundGeometry object mirrored on ``axis`` about ``mirror_point``

        Example:
            The following example mirrors a 200PFC section with a plate about the y-axis
            and the point ``(0, 0)``:

            .. plot::
                :include-source: True
                :caption: Reinforced 200PFC mirrored.

                from sectionproperties.pre.library import rectangular_section
                from sectionproperties.pre.library import channel_section

                geom1 = channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
                geom2 = rectangular_section(d=20, b=133)
                compound = geom1 + geom2.align_center(align_to=geom1).align_to(
                    other=geom1, on="top"
                )
                compound.mirror_section(axis="y", mirror_point=(0,0)).plot_geometry()
        """
        geoms_acc = [
            geom.mirror_section(axis=axis, mirror_point=mirror_point)
            for geom in self.geoms
        ]

        return CompoundGeometry(geoms=geoms_acc)

    def align_center(
        self,
        align_to: Geometry | tuple[float, float] | None = None,
    ) -> CompoundGeometry:
        """Aligns the CompoundGeometry to a center point.

        Returns a new CompoundGeometry object, translated in both ``x`` and ``y``, so
        that the center-point of the new object's material-weighted centroid will be
        aligned with centroid of the object in ``align_to``. If ``align_to`` is an
        ``x``, ``y`` coordinate, then the centroid will be aligned to the coordinate.
        If ``align_to`` is None then the new object will be aligned with its centroid at
        the origin.

        .. note::

            The material-weighted centroid refers to when individual geometries within
            the ``CompoundGeometry`` object have been assigned differing materials. The
            centroid of the compound geometry is calculated by using the elastic modulus
            of each geometry's assigned material.

        Args:
            align_to: Another Geometry to align to, an (``x``, ``y``) coordinate, or
                ``None``. Defaults to ``None``.

        Raises:
            ValueError: ``align_to`` is not valid

        Returns:
            CompoundGeometry object translated to new alignment

        Example:
            The following example creates a rectanglular steel-concrete composite
            section and uses the ``align_center()`` method to place the composite
            centroid at the origin.

            .. plot::
                :include-source: True
                :caption: Reinforced 200PFC mirrored.

                from sectionproperties.pre import Material
                from sectionproperties.pre.library import rectangular_section
                from sectionproperties.analysis import Section

                steel = Material(
                    name="Steel",
                    elastic_modulus=200e3,
                    poissons_ratio=0.3,
                    density=7.85e-6,
                    yield_strength=250,
                    color="grey",
                )
                concrete = Material(
                    name="Concrete",
                    elastic_modulus=30.1e3,
                    poissons_ratio=0.2,
                    density=2.4e-6,
                    yield_strength=32,
                    color="lightgrey",
                )

                geom_steel = rectangular_section(d=50, b=50, material=steel)
                geom_timber = rectangular_section(d=50, b=50, material=concrete)
                geom = geom_timber.align_to(geom_steel, on="right") + geom_steel
                geom = geom.align_center()
                geom.create_mesh(mesh_sizes=[10, 5])
                Section(geometry=geom).plot_mesh()
        """
        ea_sum = sum(
            [
                geom.material.elastic_modulus * geom.calculate_area()
                for geom in self.geoms
            ]
        )

        cx_ea_acc = 0.0
        cy_ea_acc = 0.0

        for geom in self.geoms:
            e = geom.material.elastic_modulus
            a = geom.calculate_area()
            ea = e * a
            cx, cy = list(geom.geom.centroid.coords[0])
            cx_ea_acc += cx * ea
            cy_ea_acc += cy * ea

        weighted_cx = cx_ea_acc / ea_sum
        weighted_cy = cy_ea_acc / ea_sum

        if align_to is None:
            shift_x, shift_y = (
                round(-weighted_cx, self.tol),
                round(-weighted_cy, self.tol),
            )

        elif isinstance(align_to, Geometry):
            align_cx, align_cy = next(iter(align_to.geom.centroid.coords))
            shift_x = round(align_cx - weighted_cx, self.tol)
            shift_y = round(align_cy - weighted_cy, self.tol)

        else:
            try:
                point_x, point_y = align_to
                shift_x = round(point_x - weighted_cx, self.tol)
                shift_y = round(point_y - weighted_cy, self.tol)
            except (TypeError, ValueError) as e:
                msg = "align_to must be either a Geometry object or an x, y "
                msg += f"coordinate, not {align_to}."
                raise ValueError(msg) from e

        return self.shift_section(x_offset=shift_x, y_offset=shift_y)

    def split_section(
        self,
        point_i: tuple[float, float],
        point_j: tuple[float, float] | None = None,
        vector: tuple[float, float] | None | npt.NDArray[np.float64] = None,
    ) -> tuple[list[Geometry], list[Geometry]]:
        """Splits CompoundGeometry about a line.

        Splits, or bisects, the geometry about an infinite line, as defined by two
        points on the line or by one point on the line and a vector. Either ``point_j``
        or ``vector`` must be given. If ``point_j`` is given, ``vector`` is ignored.

        Returns a tuple of two lists each containing new Geometry instances representing
        the ``top`` and ``bottom`` portions, respectively, of the bisected geometry.

        If the line is a vertical line then the ``right`` and ``left`` portions,
        respectively, are returned.

        Args:
            point_i: A tuple of (``x``, ``y``) coordinates to define a first point on
                the line
            point_j: A tuple of (``x``, ``y``) coordinates to define a second point
                on the line. Defaults to ``None``.
            vector: A tuple or numpy array of (``x``, ``y``) components to define the
                line direction. Defaults to ``None``.

        Returns:
            A tuple of lists containing Geometry objects that are bisected about
            the infinite line defined by the two given points. The first item in the
            tuple represents the geometries on the ``"top"`` of the line (or to the
            ``"right"`` of the line, if vertical) and the second item represents the
            geometries to the ``"bottom"`` of the line (or to the ``"left"`` of the
            line, if vertical).
        """
        top_geoms_acc: list[Geometry] = []
        bottom_geoms_acc: list[Geometry] = []

        for geom in self.geoms:
            top_geoms, bottom_geoms = geom.split_section(
                point_i=point_i, point_j=point_j, vector=vector
            )
            top_geoms_acc += top_geoms
            bottom_geoms_acc += bottom_geoms

        return (top_geoms_acc, bottom_geoms_acc)

    def offset_perimeter(
        self,
        amount: float = 0.0,
        where: str = "exterior",
        resolution: int = 12,
    ) -> CompoundGeometry:
        """Dilates/erodes the perimeter of a CompoundGeometry object by an amount.

        Args:
            amount: Distance to offset the section by. A negative value "erodes" the
                section. A positive value "dilates" the section. Defaults to ``0.0``.
            where: One of either ``"exterior"``, ``"interior"``, or ``"all"`` to specify
                which edges of the geometry to offset. If geometry has no interiors,
                then this parameter has no effect. Defaults to ``"exterior"``.
            resolution: Number of segments used to approximate a quarter circle around
                a point. Defaults to ``12``.

        Returns:
            CompoundGeometry object translated to new alignment

        Example:
            The following example erodes a 200UB25 with a 12 plate stiffener section by
            2 mm:

            .. plot::
                :include-source: True
                :caption: Eroded 200UB25 with stiffener.

                from sectionproperties.pre.library import rectangular_section, i_section

                geom1 = i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
                geom2 = rectangular_section(d=12, b=133)
                compound = geom1 + geom2.align_center(align_to=geom1).align_to(
                    other=geom1, on="top"
                )
                compound.offset_perimeter(amount=-2).plot_geometry()

        .. note::
            If performing a positive offset on a ``CompoundGeometry`` with multiple
            materials, ensure that the materials propagate as desired by performing a
            ``.plot_mesh()`` prior to performing any analysis.
        """
        if amount < 0:  # Eroding condition
            unionized_poly = unary_union([geom.geom for geom in self.geoms])

            if isinstance(unionized_poly, Polygon):
                offset_geom = Geometry(geom=unionized_poly).offset_perimeter(
                    amount=amount,
                    where=where,
                    resolution=resolution,
                )
            else:
                msg = "Geometry generation failed."
                raise RuntimeError(msg)

            # Using the offset_geom as a "mask"
            geoms_acc: list[Geometry] = []

            for geom in self.geoms:
                # Use symmetric intersection to find the region of the original
                # that intersects with the eroded unionized shape
                intersection_geom = geom & offset_geom

                if not intersection_geom.geom.is_empty:
                    geoms_acc.append(intersection_geom)

            return CompoundGeometry(geoms=geoms_acc)
        elif amount > 0:  # Ballooning condition
            # The algorithm used in the compound_dilation function cannot
            # currently handle re-entrant corners between different
            # geometries (a re-entrant corner in a single geometry is fine).
            # Re-entrant corners will require the creation of a new
            # "interface line" in the overlapping region created during
            # the dilation. I have a thought on how to do this but I just
            # have not gotten to it yet (note for later: it's like rotating
            # the overlap region between shear wall corners except cutting
            # across it from the shared vertex to the opposite vertex)
            # connorferster 2024-07-18

            return compound_dilation(self.geoms, offset=amount)
        else:
            return self

    def compile_geometry(self) -> None:
        """Creates geometry data from shapely data.

        Converts the ``shapely.Polygon`` objects stored in ``self.geoms`` into lists of
        points, facets, control_points, and hole points.
        """
        self.points: list[tuple[float, float]] = []
        self.facets: list[tuple[int, int]] = []
        self.control_points: list[tuple[float, float]] = []
        self.holes: list[tuple[float, float]] = []

        # loop through all sections
        for geom in self.geoms:
            if not all([geom.points, geom.facets, geom.control_points]):
                geom.create_facets_and_control_points()  # If not previously done

            # add points and skip duplicate points
            for point in geom.points:
                if point not in self.points:
                    self.points.append(point)

            # The facet numbering from the constituent Polygon is no longer valid
            # in the MultiPolygon.
            # We need to map the facets from the original Polygon points to the new
            # collected MultiPolygon points, which are in a different order.
            # Because points are not in their "original" order, we have to find the
            # matching point in the new self.points list and map the facet from the old
            # points list to the new self.points list. In other words, we need to change
            # the facet numbering so that connectivity is preserved but there are no
            # duplicate points.

            for facet in geom.facets:
                i_pnt, j_pnt = geom.points[facet[0]], geom.points[facet[1]]
                i_pnt_idx = self.points.index(i_pnt)  # using <list>.index method
                j_pnt_idx = self.points.index(j_pnt)  # using <list>.index method
                self.facets.append((i_pnt_idx, j_pnt_idx))

            # add control points
            for control_point in geom.control_points:
                self.control_points.append(control_point)

        # Determine if new holes have been created or if existing
        # holes have been destroyed (or "filled in").
        resultant_holes: list[tuple[float, float]] = []
        unionized_poly = unary_union([geom.geom for geom in self.geoms])
        buffer_amount = unionized_poly.area * SCALE_CONSTANT
        unionized_poly = unionized_poly.buffer(buffer_amount)

        if isinstance(unionized_poly, MultiPolygon):
            for poly in unionized_poly.geoms:
                for interior in poly.interiors:
                    rp: Point = Polygon(interior).representative_point()
                    coords = rp.coords[0]
                    resultant_holes.append(coords)  # pyright: ignore [reportArgumentType]

        elif isinstance(unionized_poly, Polygon):  # pyright: ignore [reportUnnecessaryIsInstance]
            if Geometry(unionized_poly).holes:
                resultant_holes = Geometry(geom=unionized_poly).holes
            else:  # Holes have been destroyed through operations
                resultant_holes = []

        self.holes = resultant_holes

    def calculate_perimeter(self) -> float:
        """Returns the length of the exterior perimeter of the CompoundGeometry.

        If the CompoundGeometry includes disjoint geometries then the perimeter cannot
        be calculated and the method returns -1.0.

        Returns:
            CompoundGeometry perimeter
        """
        unionized_poly = unary_union([geom.geom for geom in self.geoms])

        if isinstance(unionized_poly, Polygon):
            return float(unionized_poly.exterior.length)

        return -1.0


def load_dxf(
    dxf_filepath: str | pathlib.Path,
    spline_delta: float,
    degrees_per_segment: float,
) -> Geometry | CompoundGeometry:
    """Import any-old-shape in ``.dxf`` format for analysis.

    Args:
        dxf_filepath: Path to ``.dxf`` file to load
        spline_delta: Splines are not supported in ``shapely``, so they are approximated
            as polylines, this argument affects the spline sampling rate
        degrees_per_segment: The number of degrees discretised as a single line segment

    Raises:
        ImportError: If ``cad_to_shapely`` is not installed. To enable dxf features use
            ``pip install sectionproperties[dxf]``.
        RuntimeError: No polygon objects are found in the file
        ValueError: The filepath does not exist

    Returns:
        Geometry or CompoundGeometry loaded from ``.dxf`` file
    """
    try:
        import cad_to_shapely as c2s
    except ImportError as e:
        print(e)
        msg = "To use 'from_dxf(...)' you need to 'pip install cad_to_shapely'"
        raise ImportError(msg) from e

    if isinstance(dxf_filepath, str):
        dxf_filepath = pathlib.Path(dxf_filepath)

    if not dxf_filepath.exists():
        msg = f"The filepath does not exist: {dxf_filepath}"
        raise ValueError(msg)

    my_dxf = c2s.dxf.DxfImporter(str(dxf_filepath))
    my_dxf.process(spline_delta=spline_delta, degrees_per_segment=degrees_per_segment)
    my_dxf.cleanup()

    polygons = my_dxf.polygons
    new_polygons = c2s.utils.filter_polygons(polygons)

    # ensure list length > 0
    if len(new_polygons) == 0:
        msg = f"No shapely.Polygon objects found in file: {dxf_filepath}"
        raise RuntimeError(msg)

    # ensure only Polygons are generated
    for poly in new_polygons:
        if not isinstance(poly, Polygon):  # pyright: ignore [reportUnnecessaryIsInstance]
            msg = f"Not all objects found in file: {dxf_filepath} are Polygons"
            raise RuntimeError(msg)

    if len(new_polygons) == 1:
        return Geometry(new_polygons[0])
    else:
        return CompoundGeometry(MultiPolygon(new_polygons))


def create_facets(
    points_list: list[tuple[float, float]],
    connect_back: bool = False,
    offset: int = 0,
) -> list[tuple[int, int]]:
    """Create a list of facets.

    Returns a list of integer tuples representing the "facets" connecting the list
    of coordinates in ``points_list``. It is assumed that ``points_list`` coordinates
    are already in their order of connectivity.

    Args:
        points_list: List of coordinates
        connect_back: If True, then the last facet pair will be
            ``[len(points_list), offset]``. Defaults to ``False``.
        offset: An integer representing the value that the facets should begin
            incrementing from. Defaults to ``0``.

    Returns:
        List of facets
    """
    idx_peeker = more_itertools.peekable(
        [idx + offset for idx, _ in enumerate(points_list)]
    )
    facets = [(item, idx_peeker.peek(offset)) for item in idx_peeker]

    if connect_back:
        return facets

    return facets[:-1]


def create_exterior_points(
    shape: Polygon,
) -> list[tuple[float, float]]:
    """Creates a list of exterior points of a Polygon.

    Return a list of points tuples representing ``x``, ``y`` pairs of the exterior
    perimeter of ``polygon``.

    Args:
        shape: Shape to return exterior points

    Returns:
        List of exterior points
    """
    return cast(
        list[tuple[float, float]], [tuple(coord) for coord in shape.exterior.coords]
    )


def create_interior_points(
    lr: LinearRing,
) -> list[tuple[float, float]]:
    """Creates a list of inerior points of a Polygon.

    Return a list of point tuples representing ``x``, ``y`` pairs of the interior
    perimeter of ``lr``.

    Args:
        lr: LinearRing to return exterior points

    Returns:
        List of interior points
    """
    return cast(list[tuple[float, float]], [tuple(coord) for coord in lr.coords])


def create_points_and_facets(
    shape: Polygon,
    tol: int = 12,
) -> tuple[list[tuple[float, float]], list[tuple[int, int]]]:
    """Create the list of points and facets given a shapely Polygon.

    Args:
        shape: Polygon to create points and facets from
        tol: Number of decimal places to round the geometry vertices to. Defaults to
            ``12``,

    Returns:
        List of points (``x``, ``y``) and list of facets (``f1``, ``f2``)
    """
    master_count = 0
    points: list[tuple[float, float]] = []
    facets: list[tuple[int, int]] = []

    # Shape perimeter, note the last point == first point (shapely)
    for coords in list(shape.exterior.coords[:-1]):
        x, y = list(coords)
        x, y = round(x, tol), round(y, tol)
        points.append((x, y))
        master_count += 1

    facets += create_facets(points_list=points, connect_back=True)
    exterior_count = master_count

    # Holes
    for idx, hole in enumerate(shape.interiors):
        break_count = master_count
        int_points: list[tuple[float, float]] = []

        for coords in hole.coords[:-1]:  # The last point == first point (shapely)
            x, y = list(coords)
            x, y = round(x, tol), round(y, tol)
            int_points.append((x, y))
            master_count += 1

        # (idx > 0), (idx < 1) are like a 'step functions'
        offset = break_count * (idx > 0) + exterior_count * (idx < 1)
        facets += create_facets(
            points_list=int_points,
            connect_back=True,
            offset=offset,
        )
        points += int_points

    return points, facets


def buffer_polygon(
    polygon: Polygon,
    amount: float,
    resolution: int,
) -> Polygon | MultiPolygon:
    """Buffers (erode/corrode) a polygon by amount given a resolution.

    Args:
        polygon: Polygon to buffer
        amount: Amount to buffer
        resolution: Number of segments used to approximate a quarter circle around a
            point

    Returns:
        Buffered polygon
    """
    buffered_polygon = polygon.buffer(distance=amount, resolution=resolution)

    if isinstance(buffered_polygon, GeometryCollection):
        remaining_polygons = [
            item
            for item in buffered_polygon.geoms  # pyright: ignore [reportUnknownVariableType]
            if isinstance(item, Polygon)
        ]

        if len(remaining_polygons) == 1:
            return remaining_polygons[0]
        else:
            return MultiPolygon(remaining_polygons)

    elif isinstance(buffered_polygon, MultiPolygon | Polygon):  # pyright: ignore [reportUnnecessaryIsInstance]
        return buffered_polygon
    else:
        return Polygon()


def filter_non_polygons(
    input_geom: GeometryCollection | LineString | Point | Polygon | MultiPolygon,
) -> Polygon | MultiPolygon:
    """Filters shapely geometries to return only polygons.

    Returns a Polygon or a MultiPolygon representing any such Polygon on MultiPolygon
    that may exist in the ``input_geom``. If ``input_geom`` is a LineString or Point, an
    empty Polygon is returned.

    Args:
        input_geom: Shapely geometry to filter

    Returns:
        Filtered polygon
    """
    if isinstance(input_geom, Polygon | MultiPolygon):
        return input_geom
    elif isinstance(input_geom, GeometryCollection):
        acc = [
            item
            for item in input_geom.geoms
            if isinstance(item, MultiPolygon | Polygon)
        ]

        if len(acc) == 0:
            return Polygon()
        elif len(acc) == 1:
            return acc[0]
        else:
            return MultiPolygon(acc)  # pyright: ignore [reportArgumentType]
    else:
        return Polygon()


def round_polygon_vertices(
    poly: Polygon,
    tol: int,
) -> Polygon:
    """Returns a Polygon with all of its vertices rounded by ``tol``.

    Args:
        poly: Polygon to round
        tol: Number of decimals to round vertices to

    Returns:
        Polygon with rounded vertices
    """
    rounded_exterior = np.round(poly.exterior.coords, tol)
    rounded_interiors = [np.round(interior.coords, tol) for interior in poly.interiors]

    if not rounded_exterior.any():
        return Polygon()

    return Polygon(rounded_exterior, rounded_interiors)


def check_geometry_overlaps(
    lop: list[Polygon],
) -> bool:
    """Checks for overlaps in geometry.

    Returns True if any of the Polygon in the list of Polygons, ``lop``, are
    overlapping with any other Polygons in ``lop``. Returns False if all Polygon are
    disjoint.

    Args:
        lop: List of polygons

    Returns:
        Whether or not there is overlapping geometry
    """
    union_area = unary_union(lop).area
    sum_polygons = sum([poly.area for poly in lop])
    return not math.isclose(union_area, sum_polygons)


def compound_dilation(
    geoms: list[Geometry],
    offset: float,
) -> CompoundGeometry:
    """Returns a CompoundGeometry representing the input Geometries, dilated.

    Args:
        geoms: List of Geometry objects
        offset: A positive ``float`` or ``int``

    Raises:
        RuntimeError: If geometry generation fails

    Returns:
        The geometries dilated by ``offset``
    """
    polys = [geom.geom for geom in geoms]
    geom_network = build_geometry_network(polys)
    acc: list[Geometry] = []

    for poly_idx, connectivity in geom_network.items():
        poly_orig = polys[poly_idx]
        poly_orig_exterior = poly_orig.exterior
        connected_polys = [polys[idx].exterior for idx in connectivity]
        mucky_shared_paths = shared_paths(poly_orig_exterior, connected_polys)
        shared_path_geometries = MultiLineString(
            extract_shared_paths(mucky_shared_paths)  # pyright: ignore [reportArgumentType]
        )
        source = line_merge(poly_orig_exterior - shared_path_geometries)
        buff = source.buffer(offset, cap_style="flat")
        union = poly_orig | buff

        if isinstance(union, Polygon):
            new = Geometry(union, material=geoms[poly_idx].material)
        else:
            msg = "Geometry generation failed."
            raise RuntimeError(msg)

        acc.append(new)

    return CompoundGeometry(acc)


def check_geometry_disjoint(
    lop: list[Polygon],
) -> bool:
    """Returns True if any polygons in ``lop`` are disjoint, otherwise returns False.

    Args:
        lop: List of polygons

    Returns:
        Whether or not there is disjoint geometry
    """
    # Build polygon connectivity network
    network = build_geometry_network(lop)

    def walk_network(
        node: int,
        network: dict[int, set[int]],
        nodes_visited: list[int],
    ) -> list[int]:
        """Walks the network modifying 'nodes_visited' as it walks.

        Args:
            node: Initial node index
            network: Dictionary describing the node network
            nodes_visited: Initial list of nodes visited

        Returns:
            Updated node indices of the nodes that have been visisted
        """
        connections = network.get(node, set())

        for connection in connections:
            if connection in nodes_visited:
                continue
            else:
                nodes_visited.append(connection)
                walk_network(connection, network, nodes_visited)

        return nodes_visited

    # Traverse polygon connectivity network
    nodes_visited = [0]
    walk_network(0, network, nodes_visited)
    return set(nodes_visited) != set(network.keys())


def build_geometry_network(
    lop: list[Polygon],
) -> dict[int, set[int]]:
    """Builds a geometry connectivity graph.

    Returns a graph describing the connectivity of each polygon to each other polygon in
    ``lop``. The keys are the indexes of the polygons in ``lop`` and the values are a
    set of indexes that the key is connected to.

    Args:
        lop: List of Polygon

    Returns:
        A dictionary describing the connectivity graph of the polygons
    """
    network: dict[int, set[int]] = {}

    for idx_i, poly1 in enumerate(lop):
        for idx_j, poly2 in enumerate(lop):
            if idx_i != idx_j:
                connectivity = network.get(idx_i, set())

                if poly1.intersection(poly2):
                    connectivity.add(idx_j)

                network[idx_i] = connectivity

    return network


def extract_shared_paths(
    arr_of_geom_coll: list[GeometryCollection],
) -> list[LineString]:
    """Extracts a list of LineStrings exported by the shapely ``shared_paths`` method.

    Args:
        arr_of_geom_coll: An array of geometry collections

    Raises:
        RuntimeError: If linestrings cannot be extracted

    Returns:
        List of LineStrings
    """
    acc: list[LineString] = []

    for geom_coll in arr_of_geom_coll:
        for mls in geom_coll.geoms:
            if mls.is_empty:
                continue

            ls = line_merge(mls)

            if isinstance(ls, LineString):
                acc.append(ls)
            else:
                msg = "Extract shared paths failed."
                raise RuntimeError(msg)

    return acc