Repository: ross-rotordynamics/ross
Pull request #528 (2020-04-20 13:02:05)

Plotting statistics for random variables
- Add `.plot_random_var()`: a method to plot histograms with the Frequency or the PDF for random variables.
- Add tests for `.plot_random_var()`

Example:
```
kxx = np.random.uniform(1e6, 2e6, 1000)
cxx = 1e3
elm = ST_BearingElement(n=1, kxx=kxx, cxx=cxx, is_random=["kxx"])
show(elm.plot_random_var(["kxx"]))
show(elm.plot_random_var(["kxx"], density=False))
```
The pdf estimation use the "kernel density estimation" tool from Scipy.

![image](https://user-images.githubusercontent.com/45969994/79754083-539ed680-82ed-11ea-83ae-64927606109b.png)

![image](https://user-images.githubusercontent.com/45969994/79754259-a082ad00-82ed-11ea-8ebc-abea807f1419.png)



Patch:
diff --git a/.gitignore b/.gitignore
index 4e8193c7..81034dd3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -111,23 +111,8 @@ 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/download.rst b/docs/download.rst
index 0e46bf69..a31b1dd5 100644
--- a/docs/download.rst
+++ b/docs/download.rst
@@ -12,4 +12,3 @@ Here you can download the tutorial and examples notebooks.
 | :download:`Example 6 </examples/example_6.ipynb>`
 | :download:`Example 7 </examples/example_7.ipynb>`
 | :download:`Example 8 </examples/example_8.ipynb>`
-| :download:`Stochastic ROSS - Tutorial </examples/example_stochastic_ross.ipynb>`
diff --git a/docs/examples.rst b/docs/examples.rst
index 642ab9c2..717b8e93 100644
--- a/docs/examples.rst
+++ b/docs/examples.rst
@@ -11,5 +11,4 @@ Examples
    examples/example_5
    examples/example_6
    examples/example_7
-   examples/example_8
-   examples/example_stochastic_ross
\ No newline at end of file
+   examples/example_8
\ No newline at end of file
diff --git a/docs/examples/example_stochastic_ross.ipynb b/docs/examples/tutorial_stochastic_ross.ipynb
similarity index 100%
rename from docs/examples/example_stochastic_ross.ipynb
rename to docs/examples/tutorial_stochastic_ross.ipynb
diff --git a/docs/refs.bib b/docs/refs.bib
index 1290c2b4..6d2cd068 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -27,13 +27,6 @@
   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},
@@ -634,10 +627,3 @@ 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 95fdafc1..f1c7b782 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,3 @@
-
 numpy
 scipy
 matplotlib
@@ -6,4 +5,4 @@ toml
 pandas>=0.23
 bokeh
 xlrd
-pint
+pint
\ No newline at end of file
diff --git a/ross/__init__.py b/ross/__init__.py
index 3093ba63..97d27422 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -1,4 +1,4 @@
-__version__ = "0.3.3"
+__version__ = "0.3.2"
 
 from .api_report import *
 from .bearing_seal_element import *
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 33524cac..25981fb4 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -21,7 +21,6 @@ __all__ = [
     "SealElement",
     "BallBearingElement",
     "RollerBearingElement",
-    "BearingElement6DoF",
     "MagneticBearingElement",
 ]
 # fmt: on
@@ -419,7 +418,7 @@ class BearingElement(Element):
         return dict(x_0=0, y_0=1)
 
     def dof_global_index(self):
-        """Get the global index for an element specific degree of freedom."""
+        """Get the global index for a element specific degree of freedom."""
         global_index = super().dof_global_index()
 
         if self.n_link is not None:
@@ -1454,357 +1453,6 @@ 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
@@ -1857,33 +1505,3 @@ 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 0363e8df..eb869324 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", "DiskElement6DoF"]
+__all__ = ["DiskElement"]
 bokeh_colors = bp.RdGy[11]
 
 
@@ -557,227 +557,3 @@ 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 5eec886b..d305f3ad 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", "ShaftElement6DoF"]
+__all__ = ["ShaftElement"]
 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,953 +1191,3 @@ 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
diff --git a/ross/stochastic/st_bearing_seal_element.py b/ross/stochastic/st_bearing_seal_element.py
index e4130d9d..e9bbc115 100644
--- a/ross/stochastic/st_bearing_seal_element.py
+++ b/ross/stochastic/st_bearing_seal_element.py
@@ -9,10 +9,11 @@ from ross.bearing_seal_element import BearingElement
 from ross.fluid_flow import fluid_flow as flow
 from ross.fluid_flow.fluid_flow_coefficients import (
     calculate_damping_matrix, calculate_stiffness_matrix)
+from ross.stochastic.st_results_elements import plot_histogram
 
 # fmt: on
 
-__all__ = ["ST_BearingElement"]
+__all__ = ["ST_BearingElement", "st_bearing_example"]
 
 
 class ST_BearingElement:
@@ -88,7 +89,7 @@ class ST_BearingElement:
     ...                              cxx=cxx,
     ...                              is_random = ["kxx", "cxx"],
     ...                              )
