geometry.py

from __future__ import annotations
from ntpath import join
from typing import List, Optional, Union, Tuple, List, Any

import copy
import math
import pathlib
import more_itertools
import numpy as np
from shapely.geometry import (
    Polygon,
    MultiPolygon,
    LineString,
    LinearRing,
    Point,
    box,
    GeometryCollection,
)
from shapely.ops import split, unary_union
import shapely
import matplotlib
import matplotlib.pyplot as plt
import sectionproperties.pre.pre as pre
import sectionproperties.pre.bisect_section as bisect
import sectionproperties.post.post as post


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.

    :cvar geom: a Polygon object that defines the geometry
    :vartype geom: :class:`shapely.geometry.Polygon`

    :cvar material: Optional, a Material to associate with this geometry
    :vartype material: Optional[:class:`~sectionproperties.pre.Material`]

    :cvar control_point: Optional, 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 :func:`shapely.geometry.Polygon.representative_point`.

    :cvar tol: Optional, default is 12. 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.
    """

    def __init__(
        self,
        geom: shapely.geometry.Polygon,
        material: pre.Material = pre.DEFAULT_MATERIAL,
        control_points: Optional[Union[Point, List[float, float]]] = None,
        tol=12,
    ):
        """Inits the Geometry class."""
        if isinstance(geom, MultiPolygon):
            raise ValueError(f"Use CompoundGeometry(...) for a MultiPolygon object.")
        if not isinstance(geom, Polygon):
            raise ValueError(
                f"Argument is not a valid shapely.geometry.Polygon object: {repr(geom)}"
            )
        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  # Represents 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 = []
        self.shift = []
        self.points = []
        self.facets = []
        self.holes = []
        self.perimeter = []
        self._recovery_points = []

        self.compile_geometry()

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

    def assign_control_point(self, control_point: List[float, float]):
        """
        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 shapely.geometry.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.

        :cvar control_points: An *(x, y)* coordinate that describes the distinct, contiguous,
            region of a single material within the geometry.
            Exactly one point is required for each geometry with a distinct material.
        :vartype control_point: list[float, float]

        """
        return Geometry(self.geom, self.material, control_point)

    @staticmethod
    def from_points(
        points: List[List[float]],
        facets: List[List[int]],
        control_points: List[List[float]],
        holes: Optional[List[List[float]]] = None,
        material: Optional[pre.Material] = pre.DEFAULT_MATERIAL,
    ):
        """
        An interface for the creation of Geometry objects through the definition of points,
        facets, and holes.

        :cvar points: List of points *(x, y)* defining the vertices of the section geometry.
            If facets are not provided, it is a assumed the that the list of points are ordered
            around the perimeter, either clockwise or anti-clockwise.
        :vartype points: list[list[float, float]]
        :cvar facets: A list of *(start, end)* indexes 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.
        :vartype facets: list[list[int, int]]
        :cvar 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 CompoundGeometry.from_points()
        :vartype control_point: list[float, float]
        :cvar holes: Optional. 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.
        :vartype holes: list[list[float, float]]
        :cvar material: Optional. A :class:`~sectionproperties.pre.pre.Material` object
            that is to be assigned. If not given, then the
            :class:`~sectionproperties.pre.pre.DEFAULT_MATERIAL` will be used.
        :vartype materials: :class:`~sectionproperties.pre.pre.Material`
        """
        if len(control_points) != 1:
            raise ValueError(
                "Control points for Geometry instances must have exactly "
                "one x, y coordinate and entered as a list of list of float, e.g. [[0.1, 3.4]]."
                "CompoundGeometry.from_points() can accept multiple control points\n"
                f"Control points received: {control_points}"
            )
        if holes is None:
            holes = []
        prev_facet = []

        # Initialize the total number of accumulators needed
        # Always an exterior, plus, a separate accumulator for each interior region
        exterior = []
        interiors = [[] 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
                active_list.append(points[i_idx])
                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 (
                i_idx != prev_j_idx and holes
            ):  # If there is a break in the chain of edges...
                if (
                    active_list == exterior
                ):  # ...and we are still accumulating on the exterior...
                    active_list = interiors[
                        interior_counter
                    ]  # ... then move to the interior accumulator
                else:  # ...or if we are already in the interior accumulator...
                    interior_counter += (
                        1  # ...then start the next interior accumulator for a new hole.
                    )
                    active_list = interiors[interior_counter]
                active_list.append(points[i_idx])
            else:
                active_list.append(
                    points[i_idx]
                )  # Only need i_idx b/c shapely auto-closes polygons
            prev_facet = facet

        exterior_geometry = Polygon(exterior)
        interior_polys = [Polygon(interior) for interior in interiors]
        interior_geometry = MultiPolygon(interior_polys)
        geometry = Geometry(
            exterior_geometry - interior_geometry, material, control_points
        )
        return geometry

    @staticmethod
    def from_dxf(
        dxf_filepath: Union[str, pathlib.Path],
    ) -> Union[Geometry, CompoundGeometry]:
        """
        An interface for the creation of Geometry objects from CAD .dxf files.

        :cvar dxf_filepath: A path-like object for the dxf file
        :vartype dxf_filepath: Union[str, pathlib.Path]
        """
        return load_dxf(dxf_filepath)

    @classmethod
    def from_3dm(cls, filepath: Union[str, pathlib.Path], **kwargs) -> Geometry:
        """Class method to create a `Geometry` from the objects in a Rhino `.3dm` file.

        :param filepath:
            File path to the rhino `.3dm` file.
        :type filepath: Union[str, pathlib.Path]
        :param kwargs:
            See below.
        :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.
        :return:
            A Geometry object.
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

        :Keyword Arguments:
            * *refine_num* (``int, optional``) --
                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* (``numpy.ndarray, optional``) --
                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* (``numpy.ndarray, optional``) --
                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* (``float, optional``) --
                The distance to the Shapely plane.
                Default is 0.
            * *project* (``boolean, optional``) --
                Controls if the breps are projected onto the plane in the direction of the Shapley plane's normal.
                Default is True.
            * *parallel* (``boolean, optional``) --
                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.
        """
        try:
            import sectionproperties.pre.rhino as rhino_importer  # type: ignore
        except ImportError as e:
            print(e)
            print(
                "There is something wrong with your rhino library installation. "
                "Please report this error at https://github.com/robbievanleeuwen/section-properties/issues"
            )
            return

        geom = None
        list_poly = rhino_importer.load_3dm(filepath, **kwargs)
        if len(list_poly) == 1:
            geom = cls(geom=list_poly[0])
        else:
            raise RuntimeError(
                f"Multiple surfaces extracted from the file. "
                f"Either use CompoundGeometry or extract individual surfaces manually via pre.rhino."
            )
        return geom

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

        :param r3dm_brep:
            A Rhino3dm.Brep encoded as a string.
        :type r3dm_brep: str
        :param kwargs:
            See below.
        :return:
            A Geometry object found in the encoded string.
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

        :Keyword Arguments:
            * *refine_num* (``int, optional``) --
                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* (``numpy.ndarray, optional``) --
                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* (``numpy.ndarray, optional``) --
                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* (``float, optional``) --
                The distance to the Shapely plane.
                Default is 0.
            * *project* (``boolean, optional``) --
                Controls if the breps are projected onto the plane in the direction of the Shapley plane's normal.
                Default is True.
            * *parallel* (``boolean, optional``) --
                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.
        """
        try:
            import sectionproperties.pre.rhino as rhino_importer  # type: ignore
        except ImportError as e:
            print(e)
            print(
                "There is something wrong with your rhino library installation. "
                "Please report this error at https://github.com/robbievanleeuwen/section-properties/issues"
            )
            return
        return cls(geom=rhino_importer.load_brep_encoding(r3dm_brep, **kwargs)[0])

    def create_facets_and_control_points(self):
        self.perimeter = None
        self.perimeter = list(range(len(self.geom.exterior.coords)))
        self.holes = []
        self.points = []
        self.facets = []
        self.points, self.facets = create_points_and_facets(self.geom, self.tol)
        if not self.assigned_control_point:
            self.control_points = tuple(self.geom.representative_point().coords)
        else:
            self.control_points = tuple(self.assigned_control_point.coords)

        for hole in self.geom.interiors:
            hole_polygon = Polygon(hole)
            self.holes += tuple(hole_polygon.representative_point().coords)
        return

    def compile_geometry(self):  # Alias
        """
        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: Union[float, List[float]], coarse: bool = False):
        """Creates a quadratic triangular mesh from the Geometry object.

        :param mesh_sizes: A float describing the maximum mesh element area to be used
            within the Geometry-object finite-element mesh.
        :type mesh_sizes: Union[float, List[float]]
        :param bool coarse: If set to True, will create a coarse mesh (no area or
            quality constraints)

        :return: Geometry-object with mesh data stored in .mesh attribute. Returned
            Geometry-object is self, not a new instance.
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

        The following example creates a circular cross-section with a diameter of 50 with 64
        points, and generates a mesh with a maximum triangular area of 2.5::

            import sectionproperties.pre.library.primitive_sections as primitive_sections

            geometry = primitive_sections.circular_section(d=50, n=64)
            geometry = geometry.create_mesh(mesh_sizes=2.5)

        ..  figure:: ../images/sections/circle_mesh.png
            :align: center
            :scale: 75 %

            Mesh generated from the above geometry.
        """
        if isinstance(mesh_sizes, (list, tuple)) and len(mesh_sizes) == 1:
            mesh_size = mesh_sizes[0]
        elif isinstance(mesh_sizes, (float, int)):
            mesh_size = [mesh_sizes]
        else:
            raise ValueError(
                "Argument 'mesh_sizes' for a Geometry must be either "
                f"a float, or a list of float with length of 1, not {mesh_sizes}."
            )
        self.mesh = pre.create_mesh(
            self.points, self.facets, self.holes, self.control_points, mesh_size, coarse
        )
        return self

    def align_to(
        self,
        other: Union[Geometry, Tuple[float, float]],
        on: str,
        inner: bool = False,
    ) -> Geometry:
        """
        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.

        :param other: Either another Geometry or a tuple representing an
            *(x,y)* coordinate point that 'self' should align to.
        :type other: Union[:class:`~sectionproperties.pre.geometry.Geometry`, Tuple[float, float]]

        :param on: A str of either "left", "right", "bottom", or "top" indicating which
            side of 'other' that self should be aligned to.

        :param inner: Default False. 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.
        :type inner: bool

        :return: Geometry object translated to alignment location
        :rtype: :class:`~sectionproperties.pre.geometry.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]
        align_to_geometry = isinstance(other, (Geometry, CompoundGeometry))
        if align_to_geometry:
            align_to_idx = other_as_geom_map[on]
        else:
            align_to_idx = other_as_point_map[on]

        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]
        if align_to_geometry:
            other_extents = other.calculate_extents()
        else:
            other_extents = other
        align_to_coord = other_extents[align_to_idx]

        offset = align_to_coord - self_align_coord
        # offset = align_to_coord - self_align_coord
        arg = "x_offset"
        if on in ["top", "bottom"]:
            arg = "y_offset"
        kwargs = {arg: offset}
        new_geom = self.shift_section(**kwargs)

        return new_geom

    def align_center(
        self, align_to: Optional[Union[Geometry, Tuple[float, float]]] = None
    ):
        """
        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.

        :param align_to: Another Geometry to align to or None (default is None)
        :type align_to: Optional[Union[:class:`~sectionproperties.pre.geometry.Geometry`, Tuple[float, float]]]

        :return: Geometry object translated to new alignment
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`
        """
        cx, cy = list(self.geom.centroid.coords)[0]
        # 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 = list(align_to.geom.centroid.coords)[0]
            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,
            ):  # align_to not subscriptable, incorrect length, etc.
                raise ValueError(
                    f"align_to must be either a Geometry object or an x, y coordinate, not {align_to}."
                )
        new_geom = self.shift_section(x_offset=shift_x, y_offset=shift_y)
        return new_geom

    def shift_section(
        self,
        x_offset=0.0,
        y_offset=0.0,
    ):
        """
        Returns a new Geometry object translated by 'x_offset' and 'y_offset'.

        :param x_offset: Distance in x-direction by which to shift the geometry.
        :type x_offset: float
        :param y_offset: Distance in y-direction by which to shift the geometry.
        :type y_offset: float

        :return: New Geometry-object shifted by 'x_offset' and 'y_offset'
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`
        """
        # Move assigned control point
        new_ctrl_point = None
        if self.assigned_control_point:
            new_ctrl_point = shapely.affinity.translate(
                self.assigned_control_point, x_offset, y_offset
            ).coords[0]
        new_geom = Geometry(
            shapely.affinity.translate(self.geom, x_offset, y_offset),
            self.material,
            new_ctrl_point,
        )
        return new_geom

    def rotate_section(
        self,
        angle: float,
        rot_point: Union[List[float], str] = "center",
        use_radians: bool = False,
    ):
        """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.

        :param float angle: Angle (degrees by default) by which to rotate the section. A positive angle leads
            to a counter-clockwise rotation.
        :param rot_point: Optional. Point *(x, y)* about which to rotate the section. If not provided, will rotate
            about the center of the geometry's bounding box. Default = 'center'.
        :type rot_point: list[float, float]
        :param use_radians: Boolean to indicate whether 'angle' is in degrees or radians. If True, 'angle' is interpreted as radians.

        :return: New Geometry-object rotated by 'angle' about 'rot_point'
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

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

            import sectionproperties.pre.library.steel_sections as steel_sections

            geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
            new_geometry = geometry.rotate_section(angle=-30)
        """
        new_ctrl_point = None
        if self.assigned_control_point:
            rotate_point = rot_point
            if rot_point == "center":
                rotate_point = box(*self.geom.bounds).centroid
            new_ctrl_point = shapely.affinity.rotate(
                self.assigned_control_point, angle, rotate_point, use_radians
            ).coords[0]
        new_geom = Geometry(
            shapely.affinity.rotate(self.geom, angle, rot_point, use_radians),
            self.material,
            new_ctrl_point,
        )
        return new_geom

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

        :param string axis: Axis about which to mirror the geometry, *'x'* or *'y'*
        :param mirror_point: Point about which to mirror the geometry *(x, y)*.
            If no point is provided, mirrors the geometry about the centroid of the shape's bounding box.
            Default = 'center'.
        :type mirror_point: Union[list[float, float], str]

        :return: New Geometry-object mirrored on 'axis' about 'mirror_point'
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

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

            import sectionproperties.pre.library.steel_sections as steel_sections

            geometry = steel_sections.channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
            new_geometry = geometry.mirror_section(axis='y', mirror_point=[0, 0])
        """
        x_mirror = 1
        y_mirror = 1
        if mirror_point != "center":
            x, y = mirror_point
            mirror_point = (x, y, 0)
        if axis == "x":
            x_mirror = -x_mirror
        elif axis == "y":
            y_mirror = -y_mirror
        mirrored_geom = shapely.affinity.scale(
            self.geom, xfact=y_mirror, yfact=x_mirror, zfact=1.0, origin=mirror_point
        )

        new_ctrl_point = None
        if self.assigned_control_point:
            new_ctrl_point = shapely.affinity.scale(
                self.assigned_control_point,
                xfact=y_mirror,
                yfact=x_mirror,
                zfact=1.0,
                origin=mirror_point,
            ).coords[0]
        new_geom = Geometry(mirrored_geom, self.material, new_ctrl_point)
        return new_geom

    def split_section(
        self,
        point_i: Tuple[float, float],
        point_j: Optional[Tuple[float, float]] = None,
        vector: Union[Optional[Tuple[float, float]], np.ndarray] = None,
    ) -> Tuple[List[Geometry], List[Geometry]]:
        """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.

        :param point_i: A tuple of *(x, y)* coordinates to define a first point on the line
        :type point_i: Tuple[float, float]

        :param point_j: Optional. A tuple of *(x, y)* coordinates to define a second point on the line
        :type point_j: Tuple[float, float]

        :param vector: Optional. A tuple or numpy ndarray of *(x, y)* components to define the line direction.
        :type vector: Union[Tuple[float, float], numpy.ndarray]

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

        :rtype: Tuple[List[Geometry], List[Geometry]]

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

            import sectionproperties.pre.library.steel_sections as steel_sections
            from shapely.geometry import LineString

            geometry = steel_sections.channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
            right_geom, left_geom = geometry.split_section((0, 0), (0, 1))
        """
        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)
        elif not point_j and not vector:
            raise ValueError(
                "Either a second point or a vector must be given to define the line."
            )
        bounds = self.calculate_extents()
        line_segment = bisect.create_line_segment(point_i, vector, bounds)
        top_right_polys, bottom_left_polys = bisect.group_top_and_bottom_polys(
            split(self.geom, line_segment), line_segment
        )

        # Create new Geometry instances 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, where: str = "exterior", resolution: float = 12
    ):
        """Dilates or erodes the section perimeter by a discrete amount.

        :param amount: Distance to offset the section by. A -ve value "erodes" the section.
            A +ve value "dilates" the section.
        :type amount: float
        :param 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.
            Default is "exterior".
        :type where: str
        :param resolution: Number of segments used to approximate a quarter circle around a point
        :type resolution: float

        :return: Geometry object translated to new alignment
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

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

            import sectionproperties.pre.library.steel_sections as steel_sections

            geometry = sections.channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
            new_geometry = geometry.offset_perimeter(amount=-2)
        """
        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(
                    [Geometry(poly, self.material) for poly in exterior_polygon]
                )

            # 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(exterior_polygon, self.material, self.control_points[0])
            return Geometry(exterior_polygon, self.material)

        elif not self.geom.interiors and where == "interior":
            raise ValueError(
                "Cannot buffer interior of Geometry object if it has no holes."
            )

        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(
                    [Geometry(poly, self.material) for poly in buffered_exterior.geoms]
                )

            # 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(
                    buffered_exterior, self.material, self.control_points[0]
                )
            return Geometry(buffered_exterior, self.material)

        elif where == "all":
            buffered_geom = buffer_polygon(self.geom, amount, resolution)
            if isinstance(buffered_geom, MultiPolygon):
                compound_geom = CompoundGeometry(
                    [Geometry(poly, self.material) for poly in buffered_geom]
                )
                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
            ):
                single_geom = Geometry(
                    buffered_geom, self.material, self.control_points[0]
                )
                return single_geom
            return Geometry(buffered_geom, self.material)

    def shift_points(
        self,
        point_idxs: Union[int, List[int]],
        dx: float = 0,
        dy: float = 0,
        abs_x: Optional[float] = None,
        abs_y: Optional[float] = None,
    ) -> 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.

        :param point_idxs: An integer representing an index location or a list of integer
            index locations.
        :type point_idxs: Union[int, List[int]]
        :param dx: The number of units in the x-direction to shift the point(s) by
        :type dx: float
        :param dy: The number of units in the y-direction to shift the point(s) by
        :type dy: float
        :param 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.
        :type abs_x: Optional[float]
        :param 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.
        :type abs_y: Optional[float]

        :return: Geometry object with selected points translated to the new location.
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

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

            import sectionproperties.pre.library.primitive_sections as primitive_sections

            geometry = primitive_sections.rectangular_section(d=200, b=150)

            # Using relative shifting
            one_pt_shifted_geom = geometry.shift_points(point_idxs=1, dx=50)

            # Using absolute relocation
            both_pts_shift_geom = one_pt_shift_geom.shift_points(point_idxs=2, abs_x=200)
        """
        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(
            current_points,
            current_facets,
            current_control_points,
            current_holes,
            material=current_material,
        )
        if self.assigned_control_point and new_geom.geom.contains(
            self.assigned_control_point
        ):
            new_geom = Geometry.from_points(
                current_points,
                current_facets,
                current_control_points[0],
                current_holes,
                current_material,
            )
        return new_geom

    def plot_geometry(
        self,
        labels=["control_points"],
        title="Cross-Section Geometry",
        cp=True,
        legend=True,
        **kwargs,
    ):
        """Plots the geometry defined by the input section.

        :param labels: A list of str which indicate which labels to plot. Can be one
            or a combination of "points", "facets", "control_points", or an empty list
            to indicate no labels. Default is ["control_points"]
        :type labels: list[str]
        :param string title: Plot title
        :param bool cp: If set to True, plots the control points
        :param bool legend: If set to True, plots the legend
        :param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`

        :return: Matplotlib axes object
        :rtype: :class:`matplotlib.axes`

        The following example creates a CHS discretised with 64 points, with a diameter of 48 and
        thickness of 3.2, and plots the geometry::

            import sectionproperties.pre.library.steel_sections as steel_sections

            geometry = steel_sections.circular_hollow_section(d=48, t=3.2, n=64)
            geometry.plot_geometry()

        ..  figure:: ../images/sections/chs_geometry.png
            :align: center
            :scale: 75 %

            Geometry generated by the above example.
        """

        # create plot and setup the plot
        with post.plotting_context(title=title, **kwargs) as (fig, ax):
            # plot the points and facets
            for (i, f) in enumerate(self.facets):
                if i == 0:
                    label = "Points & Facets"
                else:
                    label = None

                ax.plot(
                    [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):
                if i == 0:
                    label = "Holes"
                else:
                    label = None

                ax.plot(h[0], h[1], "rx", markersize=5, markeredgewidth=1, label=label)

            if cp:
                # plot the control points
                for (i, cp) in enumerate(self.control_points):
                    if i == 0:
                        label = "Control Points"
                    else:
                        label = None

                    ax.plot(cp[0], cp[1], "bo", markersize=5, label=label)

            # 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")

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

                # 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")

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

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

        return ax

    def calculate_extents(self):
        """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.

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

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

        :return: Geometry perimeter.
        :rtype: float
        """
        perimeter = self.geom.exterior.length
        return perimeter

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

        :return: Geometry area.
        :rtype: float
        """
        area = self.geom.area
        return area

    def calculate_centroid(self):
        """Calculates the centroid of the geometry as a tuple
        of *(x,y)* coordinates.

        :return: Geometry centroid.
        :rtype: Tuple[float, float]
        """
        cx, cy = self.geom.centroid.coords[0]  # Nested list
        return (cx, cy)

    @property
    def recovery_points(self):
        """
        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.
        """
        return self._recovery_points

    @recovery_points.setter
    def recovery_points(
        self, new_points: Union[List[list], List[tuple], List[Point]]
    ) -> list:
        self._recovery_points = new_points

    @recovery_points.getter
    def recovery_points(
        self, new_points: Union[List[list], List[tuple], List[Point]]
    ) -> list:
        return self._recovery_points

    def __or__(self, other):
        """
        Perform union on Geometry objects with the | operator
        """
        material = self.material or other.material
        try:
            new_polygon = filter_non_polygons(self.geom | other.geom)
            if isinstance(new_polygon, MultiPolygon):

                return CompoundGeometry(
                    [Geometry(polygon, material) for polygon in new_polygon.geoms]
                )
            return Geometry(new_polygon, material, self.control_points[0])

        except:
            raise ValueError(
                f"Cannot perform 'union' on these two objects: {self} | {other}"
            )

    def __xor__(self, other):
        """
        Perform symmetric difference on Geometry objects with the ^ operator
        """
        material = self.material or other.material
        try:
            new_polygon = filter_non_polygons(self.geom ^ other.geom)
            if isinstance(new_polygon, MultiPolygon):
                return CompoundGeometry(
                    [Geometry(polygon, material) for polygon in new_polygon.geoms]
                )
            if isinstance(new_polygon, GeometryCollection):
                return None
            return Geometry(new_polygon, material)
        except:
            raise ValueError(
                f"Cannot perform 'symmetric difference' on these two objects: {self} ^ {other}"
            )

    def __sub__(self, other):
        """
        Perform difference on Geometry objects with the - operator
        """
        if isinstance(self, CompoundGeometry):
            subs_geom_acc = []
            for geom in self.geoms:
                new_geom = geom - other
                subs_geom_acc.append(new_geom)
            return CompoundGeometry(subs_geom_acc)
        else:
            material = self.material or pre.DEFAULT_MATERIAL
            try:
                new_polygon = filter_non_polygons(self.geom - other.geom)
                if isinstance(new_polygon, GeometryCollection):  # Non-polygon results
                    return None
                elif isinstance(new_polygon, MultiPolygon):
                    return CompoundGeometry(
                        [Geometry(polygon, material) for polygon in new_polygon.geoms]
                    )
                elif isinstance(new_polygon, Polygon):
                    if self.assigned_control_point and new_polygon.contains(
                        self.assigned_control_point
                    ):
                        return Geometry(
                            new_polygon, material, self.assigned_control_point
                        )
                    else:
                        return Geometry(new_polygon, material)
            except:
                raise ValueError(
                    f"Cannot perform 'difference' on these two objects: {self} - {other}"
                )

        # material = self.material or other.material
        # try:
        #     new_polygon = filter_non_polygons(self.geom - other.geom)
        #     if isinstance(new_polygon, MultiPolygon):
        #         return CompoundGeometry(
        #             [Geometry(polygon, material) for polygon in new_polygon.geoms]
        #         )
        #     if isinstance(new_polygon, GeometryCollection):
        #         return None
        #     # Check to see if assigned_control_point is still valid
        #     if self.assigned_control_point and new_polygon.contains(
        #         self.assigned_control_point
        #     ):
        #         return Geometry(new_polygon, material, self.control_points[0])
        #     return Geometry(new_polygon, material)
        # except:
        #     raise ValueError(
        #         f"Cannot perform 'difference' on these two objects: {self} - {other}"
        #     )

    def __add__(self, other):
        """
        Combine Geometry objects into a CompoundGeometry using the + operator
        """
        try:
            return CompoundGeometry([self, other])
        except:
            raise ValueError(
                f"Cannot create new CompoundGeometry with these objects: {self} + {other}"
            )

    def __and__(self, other):
        """
        Perform intersection on Geometry objects with the & operator
        """
        material = self.material or other.material
        try:
            new_polygon = filter_non_polygons(self.geom & other.geom)
            if isinstance(new_polygon, MultiPolygon) and len(new_polygon.geoms) > 1:
                return CompoundGeometry(
                    [Geometry(polygon, material) for polygon in new_polygon.geoms]
                )
            return Geometry(new_polygon, material)
        except:
            raise ValueError(
                f"Cannot perform 'intersection' on these two Geometry instances: {self} & {other}"
            )


