diff --git a/.github/release-drafter.yml b/.github/release-drafter.yml
index a779ce3e..0d8bd482 100644
--- a/.github/release-drafter.yml
+++ b/.github/release-drafter.yml
@@ -53,7 +53,10 @@ version-resolver:
 # Custom text at start of release
 header: >
 
-  Insert header here...
+  This release introduces the calculation of the yield moment when material properties
+  are applied to the section. The yield moment is calculated as part of the geometric
+  analysis with `calculate_geometric_properties()` and can be retrieved with the
+  `get_my()` or `get_my_p()` methods.
 
 template: |
 
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 3b01800e..5227b5bd 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -8,7 +8,7 @@ on:
     types: [opened, reopened, synchronize]
 
 env:
-  UV_VERSION: "0.4.27"
+  UV_VERSION: "0.4.29"
   DEFAULT_PYTHON_VERSION: "3.12"
 
 concurrency:
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 6838a01a..74d539a0 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -6,7 +6,7 @@ on:
       - master
 
 env:
-  UV_VERSION: "0.4.27"
+  UV_VERSION: "0.4.29"
   DEFAULT_PYTHON_VERSION: "3.12"
 
 jobs:
diff --git a/docs/user_guide.rst b/docs/user_guide.rst
index f2637ab5..0cba972a 100644
--- a/docs/user_guide.rst
+++ b/docs/user_guide.rst
@@ -52,7 +52,7 @@ Cross-Section Analysis
 
   * ☑ Second moments of area
   * ☑ Elastic section moduli
-  * ☐ Yield moment
+  * ☑ Yield moment
   * ☑ Radii of gyration
   * ☑ Plastic centroid
   * ☑ Plastic section moduli
@@ -62,7 +62,7 @@ Cross-Section Analysis
 
   * ☑ Second moments of area
   * ☑ Elastic section moduli
-  * ☐ Yield moment
+  * ☑ Yield moment
   * ☑ Radii of gyration
   * ☑ Plastic centroid
   * ☑ Plastic section moduli
diff --git a/docs/user_guide/results.rst b/docs/user_guide/results.rst
index 035e31f9..20cbbc1a 100644
--- a/docs/user_guide/results.rst
+++ b/docs/user_guide/results.rst
@@ -99,10 +99,12 @@ Geometric Analysis
     ~sectionproperties.analysis.section.Section.get_c
     ~sectionproperties.analysis.section.Section.get_eic
     ~sectionproperties.analysis.section.Section.get_ez
+    ~sectionproperties.analysis.section.Section.get_my
     ~sectionproperties.analysis.section.Section.get_rc
     ~sectionproperties.analysis.section.Section.get_eip
     ~sectionproperties.analysis.section.Section.get_phi
     ~sectionproperties.analysis.section.Section.get_ezp
+    ~sectionproperties.analysis.section.Section.get_my_p
     ~sectionproperties.analysis.section.Section.get_rp
     ~sectionproperties.analysis.section.Section.get_nu_eff
     ~sectionproperties.analysis.section.Section.get_e_eff
@@ -134,6 +136,8 @@ Plastic Analysis
     ~sectionproperties.analysis.section.Section.get_pc_p
     ~sectionproperties.analysis.section.Section.get_mp
     ~sectionproperties.analysis.section.Section.get_mp_p
+    ~sectionproperties.analysis.section.Section.get_sf
+    ~sectionproperties.analysis.section.Section.get_sf_p
 
 .. _label-material-affects-results:
 
diff --git a/docs/user_guide/theory.rst b/docs/user_guide/theory.rst
index bdcd084a..d78764e3 100644
--- a/docs/user_guide/theory.rst
+++ b/docs/user_guide/theory.rst
@@ -420,6 +420,20 @@ extreme (min. and max.) coordinates of the cross-section in the x and y-directio
   Z_{yy}^+ = \frac{I_{\overline{yy}}}{x_{max} - x_c} \\
   Z_{yy}^- = \frac{I_{\overline{yy}}}{x_c - x_{min}} \\
 
+Yield Moments
+~~~~~~~~~~~~~
+
+The yield moment is defined as the lowest bending moment that causes any point within
+cross-section to reach the yield strength. Note that this implementation is purely
+linear-elastic i.e. uses the linear-elastic modulus and bi-directional yield strength
+only.
+
+``sectionproperties`` applies a unit bending moment, about each axis separately, and
+determines the yield index for each point within the mesh. The yield index is defined as
+the stress divided by the material yield strength. Through this method, a critical yield
+index is determined (i.e. point which will yield first under bending) and the yield
+moment calculated as the inverse of the critical yield index.
+
 .. _label-theory-plastic-section-moduli:
 
 Plastic Section Moduli