-    >>> len(list(elms.__iter__()))
+    >>> len(list(iter(elms)))
     10
 
     # Uncertanties on bearing coefficients varying with frequency
@@ -105,7 +106,7 @@ class ST_BearingElement:
     ...                              frequency=frequency,
     ...                              is_random = ["kxx", "cxx"],
     ...                              )
-    >>> len(list(elms.__iter__()))
+    >>> len(list(iter(elms)))
     5
     """
 
@@ -163,6 +164,13 @@ class ST_BearingElement:
         Returns
         -------
         An iterator over random bearing elements.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> bearing = srs.st_bearing_example()
+        >>> len(list(iter(bearing)))
+        2
         """
         return iter(self.random_var(self.is_random, self.attribute_dict))
 
@@ -214,7 +222,7 @@ class ST_BearingElement:
         key : str
             A class parameter as string.
         value : The corresponding value for the attrbiute_dict's key.
-            ***check the correct type for each key in ST_ShaftElement
+            ***check the correct type for each key in ST_BearingElement
             docstring.
 
         Raises
@@ -246,14 +254,14 @@ class ST_BearingElement:
         """Generate a list of objects as random attributes.
 
         This function creates a list of objects with random values for selected
-        attributes from BearingElement.
+        attributes from ross.BearingElement.
 
         Parameters
         ----------
         is_random : list
             List of the object attributes to become stochastic.
         *args : dict
-            Dictionary instanciating the ShaftElement class.
+            Dictionary instanciating the ross.BearingElement class.
             The attributes that are supposed to be stochastic should be
             set as lists of random variables.
 
@@ -289,6 +297,51 @@ class ST_BearingElement:
 
         return f_list
 
+    def plot_random_var(self, var_list=[], **kwargs):
+        """Plot histogram and the PDF.
+
+        This function creates a histogram to display the random variable
+        distribution.
+
+        Parameters
+        ----------
+        var_list : list, optional
+            List of random variables, in string format, to plot.
+        **kwargs : optional
+            Additional key word arguments can be passed to change
+            the numpy.histogram (e.g. density=True, bins=11, ...)
+
+        Returns
+        -------
+        grid_plot : bokeh row
+            A row with the histogram plots.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_bearing_example()
+        >>> elm.plot_random_var(["kxx"]) # doctest: +ELLIPSIS
+        Row...
+        """
+        label = dict(
+            kxx="Direct stiffness in the X direction",
+            kxy="Cross coupled stiffness in the X direction",
+            kyx="Cross coupled stiffness in the Y direction",
+            kyy="Direct stiffness in the Y direction",
+            cxx="Direct damping in the X direction",
+            cxy="Cross coupled damping in the X direction",
+            cyx="Cross coupled damping in the Y direction",
+            cyy="Direct damping in the y direction",
+        )
+        if not all(var in self.is_random for var in var_list):
+            raise ValueError(
+                "Not random variable in var_list. Select variables from {}".format(
+                    self.is_random
+                )
+            )
+
+        return plot_histogram(self.attribute_dict, label, var_list, **kwargs)
+
     @classmethod
     def from_fluid_flow(
         cls,
@@ -414,7 +467,7 @@ class ST_BearingElement:
         ...     radius_stator, visc, rho,
         ...     eccentricity=eccentricity, is_random=["visc"]
         ... )
-        >>> len(list(elms.__iter__()))
+        >>> len(list(iter(elms)))
         5
         """
         attribute_dict = locals()
@@ -479,3 +532,28 @@ class ST_BearingElement:
             scale_factor=scale_factor,
             is_random=list(args_dict.keys()),
         )
