Repository: robbievanleeuwen/section-properties (language: python)

Issue:
Compound from points (#113)
In response to discussions on #113. The following changes have been made:

* `CompoundGeometry.from_points()` has an updated algorithm
* `Geometry.from_points()` no longer has optional arguments (other than `holes`)
* The order of arguments in both `.from_points()` methods have been updated to the following:
  * `points`, `facets`, `control_points`, `holes = None`
  * Doc strings have been updated to reflect this
* Tests in test_sections.py have been updated to reflect the above
* An additional test was created in test_sections.py: `test_nested_compound_geometry_from_points()`
* An additional example was created for advanced_geom.rst to describe the modelling of nested geometries
* Black formatting was performed which changed the formatting in test_sections.py and test_validation.py

All tests pass
Implement calculation of cross-section mass
Summary of changes:

- Material properties now require `density` variable
- `calculate_geometric_properties()` now calculates cross-section mass
- Added `get_mass()` method to `Section` class
- `display_results()` reports mass for composite cross-sections
- Updated documentation to suit new methods
- Updated documentation where materials are mentioned to include the `density` variable
- Updated tests

See #128.
Add github actions
Begins to address #80. Not diving in too deep yet, have added the following github actions as a starter:

- Lint - lints with black and prepares a PR if there are changes that need to be made
- ~~Test installation - ensures the package can be installed on all operating systems~~
- Run tests - runs all the tests for *sectionproperties* on all operating systems
- Code coverage - generates a codecov report, encouraging further tests to be added


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/black.yml b/.github/workflows/black.yml
new file mode 100644
index 00000000..ddb30a72
--- /dev/null
+++ b/.github/workflows/black.yml
@@ -0,0 +1,13 @@
+name: Lint with Black
+
+on: [push, pull_request]
+
+jobs:
+  lint:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: psf/black@stable
+        with:
+          options: "--check --verbose" 
diff --git a/.github/workflows/code_coverage.yml b/.github/workflows/code_coverage.yml
new file mode 100644
index 00000000..161ed36b
--- /dev/null
+++ b/.github/workflows/code_coverage.yml
@@ -0,0 +1,38 @@
+name: Code Coverage
+
+on:
+  push:
+    branches: [ master, v2 ]
+  pull_request:
+    branches: [ master, v2 ]
+  workflow_dispatch:
+
+jobs:
+  tests:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: [ '3.8' ]
+
+    name: Code coverage with python ${{ matrix.python-version }}
+    steps:
+    - uses: actions/checkout@v2
+
+    - name: Setup python ${{ matrix.python-version }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: ${{ matrix.python-version }}
+
+    - name: Update pip
+      run: pip install --upgrade pip
+
+    - name: Install package
+      run: python -m pip install -e .
+
+    - name: Generate report
+      run: |
+        pip install pytest-cov
+        pytest --cov=./ --cov-report=xml
+
+    - name: Upload coverage to Codecov
+      uses: codecov/codecov-action@v2.1.0
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
new file mode 100644
index 00000000..3b67a640
--- /dev/null
+++ b/.readthedocs.yaml
@@ -0,0 +1,22 @@
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-20.04
+  tools:
+    python: "3.8"
+  apt_packages:
+    - cmake
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+   configuration: docs/source/conf.py
+
+# Default is HTML
+formats: []
+
+# Optionally declare the Python requirements required to build your docs
+python:
+   install:
+   - requirements: docs/source/requirements_docs.txt
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index e9fc34d3..00000000
--- a/.travis.yml
+++ /dev/null
@@ -1,15 +0,0 @@
-language: python
-
-before_install:
-  - pip install pybind11
-
-python:
-  - "3.8"
-
-install:
-  - pip install --upgrade pip
-  - python setup.py install
-
-script:
-  - python ./sectionproperties/tests/test_rectangle.py
-  - python ./sectionproperties/tests/test_validation.py
diff --git a/MANIFEST.in b/MANIFEST.in
deleted file mode 100644
index 18e06834..00000000
--- a/MANIFEST.in
+++ /dev/null
@@ -1,1 +0,0 @@
-include CHANGELOG.rst
diff --git a/README.md b/README.md
index 90745345..39b1a0d5 100644
--- a/README.md
+++ b/README.md
@@ -10,25 +10,6 @@ A python package for the analysis of arbitrary cross-sections using the finite e
 
 For more detailed installation instructions, refer to the [documentation](https://sectionproperties.readthedocs.io/).
 
-### UNIX (MacOS/Linux):
-
-```
-$ pip install sectionproperties
-```
-
-### Windows
-
-Install *meshpy* by downloading the appropriate [installation wheel](https://www.lfd.uci.edu/~gohlke/pythonlibs/#meshpy).
-
-Navigate to the location of the downloaded wheel and install using pip:
-
-```
-$ cd Downloads
-$ pip install MeshPy‑2018.2.1‑cp36‑cp36m‑win_amd64.whl
-```
-
-Once *meshpy* has been installed, *sectionproperties* can be installed:
-
 ```
 $ pip install sectionproperties
 ```
diff --git a/codecov.yml b/codecov.yml
new file mode 100644
index 00000000..bfa8c68c
--- /dev/null
+++ b/codecov.yml
@@ -0,0 +1,17 @@
+codecov:
+  require_ci_to_pass: yes
+
+coverage:
+  precision: 2
+  round: down
+  range: "70...100"
+
+parsers:
+  gcov:
+    branch_detection:
+      conditional: yes
+      loop: yes
+      method: no
+      macro: no
+
+comment: false
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 0464e4cc..01b2ea3d 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -7,6 +7,10 @@
 # http://www.sphinx-doc.org/en/master/config
 
 from sphinx_gallery.sorting import FileNameSortKey
+import locale
+
+
+locale.setlocale(locale.LC_ALL, '')
 
 # -- Path setup --------------------------------------------------------------
 
@@ -16,20 +20,22 @@
 
 import os
 import sys
-sys.path.insert(0, os.path.abspath('../../'))
+
+sys.path.insert(0, os.path.abspath("../../"))
 from sectionproperties import __version__ as ver
-autodoc_mock_imports = ['rhino_shapely_interop', 'mpl_toolkits']
+
+autodoc_mock_imports = ["rhino_shapely_interop", "cad_to_shapely"]
 
 # -- Project information -----------------------------------------------------
 
-project = 'sectionproperties'
-copyright = '2021, Robbie van Leeuwen'
-author = 'Robbie van Leeuwen'
+project = "sectionproperties"
+copyright = "2021, Robbie van Leeuwen"
+author = "Robbie van Leeuwen"
 
 # The short Major.Minor.Build version
-_v = ver.split('.')
-_build = ''.join([c for c in _v[2] if c.isdigit()])
-version = _v[0] + '.' + _v[1] + '.' + _build
+_v = ver.split(".")
+_build = "".join([c for c in _v[2] if c.isdigit()])
+version = _v[0] + "." + _v[1] + "." + _build
 # The full version, including alpha/beta/rc tags
 release = ver
 
@@ -51,23 +57,23 @@
 ]
 
 # Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
+templates_path = ["_templates"]
 
 # The suffix(es) of source filenames.
 # You can specify multiple suffix as a list of string:
 #
 # source_suffix = ['.rst', '.md']
-source_suffix = '.rst'
+source_suffix = ".rst"
 
 # The master toctree document.
-master_doc = 'index'
+master_doc = "index"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
 #
 # This is also used if you do content translation via gettext catalogs.
 # Usually you set "language" from the command line for these cases.
-language = 'en'
+language = "en"
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
@@ -75,7 +81,7 @@
 exclude_patterns = []
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = 'sphinx'
+pygments_style = "sphinx"
 
 
 # -- Options for HTML output -------------------------------------------------
@@ -83,7 +89,7 @@
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = 'sphinx_rtd_theme'
+html_theme = "sphinx_rtd_theme"
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
@@ -110,7 +116,7 @@
 # -- Options for HTMLHelp output ---------------------------------------------
 
 # Output file base name for HTML help builder.
-htmlhelp_basename = 'sectionpropertiesdoc'
+htmlhelp_basename = "sectionpropertiesdoc"
 
 
 # -- Options for LaTeX output ------------------------------------------------
@@ -119,15 +125,12 @@
     # The paper size ('letterpaper' or 'a4paper').
     #
     # 'papersize': 'letterpaper',
-
     # The font size ('10pt', '11pt' or '12pt').
     #
     # 'pointsize': '10pt',
-
     # Additional stuff for the LaTeX preamble.
     #
     # 'preamble': '',
-
     # Latex figure (float) alignment
     #
     # 'figure_align': 'htbp',
@@ -137,8 +140,13 @@
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
-    (master_doc, 'sectionproperties.tex', 'sectionproperties Documentation',
-     'Robbie van Leeuwen', 'manual'),
+    (
+        master_doc,
+        "sectionproperties.tex",
+        "sectionproperties Documentation",
+        "Robbie van Leeuwen",
+        "manual",
+    ),
 ]
 
 
@@ -147,8 +155,7 @@
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 man_pages = [
-    (master_doc, 'sectionproperties', 'sectionproperties Documentation',
-     [author], 1)
+    (master_doc, "sectionproperties", "sectionproperties Documentation", [author], 1)
 ]
 
 
@@ -158,9 +165,15 @@
 # (source start file, target name, title, author,
 #  dir menu entry, description, category)
 texinfo_documents = [
-    (master_doc, 'sectionproperties', 'sectionproperties Documentation',
-     author, 'sectionproperties', 'One line description of project.',
-     'Miscellaneous'),
+    (
+        master_doc,
+        "sectionproperties",
+        "sectionproperties Documentation",
+        author,
+        "sectionproperties",
+        "One line description of project.",
+        "Miscellaneous",
+    ),
 ]
 
 
diff --git a/docs/source/images/examples/shapely_z_ordering.png b/docs/source/images/examples/shapely_z_order.png
similarity index 100%
rename from docs/source/images/examples/shapely_z_ordering.png
rename to docs/source/images/examples/shapely_z_order.png
diff --git a/docs/source/requirements_docs.txt b/docs/source/requirements_docs.txt
index 102b8e20..454455ae 100644
--- a/docs/source/requirements_docs.txt
+++ b/docs/source/requirements_docs.txt
@@ -3,8 +3,8 @@ numpy
 shapely
 matplotlib
 more_itertools
-pybind11
-meshpy
+triangle
 Sphinx
 sphinx-rtd-theme
 sphinx-gallery
+rhino-shapley-interop
diff --git a/docs/source/rst/advanced_geom.rst b/docs/source/rst/advanced_geom.rst
index 25d6d56c..8f017f45 100644
--- a/docs/source/rst/advanced_geom.rst
+++ b/docs/source/rst/advanced_geom.rst
@@ -232,8 +232,8 @@ Another example (but with nested geometries)
 --------------------------------------------
 
 This example demonstrates creating nested geometries using two different approaches.
-These approaches reflect the differences between how shapely (geometry pre-processor) "perceives" geometry 
-and how meshpy/Triangle (meshing tool) "perceives" geometry and how the modeller might adapt their
+These approaches reflect the differences between how shapely (geometry pre-processor) "perceives" geometry
+and how Triangle (meshing tool) "perceives" geometry and how the modeller might adapt their
 input style depending on the situation.
 
 The nested geometry we are trying to create looks as follows:
@@ -252,7 +252,7 @@ In creating this geometry consider the following:
     :align: center
     :scale: 100 %
 
-- meshpy/Triangle does not have a concept of "z-ordering" so there is only a single plane which
+- Triangle does not have a concept of "z-ordering" so there is only a single plane which
   may have regions of different materials (specified with control points). When a hole is created
   in the plane, it "punches" through "all" polygons in the plane.
 
@@ -260,7 +260,7 @@ In creating this geometry consider the following:
     :align: center
     :scale: 100 %
 
-To create the nested geometry using shapely, the code would be as follows:
+To create the nested geometry using shapely, the code would be as follows::
 
     mat1 = Material(name="Material 1", elastic_modulus=100, poissons_ratio=0.3, yield_strength=10, color="yellow")
     mat2 = Material(name="Material 2", elastic_modulus=100, poissons_ratio=0.3, yield_strength=10, color="orange")
@@ -282,7 +282,7 @@ To create the nested geometry using shapely, the code would be as follows:
     :align: center
     :scale: 100 %
 
-To create the nested geometry using the Triangle interface, the code would be as follows:
+To create the nested geometry using the Triangle interface, the code would be as follows::
 
     points = [ # Points for four squares are created
         [-50.0, 50.0], # Square 1
diff --git a/docs/source/rst/geom_mesh.rst b/docs/source/rst/geom_mesh.rst
index 7099496a..3e422b33 100644
--- a/docs/source/rst/geom_mesh.rst
+++ b/docs/source/rst/geom_mesh.rst
@@ -142,7 +142,7 @@ file (i.e. the number of contiguous regions).
 Geometry from Rhino
 -------------------
 
-Geometries can now be created from .3dm files and BREP encodings. 
+Geometries can now be created from .3dm files and BREP encodings.
 Various limitations and assumptions need to be acknowledged:
 
 * sectional analysis is based in 2d and Rhino is a 3d environment.
@@ -152,7 +152,7 @@ Various limitations and assumptions need to be acknowledged:
 
 See the keyword arguments below that are used to search and simplify the Rhino geometry.
 
-Rhino files are read via the class methods :attr:`sectionproperties.pre.sections.Geometry.from_3dm()` and 
+Rhino files are read via the class methods :attr:`sectionproperties.pre.sections.Geometry.from_3dm()` and
 :attr:`sectionproperties.pre.sections.CompoundGeometry.from_3dm()`.
 Each class method returns the respective objects.
 
@@ -168,7 +168,7 @@ Geometry objects can also be created from encodings of Rhino BREP.
    :noindex:
 
 More advanced filtering can be achieved by working with the Shapely geometries directly.
-These can be accessed by :attr:`sectionproperties.pre.rhino.load_3dm()` and 
+These can be accessed by :attr:`sectionproperties.pre.rhino.load_3dm()` and
 :attr:`sectionproperties.pre.rhino.load_rhino_brep_encoding()`.
 
 ..  autofunction:: sectionproperties.pre.rhino.load_3dm
@@ -349,6 +349,9 @@ Offsetting the Perimeter
   ..  autofunction:: sectionproperties.pre.sections.Geometry.offset_perimeter
       :noindex:
 
+  ..  autofunction:: sectionproperties.pre.sections.CompoundGeometry.offset_perimeter
+      :noindex:
+
 
 Visualising the Geometry
 ========================
diff --git a/docs/source/rst/installation.rst b/docs/source/rst/installation.rst
index db6497e9..beb07ed5 100644
--- a/docs/source/rst/installation.rst
+++ b/docs/source/rst/installation.rst
@@ -2,48 +2,17 @@ Installation
 ============
 
 These instructions will get you a copy of *sectionproperties* up and running on
-your local machine. You will need a working copy of python 3.8 on your machine.
+your local machine. You will need a working copy of python on your machine.
 
 Installing *sectionproperties*
 ------------------------------
 
-*sectionproperties* uses *meshpy* to efficiently generate a conforming triangular
-mesh in order to perform a finite element analysis of the structural cross-section.
-The installation procedure for *meshpy* depends on your local machine.
+*sectionproperties* uses `shapely <https://github.com/shapely/shapely>`_ to prepare the
+cross-section geometry and `triangle <https://github.com/drufat/triangle>`_ to efficiently
+generate a conforming triangular mesh in order to perform a finite element analysis of the
+structural cross-section.
 
-UNIX (MacOS/Linux)
-^^^^^^^^^^^^^^^^^^
-
-*sectionproperties* and all of its dependencies can be installed through the
-python package index::
-
-  $ pip install sectionproperties
-
-If you have any issues installing *meshpy*, refer to the installation instructions
-on the `MeshPy GitHub page <https://github.com/inducer/meshpy>`_ or the
-`MeshPy documentation <https://documen.tician.de/meshpy/installation.html>`_.
-
-If you have any issues installing *shapely*, refer to the installation instructions
-on the `Shapely GitHub page <https://github.com/Toblerity/Shapely>`_ or the
-`Shapely documentation <https://shapely.readthedocs.io/en/stable/manual.html>`_.
-
-Windows
-^^^^^^^
-
-If you use conda, you can install *meshpy* and *shapely* directly::
-
-  $ conda install -c conda-forge meshpy shapely
-
-Alternatively, you can install both *meshpy* and shapely by downloading the appropriate
-`installation wheels <https://www.lfd.uci.edu/~gohlke/pythonlibs>`_.
-
-Navigate to the location of the downloaded wheel and install using pip::
-
-  $ cd Downloads
-  $ pip install <wheel file names go here>
-
-Once *meshpy* and *shapely* have been installed, the rest of the *sectionproperties* package can
-be installed using the python package index::
+*sectionproperties* and all of its dependencies can be installed through the python package index::
 
   $ pip install sectionproperties
 
diff --git a/docs/source/rst/theory.rst b/docs/source/rst/theory.rst
index 542f7744..1d102b05 100644
--- a/docs/source/rst/theory.rst
+++ b/docs/source/rst/theory.rst
@@ -33,9 +33,9 @@ Mesh Generation
 ---------------
 
 The cross-section is meshed using quadratic superparametric triangular elements (Tri6) using the
-`meshpy <https://github.com/inducer/meshpy>`_ library for Python. Superparametric quadratic
+`triangle <https://github.com/drufat/triangle>`_ library for Python. Superparametric quadratic
 elements are defined as having straight edges and mid-nodes located at the mid-point between
-adjacent corner nodes. *meshpy* implements
+adjacent corner nodes. *triangle* implements
 `Triangle <https://www.cs.cmu.edu/~quake/triangle.html>`_, which is a two dimensional quality mesh
 generator and Delaunay triangulator written by Jonathan Shewchuk in C.
 
diff --git a/examples/00-basic/03_custom.py b/examples/00-basic/03_custom.py
index 12a1ec77..10c3780e 100644
--- a/examples/00-basic/03_custom.py
+++ b/examples/00-basic/03_custom.py
@@ -44,7 +44,7 @@
 # %%
 # Because we have two separate geometry regions (as indicated by our control_points)
 # we create a CompoundGeometry from points
-geometry = sections.CompoundGeometry.from_points(points, facets, holes, control_points)
+geometry = sections.CompoundGeometry.from_points(points, facets, control_points, holes)
 
 # %%
 # Create the mesh and section. For the mesh, use a smaller refinement for the angle region.
diff --git a/examples/00-basic/04_merged.py b/examples/00-basic/04_merged.py
index 05aa74fa..73b724fe 100644
--- a/examples/00-basic/04_merged.py
+++ b/examples/00-basic/04_merged.py
@@ -32,11 +32,11 @@
 rhs = sections.rectangular_hollow_section(d=150, b=100, t=6, r_out=15, n_r=8)
 
 # %%
-# Create a triangular section from points
-# only the points are needed because they are in sequential order
-# and represent only a single, contiguous region
+# Create a triangular section from points, facets, and control points
 points = [[0, 0], [50, 0], [25, 50]]
-triangle = sections.Geometry.from_points(points)
+facets = [[0, 1], [1, 2], [2, 0]]
+control_points = [[25,25]]
+triangle = sections.Geometry.from_points(points, facets, control_points)
 triangle = triangle.align_center(rhs).align_to(rhs, on="top")
 
 # %%
diff --git a/examples/00-basic/07_composite.py b/examples/00-basic/07_composite.py
index f3ac2661..6cb8cb11 100644
--- a/examples/00-basic/07_composite.py
+++ b/examples/00-basic/07_composite.py
@@ -26,9 +26,9 @@
 # %%
 # Create material properties
 steel = Material(name='Steel', elastic_modulus=200e3, poissons_ratio=0.3,
-                 yield_strength=500, color='grey')
+                 yield_strength=500, density=8.05e-6 ,color='grey')
 timber = Material(name='Timber', elastic_modulus=8e3, poissons_ratio=0.35,
-                  yield_strength=20, color='burlywood')
+                  yield_strength=20, density=0.78e-6, color='burlywood')
 
 # %%
 # Create 310UB40.4
diff --git a/sectionproperties/analysis/cross_section.py b/sectionproperties/analysis/cross_section.py
index 669d2aea..13736df0 100644
--- a/sectionproperties/analysis/cross_section.py
+++ b/sectionproperties/analysis/cross_section.py
@@ -50,8 +50,8 @@ class Section:
     :cvar int num_nodes: Number of nodes in the finite element mesh
     :cvar geometry: Cross-section geometry object used to generate the mesh
     :vartype geometry: :class:`~sectionproperties.pre.sections.Geometry`
-    :cvar mesh: Mesh object returned by meshpy
-    :vartype mesh: :class:`meshpy.triangle.MeshInfo`
+    :cvar mesh: Mesh dict returned by triangle
+    :vartype mesh: dict(mesh)
     :cvar mesh_nodes: Array of node coordinates from the mesh
     :vartype mesh_nodes: :class:`numpy.ndarray`
     :cvar mesh_elements: Array of connectivities from the mesh
@@ -95,9 +95,9 @@ def init():
                 self.materials = [self.geometry.material]
 
             # extract mesh data
-            nodes = np.array(mesh.points, dtype=np.dtype(float))
-            elements = np.array(mesh.elements, dtype=np.dtype(int))
-            attributes = np.array(mesh.element_attributes, dtype=np.dtype(int))
+            nodes = np.array(mesh["vertices"], dtype=np.dtype(float))
+            elements = np.array(mesh["triangles"], dtype=np.dtype(int))
+            attributes = np.array(mesh["triangle_attributes"].T[0], dtype=np.dtype(int))
 
             # swap mid-node order to retain node ordering consistency
             elements[:, [3, 4, 5]] = elements[:, [5, 3, 4]]
@@ -235,7 +235,17 @@ def calculate_geom():
 
             # calculate global geometric properties
             for el in self.elements:
-                (area, qx, qy, ixx_g, iyy_g, ixy_g, e, g, rho) = el.geometric_properties()
+                (
+                    area,
+                    qx,
+                    qy,
+                    ixx_g,
+                    iyy_g,
+                    ixy_g,
+                    e,
+                    g,
+                    rho,
+                ) = el.geometric_properties()
                 self.section_props.area += area
                 self.section_props.mass += area * rho
                 self.section_props.ea += area * e
@@ -855,12 +865,16 @@ def calculate_stress(self, N=0, Vx=0, Vy=0, Mxx=0, Myy=0, M11=0, M22=0, Mzz=0):
         """
 
         # check that a geometric and warping analysis has been performed
-        if None in [
-            self.section_props.area,
-            self.section_props.ixx_c,
-            self.section_props.cx,
-            self.section_props.j,
-        ] and self.section_props.omega is None:
+        if (
+            None
+            in [
+                self.section_props.area,
+                self.section_props.ixx_c,
+                self.section_props.cx,
+                self.section_props.j,
+            ]
+            and self.section_props.omega is None
+        ):
             err = "Perform a geometric and warping analysis before carrying out a stress analysis."
             raise RuntimeError(err)
 
@@ -1056,7 +1070,16 @@ def assemble_torsion(self, lg=True):
 
         return (csc_matrix(k), csc_matrix(k_lg), f_torsion)
 
-    def plot_mesh(self, ax=None, pause=True, alpha=0.5, materials=True, mask=None, size=500, dpi=96):
+    def plot_mesh(
+        self,
+        ax=None,
+        pause=True,
+        alpha=0.5,
+        materials=True,
+        mask=None,
+        size=500,
+        dpi=96,
+    ):
         """Plots the finite element mesh. If no axes object is supplied a new figure and axis is
         created.
 
@@ -1142,9 +1165,7 @@ def plot_mesh(self, ax=None, pause=True, alpha=0.5, materials=True, mask=None, s
                 if i == 0 or material not in self.materials[0:i]:
                     # add the material colour and name to the legend list
                     patch = mpatches.Patch(color=material.color, label=material.name)
-                    legend_labels.append(
-                        patch
-                    )
+                    legend_labels.append(patch)
 
             cmap = ListedColormap(color_array)  # custom colormap
 
@@ -1918,8 +1939,8 @@ class PlasticSection:
     :vartype materials: list[:class:`~sectionproperties.pre.pre.Material`]
     :cvar bool debug: If set to True, the geometry is plotted each time a new mesh is generated by
         the plastic centroid algorithm.
-    :cvar mesh: Mesh object returned by meshpy
-    :vartype mesh: :class:`meshpy.triangle.MeshInfo`
+    :cvar mesh: Mesh dict returned by triangle
+    :vartype mesh: dict(mesh)
     :cvar mesh_nodes: Array of node coordinates from the mesh
     :vartype mesh_nodes: :class:`numpy.ndarray`
     :cvar mesh_elements: Array of connectivities from the mesh
@@ -2325,7 +2346,9 @@ def __init__(self, cross_section):
         # can be saved to a new material group
         self.material_groups = copy.deepcopy(cross_section.material_groups)
 
-    def plot_stress_contour(self, sigs, title, pause, cmap, normalize=True, size=500, dpi=96):
+    def plot_stress_contour(
+        self, sigs, title, pause, cmap, normalize=True, size=500, dpi=96
+    ):
         """Plots filled stress contours over the finite element mesh.
 
         :param sigs: List of nodal stress values for each material
@@ -2400,7 +2423,9 @@ def plot_stress_contour(self, sigs, title, pause, cmap, normalize=True, size=500
 
         return (fig, ax)
 
-    def plot_stress_vector(self, sigxs, sigys, title, pause, cmap, normalize=False, size=500, dpi=96):
+    def plot_stress_vector(
+        self, sigxs, sigys, title, pause, cmap, normalize=False, size=500, dpi=96
+    ):
         """Plots stress vectors over the finite element mesh.
 
         :param sigxs: List of x-components of the nodal stress values for each material
@@ -2438,7 +2463,6 @@ def plot_stress_vector(self, sigxs, sigys, title, pause, cmap, normalize=False,
         if normalize:
             norm = CenteredNorm()
 
-
         # plot the vectors
         for (i, sigx) in enumerate(sigxs):
             sigy = sigys[i]
@@ -2634,7 +2658,9 @@ def get_stress(self):
 
         return stress
 
-    def plot_stress_n_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_n_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the normal stress :math:`\sigma_{zz,N}` resulting from the
         axial load :math:`N`.
 
@@ -2678,9 +2704,13 @@ def plot_stress_n_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz_n)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_mxx_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_mxx_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the normal stress :math:`\sigma_{zz,Mxx}` resulting from the
         bending moment :math:`M_{xx}`.
 
@@ -2724,9 +2754,13 @@ def plot_stress_mxx_zz(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz_mxx)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_myy_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_myy_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the normal stress :math:`\sigma_{zz,Myy}` resulting from the
         bending moment :math:`M_{yy}`.
 
@@ -2770,9 +2804,13 @@ def plot_stress_myy_zz(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz_myy)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_m11_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_m11_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the normal stress :math:`\sigma_{zz,M11}` resulting from the
         bending moment :math:`M_{11}`.
 
@@ -2816,9 +2854,13 @@ def plot_stress_m11_zz(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz_m11)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_m22_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_m22_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the normal stress :math:`\sigma_{zz,M22}` resulting from the
         bending moment :math:`M_{22}`.
 
@@ -2862,9 +2904,13 @@ def plot_stress_m22_zz(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz_m22)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_m_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_m_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the normal stress :math:`\sigma_{zz,\Sigma M}` resulting from
         all bending moments :math:`M_{xx} + M_{yy} + M_{11} + M_{22}`.
 
@@ -2909,9 +2955,13 @@ def plot_stress_m_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz_m)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_mzz_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_mzz_zx(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx,Mzz}`
         resulting from the torsion moment :math:`M_{zz}`.
 
@@ -2955,9 +3005,13 @@ def plot_stress_mzz_zx(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zx_mzz)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_mzz_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_mzz_zy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy,Mzz}`
         resulting from the torsion moment :math:`M_{zz}`.
 
@@ -3001,9 +3055,13 @@ def plot_stress_mzz_zy(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zy_mzz)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_mzz_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_mzz_zxy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy,Mzz}` resulting
         from the torsion moment :math:`M_{zz}`.
 
@@ -3047,9 +3105,13 @@ def plot_stress_mzz_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zxy_mzz)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_vector_mzz_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96):
+    def plot_vector_mzz_zxy(
+        self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96
+    ):
         """Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy,Mzz}` resulting
         from the torsion moment :math:`M_{zz}`.
 