diff --git a/pyproject.toml b/pyproject.toml
index 5aca03fa..6ce3fbb2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "sectionproperties"
-version = "3.5.0"
+version = "3.6.0"
 description = "A python package for the analysis of arbitrary cross-sections using the finite element method."
 readme = "README.md"
 license = {file = "LICENSE"}
@@ -74,14 +74,14 @@ docs = [
 ]
 test = [
     "pytest==8.3.3",
-    "pytest-benchmark[histogram]==4.0.0",
+    "pytest-benchmark[histogram]==5.1.0",
     "pytest-check==2.4.1",
     "coverage[toml]==7.6.4",
 ]
 lint = [
     "pre-commit==4.0.1",
     "ruff==0.7.1",
-    "pyright==1.1.386",
+    "pyright==1.1.387",
 ]
 
 [tool.uv]
diff --git a/src/sectionproperties/analysis/plastic_section.py b/src/sectionproperties/analysis/plastic_section.py
index 23767d1f..9a304730 100644
--- a/src/sectionproperties/analysis/plastic_section.py
+++ b/src/sectionproperties/analysis/plastic_section.py
@@ -204,46 +204,51 @@ def calculate_plastic_properties(
         if progress and task is not None:
             progress.update(task_id=task, advance=1)
 
-        # if there are no materials specified, calculate shape factors
+        # calculate shape factors
+        props = section.section_props
+
+        if (
+            props.sxx is None
+            or props.syy is None
+            or props.s11 is None
+            or props.s22 is None
+        ):
+            msg = "Property save failed."
+            raise RuntimeError(msg)
+
+        # for a geometric only analysis sf = s / z
         if list(set(section.materials)) == [pre.DEFAULT_MATERIAL]:
             if (
-                section.section_props.zxx_plus
-                and section.section_props.zxx_minus
-                and section.section_props.zyy_plus
-                and section.section_props.zyy_minus
+                props.zxx_plus
+                and props.zxx_minus
+                and props.zyy_plus
+                and props.zyy_minus
             ):
-                section.section_props.sf_xx_plus = (
-                    section.section_props.sxx / section.section_props.zxx_plus
-                )
-                section.section_props.sf_xx_minus = (
-                    section.section_props.sxx / section.section_props.zxx_minus
-                )
-                section.section_props.sf_yy_plus = (
-                    section.section_props.syy / section.section_props.zyy_plus
-                )
-                section.section_props.sf_yy_minus = (
-                    section.section_props.syy / section.section_props.zyy_minus
-                )
+                props.sf_xx_plus = props.sxx / props.zxx_plus
+                props.sf_xx_minus = props.sxx / props.zxx_minus
+                props.sf_yy_plus = props.syy / props.zyy_plus
+                props.sf_yy_minus = props.syy / props.zyy_minus
             if (
-                section.section_props.s11
-                and section.section_props.s22
-                and section.section_props.z11_plus
-                and section.section_props.z11_minus
-                and section.section_props.z22_plus
-                and section.section_props.z22_minus
+                props.z11_plus
+                and props.z11_minus
+                and props.z22_plus
+                and props.z22_minus
             ):
-                section.section_props.sf_11_plus = (
-                    section.section_props.s11 / section.section_props.z11_plus
-                )
-                section.section_props.sf_11_minus = (
-                    section.section_props.s11 / section.section_props.z11_minus
-                )
-                section.section_props.sf_22_plus = (
-                    section.section_props.s22 / section.section_props.z22_plus
-                )
-                section.section_props.sf_22_minus = (
-                    section.section_props.s22 / section.section_props.z22_minus
-                )
+                props.sf_11_plus = props.s11 / props.z11_plus
+                props.sf_11_minus = props.s11 / props.z11_minus
+                props.sf_22_plus = props.s22 / props.z22_plus
+                props.sf_22_minus = props.s22 / props.z22_minus
+        else:
+            if props.my_xx and props.my_yy:
+                props.sf_xx_plus = props.sxx / props.my_xx
+                props.sf_xx_minus = props.sf_xx_plus
+                props.sf_yy_plus = props.syy / props.my_yy
+                props.sf_yy_minus = props.sf_yy_plus
+            if props.my_11 and props.my_22:
+                props.sf_11_plus = props.s11 / props.my_11
+                props.sf_11_minus = props.sf_11_plus
+                props.sf_22_plus = props.s22 / props.my_22
+                props.sf_22_minus = props.sf_22_plus
 
         if progress and task is not None:
             progress.update(
diff --git a/src/sectionproperties/analysis/section.py b/src/sectionproperties/analysis/section.py
index 0e680333..ad2034f6 100644
--- a/src/sectionproperties/analysis/section.py
+++ b/src/sectionproperties/analysis/section.py
@@ -202,6 +202,7 @@ def calculate_geometric_properties(self) -> None:
           - Centroidal section moduli
           - Radii of gyration
           - Principal axis properties
+          - Yield moments (composite only)
         """
 
         def calculate_geom(progress: Progress | None = None) -> None:
@@ -259,11 +260,87 @@ def calculate_geom(progress: Progress | None = None) -> None:
             )
             self.section_props.e_eff = self.section_props.ea / self.section_props.area
             self.section_props.g_eff = self.section_props.ga / self.section_props.area
+
+            # calculate derived properties
             self.section_props.calculate_elastic_centroid()
             self.section_props.calculate_centroidal_properties(
                 node_list=self.mesh["vertices"]
             )
 
+            # calculate yield moments
+            self.section_props.my_xx = 0.0
+            self.section_props.my_yy = 0.0
+            self.section_props.my_11 = 0.0
+            self.section_props.my_22 = 0.0
+
+            # calculate yield moments:
+            # 1) loop through each material group and through each element in the group
+            # 2) for each point, calculate the bending stress from a unit bending moment
+            # in each direction (mxx, myy, m11, m22)
+            # 3) from this, calculate the yield index for each point
+            # 4) get the largest yield index and scale the bending moment such that the
+            # yield index is 1
+
+            # initialise the yield indexes
+            yield_index = {
+                "mxx": 0.0,
+                "myy": 0.0,
+                "m11": 0.0,
+                "m22": 0.0,
+            }
+
+            # get useful section properties
+            cx = self.section_props.cx
+            cy = self.section_props.cy
+            phi = self.section_props.phi
+            ixx = self.section_props.ixx_c
+            iyy = self.section_props.iyy_c
+            ixy = self.section_props.ixy_c
+            i11 = self.section_props.i11_c
+            i22 = self.section_props.i22_c
+
+            if ixx is None or iyy is None or ixy is None or i11 is None or i22 is None:
+                msg = "Section properties failed to save."
+                raise RuntimeError(msg)
+
+            # loop through each material group
+            for group in self.material_groups:
+                em = group.material.elastic_modulus
+                fy = group.material.yield_strength
+
+                # loop through each element in the material group
+                for el in group.elements:
+                    # loop through each node in the element
+                    for coord in el.coords.transpose():
+                        # calculate coordinates wrt centroidal & principal axes
+                        x = coord[0] - cx
+                        y = coord[1] - cy
+                        x11, y22 = fea.principal_coordinate(phi=phi, x=x, y=y)
+
+                        # calculate bending stresses due to unit moments
+                        sig_mxx = em * (
+                            -ixy / (ixx * iyy - ixy**2) * x
+                            + iyy / (ixx * iyy - ixy**2) * y
+                        )
+                        sig_myy = em * (
+                            -ixx / (ixx * iyy - ixy**2) * x
+                            + ixy / (ixx * iyy - ixy**2) * y
+                        )
+                        sig_m11 = em / i11 * y22
+                        sig_m22 = -em / i22 * x11
+
+                        # update yield indexes
+                        yield_index["mxx"] = max(yield_index["mxx"], abs(sig_mxx / fy))
+                        yield_index["myy"] = max(yield_index["myy"], abs(sig_myy / fy))
+                        yield_index["m11"] = max(yield_index["m11"], abs(sig_m11 / fy))
+                        yield_index["m22"] = max(yield_index["m22"], abs(sig_m22 / fy))
+
+            # calculate yield moments
+            self.section_props.my_xx = 1.0 / yield_index["mxx"]
+            self.section_props.my_yy = 1.0 / yield_index["myy"]
+            self.section_props.my_11 = 1.0 / yield_index["m11"]
+            self.section_props.my_22 = 1.0 / yield_index["m22"]
+
             if progress and task is not None:
                 msg = "[bold green]:white_check_mark: Geometric analysis complete"
                 progress.update(task_id=task, description=msg)
@@ -1107,7 +1184,8 @@ def calculate_plastic_properties(
         Note:
             If materials are specified, the values calculated for the plastic section
             moduli are displayed as plastic moments (i.e :math:`M_p = f_y S`) and the
-            shape factors are not calculated.
+            shape factors are calculated as the ratio between the plastic moment and the
+            yield moment.
 
         Warning:
             The geometric properties must be calculated prior to the calculation of the
@@ -1120,7 +1198,7 @@ def calculate_plastic_properties(
 
           - Plastic centroids (centroidal and principal axes)
           - Plastic section moduli (centroidal and principal axes)
-          - Shape factors, non-composite only (centroidal and principal axe)
+          - Shape factors (centroidal and principal axes)
         """
         # check that a geometric analysis has been performed
         if self.section_props.cx is None:
@@ -2240,6 +2318,32 @@ def get_ez(
             self.section_props.zyy_minus / e_ref,
         )
 
+    def get_my(self) -> tuple[float, float]:
+        """Returns the yield moment for bending about the centroidal axis.
+
+        This is a composite only property, as such this can only be returned if material
+        properties have been applied to the cross-section.
+
+        Returns:
+            Yield moment for bending about the centroidal ``x`` and ``y`` axes
+            (``my_xx``, ``my_yy``)
+
+        Raises:
+            RuntimeError: If material properties have *not* been applied
+            RuntimeError: If a geometric analysis has not been performed
+        """
+        if not self.is_composite():
+            msg = "Attempting to get a composite only property for a geometric analysis"
+            msg += " (material properties have not been applied). Consider using"
+            msg += " get_z()."
+            raise RuntimeError(msg)
+
+        if self.section_props.my_xx is None or self.section_props.my_yy is None:
+            msg = "Conduct a geometric analysis."
+            raise RuntimeError(msg)
+
+        return (self.section_props.my_xx, self.section_props.my_yy)
+
     def get_rc(self) -> tuple[float, float]:
         """Returns the cross-section centroidal radii of gyration.
 
@@ -2418,6 +2522,32 @@ def get_ezp(
             self.section_props.z22_minus / e_ref,
         )
 
+    def get_my_p(self) -> tuple[float, float]:
+        """Returns the yield moment for bending about the principal axis.
+
+        This is a composite only property, as such this can only be returned if material
+        properties have been applied to the cross-section.
+
+        Returns:
+            Yield moment for bending about the principal ``11`` and ``22`` axes
+            (``my_11``, ``my_22``)
+
+        Raises:
+            RuntimeError: If material properties have *not* been applied
+            RuntimeError: If a geometric analysis has not been performed
+        """
+        if not self.is_composite():
+            msg = "Attempting to get a composite only property for a geometric analysis"
+            msg += " (material properties have not been applied). Consider using"
+            msg += " get_zp()."
+            raise RuntimeError(msg)
+
+        if self.section_props.my_11 is None or self.section_props.my_22 is None:
+            msg = "Conduct a geometric analysis."
+            raise RuntimeError(msg)
+
+        return (self.section_props.my_11, self.section_props.my_22)
+
     def get_rp(self) -> tuple[float, float]:
         """Returns the cross-section principal radii of gyration.
 
@@ -3038,22 +3168,25 @@ def get_mp_p(self) -> tuple[float, float]:
     def get_sf(self) -> tuple[float, float, float, float]:
         """Returns the cross-section centroidal axis shape factors.
 
-        This is a geometric only property, as such this can only be returned if material
-        properties have *not* been applied to the cross-section.
+        For a geometric-only analysis, the shape factor is defined as the ratio between
+        the plastic section modulus and the elastic section modulus. For a composite
+        analysis the shape factors is defined as the ratio between the plastic moment
+        and the yield moment.
+
+        .. note::
+            For composite analyses, the ``plus`` values will always equal the ``minus``
+            values as the yield moment occurs when the first fibre reaches its yield
+            strength. For geometric-only analyses, the elastic section moduli is
+            calculated with respect to the top and bottom fibres (i.e. ``plus`` and
+            ``minus``).
 
         Returns:
             Centroidal axis shape factors with respect to the top and bottom fibres
             (``sf_xx_plus``, ``sf_xx_minus``, ``sf_yy_plus``, ``sf_yy_minus``)
 
         Raises:
-            RuntimeError: If material properties have been applied
             RuntimeError: If a plastic analysis has not been performed
         """
-        if self.is_composite():
-            msg = "Attempting to get a geometric only property for a composite analysis"
-            msg += " (material properties have been applied)."
-            raise RuntimeError(msg)
-
         if (
             self.section_props.sf_xx_plus is None
             or self.section_props.sf_xx_minus is None
@@ -3073,22 +3206,25 @@ def get_sf(self) -> tuple[float, float, float, float]:
     def get_sf_p(self) -> tuple[float, float, float, float]:
         """Returns the cross-section principal axis shape factors.
 
-        This is a geometric only property, as such this can only be returned if material
-        properties have *not* been applied to the cross-section.
+        For a geometric-only analysis, the shape factor is defined as the ratio between
+        the plastic section modulus and the elastic section modulus. For a composite
+        analysis the shape factors is defined as the ratio between the plastic moment
+        and the yield moment.
+
+        .. note::
+            For composite analyses, the ``plus`` values will always equal the ``minus``
+            values as the yield moment occurs when the first fibre reaches its yield
+            strength. For geometric-only analyses, the elastic section moduli is
+            calculated with respect to the top and bottom fibres (i.e. ``plus`` and
+            ``minus``).
 
         Returns:
             Principal bending axis shape factors with respect to the top and bottom
             fibres (``sf_11_plus``, ``sf_11_minus``, ``sf_22_plus``, ``sf_22_minus``)
 
         Raises:
-            RuntimeError: If material properties have been applied
             RuntimeError: If a plastic analysis has not been performed
         """
-        if self.is_composite():
-            msg = "Attempting to get a geometric only property for a composite analysis"
-            msg += " (material properties have been applied)."
-            raise RuntimeError(msg)
-
         if (
             self.section_props.sf_11_plus is None
             or self.section_props.sf_11_minus is None
diff --git a/src/sectionproperties/post/post.py b/src/sectionproperties/post/post.py
index f15896b5..5dd6ecad 100644
--- a/src/sectionproperties/post/post.py
+++ b/src/sectionproperties/post/post.py
@@ -72,6 +72,10 @@ class SectionProperties:
             negative extreme value of the 11-axis
         r11_c: Radius of gyration about the principal 11-axis.
         r22_c: Radius of gyration about the principal 22-axis.
+        my_xx: Yield moment about the x-axis
+        my_yy: Yield moment about the y-axis
+        my_11: Yield moment about the 11-axis
+        my_22: Yield moment about the 22-axis
         j: Torsion constant
         omega: Warping function
         psi_shear: Psi shear function
@@ -165,6 +169,10 @@ class SectionProperties:
     r11_c: float | None = None
     r22_c: float | None = None
     j: float | None = None
+    my_xx: float | None = None
+    my_yy: float | None = None
+    my_11: float | None = None
+    my_22: float | None = None
     omega: npt.NDArray[np.float64] | None = None
     psi_shear: npt.NDArray[np.float64] | None = None
     phi_shear: npt.NDArray[np.float64] | None = None
@@ -278,7 +286,7 @@ def calculate_centroidal_properties(
             else:
                 self.phi = np.arctan2(self.ixx_c - self.i11_c, self.ixy_c) * 180 / np.pi
 
-            # initialise min, max variables TODO: check for `if xxx:` where xxx is float
+            # initialise min, max variables
             if self.phi is not None:
                 x1, y2 = fea.principal_coordinate(
                     phi=self.phi,
@@ -662,6 +670,15 @@ def print_results(
         except RuntimeError:
             pass
 
+    # print cross-section my
+    if is_composite:
+        try:
+            my_xx, my_yy = section.get_my()
+            table.add_row("my_xx", f"{my_xx:>{fmt}}")
+            table.add_row("my_yy", f"{my_yy:>{fmt}}")
+        except RuntimeError:
+            pass
+
     # print cross-section rc
     try:
         rx, ry = section.get_rc()
@@ -713,6 +730,15 @@ def print_results(
         except RuntimeError:
             pass
 
+    # print cross-section my_p
+    if is_composite:
+        try:
+            my_11, my_22 = section.get_my_p()
+            table.add_row("my_11", f"{my_11:>{fmt}}")
+            table.add_row("my_22", f"{my_22:>{fmt}}")
+        except RuntimeError:
+            pass
+
     # print cross-section rp
     try:
         r11, r22 = section.get_rp()
@@ -887,26 +913,24 @@ def print_results(
             pass
 
     # print cross-section sf
-    if not is_composite:
-        try:
-            sf_xx_plus, sf_xx_minus, sf_yy_plus, sf_yy_minus = section.get_sf()
-            table.add_row("sf_xx+", f"{sf_xx_plus:>{fmt}}")
-            table.add_row("sf_xx-", f"{sf_xx_minus:>{fmt}}")
-            table.add_row("sf_yy+", f"{sf_yy_plus:>{fmt}}")
-            table.add_row("sf_yy-", f"{sf_yy_minus:>{fmt}}")
-        except RuntimeError:
-            pass
+    try:
+        sf_xx_plus, sf_xx_minus, sf_yy_plus, sf_yy_minus = section.get_sf()
+        table.add_row("sf_xx+", f"{sf_xx_plus:>{fmt}}")
+        table.add_row("sf_xx-", f"{sf_xx_minus:>{fmt}}")
+        table.add_row("sf_yy+", f"{sf_yy_plus:>{fmt}}")
+        table.add_row("sf_yy-", f"{sf_yy_minus:>{fmt}}")
+    except RuntimeError:
+        pass
 
     # print cross-section sf_p
-    if not is_composite:
-        try:
-            sf_11_plus, sf_11_minus, sf_22_plus, sf_22_minus = section.get_sf_p()
-            table.add_row("sf_11+", f"{sf_11_plus:>{fmt}}")
-            table.add_row("sf_11-", f"{sf_11_minus:>{fmt}}")
-            table.add_row("sf_22+", f"{sf_22_plus:>{fmt}}")
-            table.add_row("sf_22-", f"{sf_22_minus:>{fmt}}")
-        except RuntimeError:
-            pass
+    try:
+        sf_11_plus, sf_11_minus, sf_22_plus, sf_22_minus = section.get_sf_p()
+        table.add_row("sf_11+", f"{sf_11_plus:>{fmt}}")
+        table.add_row("sf_11-", f"{sf_11_minus:>{fmt}}")
+        table.add_row("sf_22+", f"{sf_22_plus:>{fmt}}")
+        table.add_row("sf_22-", f"{sf_22_minus:>{fmt}}")
+    except RuntimeError:
+        pass
 
     console = Console()
     console.print(table)
diff --git a/tests/analysis/test_yield_moment.py b/tests/analysis/test_yield_moment.py
new file mode 100644
index 00000000..86479f72
--- /dev/null
+++ b/tests/analysis/test_yield_moment.py
@@ -0,0 +1,279 @@
+"""Tests for the yield moment calculation."""
+
+import numpy as np
+import pytest
+
+from sectionproperties.analysis import Section
+from sectionproperties.post.stress_post import StressPost
+from sectionproperties.pre import Material
+from sectionproperties.pre.library import i_section, rectangular_section
+
+STEEL = Material(
+    name="Steel",
+    elastic_modulus=200e3,
+    poissons_ratio=0.3,
+    yield_strength=500,
+    density=7.85e-6,
+    color="grey",
+)
+
+
+def test_get_without_analysis():
+    """Test for raising an error if a geometric analysis has not been performed."""
+    geom = rectangular_section(d=50, b=50, material=STEEL)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+
+    with pytest.raises(RuntimeError, match="Conduct a geometric analysis."):
+        sec.get_my()
+
+    with pytest.raises(RuntimeError, match="Conduct a geometric analysis."):
+        sec.get_my_p()
+
+
+def test_non_composite():
+    """Test for raising an error for non-composite analyses."""
+    geom = rectangular_section(d=50, b=50)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    with pytest.raises(RuntimeError, match="Attempting to get a composite only"):
+        sec.get_my()
+
+    with pytest.raises(RuntimeError, match="Attempting to get a composite only"):
+        sec.get_my_p()
+
+
+def test_rectangle():
+    """Test the yield moment of a simple rectangle."""
+    geom = rectangular_section(d=100, b=50, material=STEEL)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    my_xx_calc = 50 * 100 * 100 / 6 * 500
+    my_yy_calc = 100 * 50 * 50 / 6 * 500
+
+    # compare with theoretical values
+    assert my_xx == pytest.approx(my_xx_calc)
+    assert my_yy == pytest.approx(my_yy_calc)
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+    # check shape factors
+    sec.calculate_plastic_properties()
+    sf_xx, _, sf_yy, _ = sec.get_sf()
+    sf_11, _, sf_22, _ = sec.get_sf_p()
+    assert sf_xx == pytest.approx(1.5)
+    assert sf_yy == pytest.approx(1.5)
+    assert sf_11 == pytest.approx(1.5)
+    assert sf_22 == pytest.approx(1.5)
+
+
+def test_rectangle_rotated():
+    """Test the yield moment of a simple rotated rectangle."""
+    geom = rectangular_section(d=100, b=50, material=STEEL)
+    geom.rotate_section(angle=45.0)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    my_11, my_22 = sec.get_my()
+    my_xx_calc = 50 * 100 * 100 / 6 * 500
+    my_yy_calc = 100 * 50 * 50 / 6 * 500
+
+    # compare with theoretical values
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+    # check shape factors
+    sec.calculate_plastic_properties()
+    sf_11, _, sf_22, _ = sec.get_sf_p()
+    assert sf_11 == pytest.approx(1.5)
+    assert sf_22 == pytest.approx(1.5)
+
+
+def test_isection():
+    """Test the yield moment of an isection."""
+    geom = i_section(d=200, b=100, t_f=10, t_w=5, r=12, n_r=8, material=STEEL)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    ezxx, _, ezyy, _ = sec.get_ez()
+    my_xx_calc = ezxx / 200e3 * 500
+    my_yy_calc = ezyy / 200e3 * 500
+
+    # compare with theoretical values
+    assert my_xx == pytest.approx(my_xx_calc)
+    assert my_yy == pytest.approx(my_yy_calc)
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+    # check that shape factors with a geometric-only analysis match the composite
+    sec.calculate_plastic_properties()
+    sf_xx_c, _, sf_yy_c, _ = sec.get_sf()
+
+    geom = i_section(d=200, b=100, t_f=10, t_w=5, r=12, n_r=8)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+    sec.calculate_plastic_properties()
+
+    sf_xx_plus, sf_xx_minus, sf_yy_plus, sf_yy_minus = sec.get_sf()
+    assert sf_xx_c == pytest.approx(sf_xx_plus)
+    assert sf_xx_c == pytest.approx(sf_xx_minus)
+    assert sf_yy_c == pytest.approx(sf_yy_plus)
+    assert sf_yy_c == pytest.approx(sf_yy_minus)
+
+
+def test_rectangle_composite():
+    """Test the yield moment of a composite rectangular section."""
+    mat1 = Material("a", 2, 0, 2, 1, "b")
+    mat2 = Material("b", 1, 0, 1, 1, "r")
+
+    rect1 = rectangular_section(d=20, b=20, material=mat1)
+    rect2 = rectangular_section(d=20, b=20, material=mat2).align_to(rect1, "top")
+    rect3 = rectangular_section(d=20, b=20, material=mat1).align_to(rect2, "top")
+    geom = rect1 + rect2 + rect3
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geom)
+    sec.calculate_geometric_properties()
+
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    eixx_calc = (20 * 20**3 / 12) + 2 * 2 * ((20 * 20**3 / 12) + (20 * 20 * 20**2))
+    eiyy_calc = 5 * (20 * 20**3 / 12)
+    # myield = f * Ieff / y
+    my_xx_calc = 2 * (eixx_calc / 2) / 30
+    my_yy_calc = 1 * eiyy_calc / 10
+
+    # compare with theoretical values
+    assert my_xx == pytest.approx(my_xx_calc)
+    assert my_yy == pytest.approx(my_yy_calc)
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_composite_example():
+    """Tests the composite example from the docs, compares to stress analysis."""
+    timber = Material(
+        name="Timber",
+        elastic_modulus=8e3,
+        poissons_ratio=0.35,
+        yield_strength=20,
+        density=0.78e-6,
+        color="burlywood",
+    )
+
+    ub = i_section(d=304, b=165, t_f=10.2, t_w=6.1, r=11.4, n_r=8, material=STEEL)
+    panel = rectangular_section(d=100, b=600, material=timber)
+    panel = panel.align_center(align_to=ub).align_to(other=ub, on="top")
+    geom = ub + panel
+    geom.create_mesh(mesh_sizes=[10, 500])
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    # compare with stress analysis
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_yield_internal():
+    """Test where the internal material yields first."""
+    mat1 = Material(
+        name="Material_1",
+        elastic_modulus=0.75,
+        poissons_ratio=0,
+        yield_strength=1,
+        density=1,
+        color="gold",
+    )
+    mat2 = Material(
+        name="Material_2",
+        elastic_modulus=1,
+        poissons_ratio=0,
+        yield_strength=1,
+        density=1,
+        color="blue",
+    )
+    mat3 = Material(
+        name="Material 3",
+        elastic_modulus=3,
+        poissons_ratio=0,
+        yield_strength=1,
+        density=1,
+        color="red",
+    )
+
+    sq1 = rectangular_section(d=100, b=100, material=mat1).align_center()
+    sq2 = rectangular_section(d=75, b=75, material=mat2).align_center()
+    sq3 = rectangular_section(d=50, b=50, material=mat3).align_center()
+    hole = rectangular_section(d=25, b=25).align_center()
+    compound = (sq1 - sq2) + (sq2 - sq3) + (sq3 - hole)
+    compound.create_mesh(10)
+    sec = Section(compound)
+    sec.calculate_geometric_properties()
+
+    # compare with stress analysis
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def check_yield_index(stress: StressPost, key: str) -> float:
+    """Returns the largest yield index.
+
+    Given the output from StressPost object and a dict key representing the type of
+    stress, returns the largest yield index.
+    """
+    yield_index = 0.0
+
+    for idx, mat_dict in enumerate(stress.get_stress()):
+        fy = stress.material_groups[idx].material.yield_strength
+        sigs = mat_dict[key]
+        sig_max = max(np.absolute(sigs))
+
+        yield_index = max(yield_index, sig_max / fy)
+
+    return yield_index
diff --git a/tests/post/test_get_results.py b/tests/post/test_get_results.py
index d397e4a0..ca98c60e 100644
--- a/tests/post/test_get_results.py
+++ b/tests/post/test_get_results.py
@@ -798,9 +798,6 @@ def test_get_plastic(sec_no_mat_and_mat):
     assert sf_yy_plus == pytest.approx(1.5)
     assert sf_yy_minus == pytest.approx(1.5)
 
-    with pytest.raises(RuntimeError):
-        rect_mat.get_sf()
-
     # check sf_p
     sf_11_plus, sf_11_minus, sf_22_plus, sf_22_minus = rect_no_mat.get_sf_p()
     assert sf_11_plus == pytest.approx(1.5)
@@ -808,9 +805,6 @@ def test_get_plastic(sec_no_mat_and_mat):
     assert sf_22_plus == pytest.approx(1.5)
     assert sf_22_minus == pytest.approx(1.5)
 
-    with pytest.raises(RuntimeError):
-        rect_mat.get_sf_p()
-
 
 def test_get_effective_material():
     """Tests effective material properties."""
diff --git a/uv.lock b/uv.lock
index a1a0cf41..07587798 100644
--- a/uv.lock
+++ b/uv.lock
@@ -1707,15 +1707,15 @@ wheels = [
 
 [[package]]
 name = "pyright"
-version = "1.1.386"
+version = "1.1.387"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "nodeenv" },
     { name = "typing-extensions" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/92/50/1a57054b5585fa72a93a6244c1b4b5639f8f7a1cc60b2e807cc67da8f0bc/pyright-1.1.386.tar.gz", hash = "sha256:8e9975e34948ba5f8e07792a9c9d2bdceb2c6c0b61742b068d2229ca2bc4a9d9", size = 21949 }
+sdist = { url = "https://files.pythonhosted.org/packages/c2/32/e7187478d3105d6d7edc9b754d56472ee06557c25cc404911288fee1796a/pyright-1.1.387.tar.gz", hash = "sha256:577de60224f7fe36505d5b181231e3a395d427b7873be0bbcaa962a29ea93a60", size = 21939 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/cc/68/47fd6b3ffa27c99d7e0c866c618f07784b8806712059049daa492ca7e526/pyright-1.1.386-py3-none-any.whl", hash = "sha256:7071ac495593b2258ccdbbf495f1a5c0e5f27951f6b429bed4e8b296eb5cd21d", size = 18577 },
+    { url = "https://files.pythonhosted.org/packages/a0/18/c497df36641b0572f5bd59ae147b08ccaa6b8086397d50e1af97cc2ddcf6/pyright-1.1.387-py3-none-any.whl", hash = "sha256:6a1f495a261a72e12ad17e20d1ae3df4511223c773b19407cfa006229b1b08a5", size = 18577 },
 ]
 
 [[package]]
@@ -1737,21 +1737,22 @@ wheels = [
 
 [[package]]
 name = "pytest-benchmark"
-version = "4.0.0"
+version = "5.1.0"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "py-cpuinfo" },
     { name = "pytest" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/28/08/e6b0067efa9a1f2a1eb3043ecd8a0c48bfeb60d3255006dcc829d72d5da2/pytest-benchmark-4.0.0.tar.gz", hash = "sha256:fb0785b83efe599a6a956361c0691ae1dbb5318018561af10f3e915caa0048d1", size = 334641 }
+sdist = { url = "https://files.pythonhosted.org/packages/39/d0/a8bd08d641b393db3be3819b03e2d9bb8760ca8479080a26a5f6e540e99c/pytest-benchmark-5.1.0.tar.gz", hash = "sha256:9ea661cdc292e8231f7cd4c10b0319e56a2118e2c09d9f50e1b3d150d2aca105", size = 337810 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/4d/a1/3b70862b5b3f830f0422844f25a823d0470739d994466be9dbbbb414d85a/pytest_benchmark-4.0.0-py3-none-any.whl", hash = "sha256:fdb7db64e31c8b277dff9850d2a2556d8b60bcb0ea6524e36e28ffd7c87f71d6", size = 43951 },
+    { url = "https://files.pythonhosted.org/packages/9e/d6/b41653199ea09d5969d4e385df9bbfd9a100f28ca7e824ce7c0a016e3053/pytest_benchmark-5.1.0-py3-none-any.whl", hash = "sha256:922de2dfa3033c227c96da942d1878191afa135a29485fb942e85dff1c592c89", size = 44259 },
 ]
 
 [package.optional-dependencies]
 histogram = [
     { name = "pygal" },
     { name = "pygaljs" },
+    { name = "setuptools" },
 ]
 
 [[package]]
@@ -2151,7 +2152,7 @@ rhino = [
     { name = "rhino3dm" },
 ]
 
-[package.dependency-groups]
+[package.dev-dependencies]
 docs = [
     { name = "furo" },
     { name = "ipykernel" },
@@ -2191,7 +2192,7 @@ requires-dist = [
     { name = "shapely", specifier = "~=2.0.6" },
 ]
 
-[package.metadata.dependency-groups]
+[package.metadata.requires-dev]
 docs = [
     { name = "furo", specifier = "==2024.8.6" },
     { name = "ipykernel", specifier = "==6.29.5" },
@@ -2205,13 +2206,13 @@ docs = [
 ]
 lint = [
     { name = "pre-commit", specifier = "==4.0.1" },
-    { name = "pyright", specifier = "==1.1.386" },
+    { name = "pyright", specifier = "==1.1.387" },
     { name = "ruff", specifier = "==0.7.1" },
 ]
 test = [
     { name = "coverage", extras = ["toml"], specifier = "==7.6.4" },
     { name = "pytest", specifier = "==8.3.3" },
-    { name = "pytest-benchmark", extras = ["histogram"], specifier = "==4.0.0" },
+    { name = "pytest-benchmark", extras = ["histogram"], specifier = "==5.1.0" },
     { name = "pytest-check", specifier = "==2.4.1" },
 ]