+
+
+def st_bearing_example():
+    """Return an instance of a simple random bearing.
+
+    The purpose is to make available a simple model so that doctest can be
+    written using it.
+
+    Returns
+    -------
+    elm : ross.stochastic.ST_BearingElement
+        An instance of a random bearing element object.
+
+    Examples
+    --------
+    >>> import ross.stochastic as srs
+    >>> elm = srs.st_bearing_example()
+    >>> len(list(iter(elm)))
+    2
+    """
+    kxx = [1e6, 2e6]
+    cxx = [1e3, 2e3]
+    elm = ST_BearingElement(n=1, kxx=kxx, cxx=cxx, is_random=["kxx", "cxx"])
+
+    return elm
diff --git a/ross/stochastic/st_disk_element.py b/ross/stochastic/st_disk_element.py
index 1fe31d68..e3b9f9cd 100644
--- a/ross/stochastic/st_disk_element.py
+++ b/ross/stochastic/st_disk_element.py
@@ -1,12 +1,18 @@
+"""Disk element module for STOCHASTIC ROSS.
+
+This module creates an instance of random disk element for stochastic
+analysis.
+"""
 import bokeh.palettes as bp
 import numpy as np
 
 from ross.disk_element import DiskElement
 from ross.stochastic.st_materials import ST_Material
+from ross.stochastic.st_results_elements import plot_histogram
 
 bokeh_colors = bp.RdGy[11]
 
-__all__ = ["ST_DiskElement"]
+__all__ = ["ST_DiskElement", "st_disk_example"]
 
 
 class ST_DiskElement:
@@ -48,7 +54,7 @@ class ST_DiskElement:
     ...                           Ip=np.random.uniform(0.15, 0.25, 5),
     ...                           is_random=["Id", "Ip"],
     ...                           )
-    >>> len(list(elms.__iter__()))
+    >>> len(list(iter(elms)))
     5
     """
 
@@ -66,6 +72,13 @@ class ST_DiskElement:
         Returns
         -------
         An iterator over random disk elements.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_disk_example()
+        >>> len(list(iter(elm)))
+        2
         """
         return iter(self.random_var(self.is_random, self.attribute_dict))
 
@@ -115,7 +128,7 @@ class ST_DiskElement:
         key : str
             A class parameter as string.
         value : The corresponding value for the attrbiute_dict's key.
-            ***check the correct type for each key in ST_ShaftElement
+            ***check the correct type for each key in ST_DiskElement
             docstring.
 
         Raises
@@ -145,14 +158,14 @@ class ST_DiskElement:
         """Generate a list of objects as random attributes.
 
         This function creates a list of objects with random values for selected
-        attributes from DiskElement.
+        attributes from ross.DiskElement.
 
         Parameters
         ----------
         is_random : list
             List of the object attributes to become stochastic.
         *args : dict
-            Dictionary instanciating the ShaftElement class.
+            Dictionary instanciating the ross.DiskElement class.
             The attributes that are supposed to be stochastic should be
             set as lists of random variables.
 
@@ -175,6 +188,44 @@ class ST_DiskElement:
 
         return f_list
 
+    def plot_random_var(self, var_list=[], **kwargs):
+        """Plot histogram and the PDF.
+
+        This function creates a histogram to display the random variable
+        distribution.
+
+        Parameters
+        ----------
+        var_list : list, optional
+            List of random variables, in string format, to plot.
+        **kwargs : optional
+            Additional key word arguments can be passed to change
+            the numpy.histogram (e.g. density=True, bins=11, ...)
+
+        Returns
+        -------
+        grid_plot : bokeh row
+            A row with the histogram plots.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_disk_example()
+        >>> elm.plot_random_var(["Id"]) # doctest: +ELLIPSIS
+        Row...
+        """
+        label = dict(
+            m="Mass", Id="Diametral moment of inertia", Ip="Polar moment of inertia",
+        )
+        if not all(var in self.is_random for var in var_list):
+            raise ValueError(
+                "Not random variable in var_list. Select variables from {}".format(
+                    self.is_random
+                )
+            )
+
+        return plot_histogram(self.attribute_dict, label, var_list, **kwargs)
+
     @classmethod
     def from_geometry(
         cls, n, material, width, i_d, o_d, tag=None, is_random=None,
@@ -222,7 +273,7 @@ class ST_DiskElement:
         ...                                         o_d=o_d,
         ...                                         is_random=["i_d", "o_d"],
         ...                                         )
-        >>> len(list(elms.__iter__()))
+        >>> len(list(iter(elms)))
         5
         """
         if isinstance(material, ST_Material):
@@ -261,3 +312,27 @@ class ST_DiskElement:
         is_random = ["m", "Id", "Ip"]
 
         return cls(n, m, Id, Ip, tag, is_random=is_random)