class CompoundGeometry(Geometry):
    """Class for defining a geometry of multiple distinct regions, each potentially
    having different material properties.

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

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

    :cvar geoms: either a list of Geometry objects or a :class:`shapely.geometry.MultiPolygon`
        instance.
    :vartype geoms: Union[:class:`shapely.geometry.MultiPolygon`,
        List[:class:`~sectionproperties.pre.geometry.Geometry`]]
    """

    def __init__(self, geoms: Union[MultiPolygon, List[Geometry]]):
        if isinstance(geoms, MultiPolygon):
            self.geoms = [
                Geometry(poly, material=pre.DEFAULT_MATERIAL) for poly in geoms.geoms
            ]
            self.geom = geoms
        elif isinstance(geoms, list):
            processed_geoms = []
            for item in geoms:
                if isinstance(item, CompoundGeometry):
                    # "Flatten" any CompoundGeometry objects in the 'geoms' list
                    processed_geoms += item.geoms
                elif isinstance(item, Geometry):
                    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.perimeter = []
        self.material = None
        self.compile_geometry()
        self.tol = 12
        self.mesh = None

    def _repr_svg_(self):
        """
        Returns an svg representation of the CompoundGeometry.
        Wraps :func:`shapely.geometry.MultiPolygon._repr_svg_()` by returning
        ``self.geom._repr_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 self.geom._repr_svg_()

    @staticmethod
    def from_points(
        points: List[List[float]],
        facets: List[List[int]],
        control_points: List[List[float]],
        holes: Optional[List[List[float]]] = None,
        materials: Optional[List[pre.Material]] = pre.DEFAULT_MATERIAL,
    ):
        """
        An interface for the creation of CompoundGeometry objects through the definition of points,
        facets, holes, and control_points. 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).

        :cvar points: List of points *(x, y)* defining the vertices of the section geometry.
            If facets are not provided, it is a assumed the that the list of points are ordered
            around the perimeter, either clockwise or anti-clockwise
        :vartype points: list[list[float]]
        :cvar facets: A list of *(start, end)* indexes 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.
        :vartype facets: list[list[int]]
        :cvar control_points: Optional. A list of points *(x, y)* that define non-interior 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'.
            If not given, then points will be assigned automatically using
            shapely.geometry.Polygon.representative_point()
        :vartype control_points: list[list[float]]
        :cvar holes: Optional. 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.
        :vartype holes: list[list[float]]
        :cvar materials: Optional. 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.pre.DEFAULT_MATERIAL` will be used for each region.
        :vartype materials: list[:class:`~sectionproperties.pre.pre.Material`]
        """
        if materials and not control_points:
            raise ValueError(
                "Materials cannot be assigned without control_points. "
                "Please provide corresponding control_points for each material."
            )
        if holes is None:
            holes = list()
        if materials is not pre.DEFAULT_MATERIAL:
            if len(materials) != len(control_points):
                raise ValueError(
                    f"If materials are provided, the number of materials in the list must "
                    "match the number of control_points provided.\n"
                    f"len(materials)=={len(materials)}, len(control_points)=={len(control_points)}."
                )

        # First, generate all invidual polygons from points and facets
        current_polygon_points = []
        all_polygons = []
        prev_facet = None
        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
        else:  # Use the for...else clause to add the last point and close the last polygon.
            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 = []
        interiors = []

        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):
            raise ValueError(
                f"The number of exterior regions ({len(exteriors)}) "
                f"does not match the number of control_points given ({len(control_points)})."
            )
        if not interiors:
            if materials is pre.DEFAULT_MATERIAL:
                return CompoundGeometry(
                    [
                        Geometry(
                            exterior,
                            control_points=control_points[idx],
                            material=materials,
                        )
                        for idx, exterior in enumerate(exteriors)
                    ]
                )
            else:
                return CompoundGeometry(
                    [
                        Geometry(
                            exterior,
                            control_points=control_points[idx],
                            material=materials[idx],
                        )
                        for idx, exterior in enumerate(exteriors)
                    ]
                )

        else:
            # "Punch" all holes through each exterior geometry
            punched_exteriors = []
            punched_exterior_geometries = []
            for idx, exterior in enumerate(exteriors):
                punched_exterior = exterior
                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:
                        raise ValueError(
                            f"Control points given are not contained within the geometry"
                            f" once holes are subtracted: {control_points}"
                        )
                if materials is pre.DEFAULT_MATERIAL:

                    exterior_geometry = Geometry(
                        punched_exterior,
                        control_points=exterior_control_point,
                        material=materials,
                    )
                    punched_exterior_geometries.append(exterior_geometry)

                else:

                    exterior_geometry = Geometry(
                        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: Union[str, pathlib.Path], **kwargs) -> CompoundGeometry:
        """Class method to create a `CompoundGeometry` from the objects in a Rhino `3dm` file.

        :param filepath:
            File path to the rhino `.3dm` file.
        :type filepath: Union[str, pathlib.Path]
        :param kwargs:
            See below.
        :return:
            A `CompoundGeometry` object.
        :rtype: :class:`~sectionproperties.pre.geometry.CompoundGeometry`

        :Keyword Arguments:
            * *refine_num* (``int, optional``) --
                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* (``numpy.ndarray, optional``) --
                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* (``numpy.ndarray, optional``) --
                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* (``float, optional``) --
                The distance to the Shapely plane.
                Default is 0.
            * *project* (``boolean, optional``) --
                Controls if the breps are projected onto the plane in the direction of the Shapley plane's normal.
                Default is True.
            * *parallel* (``boolean, optional``) --
                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.
        """
        try:
            import sectionproperties.pre.rhino as rhino_importer  # type: ignore
        except ImportError as e:
            print(e)
            print(
                "There is something wrong with your rhino library installation. "
                "Please report this error at https://github.com/robbievanleeuwen/section-properties/issues"
            )
            return
        list_poly = rhino_importer.load_3dm(filepath, **kwargs)
        return cls(geoms=MultiPolygon(list_poly))

    def create_mesh(self, mesh_sizes: List[float], coarse: bool = False):
        """Creates a quadratic triangular mesh from the Geometry object.

        :param mesh_size: 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 is passed, then the one
            size will be applied to all constituent Geometry meshes.
        :type mesh_sizes: List[float]
        :param bool coarse: If set to True, will create a coarse mesh (no area or
            quality constraints)

        :return: Geometry-object with mesh data stored in .mesh attribute. Returned
            Geometry-object is self, not a new instance.
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

        The following example creates a circular cross-section with a diameter of 50 with 64
        points, and generates a mesh with a maximum triangular area of 2.5::

            import sectionproperties.pre.library.primitive_sections as primitive_sections

            geometry = primitive_sections.circular_section(d=50, n=64)
            geometry = geometry.create_mesh(mesh_sizes=[2.5])

        ..  figure:: ../images/sections/circle_mesh.png
            :align: center
            :scale: 75 %

            Mesh generated from the above geometry.
        """
        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(
            self.points,
            self.facets,
            self.holes,
            self.control_points,
            mesh_sizes,
            coarse,
        )
        return self

    def shift_section(self, x_offset: float = 0, y_offset: float = 0):
        """
        Returns a new CompoundGeometry object translated by 'x_offset' and 'y_offset'.

        :param x_offset: Distance in x-direction by which to shift the geometry.
        :type x_offset: float
        :param y_offset: Distance in y-direction by which to shift the geometry.
        :type y_offset: float

        :return: CompoundGeometry object shifted by 'x_offset' and 'y_offset'
        :rtype: :class:`~sectionproperties.pre.geometry.CompoundGeometry`
        """
        geoms_acc = []
        for geom in self.geoms:
            geoms_acc.append(geom.shift_section(x_offset=x_offset, y_offset=y_offset))
        new_geom = CompoundGeometry(geoms_acc)
        return new_geom

    def rotate_section(
        self,
        angle: float,
        rot_point: Union[List[float], str] = "center",
        use_radians: bool = False,
    ):
        """Rotates the compound geometry and specified angle about a point. If the rotation point is not
        provided, rotates the section about the center of the compound geometry's bounding box.

        :param float angle: Angle (degrees by default) by which to rotate the section. A positive angle leads
            to a counter-clockwise rotation.
        :param rot_point: Optional. Point *(x, y)* about which to rotate the section. If not provided, will rotate
            about the center of the compound geometry's bounding box. Default = 'center'.
        :type rot_point: list[float, float]
        :param use_radians: Boolean to indicate whether 'angle' is in degrees or radians. If True, 'angle' is interpreted as radians.

        :return: CompoundGeometry object rotated by 'angle' about 'rot_point'
        :rtype: :class:`~sectionproperties.pre.geometry.CompoundGeometry`

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

            import sectionproperties.pre.library.steel_sections as steel_sections
            import sectionproperties.pre.library.primitive_sections as primitive_sections

            geometry_1 = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
            geometry_2 = primitive_sections.rectangular_section(d=20, b=133)
            compound = geometry_2.align_center(geometry_1).align_to(geometry_1, on="top") + geometry_1
            new_compound = compound.rotate_section(angle=-30)
        """
        geoms_acc = []
        if rot_point == "center":
            rot_point = box(*MultiPolygon(self.geom).bounds).centroid
        for geom in self.geoms:
            geoms_acc.append(geom.rotate_section(angle, rot_point, use_radians))
        new_geom = CompoundGeometry(geoms_acc)
        return new_geom

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

        :param string axis: Axis about which to mirror the geometry, *'x'* or *'y'*
        :param mirror_point: Point about which to mirror the geometry *(x, y)*.
            If no point is provided, mirrors the geometry about the centroid of the shape's bounding box.
            Default = 'center'.
        :type mirror_point: Union[list[float, float], str]

        :return: Geometry object mirrored on 'axis' about 'mirror_point'
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

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

            import sectionproperties.pre.library.steel_sections as steel_sections
            import sectionproperties.pre.library.primitive_sections as primitive_sections

            geometry_1 = steel_sections.channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
            geometry_2 = primitive_sections.rectangular_section(d=20, b=133)
            compound = geometry_2.align_center(geometry_1).align_to(geometry_1, on="top") + geometry_1
            new_compound = compound.mirror_section(axis='y', mirror_point=[0,0])
        """
        geoms_acc = []
        for geom in self.geoms:
            geoms_acc.append(geom.mirror_section(axis, mirror_point))
        new_geom = CompoundGeometry(geoms_acc)
        return new_geom

    def align_center(
        self, align_to: Optional[Union[Geometry, Tuple[float, float]]] = None
    ):
        """
        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 E modulus of each
        geometry's assigned material.

        :param align_to: Another Geometry to align to, an xy coordinate, or None (default is None)
        :type align_to: Optional[Union[:class:`~sectionproperties.pre.geometry.Geometry`, Tuple[float, float]]]

        :return: Geometry object translated to new alignment
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`
        """
        EA_sum = sum(
            [
                geom.material.elastic_modulus * geom.calculate_area()
                for geom in self.geoms
            ]
        )
        cx_EA_acc = 0
        cy_EA_acc = 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(Geometry):
            align_cx, align_cy = list(align_to.geom.centroid.coords)[0]
            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):
                raise ValueError(
                    f"align_to must be either a Geometry object or an x, y coordinate, not {align_to}."
                )

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

    def split_section(
        self,
        point_i: Tuple[float, float],
        point_j: Optional[Tuple[float, float]] = None,
        vector: Union[Optional[Tuple[float, float]], np.ndarray] = None,
    ) -> Tuple[List[Geometry], List[Geometry]]:
        """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.

        :param point_i: A tuple of *(x, y)* coordinates to define a first point on the line
        :type point_i: Tuple[float, float]

        :param point_j: Optional. A tuple of *(x, y)* coordinates to define a second point on the line
        :type point_j: Tuple[float, float]

        :param vector: Optional. A tuple or numpy ndarray of *(x, y)* components to define the line direction.
        :type vector: Union[Tuple[float, float], numpy.ndarray]

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

        :rtype: Tuple[List[Geometry], List[Geometry]]

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

            import sectionproperties.pre.library.steel_sections as steel_sections
            from shapely.geometry import LineString

            geometry = steel_sections.channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
            right_geom, left_geom = geometry.split_section((0, 0), (0, 1))
        """
        top_geoms_acc = []
        bottom_geoms_acc = []
        for geom in self.geoms:
            top_geoms, bottom_geoms = geom.split_section(point_i, point_j, 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, where="exterior", resolution: float = 12
    ):
        """Dilates or erodes the perimeter of a CompoundGeometry object by a discrete amount.

        :param amount: Distance to offset the section by. A -ve value "erodes" the section. A +ve
            value "dilates" the section.
        :type amount: float
        :param 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.
            Default is "exterior".
        :type where: str
        :param resolution: Number of segments used to approximate a quarter circle around a point
        :type resolution: float
        :return: Geometry object translated to new alignment
        :rtype: :class:`~sectionproperties.pre.geometry.Geometry`

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

            import sectionproperties.pre.library.steel_sections as steel_sections
            import sectionproperties.pre.library.primitive_sections as primitive_sections

            geometry_1 = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
            geometry_2 = primitive_sections.rectangular_section(d=12, b=133)
            compound = geometry_2.align_center(geometry_1).align_to(geometry_1, on="top") + geometry_1
            new_geometry = compound.offset_perimeter(amount=-2)

        .. 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])
            offset_geom = Geometry(unionized_poly).offset_perimeter(
                amount, where, resolution
            )
            # Using the offset_geom as a "mask"
            geoms_acc = []
            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)
            new_geom = CompoundGeometry(geoms_acc)
            return new_geom

        elif amount > 0:  # Ballooning condition
            # This produces predictable results up to a point.
            # That point is when the offset is so great it exceeds the thickness
            # of the material at an interface of two materials.
            # e.g. A 50 deep plate on top of the top flange of an I Section with a flange depth of 10
            # When the offset exceeds 10 (the depth of the flange at the intersection), the meshed
            # material regions will become unpredictable.

            geoms_acc = []
            for i_idx, geom in enumerate(self.geoms):
                # Offset each geom...
                offset_geom = geom.offset_perimeter(amount, where, resolution)
                for j_idx, orig_geom in enumerate(self.geoms):
                    if i_idx != j_idx:
                        # ... then remove the parts that intersect with the other
                        # constituents of the compound geometry (because they are
                        # occupying that space already)
                        offset_geom = offset_geom - orig_geom
                if not offset_geom.geom.is_empty:
                    geoms_acc.append(offset_geom)
            new_geom = CompoundGeometry(geoms_acc)
            return new_geom
        else:
            return self

    def compile_geometry(self):
        """
        Converts the shapely.geometry.Polygon objects stored in self.geoms into lists of
        points, facets, control_points, and hole points.
        """
        self.points = []
        self.facets = []
        self.control_points = []
        self.holes = []

        # 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 list(point) not in self.points:
                    self.points.append(list(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 holes
            existing_geom_holes = []
            for hole in geom.holes:
                existing_geom_holes.append(tuple(hole))

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

        # Determine if new holes have been created or if existing
        # holes have been destroyed (or "filled in").
        resultant_holes = []
        unionized_poly = unary_union([geom.geom for geom in self.geoms])
        if isinstance(unionized_poly, MultiPolygon):
            for poly in unionized_poly.geoms:
                for interior in poly.interiors:
                    resultant_holes.append(
                        tuple(Polygon(interior).representative_point().coords[0])
                    )

        elif isinstance(unionized_poly, Polygon):
            if Geometry(unionized_poly).holes:
                resultant_holes = Geometry(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
        """
        unionized_poly = unary_union([geom.geom for geom in self.geoms])
        if isinstance(unionized_poly, MultiPolygon):
            return -1.0
        return unionized_poly.exterior.length


### Helper functions for Geometry


def load_dxf(dxf_filepath: pathlib.Path):
    """
    Import any-old-shape in dxf format for analysis.
    Code by aegis1980 and connorferster
    """
    c2s = None
    try:
        import cad_to_shapely as c2s  # type: ignore
    except ImportError as e:
        print(e)
        print("To use 'from_dxf(...)' you need to 'pip install cad_to_shapely'")
        return

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

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

    my_dxf = c2s.dxf.DxfImporter(dxf_filepath)
    my_dxf.process()
    my_dxf.cleanup()

    polygons = my_dxf.polygons
    new_polygons = c2s.utils.find_holes(polygons)
    if isinstance(new_polygons, MultiPolygon):
        return CompoundGeometry(new_polygons)
    elif isinstance(new_polygons, Polygon):
        return Geometry(new_polygons)
    else:
        print(f"No shapely.Polygon objects found in file: {dxf_filepath}")


def create_facets(
    points_list: list, connect_back: bool = False, offset: int = 0
) -> list:
    """
    Returns a list of lists of integers representing the "facets" connecting
    the list of coordinates in 'loc'. It is assumed that 'loc' coordinates are
    already in their order of connectivity.

    'loc': a list of coordinates
    'connect_back': if True, then the last facet pair will be [len(loc), offset]
    'offset': an integer representing the value that the facets should begin incrementing from.
    """
    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:
    """
    Return a list of lists representing x,y pairs of the exterior
    perimeter of `polygon`.
    """
    acc = [list(coord) for coord in shape.exterior.coords]
    return acc


def create_interior_points(lr: LinearRing) -> list:
    """
    Return a list of lists representing x,y pairs of the exterior
    perimeter of `polygon`.
    """
    acc = [list(coord) for coord in lr.coords]
    return acc


def create_points_and_facets(shape: Polygon, tol=12) -> tuple:
    """
    Return a list of lists representing x,y pairs of the exterior
    perimeter of `polygon`.
    """
    master_count = 0
    points = []
    facets = []

    # Shape perimeter
    for coords in list(
        shape.exterior.coords[:-1]
    ):  # The last point == first point (shapely)
        x, y = list(coords)
        x, y = round(x, tol), round(y, tol)
        points.append([x, y])
        master_count += 1
    facets += create_facets(points, connect_back=True)
    exterior_count = master_count

    # Holes
    for idx, hole in enumerate(shape.interiors):
        break_count = master_count
        int_points = []
        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

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

    return points, facets


def buffer_polygon(polygon: Polygon, amount: float, resolution: int):
    buffered_polygon = polygon.buffer(
        distance=amount, join_style=1, resolution=resolution
    )
    if isinstance(buffered_polygon, GeometryCollection):
        remaining_polygons = []
        for item in buffered_polygon.geoms:
            if isinstance(item, Polygon):
                remaining_polygons.append(Polygon)
        if len(remaining_polygons) == 1:
            return remaining_polygons[0]
        else:
            return MultiPolygon(remaining_polygons)
    elif isinstance(buffered_polygon, MultiPolygon) or isinstance(
        buffered_polygon, Polygon
    ):
        return buffered_polygon
    else:
        return Polygon()


def filter_non_polygons(
    input_geom: Union[GeometryCollection, LineString, Point, Polygon, MultiPolygon]
) -> Union[Polygon, MultiPolygon]:
    """
    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.
    """
    if isinstance(input_geom, (Polygon, MultiPolygon)):
        return input_geom
    elif isinstance(input_geom, GeometryCollection):
        acc = []
        for item in input_geom.geoms:
            if isinstance(item, MultiPolygon):
                acc.append(item)
            elif isinstance(item, Polygon):
                acc.append(item)
        if len(acc) == 0:
            return Polygon()
        elif len(acc) == 1:
            return acc[0]
        else:
            return MultiPolygon(acc)
    elif isinstance(input_geom, (Point, LineString)):
        return Polygon()


def round_polygon_vertices(poly: Polygon, tol: int) -> Polygon:
    """
    Returns a Polygon representing 'poly' with all of its vertices rounded by 'tol'.
    """
    rounded_exterior = np.round(poly.exterior.coords, tol)
    rounded_interiors = []
    for interior in poly.interiors:
        rounded_interiors.append(np.round(interior.coords, tol))
    # print("Exterior: ", rounded_exterior)
    # print("Interior: ", rounded_interiors)
    if not rounded_exterior.any():
        return Polygon()
    return Polygon(rounded_exterior, rounded_interiors)


def check_geometry_overlaps(lop: List[Polygon]) -> bool:
    """
    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.
    """
    union_area = unary_union(lop).area
    sum_polygons = sum([poly.area for poly in lop])
    return not math.isclose(union_area, sum_polygons)


def check_geometry_disjoint(lop: List[Polygon]) -> bool:
    """
    Returns True if any polygons in 'lop' are disjoint. Returns
    False, otherwise.
    """
    # Build polygon connectivity network
    network = {}
    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

    def walk_network(node: int, network: dict, nodes_visited: list[int]) -> list[int]:
        """
        Walks the network modifying 'nodes_visited' as it walks.
        """
        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())