Repository: ross-rotordynamics/ross
Pull request #502 (2020-03-26 12:29:34)

Updating to match master current build


Patch:
diff --git a/.gitignore b/.gitignore
index 81034dd3..4e8193c7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -111,8 +111,23 @@ Thumbs.db
 notebooks/
 
 # docs #
+########
 /docs/examples/*_files
 /docs/generated/
 /paper/whedon/
 /paper/paper.pdf
 /paper/Makefile
+
+### VirtualEnv ###
+##################
+.Python
+pyvenv.cfg
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+pip-selfcheck.json
+
diff --git a/docs/refs.bib b/docs/refs.bib
index 6d2cd068..1290c2b4 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -27,6 +27,13 @@
   publisher={Cambridge university press}
 }
 
+@book{lalanne2001,
+  title={Rotordynamics Prediction in Engineering},
+  author={M Lalanne and G Ferraris},
+  year={1998},
+  publisher={John Wiley and Sons
+}
+
 @article{zang2005review,
   title={A review of robust optimal design and its application in dynamics},
   author={Zang, C and Friswell, MI and Mottershead, JE},
@@ -627,3 +634,10 @@ year={2011}
 	year={2012}
 }
 
+@article{Hutchingson2001,
+	title={Shear Coefficients for Timoshenko Beam Theory},
+	author={Hutchingson, J. R.},
+	journal={J. Appl. Mech.},
+	volume={68},
+	year={2001}
+}
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index f1c7b782..95fdafc1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
+
 numpy
 scipy
 matplotlib
@@ -5,4 +6,4 @@ toml
 pandas>=0.23
 bokeh
 xlrd
-pint
\ No newline at end of file
+pint
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 25981fb4..33524cac 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -21,6 +21,7 @@ __all__ = [
     "SealElement",
     "BallBearingElement",
     "RollerBearingElement",
+    "BearingElement6DoF",
     "MagneticBearingElement",
 ]
 # fmt: on
@@ -418,7 +419,7 @@ class BearingElement(Element):
         return dict(x_0=0, y_0=1)
 
     def dof_global_index(self):
-        """Get the global index for a element specific degree of freedom."""
+        """Get the global index for an element specific degree of freedom."""
         global_index = super().dof_global_index()
 
         if self.n_link is not None:
@@ -1453,6 +1454,357 @@ class MagneticBearingElement(BearingElement):
         )
 
 
+class BearingElement6DoF(BearingElement):
+    """A generalistic 6 DoF bearing element. This class will create a bearing
+    element based on the user supplied stiffness and damping coefficients. These
+    are determined alternatively, via purposefully built codes.
+
+    Parameters
+    ----------
+    kxx: float, optional
+        Direct stiffness in the x direction.
+        Default is None
+    kxy: float, optional
+        Cross stiffness between xy directions.
+        Default is None
+    kyx: float, optional
+        Cross stiffness between yx directions.
+        Default is None
+    kyy: float, optional
+        Direct stiffness in the y direction.
+        Default is None
+    kzz: float, optional
+        Direct stiffness in the z direction.
+        Default is None
+    cxx: float, optional
+        Direct damping in the x direction.
+        Default is None
+    cxy: float, optional
+        Cross damping between xy directions.
+        Default is None
+    cyx: float, optional
+        Cross damping between yx directions.
+        Default is None
+    cyy: float, optional
+        Direct damping in the y direction.
+        Default is None
+    czz: float, optional
+        Direct damping in the z direction.
+        Default is None
+    tag : str, optional
+        A tag to name the element
+        Default is None
+    Examples
+    --------
+    >>> n = 0
+    >>> kxx = 1.0e7
+    >>> kyy = 1.5e7
+    >>> kzz = 5.0e5
+    >>> bearing = BearingElement6DoF(n=n, kxx=kxx, kyy=kyy, kzz=kzz,
+    ...                              cxx=0, cyy=0)
+    >>> bearing.K(0)
+    array([[10000000.,        0.,        0.],
+           [       0., 15000000.,        0.],
+           [       0.,        0.,   500000.]])
+    """
+
+    def __init__(
+        self,
+        n,
+        kxx,
+        kyy,
+        cxx,
+        cyy,
+        kxy=0.0,
+        kyx=0.0,
+        kzz=0.0,
+        cxy=0.0,
+        cyx=0.0,
+        czz=0.0,
+        tag=None,
+        scale_factor=1,
+        color="#355d7a",
+    ):
+        super().__init__(
+            n=n,
+            kxx=kxx,
+            cxx=cxx,
+            kyy=kyy,
+            kxy=kxy,
+            kyx=kyx,
+            cyy=cyy,
+            cxy=cxy,
+            cyx=cyx,
+            frequency=None,
+            tag=tag,
+            n_link=None,
+            scale_factor=scale_factor,
+            color=color,
+        )
+
+        new_args = ["kzz", "czz"]
+
+        args_dict = locals()
+        coefficients = {}
+
+        if kzz is None:
+            args_dict["kzz"] = (
+                kxx * 0.6
+            )  # NSK manufacturer sugestion for deep groove ball bearings
+        if czz is None:
+            args_dict["czz"] = cxx
+
+        for arg in new_args:
+            if arg[0] == "k":
+                coefficients[arg] = _Stiffness_Coefficient(
+                    coefficient=args_dict[arg], w=None
+                )
+            else:
+                coefficients[arg] = _Damping_Coefficient(args_dict[arg], None)
+
+        coefficients_len = [len(v.coefficient) for v in coefficients.values()]
+
+        for c in coefficients_len:
+            if c != 1:
+                raise ValueError(
+                    "Arguments (coefficients and frequency)"
+                    " must have the same dimension"
+                )
+
+        for k, v in coefficients.items():
+            setattr(self, k, v)
+
+    def __repr__(self):
+        """This function returns a string representation of a bearing element.
+        Parameters
+        ----------
+
+        Returns
+        -------
+        A string representation of a bearing object.
+        Examples
+        --------
+        >>> bearing = bearing_example()
+        >>> bearing # doctest: +ELLIPSIS
+        BearingElement(n=0, n_link=None,
+         kxx=[...
+        """
+        return (
+            f"{self.__class__.__name__}"
+            f"(n={self.n}, n_link={self.n_link},\n"
+            f" kxx={self.kxx}, kxy={self.kxy},\n"
+            f" kyx={self.kyx}, kyy={self.kyy},\n"
+            f" kzz={self.kzz}, cxx={self.cxx},\n"
+            f" cxy={self.cxy}, cyx={self.cyx},\n"
+            f" cyy={self.cyy}, czz={self.czz},\n"
+            f" frequency={self.frequency}, tag={self.tag!r})"
+        )
+
+    def __eq__(self, other):
+        """This function allows bearing elements to be compared.
+        Parameters
+        ----------
+        other: object
+            The second object to be compared with.
+
+        Returns
+        -------
+        bool
+            True if the comparison is true; False otherwise.
+        Examples
+        --------
+        >>> bearing1 = bearing_example()
+        >>> bearing2 = bearing_example()
+        >>> bearing1 == bearing2
+        True
+        """
+        compared_attributes = [
+            "kxx",
+            "kyy",
+            "kxy",
+            "kyx",
+            "cxx",
+            "cyy",
+            "cxy",
+            "cyx",
+            "kzz",
+            "czz",
+            "frequency",
+            "n",
+        ]
+        if isinstance(other, self.__class__):
+            return all(
+                (
+                    np.array(getattr(self, attr)).all()
+                    == np.array(getattr(other, attr)).all()
+                    for attr in compared_attributes
+                )
+            )
+        return False
+
+    def save(self, file_name=Path(os.getcwd())):
+        """Saves a bearing element in a toml format. It works as an auxiliary
+        function of the save function in the Rotor class.
+
+        Parameters
+        ----------
+        file_name: string
+            The name of the file the bearing element will be saved in.
+
+        Returns
+        -------
+        None
+
+        Examples
+        --------
+        >>> bearing = bearing_example()
+        >>> bearing.save(Path(os.getcwd()))
+        """
+        data = self.get_data(Path(file_name) / "BearingElement6DoF.toml")
+
+        if type(self.frequency) == np.ndarray:
+            try:
+                self.frequency[0]
+                frequency = list(self.frequency)
+            except IndexError:
+                frequency = None
+
+        data["BearingElement6DoF"][str(self.n)] = {
+            "n": self.n,
+            "kxx": self.kxx.coefficient,
+            "cxx": self.cxx.coefficient,
+            "kyy": self.kyy.coefficient,
+            "kxy": self.kxy.coefficient,
+            "kyx": self.kyx.coefficient,
+            "kzz": self.kzz.coefficient,
+            "cyy": self.cyy.coefficient,
+            "cxy": self.cxy.coefficient,
+            "cyx": self.cyx.coefficient,
+            "czz": self.czz.coefficient,
+            "frequency": frequency,
+            "tag": self.tag,
+        }
+        self.dump_data(data, Path(file_name) / "BearingElement6DoF.toml")
+
+    @staticmethod
+    def load(file_name=""):
+        """Loads a list of bearing elements saved in a toml format.
+
+        Parameters
+        ----------
+        file_name: string
+            The name of the file of the bearing element to be loaded.
+
+        Returns
+        -------
+        A list of bearing elements.
+
+        Examples
+        --------
+
+        """
+        bearing_elements = []
+        bearing_elements_dict = BearingElement.get_data(
+            file_name=Path(file_name) / "BearingElement6DoF.toml"
+        )
+        for element in bearing_elements_dict["BearingElement6DoF"]:
+            bearing = BearingElement6DoF(
+                **bearing_elements_dict["BearingElement6DoF"][element]
+            )
+            data = bearing_elements_dict["BearingElement6DoF"]
+            bearing.kxx.coefficient = data[element]["kxx"]
+
+            bearing.kxy.coefficient = data[element]["kxy"]
+
+            bearing.kyx.coefficient = data[element]["kyx"]
+
+            bearing.kyy.coefficient = data[element]["kyy"]
+
+            bearing.kzz.coefficient = data[element]["kzz"]
+
+            bearing.cxx.coefficient = data[element]["cxx"]
+
+            bearing.cxy.coefficient = data[element]["cxy"]
+
+            bearing.cyx.coefficient = data[element]["cyx"]
+
+            bearing.cyy.coefficient = data[element]["cyy"]
+
+            bearing.czz.coefficient = data[element]["czz"]
+
+            bearing_elements.append(bearing)
+        return bearing_elements
+
+    def dof_mapping(self):
+        return dict(x_0=0, y_0=1, z_0=2)
+
+    def dof_global_index(self):
+        """Get the global index for an element specific degree of freedom for the 6DoF element."""
+        global_index = super().dof_global_index()
+
+        if self.n_link is not None:
+            global_index = global_index._asdict()
+            global_index[f"x_{self.n_link}"] = 6 * self.n_link
+            global_index[f"y_{self.n_link}"] = 6 * self.n_link + 1
+            global_index[f"z_{self.n_link}"] = 6 * self.n_link + 2
+            dof_tuple = namedtuple("GlobalIndex", global_index)
+            global_index = dof_tuple(**global_index)
+
+        return global_index
+
+    def K(self, frequency):
+        """Returns the stiffness matrix for a given excitation frequency.
+
+        Returns
+        -------
+        K: np.ndarray
+            A 3x3 matrix of floats containing the kxx, kxy, kyx, kyy and kzz values.
+
+        Examples
+        --------
+        >>> bearing = bearing_6dof_example()
+        >>> bearing.K(0)
+        array([[1000000.,       0.,       0.],
+               [      0.,  800000.,       0.],
+               [      0.,       0.,  100000.]])
+        """
+        kxx = self.kxx.interpolated(frequency)
+        kyy = self.kyy.interpolated(frequency)
+        kxy = self.kxy.interpolated(frequency)
+        kyx = self.kyx.interpolated(frequency)
+        kzz = self.kzz.interpolated(frequency)
+
+        K = np.array([[kxx, kxy, 0], [kyx, kyy, 0], [0, 0, kzz]])
+
+        return K
+
+    def C(self, frequency):
+        """Returns the damping matrix for a given excitation frequency.
+
+        Returns
+        -------
+        C: np.ndarray
+            A 3x3 matrix of floats containing the cxx, cxy, cyx, cyy, and czz values.
+
+        Examples
+        --------
+        >>> bearing = bearing_6dof_example()
+        >>> bearing.C(0)
+        array([[200.,   0.,   0.],
+               [  0., 150.,   0.],
+               [  0.,   0.,  50.]])
+        """
+        cxx = self.cxx.interpolated(frequency)
+        cyy = self.cyy.interpolated(frequency)
+        cxy = self.cxy.interpolated(frequency)
+        cyx = self.cyx.interpolated(frequency)
+        czz = self.czz.interpolated(frequency)
+
+        C = np.array([[cxx, cxy, 0], [cyx, cyy, 0], [0, 0, czz]])
+
+        return C
+
+
 def bearing_example():
     """This function returns an instance of a simple bearing.
     The purpose is to make available a simple model
@@ -1505,3 +1857,33 @@ def seal_example():
     w = np.linspace(0, 200, 11)
     seal = SealElement(n=0, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, frequency=w)
     return seal
+
+
+def bearing_6dof_example():
+    """This function returns an instance of a simple bearing.
+    The purpose is to make available a simple model
+    so that doctest can be written using it.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    An instance of a bearing object.
+
+    Examples
+    --------
+    >>> bearing = bearing_example()
+    >>> bearing.frequency[0]
+    0.0
+    """
+    kxx = 1e6
+    kyy = 0.8e6
+    kzz = 1e5
+    cxx = 2e2
+    cyy = 1.5e2
+    czz = 0.5e2
+    bearing = BearingElement6DoF(
+        n=0, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, kzz=kzz, czz=czz
+    )
+    return bearing
diff --git a/ross/disk_element.py b/ross/disk_element.py
index eb869324..0363e8df 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -10,7 +10,7 @@ from bokeh.models import ColumnDataSource, HoverTool
 from ross.element import Element
 from ross.utils import read_table_file
 
-__all__ = ["DiskElement"]
+__all__ = ["DiskElement", "DiskElement6DoF"]
 bokeh_colors = bp.RdGy[11]
 
 
@@ -557,3 +557,227 @@ def disk_example():
     """
     disk = DiskElement(0, 32.589_727_65, 0.178_089_28, 0.329_563_62)
     return disk
+
+
+class DiskElement6DoF(DiskElement):
+    """A disk element for 6 DoFs.
+
+     This class will create a disk element with 6 DoF from input data of inertia and mass.
+
+     Parameters
+     ----------
+     n: int
+         Node in which the disk will be inserted.
+     m : float
+         Mass of the disk element.
+     Id : float
+         Diametral moment of inertia.
+     Ip : float
+         Polar moment of inertia
+     tag : str, optional
+         A tag to name the element
+         Default is None
+
+     Examples
+     --------
+     >>> disk = DiskElement6DoF(n=0, m=32, Id=0.2, Ip=0.3)
+     >>> disk.Ip
+     0.3
+     """
+
+    def dof_mapping(self):
+        """6DoFs degrees of freedom mapping.
+
+        Returns a dictionary with a mapping between degree of freedom and its
+        index.
+
+        Returns
+        -------
+        dof_mapping: dict
+            A dictionary containing the degrees of freedom and their indexes.
+
+        Examples
+        --------
+        The numbering of the degrees of freedom for each node. 
+        
+        Being the following their ordering for a node:
+
+        x_0,u_0 - horizontal translation, 
+        y_0,v_0 - vertical translation, 
+        z_0,w_0 - axial translation, 
+        theta_0 - rotation around horizontal, 
+        psi_0   - rotation around vertical, 
+        phi_0   - torsion around axial,
+        """
+        return dict(u_0=0, v_0=1, w_0=2, theta_0=3, psi_0=4, phi_0=5,)
+
+    def M(self):
+        """6DoFs mass matrix.
+
+        This method will return the mass matrix for an instance of a disk
+        element with 6DoFs.
+
+        Returns
+        -------
+        Mass matrix for the 6DoFs disk element.
+
+        Examples
+        --------
+        >>> disk = DiskElement6DoF(0, 32.58972765, 0.17808928, 0.32956362)
+        >>> disk.M().round(2)
+        array([[32.59,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
+               [ 0.  , 32.59,  0.  ,  0.  ,  0.  ,  0.  ],
+               [ 0.  ,  0.  , 32.59,  0.  ,  0.  ,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.18,  0.  ,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.18,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.33]])
+        """
+        m = self.m
+        Id = self.Id
+        Ip = self.Ip
+        # fmt: off
+        M = np.array([
+            [  m,  0,  0,  0,  0,  0],
+            [  0,  m,  0,  0,  0,  0],
+            [  0,  0,  m,  0,  0,  0],
+            [  0,  0,  0, Id,  0,  0],
+            [  0,  0,  0,  0, Id,  0],
+            [  0,  0,  0,  0,  0, Ip]])
+        # fmt: on
+        return M
+
+    def K(self):
+        """6DoFs stiffness matrix.
+
+        This method will return the stiffness matrix for an instance of a disk
+        element with 6DoFs.
+
+        Returns
+        -------
+        K: np.ndarray
+            A matrix of floats containing the values of the stiffness matrix.
+
+        Examples
+        --------
+        >>> disk = disk_example_6dof()
+        >>> disk.K().round(2)
+        array([[0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
+               [0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
+               [0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
+               [0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
+               [0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
+               [0.  , 0.  , 0.  , 0.33, 0.  , 0.  ]])
+        """
+        Ip = self.Ip
+        # fmt: off
+        K = np.array([
+            [ 0, 0, 0,  0, 0, 0],
+            [ 0, 0, 0,  0, 0, 0],
+            [ 0, 0, 0,  0, 0, 0],
+            [ 0, 0, 0,  0, 0, 0],
+            [ 0, 0, 0,  0, 0, 0],
+            [ 0, 0, 0, Ip, 0, 0]])
+        # fmt: on
+        return K
+
+    def C(self):
+        """6DoFs damping matrix.
+
+        Returns
+        -------
+        C: np.ndarray
+            A matrix of floats containing the values of the damping matrix.
+
+        Examples
+        --------
+        >>> disk = disk_example_6dof()
+        >>> disk.C()
+        array([[0., 0., 0., 0., 0., 0.],
+               [0., 0., 0., 0., 0., 0.],
+               [0., 0., 0., 0., 0., 0.],
+               [0., 0., 0., 0., 0., 0.],
+               [0., 0., 0., 0., 0., 0.],
+               [0., 0., 0., 0., 0., 0.]])
+        """
+        # fmt: off
+        C = np.zeros((6, 6))
+        # fmt: on
+        return C
+
+    def G(self):
+        """6DoFs gyroscopic matrix.
+
+        This method will return the gyroscopic matrix for an instance of a disk
+        element.
+
+        Returns
+        -------
+        G: np.ndarray
+            Gyroscopic matrix for the disk element.
+
+        Examples
+        --------
+        >>> disk = DiskElement6DoF(0, 32.58972765, 0.17808928, 0.32956362)
+        >>> disk.G().round(2)
+        array([[ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.33,  0.  ],
+               [ 0.  ,  0.  ,  0.  , -0.33,  0.  ,  0.  ],
+               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ]])
+        """
+        Ip = self.Ip
+        # fmt: off
+        G = np.array([
+            [ 0, 0, 0,   0,  0, 0],
+            [ 0, 0, 0,   0,  0, 0],
+            [ 0, 0, 0,   0,  0, 0],
+            [ 0, 0, 0,   0, Ip, 0],
+            [ 0, 0, 0, -Ip,  0, 0],
+            [ 0, 0, 0,   0,  0, 0]])
+        # fmt: on
+        return G
+
+
+def disk_example():
+    """This function returns an instance of a simple disk.
+    The purpose is to make available a simple model
+    so that doctest can be written using it.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    An instance of a disk object.
+
+    Examples
+    --------
+    >>> disk = disk_example()
+    >>> disk.Ip
+    0.32956362
+    """
+    disk = DiskElement(0, 32.589_727_65, 0.178_089_28, 0.329_563_62)
+    return disk
+
+
+def disk_example_6dof():
+    """This function returns an instance of a simple disk.
+    The purpose is to make available a simple 6DoFs model
+    so that doctest can be written using it.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    An instance of a disk object.
+
+    Examples
+    --------
+    >>> disk = disk_example()
+    >>> disk.Ip
+    0.32956362
+    """
+    disk = DiskElement6DoF(0, 32.589_727_65, 0.178_089_28, 0.329_563_62)
+    return disk
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index d305f3ad..5eec886b 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -12,7 +12,7 @@ from ross.materials import Material, steel
 from ross.units import check_units
 from ross.utils import read_table_file
 
-__all__ = ["ShaftElement"]
+__all__ = ["ShaftElement", "ShaftElement6DoF"]
 bokeh_colors = bp.RdGy[11]
 
 
@@ -25,7 +25,7 @@ class ShaftElement(Element):
 
     .. math::
 
-        [x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]^T
+        [x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]**T
     Where :math:`\alpha_0` and :math:`\alpha_1` are the bending on the yz plane
     and :math:`\beta_0` and :math:`\beta_1` are the bending on the xz plane.
 
@@ -246,7 +246,7 @@ class ShaftElement(Element):
         self.beam_cg = L * c1 / (4 * c2)
         self.axial_cg_pos = None
 
-        # Slenderness ratio of beam elements (G*A*L^2) / (E*I)
+        # Slenderness ratio of beam elements (G*A*L**2) / (E*I)
         sld = (self.material.G_s * self.A * self.L ** 2) / (self.material.E * Ie)
         self.slenderness_ratio = sld
 
@@ -1191,3 +1191,953 @@ class ShaftElement(Element):
         ]
 
         return elements
+
+
+class ShaftElement6DoF(ShaftElement):
+
+    r"""A 6 Degrees of Freedom shaft element.
+    =========================================
+    This class will create a shaft element that takes into
+    account shear stress, rotary inertia and gyroscopic effects.
+    The matrices will be defined considering the following local
+    coordinate vector:
+
+    :math:`[u_0, v_0, w_0, \theta_0, \psi_0, \phi_0, u_1, v_1, w_1, \theta_1, \psi_1, \phi_1]**T`
+    
+    Being the following their ordering for an element:
+
+        :math:`x_0,u_0`  - horizontal translation;
+        :math:`y_0,v_0`  - vertical translation;
+        :math:`z_0,w_0`  - axial translation;
+        :math:`\theta_0` - rotation around horizontal, bending on the yz plane;
+        :math:`\psi_0`   - rotation around vertical, bending on the xz plane;
+        :math:`\phi_0`   - torsion around axial, z direction.
+
+    Parameters
+    ----------
+    L : float, pint.Quantity
+        Element length.
+    idl : float, pint.Quantity
+        Inner diameter of the element at the left node.
+    odl : float, pint.Quantity
+        Outer diameter of the element at the left node.
+    idr : float, pint.Quantity, optional
+        Inner diameter of the element at the right node;
+        Default is equal to idl value for cylindrical element.
+    odr : float, pint.Quantity, optional
+        Outer diameter of the element at the right node;
+        Default is equal to odl value for cylindrical element.
+    material : ross.material
+        Shaft material.
+    alpha : float, optional
+        Proportional damping coefficient, associated to the element Mass matrix
+    beta : float, optional
+        Proportional damping coefficient, associated to the element Stiffness matrix
+    n : int, optional
+        Element number, coincident with it's first node.
+        If not given, it will be set when the rotor is assembled
+        according to the element's position in the list supplied to
+        the rotor constructor.
+    axial_force : float, optional
+        Axial force;
+        Default is zero.
+    torque : float, optional
+        Torque moment;
+        Default is zero.
+    shear_effects : bool, optional
+        Determine if shear effects are taken into account;
+        Default is True.
+    rotary_inertia : bool, optional
+        Determine if rotary_inertia effects are taken into account;
+        Default is True.
+    gyroscopic : bool, optional
+        Determine if gyroscopic effects are taken into account;
+        Default is True.
+    tag : str, optional
+        Element tag;
+        Default is None.
+
+    Returns
+    -------
+    A 6 degrees of freedom shaft element, with available gyroscopic, shear and rotary inertia effects.
+
+    Attributes
+    ----------
+    Poisson : float
+        Poisson coefficient for the element.
+    kappa : float
+        Shear coefficient for the element, determined from :cite:`Hutchingson2001` formulation.
+
+    References
+    ----------
+    .. bibliography:: ../../../docs/refs.bib
+
+    Examples
+    --------
+    >>> from ross.materials import steel
+    >>> eixo1 = ShaftElement6DoF(L=0.5, idl=0.0, odl=0.01, idr=0.0, odr=0.01, material=steel, n=0, axial_force=10, torque=30)
+    >>> Euler_Bernoulli_Element = ShaftElement6DoF(
+    ...                         material=steel, L=0.5, idl=0.05, odl=0.1,
+    ...                         idr=0.05, odr=0.15,
+    ...                         alpha=0.01, beta=100,
+    ...                         rotary_inertia=False,
+    ...                         shear_effects=False)
+    >>> Euler_Bernoulli_Element.kappa
+    0.7099387976608923
+    """
+
+    @check_units
+    def __init__(
+        self,
+        L,
+        idl,
+        odl,
+        idr=None,
+        odr=None,
+        material=None,
+        n=None,
+        axial_force=0,
+        torque=0,
+        shear_effects=True,
+        rotary_inertia=True,
+        gyroscopic=True,
+        alpha=0,
+        beta=0,
+        tag=None,
+    ):
+
+        if idr is None:
+            idr = idl
+        if odr is None:
+            odr = odl
+
+        if material is None:
+            raise AttributeError("Material is not defined.")
+
+        if type(material) is str:
+            # os.chdir(Path(os.path.dirname(ross.__file__)))
+            self.material = Material.use_material(material)
+        else:
+            self.material = material
+
+        self.shear_effects = shear_effects
+        self.rotary_inertia = rotary_inertia
+        self.gyroscopic = gyroscopic
+        self.axial_force = axial_force
+        self.torque = torque
+        self._n = n
+        self.n_l = n
+        self.n_r = None
+        if n is not None:
+            self.n_r = n + 1
+
+        self.tag = tag
+
+        self.L = float(L)
+        self.o_d = (float(odl) + float(odr)) / 2
+        self.i_d = (float(idl) + float(idr)) / 2
+        self.idl = float(idl)
+        self.odl = float(odl)
+        self.idr = float(idr)
+        self.odr = float(odr)
+        self.color = self.material.color
+
+        self.alpha = float(alpha)
+        self.beta = float(beta)
+
+        # Timoshenko kappa factor determination, based on the diameters relation
+        if self.__is_circular():
+            kappa = (6 * (1 + self.material.Poisson) ** 2) / (
+                7 + 12 * self.material.Poisson + 4 * self.material.Poisson ** 2
+            )
+        elif self.__is_thickwall():
+            a = (self.idl + self.idr) / 2
+            b = (self.odl + self.odr) / 2
+            v = self.material.Poisson
+            kappa = (6 * (a ** 2 + b ** 2) ** 2 * (1 + v) ** 2) / (
+                7 * a ** 4
+                + 34 * a ** 2 * b ** 2
+                + 7 * b ** 4
+                + v * (12 * a ** 4 + 48 * a ** 2 * b ** 2 + 12 * b ** 4)
+                + v ** 2 * (4 * a ** 4 + 16 * a ** 2 * b ** 2 + 4 * b ** 4)
+            )
+        else:
+            kappa = (1 + self.material.Poisson) / (2 + self.material.Poisson)
+
+        self.kappa = kappa
+
+    def __is_circular(self):
+        return self.idl == 0 and self.idr == 0
+
+    def __is_thickwall(self):
+        p = (((self.odl + self.odr) / 2) - ((self.idl + self.idr) / 2)) / (
+            (self.odl + self.odr) / 2
+        )
+        return p >= 0.2
+
+    def save(self, file_name=Path(os.getcwd())):
+        """Save shaft elements to toml file.
+
+        Parameters
+        ----------
+        file_name : str
+
+        Examples
+        --------
+
+        """
+
+        data = self.get_data(Path(file_name) / "ShaftElement6DoF.toml")
+        data["ShaftElement6DoF"][str(self.n)] = {
+            "L": self.L,
+            "idl": self.idl,
+            "odl": self.odl,
+            "idr": self.idr,
+            "odr": self.odr,
+            "alpha": self.alpha,
+            "beta": self.beta,
+            "material": self.material.name,
+            "n": self.n,
+            "axial_force": self.axial_force,
+            "torque": self.torque,
+            "shear_effects": self.shear_effects,
+            "rotary_inertia": self.rotary_inertia,
+            "gyroscopic": self.gyroscopic,
+        }
+        self.dump_data(data, Path(file_name) / "ShaftElement6DoF.toml")
+
+    @staticmethod
+    def load(file_name="ShaftElement6DoF"):
+        """Load previously saved shaft elements from toml file.
+
+        Parameters
+        ----------
+        file_name : str, optional
+
+        Examples
+        --------
+
+        """
+        shaft_elements = []
+        with open("ShaftElement6DoF.toml", "r") as f:
+            shaft_elements_dict = toml.load(f)
+            for element in shaft_elements_dict["ShaftElement6DoF"]:
+                shaft_elements.append(
+                    ShaftElement6DoF(**shaft_elements_dict["ShaftElement6DoF"][element])
+                )
+        return shaft_elements
+
+    @classmethod
+    def from_table(cls, file, sheet_type="Simple", sheet_name=0):
+        """Instantiate one or more shafts using inputs from an Excel table.
+
+        A header with the names of the columns is required. These names should
+        match the names expected by the routine (usually the names of the
+        parameters, but also similar ones). The program will read every row
+        bellow the header until they end or it reaches a NaN.
+
+        Parameters
+        ----------
+        file: str
+            Path to the file containing the shaft parameters.
+        sheet_type: str, optional
+            Describes the kind of sheet the function should expect:
+                Simple: The input table should specify only the number of the materials to be used.
+                They must be saved prior to calling the method.
+                Model: The materials parameters must be passed along with the shaft parameters. Each
+                material must have an id number and each shaft must reference one of the materials ids.
+        sheet_name: int or str, optional
+            Position of the sheet in the file (starting from 0) or its name. If none is passed, it is
+            assumed to be the first sheet in the file.
+
+        Returns
+        -------
+        shaft: list
+            A list of shaft objects.
+        """
+        parameters = read_table_file(
+            file, "shaft", sheet_name=sheet_name, sheet_type=sheet_type
+        )
+        list_of_shafts = []
+        if sheet_type == "Model":
+            new_materials = {}
+            for i in range(0, len(parameters["matno"])):
+                new_material = Material(
+                    name="shaft_mat_" + str(parameters["matno"][i]),
+                    rho=parameters["rhoa"][i],
+                    E=parameters["ea"][i],
+                    G_s=parameters["ga"][i],
+                )
+                new_materials["shaft_mat_" + str(parameters["matno"][i])] = new_material
+            for i in range(0, len(parameters["L"])):
+                list_of_shafts.append(
+                    cls(
+                        L=parameters["L"][i],
+                        i_d=parameters["i_d"][i],
+                        o_d=parameters["o_d"][i],
+                        alpha=parameters["alpha"][i],
+                        beta=parameters["beta"][i],
+                        material=new_materials[parameters["material"][i]],
+                        n=parameters["n"][i],
+                        axial_force=parameters["axial_force"][i],
+                        torque=parameters["torque"][i],
+                        gyroscopic=parameters["gyroscopic"][i],
+                    )
+                )
+        elif sheet_type == "Simple":
+            for i in range(0, len(parameters["L"])):
+                list_of_shafts.append(
+                    cls(
+                        L=parameters["L"][i],
+                        i_d=parameters["i_d"][i],
+                        o_d=parameters["o_d"][i],
+                        alpha=parameters["alpha"][i],
+                        beta=parameters["beta"][i],
+                        material=parameters["material"][i],
+                        n=parameters["n"][i],
+                        axial_force=parameters["axial_force"][i],
+                        torque=parameters["torque"][i],
+                        gyroscopic=parameters["gyroscopic"][i],
+                    )
+                )
+        return list_of_shafts
+
+    def dof_mapping(self):
+        """
+        Method to map the element's 12 degrees of freedom
+
+        Parameters
+        ----------
+        Returns
+        -------
+        The numbering of the degrees of freedom for each element node. 
+        
+        Being the following their ordering for an element:
+
+        x_0,u_0 - horizontal translation, 
+        y_0,v_0 - vertical translation, 
+        z_0,w_0 - axial translation, 
+        theta_0 - rotation around horizontal, 
+        psi_0   - rotation around vertical, 
+        phi_0   - torsion around axial
+        x_1,u_1 - horizontal translation, 
+        y_1,v_1 - vertical translation, 
+        z_1,w_1 - axial translation, 
+        theta_1 - rotation around horizontal, 
+        psi_1   - rotation around vertical, 
+        phi_1   - torsion around axial
+        """
+        return dict(
+            u_0=0,
+            v_0=1,
+            w_0=2,
+            theta_0=3,
+            psi_0=4,
+            phi_0=5,
+            u_1=6,
+            v_1=7,
+            w_1=8,
+            theta_1=9,
+            psi_1=10,
+            phi_1=11,
+        )
+
+    def __repr__(self):
+        """This function returns a string representation of a 6 DoF shaft element.
+
+        Parameters
+        ----------
+        Returns
+        -------
+        A string representation of a 6 DoF shaft object.
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> shaft1 = ShaftElement(
+        ...        L=0.25, idl=0, idr=0, odl=0.05, odr=0.08,
+        ...        material=steel, rotary_inertia=True, shear_effects=True
+        ... )
+        >>> shaft1 # doctest: +ELLIPSIS
+        ShaftElement(L=0.25, idl=0.0...
+        """
+        return (
+            f"{self.__class__.__name__}"
+            f"(L={self.L:{0}.{5}}, idl={self.idl:{0}.{5}}, "
+            f"idr={self.idr:{0}.{5}}, odl={self.odl:{0}.{5}},  "
+            f"odr={self.odr:{0}.{5}}, material={self.material.name!r}, "
+            f"alpha={self.alpha:{0}.{5}}, beta={self.beta:{0}.{5}},  "
+            f"n={self.n})"
+        )
+
+    def __str__(self):
+        """
+        Method to convert object into string
+        
+        Parameters
+        ----------
+        Returns
+        -------
+        The object's parameters translated to strings
+        """
+        return (
+            f"\nElem. N:    {self.n}"
+            f"\nLenght:     {self.L:{10}.{5}}"
+            f"\nLeft Int. Diam.: {self.idl:{10}.{5}}"
+            f"\nLeft Out. Diam.: {self.odl:{10}.{5}}"
+            f"\nRight Int. Diam.: {self.idr:{10}.{5}}"
+            f"\nRight Out. Diam.: {self.odr:{10}.{5}}"
+            f"\nAlpha damp. factor: {self.alpha:{10}.{5}}"
+            f"\nBeta damp. factor: {self.beta:{10}.{5}}"
+            f'\n{35*"-"}'
+            f"\n{self.material}"
+            f"\n"
+        )
+
+    def M(self):
+        r"""Mass matrix for an instance of a 6 DoF shaft element.
+
+        Returns
+        -------
+        M: np.ndarray
+            Mass matrix for the 6 DoF shaft element.
+
+        Examples
+        --------
+        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
+        >>> Timoshenko_Element.M().shape
+        (12, 12)
+        """
+
+        # temporary material and geometrical constants
+        L = self.L
+        tempG = self.material.E / (2 * (1 + self.material.Poisson))
+        tempS = np.pi * (
+            ((self.odr / 2) ** 2 + (self.odl / 2) ** 2) / 2
+            - ((self.idr / 2) ** 2 + (self.idl / 2) ** 2) / 2
+        )
+        tempI = (
+            np.pi
+            / 4
+            * (
+                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
+                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
+            )
+        )
+        tempJ = (
+            np.pi
+            / 2
+            * (
+                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
+                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
+            )
+        )
+
+        # element level matrix declaration
+
+        aux1 = self.material.rho * tempS * L / 420
+        aux2 = self.material.rho * tempJ * L / 6
+        # fmt: off
+
+        # Standard mass matrix
+        M = aux1 * np.array([
+            [   156,     0, 0,      0,  -22*L, 0,    54,     0, 0,      0,   13*L, 0],
+            [     0,   156, 0,   22*L,      0, 0,     0,    54, 0,  -13*L,      0, 0],
+            [     0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
+            [     0,  22*L, 0, 4*L**2,      0, 0,     0,  13*L, 0,-3*L**2,      0, 0],
+            [ -22*L,     0, 0,      0, 4*L**2, 0, -13*L,     0, 0,      0,-3*L**2, 0],
+            [     0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
+            [    54,     0, 0,      0,  -13*L, 0,   156,     0, 0,      0,   22*L, 0],
+            [     0,    54, 0,   13*L,      0, 0,     0,   156, 0,  -22*L,      0, 0],
+            [     0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
+            [     0, -13*L, 0,-3*L**2,      0, 0,     0, -22*L, 0, 4*L**2,      0, 0],
+            [  13*L,     0, 0,      0,-3*L**2, 0,  22*L,     0, 0,      0, 4*L**2, 0],
+            [     0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
+        ])
+
+        # Secondary inertias mass matrix
+        Ms = self.material.rho * tempI / (30 * L) * np.array([
+            [   36,   0, 0,     0,  -3*L, 0, -36,    0, 0,     0,  -3*L, 0],
+            [    0,  36, 0,   3*L,     0, 0,   0,  -36, 0,   3*L,     0, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+            [    0, 3*L, 0,4*L**2,     0, 0,   0, -3*L, 0, -L**2,     0, 0],
+            [ -3*L,   0, 0,     0,4*L**2, 0, 3*L,    0, 0,     0, -L**2, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+            [  -36,   0, 0,     0,   3*L, 0,  36,    0, 0,     0,   3*L, 0],
+            [    0, -36, 0,  -3*L,     0, 0,   0,   36, 0,  -3*L,     0, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+            [    0, 3*L, 0, -L**2,     0, 0,   0, -3*L, 0,4*L**2,     0, 0],
+            [ -3*L,   0, 0,     0, -L**2, 0, 3*L,    0, 0,     0,4*L**2, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+        ])
+
+        # Axial terms inertia matrix
+        Ma = self.material.rho * tempS * L / 6 * np.array([
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        ])
+
+        # Torsional terms inertias matrix
+        Mr = self.material.rho * tempI * L / 6 * np.array([
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2],
+        ])
+
+        # fmt: on
+
+        M = M + Ms + Ma + Mr
+
+        return M
+
+    def K(self):
+        r"""Stiffness matrix for an instance of a 6 DoF shaft element.
+
+        Returns
+        -------
+        K: np.ndarray
+            Omega independent stiffness matrix for the 6 DoF shaft element.
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
+        >>> Timoshenko_Element.K().shape
+        (12, 12)
+        """
+
+        # Axial force and torque applied to the element.
+        Fa = self.axial_force
+        T = self.torque
+
+        # temporary material and geometrical constants, determined as mean values
+        # from the left and right radii of the taperad shaft
+        L = self.L
+        tempS = np.pi * (
+            ((self.odr / 2) ** 2 + (self.odl / 2) ** 2) / 2
+            - ((self.idr / 2) ** 2 + (self.idl / 2) ** 2) / 2
+        )
+        tempI = (
+            np.pi
+            / 4
+            * (
+                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
+                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
+            )
+        )
+        tempJ = (
+            np.pi
+            / 2
+            * (
+                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
+                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
+            )
+        )
+
+        # temporary variables
+        A = (
+            12
+            * self.material.E
+            * tempI
+            / (self.material.G_s * self.kappa * tempS * L ** 2)
+        )
+
+        # auxiliary variables
+        a1 = self.material.E * tempI / ((1 + A) * L ** 3)
+        a2 = self.material.G_s * tempJ / L
+        a3 = self.material.E * tempS / L
+
+        # fmt: off
+        # pure stiffness matrix [Kc], added to the axial loads stiffness matrix [Ka],
+        # torsional stiffnesses matrix [Kr] and Tinoshenko shear compensation [Ks].
+        Kc_plus = a1 * np.array([
+            [   12,   0,      0,            0,         -6*L,      0, -12,    0,      0,            0,         -6*L,      0],
+            [    0,  12,      0,          6*L,            0,      0,   0,  -12,      0,          6*L,            0,      0],
+            [    0,   0,  a3/a1,            0,            0,      0,   0,    0, -a3/a1,            0,            0,      0],
+            [    0, 6*L,      0, L**2*(A + 4),            0,      0,   0, -6*L,      0,-L**2*(A - 2),            0,      0],
+            [ -6*L,   0,      0,            0, L**2*(A + 4),      0, 6*L,    0,      0,            0,-L**2*(A - 2),      0],
+            [    0,   0,      0,            0,            0,  a2/a1,   0,    0,      0,            0,            0, -a2/a1],
+            [  -12,   0,      0,            0,          6*L,      0,  12,    0,      0,            0,          6*L,      0],
+            [    0, -12,      0,         -6*L,            0,      0,   0,   12,      0,         -6*L,            0,      0],
+            [    0,   0, -a3/a1,            0,            0,      0,   0,    0,  a3/a1,            0,            0,      0],
+            [    0, 6*L,      0,-L**2*(A - 2),            0,      0,   0, -6*L,      0, L**2*(A + 4),            0,      0],
+            [ -6*L,   0,      0,            0,-L**2*(A - 2),      0, 6*L,    0,      0,            0, L**2*(A + 4),      0],
+            [    0,   0,      0,            0,            0, -a2/a1,   0,    0,      0,            0,            0,  a2/a1],
+        ])
+
+        # stiffness matrix due to axial loading influence
+        Kf = Fa / (30 * L) * np.array([
+            [   36,   0, 0,     0,  -3*L, 0, -36,    0, 0,     0,  -3*L, 0],
+            [    0,  36, 0,   3*L,     0, 0,   0,  -36, 0,   3*L,     0, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+            [    0, 3*L, 0,4*L**2,     0, 0,   0, -3*L, 0, -L**2,     0, 0],
+            [ -3*L,   0, 0,     0,4*L**2, 0, 3*L,    0, 0,     0, -L**2, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+            [  -36,   0, 0,     0,   3*L, 0,  36,    0, 0,     0,   3*L, 0],
+            [    0, -36, 0,  -3*L,     0, 0,   0,   36, 0,  -3*L,     0, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+            [    0, 3*L, 0, -L**2,     0, 0,   0, -3*L, 0,4*L**2,     0, 0],
+            [ -3*L,   0, 0,     0, -L**2, 0, 3*L,    0, 0,     0,4*L**2, 0],
+            [    0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
+        ])
+
+        # stiffness matrix due to torque loading influence
+        Kt = T * np.array([
+            [    0,    0, 0, -1/L,    0, 0,    0,    0, 0,  1/L,    0, 0],
+            [    0,    0, 0,    0, -1/L, 0,    0,    0, 0,    0,  1/L, 0],
+            [    0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+            [ -1/L,    0, 0,    0,  1/2, 0,  1/L,    0, 0,    0,  1/2, 0],
+            [    0, -1/L, 0, -1/2,    0, 0,    0,  1/L, 0, -1/2,    0, 0],
+            [    0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+            [    0,    0, 0,  1/L,    0, 0,    0,    0, 0, -1/L,    0, 0],
+            [    0,    0, 0,    0,  1/L, 0,    0,    0, 0,    0, -1/L, 0],
+            [    0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+            [  1/L,    0, 0,    0, -1/2, 0, -1/L,    0, 0,    0, -1/2, 0],
+            [    0,  1/L, 0,  1/2,    0, 0,    0, -1/L, 0,  1/2,    0, 0],
+            [    0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+        ])
+
+        # fmt: on
+        # Dynamic stiffness matrix is added independently in "def Kst"
+        # Kst = self.material.rho*tempI/(15*L)*np.array(12,12)
+
+        K = Kc_plus + Kf + Kt
+
+        return K
+
+    def Kst(self):
+        r"""Dynamic stiffness matrix for an instance of a 6 DoF shaft element.
+
+        Returns
+        -------
+        Kst: np.ndarray
+            Dynamic stiffness matrix for the 6 DoF shaft element. This is
+            directly dependent on the rotation speed Omega. It needs to be
+            multiplied by the adequate Omega value when used in time depen-
+            dent analyses. The matrix multiplier term is: 
+            
+            [(Iz*Omega*rho)/(15*L)] * [Kst]
+
+            and here the Omega value has been suppressed and must be added
+            in the adequate analyses.
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
+        >>> Timoshenko_Element.Kst().shape
+        (12, 12)
+        """
+        # temporary material and geometrical constants, determined as mean values
+        # from the left and right radii of the taperad shaft
+        L = self.L
+        tempI = (
+            np.pi
+            / 4
+            * (
+                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
+                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
+            )
+        )
+
+        # fmt: off
+        # dynamic stiffening matrix
+        Kst = self.material.rho * tempI / (15 * L) * np.array([
+            [ 0, -36, 0,  -3*L, 0, 0, 0,   36, 0,  -3*L, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0, 3*L, 0,4*L**2, 0, 0, 0, -3*L, 0, -L**2, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0,  36, 0,   3*L, 0, 0, 0,  -36, 0,   3*L, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+            [ 0, 3*L, 0, -L**2, 0, 0, 0, -3*L, 0,4*L**2, 0, 0],
+            [ 0,   0, 0,     0, 0, 0, 0,    0, 0,     0, 0, 0],
+        ])
+        # fmt: on
+
+        return Kst
+
+    def C(self):
+        r"""Proportional damping matrix for an instance of a 6 DoF shaft element.
+
+        Returns
+        -------
+        C: np.ndarray
+            Proportional damping matrix for the 6 DoF shaft element.
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> Timoshenko_Element = ShaftElement6DoF(L=0.25, idl=0, odl=0.05, material=steel)
+        >>> Timoshenko_Element.C().shape
+        (12, 12)
+        """
+
+        # proportional damping matrix
+        C = self.alpha * self.M() + self.beta * self.K()
+
+        return C
+
+    def G(self):
+        r"""Gyroscopic matrix for an instance of a 6 DoFs shaft element.
+
+        Returns
+        -------
+        G: np.ndarray
+            Gyroscopic matrix for the 6 DoF shaft element. Similar to the Kst
+            stiffness matrix, this Gyro matrix is also multiplied by the value
+            of the rotating speed Omega. It is omitted from this and must be 
+            added in the respective analyses.
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
+        >>> Timoshenko_Element.G().shape
+        (12, 12)
+        """
+
+        if self.gyroscopic:
+            # temporary material and geometrical constants, determined as mean values
+            # from the left and right radii of the tapered shaft
+            L = self.L
+            tempI = (
+                np.pi
+                / 4
+                * (
+                    ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
+                    - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
+                )
+            )
+
+            # fmt: off
+            # Gyroscopic effect matrix
+            G = (self.material.rho * tempI / (15 * L)) * np.array([
+                [   0, -36, 0,  -3*L,      0, 0,    0,   36, 0,  -3*L,      0, 0],
+                [  36,   0, 0,     0,   -3*L, 0,  -36,    0, 0,     0,   -3*L, 0],
+                [   0,   0, 0,     0,      0, 0,    0,    0, 0,     0,      0, 0],
+                [ 3*L,   0, 0,     0,-4*L**2, 0, -3*L,    0, 0,     0,   L**2, 0],
+                [   0, 3*L, 0,4*L**2,      0, 0,    0, -3*L, 0, -L**2,      0, 0],
+                [   0,   0, 0,     0,      0, 0,    0,    0, 0,     0,      0, 0],
+                [   0,  36, 0,   3*L,      0, 0,    0,  -36, 0,   3*L,      0, 0],
+                [ -36,   0, 0,     0,    3*L, 0,   36,    0, 0,     0,    3*L, 0],
+                [   0,   0, 0,     0,      0, 0,    0,    0, 0,     0,      0, 0],
+                [ 3*L,   0, 0,     0,   L**2, 0, -3*L,    0, 0,     0,-4*L**2, 0],
+                [   0, 3*L, 0, -L**2,      0, 0,    0, -3*L, 0,4*L**2,      0, 0],
+                [   0,   0, 0,     0,      0, 0,    0,    0, 0,     0,      0, 0],
+            ])
+            # fmt: on
+
+        return G
+
+    def bokeh_patch(self, position, check_sld, bk_ax):
+        """Shaft element patch.
+        Patch that will be used to draw the shaft element.
+        Parameters
+        ----------
+        position : float
+            Position in which the patch will be drawn.
+        check_sld : bool
+            If True, HoverTool displays only the slenderness ratio and color
+            the elements in yellow if slenderness ratio < 1.6
+        bk_ax : bokeh plotting axes, optional
+            Axes in which the plot will be drawn.
+        Returns
+        -------
+        hover : Bokeh HoverTool
+            Bokeh HoverTool axes
+        """
+        if check_sld is True and self.slenderness_ratio < 1.6:
+            bk_color = "yellow"
+            legend = "Shaft - Slenderness Ratio < 1.6"
+        else:
+            bk_color = self.material.color
+            legend = "Shaft"
+
+        # bokeh plot - plot the shaft
+        z_upper = [position, position, position + self.L, position + self.L]
+        y_upper = [self.idl / 2, self.odl / 2, self.odr / 2, self.idr / 2]
+        z_lower = [position, position, position + self.L, position + self.L]
+        y_lower = [-self.idl / 2, -self.odl / 2, -self.odr / 2, -self.idr / 2]
+
+        source = ColumnDataSource(
+            dict(
+                z_l=[z_lower],
+                y_l=[y_lower],
+                z_u=[z_upper],
+                y_u=[y_upper],
+                sld=[self.slenderness_ratio],
+                elnum=[self.n],
+                out_d_l=[self.odl],
+                out_d_r=[self.odr],
+                in_d_l=[self.idl],
+                in_d_r=[self.idr],
+                alpha_factor=[self.alpha],
+                beta_factor=[self.beta],
+                length=[self.L],
+                mat=[self.material.name],
+            )
+        )
+
+        bk_ax.patches(
+            xs="z_u",
+            ys="y_u",
+            source=source,
+            line_color=bokeh_colors[0],
+            line_width=1,
+            fill_alpha=0.5,
+            fill_color=bk_color,
+            legend_label=legend,
+            name="u_shaft",
+        )
+        bk_ax.patches(
+            xs="z_l",
+            ys="y_l",
+            source=source,
+            line_color=bokeh_colors[0],
+            line_width=1,
+            fill_alpha=0.5,
+            fill_color=bk_color,
+            legend_label=legend,
+            name="l_shaft",
+        )
+
+        hover = HoverTool(names=["l_shaft", "u_shaft"])
+        if check_sld:
+            hover.tooltips = [
+                ("Element Number :", "@elnum"),
+                ("Slenderness Ratio :", "@sld"),
+            ]
+        else:
+            hover.tooltips = [
+                ("Element Number :", "@elnum"),
+                ("Left Outer Diameter :", "@out_d_l"),
+                ("Left Inner Diameter :", "@in_d_l"),
+                ("Right Outer Diameter :", "@out_d_r"),
+                ("Right Inner Diameter :", "@in_d_r"),
+                ("Alpha Damp. Factor :", "@alpha_factor"),
+                ("Beta Damp. Factor :", "@beta_factor"),
+                ("Element Length :", "@length"),
+                ("Material :", "@mat"),
+            ]
+        hover.mode = "mouse"
+
+        return hover
+
+    @classmethod
+    def section(
+        cls,
+        L,
+        ne,
+        s_idl,
+        s_odl,
+        s_idr=None,
+        s_odr=None,
+        alpha=0,
+        beta=0,
+        material=None,
+        n=None,
+        gyroscopic=True,
+    ):
+        """Shaft section constructor.
+
+        This method will create a shaft section with length 'L' divided into
+        'ne' elements.
+
+        Parameters
+        ----------
+        odl : float, pint.Quantity
+            Outer diameter of the element at the left node.
+        idr : float, pint.Quantity, optional
+            Inner diameter of the element at the right node;
+            Default is equal to idl value for cylindrical element.
+        odr : float, pint.Quantity, optional
+            Outer diameter of the element at the right node;
+            Default is equal to odl value for cylindrical element.
+        E : float
+            Young's modulus.
+        G_s : float
+            Shear modulus.
+        alpha : float
+            Proportional damping coefficient, associated to the element Mass matrix
+        beta : float
+            Proportional damping coefficient, associated to the element Stiffness matrix
+        material : ross.material
+            Shaft material.
+        n : int, optional
+            Element number (coincident with it's first node).
+            If not given, it will be set when the rotor is assembled
+            according to the element's position in the list supplied to
+            the rotor constructor.
+        axial_force : float
+            Axial force.
+        torque : float
+            Torque.
+        gyroscopic : bool
+            Determine if gyroscopic effects are taken into account.
+            Default is False.
+
+        Returns
+        -------
+        elements: list
+            List with the 'ne' shaft elements.
+
+        Examples
+        --------
+        >>> # shaft material
+        >>> from ross.materials import steel
+        >>> # shaft inner and outer diameters
+        >>> s_idl = 0
+        >>> s_odl = 0.01585
+        >>> sec = ShaftElement.section(247.65e-3, 4, 0, 15.8e-3, material=steel)
+        >>> len(sec)
+        4
+        >>> sec[0].i_d
+        0.0
+        """
+        if s_idr is None:
+            s_idr = s_idl
+        if s_odr is None:
+            s_odr = s_odl
+
+        le = L / ne
+
+        elements = [
+            cls(
+                le,
+                (s_idr - s_idl) * i * le / L + s_idl,
+                (s_odr - s_odl) * i * le / L + s_odl,
+                (s_idr - s_idl) * (i + 1) * le / L + s_idl,
+                (s_odr - s_odl) * (i + 1) * le / L + s_odl,
+                alpha,
+                beta,
+                material,
+                n,
+                gyroscopic,
+            )
+            for i in range(ne)
+        ]
+
+        return elements


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index ccc6e6e4..f3e684fe 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -6,6 +6,7 @@ import pytest
 from numpy.testing import assert_allclose
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
+                                       BearingElement6DoF,
                                        RollerBearingElement)
 
 # fmt: on
@@ -287,3 +288,42 @@ def test_roller_bearing_element():
     assert_allclose(rollerbearing.K(0), K)
     assert_allclose(rollerbearing.C(0), C)
     assert_allclose(rollerbearing.G(), G)
+
+
+def test_bearing_6dof():
+    bearing_6dof = BearingElement6DoF(
+        n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
+    )
+    # fmt: off
+    K = np.array(
+        [[1000000.,       0.,       0.],
+         [      0.,  800000.,       0.],
+         [      0.,       0.,  100000.]])
+    C = np.array(
+        [[200.,   0.,   0.],
+         [  0., 150.,   0.],
+         [  0.,   0.,  50.]])
+    M = np.zeros((3, 3))
+    G = np.zeros((3, 3))
+    # fmt: on
+
+    assert_allclose(bearing_6dof.K(0), K, rtol=1e-3)
+    assert_allclose(bearing_6dof.C(0), C, rtol=1e-3)
+    assert_allclose(bearing_6dof.M(), M, rtol=1e-3)
+    assert_allclose(bearing_6dof.G(), G, rtol=1e-3)
+
+
+def test_bearing_6dof_equality():
+    bearing_6dof_0 = BearingElement6DoF(
+        n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
+    )
+    bearing_6dof_1 = BearingElement6DoF(
+        n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
+    )
+    bearing_6dof_2 = BearingElement6DoF(
+        n=0, kxx=2e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
+    )
+
+    assert bearing_6dof_0 == bearing_6dof_1
+    assert not bearing_6dof_1 == bearing_6dof_2
+    assert not bearing_6dof_0 == bearing_6dof_2
diff --git a/ross/tests/test_disk_element.py b/ross/tests/test_disk_element.py
index c1dcc4a1..c46a57aa 100644
--- a/ross/tests/test_disk_element.py
+++ b/ross/tests/test_disk_element.py
@@ -4,7 +4,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_allclose, assert_almost_equal
 
-from ross.disk_element import DiskElement
+from ross.disk_element import DiskElement, DiskElement6DoF
 from ross.materials import steel
 
 
@@ -75,3 +75,85 @@ def test_from_table():
         assert_allclose(disks[1].m, 6.90999178227835)
         assert_allclose(disks[1].Ip, 0.0469996988106328, atol=1.6e-07)
         assert_allclose(disks[1].Id, 0.0249998397928898, atol=1.6e-07)
+
+
+@pytest.fixture
+def disk_6dof():
+    return DiskElement6DoF(0, 32.58973, 0.17809, 0.32956)
+
+
+def test_index_6dof(disk_6dof):
+    assert disk_6dof.dof_local_index().u_0 == 0
+    assert disk_6dof.dof_local_index().v_0 == 1
+    assert disk_6dof.dof_local_index().w_0 == 2
+    assert disk_6dof.dof_local_index().theta_0 == 3
+    assert disk_6dof.dof_local_index().psi_0 == 4
+    assert disk_6dof.dof_local_index().phi_0 == 5
+
+
+def test_mass_matrix_disk_6dof(disk_6dof):
+    # fmt: off
+    Md1 = np.array([[32.58973, 0.,       0.,       0.,      0.,      0.   ],
+                    [0.,       32.58973, 0.,       0.,      0.,      0.   ],
+                    [0.,       0.,       32.58973, 0.,      0.,      0.   ],
+                    [0.,       0.,       0.,       0.17809, 0.,      0.   ],
+                    [0.,       0.,       0.,       0.,      0.17809, 0.   ],
+                    [0.,       0.,       0.,       0,       0.,      0.32956]])
+    # fmt: on
+    assert_almost_equal(disk_6dof.M(), Md1, decimal=5)
+
+
+def test_gyroscopic_matrix_disk_6dof(disk_6dof):
+    # fmt: off
+    Gd1 = np.array([
+            [ 0, 0, 0,  0,         0,       0],
+            [ 0, 0, 0,  0,         0,       0],
+            [ 0, 0, 0,  0,         0,       0],
+            [ 0, 0, 0,  0,         0.32956, 0],
+            [ 0, 0, 0, -0.32956,   0,       0],
+            [ 0, 0, 0,  0,         0,       0]])
+    # fmt: on
+
+    assert_almost_equal(disk_6dof.G(), Gd1, decimal=5)
+
+
+@pytest.fixture
+def disk_from_geometry_6dof():
+    return DiskElement6DoF.from_geometry(0, steel, 0.07, 0.05, 0.28)
+
+
+def test_mass_matrix_disk1_6dof(disk_from_geometry_6dof):
+    # fmt: off
+    Md1 = np.array([[32.58973, 0.,       0.,       0.,      0.,      0.   ],
+                    [0.,       32.58973, 0.,       0.,      0.,      0.   ],
+                    [0.,       0.,       32.58973, 0.,      0.,      0.   ],
+                    [0.,       0.,       0.,       0.17809, 0.,      0.   ],
+                    [0.,       0.,       0.,       0.,      0.17809, 0.   ],
+                    [0.,       0.,       0.,       0,       0.,      0.32956]])
+    # fmt: on
+
+    M_class = disk_from_geometry_6dof.M()
+    assert_almost_equal(M_class, Md1, decimal=5)
+
+
+def test_gyroscopic_matrix_disk1_6dof(disk_from_geometry_6dof):
+    # fmt: off
+    Gd1 = np.array([
+            [ 0, 0, 0,  0,         0,       0],
+            [ 0, 0, 0,  0,         0,       0],
+            [ 0, 0, 0,  0,         0,       0],
+            [ 0, 0, 0,  0,         0.32956, 0],
+            [ 0, 0, 0, -0.32956,   0,       0],
+            [ 0, 0, 0,  0,         0,       0]])
+    # fmt: on
+    assert_almost_equal(disk_from_geometry_6dof.G(), Gd1, decimal=5)
+
+
+def test_from_table_6dof():
+    for file_name in ["/data/shaft_us.xls", "/data/shaft_si.xls"]:
+        file_name = os.path.dirname(os.path.realpath(__file__)) + file_name
+        disks = DiskElement6DoF.from_table(file_name, sheet_name="More")
+
+        assert_allclose(disks[1].m, 6.90999178227835)
+        assert_allclose(disks[1].Ip, 0.0469996988106328, atol=1.6e-07)
+        assert_allclose(disks[1].Id, 0.0249998397928898, atol=1.6e-07)
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index 75c611a0..22392e23 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -5,7 +5,7 @@ import pytest
 from numpy.testing import assert_allclose, assert_almost_equal
 
 from ross.materials import steel
-from ross.shaft_element import ShaftElement
+from ross.shaft_element import ShaftElement, ShaftElement6DoF
 
 test_dir = os.path.dirname(__file__)
 
@@ -306,3 +306,121 @@ def test_match_gyroscopic_matrix(tap2, tim2):
     G_tap = tap2.G()
     G_tim = tim2.G()
     assert_almost_equal(G_tap, G_tim, decimal=5)
+
+
+@pytest.fixture
+def s6_eb():
+
+    #  Euler-Bernoulli element
+    """
+    L = 0.1
+    odr = 25.4 * 3 / 4
+    odl = 25.4 * 3 / 4
+    idr = 25.4 * 1 / 2
+    idl = 25.4 * 1 / 2
+    axial_force = 1
+    torque = 1
+    alpha = 2.7
+    beta = 4.8 * 10 ** (-6)
+    """
+
+    return ShaftElement6DoF(
+        L=0.1,
+        idl=25.4 * 1 / 2,
+        odl=25.4 * 3 / 4,
+        idr=25.4 * 1 / 2,
+        odr=25.4 * 3 / 4,
+        material=steel,
+        axial_force=1,
+        torque=1,
+        alpha=2.7,
+        beta=4.8 * 10 ** (-6),
+    )
+
+
+def test_s6_parameters_eb(s6_eb):
+    assert s6_eb.L == 0.1
+    assert s6_eb.idl == 25.4 * 1 / 2
+    assert s6_eb.odl == 25.4 * 3 / 4
+    assert s6_eb.idr == 25.4 * 1 / 2
+    assert s6_eb.odr == 25.4 * 3 / 4
+    assert s6_eb.material.E == 211e9
+    assert s6_eb.material.G_s == 81.2e9
+    assert s6_eb.material.rho == 7810
+
+
+def test_s6_M_matrix(s6_eb):
+    # fmt: off
+    M_mat_6DoF = np.array(
+        [[  4.8624e+08,            0,            0,            0,  -4.0523e+06,            0,  -4.8618e+08,            0,            0,            0,  -4.0512e+06,            0],
+         [           0,   4.8624e+08,            0,   4.0523e+06,            0,            0,            0,  -4.8618e+08,            0,   4.0512e+06,            0,            0],
+         [           0,            0,   4.1223e+04,            0,            0,            0,            0,            0,   2.0611e+04,            0,            0,            0],
+         [           0,   4.0523e+06,            0,   5.4023e+05,            0,            0,            0,  -4.0512e+06,            0,  -1.3506e+05,            0,            0],
+         [ -4.0523e+06,            0,            0,            0,   5.4023e+05,            0,   4.0512e+06,            0,            0,            0,  -1.3506e+05,            0],
+         [           0,            0,            0,            0,            0,   1.3505e+06,            0,            0,            0,            0,            0,   6.7527e+05],
+         [ -4.8618e+08,            0,            0,            0,   4.0512e+06,            0,   4.8624e+08,            0,            0,            0,   4.0523e+06,            0],
+         [           0,  -4.8618e+08,            0,  -4.0512e+06,            0,            0,            0,   4.8624e+08,            0,  -4.0523e+06,            0,            0],
+         [           0,            0,   2.0611e+04,            0,            0,            0,            0,            0,   4.1223e+04,            0,            0,            0],
+         [           0,   4.0512e+06,            0,  -1.3506e+05,            0,            0,            0,  -4.0523e+06,            0,   5.4023e+05,            0,            0],
+         [ -4.0512e+06,            0,            0,            0,  -1.3506e+05,            0,   4.0523e+06,            0,            0,            0,   5.4023e+05,            0],
+         [           0,            0,            0,            0,            0,   6.7527e+05,            0,            0,            0,            0,            0,   1.3505e+06]])
+    # fmt: on
+    assert_allclose(s6_eb.M(), M_mat_6DoF, rtol=1e-3)
+
+
+def test_s6_G_matrix(s6_eb):
+    # fmt: off
+    G_mat_6DoF = 0.01 * np.array(
+        [[           0,  -9.7239e+10,            0,  -8.1032e+08,            0,            0,            0,   9.7239e+10,            0,  -8.1032e+08,            0,            0],
+         [  9.7239e+10,            0,            0,            0,  -8.1032e+08,            0,  -9.7239e+10,            0,            0,            0,  -8.1032e+08,            0],
+         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+         [  8.1032e+08,            0,            0,            0,  -1.0804e+08,            0,  -8.1032e+08,            0,            0,            0,   2.7011e+07,            0],
+         [           0,   8.1032e+08,            0,   1.0804e+08,            0,            0,            0,  -8.1032e+08,            0,  -2.7011e+07,            0,            0],
+         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+         [           0,   9.7239e+10,            0,   8.1032e+08,            0,            0,            0,  -9.7239e+10,            0,   8.1032e+08,            0,            0],
+         [ -9.7239e+10,            0,            0,            0,   8.1032e+08,            0,   9.7239e+10,            0,            0,            0,   8.1032e+08,            0],
+         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+         [  8.1032e+08,            0,            0,            0,   2.7011e+07,            0,  -8.1032e+08,            0,            0,            0,  -1.0804e+08,            0],
+         [           0,   8.1032e+08,            0,  -2.7011e+07,            0,            0,            0,  -8.1032e+08,            0,   1.0804e+08,            0,            0],
+         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0]])
+    # fmt: on
+    assert_allclose(s6_eb.G(), G_mat_6DoF, rtol=1e-3)
+
+
+def test_s6_K_matrix(s6_eb):
+    # fmt: off
+    K_mat_6DoF = np.array(
+        [[ 7.7093e+13,            0,            0,  -1.0000e+01,  -3.8547e+12,            0,  -7.7093e+13,            0,            0,   1.0000e+01,  -3.8547e+12,            0],
+         [          0,   7.7093e+13,            0,   3.8547e+12,  -1.0000e+01,            0,            0,  -7.7093e+13,            0,   3.8547e+12,   1.0000e+01,            0],
+         [          0,            0,   3.3411e+14,            0,            0,            0,            0,            0,  -3.3411e+14,            0,            0,            0],
+         [-1.0000e+01,   3.8547e+12,            0,   1.0946e+16,   5.0000e-01,            0,   1.0000e+01,  -3.8547e+12,            0,  -1.0946e+16,   5.0000e-01,            0],
+         [-3.8547e+12,  -1.0000e+01,            0,  -5.0000e-01,   1.0946e+16,            0,   3.8547e+12,   1.0000e+01,            0,  -5.0000e-01,  -1.0946e+16,            0],
+         [          0,            0,            0,            0,            0,   8.4249e+15,            0,            0,            0,            0,            0,  -8.4249e+15],
+         [-7.7093e+13,            0,            0,   1.0000e+01,   3.8547e+12,            0,   7.7093e+13,            0,            0,  -1.0000e+01,   3.8547e+12,            0],
+         [          0,  -7.7093e+13,            0,  -3.8547e+12,   1.0000e+01,            0,            0,   7.7093e+13,            0,  -3.8547e+12,  -1.0000e+01,            0],
+         [          0,            0,  -3.3411e+14,            0,            0,            0,            0,            0,   3.3411e+14,            0,            0,            0],
+         [ 1.0000e+01,   3.8547e+12,            0,  -1.0946e+16,  -5.0000e-01,            0,  -1.0000e+01,  -3.8547e+12,            0,   1.0946e+16,  -5.0000e-01,            0],
+         [-3.8547e+12,   1.0000e+01,            0,   5.0000e-01,  -1.0946e+16,            0,   3.8547e+12,  -1.0000e+01,            0,   5.0000e-01,   1.0946e+16,            0],
+         [          0,            0,            0,            0,            0,  -8.4249e+15,            0,            0,            0,            0,            0,   8.4249e+15]])
+
+    # fmt: on
+    assert_allclose(s6_eb.K(), K_mat_6DoF, rtol=1e-3)
+
+
+def test_s6_Kst_matrix(s6_eb):
+    # fmt: off
+    Kst_mat_6DoF = 0.01 * np.array([
+        [          0,  -9.7239e+10,            0,  -8.1032e+08,            0,            0,            0,   9.7239e+10,            0,  -8.1032e+08,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,   8.1032e+08,            0,   1.0804e+08,            0,            0,            0,  -8.1032e+08,            0,  -2.7011e+07,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,   9.7239e+10,            0,   8.1032e+08,            0,            0,            0,  -9.7239e+10,            0,   8.1032e+08,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
+        [          0,   8.1032e+08,            0,  -2.7011e+07,            0,            0,            0,  -8.1032e+08,            0,   1.0804e+08,            0,            0],
+        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0]])
+    # fmt: on
+    assert_allclose(s6_eb.Kst(), Kst_mat_6DoF, rtol=1e-3)