+
+
+def st_disk_example():
+    """Return an instance of a simple random disk.
+
+    The purpose is to make available a simple model so that doctest can be
+    written using it.
+
+    Returns
+    -------
+    elm : ross.stochastic.ST_DiskElement
+        An instance of a random disk element object.
+
+    Examples
+    --------
+    >>> import ross.stochastic as srs
+    >>> elm = srs.st_disk_example()
+    >>> len(list(iter(elm)))
+    2
+    """
+    elm = ST_DiskElement(
+        n=1, m=[30, 40], Id=[0.2, 0.3], Ip=[0.5, 0.7], is_random=["m", "Id", "Ip"],
+    )
+    return elm
diff --git a/ross/stochastic/st_materials.py b/ross/stochastic/st_materials.py
index a1d1386f..1e82bdfe 100644
--- a/ross/stochastic/st_materials.py
+++ b/ross/stochastic/st_materials.py
@@ -8,6 +8,7 @@ from collections.abc import Iterable
 import numpy as np
 
 from ross.materials import Material
+from ross.stochastic.st_results_elements import plot_histogram
 
 __all__ = ["ST_Material"]
 
@@ -55,7 +56,7 @@ class ST_Material:
     >>> import ross.stochastic as srs
     >>> E = np.random.uniform(208e9, 211e9, 5)
     >>> st_steel = srs.ST_Material(name="Steel", rho=7810, E=E, G_s=81.2e9)
-    >>> len(list(st_steel.__iter__()))
+    >>> len(list(iter(st_steel)))
     5
     """
 
@@ -89,19 +90,29 @@ class ST_Material:
         if type(Poisson) == list:
             Poisson = np.asarray(Poisson)
 
-        self.name = name
-        self.rho = rho
-        self.E = E
-        self.G_s = G_s
-        self.Poisson = Poisson
-        self.color = color
-
         attribute_dict = dict(
             name=name, rho=rho, E=E, G_s=G_s, Poisson=Poisson, color=color,
         )
         self.is_random = is_random
         self.attribute_dict = attribute_dict
 
+    def __iter__(self):
+        """Return an iterator for the container.
+
+        Returns
+        -------
+        An iterator over random material properties.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> E = np.random.uniform(208e9, 211e9, 5)
+        >>> st_steel = srs.ST_Material(name="Steel", rho=7810, E=E, G_s=81.2e9)
+        >>> len(list(iter(st_steel)))
+        5
+        """
+        return iter(self.random_var(self.is_random, self.attribute_dict))
+
     def __getitem__(self, key):
         """Return the value for a given key from attribute_dict.
 
@@ -143,7 +154,7 @@ class ST_Material:
         key : str
             A class parameter as string.
         value : The corresponding value for the attrbiute_dict's key.
-            ***check the correct type for each key in ST_ShaftElement
+            ***check the correct type for each key in ST_Material
             docstring.
 
         Raises
@@ -168,14 +179,14 @@ class ST_Material:
         """Generate a list of objects as random attributes.
 
         This function creates a list of objects with random values for selected
-        attributes from ShaftElement.
+        attributes from ross.Material.
 
         Parameters
         ----------
         is_random : list
             List of the object attributes to become stochastic.
         *args : dict
-            Dictionary instanciating the ShaftElement class.
+            Dictionary instanciating the ross.Material class.
             The attributes that are supposed to be stochastic should be
             set as lists of random variables.
 
@@ -198,11 +209,46 @@ class ST_Material:
 
         return f_list
 
-    def __iter__(self):
-        """Return an iterator for the container.
+    def plot_random_var(self, var_list=[], **kwargs):
+        """Plot histogram and the PDF.
+
+        This function creates a histogram to display the random variable
+        distribution.
+
+        Parameters
+        ----------
+        var_list : list, optional
+            List of random variables, in string format, to plot.
+        **kwargs : optional
+            Additional key word arguments can be passed to change
+            the numpy.histogram (e.g. density=True, bins=11, ...)
 
         Returns
         -------
-        An iterator over random material properties.
+        grid_plot : bokeh row
+            A row with the histogram plots.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> E = np.random.uniform(208e9, 211e9, 5)
+        >>> st_steel = ST_Material(name="Steel", rho=7810, E=E, G_s=81.2e9)
+        >>> st_steel.plot_random_var(["E"]) # doctest: +ELLIPSIS
+        Row...
         """
-        return iter(self.random_var(self.is_random, self.attribute_dict))
+        label = dict(
+            E="Young's Modulus",
+            G_s="Shear Modulus",
+            Poisson="Poisson coefficient",
+            rho="density",
+        )
+        is_random = self.is_random
+
+        if not all(var in self.is_random for var in var_list):
+            raise ValueError(
+                "Not random variable in var_list. Select variables from {}".format(
+                    is_random
+                )
+            )
+
+        return plot_histogram(self.attribute_dict, label, var_list, **kwargs)
diff --git a/ross/stochastic/st_point_mass.py b/ross/stochastic/st_point_mass.py
index 19fc284a..9d0d6e46 100644
--- a/ross/stochastic/st_point_mass.py
+++ b/ross/stochastic/st_point_mass.py
@@ -1,4 +1,12 @@
+"""Point mass element module for STOCHASTIC ROSS.
+
+This module creates an instance of random point mass for stochastic
+analysis.
+"""
 from ross.point_mass import PointMass