@@ -3095,9 +3157,13 @@ def plot_vector_mzz_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=5
             sigxs.append(group.stress_result.sig_zx_mzz)
             sigys.append(group.stress_result.sig_zy_mzz)
 
-        return self.plot_stress_vector(sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_vector(
+            sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vx_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vx_zx(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx,Vx}`
         resulting from the shear force :math:`V_{x}`.
 
@@ -3141,9 +3207,13 @@ def plot_stress_vx_zx(self, pause=True, cmap="coolwarm", normalize=True, size=50
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zx_vx)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vx_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vx_zy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy,Vx}`
         resulting from the shear force :math:`V_{x}`.
 
@@ -3187,9 +3257,13 @@ def plot_stress_vx_zy(self, pause=True, cmap="coolwarm", normalize=True, size=50
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zy_vx)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vx_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vx_zxy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy,Vx}` resulting
         from the shear force :math:`V_{x}`.
 
@@ -3233,9 +3307,13 @@ def plot_stress_vx_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zxy_vx)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_vector_vx_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96):
+    def plot_vector_vx_zxy(
+        self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96
+    ):
         """Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy,Vx}` resulting
         from the shear force :math:`V_{x}`.
 
@@ -3281,9 +3359,13 @@ def plot_vector_vx_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=50
             sigxs.append(group.stress_result.sig_zx_vx)
             sigys.append(group.stress_result.sig_zy_vx)
 
-        return self.plot_stress_vector(sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_vector(
+            sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vy_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vy_zx(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx,Vy}`
         resulting from the shear force :math:`V_{y}`.
 
@@ -3327,9 +3409,13 @@ def plot_stress_vy_zx(self, pause=True, cmap="coolwarm", normalize=True, size=50
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zx_vy)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vy_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vy_zy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy,Vy}`
         resulting from the shear force :math:`V_{y}`.
 
@@ -3373,9 +3459,13 @@ def plot_stress_vy_zy(self, pause=True, cmap="coolwarm", normalize=True, size=50
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zy_vy)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vy_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vy_zxy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy,Vy}` resulting
         from the shear force :math:`V_{y}`.
 
@@ -3419,9 +3509,13 @@ def plot_stress_vy_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=5
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zxy_vy)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_vector_vy_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96):
+    def plot_vector_vy_zxy(
+        self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96
+    ):
         """Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy,Vy}` resulting
         from the shear force :math:`V_{y}`.
 
@@ -3467,9 +3561,13 @@ def plot_vector_vy_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=50
             sigxs.append(group.stress_result.sig_zx_vy)
             sigys.append(group.stress_result.sig_zy_vy)
 
-        return self.plot_stress_vector(sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_vector(
+            sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_v_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_v_zx(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *x*-component of the shear stress
         :math:`\sigma_{zx,\Sigma V}` resulting from the sum of the applied shear forces
         :math:`V_{x} + V_{y}`.
@@ -3515,9 +3613,13 @@ def plot_stress_v_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zx_v)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_v_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_v_zy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *y*-component of the shear stress
         :math:`\sigma_{zy,\Sigma V}` resulting from the sum of the applied shear forces
         :math:`V_{x} + V_{y}`.
@@ -3563,9 +3665,13 @@ def plot_stress_v_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zy_v)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_v_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_v_zxy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the resultant shear stress
         :math:`\sigma_{zxy,\Sigma V}` resulting from the sum of the applied shear forces
         :math:`V_{x} + V_{y}`.
@@ -3611,9 +3717,13 @@ def plot_stress_v_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=50
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zxy_v)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_vector_v_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96):
+    def plot_vector_v_zxy(
+        self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96
+    ):
         """Produces a vector plot of the resultant shear stress
         :math:`\sigma_{zxy,\Sigma V}` resulting from the sum of the  applied shear forces
         :math:`V_{x} + V_{y}`.
@@ -3661,9 +3771,13 @@ def plot_vector_v_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500
             sigxs.append(group.stress_result.sig_zx_v)
             sigys.append(group.stress_result.sig_zy_v)
 
-        return self.plot_stress_vector(sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_vector(
+            sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_zz(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the combined normal stress :math:`\sigma_{zz}` resulting from
         all actions.
 
@@ -3708,9 +3822,13 @@ def plot_stress_zz(self, pause=True, cmap="coolwarm", normalize=True, size=500,
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zz)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_zx(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx}`
         resulting from all actions.
 
@@ -3755,9 +3873,13 @@ def plot_stress_zx(self, pause=True, cmap="coolwarm", normalize=True, size=500,
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zx)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_zy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy}`
         resulting from all actions.
 
@@ -3802,9 +3924,13 @@ def plot_stress_zy(self, pause=True, cmap="coolwarm", normalize=True, size=500,
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zy)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_zxy(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy}` resulting
         from all actions.
 
@@ -3849,9 +3975,13 @@ def plot_stress_zxy(self, pause=True, cmap="coolwarm", normalize=True, size=500,
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_zxy)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_vector_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96):
+    def plot_vector_zxy(
+        self, pause=True, cmap="YlOrBr", normalize=False, size=500, dpi=96
+    ):
         """Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy}` resulting
         from all actions.
 
@@ -3898,9 +4028,13 @@ def plot_vector_zxy(self, pause=True, cmap="YlOrBr", normalize=False, size=500,
             sigxs.append(group.stress_result.sig_zx)
             sigys.append(group.stress_result.sig_zy)
 
-        return self.plot_stress_vector(sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_vector(
+            sigxs, sigys, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_1(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_1(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the major principal stress :math:`\sigma_{1}` resulting from all
         actions.
 
@@ -3944,15 +4078,19 @@ def plot_stress_1(self, pause=True, cmap="coolwarm", normalize=True, size=500, d
             Contour plot of the major principal stress.
         """
 
-        title = 'Stress Contour Plot - $\sigma_{1}$'
+        title = "Stress Contour Plot - $\sigma_{1}$"
         sigs = []
 
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_1)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_3(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_3(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the Minor principal stress :math:`\sigma_{3}` resulting from all
         actions.
 
@@ -3996,15 +4134,19 @@ def plot_stress_3(self, pause=True, cmap="coolwarm", normalize=True, size=500, d
             Contour plot of the minor principal stress.
         """
 
-        title = 'Stress Contour Plot - $\sigma_{3}$'
+        title = "Stress Contour Plot - $\sigma_{3}$"
         sigs = []
 
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_3)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
-    def plot_stress_vm(self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96):
+    def plot_stress_vm(
+        self, pause=True, cmap="coolwarm", normalize=True, size=500, dpi=96
+    ):
         """Produces a contour plot of the von Mises stress :math:`\sigma_{vM}` resulting from all
         actions.
 
@@ -4059,8 +4201,9 @@ def plot_stress_vm(self, pause=True, cmap="coolwarm", normalize=True, size=500,
         for group in self.material_groups:
             sigs.append(group.stress_result.sig_vm)
 
-        return self.plot_stress_contour(sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi)
-
+        return self.plot_stress_contour(
+            sigs, title, pause, cmap, normalize=normalize, size=size, dpi=dpi
+        )
 
     def plot_mohrs_circles(self, x, y, pause=True, cmap="coolwarm", size=500, dpi=96):
         """Plots Mohr's Circles of the 3D stress state at position x,y
@@ -4110,7 +4253,7 @@ def plot_mohrs_circles(self, x, y, pause=True, cmap="coolwarm", size=500, dpi=96
 
         """
 
-        pt = (x,y)
+        pt = (x, y)
         nodes = self.cross_section.mesh_nodes
         ele = self.cross_section.mesh_elements
         triang = tri.Triangulation(nodes[:, 0], nodes[:, 1], ele[:, 0:3])
@@ -4148,50 +4291,76 @@ def plot_mohrs_circles(self, x, y, pause=True, cmap="coolwarm", size=500, dpi=96
         sigma_xx = 0
         sigma_yy = 0
         tau_xy = 0
-        sigma = np.array([[sigma_xx, tau_xy, tau_xz],
-                          [tau_xy, sigma_yy, tau_yz],
-                          [tau_xz, tau_yz, sigma_zz]])
+        sigma = np.array(
+            [
+                [sigma_xx, tau_xy, tau_xz],
+                [tau_xy, sigma_yy, tau_yz],
+                [tau_xz, tau_yz, sigma_zz],
+            ]
+        )
 
         # Solve for the principal stresses using the general approach
         s, n = np.linalg.eig(sigma)
         sigma_3, sigma_2, sigma_1 = np.sort(s)
         s_l = s.tolist()
-        sigma_1i, sigma_2i, sigma_3i = s_l.index(sigma_1), s_l.index(sigma_2), s_l.index(sigma_3)
-        n1, n2, n3 = n[:,sigma_1i], n[:,sigma_2i], n[:,sigma_3i]
+        sigma_1i, sigma_2i, sigma_3i = (
+            s_l.index(sigma_1),
+            s_l.index(sigma_2),
+            s_l.index(sigma_3),
+        )
+        n1, n2, n3 = n[:, sigma_1i], n[:, sigma_2i], n[:, sigma_3i]
 
         # The tractions on each plane in cartesian coords wrt principal axes
         n_inv = np.linalg.inv(n)
         tractions = []
         for col in range(3):
-            ss = n_inv[:,col].T @ np.diag(s) @ n_inv[:,col]
-            ts = np.sqrt(np.linalg.norm(np.diag(s) @ n_inv[:,col])**2 - ss**2)
-            tractions.append((ss,ts))
+            ss = n_inv[:, col].T @ np.diag(s) @ n_inv[:, col]
+            ts = np.sqrt(np.linalg.norm(np.diag(s) @ n_inv[:, col]) ** 2 - ss ** 2)
+            tractions.append((ss, ts))
 
-        def plot_circle(ax,c,R,col,label=None,fill=None):
+        def plot_circle(ax, c, R, col, label=None, fill=None):
             circ = plt.Circle(c, R, fill=fill, ec=col, label=label)
             ax.add_patch(circ)
             ax.set_aspect(1)
             ax.autoscale_view()
 
         fig, ax = plt.subplots()
-        plot_circle(ax, (0.5*(sigma_2+sigma_3),0), 0.5*(sigma_2-sigma_3),'r', r"C1: ($\sigma_2$, $\sigma_3$)")
-        plot_circle(ax, (0.5*(sigma_1+sigma_3),0), 0.5*(sigma_1-sigma_3),'b', r"C2: ($\sigma_1$, $\sigma_3$)")
-        plot_circle(ax, (0.5*(sigma_1+sigma_2),0), 0.5*(sigma_1-sigma_2),'k', r"C3: ($\sigma_1$, $\sigma_2$)")
-        for i,plane,col in zip(range(3),['X','Y','Z'],['r','b','k']):
-            ax.plot(*tractions[i],f'{col}.', label=rf"{plane}-face")
+        plot_circle(
+            ax,
+            (0.5 * (sigma_2 + sigma_3), 0),
+            0.5 * (sigma_2 - sigma_3),
+            "r",
+            r"C1: ($\sigma_2$, $\sigma_3$)",
+        )
+        plot_circle(
+            ax,
+            (0.5 * (sigma_1 + sigma_3), 0),
+            0.5 * (sigma_1 - sigma_3),
+            "b",
+            r"C2: ($\sigma_1$, $\sigma_3$)",
+        )
+        plot_circle(
+            ax,
+            (0.5 * (sigma_1 + sigma_2), 0),
+            0.5 * (sigma_1 - sigma_2),
+            "k",
+            r"C3: ($\sigma_1$, $\sigma_2$)",
+        )
+        for i, plane, col in zip(range(3), ["X", "Y", "Z"], ["r", "b", "k"]):
+            ax.plot(*tractions[i], f"{col}.", label=rf"{plane}-face")
 
         ax.set_axisbelow(True)
-        ax.grid(which='both')
-        ax.legend(loc='upper left', bbox_to_anchor=(1.0, 1.0))
+        ax.grid(which="both")
+        ax.legend(loc="upper left", bbox_to_anchor=(1.0, 1.0))
 
-        ax.spines['top'].set_visible(False)
-        ax.spines['right'].set_visible(False)
+        ax.spines["top"].set_visible(False)
+        ax.spines["right"].set_visible(False)
 
-        ax.set_ylabel(r'Shear stress $\tau$ (MPa)')
-        ax.set_xlabel(r'Direct stress $\sigma$ (MPa)')
+        ax.set_ylabel(r"Shear stress $\tau$ (MPa)")
+        ax.set_xlabel(r"Direct stress $\sigma$ (MPa)")
 
-        ax.xaxis.set_tick_params(bottom=True, top=False, direction='inout')
-        ax.yaxis.set_tick_params(left=True, right=False, direction='inout')
+        ax.xaxis.set_tick_params(bottom=True, top=False, direction="inout")
+        ax.yaxis.set_tick_params(left=True, right=False, direction="inout")
 
         #
         # The following is just to get the labels positioned outside the axes
@@ -4200,14 +4369,20 @@ def plot_circle(ax,c,R,col,label=None,fill=None):
         x_lbl_pos = ax.xaxis.label.get_position()
         y_lbl_pos = ax.yaxis.label.get_position()
         # Make spines pass through zero of the other axis
-        ax.spines['bottom'].set_position('zero')
-        ax.spines['left'].set_position('zero')
+        ax.spines["bottom"].set_position("zero")
+        ax.spines["left"].set_position("zero")
         # Now set the coords
         ax.xaxis.set_label_coords(*x_lbl_pos)
         ax.yaxis.set_label_coords(*y_lbl_pos)
 
         fig.subplots_adjust(right=0.8)
-        post.finish_plot(ax, pause, f"Mohr's Circles for 3D Stress State at {(*pt,)}", size=size, dpi=dpi)
+        post.finish_plot(
+            ax,
+            pause,
+            f"Mohr's Circles for 3D Stress State at {(*pt,)}",
+            size=size,
+            dpi=dpi,
+        )
         return (fig, ax)
 
 
@@ -4375,8 +4550,12 @@ def calculate_combined_stresses(self):
         self.sig_zx = self.sig_zx_mzz + self.sig_zx_v
         self.sig_zy = self.sig_zy_mzz + self.sig_zy_v
         self.sig_zxy = (self.sig_zx ** 2 + self.sig_zy ** 2) ** 0.5
-        self.sig_1 = self.sig_zz / 2 + np.sqrt( (self.sig_zz / 2) ** 2 + self.sig_zxy ** 2)
-        self.sig_3 = self.sig_zz / 2 - np.sqrt( (self.sig_zz / 2) ** 2 + self.sig_zxy ** 2)
+        self.sig_1 = self.sig_zz / 2 + np.sqrt(
+            (self.sig_zz / 2) ** 2 + self.sig_zxy ** 2
+        )
+        self.sig_3 = self.sig_zz / 2 - np.sqrt(
+            (self.sig_zz / 2) ** 2 + self.sig_zxy ** 2
+        )
         self.sig_vm = (self.sig_zz ** 2 + 3 * self.sig_zxy ** 2) ** 0.5
 
 
@@ -4588,7 +4767,7 @@ def calculate_centroidal_properties(self, mesh):
         self.ixy_c = self.ixy_g - self.qx * self.qy / self.ea
 
         # calculate section moduli about the centroidal xy axis
-        nodes = np.array(mesh.points)
+        nodes = np.array(mesh["vertices"])
         xmax = nodes[:, 0].max()
         xmin = nodes[:, 0].min()
         ymax = nodes[:, 1].max()
diff --git a/sectionproperties/analysis/fea.py b/sectionproperties/analysis/fea.py
index 9dd75bf5..b02baf9d 100644
--- a/sectionproperties/analysis/fea.py
+++ b/sectionproperties/analysis/fea.py
@@ -88,7 +88,7 @@ def geometric_properties(self):
             ixy,
             self.material.elastic_modulus,
             self.material.shear_modulus,
-            self.material.density
+            self.material.density,
         )
 
     def torsion_properties(self):
diff --git a/sectionproperties/post/post.py b/sectionproperties/post/post.py
index 42235841..56fd2fdb 100644
--- a/sectionproperties/post/post.py
+++ b/sectionproperties/post/post.py
@@ -28,6 +28,8 @@ def finish_plot(ax, pause, title="", size=500, dpi=96):
     :param str title: Plot title
     """
     ax.set_aspect("equal", anchor="C")
+    ax.set_title(title, fontdict={"fontsize": 4 * size / dpi})
+    ax.tick_params(labelsize=2.5 * size / dpi)
     plt.tight_layout()
 
     fig = plt.gcf()
@@ -41,8 +43,6 @@ def finish_plot(ax, pause, title="", size=500, dpi=96):
 
     fig.set_size_inches(new_size)
     fig.set_dpi(dpi)
-    ax.set_title(title, fontdict={"fontsize": 4*size/dpi})
-    ax.tick_params(labelsize=2.5*size/dpi)
 
     plt.show()
 
diff --git a/sectionproperties/pre/bisect_section.py b/sectionproperties/pre/bisect_section.py
index b63b0de1..a221230e 100644
--- a/sectionproperties/pre/bisect_section.py
+++ b/sectionproperties/pre/bisect_section.py
@@ -35,7 +35,8 @@ def create_line_segment(
 
 
 def group_top_and_bottom_polys(
-    polys: shapely.geometry.GeometryCollection, line: shapely.geometry.LineString,
+    polys: shapely.geometry.GeometryCollection,
+    line: shapely.geometry.LineString,
 ) -> Tuple[list, list]:
     """
     Returns tuple of two lists representing the list of Polygons in 'polys' on the "top" side of 'line' and the
@@ -66,7 +67,9 @@ def group_top_and_bottom_polys(
     return top_acc, bot_acc
 
 
-def line_mx_plus_b(line: shapely.geometry.LineString,) -> Tuple[float, float]:
+def line_mx_plus_b(
+    line: shapely.geometry.LineString,
+) -> Tuple[float, float]:
     """
     Returns a tuple representing the values of "m" and "b" from the definition of 'line' as "y = mx + b".
     """
@@ -84,7 +87,8 @@ def line_mx_plus_b(line: shapely.geometry.LineString,) -> Tuple[float, float]:
 
 
 def perp_mx_plus_b(
-    m_slope: float, point_on_line: Tuple[float, float],
+    m_slope: float,
+    point_on_line: Tuple[float, float],
 ) -> Tuple[float, float]:
     """
     Returns a tuple representing the values of "m" and "b" from
@@ -98,7 +102,10 @@ def perp_mx_plus_b(
 
 
 def line_intersection(
-    m_1: float, b_1: float, m_2: float, b_2: float,
+    m_1: float,
+    b_1: float,
+    m_2: float,
+    b_2: float,
 ) -> Optional[float]:
     """
     Returns a float representing the x-ordinate of the intersection
@@ -113,7 +120,9 @@ def line_intersection(
     return x
 
 
-def sum_poly_areas(lop: List[shapely.geometry.Polygon],) -> float:
+def sum_poly_areas(
+    lop: List[shapely.geometry.Polygon],
+) -> float:
     """
     Returns a float representing the total area of all polygons
     in 'lop', the list of polygons.
diff --git a/sectionproperties/pre/dxf.py b/sectionproperties/pre/dxf.py
index 9b7a4502..50df83a1 100644
--- a/sectionproperties/pre/dxf.py
+++ b/sectionproperties/pre/dxf.py
@@ -6,9 +6,9 @@
 
 def load_dxf(dxf_filepath: pathlib.Path):
     """
-            Import any-old-shape in dxf format for analysis.
-            Code by aegis1980 and connorferster
-        """
+    Import any-old-shape in dxf format for analysis.
+    Code by aegis1980 and connorferster
+    """
     if not dxf_filepath.exists():
         raise ValueError(f"The filepath does not exist: {dxf_filepath}")
 
diff --git a/sectionproperties/pre/nastran_sections.py b/sectionproperties/pre/nastran_sections.py
index d618aa5a..023e362a 100644
--- a/sectionproperties/pre/nastran_sections.py
+++ b/sectionproperties/pre/nastran_sections.py
@@ -5,7 +5,9 @@
 from sectionproperties.pre.sections import draw_radius
 
 
-def nastran_bar(DIM1: float, DIM2: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_bar(
+    DIM1: float, DIM2: float, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a BAR section with the center at the origin *(0, 0)*, with two parameters
     defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ [5]_ for definition of
     parameters. Added by JohnDN90.
@@ -49,8 +51,14 @@ def nastran_bar(DIM1: float, DIM2: float, material: pre.Material = pre.DEFAULT_M
     return geometry
 
 
-def nastran_box(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs a BOX section with the center at the origin *(0, 0)*, with four parameters
+def nastran_box(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
+    """Constructs a BOX section with the center at the origin *(0, 0)*, with four parameters
     defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ [5]_ for definition of
     parameters. Added by JohnDN90.
 
@@ -104,7 +112,13 @@ def nastran_box(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pr
 
 
 def nastran_box1(
-    DIM1: float, DIM2: float, DIM3: float, DIM4: float, DIM5: float, DIM6: float, material: pre.Material = pre.DEFAULT_MATERIAL
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    DIM5: float,
+    DIM6: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
 ) -> Geometry:
     """Constructs a BOX1 section with the center at the origin *(0, 0)*, with six parameters
     defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more details. Added by
@@ -166,8 +180,14 @@ def nastran_box1(
     return geometry
 
 
-def nastran_chan(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs a CHAN (C-Channel) section with the web's middle center at the origin *(0, 0)*,
+def nastran_chan(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
+    """Constructs a CHAN (C-Channel) section with the web's middle center at the origin *(0, 0)*,
     with four parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for
     more details. Added by JohnDN90.
 
@@ -224,8 +244,14 @@ def nastran_chan(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: p
     return geometry
 
 
-def nastran_chan1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs a CHAN1 (C-Channel) section with the web's middle center at the origin *(0, 0)*,
+def nastran_chan1(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
+    """Constructs a CHAN1 (C-Channel) section with the web's middle center at the origin *(0, 0)*,
     with four parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for
     more details. Added by JohnDN90.
 
@@ -280,8 +306,14 @@ def nastran_chan1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material:
     return geometry
 
 
-def nastran_chan2(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs a CHAN2 (C-Channel) section with the bottom web's middle center at the origin
+def nastran_chan2(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
+    """Constructs a CHAN2 (C-Channel) section with the bottom web's middle center at the origin
     *(0, 0)*, with four parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_
     [4]_ for more details. Added by JohnDN90.
 
@@ -336,8 +368,14 @@ def nastran_chan2(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material:
     return geometry
 
 
-def nastran_cross(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs Nastran's cruciform/cross section with the intersection's middle center at the
+def nastran_cross(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
+    """Constructs Nastran's cruciform/cross section with the intersection's middle center at the
     origin *(0, 0)*, with four parameters defining dimensions. See Nastran documentation [1]_ [2]_
     [3]_ [4]_ for more details. Added by JohnDN90.
 
@@ -405,9 +443,9 @@ def nastran_fcross(
     DIM6: float,
     DIM7: float,
     DIM8: float,
-    material: pre.Material = pre.DEFAULT_MATERIAL
+    material: pre.Material = pre.DEFAULT_MATERIAL,
 ) -> Geometry:
-    """ Constructs a flanged cruciform/cross section with the intersection's middle center at the
+    """Constructs a flanged cruciform/cross section with the intersection's middle center at the
     origin *(0, 0)*, with eight parameters defining dimensions. Added by JohnDN90.
 
     :param float DIM1: Depth (y) of flanged cruciform
@@ -502,9 +540,9 @@ def nastran_dbox(
     DIM8: float,
     DIM9: float,
     DIM10: float,
-    material: pre.Material = pre.DEFAULT_MATERIAL
+    material: pre.Material = pre.DEFAULT_MATERIAL,
 ) -> Geometry:
-    """ Constructs a DBOX section with the center at the origin *(0, 0)*, with ten parameters
+    """Constructs a DBOX section with the center at the origin *(0, 0)*, with ten parameters
     defining dimensions. See MSC Nastran documentation [1]_ for more details. Added by JohnDN90.
 
     :param float DIM1: Width (x) of the DBOX-section
@@ -581,9 +619,15 @@ def nastran_dbox(
 
 
 def nastran_gbox(
-    DIM1: float, DIM2: float, DIM3: float, DIM4: float, DIM5: float, DIM6: float, material: pre.Material = pre.DEFAULT_MATERIAL
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    DIM5: float,
+    DIM6: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
 ) -> Geometry:
-    """ Constructs a GBOX section with the center at the origin *(0, 0)*, with six parameters
+    """Constructs a GBOX section with the center at the origin *(0, 0)*, with six parameters
     defining dimensions. See ASTROS documentation [5]_ for more details. Added by JohnDN90.
 
     :param float DIM1: Width (x) of the GBOX-section
@@ -653,7 +697,13 @@ def nastran_gbox(
     return geometry
 
 
-def nastran_h(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_h(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a H section with the middle web's middle center at the origin *(0, 0)*, with four
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more details.
     Added by JohnDN90.
@@ -714,7 +764,13 @@ def nastran_h(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.
     return geometry
 
 
-def nastran_hat(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_hat(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a Hat section with the top most section's middle center at the origin *(0, 0)*,
     with four parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for
     more details. Note that HAT in ASTROS is actually HAT1 in this code. Added by JohnDN90.
@@ -774,8 +830,15 @@ def nastran_hat(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pr
     return geometry
 
 
-def nastran_hat1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, DIM5: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs a HAT1 section with the bottom plate's bottom center at the origin *(0, 0)*,
+def nastran_hat1(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    DIM5: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
+    """Constructs a HAT1 section with the bottom plate's bottom center at the origin *(0, 0)*,
     with five parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [5]_ for
     definition of parameters. Note that in ASTROS, HAT1 is called HAT. Added by JohnDN90.
 
@@ -821,7 +884,9 @@ def nastran_hat1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, DIM5: float
 
     hat = nastran_hat(DIM1=d1, DIM2=d2, DIM3=d3, DIM4=d4)
     # Merge the two sections into one geometry
-    geometry = hat.align_center(bottom_plate).align_to(bottom_plate, on="top") + bottom_plate
+    geometry = (
+        hat.align_center(bottom_plate).align_to(bottom_plate, on="top") + bottom_plate
+    )
 
     C = (-0.5 * DIM1, 0)
     D = (0.5 * DIM1, 0)
@@ -834,8 +899,10 @@ def nastran_hat1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, DIM5: float
     return geometry
 
 
-def nastran_hexa(DIM1: float, DIM2: float, DIM3: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
-    """ Constructs a HEXA (hexagon) section with the center at the origin *(0, 0)*, with three
+def nastran_hexa(
+    DIM1: float, DIM2: float, DIM3: float, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
+    """Constructs a HEXA (hexagon) section with the center at the origin *(0, 0)*, with three
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more details.
     Added by JohnDN90.
 
@@ -888,7 +955,13 @@ def nastran_hexa(DIM1: float, DIM2: float, DIM3: float, material: pre.Material =
 
 
 def nastran_i(
-    DIM1: float, DIM2: float, DIM3: float, DIM4: float, DIM5: float, DIM6: float, material: pre.Material = pre.DEFAULT_MATERIAL
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    DIM5: float,
+    DIM6: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
 ) -> Geometry:
     """Constructs Nastran's I section with the bottom flange's middle center at the origin
     *(0, 0)*, with six parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_
@@ -957,7 +1030,13 @@ def nastran_i(
     return geometry
 
 
-def nastran_i1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_i1(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a I1 section with the web's middle center at the origin *(0, 0)*, with four
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more details.
     Added by JohnDN90.
@@ -1019,7 +1098,13 @@ def nastran_i1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre
     return geometry
 
 
-def nastran_l(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_l(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a L section with the intersection's center at the origin *(0, 0)*, with four
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ for more details.
     Added by JohnDN90.
@@ -1068,7 +1153,9 @@ def nastran_l(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.
     return geometry
 
 
-def nastran_rod(DIM1: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_rod(
+    DIM1: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a circular rod section with the center at the origin *(0, 0)*, with one parameter
     defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more details. Added by
     JohnDN90.
@@ -1123,7 +1210,13 @@ def nastran_rod(DIM1: float, n: int, material: pre.Material = pre.DEFAULT_MATERI
     return geometry
 
 
-def nastran_tee(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_tee(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a T section with the top flange's middle center at the origin *(0, 0)*, with four
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ [5]_ for more
     details. Added by JohnDN90.
@@ -1196,7 +1289,13 @@ def nastran_tee(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pr
     return geometry
 
 
-def nastran_tee1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_tee1(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a T1 section with the right flange's middle center at the origin *(0, 0)*, with
     four parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more
     details. Added by JohnDN90.
@@ -1254,7 +1353,13 @@ def nastran_tee1(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: p
     return geometry
 
 
-def nastran_tee2(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_tee2(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a T2 section with the bottom flange's middle center at the origin *(0, 0)*, with
     four parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more
     details. Added by JohnDN90.
@@ -1312,7 +1417,9 @@ def nastran_tee2(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: p
     return geometry
 
 
-def nastran_tube(DIM1: float, DIM2: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_tube(
+    DIM1: float, DIM2: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a circular tube section with the center at the origin *(0, 0)*, with two
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more
     details. Added by JohnDN90.
@@ -1378,7 +1485,9 @@ def nastran_tube(DIM1: float, DIM2: float, n: int, material: pre.Material = pre.
     return geometry
 
 
-def nastran_tube2(DIM1: float, DIM2: float, n: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_tube2(
+    DIM1: float, DIM2: float, n: float, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a circular TUBE2 section with the center at the origin *(0, 0)*, with two
     parameters defining dimensions. See MSC Nastran documentation [1]_ for more details. Added by
     JohnDN90.
@@ -1445,7 +1554,13 @@ def nastran_tube2(DIM1: float, DIM2: float, n: float, material: pre.Material = p
         return geometry
 
 
-def nastran_zed(DIM1: float, DIM2: float, DIM3: float, DIM4: float, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def nastran_zed(
+    DIM1: float,
+    DIM2: float,
+    DIM3: float,
+    DIM4: float,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a Z section with the web's middle center at the origin *(0, 0)*, with four
     parameters defining dimensions. See Nastran documentation [1]_ [2]_ [3]_ [4]_ for more details.
     Added by JohnDN90.
diff --git a/sectionproperties/pre/pre.py b/sectionproperties/pre/pre.py
index 92a69c1e..aca63a7a 100644
--- a/sectionproperties/pre/pre.py
+++ b/sectionproperties/pre/pre.py
@@ -1,8 +1,7 @@
 from typing import Union, List
 from dataclasses import dataclass
 import numpy as np
-import meshpy.triangle as triangle
-import meshpy
+import triangle
 
 
 class GeometryError(Exception):
@@ -76,9 +75,8 @@ def create_mesh(
     holes: List[List[float]],
     control_points: List[List[float]],
     mesh_sizes: Union[List[float], float],
-    atol=1.0e-8,
-) -> meshpy.triangle.MeshInfo:
-    """Creates a quadratic triangular mesh using the meshpy module, which utilises the code
+):
+    """Creates a quadratic triangular mesh using the triangle module, which utilises the code
     'Triangle', by Jonathan Shewchuk.
 
     :param points: List of points *(x, y)* defining the vertices of the cross-section
@@ -93,35 +91,29 @@ def create_mesh(
     :type control_points: list[list[float, float]]
     :param mesh_sizes: List of maximum element areas for each region defined by a control point
     :type mesh_sizes: list[float]
-    :param atol: minimum permissable point distance from any section facet
-    :type atol: float
 
-    :return: Object containing generated mesh data
-    :rtype: :class:`meshpy.triangle.MeshInfo`
+    :return: Dictionary containing mesh data
+    :rtype: dict()
     """
-    # check_geometry(points, facets, holes, control_points, atol=atol)
     if not isinstance(mesh_sizes, list):
         mesh_sizes = [mesh_sizes]
-    mesh = triangle.MeshInfo()  # create mesh info object
-    mesh.set_points(points)  # set points
-    mesh.set_facets(facets)  # set facets
-    mesh.set_holes(holes)  # set holes
 
-    # set regions
-    mesh.regions.resize(len(control_points))  # resize regions list
-    region_id = 0  # initialise region ID variable
+    tri = {}  # create tri dictionary
+    tri["vertices"] = points  # set point
+    tri["segments"] = facets  # set facets
+
+    if holes:
+        tri["holes"] = holes  # set holes
+
+    # prepare regions
+    regions = []
 
     for (i, cp) in enumerate(control_points):
-        mesh.regions[i] = [cp[0], cp[1], region_id, mesh_sizes[i]]
-        region_id += 1
-
-    mesh = triangle.build(
-        mesh,
-        min_angle=30,
-        mesh_order=2,
-        quality_meshing=True,
-        attributes=True,
-        volume_constraints=True,
-    )
+        regions.append([cp[0], cp[1], i, mesh_sizes[i]])
+
+    tri["regions"] = regions  # set regions
+
+    # generate mesh
+    mesh = triangle.triangulate(tri, "pq30Aao2")
 
     return mesh
diff --git a/sectionproperties/pre/rhino.py b/sectionproperties/pre/rhino.py
index 5c472e4d..a1043d83 100644
--- a/sectionproperties/pre/rhino.py
+++ b/sectionproperties/pre/rhino.py
@@ -3,19 +3,20 @@
 from shapely.geometry.polygon import Polygon
 from rhino_shapely_interop.importers import RhImporter
 
+
 def load_3dm(r3dm_filepath: Union[pathlib.Path, str], **kwargs) -> List[Polygon]:
     """Load a Rhino `.3dm` file and import the single surface planer breps.
 
-    :param r3dm_filepath: 
+    :param r3dm_filepath:
         File path to the rhino `.3dm` file.
     :type r3dm_filepath: pathlib.Path or string
     :param \**kwargs:
         See below.
-    :raises RuntimeError: 
-        A RuntimeError is raised if no polygons are found in the file. 
-        This is dependent on the keyword arguments. 
-        Try adjusting the keyword arguments if this error is raised. 
-    :return: 
+    :raises RuntimeError:
+        A RuntimeError is raised if no polygons are found in the file.
+        This is dependent on the keyword arguments.
+        Try adjusting the keyword arguments if this error is raised.
+    :return:
         List of Polygons found in the file.
     :rtype: List[shapely.geometry.Polygon]
 
@@ -23,7 +24,7 @@ def load_3dm(r3dm_filepath: Union[pathlib.Path, str], **kwargs) -> List[Polygon]
         * *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. 
+            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.
@@ -48,27 +49,28 @@ def load_3dm(r3dm_filepath: Union[pathlib.Path, str], **kwargs) -> List[Polygon]
     """
     rhi = RhImporter.from_file(str(r3dm_filepath))
     list_polygons = list(rhi.get_planer_brep(**kwargs))
-    if len(list_polygons)==0:
+    if len(list_polygons) == 0:
         raise RuntimeError(
             f"No shapely.Polygon objects found. "
             f"Consider adjusting the keyword arguments. "
             f"File name: {r3dm_filepath}. "
-            )
+        )
     return list_polygons
 
+
 def load_brep_encoding(brep: str, **kwargs) -> Polygon:
     """Load an encoded single surface planer brep.
 
-    :param brep: 
+    :param brep:
         Rhino3dm.Brep encoded as a string.
     :type brep: str
     :param \**kwargs:
         See below.
-    :raises RuntimeError: 
-        A RuntimeError is raised if no polygons are found in the encoding. 
-        This is dependent on the keyword arguments. 
-        Try adjusting the keyword arguments if this error is raised. 
-    :return: 
+    :raises RuntimeError:
+        A RuntimeError is raised if no polygons are found in the encoding.
+        This is dependent on the keyword arguments.
+        Try adjusting the keyword arguments if this error is raised.
+    :return:
         The Polygons found in the encoding string.
     :rtype: shapely.geometry.Polygon
 
@@ -76,7 +78,7 @@ def load_brep_encoding(brep: str, **kwargs) -> Polygon:
         * *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. 
+            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.
@@ -98,9 +100,9 @@ def load_brep_encoding(brep: str, **kwargs) -> Polygon:
             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.
-    """    
+    """
     rhi = RhImporter.from_serialzed_brep(brep)
     geom = list(rhi.get_planer_brep(**kwargs))
-    if len(geom)==0:
+    if len(geom) == 0:
         raise RuntimeError(f"No shapely.Polygon objects found for encoded object")
-    return geom
\ No newline at end of file
+    return geom
diff --git a/sectionproperties/pre/sections.py b/sectionproperties/pre/sections.py
index 3c1b1656..197ceffa 100644
--- a/sectionproperties/pre/sections.py
+++ b/sectionproperties/pre/sections.py
@@ -1,5 +1,6 @@
 from __future__ import annotations
-from typing import List, Optional, Union, Tuple
+from ntpath import join
+from typing import List, Optional, Union, Tuple, Any
 
 import copy
 import pathlib
@@ -8,9 +9,11 @@
 from shapely.geometry import (
     Polygon,
     MultiPolygon,
+    LineString,
     LinearRing,
     Point,
     box,
+    GeometryCollection,
 )
 from shapely.ops import split, unary_union
 import shapely
@@ -21,6 +24,7 @@
 import sectionproperties.post.post as post
 import sectionproperties.pre.rhino as rhino_importer
 
+
 class Geometry:
     """Class for defining the geometry of a contiguous section of a single material.
 
@@ -47,13 +51,12 @@ def __init__(
         material: pre.Material = pre.DEFAULT_MATERIAL,
         control_points: Optional[List[float, float]] = None,
     ):
-        """Inits the Geometry class.
-        """
+        """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: {geom}"
+                f"Argument is not a valid shapely.geometry.Polygon object: {repr(geom)}"
             )
         self.assigned_control_point = None
         if control_points is not None and len(control_points) == 2:
@@ -67,7 +70,9 @@ def __init__(
         self.holes = []
         self.perimeter = []
         self._recovery_points = []
-        self.tol = 12 # Represents num of decimal places of precision for point locations
+        self.tol = (
+            12  # Represents num of decimal places of precision for point locations
+        )
         self.compile_geometry()
 
     def _repr_svg_(self):
@@ -143,7 +148,7 @@ def from_points(
                 "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 = []
@@ -189,11 +194,15 @@ def from_points(
         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)
+        geometry = Geometry(
+            exterior_geometry - interior_geometry, material, control_points
+        )
         return geometry
 
     @staticmethod
-    def from_dxf(dxf_filepath: Union[str, pathlib.Path]) -> Union[Geometry, CompoundGeometry]:
+    def from_dxf(
+        dxf_filepath: Union[str, pathlib.Path]
+    ) -> Union[Geometry, CompoundGeometry]:
         """
         An interface for the creation of Geometry objects from CAD .dxf files.
 
@@ -247,8 +256,8 @@ def from_3dm(cls, filepath: Union[str, pathlib.Path], **kwargs) -> Geometry:
                 Default is False.
         """
         geom = None
-        list_poly = rhino_importer.load_3dm(filepath,**kwargs)
-        if (len(list_poly)==1):
+        list_poly = rhino_importer.load_3dm(filepath, **kwargs)
+        if len(list_poly) == 1:
             geom = cls(geom=list_poly[0])
         else:
             raise RuntimeError(
@@ -297,7 +306,7 @@ def from_rhino_encoding(cls, r3dm_brep: str, **kwargs) -> Geometry:
                 If true, all non parallel surfaces are filtered out.
                 Default is False.
         """
-        return cls(geom = rhino_importer.load_brep_encoding(r3dm_brep, **kwargs)[0])
+        return cls(geom=rhino_importer.load_brep_encoding(r3dm_brep, **kwargs)[0])
 
     def create_facets_and_control_points(self):
         self.perimeter = None
@@ -323,7 +332,6 @@ def compile_geometry(self):  # Alias
         """
         self.create_facets_and_control_points()
 
-
     def create_mesh(self, mesh_sizes: Union[float, List[float]]):
         """Creates a quadratic triangular mesh from the Geometry object.
 
@@ -364,7 +372,10 @@ def create_mesh(self, mesh_sizes: Union[float, List[float]]):
         return self
 
     def align_to(
-        self, other: Union[Geometry, Tuple[float, float]], on: str, inner: bool = False,
+        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
@@ -470,7 +481,9 @@ def align_center(self, align_to: Optional[Geometry] = None):
         return new_geom
 
     def shift_section(
-        self, x_offset=0.0, y_offset=0.0,
+        self,
+        x_offset=0.0,
+        y_offset=0.0,
     ):
         """
         Returns a new Geometry object translated by 'x_offset' and 'y_offset'.
@@ -486,9 +499,13 @@ def shift_section(
         # 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_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
+            shapely.affinity.translate(self.geom, x_offset, y_offset),
+            self.material,
+            new_ctrl_point,
         )
         return new_geom
 
@@ -525,10 +542,11 @@ def rotate_section(
                 rotate_point = box(*self.geom.bounds).centroid
             new_ctrl_point = shapely.affinity.rotate(
                 self.assigned_control_point, angle, rotate_point, use_radians
-                ).coords[0]
+            ).coords[0]
         new_geom = Geometry(
             shapely.affinity.rotate(self.geom, angle, rot_point, use_radians),
-            self.material, new_ctrl_point
+            self.material,
+            new_ctrl_point,
         )
         return new_geom
 
@@ -558,10 +576,10 @@ def mirror_section(
         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
+        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
         )
@@ -569,7 +587,11 @@ def mirror_section(
         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
+                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
@@ -616,7 +638,12 @@ def split_section(
             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],])
+            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:
@@ -637,7 +664,9 @@ def split_section(
 
         return (top_right_geoms, bottom_left_geoms)
 
-    def offset_perimeter(self, amount: float = 0, where: str = "exterior", resolution: float = 12):
+    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.
@@ -653,47 +682,59 @@ def offset_perimeter(self, amount: float = 0, where: str = "exterior", resolutio
         :return: Geometry object translated to new alignment
         :rtype: :class:`sections.pre.sections.Geometry`
 
-        The following example erodes a 200PFC section by 3::
+        The following example erodes a 200PFC section by 2 mm::
 
             import sectionproperties.pre.sections as 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=-3)
+            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 = Polygon(interior).buffer(distance=amount, join_style=1, resolution=resolution)
+                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])
+                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):
+            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.")
+            raise ValueError(
+                "Cannot buffer interior of Geometry object if it has no holes."
+            )
 
         elif where == "exterior":
             exterior_polygon = Polygon(self.geom.exterior)
-            buffered_exterior = exterior_polygon.buffer(distance=amount, join_style=1, resolution=resolution)
+            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])
+                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])
+            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 = self.geom.buffer(
-                distance=amount, join_style=1, resolution=resolution
-            )
+            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]
@@ -701,8 +742,12 @@ def offset_perimeter(self, amount: float = 0, where: str = "exterior", resolutio
                 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])
+            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)
 
@@ -758,7 +803,7 @@ def shift_points(
         current_points = copy.copy(self.points)
         current_facets = copy.copy(self.facets)
         current_holes = copy.copy(self.holes)
-        current_materials = copy.copy(self.materials)
+        current_material = copy.copy(self.material)
         current_control_points = copy.copy(self.control_points)
 
         if isinstance(point_idxs, int):
@@ -771,15 +816,32 @@ def shift_points(
             current_points[point_idx] = (new_x, new_y)
 
         new_geom = Geometry.from_points(
-            current_points, current_facets, current_holes, control_points=None, materials=current_materials
+            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):
+        if self.assigned_control_point and new_geom.geom.contains(
+            self.assigned_control_point
+        ):
             new_geom = Geometry.from_points(
-                current_points, current_facets, current_holes, current_control_points[0], current_materials
+                current_points,
+                current_facets,
+                current_control_points[0],
+                current_holes,
+                current_material,
             )
         return new_geom
 
-    def plot_geometry(self, ax: Optional[matplotlib.axes.Axes]=None, labels=["control_points"], size=500, dpi=96, pause=True):
+    def plot_geometry(
+        self,
+        ax: Optional[matplotlib.axes.Axes] = None,
+        labels=["control_points"],
+        size=500,
+        dpi=96,
+        pause=True,
+    ):
         """Plots the geometry defined by the input section. If no axes object is supplied a new
         figure and axis is created. Allows for 'size' and 'dpi' arguments to be given however the
         function dynamically sizes nodes and lineweights so that they appear appropriately sized at
@@ -818,7 +880,6 @@ def plot_geometry(self, ax: Optional[matplotlib.axes.Axes]=None, labels=["contro
             Geometry generated by the above example.
         """
         # if no axes object is supplied, create and setup the plot
-        perimeter = None  # To be implemented for CompoundGeometry at a later date
         fig = None
 
         if ax is None:
@@ -829,13 +890,7 @@ def plot_geometry(self, ax: Optional[matplotlib.axes.Axes]=None, labels=["contro
             ax_supplied = True
 
         for (i, f) in enumerate(self.facets):
-            if perimeter:
-                if i in self.perimeter:
-                    linewidth = 0.3*size/dpi
-                else:
-                    linewidth = 0.15*size/dpi
-            else:
-                linewidth = 0.15*size/dpi
+            linewidth = 0.15 * size / dpi
 
             # plot the points and facets
             if i == 0:
@@ -843,7 +898,7 @@ def plot_geometry(self, ax: Optional[matplotlib.axes.Axes]=None, labels=["contro
                     [self.points[f[0]][0], self.points[f[1]][0]],
                     [self.points[f[0]][1], self.points[f[1]][1]],
                     "ko-",
-                    markersize=0.15*size/dpi,
+                    markersize=0.15 * size / dpi,
                     linewidth=linewidth,
                     label="Points & Facets",
                 )
@@ -852,23 +907,42 @@ def plot_geometry(self, ax: Optional[matplotlib.axes.Axes]=None, labels=["contro
                     [self.points[f[0]][0], self.points[f[1]][0]],
                     [self.points[f[0]][1], self.points[f[1]][1]],
                     "ko-",
-                    markersize=0.15*size/dpi,
+                    markersize=0.15 * size / dpi,
                     linewidth=linewidth,
                 )
 
         for (i, h) in enumerate(self.holes):
             # plot the holes
             if i == 0:
-                ax.plot(h[0], h[1], "rx", markersize=1*size/dpi, markeredgewidth=0.2*size/dpi, label="Holes")
+                ax.plot(
+                    h[0],
+                    h[1],
+                    "rx",
+                    markersize=1 * size / dpi,
+                    markeredgewidth=0.2 * size / dpi,
+                    label="Holes",
+                )
             else:
-                ax.plot(h[0], h[1], "rx", markersize=1*size/dpi, markeredgewidth=0.2*size/dpi)
+                ax.plot(
+                    h[0],
+                    h[1],
+                    "rx",
+                    markersize=1 * size / dpi,
+                    markeredgewidth=0.2 * size / dpi,
+                )
 
         for (i, cp) in enumerate(self.control_points):
             # plot the control points
             if i == 0:
-                ax.plot(cp[0], cp[1], "bo", markersize=1*size/dpi, label="Control Points")
+                ax.plot(
+                    cp[0],
+                    cp[1],
+                    "bo",
+                    markersize=1 * size / dpi,
+                    label="Control Points",
+                )
             else:
-                ax.plot(cp[0], cp[1], "bo", markersize=1*size/dpi)
+                ax.plot(cp[0], cp[1], "bo", markersize=1 * size / dpi)
 
         # display the labels
         for label in labels:
@@ -891,11 +965,13 @@ def plot_geometry(self, ax: Optional[matplotlib.axes.Axes]=None, labels=["contro
                     ax.annotate(str(i), xy=xy, color="b")
 
         # display the legend
-        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+        ax.legend(loc="center left", bbox_to_anchor=(1, 0.5))
 
         # if no axes object is supplied, finish the plot
         if not ax_supplied:
-            post.finish_plot(ax, pause, title="Cross-Section Geometry", size=size, dpi=dpi)
+            post.finish_plot(
+                ax, pause, title="Cross-Section Geometry", size=size, dpi=dpi
+            )
             return (fig, ax)
         return (fig, ax)
 
@@ -934,7 +1010,7 @@ def calculate_centroid(self):
         :return: Geometry centroid.
         :rtype: Tuple[float, float]
         """
-        cx, cy = self.geom.centroid.coords[0] # Nested list
+        cx, cy = self.geom.centroid.coords[0]  # Nested list
         return (cx, cy)
 
     @property
@@ -964,13 +1040,14 @@ def __or__(self, other):
         """
         material = self.material or other.material
         try:
-            new_polygon = self.geom | other.geom
+            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}"
@@ -982,11 +1059,13 @@ def __xor__(self, other):
         """
         material = self.material or other.material
         try:
-            new_polygon = self.geom ^ other.geom
+            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(
@@ -999,13 +1078,17 @@ def __sub__(self, other):
         """
         material = self.material or other.material
         try:
-            new_polygon = self.geom - other.geom
+            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):
+            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:
@@ -1030,8 +1113,8 @@ def __and__(self, other):
         """
         material = self.material or other.material
         try:
-            new_polygon = self.geom & other.geom
-            if isinstance(new_polygon, MultiPolygon):
+            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]
                 )
@@ -1040,8 +1123,11 @@ def __and__(self, other):
             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.
@@ -1059,9 +1145,12 @@ class CompoundGeometry(Geometry):
         instance.
     :vartype geoms: Union[shapely.geometry.MultiPolygon, List[sectionproperties.pre.sections.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.geoms = [
+                Geometry(poly, material=pre.DEFAULT_MATERIAL) for poly in geoms.geoms
+            ]
             self.geom = geoms
         elif isinstance(geoms, list):
             processed_geoms = []
@@ -1096,7 +1185,6 @@ def _repr_svg_(self):
         print(f"Materials incl.: {list(set(materials_list))}")
         return self.geom._repr_svg_()
 
-
     @staticmethod
     def from_points(
         points: List[List[float]],
@@ -1140,7 +1228,7 @@ def from_points(
             raise ValueError(
                 "Materials cannot be assigned without control_points. "
                 "Please provide corresponding control_points for each material."
-                )
+            )
         if holes is None:
             holes = list()
 
@@ -1167,7 +1255,7 @@ def from_points(
                     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.
+        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))
 
@@ -1176,28 +1264,34 @@ def from_points(
         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))]
-            )
+            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:
             return CompoundGeometry(
-                [Geometry(exterior, control_points=control_points[idx]) for idx, exterior in enumerate(exteriors)]
-                )
+                [
+                    Geometry(exterior, control_points=control_points[idx])
+                    for idx, exterior in enumerate(exteriors)
+                ]
+            )
         else:
             # "Punch" all holes through each exterior geometry
             punched_exteriors = []
@@ -1208,15 +1302,18 @@ def from_points(
                     punched_exterior = punched_exterior - interior
                     try:
                         exterior_control_point = next(
-                            control_point for control_point in control_points
+                            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}"
-                            )
-                exterior_geometry = Geometry(punched_exterior, control_points=exterior_control_point)
+                        )
+                exterior_geometry = Geometry(
+                    punched_exterior, control_points=exterior_control_point
+                )
                 punched_exterior_geometries.append(exterior_geometry)
 
         return CompoundGeometry(punched_exterior_geometries)
@@ -1262,7 +1359,7 @@ def from_3dm(cls, filepath: Union[str, pathlib.Path], **kwargs) -> CompoundGeome
                 Default is False.
         """
         list_poly = rhino_importer.load_3dm(filepath, **kwargs)
-        return cls(geoms = MultiPolygon(list_poly))
+        return cls(geoms=MultiPolygon(list_poly))
 
     def create_mesh(self, mesh_sizes: List[float]):
         """Creates a quadratic triangular mesh from the Geometry object.
@@ -1347,7 +1444,7 @@ def rotate_section(
             new_compound = compound.rotate_section(angle=-30)
         """
         geoms_acc = []
-        if rot_point =='center':
+        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))
@@ -1383,7 +1480,6 @@ def mirror_section(
         new_geom = CompoundGeometry(geoms_acc)
         return new_geom
 
-
     def align_center(self, align_to: Optional[Geometry] = None):
         """
         Returns a new CompoundGeometry object, translated in both x and y, so that the
@@ -1402,21 +1498,29 @@ def align_center(self, align_to: Optional[Geometry] = None):
         :return: Geometry object translated to new alignment
         :rtype: :class:`sections.pre.sections.Geometry`
         """
-        EA_sum = sum([geom.material.elastic_modulus * geom.calculate_area() for geom in self.geoms])
+        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
+            EA = E * A
             cx, cy = list(geom.geom.centroid.coords[0])
-            cx_EA_acc += cx*EA
-            cy_EA_acc += cy*EA
+            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)
+            shift_x, shift_y = (
+                round(-weighted_cx, self.tol),
+                round(-weighted_cy, self.tol),
+            )
 
         else:
             align_cx, align_cy = list(align_to.geom.centroid.coords)[0]
@@ -1425,7 +1529,6 @@ def align_center(self, align_to: Optional[Geometry] = None):
         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],
@@ -1475,40 +1578,78 @@ def split_section(
             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 perimeter of the individual geometries within the CompoundGeometry
-        object by a discrete amount. Note, because the individual geometries have their own
-        perimeters offset independently, sections don't "stick" as though they were a joined section.
-        Any aligned geometries within the CompoundGeometry will need to be re-aligned after
-        the perimeter offset.
+    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 resolution: Number of segments used to approximate a quarter circle around a point
-        :type resolution: 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:`sections.pre.sections.Geometry`
 
-        The following example erodes a 200PFC section by 3::
+        The following example erodes a 200UB25 with a 12 plate stiffener section by 2 mm::
 
             import sectionproperties.pre.sections as sections
 
             geometry_1 = sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
-            geometry_2 = sections.rectangular_section(d=20, b=133)
+            geometry_2 = 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_section_perimeter(amount=-3)
+            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.
         """
-        geoms_acc = []
-        for geom in self.geoms:
-            geoms_acc.append(geom.offset_perimeter(amount, where, resolution))
-        new_geom = CompoundGeometry(geoms_acc)
-        return new_geom
+        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):
         """
@@ -1586,6 +1727,7 @@ def calculate_perimeter(self) -> float:
 
 ### Helper functions for Geometry
 
+
 def load_dxf(dxf_filepath: pathlib.Path):
     """
     Import any-old-shape in dxf format for analysis.
@@ -1694,6 +1836,54 @@ def create_points_and_facets(shape: Polygon, tol=12) -> tuple:
     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 draw_radius(
     pt: list, r: float, theta: float, n, ccw: bool = True
 ):  # Changed 'anti' to ccw to match shapely
@@ -1728,7 +1918,9 @@ def draw_radius(
     return points
 
 
-def rectangular_section(b, d, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def rectangular_section(
+    b, d, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a rectangular section with the bottom left corner at the origin *(0, 0)*, with
     depth *d* and width *b*.
 
@@ -1761,7 +1953,9 @@ def rectangular_section(b, d, material: pre.Material = pre.DEFAULT_MATERIAL) ->
     return Geometry(rectangle, material)
 
 
-def circular_section(d: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def circular_section(
+    d: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a solid circle centered at the origin *(0, 0)* with diameter *d* and using *n*
     points to construct the circle.
 
@@ -1807,7 +2001,9 @@ def circular_section(d: float, n: int, material: pre.Material = pre.DEFAULT_MATE
     return Geometry(circle, material)
 
 
-def circular_hollow_section(d: float, t: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def circular_hollow_section(
+    d: float, t: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a circular hollow section (CHS) centered at the origin *(0, 0)*, with diameter *d* and
     thickness *t*, using *n* points to construct the inner and outer circles.
 
@@ -1858,7 +2054,9 @@ def circular_hollow_section(d: float, t: float, n: int, material: pre.Material =
     return Geometry(outer_circle - inner_circle, material)
 
 
-def elliptical_section(d_y: float, d_x: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def elliptical_section(
+    d_y: float, d_x: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a solid ellipse centered at the origin *(0, 0)* with vertical diameter *d_y* and
     horizontal diameter *d_x*, using *n* points to construct the ellipse.
 
@@ -1906,7 +2104,13 @@ def elliptical_section(d_y: float, d_x: float, n: int, material: pre.Material =
     return Geometry(ellipse, material)
 
 
-def elliptical_hollow_section(d_y: float, d_x: float, t: float, n: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def elliptical_hollow_section(
+    d_y: float,
+    d_x: float,
+    t: float,
+    n: int,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs an elliptical hollow section (EHS) centered at the origin *(0, 0)*, with outer vertical
     diameter *d_y*, outer horizontal diameter *d_x*, and thickness *t*, using *n* points to
     construct the inner and outer ellipses.
@@ -1960,7 +2164,14 @@ def elliptical_hollow_section(d_y: float, d_x: float, t: float, n: int, material
     return Geometry(outer - inner, material)
 
 
-def rectangular_hollow_section(b: float, d: float, t: float, r_out: float, n_r: int, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def rectangular_hollow_section(
+    b: float,
+    d: float,
+    t: float,
+    r_out: float,
+    n_r: int,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
+) -> Geometry:
     """Constructs a rectangular hollow section (RHS) centered at *(b/2, d/2)*, with depth *d*, width *b*,
     thickness *t* and outer radius *r_out*, using *n_r* points to construct the inner and outer
     radii. If the outer radius is less than the thickness of the RHS, the inner radius is set to
@@ -2015,7 +2226,13 @@ def rectangular_hollow_section(b: float, d: float, t: float, r_out: float, n_r:
 
 
 def i_section(
-    d: float, b: float, t_f: float, t_w: float, r: float, n_r: int, material: pre.Material = pre.DEFAULT_MATERIAL
+    d: float,
+    b: float,
+    t_f: float,
+    t_w: float,
+    r: float,
+    n_r: int,
+    material: pre.Material = pre.DEFAULT_MATERIAL,
 ) -> Geometry:  # More specific description and less ambiguous? e.g. not an "S" section.
     """Constructs an I-section centered at *(b/2, d/2)*, with depth *d*, width *b*, flange
     thickness *t_f*, web thickness *t_w*, and root radius *r*, using *n_r* points to construct the
@@ -2085,7 +2302,9 @@ def i_section(
     return Geometry(i_section, material)
 
 
-def mono_i_section(d, b_t, b_b, t_fb, t_ft, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def mono_i_section(
+    d, b_t, b_b, t_fb, t_ft, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a monosymmetric I-section centered at *(max(b_t, b_b)/2, d/2)*, with depth *d*,
     top flange width *b_t*, bottom flange width *b_b*, top flange thickness *t_ft*, top flange
     thickness *t_fb*, web thickness *t_w*, and root radius *r*, using *n_r* points to construct the
@@ -2163,7 +2382,9 @@ def mono_i_section(d, b_t, b_b, t_fb, t_ft, t_w, r, n_r, material: pre.Material
     return Geometry(polygon, material)
 
 
-def tapered_flange_i_section(d, b, t_f, t_w, r_r, r_f, alpha, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def tapered_flange_i_section(
+    d, b, t_f, t_w, r_r, r_f, alpha, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a Tapered Flange I-section centered at *(b/2, d/2)*, with depth *d*, width *b*,
     mid-flange thickness *t_f*, web thickness *t_w*, root radius *r_r*, flange radius *r_f* and
     flange angle *alpha*, using *n_r* points to construct the radii.
@@ -2355,7 +2576,9 @@ def tapered_flange_i_section(d, b, t_f, t_w, r_r, r_f, alpha, n_r, material: pre
     return Geometry(polygon, material)
 
 
-def channel_section(d, b, t_f, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def channel_section(
+    d, b, t_f, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a parallel-flange channel (PFC) section with the bottom left corner at the origin *(0, 0)*, with depth *d*,
     width *b*, flange thickness *t_f*, web  thickness *t_w* and root radius *r*, using *n_r* points
     to construct the root radius.
@@ -2415,7 +2638,9 @@ def channel_section(d, b, t_f, t_w, r, n_r, material: pre.Material = pre.DEFAULT
     return Geometry(polygon, material)
 
 
-def tapered_flange_channel(d, b, t_f, t_w, r_r, r_f, alpha, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def tapered_flange_channel(
+    d, b, t_f, t_w, r_r, r_f, alpha, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a Tapered Flange Channel section with the bottom left corner at the origin
     *(0, 0)*, with depth *d*, width *b*, mid-flange thickness *t_f*, web thickness *t_w*, root
     radius *r_r*, flange radius *r_f* and flange angle *alpha*, using *n_r* points to construct the
@@ -2544,7 +2769,9 @@ def tapered_flange_channel(d, b, t_f, t_w, r_r, r_f, alpha, n_r, material: pre.M
     return Geometry(polygon, material)
 
 
-def tee_section(d, b, t_f, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def tee_section(
+    d, b, t_f, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a Tee section with the top left corner at *(0, d)*, with depth *d*, width *b*,
     flange thickness *t_f*, web thickness *t_w* and root radius *r*, using *n_r* points to
     construct the root radius.
@@ -2601,7 +2828,9 @@ def tee_section(d, b, t_f, t_w, r, n_r, material: pre.Material = pre.DEFAULT_MAT
     return Geometry(polygon, material)
 
 
-def angle_section(d, b, t, r_r, r_t, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def angle_section(
+    d, b, t, r_r, r_t, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs an angle section with the bottom left corner at the origin *(0, 0)*, with depth
     *d*, width *b*, thickness *t*, root radius *r_r* and toe radius *r_t*, using *n_r* points to
     construct the radii.
@@ -2663,7 +2892,9 @@ def angle_section(d, b, t, r_r, r_t, n_r, material: pre.Material = pre.DEFAULT_M
     return Geometry(polygon, material)
 
 
-def cee_section(d, b, l, t, r_out, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def cee_section(
+    d, b, l, t, r_out, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a Cee section (typical of cold-formed steel) with the bottom left corner at the
     origin *(0, 0)*, with depth *d*, width *b*, lip *l*, thickness *t* and outer radius *r_out*,
     using *n_r* points to construct the radius. If the outer radius is less than the thickness
@@ -2744,7 +2975,9 @@ def cee_section(d, b, l, t, r_out, n_r, material: pre.Material = pre.DEFAULT_MAT
     return Geometry(polygon, material)
 
 
-def zed_section(d, b_l, b_r, l, t, r_out, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def zed_section(
+    d, b_l, b_r, l, t, r_out, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a zed section with the bottom left corner at the origin *(0, 0)*, with depth *d*,
     left flange width *b_l*, right flange width *b_r*, lip *l*, thickness *t* and outer radius
     *r_out*, using *n_r* points to construct the radius. If the outer radius is less than the
@@ -2827,7 +3060,9 @@ def zed_section(d, b_l, b_r, l, t, r_out, n_r, material: pre.Material = pre.DEFA
 zee_section = zed_section  # An alias for our American friends (and friends who use "American English")
 
 
-def cruciform_section(d, b, t, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def cruciform_section(
+    d, b, t, r, n_r, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a cruciform section centered at the origin *(0, 0)*, with depth *d*, width *b*,
     thickness *t* and root radius *r*, using *n_r* points to construct the root radius.
 
@@ -2894,7 +3129,9 @@ def cruciform_section(d, b, t, r, n_r, material: pre.Material = pre.DEFAULT_MATE
     return Geometry(polygon, material)
 
 
-def polygon_hollow_section(d, t, n_sides, r_in=0, n_r=1, rot=0, material: pre.Material = pre.DEFAULT_MATERIAL) -> Geometry:
+def polygon_hollow_section(
+    d, t, n_sides, r_in=0, n_r=1, rot=0, material: pre.Material = pre.DEFAULT_MATERIAL
+) -> Geometry:
     """Constructs a regular hollow polygon section centered at *(0, 0)*, with a pitch circle
     diameter of bounding polygon *d*, thickness *t*, number of sides *n_sides* and an optional
     inner radius *r_in*, using *n_r* points to construct the inner and outer radii (if radii is
@@ -3033,7 +3270,9 @@ def rotate(point, angle: float):
     return [new_x, new_y]
 
 
-def box_girder_section(d, b_t, b_b, t_ft, t_fb, t_w, material: pre.Material = pre.DEFAULT_MATERIAL):
+def box_girder_section(
+    d, b_t, b_b, t_ft, t_fb, t_w, material: pre.Material = pre.DEFAULT_MATERIAL
+):
     """Constructs a box girder section centered at at *(max(b_t, b_b)/2, d/2)*, with depth *d*, top
     width *b_t*, bottom width *b_b*, top flange thickness *t_ft*, bottom flange thickness *t_fb*
     and web thickness *t_w*.
diff --git a/setup.py b/setup.py
index 01343fc4..9f2b0ea9 100644
--- a/setup.py
+++ b/setup.py
@@ -4,48 +4,56 @@
 
 
 def description_text():
-    description = 'A python package for the analysis of arbitrary cross-sections using the'
-    description += ' finite element method.'
+    description = (
+        "A python package for the analysis of arbitrary cross-sections using the"
+    )
+    description += " finite element method."
 
     return description
 
 
 def readme():
-    with open('README.md') as f:
+    with open("README.md") as f:
         return f.read()
 
 
 if sys.version_info[0] < 3 or sys.version_info[0] == 3 and sys.version_info[1] < 7:
-    sys.exit('Sorry, Python < 3.7 is not supported')
+    sys.exit("Sorry, Python < 3.7 is not supported")
 
 install_requires = [
-    'numpy', 'scipy', 'pytest_check', 'more_itertools', 'matplotlib', 'shapely', 'pybind11',
-    'meshpy', 'rhino-shapley-interop'
+    "numpy",
+    "scipy",
+    "pytest_check",
+    "more_itertools",
+    "matplotlib",
+    "shapely",
+    "triangle",
+    "rhino-shapley-interop",
+    "cad_to_shapely",
 ]
 
 setup(
-    name='sectionproperties',
+    name="sectionproperties",
     version=version,
     description=description_text(),
     long_description=readme(),
-    long_description_content_type='text/markdown',
+    long_description_content_type="text/markdown",
     classifiers=[
-        'Development Status :: 5 - Production/Stable',
-        'Environment :: Console',
-        'Intended Audience :: Science/Research',
-        'License :: OSI Approved :: MIT License',
-        'Programming Language :: Python :: 3 :: Only',
-        'Topic :: Scientific/Engineering',
+        "Development Status :: 5 - Production/Stable",
+        "Environment :: Console",
+        "Intended Audience :: Science/Research",
+        "License :: OSI Approved :: MIT License",
+        "Programming Language :: Python :: 3 :: Only",
+        "Topic :: Scientific/Engineering",
     ],
-    url='https://github.com/robbievanleeuwen/section-properties',
-    author='Robbie van Leeuwen',
-    author_email='robbie.vanleeuwen@gmail.com',
-    license='MIT',
+    url="https://github.com/robbievanleeuwen/section-properties",
+    author="Robbie van Leeuwen",
+    author_email="robbie.vanleeuwen@gmail.com",
+    license="MIT",
     packages=[
-        'sectionproperties', 'sectionproperties.analysis', 'sectionproperties.post',
-        'sectionproperties.pre', 'sectionproperties.tests'
+        "sectionproperties",
     ],
     install_requires=install_requires,
     include_package_data=True,
-    zip_safe=False
+    zip_safe=False,
 )
diff --git a/todo.md b/todo.md
deleted file mode 100755
index 7d2a9128..00000000
--- a/todo.md
+++ /dev/null
@@ -1,3 +0,0 @@
-# TODO LIST:
-- [ ] add lots of tests!
-- [ ] concrete reinforcement - create circle with slightly larger diameter but correct area??


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 00000000..fd23a827
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,30 @@
+name: Run Tests
+
+on: [workflow_dispatch, push, pull_request]
+
+jobs:
+  tests:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest, windows-latest]
+        python-version: ['3.7', '3.8', '3.9']
+
+    name: Runs tests with python ${{ matrix.python-version }} using ${{ matrix.os }}
+    steps:
+    - uses: actions/checkout@v2
+
+    - name: Setup python ${{ matrix.python-version }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: ${{ matrix.python-version }}
+
+    - name: Update pip
+      run: pip install --upgrade pip
+
+    - name: Install package
+      run: python -m pip install -e .
+
+    - name: Run package tests
+      run: python -m pytest
diff --git a/sectionproperties/tests/test_offset.py b/sectionproperties/tests/test_offset.py
new file mode 100644
index 00000000..b4211a41
--- /dev/null
+++ b/sectionproperties/tests/test_offset.py
@@ -0,0 +1,195 @@
+import pytest_check as check
+import numpy as np
+import sectionproperties.pre.sections as sections
+from sectionproperties.analysis.cross_section import Section
+from shapely.geometry import Polygon
+from sectionproperties.pre.sections import Geometry
+
+
+r_tol = 1e-3
+
+
+def test_rectangular_offset():
+    # exterior negative offset
+    rect = sections.rectangular_section(d=500, b=300)
+    rect = rect.offset_perimeter(amount=-10, where="exterior")
+    rect.create_mesh([200])
+    section = Section(rect)
+    section.calculate_geometric_properties()
+    area = 480 * 280
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+    # exterior positive offset
+    rect = sections.rectangular_section(d=500, b=300)
+    rect = rect.offset_perimeter(amount=10, where="exterior")
+    rect.create_mesh([200])
+    section = Section(rect)
+    section.calculate_geometric_properties()
+    area = 520 * 320 - (20 * 20 - np.pi * 10 * 10)
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+
+def test_box_offset():
+    # exterior negative offset
+    box = sections.rectangular_hollow_section(d=200, b=100, t=10, r_out=0, n_r=1)
+    box = box.offset_perimeter(amount=-5, where="exterior")
+    box.create_mesh([50])
+    section = Section(box)
+    section.calculate_geometric_properties()
+    area = 190 * 90 - 180 * 80
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+    # exterior positve offset
+    box = sections.rectangular_hollow_section(d=200, b=100, t=10, r_out=0, n_r=1)
+    box = box.offset_perimeter(amount=5, where="exterior")
+    box.create_mesh([50])
+    section = Section(box)
+    section.calculate_geometric_properties()
+    area = 210 * 110 - (10 * 10 - np.pi * 5 * 5) - 180 * 80
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+    # interior negative offset
+    box = sections.rectangular_hollow_section(d=200, b=100, t=10, r_out=0, n_r=1)
+    box = box.offset_perimeter(amount=-5, where="interior")
+    box.create_mesh([50])
+    section = Section(box)
+    section.calculate_geometric_properties()
+    area = 200 * 100 - 170 * 70
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+    # interior positive offset
+    box = sections.rectangular_hollow_section(d=200, b=100, t=10, r_out=0, n_r=1)
+    box = box.offset_perimeter(amount=5, where="interior")
+    box.create_mesh([50])
+    section = Section(box)
+    section.calculate_geometric_properties()
+    area = 200 * 100 - 190 * 90 + (10 * 10 - np.pi * 5 * 5)
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+    # all negative offset
+    box = sections.rectangular_hollow_section(d=200, b=100, t=10, r_out=0, n_r=1)
+    box = box.offset_perimeter(amount=-2.5, where="all")
+    box.create_mesh([50])
+    section = Section(box)
+    section.calculate_geometric_properties()
+    area = 195 * 95 - 185 * 85 + (5 * 5 - np.pi * 2.5 * 2.5)
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+    # all positive offset
+    box = sections.rectangular_hollow_section(d=200, b=100, t=10, r_out=0, n_r=1)
+    box = box.offset_perimeter(amount=5, where="all")
+    box.create_mesh([50])
+    section = Section(box)
+    section.calculate_geometric_properties()
+    area = 210 * 110 - (10 * 10 - np.pi * 5 * 5) - 170 * 70
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+
+def test_compound_rectangular_offset():
+    rect1 = sections.rectangular_section(d=50, b=50)
+    rect2 = sections.rectangular_section(d=50, b=50).align_to(rect1, "right")
+    geom = rect1 + rect2
+    geom = geom.offset_perimeter(amount=-5, where="exterior")
+    geom.create_mesh([50])
+    section = Section(geom)
+    section.calculate_geometric_properties()
+    area = 90 * 40
+    check.almost_equal(section.get_area(), area, rel=r_tol)
+
+
+def test_compound_rectangular_isection_offset_corrode():
+    d = 300
+    b = 150
+    tf = 10
+    tw = 8
+    r = 12
+    b_p = 250
+    t_p = 16
+    ub = sections.i_section(d=d, b=b, t_f=tf, t_w=tw, r=r, n_r=16)
+    plate = (
+        sections.rectangular_section(b=b_p, d=t_p)
+        .align_center(ub)
+        .align_to(ub, on="top")
+    )
+    geom_test = ub + plate
+    geom_test = geom_test.offset_perimeter(amount=-2, where="exterior")
+    geom_test.create_mesh([100])
+    section_test = Section(geom_test)
+    section_test.calculate_geometric_properties()
+
+    ub_corroded = sections.mono_i_section(
+        d=298, b_t=146, b_b=146, t_ft=8, t_fb=6, t_w=4, r=14, n_r=16
+    )
+    plate_corroded1 = (
+        sections.rectangular_section(b=146, d=2)
+        .align_center(ub_corroded)
+        .align_to(ub_corroded, "top")
+    )
+    plate_corroded2 = (
+        sections.rectangular_section(b=246, d=12)
+        .align_center(ub_corroded)
+        .align_to(plate_corroded1, "top")
+    )
+    rad_l = (
+        draw_radius(2, 8)
+        .align_to(plate_corroded1, "left")
+        .align_to(plate_corroded2, "bottom")
+    )
+    rad_r = (
+        draw_radius(2, 8)
+        .mirror_section("y", [2, 0])
+        .align_to(plate_corroded1, "right")
+        .align_to(plate_corroded2, "bottom")
+    )
+    geom_corroded = ub_corroded + plate_corroded1 + plate_corroded2 + rad_l + rad_r
+    geom_corroded.create_mesh([100])
+    section_corroded = Section(geom_corroded)
+    section_corroded.calculate_geometric_properties()
+
+    check.almost_equal(section_test.get_area(), section_corroded.get_area(), rel=r_tol)
+
+
+def test_compound_stiffened_isection():
+    """
+    Tests that plates 1 and 2 can be eroded to nothing and a valid Section can
+    still be generated without errors.
+    """
+    uc = sections.i_section(d=400, b=400, t_f=25, t_w=25, r=30, n_r=8)
+    plate1 = (
+        sections.rectangular_section(b=500, d=10).align_center(uc).align_to(uc, "top")
+    )
+    plate2 = (
+        sections.rectangular_section(b=500, d=10)
+        .align_center(uc)
+        .align_to(uc, "bottom")
+    )
+    geom = uc + plate1 + plate2
+
+    new_geom = geom.offset_perimeter(-9)
+    new_geom.create_mesh([100])
+    section = Section(new_geom)
+
+    new_geom = geom.offset_perimeter(-10)
+    new_geom.create_mesh([100])
+    section = Section(new_geom)
+
+    new_geom = geom.offset_perimeter(-11)
+    new_geom.create_mesh([100])
+    section = Section(new_geom)
+
+
+def draw_radius(r, n):
+    points = []
+
+    # calculate radius of points
+    for i in range(n):
+        # determine angle
+        t = i * 1.0 / max(1, n - 1) * np.pi * 0.5
+
+        x = r * np.cos(t)
+        y = r * np.sin(t)
+        points.append([x, y])
+
+    points.append([r, r])
+
+    return Geometry(Polygon(points))
diff --git a/sectionproperties/tests/test_perimeter.py b/sectionproperties/tests/test_perimeter.py
new file mode 100644
index 00000000..8426240d
--- /dev/null
+++ b/sectionproperties/tests/test_perimeter.py
@@ -0,0 +1,151 @@
+import pytest_check as check
+import numpy as np
+import sectionproperties.pre.sections as sections
+from sectionproperties.analysis.cross_section import Section
+
+
+r_tol = 1e-3
+
+
+def test_rectangular_perimeter():
+    rect = sections.rectangular_section(d=500, b=300)
+    rect.create_mesh([200])
+    section = Section(rect)
+    section.calculate_geometric_properties()
+    assert section.get_perimeter() == 2 * (500 + 300)
+
+
+def test_i_section():
+    i_section = sections.i_section(d=308, b=305, t_f=15.4, t_w=9.9, r=16.5, n_r=16)
+    i_section.create_mesh([100])
+    section = Section(i_section)
+    section.calculate_geometric_properties()
+    perim = (
+        (2 * 305)
+        + (4 * 15.4)
+        + 2 * (305 - 9.9 - 2 * 16.5)
+        + (2 * np.pi * 16.5)
+        + 2 * (308 - 2 * 15.4 - 2 * 16.5)
+    )
+    check.almost_equal(section.get_perimeter(), perim, rel=r_tol)
+
+
+def test_box_girder_perimeter():
+    box_girder = sections.box_girder_section(
+        d=400, b_t=700, b_b=100, t_ft=20, t_fb=20, t_w=12
+    )
+    box_girder.create_mesh([100])
+    section = Section(box_girder)
+    section.calculate_geometric_properties()
+    assert section.get_perimeter() == 700 + 100 + 2 * 500
+
+
+def test_custom_geometry_perimeter():
+    points = [[0, 0], [5, 0], [11, 8], [3, 2], [0, 2]]
+    facets = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 0]]
+    control_points = [[5, 5]]
+    custom = sections.Geometry.from_points(points, facets, control_points, holes=None)
+    custom.create_mesh([100])
+    section = Section(custom)
+    section.calculate_geometric_properties()
+    assert section.get_perimeter() == 5 + 2 + 10 + 10 + 3
+
+
+def test_compound_rectangular_perimeter():
+    rect1 = sections.rectangular_section(d=100, b=100)
+    rect2 = sections.rectangular_section(d=50, b=50).align_to(rect1, on="top")
+    rect3 = sections.rectangular_section(d=50, b=100).align_to(rect1, on="right")
+    rect4 = (
+        sections.rectangular_section(d=50, b=50)
+        .align_to(rect3, on="bottom")
+        .align_to(rect3, on="right", inner=True)
+    )
+    geom = rect1 + rect2 + rect3 + rect4
+    geom.create_mesh([100])
+    section = Section(geom)
+    section.calculate_geometric_properties()
+    assert section.get_perimeter() == (
+        150 + 50 + 50 + 100 + 100 + 50 + 50 + 50 + 50 + 150
+    )
+
+
+def test_compound_rectangular_isection_perimeter1():
+    d = 300
+    b = 150
+    tf = 10
+    tw = 6
+    r = 12
+    b_p = 250
+    t_p = 16
+    ub = sections.i_section(d=d, b=b, t_f=tf, t_w=tw, r=r, n_r=16)
+    plate = (
+        sections.rectangular_section(b=b_p, d=t_p)
+        .align_center(ub)
+        .align_to(ub, on="top")
+    )
+    geom = ub + plate
+    geom.create_mesh([100])
+    section = Section(geom)
+    section.calculate_geometric_properties()
+    perim = (
+        b
+        + (4 * tf)
+        + 2 * (b - tw - 2 * r)
+        + (2 * np.pi * r)
+        + 2 * (d - 2 * tf - 2 * r)
+        + (b_p - b)
+        + (2 * t_p)
+        + b_p
+    )
+    check.almost_equal(section.get_perimeter(), perim, rel=r_tol)
+
+
+def test_compound_rectangular_isection_perimeter2():
+    i_section = sections.i_section(d=308, b=305, t_f=15.4, t_w=9.9, r=16.5, n_r=16)
+    rect1 = (
+        sections.rectangular_section(d=330, b=16)
+        .align_center(i_section)
+        .align_to(i_section, "left")
+    )
+    rect2 = (
+        sections.rectangular_section(d=330, b=16)
+        .align_center(i_section)
+        .align_to(i_section, "right")
+    )
+    geom = i_section + rect1 + rect2
+    geom.create_mesh([100])
+    section = Section(geom)
+    section.calculate_geometric_properties()
+    assert section.get_perimeter() == 2 * 330 + 4 * 16 + 2 * 305 + 2 * (330 - 308)
+
+
+def test_compound_rhs_isection_perimeter():
+    d = 200
+    b = 150
+    t = 9
+    r = 15
+    b_p = 250
+    t_p = 16
+    rhs = sections.rectangular_hollow_section(d=d, b=b, t=t, r_out=r, n_r=16)
+    plate1 = (
+        sections.rectangular_section(b=b_p, d=t_p)
+        .align_center(rhs)
+        .align_to(rhs, on="top")
+    )
+    plate2 = (
+        sections.rectangular_section(b=b_p, d=t_p)
+        .align_center(rhs)
+        .align_to(rhs, on="bottom")
+    )
+    geom = rhs + plate1 + plate2
+    geom.create_mesh([100])
+    section = Section(geom)
+    section.calculate_geometric_properties()
+    perim = (
+        (2 * b_p)
+        + (4 * t_p)
+        + 2 * (b_p - b + 2 * r)
+        + (2 * np.pi * r)
+        + 2 * (d - 2 * r)
+    )
+    check.almost_equal(section.get_perimeter(), perim, rel=r_tol)
diff --git a/sectionproperties/tests/test_sections.py b/sectionproperties/tests/test_sections.py
index fe954680..4bc2c3c4 100644
--- a/sectionproperties/tests/test_sections.py
+++ b/sectionproperties/tests/test_sections.py
@@ -2,10 +2,17 @@
 import pytest
 
 from sectionproperties.pre.sections import *
+from sectionproperties.pre.nastran_sections import *
 from sectionproperties.analysis.cross_section import Section
 from sectionproperties.pre.pre import Material
 from sectionproperties.pre.rhino import load_3dm, load_brep_encoding
-from shapely.geometry import Polygon, MultiPolygon
+from shapely.geometry import (
+    Polygon,
+    MultiPolygon,
+    LineString,
+    Point,
+    GeometryCollection,
+)
 from shapely import wkt
 import json
 
@@ -77,11 +84,13 @@ def test_geometry_from_points():
     ]
     control_points = [[0, 0]]
     holes = [[0, 6], [0, -6]]
-    new_geom = Geometry.from_points(points=points, facets=facets, control_points=control_points, holes=holes)
-    assert (
-        new_geom.geom.wkt
-        == "POLYGON ((6 10, 6 -10, -6 -10, -6 10, 6 10), (-4 4, 4 4, 4 8, -4 8, -4 4), (4 -8, 4 -4, -4 -4, -4 -8, 4 -8))"
-    )  # Note, the order of point coordinates can change based on geometry operations, e.g. +/- and how many; seems to be a shapely issue
+    new_geom = Geometry.from_points(
+        points=points, facets=facets, control_points=control_points, holes=holes
+    )
+    wkt_test_geom = shapely.wkt.loads(
+        "POLYGON ((6 10, 6 -10, -6 -10, -6 10, 6 10), (-4 4, 4 4, 4 8, -4 8, -4 4), (4 -8, 4 -4, -4 -4, -4 -8, 4 -8))"
+    )
+    assert (new_geom.geom - wkt_test_geom) == Polygon()
 
 
 def test_compound_geometry_from_points():
@@ -120,10 +129,10 @@ def test_compound_geometry_from_points():
     ]
     control_points = [[0, 0], [0, -2 * a - t / 2]]
     new_geom = CompoundGeometry.from_points(points, facets, control_points)
-    assert (
-        new_geom.geom.wkt
-        == "MULTIPOLYGON (((-0.05 -2, 0.05 -2, 0.05 -0.05, 1 -0.05, 1 0.05, -0.05 0.05, -0.05 -2)), ((-1 -2, 1 -2, 1 -2.1, -1 -2.1, -1 -2)))"
+    wkt_test_geom = shapely.wkt.loads(
+        "MULTIPOLYGON (((-0.05 -2, 0.05 -2, 0.05 -0.05, 1 -0.05, 1 0.05, -0.05 0.05, -0.05 -2)), ((-1 -2, 1 -2, 1 -2.1, -1 -2.1, -1 -2)))"
     )
+    assert (new_geom.geom - wkt_test_geom) == Polygon()
 
 
 def test_nested_compound_geometry_from_points():
@@ -172,10 +181,11 @@ def test_nested_compound_geometry_from_points():
     nested_compound = CompoundGeometry.from_points(
         points=points, facets=facets, control_points=control_points, holes=holes
     )
-    assert (
-        nested_compound.geom.wkt
-        == "MULTIPOLYGON (((50 50, 50 -50, -50 -50, -50 50, 50 50), (12.5 12.5, -12.5 12.5, -12.5 -12.5, 12.5 -12.5, 12.5 12.5)), ((-37.5 -37.5, -37.5 37.5, 37.5 37.5, 37.5 -37.5, -37.5 -37.5), (12.5 12.5, -12.5 12.5, -12.5 -12.5, 12.5 -12.5, 12.5 12.5)), ((-25 -25, -25 25, 25 25, 25 -25, -25 -25), (12.5 12.5, -12.5 12.5, -12.5 -12.5, 12.5 -12.5, 12.5 12.5)))"
+    wkt_test_geom = shapely.wkt.loads(
+        "MULTIPOLYGON (((50 50, 50 -50, -50 -50, -50 50, 50 50), (12.5 12.5, -12.5 12.5, -12.5 -12.5, 12.5 -12.5, 12.5 12.5)), ((-37.5 -37.5, -37.5 37.5, 37.5 37.5, 37.5 -37.5, -37.5 -37.5), (12.5 12.5, -12.5 12.5, -12.5 -12.5, 12.5 -12.5, 12.5 12.5)), ((-25 -25, -25 25, 25 25, 25 -25, -25 -25), (12.5 12.5, -12.5 12.5, -12.5 -12.5, 12.5 -12.5, 12.5 12.5)))"
     )
+    assert (nested_compound.geom - wkt_test_geom) == Polygon()
+
     assert nested_compound.control_points == [
         (-43.75, 0.0),
         (-31.25, 0.0),
@@ -207,10 +217,22 @@ def test_plastic_centroid():
     # correct "center" of the original section which is affected by EA of each
     # of the constituent geometries.
 
-    steel = Material(name='Steel', elastic_modulus=200e3, poissons_ratio=0.3, density=7.85e-6,
-                    yield_strength=500, color='grey')
-    timber = Material(name='Timber', elastic_modulus=5e3, poissons_ratio=0.35, density=6.5e-7,
-                    yield_strength=20, color='burlywood')
+    steel = Material(
+        name="Steel",
+        elastic_modulus=200e3,
+        poissons_ratio=0.3,
+        density=7.85e-6,
+        yield_strength=500,
+        color="grey",
+    )
+    timber = Material(
+        name="Timber",
+        elastic_modulus=5e3,
+        poissons_ratio=0.35,
+        density=6.5e-7,
+        yield_strength=20,
+        color="burlywood",
+    )
 
     # create 310UB40.4
     ub = i_section(d=304, b=165, t_f=10.2, t_w=6.1, r=11.4, n_r=8, material=steel)
@@ -278,4 +300,48 @@ def test_geometry_from_3dm_encode():
         brep_encoded = json.load(file)
     exp = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
     test = Geometry.from_rhino_encoding(brep_encoded)
-    assert (test.geom-exp).is_empty
+    assert (test.geom - exp).is_empty
+
+
+def test_shift_points():
+    assymetrical_chan = nastran_chan(75, 200, 8, 16).shift_points(1, dy=-10)
+    assert (
+        assymetrical_chan.geom.wkt
+        == "POLYGON ((0 0, 75 -10, 75 16, 8 16, 8 184, 75 184, 75 200, 0 200, 0 0))"
+    )
+
+
+def test_mirror_section():
+    assymetrical_chan = nastran_chan(75, 200, 8, 16).shift_points(1, dy=-10)
+    assert (
+        assymetrical_chan.mirror_section(axis="x").geom.wkt
+        == "POLYGON ((0 190, 75 200, 75 174, 8 174, 8 6, 75 6, 75 -10, 0 -10, 0 190))"
+    )
+    assert (
+        assymetrical_chan.mirror_section(axis="y").geom.wkt
+        == "POLYGON ((75 0, 0 -10, 0 16, 67 16, 67 184, 0 184, 0 200, 75 200, 75 0))"
+    )
+    assert (
+        assymetrical_chan.mirror_section(axis="y", mirror_point=[50, 50]).geom.wkt
+        == "POLYGON ((100 0, 25 -10, 25 16, 92 16, 92 184, 25 184, 25 200, 100 200, 100 0))"
+    )
+    assert (
+        assymetrical_chan.mirror_section(axis="x", mirror_point=[50, 50]).geom.wkt
+        == "POLYGON ((0 100, 75 110, 75 84, 8 84, 8 -84, 75 -84, 75 -100, 0 -100, 0 100))"
+    )
+
+
+def test_filter_non_polygons():
+    point1 = Point([0, 0])
+    point2 = Point([1, 1])
+    point3 = Point([1, 0])
+    line = LineString([point1, point2])
+    poly = Polygon([point1, point2, point3])
+    multi_poly = MultiPolygon([poly, poly])
+    collection = GeometryCollection([poly, point1, line])
+    out = filter_non_polygons(collection)
+    assert filter_non_polygons(poly) == poly
+    assert filter_non_polygons(multi_poly) == multi_poly
+    assert filter_non_polygons(point1) == Polygon()
+    assert filter_non_polygons(line) == Polygon()
+    assert filter_non_polygons(collection) == poly