diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 00000000..7a1ad43d
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,6 @@
+repos:
+  - repo: https://github.com/asottile/pyupgrade
+    rev: v3.21.2
+    hooks:
+      - id: pyupgrade
+        args: [--py310-plus]
diff --git a/Examples/Beam on Elastic Foundation.py b/Examples/Beam on Elastic Foundation.py
index a8db9562..966280cb 100644
--- a/Examples/Beam on Elastic Foundation.py	
+++ b/Examples/Beam on Elastic Foundation.py	
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 This example shows how to generate a beam on elastic foundation by using spring supports. All units
 in this model are expressed in terms of kips (force) and inches (length).
diff --git a/Pynite/Analysis.py b/Pynite/Analysis.py
index e96b761a..b279b653 100644
--- a/Pynite/Analysis.py
+++ b/Pynite/Analysis.py
@@ -75,7 +75,7 @@ def _prepare_model(model: FEModel3D, n_modes: int = 0) -> None:
     _renumber(model)
 
 
-def _identify_combos(model: FEModel3D, combo_tags: List[str] | None = None) -> List[LoadCombo]:
+def _identify_combos(model: FEModel3D, combo_tags: list[str] | None = None) -> list[LoadCombo]:
     """Returns a list of load combinations that are to be run based on tags given by the user.
 
     :param model: The model being analyzed.
@@ -164,7 +164,7 @@ def _check_stability(model: FEModel3D, K: NDArray[float64]) -> None:
     return
 
 
-def _PDelta(model: FEModel3D, combo_name: str, P1: NDArray[float64], FER1: NDArray[float64], D1_indices: List[int], D2_indices: List[int], D2: NDArray[float64], log: bool = True, sparse: bool = True, check_stability: bool = False, max_iter: int = 30) -> None:
+def _PDelta(model: FEModel3D, combo_name: str, P1: NDArray[float64], FER1: NDArray[float64], D1_indices: list[int], D2_indices: list[int], D2: NDArray[float64], log: bool = True, sparse: bool = True, check_stability: bool = False, max_iter: int = 30) -> None:
     """Performs second order (P-Delta) analysis. This type of analysis is appropriate for most models using beams, columns and braces. Second order analysis is usually required by material-specific codes. Models with slender members and/or members with combined bending and axial loads will generally have more significant P-Delta effects. P-Delta effects in plates/quads are not considered by Pynite at this time.
 
     :param model: The finite element model to be solved.