+from ross.stochastic.st_results_elements import plot_histogram
+
+__all__ = ["ST_PointMass", "st_pointmass_example"]
 
 
 class ST_PointMass:
@@ -41,7 +49,7 @@ class ST_PointMass:
     ...                         my=np.random.uniform(2.0, 2.5, 5),
     ...                         is_random=["mx", "my"],
     ...                         )
-    >>> len(list(elms.__iter__()))
+    >>> len(list(iter(elms)))
     5
     """
 
@@ -59,6 +67,13 @@ class ST_PointMass:
         Returns
         -------
         An iterator over random point mass elements.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_pointmass_example()
+        >>> len(list(iter(elm)))
+        2
         """
         return iter(self.random_var(self.is_random, self.attribute_dict))
 
@@ -107,7 +122,7 @@ class ST_PointMass:
         key : str
             A class parameter as string.
         value : The corresponding value for the attrbiute_dict's key.
-            ***check the correct type for each key in ST_ShaftElement
+            ***check the correct type for each key in ST_PointMass
             docstring.
 
         Raises
@@ -136,14 +151,14 @@ class ST_PointMass:
         """Generate a list of objects as random attributes.
 
         This function creates a list of objects with random values for selected
-        attributes from PointMass.
+        attributes from ross.PointMass.
 
         Parameters
         ----------
         is_random : list
             List of the object attributes to become stochastic.
         *args : dict
-            Dictionary instanciating the PointMass class.
+            Dictionary instanciating the ross.PointMass class.
             The attributes that are supposed to be stochastic should be
             set as lists of random variables.
 
@@ -165,3 +180,65 @@ class ST_PointMass:
         f_list = (PointMass(*arg) for arg in new_args)
 
         return f_list
+
+    def plot_random_var(self, var_list=[], **kwargs):
+        """Plot histogram and the PDF.
+
+        This function creates a histogram to display the random variable
+        distribution.
+
+        Parameters
+        ----------
+        var_list : list, optional
+            List of random variables, in string format, to plot.
+        **kwargs : optional
+            Additional key word arguments can be passed to change
+            the numpy.histogram (e.g. density=True, bins=11, ...)
+
+        Returns
+        -------
+        grid_plot : bokeh row
+            A row with the histogram plots.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_pointmass_example()
+        >>> elm.plot_random_var(["mx"]) # doctest: +ELLIPSIS
+        Row...
+        """
+        label = dict(
+            mx="Mass on the X direction", my="Mass on the Y direction", m="Mass",
+        )
+        if not all(var in self.is_random for var in var_list):
+            raise ValueError(
+                "Not random variable in var_list. Select variables from {}".format(
+                    self.is_random
+                )
+            )
+
+        return plot_histogram(self.attribute_dict, label, var_list, **kwargs)
+
+
+def st_pointmass_example():
+    """Return an instance of a simple random point mass.
+
+    The purpose is to make available a simple model so that doctest can be
+    written using it.
+
+    Returns
+    -------
+    elm : ross.stochastic.ST_PointMass
+        An instance of a random point mass element object.
+
+    Examples
+    --------
+    >>> import ross.stochastic as srs
+    >>> elm = srs.st_pointmass_example()
+    >>> len(list(iter(elm)))
+    2
+    """
+    mx = [2.0, 2.5]
+    my = [3.0, 3.5]
+    elm = ST_PointMass(n=1, mx=mx, my=my, is_random=["mx", "my"])
+    return elm
diff --git a/ross/stochastic/st_results_elements.py b/ross/stochastic/st_results_elements.py
new file mode 100644
index 00000000..537849a1
--- /dev/null
+++ b/ross/stochastic/st_results_elements.py
@@ -0,0 +1,89 @@
+"""Plotting module for elements.
+
+This modules provides functions to plot the elements statistic data.
+"""
+import numpy as np
+from bokeh.layouts import gridplot
+from bokeh.plotting import figure
+from scipy.stats import gaussian_kde
+
+
+def plot_histogram(attribute_dict, label={}, var_list=[], **kwargs):
+    """Plot histogram and the PDF.
+
+    This function creates a histogram to display the random variable
+    distribution.
+
+    Parameters
+    ----------
+    attribute_dict : dict
+        Dictionary with element parameters.
+    label : dict
+        Dictionary with labels for each element parameter. Labels are displayed
+        on bokeh figure.
+    var_list : list, optional
+        List of random variables, in string format, to plot.
+    **kwargs : optional
+        Additional key word arguments can be passed to change
+        the numpy.histogram (e.g. density=True, bins=11, ...)
+
+    Returns
+    -------
+    grid_plot : bokeh row
+        A row with the histogram plots.
+    """
+    default_values = dict(density=True, bins=21)
+    for k, v in default_values.items():
+        kwargs.setdefault(k, v)
+
+    figures = []
+
+    for var in var_list:
+        hist, edges = np.histogram(attribute_dict[var], **kwargs)
+        if kwargs["density"] is True:
+            y_label = "PDF"
+        else:
+            y_label = "Frequency"
+
+        fig = figure(
+            width=640,
+            height=480,
+            title="Histogram - {}".format(var),
+            x_axis_label="{}".format(label[var]),
+            y_axis_label="{}".format(y_label),
+        )
+        fig.xaxis.axis_label_text_font_size = "14pt"
+        fig.yaxis.axis_label_text_font_size = "14pt"
+        fig.axis.major_label_text_font_size = "14pt"
+        fig.title.text_font_size = "14pt"
+
+        fig.quad(
+            top=hist,
+            bottom=0,
+            left=edges[:-1],
+            right=edges[1:],
+            fill_color="red",
+            line_color="white",
+            alpha=1.0,
+        )
+        if y_label == "PDF":
+            x = np.linspace(
+                min(attribute_dict[var]),
+                max(attribute_dict[var]),
+                len(attribute_dict[var]),
+            )
+            kernel = gaussian_kde(attribute_dict[var])
+            fig.line(
+                x,
+                kernel(x),
+                line_alpha=1.0,
+                line_width=3.0,
+                line_color="royalblue",
+                legend_label="pdf estimation",
+            )
+
+        figures.append(fig)
+
+    grid_plot = gridplot([figures], toolbar_location="right")
+
+    return grid_plot
diff --git a/ross/stochastic/st_shaft_element.py b/ross/stochastic/st_shaft_element.py
index fe29e317..7806eeb5 100644
--- a/ross/stochastic/st_shaft_element.py
+++ b/ross/stochastic/st_shaft_element.py
@@ -5,8 +5,9 @@ analysis.
 """
 from ross.shaft_element import ShaftElement
 from ross.stochastic.st_materials import ST_Material
+from ross.stochastic.st_results_elements import plot_histogram
 
-__all__ = ["ST_ShaftElement"]
+__all__ = ["ST_ShaftElement", "st_shaft_example"]
 
 
 class ST_ShaftElement:
@@ -71,7 +72,7 @@ class ST_ShaftElement:
     ...                            material=st_steel,
     ...                            is_random=["odl", "material"],
     ...                            )
-    >>> len(list(elms.__iter__()))
+    >>> len(list(iter(elms)))
     5
     """
 
@@ -127,6 +128,13 @@ class ST_ShaftElement:
         Returns
         -------
         An iterator over random shaft elements.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_shaft_example()
+        >>> len(list(iter(elm)))
+        2
         """
         return iter(self.random_var(self.is_random, self.attribute_dict))
 
@@ -212,14 +220,14 @@ class ST_ShaftElement:
         """Generate a list of objects as random attributes.
 
         This function creates a list of objects with random values for selected
-        attributes from ShaftElement.
+        attributes from ross.ShaftElement.
 
         Parameters
         ----------
         is_random : list
             List of the object attributes to become stochastic.
         *args : dict
-            Dictionary instanciating the ShaftElement class.
+            Dictionary instanciating the ross.ShaftElement class.
             The attributes that are supposed to be stochastic should be
             set as lists of random variables.
 
@@ -241,3 +249,75 @@ class ST_ShaftElement:
         f_list = (ShaftElement(*arg) for arg in new_args)
 
         return f_list
+
+    def plot_random_var(self, var_list=[], **kwargs):
+        """Plot histogram and the PDF.
+
+        This function creates a histogram to display the random variable
+        distribution.
+
+        Parameters
+        ----------
+        var_list : list, optional
+            List of random variables, in string format, to plot.
+        **kwargs : optional
+            Additional key word arguments can be passed to change
+            the numpy.histogram (e.g. density=True, bins=11, ...)
+
+        Returns
+        -------
+        grid_plot : bokeh row
+            A row with the histogram plots.
+
+        Examples
+        --------
+        >>> import ross.stochastic as srs
+        >>> elm = srs.st_shaft_example()
+        >>> elm.plot_random_var(["odl"]) # doctest: +ELLIPSIS
+        Row...
+        """
+        label = dict(
+            L="Length",
+            idl="Left inner diameter",
+            odl="Left outer diameter",
+            idr="Right inner diameter",
+            odr="Right outer diameter",
+        )
+        is_random = self.is_random
+        if "material" in is_random:
+            is_random.remove("material")
+
+        if not all(var in self.is_random for var in var_list):
+            raise ValueError(
+                "Not random variable in var_list. Select variables from {}".format(
+                    is_random
+                )
+            )
+
+        return plot_histogram(self.attribute_dict, label, var_list, **kwargs)
+
+
+def st_shaft_example():
+    """Return an instance of a simple random shaft element.
+
+    The purpose is to make available a simple model so that doctest can be
+    written using it.
+
+    Returns
+    -------
+    elm : ross.stochastic.ST_ShaftElement
+        An instance of a random shaft element object.
+
+    Examples
+    --------
+    >>> import ross.stochastic as srs
+    >>> elm = srs.st_shaft_example()
+    >>> len(list(iter(elm)))
+    2
+    """
+    from ross.materials import steel
+
+    elm = ST_ShaftElement(
+        L=[1.0, 1.1], idl=0.0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"],
+    )
+    return elm


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index f3e684fe..ccc6e6e4 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -6,7 +6,6 @@ import pytest
 from numpy.testing import assert_allclose
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
-                                       BearingElement6DoF,
                                        RollerBearingElement)
 
 # fmt: on
@@ -288,42 +287,3 @@ 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 c46a57aa..c1dcc4a1 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, DiskElement6DoF
+from ross.disk_element import DiskElement
 from ross.materials import steel
 
 
@@ -75,85 +75,3 @@ 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 22392e23..75c611a0 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, ShaftElement6DoF
+from ross.shaft_element import ShaftElement
 
 test_dir = os.path.dirname(__file__)
 
@@ -306,121 +306,3 @@ 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)
diff --git a/ross/tests/test_stochastic_elements.py b/ross/tests/test_stochastic_elements.py
index 92be7ec7..d12719f5 100644
--- a/ross/tests/test_stochastic_elements.py
+++ b/ross/tests/test_stochastic_elements.py
@@ -7,6 +7,7 @@ Tests for:
 """
 import numpy as np
 import pytest