@@ -309,7 +309,7 @@ def _PDelta(model: FEModel3D, combo_name: str, P1: NDArray[float64], FER1: NDArr
     model.solution = 'P-Delta'
 
 
-def _pushover_step(model: FEModel3D, combo_name: str, push_combo: str, step_num: int, P1: NDArray[float64], FER1: NDArray[float64], D1_indices: List[int], D2_indices: List[int], D2: NDArray[float64], log: bool = True, sparse: bool = True, check_stability: bool = False) -> None:
+def _pushover_step(model: FEModel3D, combo_name: str, push_combo: str, step_num: int, P1: NDArray[float64], FER1: NDArray[float64], D1_indices: list[int], D2_indices: list[int], D2: NDArray[float64], log: bool = True, sparse: bool = True, check_stability: bool = False) -> None:
 
     # Run at least one iteration
     run_step = True
@@ -434,7 +434,7 @@ def _pushover_step(model: FEModel3D, combo_name: str, push_combo: str, step_num:
     _sum_displacements(model, Delta_D1, D2, D1_indices, D2_indices, model.load_combos[combo_name])
 
 
-def _unpartition_disp(model: FEModel3D, D1: NDArray[float64], D2: NDArray[float64], D1_indices: List[int], D2_indices: List[int]) -> NDArray[float64]:
+def _unpartition_disp(model: FEModel3D, D1: NDArray[float64], D2: NDArray[float64], D1_indices: list[int], D2_indices: list[int]) -> NDArray[float64]:
     """Unpartitions displacements from the solver and returns them as a global displacement vector
 
     :param model: The finite element model being evaluated
@@ -471,7 +471,7 @@ def _unpartition_disp(model: FEModel3D, D1: NDArray[float64], D2: NDArray[float6
     return D
 
 
-def _store_displacements(model: FEModel3D, D1: NDArray[float64], D2: NDArray[float64], D1_indices: List[int], D2_indices: List[int], combo: LoadCombo) -> None:
+def _store_displacements(model: FEModel3D, D1: NDArray[float64], D2: NDArray[float64], D1_indices: list[int], D2_indices: list[int], combo: LoadCombo) -> None:
     """Stores calculated displacements from the solver into the model's displacement vector `_D` and into each node object in the model
 
     :param model: The finite element model being evaluated.
@@ -507,7 +507,7 @@ def _store_displacements(model: FEModel3D, D1: NDArray[float64], D2: NDArray[flo
             node.RZ[combo.name] = D[node.ID*6 + 5, 0]
 
 
-def _sum_displacements(model: FEModel3D, Delta_D1: NDArray[float64], Delta_D2: NDArray[float64], D1_indices: List[int], D2_indices: List[int], combo: LoadCombo) -> None:
+def _sum_displacements(model: FEModel3D, Delta_D1: NDArray[float64], Delta_D2: NDArray[float64], D1_indices: list[int], D2_indices: list[int], combo: LoadCombo) -> None:
     """Sums calculated displacements for a load step from the solver into the model's displacement vector `_D` and into each node object in the model.
 
     :param model: The finite element model being evaluated.
@@ -682,7 +682,7 @@ def _check_TC_convergence(model: FEModel3D, combo_name: str = "Combo 1", log: bo
     return convergence
 
 
-def _calc_reactions(model: FEModel3D, log: bool = False, combo_tags: List[str] | None = None) -> None:
+def _calc_reactions(model: FEModel3D, log: bool = False, combo_tags: list[str] | None = None) -> None:
     """
     Calculates reactions internally once the model is solved.
 
@@ -939,7 +939,7 @@ def _calc_reactions(model: FEModel3D, log: bool = False, combo_tags: List[str] |
                 node.RxnMZ[combo.name] -= k*RZ
 
 
-def _check_statics(model: FEModel3D, combo_tags: List[str] | None = None) -> None:
+def _check_statics(model: FEModel3D, combo_tags: list[str] | None = None) -> None:
     '''
     Checks static equilibrium and prints results to the console.
 
@@ -1023,19 +1023,19 @@ def _check_statics(model: FEModel3D, combo_tags: List[str] | None = None) -> Non
             SumRMZ += RMZ - RFX*Y + RFY*X 
 
         # Add the results to the table
-        statics_table.add_row([combo.name, '{:.3g}'.format(SumFX), '{:.3g}'.format(SumRFX),
-                                            '{:.3g}'.format(SumFY), '{:.3g}'.format(SumRFY),
-                                            '{:.3g}'.format(SumFZ), '{:.3g}'.format(SumRFZ),
-                                            '{:.3g}'.format(SumMX), '{:.3g}'.format(SumRMX),
-                                            '{:.3g}'.format(SumMY), '{:.3g}'.format(SumRMY),
-                                            '{:.3g}'.format(SumMZ), '{:.3g}'.format(SumRMZ)])
+        statics_table.add_row([combo.name, f'{SumFX:.3g}', f'{SumRFX:.3g}',
+                                            f'{SumFY:.3g}', f'{SumRFY:.3g}',
+                                            f'{SumFZ:.3g}', f'{SumRFZ:.3g}',
+                                            f'{SumMX:.3g}', f'{SumRMX:.3g}',
+                                            f'{SumMY:.3g}', f'{SumRMY:.3g}',
+                                            f'{SumMZ:.3g}', f'{SumRMZ:.3g}'])
 
     # Print the static check table
     print(statics_table)
     print('')
 
 
-def _partition_D(model: FEModel3D) -> Tuple[List[int], List[int], NDArray[float64]]:
+def _partition_D(model: FEModel3D) -> tuple[list[int], list[int], NDArray[float64]]:
     """Builds a list with known nodal displacements and with the positions in global stiffness matrix of known and unknown nodal displacements
 
     :return: A list of the global matrix indices for the unknown nodal displacements (D1_indices). A list of the global matrix indices for the known nodal displacements (D2_indices). A list of the known nodal displacements (D2).
@@ -1131,7 +1131,7 @@ def _partition_D(model: FEModel3D) -> Tuple[List[int], List[int], NDArray[float6
     return D1_indices, D2_indices, D2
 
 
-def _partition(model: FEModel3D, unp_matrix: NDArray[float64] | lil_matrix, D1_indices: List[int], D2_indices: List[int]) -> Tuple[NDArray[float64], NDArray[float64]] | Tuple[NDArray[float64], NDArray[float64], NDArray[float64], NDArray[float64]]:
+def _partition(model: FEModel3D, unp_matrix: NDArray[float64] | lil_matrix, D1_indices: list[int], D2_indices: list[int]) -> tuple[NDArray[float64], NDArray[float64]] | tuple[NDArray[float64], NDArray[float64], NDArray[float64], NDArray[float64]]:
     """Partitions a matrix (or vector) into submatrices (or subvectors) based on degree of freedom boundary conditions.
 
     :param unp_matrix: The unpartitioned matrix (or vector) to be partitioned.
diff --git a/Pynite/BeamSegZ.py b/Pynite/BeamSegZ.py
index 9b40c687..840aba97 100644
--- a/Pynite/BeamSegZ.py
+++ b/Pynite/BeamSegZ.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 Created on Mon Nov  6 20:52:31 2017
 
@@ -138,7 +137,7 @@ def axial(self, x: float) -> float:
 
         return P1 + (p2 - p1)/(2*L)*x**2 + p1*x
 
-    def Torsion(self, x: float | List[float] = 0) -> float | None | NDArray[Any]:
+    def Torsion(self, x: float | list[float] = 0) -> float | None | NDArray[Any]:
         """
         Returns the torsional moment in the segment.
         """
diff --git a/Pynite/FEModel3D.py b/Pynite/FEModel3D.py
index 5f1fed05..98e22e1c 100644
--- a/Pynite/FEModel3D.py
+++ b/Pynite/FEModel3D.py
@@ -40,28 +40,28 @@ def __init__(self) -> None:
         # the data types they store, and then those types will be removed. This will give us the
         # ability to get type-based hints when using the dictionaries.
 
-        self.nodes: Dict[str, Node3D] = {}             # A dictionary of the model's nodes
-        self.materials: Dict[str, Material] = {}       # A dictionary of the model's materials
-        self.sections: Dict[str, Section] = {}         # A dictonary of the model's cross-sections
-        self.springs: Dict[str, Spring3D] = {}         # A dictionary of the model's springs
-        self.members: Dict[str, PhysMember] = {}       # A dictionary of the model's physical members
-        self.quads: Dict[str, Quad3D] = {}             # A dictionary of the model's quadiralterals
-        self.plates: Dict[str, Plate3D] = {}           # A dictionary of the model's rectangular plates
-        self.meshes: Dict[str, Mesh] = {}              # A dictionary of the model's meshes
-        self.shear_walls: Dict[str, ShearWall] = {}    # A dictionary of the model's shear walls
-        self.mats: Dict[str, MatFoundation] = {}       # A dictionary of the model's mat foundations
-        self.load_combos: Dict[str, LoadCombo] = {}    # A dictionary of the model's load combinations
-        self._D: Dict[str, NDArray[float64]] = {}      # A dictionary of the model's nodal displacements by load combination
+        self.nodes: dict[str, Node3D] = {}             # A dictionary of the model's nodes
+        self.materials: dict[str, Material] = {}       # A dictionary of the model's materials
+        self.sections: dict[str, Section] = {}         # A dictonary of the model's cross-sections
+        self.springs: dict[str, Spring3D] = {}         # A dictionary of the model's springs
+        self.members: dict[str, PhysMember] = {}       # A dictionary of the model's physical members
+        self.quads: dict[str, Quad3D] = {}             # A dictionary of the model's quadiralterals
+        self.plates: dict[str, Plate3D] = {}           # A dictionary of the model's rectangular plates
+        self.meshes: dict[str, Mesh] = {}              # A dictionary of the model's meshes
+        self.shear_walls: dict[str, ShearWall] = {}    # A dictionary of the model's shear walls
+        self.mats: dict[str, MatFoundation] = {}       # A dictionary of the model's mat foundations
+        self.load_combos: dict[str, LoadCombo] = {}    # A dictionary of the model's load combinations
+        self._D: dict[str, NDArray[float64]] = {}      # A dictionary of the model's nodal displacements by load combination
 
         self.solution: str | None = None  # Indicates the solution type for the latest run of the model
 
     @property
-    def load_cases(self) -> List[str]:
+    def load_cases(self) -> list[str]:
         """Returns a list of all the load cases in the model (in alphabetical order).
         """
 
         # Create an empty list of load cases
-        cases: List[str] = []
+        cases: list[str] = []
 
         # Step through each node
         for node in self.nodes.values():
@@ -757,7 +757,7 @@ def add_cylinder_mesh(self, name:str, mesh_size:float, radius:float, height:floa
         # Return the mesh's name
         return name
 
-    def add_shear_wall(self, name: str, mesh_size: float, length: float, height: float, thickness: float, material_name: str, ky_mod: float = 0.35, plane: Literal['XY', 'YZ'] = 'XY', origin: List[float] = [0, 0, 0]):
+    def add_shear_wall(self, name: str, mesh_size: float, length: float, height: float, thickness: float, material_name: str, ky_mod: float = 0.35, plane: Literal['XY', 'YZ'] = 'XY', origin: list[float] = [0, 0, 0]):
         """Adds a meshed shear wall helper to the model.
 
         The shear wall utility generates a regular mesh for a rectangular wall panel and
diff --git a/Pynite/FixedEndReactions.py b/Pynite/FixedEndReactions.py
index 71e79509..8e27cadf 100644
--- a/Pynite/FixedEndReactions.py
+++ b/Pynite/FixedEndReactions.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 Created on Fri Nov  3 20:58:03 2017
 
diff --git a/Pynite/LoadCombo.py b/Pynite/LoadCombo.py
index 4ab0eee4..41f293a7 100644
--- a/Pynite/LoadCombo.py
+++ b/Pynite/LoadCombo.py
@@ -5,7 +5,7 @@ class LoadCombo():
     """A class that stores all the information necessary to define a load combination.
     """
 
-    def __init__(self, name: str, combo_tags: List[str] | None = None, factors: Dict[str, float] = {}) -> None:
+    def __init__(self, name: str, combo_tags: list[str] | None = None, factors: dict[str, float] = {}) -> None:
         """Initializes a new load combination.
 
         :param name: A unique name for the load combination.
@@ -17,8 +17,8 @@ def __init__(self, name: str, combo_tags: List[str] | None = None, factors: Dict
         """
         
         self.name: str = name             # A unique user-defined name for the load combination
-        self.combo_tags: List[str] | None = combo_tags   # Used to categorize the load combination (e.g. strength or serviceability)
-        self.factors: Dict[str, float] = factors       # A dictionary containing each load case name and associated load factor
+        self.combo_tags: list[str] | None = combo_tags   # Used to categorize the load combination (e.g. strength or serviceability)
+        self.factors: dict[str, float] = factors       # A dictionary containing each load case name and associated load factor
     
     def AddLoadCase(self, case_name: str, factor: float) -> None:
         '''
diff --git a/Pynite/Member3D.py b/Pynite/Member3D.py
index ad55e380..c0794203 100644
--- a/Pynite/Member3D.py
+++ b/Pynite/Member3D.py
@@ -89,17 +89,17 @@ def __init__(self, model: FEModel3D, name: str, i_node: Node3D,
         self.j_reversal: bool = False
 
         self.rotation: float = rotation  # Member rotation (degrees) about its local x-axis
-        self.PtLoads: List[Tuple] = []  # A list of point loads & moments applied to the element (Direction, P, x, case='Case 1') or (Direction, M, x, case='Case 1')
-        self.DistLoads: List[Tuple] = []       # A list of linear distributed loads applied to the element (Direction, w1, w2, x1, x2, case='Case 1', self_weight=False)
-        self.SegmentsZ: List[BeamSegZ] = []       # A list of mathematically continuous beam segments for z-bending
-        self.SegmentsY: List[BeamSegY] = []       # A list of mathematically continuous beam segments for y-bending
-        self.SegmentsX: List[BeamSegZ] = []       # A list of mathematically continuous beam segments for torsion
-        self.Releases: List[bool] = [False, False, False, False, False, False, False, False, False, False, False, False]
+        self.PtLoads: list[tuple] = []  # A list of point loads & moments applied to the element (Direction, P, x, case='Case 1') or (Direction, M, x, case='Case 1')
+        self.DistLoads: list[tuple] = []       # A list of linear distributed loads applied to the element (Direction, w1, w2, x1, x2, case='Case 1', self_weight=False)
+        self.SegmentsZ: list[BeamSegZ] = []       # A list of mathematically continuous beam segments for z-bending
+        self.SegmentsY: list[BeamSegY] = []       # A list of mathematically continuous beam segments for y-bending
+        self.SegmentsX: list[BeamSegZ] = []       # A list of mathematically continuous beam segments for torsion
+        self.Releases: list[bool] = [False, False, False, False, False, False, False, False, False, False, False, False]
         self.tension_only: bool = tension_only  # Indicates whether the member is tension-only
         self.comp_only: bool = comp_only  # Indicates whether the member is compression-only
 
         # Members need to track whether they are active or not for any given load combination. They may become inactive for a load combination during a tension/compression-only analysis. This dictionary will be used when the model is solved.
-        self.active: Dict[str, bool] = {}  # Key = load combo name, Value = True or False
+        self.active: dict[str, bool] = {}  # Key = load combo name, Value = True or False
 
         # The 'Member3D' object will store results for one load combination at a time. To reduce repetative calculations the '_solved_combo' variable will be used to track whether the member needs to be resegmented before running calculations for any given load combination.
         self._solved_combo: LoadCombo | None = None  # The current solved load combination
@@ -120,7 +120,7 @@ def L(self) -> float:
         return self.i_node.distance(self.j_node)
 
 # %%
-    def _partition_D(self) -> Tuple[List[int], List[int]]:
+    def _partition_D(self) -> tuple[list[int], list[int]]:
         """
         Builds lists of unreleased and released degree of freedom indices for the member.
 
@@ -1154,7 +1154,7 @@ def shear(self, Direction: Literal['Fy', 'Fz'], x: float, combo_name: str = 'Com
 
             return 0
 
-    def max_shear(self, Direction: Literal['Fy', 'Fz'], combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def max_shear(self, Direction: Literal['Fy', 'Fz'], combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the maximum shear force in the member for the specified direction
         and load combination(s).
@@ -1223,7 +1223,7 @@ def max_shear(self, Direction: Literal['Fy', 'Fz'], combo_tags: Union[str, List[
         # Return 0 if no valid combos were found
         return Vmax_global if Vmax_global is not None else 0
 
-    def min_shear(self, Direction: Literal['Fy', 'Fz'], combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def min_shear(self, Direction: Literal['Fy', 'Fz'], combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the minimum shear force in the member for the specified direction
         and load combination(s).
@@ -1436,7 +1436,7 @@ def moment(self, Direction: Literal['My', 'Mz'], x: float, combo_name: str = 'Co
 
             return 0
 
-    def max_moment(self, Direction: Literal['My', 'Mz'], combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def max_moment(self, Direction: Literal['My', 'Mz'], combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the maximum bending moment in the member for the specified direction
         and load combination(s).
@@ -1508,7 +1508,7 @@ def max_moment(self, Direction: Literal['My', 'Mz'], combo_tags: Union[str, List
         # Return 0 if no valid combos were found
         return Mmax_global if Mmax_global is not None else 0
 
-    def min_moment(self, Direction: Literal['My', 'Mz'], combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def min_moment(self, Direction: Literal['My', 'Mz'], combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the minimum bending moment in the member for the specified direction
         and load combination(s).
@@ -1620,7 +1620,7 @@ def plot_moment(self, Direction: Literal['My', 'Mz'], combo_name: str = 'Combo 1
         Member3D.__plt.title('Member ' + self.name + '\n' + combo_name)
         Member3D.__plt.show()
 
-    def moment_array(self, Direction: Literal['My', 'Mz'], n_points: int, combo_name: str = 'Combo 1', x_array: Optional[NDArray[float64]] = None) -> NDArray[float64]:
+    def moment_array(self, Direction: Literal['My', 'Mz'], n_points: int, combo_name: str = 'Combo 1', x_array: NDArray[float64] | None = None) -> NDArray[float64]:
         """
         Returns the array of the moment in the member for the given direction
 
@@ -1708,7 +1708,7 @@ def torque(self, x: float, combo_name: str = 'Combo 1') -> float:
 
             return 0
 
-    def max_torque(self, combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def max_torque(self, combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the maximum torsional moment in the member across the
         specified load combination(s).
@@ -1767,7 +1767,7 @@ def max_torque(self, combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
         return Tmax_global if Tmax_global is not None else 0
 
 
-    def min_torque(self, combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def min_torque(self, combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the minimum torsional moment in the member across the
         specified load combination(s).
@@ -1921,7 +1921,7 @@ def axial(self, x: float, combo_name: str = 'Combo 1') -> float:
 
             return 0
 
-    def max_axial(self, combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def max_axial(self, combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the maximum axial force in the member across the specified
         load combination(s).
@@ -1980,7 +1980,7 @@ def max_axial(self, combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
         # Return the global maximum, or 0 if nothing was found
         return Pmax_global if Pmax_global is not None else 0
 
-    def min_axial(self, combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def min_axial(self, combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the minimum axial force in the member across the specified
         load combination(s).
@@ -2073,7 +2073,7 @@ def plot_axial(self, combo_name: str = 'Combo 1', n_points=20) -> None:
         Member3D.__plt.title('Member ' + self.name + '\n' + combo_name)
         Member3D.__plt.show()    
 
-    def axial_array(self, n_points: int, combo_name: str = 'Combo 1', x_array: Optional[NDArray[float64]] = None) -> NDArray[float64]:
+    def axial_array(self, n_points: int, combo_name: str = 'Combo 1', x_array: NDArray[float64] | None = None) -> NDArray[float64]:
         """
         Returns the array of the axial force in the member for the given direction
         
@@ -2178,7 +2178,7 @@ def deflection(self, Direction: Literal['dx', 'dy', 'dz'], x: float, combo_name:
 
             return 0
 
-    def max_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def max_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the maximum deflection in the member across the specified
         load combination(s).
@@ -2244,7 +2244,7 @@ def max_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_tags: Union
         return dmax_global if dmax_global is not None else 0
 
 
-    def min_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_tags: Union[str, List[str]] = 'Combo 1') -> float:
+    def min_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_tags: str | list[str] = 'Combo 1') -> float:
         """
         Returns the minimum deflection in the member across the specified
         load combination(s).
@@ -2348,7 +2348,7 @@ def plot_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_name: str
         Member3D.__plt.title('Member ' + self.name + '\n' + combo_name)
         Member3D.__plt.show()
 
-    def deflection_array(self, Direction: Literal['dx', 'dy', 'dz'], n_points: int, combo_name: str = 'Combo 1', x_array: Optional[NDArray[float64]] = None) -> NDArray[float64]:
+    def deflection_array(self, Direction: Literal['dx', 'dy', 'dz'], n_points: int, combo_name: str = 'Combo 1', x_array: NDArray[float64] | None = None) -> NDArray[float64]:
         """
         Returns the array of the deflection in the member for the given direction
         
@@ -2508,7 +2508,7 @@ def plot_rel_deflection(self, Direction: Literal['dx', 'dy', 'dz'], combo_name:
         Member3D.__plt.title('Member ' + self.name + '\n' + combo_name)
         Member3D.__plt.show()
 
-    def rel_deflection_array(self, Direction: Literal['dx', 'dy', 'dz'], n_points: int, combo_name: str = 'Combo 1', x_array: Optional[NDArray[float64]] = None) -> NDArray[float64]:
+    def rel_deflection_array(self, Direction: Literal['dx', 'dy', 'dz'], n_points: int, combo_name: str = 'Combo 1', x_array: NDArray[float64] | None = None) -> NDArray[float64]:
         """
         Returns the array of the relative deflection in the member for the given direction
         
@@ -2839,7 +2839,7 @@ def _segment_member(self, combo_name='Combo 1'):
                                 SegmentsY[i].V1 += (f1[2] + f2[2])/2*(x2 - x1)
                                 SegmentsY[i].M1 += (x1 - x2)*(2*f1[2]*x1 - 3*f1[2]*x + f1[2]*x2 + f2[2]*x1 - 3*f2[2]*x + 2*f2[2]*x2)/6
 
-    def _extract_vector_results(self, segments: List, x_array: NDArray[float64], result_name: Literal['moment', 'shear', 'axial', 'torque', 'deflection', 'axial_deflection'], P_delta: bool = False) -> NDArray[float64]:
+    def _extract_vector_results(self, segments: list, x_array: NDArray[float64], result_name: Literal['moment', 'shear', 'axial', 'torque', 'deflection', 'axial_deflection'], P_delta: bool = False) -> NDArray[float64]:
         """
         Extracts result values at specified locations along a structural member using efficient, 
         vectorized evaluation of piecewise segment functions.
diff --git a/Pynite/Mesh.py b/Pynite/Mesh.py
index b85b5f12..811b7d7c 100644
--- a/Pynite/Mesh.py
+++ b/Pynite/Mesh.py
@@ -44,8 +44,8 @@ def __init__(self, thickness: float, material_name: str, model: FEModel3D, kx_mo
         self.last_node = None               # The name of the last node in the mesh
         self.start_element = start_element  # The name of the first element in the mesh
         self.last_element = None            # The name of the last element in the mesh
-        self.nodes: Dict[str, Node3D] = {}  # A dictionary containing the nodes in the mesh
-        self.elements: Dict[str, Union[Quad3D, Plate3D]] = {}  # A dictionary containing the elements in the mesh
+        self.nodes: dict[str, Node3D] = {}  # A dictionary containing the nodes in the mesh
+        self.elements: dict[str, Quad3D | Plate3D] = {}  # A dictionary containing the elements in the mesh
         self.element_type = 'Quad'          # The type of element used in the mesh
         self.is_generated = False          # A flag indicating whether the mesh has been generated
 
@@ -54,8 +54,8 @@ def _rename_duplicates(self) -> None:
         """
 
         # Initialize lists to track node and element name changes
-        revised_nodes: Dict[str, Node3D] = {}
-        revised_elements: Dict[str, Union[Quad3D, Plate3D]] = {}
+        revised_nodes: dict[str, Node3D] = {}
+        revised_elements: dict[str, Quad3D | Plate3D] = {}
 
         # Step through each node in the mesh
         for node in self.nodes.values():
@@ -660,7 +660,7 @@ def min_membrane(self, direction: str = 'Sx', combo_tags: str | list[str] = 'Com
 
 class RectangleMesh(Mesh):
 
-    def __init__(self, mesh_size: float, width: float, height: float, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1.0, ky_mod: float = 1.0, origin: List[float] = [0, 0, 0], plane: str = 'XY', x_control: List[float] | None = None, y_control: List[float] | None = None, start_node: str = 'N1', start_element: str = 'Q1', element_type: str = 'Quad') -> None:
+    def __init__(self, mesh_size: float, width: float, height: float, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1.0, ky_mod: float = 1.0, origin: list[float] = [0, 0, 0], plane: str = 'XY', x_control: list[float] | None = None, y_control: list[float] | None = None, start_node: str = 'N1', start_element: str = 'Q1', element_type: str = 'Quad') -> None:
         """
         A rectangular mesh of elements.
 
@@ -714,7 +714,7 @@ def __init__(self, mesh_size: float, width: float, height: float, thickness: flo
         else: self.y_control = y_control
 
         self.element_type = element_type
-        self.openings: Dict[str, RectOpening] = {}
+        self.openings: dict[str, RectOpening] = {}
     
     def generate(self) -> None:
 
@@ -1054,7 +1054,7 @@ class AnnulusMesh(Mesh):
     """
 
     def __init__(self, mesh_size: float, outer_radius: float, inner_radius: float, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1,
-        ky_mod: float = 1, origin: List[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
+        ky_mod: float = 1, origin: list[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
 
         """Annular (donut) mesh between inner and outer radii.
 
@@ -1179,7 +1179,7 @@ class AnnulusRingMesh(Mesh):
     """
 
     def __init__(self, outer_radius: float, inner_radius: float, num_quads: int, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1,
-                 origin: List[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
+                 origin: list[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
 
         """Single annular ring of quads between two radii.
 
@@ -1331,7 +1331,7 @@ class AnnulusTransRingMesh(Mesh):
     A mesh of quadrilaterals forming an annular ring (a donut) with the mesh transitioning to a finer on the outer edge.
     """
 
-    def __init__(self, outer_radius: float, inner_radius: float, num_inner_quads: int, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1, origin: List[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
+    def __init__(self, outer_radius: float, inner_radius: float, num_inner_quads: int, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1, origin: list[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
         """Creates an annular ring (a donut) with the mesh transitioning to a finer mesh on the outer edge
 
         :param outer_radius: The outer radius of the annular ring.
@@ -1527,7 +1527,7 @@ class FrustrumMesh(AnnulusMesh):
     """
 
     def __init__(self, mesh_size: float, large_radius: float, small_radius: float, height: float, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1,
-                 origin: List[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
+                 origin: list[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1') -> None:
         """Conical frustum mesh generated from an annulus and then tapered to height.
 
         :param mesh_size: Target element size for the base annulus.
@@ -1600,7 +1600,7 @@ def generate(self) -> None:
 #%%
 class CylinderMesh(Mesh):
 
-    def __init__(self, mesh_size: float, radius: float, height: float, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1,origin: List[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1', num_elements: int | None = None, element_type: str = 'Quad') -> None:
+    def __init__(self, mesh_size: float, radius: float, height: float, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1,origin: list[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1', num_elements: int | None = None, element_type: str = 'Quad') -> None:
 
         """Cylindrical shell mesh.
 
@@ -1775,7 +1775,7 @@ class CylinderRingMesh(Mesh):
     """
 
     def __init__(self, radius: float, height: float, num_elements: int, thickness: float, material_name: str, model: FEModel3D, kx_mod: float = 1, ky_mod: float = 1,
-                 origin: List[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1',
+                 origin: list[float] = [0, 0, 0], axis: str = 'Y', start_node: str = 'N1', start_element: str = 'Q1',
                  element_type: str = 'Quad') -> None:
 
         super().__init__(thickness, material_name, model, kx_mod, ky_mod, start_node=start_node, start_element=start_element)
@@ -1917,7 +1917,7 @@ def generate(self) -> None:
         # Flag the mesh as generated
         self.is_generated = True
         
-def check_mesh_integrity(mesh: Mesh, console_log: bool = True) -> Union[str, List[str], None]:
+def check_mesh_integrity(mesh: Mesh, console_log: bool = True) -> str | list[str] | None:
     """Runs basic integrity checks to ensure the mesh is in sync with its model. Usually you don't
     want to run this check unless the mesh has been generated since generating the mesh is what
     syncs it to the model.
diff --git a/Pynite/Node3D.py b/Pynite/Node3D.py
index b474b823..62f7e248 100644
--- a/Pynite/Node3D.py
+++ b/Pynite/Node3D.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 Created on Thu Nov  2 18:04:56 2017
 
@@ -28,29 +27,29 @@ class Node3D():
     def __init__(self, model: FEModel3D, name: str, X: float, Y: float, Z: float):
 
         self.name = name                 # A unique name for the node assigned by the user
-        self.ID: Optional[int] = None    # A unique index number for the node assigned by the program
+        self.ID: int | None = None    # A unique index number for the node assigned by the program
 
         self.X = X          # Global X coordinate
         self.Y = Y          # Global Y coordinate
         self.Z = Z          # Global Z coordinate
 
-        self.NodeLoads: List[Tuple[str, float, str]] = []  # A list of loads applied to the node (Direction, P, case) or (Direction, M, case)
+        self.NodeLoads: list[tuple[str, float, str]] = []  # A list of loads applied to the node (Direction, P, case) or (Direction, M, case)
 
         # Initialize the dictionaries of calculated node displacements
-        self.DX: Dict[str, float] = {}
-        self.DY: Dict[str, float] = {}
-        self.DZ: Dict[str, float] = {}
-        self.RX: Dict[str, float] = {}
-        self.RY: Dict[str, float] = {}
-        self.RZ: Dict[str, float] = {}
+        self.DX: dict[str, float] = {}
+        self.DY: dict[str, float] = {}
+        self.DZ: dict[str, float] = {}
+        self.RX: dict[str, float] = {}
+        self.RY: dict[str, float] = {}
+        self.RZ: dict[str, float] = {}
 
         # Initialize the dictionaries of calculated node reactions
-        self.RxnFX: Dict[str, float] = {}
-        self.RxnFY: Dict[str, float] = {}
-        self.RxnFZ: Dict[str, float] = {}
-        self.RxnMX: Dict[str, float] = {}
-        self.RxnMY: Dict[str, float] = {}
-        self.RxnMZ: Dict[str, float] = {}
+        self.RxnFX: dict[str, float] = {}
+        self.RxnFY: dict[str, float] = {}
+        self.RxnFZ: dict[str, float] = {}
+        self.RxnMX: dict[str, float] = {}
+        self.RxnMY: dict[str, float] = {}
+        self.RxnMZ: dict[str, float] = {}
 
         # Initialize all support conditions to `False`
         self.support_DX: bool = False
@@ -61,12 +60,12 @@ def __init__(self, model: FEModel3D, name: str, X: float, Y: float, Z: float):
         self.support_RZ: bool = False
 
         # Inititialize all support springs
-        self.spring_DX: List[float | str | bool | None] = [None, None, None]  # [stiffness, direction, active]
-        self.spring_DY: List[float | str | bool | None] = [None, None, None]
-        self.spring_DZ: List[float | str | bool | None] = [None, None, None]
-        self.spring_RX: List[float | str | bool | None] = [None, None, None]
-        self.spring_RY: List[float | str | bool | None] = [None, None, None]
-        self.spring_RZ: List[float | str | bool | None] = [None, None, None]
+        self.spring_DX: list[float | str | bool | None] = [None, None, None]  # [stiffness, direction, active]
+        self.spring_DY: list[float | str | bool | None] = [None, None, None]
+        self.spring_DZ: list[float | str | bool | None] = [None, None, None]
+        self.spring_RX: list[float | str | bool | None] = [None, None, None]
+        self.spring_RY: list[float | str | bool | None] = [None, None, None]
+        self.spring_RZ: list[float | str | bool | None] = [None, None, None]
 
         # Initialize all enforced displacements to `None`
         self.EnforcedDX: float | None = None
@@ -77,14 +76,14 @@ def __init__(self, model: FEModel3D, name: str, X: float, Y: float, Z: float):
         self.EnforcedRZ: float | None = None
 
         # Initialize the color contour value for the node. This will be used for contour smoothing.
-        self.contour: List[float] = []
+        self.contour: list[float] = []
 
         # The 'Node3D' object will store results for one load combination at a time. 
 
         # Adding a link to the model that Nodes belong to 
         self.model: FEModel3D = model
 
-    def distance(self, other: 'Node3D') -> float:
+    def distance(self, other: Node3D) -> float:
         """
         Returns the distance to another node.
 
diff --git a/Pynite/PhysMember.py b/Pynite/PhysMember.py
index b8beedbf..0639752e 100644
--- a/Pynite/PhysMember.py
+++ b/Pynite/PhysMember.py
@@ -29,7 +29,7 @@ def __init__(self, model: FEModel3D, name: str, i_node: Node3D, j_node: Node3D,
                  tension_only: bool = False, comp_only: bool = False) -> None:
 
         super().__init__(model, name, i_node, j_node, material_name, section_name, rotation, tension_only, comp_only)
-        self.sub_members: Dict[str, Member3D] = {}
+        self.sub_members: dict[str, Member3D] = {}
 
     def descritize(self) -> None:
         """
@@ -40,7 +40,7 @@ def descritize(self) -> None:
         self.sub_members = {}
 
         # Start a new list of nodes along the member
-        int_nodes: List[Tuple[Node3D, float]] = []
+        int_nodes: list[tuple[Node3D, float]] = []
 
         # Create a vector from the i-node to the j-node
         Xi, Yi, Zi = self.i_node.X, self.i_node.Y, self.i_node.Z
@@ -1012,7 +1012,7 @@ def deflection_array(self, Direction: Literal['dx', 'dy', 'dz'], n_points: int,
         # Return the results
         return d_array2
 
-    def find_member(self, x: float) -> Tuple[Member3D, float]:
+    def find_member(self, x: float) -> tuple[Member3D, float]:
         """
         Returns the sub-member that the physical member's local point 'x' lies on, and 'x' modified for that sub-member's local coordinate system.
         """
diff --git a/Pynite/Plate3D.py b/Pynite/Plate3D.py
index 03b3f056..93f8f739 100644
--- a/Pynite/Plate3D.py
+++ b/Pynite/Plate3D.py
@@ -46,7 +46,7 @@ def __init__(self, name: str, i_node: Node3D, j_node: Node3D, m_node: Node3D, n_
         """
 
         self.name: str = name
-        self.ID: Optional[int] = None
+        self.ID: int | None = None
         self.type: str = 'Rect'
 
         self.i_node: Node3D = i_node
@@ -59,7 +59,7 @@ def __init__(self, name: str, i_node: Node3D, j_node: Node3D, m_node: Node3D, n_
         self.kx_mod: float = kx_mod
         self.ky_mod: float = ky_mod
 
-        self.pressures: List[Tuple[float, str]] = []  # A list of surface pressures [pressure, case='Case 1']
+        self.pressures: list[tuple[float, str]] = []  # A list of surface pressures [pressure, case='Case 1']
 
         # Plates need a link to the model they belong to
         self.model: FEModel3D = model
diff --git a/Pynite/Quad3D.py b/Pynite/Quad3D.py
index d6615acf..37b0922e 100644
--- a/Pynite/Quad3D.py
+++ b/Pynite/Quad3D.py
@@ -32,7 +32,7 @@ def __init__(self, name: str, i_node: Node3D, j_node: Node3D, m_node: Node3D, n_
                  ky_mod: float = 1.0):
 
         self.name: str = name
-        self.ID: Optional[int] = None
+        self.ID: int | None = None
         self.type: str = 'Quad'
 
         self.i_node: Node3D = i_node
@@ -44,7 +44,7 @@ def __init__(self, name: str, i_node: Node3D, j_node: Node3D, m_node: Node3D, n_
         self.kx_mod: float = kx_mod
         self.ky_mod: float = ky_mod
 
-        self.pressures: List[Tuple[float, str]] = []  # A list of surface pressures [pressure, case='Case 1']
+        self.pressures: list[tuple[float, str]] = []  # A list of surface pressures [pressure, case='Case 1']
 
         # Quads need a link to the model they belong to
         self.model: FEModel3D = model
@@ -151,7 +151,7 @@ def L_k(self, k: Literal[5, 6, 7, 8]) -> float:
         else:
             raise Exception('Invalid value for k. k must be 5, 6, 7, or 8.')
 
-    def dir_cos(self, k: Literal[5, 6, 7, 8]) -> Tuple[float, float]:
+    def dir_cos(self, k: Literal[5, 6, 7, 8]) -> tuple[float, float]:
 
         L_k = self.L_k(k)
 
diff --git a/Pynite/Rendering.py b/Pynite/Rendering.py
index 4793839c..89eb4b0f 100644
--- a/Pynite/Rendering.py
+++ b/Pynite/Rendering.py
@@ -1,7 +1,8 @@
 from __future__ import annotations # Allows more recent type hints features
 from json import load
 import warnings
-from typing import TYPE_CHECKING, Callable, List, Any
+from typing import TYPE_CHECKING, List, Any
+from collections.abc import Callable
 
 from IPython.display import Image
 import numpy as np
@@ -28,7 +29,7 @@ class Renderer:
     """Used to render finite element models.
     """
 
-    scalar: Optional[str] = None
+    scalar: str | None = None
 
     def __init__(self, model: FEModel3D) -> None:
 
@@ -40,9 +41,9 @@ def __init__(self, model: FEModel3D) -> None:
         self._deformed_scale: float = 30.0
         self._render_nodes: bool = True
         self._render_loads: bool = True
-        self._color_map: Optional[str] = None
-        self._combo_name: Optional[str] = 'Combo 1'
-        self._case: Optional[str] = None
+        self._color_map: str | None = None
+        self._combo_name: str | None = 'Combo 1'
+        self._case: str | None = None
         self._labels: bool = True
         self._scalar_bar: bool = False
         self._scalar_bar_text_size: int = 24
@@ -52,7 +53,7 @@ def __init__(self, model: FEModel3D) -> None:
         # This is added because `self.update()` clears the plotter, removing user self.plotter configurations.
         # Functions in this list run after Pynite adds actors, allowing further PyVista customizations 
         # (e.g., grid, axes) before render. Each func in this list must accept a `pyvista.Plotter` argument.
-        self.post_update_callbacks: List[Callable[[pv.Plotter], None]] = []
+        self.post_update_callbacks: list[Callable[[pv.Plotter], None]] = []
 
         self.plotter: pv.Plotter = pv.Plotter()
         self.plotter.set_background('white')  # Setting background color
@@ -63,12 +64,12 @@ def __init__(self, model: FEModel3D) -> None:
         self.plotter.set_viewup((0, 1, 0))  # Set the Y axis to vertical for 3D plots
 
         # Initialize load labels
-        self._load_label_points: List[List[float]] = []
-        self._load_labels: List[Union[str, float, int]] = []
+        self._load_label_points: list[list[float]] = []
+        self._load_labels: list[str | float | int] = []
 
         # Initialize spring labels
-        self._spring_label_points: List[List[float]] = []
-        self._spring_labels: List[str] = []
+        self._spring_label_points: list[list[float]] = []
+        self._spring_labels: list[str] = []
 
     @property
     def window_width(self) -> int:
@@ -129,28 +130,28 @@ def render_loads(self, render_loads: bool) -> None:
         self._render_loads = render_loads
     
     @property
-    def color_map(self) -> Optional[str]:
+    def color_map(self) -> str | None:
         return self._color_map
     
     @color_map.setter
-    def color_map(self, color_map: Optional[str]) -> None:
+    def color_map(self, color_map: str | None) -> None:
         self._color_map = color_map
     
     @property
-    def combo_name(self) -> Optional[str]:
+    def combo_name(self) -> str | None:
         return self._combo_name
     
     @combo_name.setter
-    def combo_name(self, combo_name: Optional[str]) -> None:
+    def combo_name(self, combo_name: str | None) -> None:
         self._combo_name = combo_name
         self._case = None
 
     @property
-    def case(self) -> Optional[str]:
+    def case(self) -> str | None:
         return self._case
 
     @case.setter
-    def case(self, case: Optional[str]) -> None:
+    def case(self, case: str | None) -> None:
         self._case = case
         self._combo_name = None
 
@@ -621,7 +622,7 @@ def plot_spring(self, spring: Spring3D, color: str = 'grey', deformed: bool = Fa
         self._spring_label_points.append([(Xi + Xj) / 2, (Yi + Yj) / 2, (Zi + Zj) / 2])
 
             
-    def plot_plates(self, deformed_shape: bool, deformed_scale: float, color_map: Optional[str], combo_name: Optional[str]) -> None:
+    def plot_plates(self, deformed_shape: bool, deformed_scale: float, color_map: str | None, combo_name: str | None) -> None:
         
         # Start a list of vertices
         plate_vertices = []
@@ -776,8 +777,8 @@ def plot_deformed_member(self, member: Member3D, scale_factor: float) -> None:
                 line = pv.Line(D_plot[i], D_plot[i+1])
                 self.plotter.add_mesh(line, color='red', line_width=2)
 
-    def plot_pt_load(self, position: Tuple[float, float, float], direction: Union[Tuple[float, float, float], np.ndarray], 
-                    length: float, label_text: Optional[Union[str, float, int]] = None, color: str = 'green') -> None:
+    def plot_pt_load(self, position: tuple[float, float, float], direction: tuple[float, float, float] | np.ndarray, 
+                    length: float, label_text: str | float | int | None = None, color: str = 'green') -> None:
 
         # Create a unit vector in the direction of the 'direction' vector
         unitVector = direction/np.linalg.norm(direction)
@@ -814,9 +815,9 @@ def plot_pt_load(self, position: Tuple[float, float, float], direction: Union[Tu
         # Plot the shaft
         self.plotter.add_mesh(shaft, line_width=2, color=color)                         
 
-    def plot_dist_load(self, position1: Tuple[float, float, float], position2: Tuple[float, float, float], 
-                      direction: Union[np.ndarray, Tuple[float, float, float]], length1: float, length2: float,
-                      label_text1: Optional[Union[str, float, int]], label_text2: Optional[Union[str, float, int]], 
+    def plot_dist_load(self, position1: tuple[float, float, float], position2: tuple[float, float, float], 
+                      direction: np.ndarray | tuple[float, float, float], length1: float, length2: float,
+                      label_text1: str | float | int | None, label_text2: str | float | int | None, 
                       color: str = 'green') -> None:
 
         # Calculate the length of the distributed load
@@ -867,8 +868,8 @@ def plot_dist_load(self, position1: Tuple[float, float, float], position2: Tuple
         # Combine all geometry into a single PolyData object
         self.plotter.add_mesh(tail_line, color=color)
 
-    def plot_moment(self, center: Tuple[float, float, float], direction: Union[Tuple[float, float, float], np.ndarray], 
-                    radius: float, label_text: Optional[Union[str, float, int]] = None, color: str = 'green') -> None:
+    def plot_moment(self, center: tuple[float, float, float], direction: tuple[float, float, float] | np.ndarray, 
+                    radius: float, label_text: str | float | int | None = None, color: str = 'green') -> None:
 
         # Convert the direction vector into a unit vector
         v1 = direction/np.linalg.norm(direction)
diff --git a/Pynite/Section.py b/Pynite/Section.py
index e24d3015..909e1d5e 100644
--- a/Pynite/Section.py
+++ b/Pynite/Section.py
@@ -14,7 +14,7 @@ class Section():
 
     This class stores all properties related to the geometry of the member
     """
-    def __init__(self, model: 'FEModel3D', name: str, A: float, Iy: float, Iz: float, J: float) -> None:
+    def __init__(self, model: FEModel3D, name: str, A: float, Iy: float, Iz: float, J: float) -> None:
         """
         :param model: The finite element model to which this section belongs
         :type model: FEModel3D
@@ -29,7 +29,7 @@ def __init__(self, model: 'FEModel3D', name: str, A: float, Iy: float, Iz: float
         :param J: The torsion constant of the section
         :type J: float
         """        
-        self.model: 'FEModel3D' = model
+        self.model: FEModel3D = model
         self.name: str = name
         self.A: float = A
         self.Iy: float = Iy
@@ -81,7 +81,7 @@ def G(self, fx: float, my: float, mz: float) -> NDArray:
 
 class SteelSection(Section):
 
-    def __init__(self, model: 'FEModel3D', name: str, A: float, Iy: float, Iz: float, J: float, 
+    def __init__(self, model: FEModel3D, name: str, A: float, Iy: float, Iz: float, J: float, 
                  Zy: float, Zz: float, material_name: str) -> None:
         """
         Initialize a steel section
diff --git a/Pynite/ShearWall.py b/Pynite/ShearWall.py
index ca36259d..10845ad0 100644
--- a/Pynite/ShearWall.py
+++ b/Pynite/ShearWall.py
@@ -29,15 +29,15 @@ def __init__(self, model, name, mesh_size, length, height, thickness, material_n
         self.ky_mod = ky_mod
         self.origin = origin
         self.plane = plane
-        self._openings: List[List[str | float | None]] = []
-        self._flanges: List[List[str | float]] = []
-        self._supports: List[List[float]] = []
-        self._stories: List[List[str | float]] = []
-        self._shears: List[List[str | float]] = []
-        self._axials: List[List[str | float]] = []
-        self._materials: List[List[str | float]] = []
-        self.piers: Dict[str, Pier] = {}
-        self.coupling_beams: Dict[str, CouplingBeam] = {}
+        self._openings: list[list[str | float | None]] = []
+        self._flanges: list[list[str | float]] = []
+        self._supports: list[list[float]] = []
+        self._stories: list[list[str | float]] = []
+        self._shears: list[list[str | float]] = []
+        self._axials: list[list[str | float]] = []
+        self._materials: list[list[str | float]] = []
+        self.piers: dict[str, Pier] = {}
+        self.coupling_beams: dict[str, CouplingBeam] = {}
         self.asign_material(material_name, thickness)
 
     def asign_material(self, name: str, t: float, x_start: float | None = None, x_end: float | None = None, y_start: float | None = None, y_end: float | None = None) -> None:
@@ -96,8 +96,8 @@ def add_axial(self, story_name: str, force: float, case: str = 'Case 1') -> None
     def generate(self) -> None:
 
         # Identify mesh control points
-        x_control: List[float] = [0, self.L]
-        y_control: List[float] = [0, self.H]
+        x_control: list[float] = [0, self.L]
+        y_control: list[float] = [0, self.H]
 
         for material in self._materials:
             x_control.append(material[2])
@@ -105,7 +105,7 @@ def generate(self) -> None:
             y_control.append(material[4])
             y_control.append(material[5])
 
-        z_control: List[float] = [0]
+        z_control: list[float] = [0]
         for flg in self._flanges:
             if flg[6] == '+z': z_control.append(flg[1])
             else: z_control.append(-flg[1])
@@ -300,8 +300,8 @@ def _identify_piers(self) -> None:
         self.piers = {}
 
         # Create a list of x and y coordinates that represent the edges of the wall
-        x_vals: List[float] = [0, self.L]
-        y_vals: List[float] = [0, self.H]
+        x_vals: list[float] = [0, self.L]
+        y_vals: list[float] = [0, self.H]
 
         # Add the edges of the openings to the lists
         for opng in self._openings:
@@ -315,7 +315,7 @@ def _identify_piers(self) -> None:
         y_vals = sorted(y_vals)
 
         # Remove duplicate (or near duplicate) values
-        unique_list: List[float] = []
+        unique_list: list[float] = []
         for i in range(len(x_vals) - 1):
             # Only keep the value at `i` if it's not a duplicate or near duplicate of the next value
             if not isclose(x_vals[i], x_vals[i+1]):
@@ -341,7 +341,7 @@ def _identify_piers(self) -> None:
             self.piers['P' + str(i+1)] = Pier('P' + str(i+1), x, y, width, height, self)
 
         # Divide the strip piers further into rectanglular piers using the top and bottom of each opening as pier boundaries
-        new_piers: Dict[str, Pier] = {}
+        new_piers: dict[str, Pier] = {}
         pier_count = 1
         for pier in self.piers.values():
             for i in range(len(y_vals) - 1):
@@ -354,7 +354,7 @@ def _identify_piers(self) -> None:
         self.piers = new_piers
 
         # Delete any piers that fall within an opening
-        delete_list: List[str] = []
+        delete_list: list[str] = []
         for pier in self.piers.values():
 
             # Check if this pier is inside any of the openings
@@ -431,7 +431,7 @@ def _identify_piers(self) -> None:
                     break
 
         # Generate a list of new keys in ascending order
-        new_keys: List[str] = [f'P{i+1}' for i in range(len(self.piers))]
+        new_keys: list[str] = [f'P{i+1}' for i in range(len(self.piers))]
 
         # Replace the old dicionary with one that has updated keys
         self.piers = dict(zip(new_keys, self.piers.values()))
@@ -463,8 +463,8 @@ def _identify_coupling_beams(self) -> None:
         self.coupling_beams = {}
 
         # Create a list of x and y coordinates that represent the edges of the wall
-        x_vals: List[float] = [0, self.L]
-        y_vals: List[float] = [0, self.H]
+        x_vals: list[float] = [0, self.L]
+        y_vals: list[float] = [0, self.H]
 
         # Add the edges of the openings to the lists
         for opng in self._openings:
@@ -478,7 +478,7 @@ def _identify_coupling_beams(self) -> None:
         y_vals = sorted(y_vals)
 
         # Remove duplicate (or near duplicate) values
-        unique_list: List[float] = []
+        unique_list: list[float] = []
         for i in range(len(x_vals) - 1):
             # Only keep the value at `i` if it's not a duplicate or near duplicate of the next value
             if not isclose(x_vals[i], x_vals[i+1]):
@@ -504,7 +504,7 @@ def _identify_coupling_beams(self) -> None:
             self.coupling_beams['B' + str(i+1)] = CouplingBeam('B' + str(i+1), x, y, length, height, self)
 
         # Divide the strips further into rectanglular beams using the left and right of each opening as beam boundaries
-        new_beams: Dict[str, CouplingBeam] = {}
+        new_beams: dict[str, CouplingBeam] = {}
         beam_count = 1
         for beam in self.coupling_beams.values():
             for i in range(len(x_vals) - 1):
@@ -517,7 +517,7 @@ def _identify_coupling_beams(self) -> None:
         self.coupling_beams = new_beams
 
         # Delete any beams that fall within an opening
-        delete_list: List[str] = []
+        delete_list: list[str] = []
         for beam in self.coupling_beams.values():
            
            # Check if this beam is inside any of the openings
@@ -604,7 +604,7 @@ def _identify_coupling_beams(self) -> None:
             del self.coupling_beams[beam]
 
         # Generate a list of new keys in ascending order
-        new_keys: List[str] = [f'B{i + 1}' for i in range(len(self.coupling_beams))]
+        new_keys: list[str] = [f'B{i + 1}' for i in range(len(self.coupling_beams))]
 
         # Replace the old dicionary with one that has updated keys
         self.coupling_beams = dict(zip(new_keys, self.coupling_beams.values()))
@@ -815,7 +815,7 @@ def print_coupling_beams(self, combo_name: str = 'Combo 1') -> None:
         print('+----------------------------+')
         print(table)
 
-    def _local2global(self, x: float, y: float, z: float, plane: Literal['XY', 'XZ', 'YZ'] = 'XY') -> List[float]:
+    def _local2global(self, x: float, y: float, z: float, plane: Literal['XY', 'XZ', 'YZ'] = 'XY') -> list[float]:
 
         Xo, Yo, Zo = self.origin[0], self.origin[1], self.origin[2]
 
@@ -835,7 +835,7 @@ def _local2global(self, x: float, y: float, z: float, plane: Literal['XY', 'XZ',
         return [X, Y, Z]
 
 
-def _global2local(X: float, Y: float, Z: float, origin: List[float] = [0, 0, 0], plane: Literal['XY', 'YZ', 'XZ'] = 'XY') -> List[float]:
+def _global2local(X: float, Y: float, Z: float, origin: list[float] = [0, 0, 0], plane: Literal['XY', 'YZ', 'XZ'] = 'XY') -> list[float]:
 
     Xo, Yo, Zo = origin[0], origin[1], origin[2]
 
@@ -873,9 +873,9 @@ def __init__(self, name: str, x: float, y: float, width: float, height: float, s
         self.origin = shear_wall.origin
 
         # This list will be used by the parent shear wall to store a list of only the plates in this pier
-        self.plates: List[Quad3D] = []
+        self.plates: list[Quad3D] = []
 
-    def sum_forces(self, combo_name: str = 'Combo 1') -> Tuple[float, float, float, float]:
+    def sum_forces(self, combo_name: str = 'Combo 1') -> tuple[float, float, float, float]:
 
         # Initialize the forces in the plate
         P, M, V = 0, 0, 0
@@ -933,9 +933,9 @@ def __init__(self, name: str, x: float, y: float, length: float, height: float,
         self.plane = shear_wall.plane
         self.origin = shear_wall.origin
 
-        self.plates: List[Quad3D] = []
+        self.plates: list[Quad3D] = []
 
-    def sum_forces(self, combo_name: str = 'Combo 1') -> Tuple[float, float, float, float]:
+    def sum_forces(self, combo_name: str = 'Combo 1') -> tuple[float, float, float, float]:
 
         # Initialize plate forces to zero
         P, M, V = 0, 0, 0
diff --git a/Pynite/Spring3D.py b/Pynite/Spring3D.py
index 01312579..bc64d41b 100644
--- a/Pynite/Spring3D.py
+++ b/Pynite/Spring3D.py
@@ -23,24 +23,24 @@ class Spring3D():
 
 #%%
     def __init__(self, name: str, i_node: Node3D, j_node: Node3D, ks: float, 
-                 LoadCombos: Dict[str, LoadCombo] = {'Combo 1': LoadCombo('Combo 1', factors={'Case 1': 1.0})},
+                 LoadCombos: dict[str, LoadCombo] = {'Combo 1': LoadCombo('Combo 1', factors={'Case 1': 1.0})},
                  tension_only: bool = False, comp_only: bool = False) -> None:
         '''
         Initializes a new spring.
         '''
         self.name: str = name      # A unique name for the spring given by the user
-        self.ID: Optional[int] = None        # Unique index number for the spring assigned by the program
+        self.ID: int | None = None        # Unique index number for the spring assigned by the program
         self.i_node: Node3D = i_node  # The spring's i-node
         self.j_node: Node3D = j_node  # The spring's j-node
         self.ks: float = ks          # The spring constant (force/displacement)
-        self.load_combos: Dict[str, LoadCombo] = LoadCombos # The dictionary of load combinations in the model this spring belongs to
+        self.load_combos: dict[str, LoadCombo] = LoadCombos # The dictionary of load combinations in the model this spring belongs to
         self.tension_only: bool = tension_only # Indicates whether the spring is tension-only
         self.comp_only: bool = comp_only # Indicates whether the spring is compression-only
 
         # Springs need to track whether they are active or not for any given load combination.
         # They may become inactive for a load combination during a tension/compression-only
         # analysis. This dictionary will be used when the model is solved.
-        self.active: Dict[str, bool] = {} # Key = load combo name, Value = True or False
+        self.active: dict[str, bool] = {} # Key = load combo name, Value = True or False
 
 #%%
     def L(self) -> float:
diff --git a/Pynite/VTKWriter.py b/Pynite/VTKWriter.py
index d4848b49..8a46c606 100644
--- a/Pynite/VTKWriter.py
+++ b/Pynite/VTKWriter.py
@@ -77,7 +77,7 @@ def _write_node_data(self, path:str):
         node_names.SetName("Name")
 
         node_cells = vtk.vtkCellArray()
-        node_ids: Dict[str, int] = {}
+        node_ids: dict[str, int] = {}
         for node in self.model.nodes.values():
             point_id = points.InsertNextPoint(node.X, node.Y, node.Z)
             node_ids[node.name] = point_id
@@ -127,8 +127,8 @@ def _write_node_data(self, path:str):
                 moments.InsertTuple3(node_id, node.RxnMX[combo], node.RxnMY[combo], node.RxnMZ[combo])
 
                 # calculate the NodeLoad for each node
-                fl: Dict[str,float] = {"X":0,"Y":0,"Z":0}
-                ml: Dict[str,float] = {"X":0,"Y":0,"Z":0}
+                fl: dict[str,float] = {"X":0,"Y":0,"Z":0}
+                ml: dict[str,float] = {"X":0,"Y":0,"Z":0}
                 for (f_or_m, direction), magnitude, case in node.NodeLoads:
                     if f_or_m == "F":
                         fl[direction] += magnitude
@@ -167,7 +167,7 @@ def _write_member_data(self, path: str):
         #### CREATE LINE CELLS ####
         # each (sub)member is further subdivided into line segments
         lines = vtk.vtkCellArray()
-        submembers: List[Tuple[Tuple[float, float], Member3D, vtk.vtkLine]] = []
+        submembers: list[tuple[tuple[float, float], Member3D, vtk.vtkLine]] = []
         for member in self.model.members.values():
             if len(member.sub_members) == 0:
                 # The model has not been analyzed yet. Only add straight lines between the nodes
@@ -332,7 +332,7 @@ def _write_quad_data(self, path: str):
         quads = vtk.vtkCellArray() # this holds all cell data. A cell is in this case the vtkBiQuadraticQuads
 
         # keeps track of all vtk subquads for every Pynite Quad by name
-        quad_references: Dict[str, List[vtk.vtkBiQuadraticQuad]] = {}
+        quad_references: dict[str, list[vtk.vtkBiQuadraticQuad]] = {}
 
         # this loops through every Pynite quad and creates 4 vtkBiQuadraticQuad elements and their point data
         for quad in self.model.quads.values():
diff --git a/Pynite/Visualization.py b/Pynite/Visualization.py
index b93926d7..c9068851 100644
--- a/Pynite/Visualization.py
+++ b/Pynite/Visualization.py
@@ -1350,17 +1350,17 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
 
                 # Display the load
                 if load[0] == 'FX':
-                    ptLoad = VisPtLoad((node.X - 0.6*annotation_size*sign, node.Y, node.Z), [1, 0, 0], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad((node.X - 0.6*annotation_size*sign, node.Y, node.Z), [1, 0, 0], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'FY':
-                    ptLoad = VisPtLoad((node.X, node.Y - 0.6*annotation_size*sign, node.Z), [0, 1, 0], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad((node.X, node.Y - 0.6*annotation_size*sign, node.Z), [0, 1, 0], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'FZ':
-                    ptLoad = VisPtLoad((node.X, node.Y, node.Z - 0.6*annotation_size*sign), [0, 0, 1], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad((node.X, node.Y, node.Z - 0.6*annotation_size*sign), [0, 0, 1], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'MX':
-                    ptLoad = VisMoment((node.X, node.Y, node.Z), (1*sign, 0, 0), abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment((node.X, node.Y, node.Z), (1*sign, 0, 0), abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'MY':
-                    ptLoad = VisMoment((node.X, node.Y, node.Z), (0, 1*sign, 0), abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment((node.X, node.Y, node.Z), (0, 1*sign, 0), abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'MZ':
-                    ptLoad = VisMoment((node.X, node.Y, node.Z), (0, 0, 1*sign), abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment((node.X, node.Y, node.Z), (0, 0, 1*sign), abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
 
                 polydata.AddInputData(ptLoad.polydata.GetOutput())
                 renderer.AddActor(ptLoad.lblActor)
@@ -1391,29 +1391,29 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
 
                 # Display the load
                 if load[0] == 'Fx':
-                    ptLoad = VisPtLoad(position, dir_cos[0, :], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad(position, dir_cos[0, :], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'Fy':
-                    ptLoad = VisPtLoad(position, dir_cos[1, :], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad(position, dir_cos[1, :], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'Fz':
-                    ptLoad = VisPtLoad(position, dir_cos[2, :], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad(position, dir_cos[2, :], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'Mx':
-                    ptLoad = VisMoment(position, dir_cos[0, :]*sign, abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment(position, dir_cos[0, :]*sign, abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'My':
-                    ptLoad = VisMoment(position, dir_cos[1, :]*sign, abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment(position, dir_cos[1, :]*sign, abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'Mz':
-                    ptLoad = VisMoment(position, dir_cos[2, :]*sign, abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment(position, dir_cos[2, :]*sign, abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'FX':
-                    ptLoad = VisPtLoad(position, [1, 0, 0], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad(position, [1, 0, 0], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'FY':
-                    ptLoad = VisPtLoad(position, [0, 1, 0], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad(position, [0, 1, 0], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'FZ':
-                    ptLoad = VisPtLoad(position, [0, 0, 1], load_value/max_pt_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisPtLoad(position, [0, 0, 1], load_value/max_pt_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'MX':
-                    ptLoad = VisMoment(position, [1*sign, 0, 0], abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment(position, [1*sign, 0, 0], abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'MY':
-                    ptLoad = VisMoment(position, [0, 1*sign, 0], abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment(position, [0, 1*sign, 0], abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
                 elif load[0] == 'MZ':
-                    ptLoad = VisMoment(position, [0, 0, 1*sign], abs(load_value)/max_moment*2.5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                    ptLoad = VisMoment(position, [0, 0, 1*sign], abs(load_value)/max_moment*2.5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
 
                 polydata.AddInputData(ptLoad.polydata.GetOutput())
                 renderer.AddActor(ptLoad.lblActor)
@@ -1437,17 +1437,17 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
 
                 # Display the load
                 if load[0] == 'Fx':
-                    distLoad = VisDistLoad(position1, position2, dir_cos[0, :], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, '{:.3g}'.format(w1), '{:.3g}'.format(w2), annotation_size)
+                    distLoad = VisDistLoad(position1, position2, dir_cos[0, :], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, f'{w1:.3g}', f'{w2:.3g}', annotation_size)
                 elif load[0] == 'Fy':
-                    distLoad = VisDistLoad(position1, position2, dir_cos[1, :], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, '{:.3g}'.format(w1), '{:.3g}'.format(w2), annotation_size)
+                    distLoad = VisDistLoad(position1, position2, dir_cos[1, :], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, f'{w1:.3g}', f'{w2:.3g}', annotation_size)
                 elif load[0] == 'Fz':
-                    distLoad = VisDistLoad(position1, position2, dir_cos[2, :], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, '{:.3g}'.format(w1), '{:.3g}'.format(w2), annotation_size)
+                    distLoad = VisDistLoad(position1, position2, dir_cos[2, :], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, f'{w1:.3g}', f'{w2:.3g}', annotation_size)
                 elif load[0] == 'FX':
-                    distLoad = VisDistLoad(position1, position2, [1, 0, 0], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, '{:.3g}'.format(w1), '{:.3g}'.format(w2), annotation_size)
+                    distLoad = VisDistLoad(position1, position2, [1, 0, 0], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, f'{w1:.3g}', f'{w2:.3g}', annotation_size)
                 elif load[0] == 'FY':
-                    distLoad = VisDistLoad(position1, position2, [0, 1, 0], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, '{:.3g}'.format(w1), '{:.3g}'.format(w2), annotation_size)
+                    distLoad = VisDistLoad(position1, position2, [0, 1, 0], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, f'{w1:.3g}', f'{w2:.3g}', annotation_size)
                 elif load[0] == 'FZ':
-                    distLoad = VisDistLoad(position1, position2, [0, 0, 1], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, '{:.3g}'.format(w1), '{:.3g}'.format(w2), annotation_size)
+                    distLoad = VisDistLoad(position1, position2, [0, 0, 1], w1/max_dist_load*5*annotation_size, w2/max_dist_load*5*annotation_size, f'{w1:.3g}', f'{w2:.3g}', annotation_size)
 
                 polydata.AddInputData(distLoad.polydata.GetOutput())
                 renderer.AddActor(distLoad.lblActors[0])
@@ -1485,7 +1485,7 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
                 position3 = [plate.n_node.X, plate.n_node.Y, plate.n_node.Z]
 
                 # Create an area load and get its data
-                area_load = VisAreaLoad(position0, position1, position2, position3, dir_cos*sign, abs(load_value)/max_area_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
+                area_load = VisAreaLoad(position0, position1, position2, position3, dir_cos*sign, abs(load_value)/max_area_load*5*annotation_size, f'{load_value:.3g}', annotation_size, theme)
 
                 # Add the area load's arrows to the overall load polydata
                 polydata.AddInputData(area_load.polydata.GetOutput())
diff --git a/Testing/test_2D_frames.py b/Testing/test_2D_frames.py
index 84c48fac..bcdf275a 100644
--- a/Testing/test_2D_frames.py
+++ b/Testing/test_2D_frames.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_AISC_PDelta_benchmarks.py b/Testing/test_AISC_PDelta_benchmarks.py
index 42bff1bf..79dd6f93 100644
--- a/Testing/test_AISC_PDelta_benchmarks.py
+++ b/Testing/test_AISC_PDelta_benchmarks.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_TC_analysis.py b/Testing/test_TC_analysis.py
index ef44e617..fdff7adf 100644
--- a/Testing/test_TC_analysis.py
+++ b/Testing/test_TC_analysis.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_end_releases.py b/Testing/test_end_releases.py
index 1ba61f89..bc8a8c84 100644
--- a/Testing/test_end_releases.py
+++ b/Testing/test_end_releases.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_loads.py b/Testing/test_loads.py
index 962c0760..da20f76d 100644
--- a/Testing/test_loads.py
+++ b/Testing/test_loads.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_plates&quads.py b/Testing/test_plates&quads.py
index 5f367dc5..4a9b0cbe 100644
--- a/Testing/test_plates&quads.py
+++ b/Testing/test_plates&quads.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_reactions.py b/Testing/test_reactions.py
index 28639a0a..cd6e1f87 100644
--- a/Testing/test_reactions.py
+++ b/Testing/test_reactions.py
@@ -1,5 +1,3 @@
-# -*- coding: utf-8 -*-
-
 from Pynite import FEModel3D
 
 
diff --git a/Testing/test_reanalysis.py b/Testing/test_reanalysis.py
index 3f47911c..1a9e6008 100644
--- a/Testing/test_reanalysis.py
+++ b/Testing/test_reanalysis.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_sloped_beam.py b/Testing/test_sloped_beam.py
index f09e5c7c..2abe1d68 100644
--- a/Testing/test_sloped_beam.py
+++ b/Testing/test_sloped_beam.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_spring_support.py b/Testing/test_spring_support.py
index 84156d6d..7e35d9a9 100644
--- a/Testing/test_spring_support.py
+++ b/Testing/test_spring_support.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_springs.py b/Testing/test_springs.py
index 59d8a725..f427bf6d 100644
--- a/Testing/test_springs.py
+++ b/Testing/test_springs.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 from Pynite import FEModel3D
 import numpy as np
 
diff --git a/Testing/test_support_settlement.py b/Testing/test_support_settlement.py
index d7762c11..13a4d2de 100644
--- a/Testing/test_support_settlement.py
+++ b/Testing/test_support_settlement.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_torsion.py b/Testing/test_torsion.py
index 6f09f542..9cb2193a 100644
--- a/Testing/test_torsion.py
+++ b/Testing/test_torsion.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_unstable_structure.py b/Testing/test_unstable_structure.py
index 158be3e9..e5281b54 100644
--- a/Testing/test_unstable_structure.py
+++ b/Testing/test_unstable_structure.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 MIT License
 
diff --git a/Testing/test_vtkwriter.py b/Testing/test_vtkwriter.py
index f0d32c13..52dfe992 100644
--- a/Testing/test_vtkwriter.py
+++ b/Testing/test_vtkwriter.py
@@ -1,4 +1,3 @@
-
 import unittest
 from Pynite import FEModel3D
 from Pynite.VTKWriter import VTKWriter
diff --git a/setup.py b/setup.py
index b1369025..4aa4e3b2 100644
--- a/setup.py
+++ b/setup.py
@@ -1,6 +1,6 @@
 import setuptools
 
-with open("README.md", "r") as fh:
+with open("README.md") as fh:
     long_description = fh.read()
 
 setuptools.setup(