+from bokeh.layouts import row
 
 from ross.stochastic.st_bearing_seal_element import ST_BearingElement
 from ross.stochastic.st_disk_element import ST_DiskElement
@@ -97,8 +98,11 @@ def rand_point_mass():
     return elm
 
 
+###############################################################################
+# testing attributes and parameters
+###############################################################################
 def test_st_shaft_element(rand_shaft):
-    elm = list(rand_shaft.__iter__())
+    elm = list(iter(rand_shaft))
     assert [sh.L for sh in elm] == [1.0, 1.1]
     assert [sh.idl for sh in elm] == [0.01, 0.02]
     assert [sh.odl for sh in elm] == [0.1, 0.2]
@@ -108,7 +112,7 @@ def test_st_shaft_element(rand_shaft):
 
 
 def test_st_disk_element_from_inertia(rand_disk_from_inertia):
-    elm = list(rand_disk_from_inertia.__iter__())
+    elm = list(iter(rand_disk_from_inertia))
     assert [dk.n for dk in elm] == [1, 1]
     assert [dk.m for dk in elm] == [30, 40]
     assert [dk.Id for dk in elm] == [0.2, 0.3]
@@ -116,7 +120,7 @@ def test_st_disk_element_from_inertia(rand_disk_from_inertia):
 
 
 def test_st_disk_element_from_geometry(rand_disk_from_geometry):
-    elm = list(rand_disk_from_geometry.__iter__())
+    elm = list(iter(rand_disk_from_geometry))
     assert [dk.n for dk in elm] == [1, 1]
     assert [dk.m for dk in elm] == [37.570502893821185, 76.74810321754951]
     assert [dk.Id for dk in elm] == [0.2325457585365474, 0.8256816771154702]
@@ -124,7 +128,7 @@ def test_st_disk_element_from_geometry(rand_disk_from_geometry):
 
 
 def test_st_bearing_element_constant_coef(rand_bearing_constant_coefs):
-    elm = list(rand_bearing_constant_coefs.__iter__())
+    elm = list(iter(rand_bearing_constant_coefs))
     assert [brg.n for brg in elm] == [1, 1]
     assert [brg.kxx.coefficient for brg in elm] == [[1000000.0], [2000000.0]]
     assert [brg.kyy.coefficient for brg in elm] == [[1000000.0], [2000000.0]]
@@ -137,7 +141,7 @@ def test_st_bearing_element_constant_coef(rand_bearing_constant_coefs):
 
 
 def test_st_bearing_element_varying_coef(rand_bearing_varying_coefs):
-    elm = list(rand_bearing_varying_coefs.__iter__())
+    elm = list(iter(rand_bearing_varying_coefs))
     assert [brg.n for brg in elm] == [1, 1]
     assert [brg.kxx.coefficient for brg in elm] == [
         [1000000.0, 2000000.0],
@@ -162,7 +166,7 @@ def test_st_bearing_element_varying_coef(rand_bearing_varying_coefs):
 
 
 def test_st_point_mass(rand_point_mass):
-    elm = list(rand_point_mass.__iter__())
+    elm = list(iter(rand_point_mass))
     assert [pm.n for pm in elm] == [1, 1]
     assert [pm.mx for pm in elm] == [2.0, 2.5]
     assert [pm.my for pm in elm] == [3.0, 3.5]
@@ -171,9 +175,7 @@ def test_st_point_mass(rand_point_mass):
 ###############################################################################
 # testing error messages
 ###############################################################################
-
-
-def test_st_bearing_error_messages():
+def test_st_bearing_error_messages(rand_bearing_constant_coefs):
     kxx = [1e6, 2e6]
     cxx = [1e3, 2e3]
     freq = [500, 1000]
@@ -182,3 +184,88 @@ def test_st_bearing_error_messages():
             n=1, kxx=kxx, cxx=cxx, frequency=freq, is_random=["kxx", "cxx", "frequency"]
         )
     assert "frequency can not be a random variable" in str(ex.value)
+
+    with pytest.raises(ValueError) as ex:
+        rand_bearing_constant_coefs.plot_random_var(["kxy"])
+    assert (
+        "Not random variable in var_list. Select variables from ['kxx', 'cxx', 'kyy', 'cyy']"
+        in str(ex.value)
+    )
+
+    with pytest.raises(KeyError) as ex:
+        rand_bearing_constant_coefs["odd"] = [1, 2]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+    with pytest.raises(KeyError) as ex:
+        rand_bearing_constant_coefs["odd"]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+
+def test_st_disk_error_messages(rand_disk_from_inertia):
+    with pytest.raises(ValueError) as ex:
+        rand_disk_from_inertia.plot_random_var(["n"])
+    assert (
+        "Not random variable in var_list. Select variables from ['m', 'Id', 'Ip']"
+        in str(ex.value)
+    )
+
+    with pytest.raises(KeyError) as ex:
+        rand_disk_from_inertia["odd"] = [1, 2]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+    with pytest.raises(KeyError) as ex:
+        rand_disk_from_inertia["odd"]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+
+def test_st_shaft_error_messages(rand_shaft):
+    with pytest.raises(ValueError) as ex:
+        rand_shaft.plot_random_var(["n"])
+    assert (
+        "Not random variable in var_list. Select variables from ['L', 'idl', 'odl', 'idr', 'odr']"
+        in str(ex.value)
+    )
+
+    with pytest.raises(KeyError) as ex:
+        rand_shaft["odd"] = [1, 2]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+    with pytest.raises(KeyError) as ex:
+        rand_shaft["odd"]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+
+def test_st_point_mass_error_messages(rand_point_mass):
+    with pytest.raises(ValueError) as ex:
+        rand_point_mass.plot_random_var(["n"])
+    assert "Not random variable in var_list. Select variables from ['mx', 'my']" in str(
+        ex.value
+    )
+
+    with pytest.raises(KeyError) as ex:
+        rand_point_mass["odd"] = [1, 2]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+    with pytest.raises(KeyError) as ex:
+        rand_point_mass["odd"]
+    assert "Object does not have parameter: odd." in str(ex.value)
+
+
+###############################################################################
+# testing elements plot type
+###############################################################################
+def test_element_stats_plot(
+    rand_shaft, rand_disk_from_inertia, rand_bearing_constant_coefs, rand_point_mass
+):
+    figure_type = type(row())
+    fig = rand_shaft.plot_random_var(["L"])
+    assert type(fig) == figure_type
+
+    fig = rand_disk_from_inertia.plot_random_var(["Id"])
+    assert type(fig) == figure_type
+
+    fig = rand_bearing_constant_coefs.plot_random_var(["kxx"])
+    assert type(fig) == figure_type
+
+    fig = rand_point_mass.plot_random_var(["mx"])
+    assert type(fig) == figure_type