Repository: ross-rotordynamics/ross
Pull request #109 (2019-05-07 18:56:25)

Add tests for static analysis and convergence analysis
This Pull Request is to replace the last PR related to this issue. Since the last one was getting too difficult to solve the conflicts. 
EDIT 1: I'm also adding the Bokeh plots in this PR instead of opening a new branch / PR to avoid future conflicts.
EDIT 2: The function interp1d from module scipy.interpolate was not
being imported. This PR also corrects this problem.


Patch:
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index f9186997..f1dd8695 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -86,6 +86,8 @@ class BearingElement(Element):
     Values for each parameter will be interpolated for the speed.
     Parameters
     ----------
+    n: int
+        Node which the bearing will be located in
     kxx: float, array
         Direct stiffness in the x direction.
     cxx: float, array
@@ -112,6 +114,19 @@ class BearingElement(Element):
         Array with the speeds (rad/s).
     Examples
     --------
+    >>> # A bearing element located in the first rotor node, with these
+    >>> # following stiffness and damping coefficients and speed range from
+    >>> # 0 to 200 rad/s
+    >>> Kxx = 1e6
+    >>> Kyy = 0.8e6
+    >>> Cxx = 2e2
+    >>> Cyy = 1.5e2
+    >>> W = np.linspace(0,200,11)
+    >>> bearing0 = rs.BearingElement(n=0, kxx=Kxx, kyy=Kyy, cxx=Cxx, cyy=Cyy, w=W)
+    >>> bearing0.K(W)
+    array([[[1000000., 1000000., ... , 800000.,  800000.]]])
+    >>> bearing0.C(W)
+    array([[[200., 200., ... , 150., 150.]]])
     """
 
     def __init__(
@@ -468,4 +483,4 @@ class SealElement(BearingElement):
         bk_ax.patch(
             bk_seal_points_l[0], bk_seal_points_l[1],
             alpha=0.5, line_width=2, color=bokeh_colors[6]
-        )
\ No newline at end of file
+        )
diff --git a/ross/results.py b/ross/results.py
index 917c732c..d46f3171 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -3,6 +3,14 @@ import numpy as np
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
+import bokeh.palettes as bp
+from bokeh.transform import linear_cmap
+from bokeh.layouts import gridplot
+from bokeh.plotting import figure, output_file, show
+from bokeh.models import ColumnDataSource, ColorBar
+
+# set bokeh palette of colors
+bokeh_colors = bp.RdGy[11]
 
 
 class Results(np.ndarray):
@@ -90,7 +98,26 @@ class CampbellResults(Results):
         for k, v in default_values.items():
             kwargs.setdefault(k, v)
 
-        for mark, whirl_dir in zip(["^", "o", "v"], [0.0, 0.5, 1.0]):
+        # bokeh plot - output to static HTML file
+        output_file("Campbell_diagram.html")
+
+        log_dec_map = log_dec.flatten()
+
+        # bokeh plot - create a new plot
+        camp = figure(
+           tools="pan, box_zoom, wheel_zoom, reset, save",
+           width=900,
+           height=500,
+           title="Campbell Diagram",
+           x_axis_label="Rotor speed (rad/s)",
+           y_axis_label="Damped natural frequencies (rad/s)"
+        )
+
+        for mark, whirl_dir, legend in zip(
+                ["^", "o", "v"],
+                [0.0, 0.5, 1.0],
+                ["Foward", "Mixed", "Backward"]
+            ):
             num_frequencies = wd.shape[1]
             for i in range(num_frequencies):
                 if wn is True:
@@ -113,6 +140,61 @@ class CampbellResults(Results):
                         **kwargs,
                     )
 
+                    # Bokeh plot
+                    source = ColumnDataSource(
+                            dict(x=speed_range_i[whirl_mask],
+                                 y=w_i[whirl_mask],
+                                 color=log_dec_i[whirl_mask]
+                                 )
+                    )
+                    color_mapper = linear_cmap(
+                            field_name='color',
+                            palette=bp.magma(256),
+                            low=min(log_dec_map),
+                            high=max(log_dec_map)
+                    )
+                    camp.scatter(
+                        x='x',
+                        y='y',
+                        color=color_mapper,
+                        marker=mark,
+                        fill_alpha=1.0,
+                        size=5,
+                        muted_color=color_mapper,
+                        muted_alpha=0.2, 
+                        source=source,
+                        legend=legend
+                    )
+
+        ax.plot(speed_range[:, 0],
+                speed_range[:, 0],
+                color='k',
+                linewidth=1.5,
+                linestyle="-.",
+                alpha=0.75,
+                label="Rotor speed"
+        )
+
+        camp.line(x=speed_range[:, 0],
+                  y=speed_range[:, 0],
+                  line_width=3,
+                  color=bokeh_colors[0],
+                  line_dash="dotdash",
+                  line_alpha=0.75,
+                  legend="Rotor speed",
+                  muted_color=bokeh_colors[0],
+                  muted_alpha=0.2,
+        )
+
+        color_bar = ColorBar(
+                color_mapper=color_mapper['transform'], width=8, location=(0, 0)
+        )
+
+        camp.legend.click_policy = "mute"
+        camp.legend.location = "top_left"
+        camp.add_layout(color_bar, 'right')
+        show(camp)
+
         if len(fig.axes) == 1:
             cbar = fig.colorbar(im)
             cbar.ax.set_ylabel("log dec")
@@ -160,7 +242,9 @@ class FrequencyResponseResults(Results):
         Returns
         -------
         ax : matplotlib.axes
-            Matplotlib axes with phase plot.
+            Matplotlib axes with amplitude plot.
+        mag_plot : bokeh plot axes
+            Bokeh plot axes with amplitude plot.
         Examples
         --------
         """
@@ -178,14 +262,35 @@ class FrequencyResponseResults(Results):
 
         if units == "m":
             ax.set_ylabel("Amplitude $(m)$")
+            y_axis_label = "Amplitude (m)"
         elif units == "mic-pk-pk":
             ax.set_ylabel("Amplitude $(\mu pk-pk)$")
+            y_axis_label = ("Amplitude (\mu pk-pk)")
         else:
             ax.set_ylabel("Amplitude $(dB)$")
+            y_axis_label = ("Amplitude (dB")
 
         ax.set_xlabel("Frequency (rad/s)")
 
-        return ax
+        # bokeh plot - create a new plot
+        mag_plot = figure(
+                tools="pan, box_zoom, wheel_zoom, reset, save",
+                width=900,
+                height=400,
+                title="Frequency Response - Magnitude",
+                x_axis_label="Frequency",
+                y_axis_label=y_axis_label
+        )
+        source = ColumnDataSource(dict(x=frequency_range, y=mag[inp, out, :]))
+        mag_plot.line(x='x',
+                      y='y',
+                      source=source,
+                      line_color=bokeh_colors[0],
+                      line_alpha=1.0,
+                      line_width=3,
+                      )
+
+        return ax, mag_plot
 
     def plot_phase(self, inp, out, ax=None, **kwargs):
         """Plot frequency response.
@@ -205,10 +310,10 @@ class FrequencyResponseResults(Results):
             the plot (e.g. linestyle='--')
         Returns
         -------
-        ax0 : matplotlib.axes
-            Matplotlib axes with amplitude plot.
-        ax1 : matplotlib.axes
+        ax : matplotlib.axes
             Matplotlib axes with phase plot.
+        phase_plot : bokeh plot axes
+            Bokeh plot axes with phase plot.
         Examples
         --------
         """
@@ -227,7 +332,25 @@ class FrequencyResponseResults(Results):
         ax.set_ylabel("Phase")
         ax.set_xlabel("Frequency (rad/s)")
 
-        return ax
+        # bokeh plot - create a new plot
+        phase_plot = figure(
+                tools="pan, box_zoom, wheel_zoom, reset, save",
+                width=900,
+                height=400,
+                title="Frequency Response - Phase",
+                x_axis_label="Frequency",
+                y_axis_label="Phase"
+        )
+        source = ColumnDataSource(dict(x=frequency_range, y=phase[inp, out, :]))
+        phase_plot.line(x='x',
+                        y='y',
+                        source=source,
+                        line_color=bokeh_colors[0],
+                        line_alpha=1.0,
+                        line_width=3,
+                        )
+
+        return ax, phase_plot
 
     def plot(self, inp, out, ax0=None, ax1=None, **kwargs):
         """Plot frequency response.
@@ -239,11 +362,11 @@ class FrequencyResponseResults(Results):
             Input.
         out : int
             Output.
-        ax0 : matplotlib.axes, optional
-            Matplotlib axes where the amplitude will be plotted.
+        ax0 : matplotlib.axes, bokeh plot axes optional
+            Matplotlib and bokeh plot axes where the amplitude will be plotted.
             If None creates a new.
-        ax1 : matplotlib.axes, optional
-            Matplotlib axes where the phase will be plotted.
+        ax1 : matplotlib.axes, bokeh plot axes optional
+            Matplotlib and bokeh plot axes where the phase will be plotted.
             If None creates a new.
         kwargs : optional
             Additional key word arguments can be passed to change
@@ -254,18 +377,34 @@ class FrequencyResponseResults(Results):
             Matplotlib axes with amplitude plot.
         ax1 : matplotlib.axes
             Matplotlib axes with phase plot.
+        bk_ax0 : bokeh plot axes
+            Bokeh plot axes with amplitude plot
+        bk_ax1 : bokeh plot axes
+            Bokeh plot axes with phase plot
         Examples
         --------
         """
         if ax0 is None and ax1 is None:
             fig, (ax0, ax1) = plt.subplots(2)
 
-        ax0 = self.plot_magnitude(inp, out, ax=ax0)
-        ax1 = self.plot_phase(inp, out, ax=ax1)
+        # bokeh plot - output to static HTML file
+        output_file("plot_freq_response.html")
+
+        # matplotlib axes
+        ax0 = self.plot_magnitude(inp, out, ax=ax0)[0]
+        ax1 = self.plot_phase(inp, out, ax=ax1)[0]
+
+        # bokeh plot axes
+        bk_ax0 = self.plot_magnitude(inp, out, ax=ax0)[1]
+        bk_ax1 = self.plot_phase(inp, out, ax=ax1)[1]
 
         ax0.set_xlabel("")
 
-        return ax0, ax1
+        # show the bokeh plot results
+        grid_plots = gridplot([[bk_ax0], [bk_ax1]])
+        show(grid_plots)
+
+        return ax0, ax1, bk_ax0, bk_ax1
 
     def plot_freq_response_grid(self, outs, inps, ax=None, **kwargs):
 
@@ -559,7 +698,7 @@ class ModeShapeResults(Results):
 
         ax.set_title(
             f"$speed$ = {self.w:.1f} rad/s\n$"
-            f"\omega_d$ = {self.wd[mode]:.1f} rad/s\n"
+            f"\frequency_range_d$ = {self.wd[mode]:.1f} rad/s\n"
             f"$log dec$ = {self.log_dec[mode]:.1f}"
         )
 
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 73d38a74..69a78e26 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -6,6 +6,7 @@ import scipy.linalg as la
 import scipy.sparse.linalg as las
 import scipy.signal as signal
 import scipy.io as sio
+import scipy.interpolate as interpolate
 from copy import copy, deepcopy
 from collections import Iterable
 import shutil
@@ -279,6 +280,11 @@ class Rotor(object):
         # number of dofs
         self.ndof = 4 * max([el.n for el in shaft_elements]) + 8
 
+        #  values for static analysis will be calculated by def static
+        self.Vx = None
+        self.Bm = None
+        self.disp_y = None
+        
         #  diameter at node position
 
         print(
@@ -307,6 +313,146 @@ class Rotor(object):
             )
         self.lti = self._lti()
 
+    def convergence(self, n_eigval=0, err_max=1e-02):
+        """
+        Function to analyze the eigenvalues convergence through the number of
+        shaft elements. Every new run doubles the number os shaft elements.
+
+        -----------
+        Parameters:
+        -----------
+        n_eigval : int
+            The nth eigenvalue which the convergence analysis will run.
+            Default is 0 (the first eigenvalue).
+        err_max : float
+            Maximum allowable convergence error.
+            Default is 1e-02
+
+        Example:
+        --------
+        >>> i_d = 0
+        >>> o_d = 0.05
+        >>> n = 6
+        >>> L = [0.25 for _ in range(n)]
+        ...
+        >>> shaft_elem = [rs.ShaftElement(l, i_d, o_d, steel,
+        shear_effects=True, rotary_inertia=True, gyroscopic=True) for l in L]
+        ...
+        >>> disk0 = DiskElement.from_geometry(2, steel, 0.07, 0.05, 0.28)
+        >>> disk1 = DiskElement.from_geometry(4, steel, 0.07, 0.05, 0.35)
+        >>> bearing0 = BearingElement(0, kxx=1e6, kyy=8e5, cxx=2e3)
+        >>> bearing1 = BearingElement(6, kxx=1e6, kyy=8e5, cxx=2e3)
+        >>> rotor0 = Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+        >>> len(rotor0.shaft_elements)
+        6
+        >>> rotor0.convergence(n_eigval=0, err_max=1e-08)
+        >>> len(rotor0.shaft_elements)
+        96
+        """
+        el_num = np.array([len(self.shaft_elements)])
+        eigv_arr = np.array([])
+        error_arr = np.array([0])
+
+        self.run()
+        eigv_arr = np.append(eigv_arr, self.wn[n_eigval])
+
+        # this value is up to start the loop while
+        error = 1
+        nel_r = 2
+
+        while error > err_max:
+            shaft_elem = []
+            disk_elem = []
+            brgs_elem = []
+
+            for i, leng in enumerate(self.shaft_elements):
+                le = self.shaft_elements[i].L / nel_r
+                o_ds = self.shaft_elements[i].o_d
+                i_ds = self.shaft_elements[i].i_d
+
+                # loop to double the number of element
+                for j in range(nel_r):
+                    shaft_elem.append(
+                        ShaftElement(
+                            le,
+                            i_ds,
+                            o_ds,
+                            material=self.shaft_elements[i].material,
+                            shear_effects=True,
+                            rotary_inertia=True,
+                            gyroscopic=True,
+                        )
+                    )
+
+            for DiskEl in self.disk_elements:
+                aux_DiskEl = deepcopy(DiskEl)
+                aux_DiskEl.n = nel_r * DiskEl.n
+                disk_elem.append(aux_DiskEl)
+
+            for Brg_SealEl in self.bearing_seal_elements:
+                aux_Brg_SealEl = deepcopy(Brg_SealEl)
+                aux_Brg_SealEl.n = nel_r * Brg_SealEl.n
+                brgs_elem.append(aux_Brg_SealEl)
+
+            rotor = Rotor(
+                shaft_elem, disk_elem, brgs_elem, w=self.w, n_eigen=self.n_eigen
+            )
+            rotor.run()
+
+            eigv_arr = np.append(eigv_arr, rotor.wn[n_eigval])
+            el_num = np.append(el_num, len(shaft_elem))
+
+            error = min(eigv_arr[-1], eigv_arr[-2]) / max(eigv_arr[-1], eigv_arr[-2])
+            error = 1 - error
+
+            error_arr = np.append(error_arr, 100 * error)
+            nel_r *= 2
+
+        self.__dict__ = rotor.__dict__
+        self.error_arr = error_arr
+
+        output_file("convergence.html")
+        source = ColumnDataSource(
+            data=dict(x0=el_num[1:], y0=eigv_arr[1:], x1=el_num[1:], y1=error_arr[1:])
+        )
+
+        TOOLS = "pan,wheel_zoom,box_zoom,hover,reset,save,"
+        TOOLTIPS1 = [("Frquency:", "@y0 Hz"), ("Number of Elements", "@x0")]
+        TOOLTIPS2 = [("Relative Error:", "@y1"), ("Number of Elements", "@x1")]
+        # create a new plot and add a renderer
+        freq_arr = figure(
+            tools=TOOLS,
+            tooltips=TOOLTIPS1,
+            width=500,
+            height=500,
+            title="Frequency Evaluation",
+            x_axis_label="Numer of Elements",
+            y_axis_label="Frequency (Hz)",
+        )
+        freq_arr.line("x0", "y0", source=source, line_width=3, line_color="crimson")
+        freq_arr.circle("x0", "y0", source=source, fill_color="crimson", size=8)
+
+        # create another new plot and add a renderer
+        rel_error = figure(
+            tools=TOOLS,
+            tooltips=TOOLTIPS2,
+            width=500,
+            height=500,
+            title="Relative Error Evaluation",
+            x_axis_label="Number of Rlements",
+            y_axis_label="Relative Rrror",
+        )
+        rel_error.line(
+            "x1", "y1", source=source, line_width=3, line_color="darkslategray"
+        )
+        rel_error.circle("x1", "y1", source=source, fill_color="darkslategray", size=8)
+
+        # put the subplots in a gridplot
+        p = gridplot([[freq_arr, rel_error]])
+
+        # show the results
+        show(p)
+
     @property
     def w(self):
         return self._w
@@ -1564,6 +1710,7 @@ class Rotor(object):
         disp_y = np.array([])
         for node_dof in range(int(len(disp) / 4)):
             disp_y = np.append(disp_y, disp[4 * node_dof + 1])
+        self.disp_y = disp_y
 
         # Shearing Force
         BrgForce = [0] * len(self.nodes_pos)
@@ -1602,7 +1749,8 @@ class Rotor(object):
                 DskForce.insert(i + 1, 0)
                 SchForce.insert(i + 1, 0)
                 Vx_axis.insert(i, Vx_axis[i])
-        Vx = [x * -1 for x in Vx]
+        self.Vx = [x * -1 for x in Vx]
+        Vx = self.Vx
 
         # Bending Moment vector
         Mx = []
@@ -1622,6 +1770,7 @@ class Rotor(object):
         Bm = [0]
         for i in range(len(Mx)):
             Bm.append(Bm[i] + Mx[i])
+        self.Bm = Bm
 
         output_file("static.html")
         source = ColumnDataSource(
@@ -1630,9 +1779,9 @@ class Rotor(object):
 
         TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select,hover"
         TOOLTIPS = [
-            ("Shaft lenght:", "@x0 m"),
-            ("Underformed:", "@y1 mm"),
-            ("Displacement:", "@y0 mm"),
+            ("Shaft lenght:", "@x0"),
+            ("Underformed:", "@y1"),
+            ("Displacement:", "@y0"),
         ]
 
         # create displacement plot
@@ -1642,8 +1791,8 @@ class Rotor(object):
             width=800,
             height=400,
             title="Static Analysis",
-            x_axis_label="shaft lenght(m)",
-            y_axis_label="lateral displacement (mm)",
+            x_axis_label="shaft lenght",
+            y_axis_label="lateral displacement",
         )
 
         interpolated = interpolate.interp1d(source.data['x0'], source.data['y0'], kind='cubic')
@@ -1703,8 +1852,8 @@ class Rotor(object):
             width=800,
             height=400,
             title="Free-Body Diagram",
-            x_axis_label="shaft lenght (m)",
-            y_axis_label="Force (N)",
+            x_axis_label="shaft lenght",
+            y_axis_label="Force",
             x_range=[-0.1 * shaft_end, 1.1 * shaft_end],
             y_range=[-max(y_range) * 1.4, max(y_range) * 1.4],
         )
@@ -1814,15 +1963,15 @@ class Rotor(object):
 
         # Shearing Force Diagram plot (SF)
         source_SF = ColumnDataSource(data=dict(x=Vx_axis, y=Vx))
-        TOOLTIPS_SF = [("Shearing Force:", "@y N")]
+        TOOLTIPS_SF = [("Shearing Force:", "@y")]
         SF = figure(
             tools=TOOLS,
             tooltips=TOOLTIPS_SF,
             width=800,
             height=400,
             title="Shearing Force Diagram",
-            x_axis_label="Shaft lenght (m)",
-            y_axis_label="Force (N)",
+            x_axis_label="Shaft lenght",
+            y_axis_label="Force",
             x_range=[-0.1 * shaft_end, 1.1 * shaft_end],
         )
         SF.line("x", "y", source=source_SF, line_width=4, line_color=bokeh_colors[0])
@@ -1839,15 +1988,15 @@ class Rotor(object):
 
         # Bending Moment Diagram plot (BM)
         source_BM = ColumnDataSource(data=dict(x=self.nodes_pos, y=Bm))
-        TOOLTIPS_BM = [("Bending Moment:", "@y N.m")]
+        TOOLTIPS_BM = [("Bending Moment:", "@y")]
         BM = figure(
             tools=TOOLS,
             tooltips=TOOLTIPS_BM,
             width=800,
             height=400,
             title="Bending Moment Diagram",
-            x_axis_label="Shaft lenght (m)",
-            y_axis_label="Bending Moment (N.m)",
+            x_axis_label="Shaft lenght",
+            y_axis_label="Bending Moment",
             x_range=[-0.1 * shaft_end, 1.1 * shaft_end],
         )
         i=0
@@ -1882,6 +2031,7 @@ class Rotor(object):
         grid_plots = gridplot([[FBD, SF], [disp_graph, BM]])
         show(grid_plots)
 
+
     @classmethod
     def from_section(
         cls,
@@ -1893,18 +2043,14 @@ class Rotor(object):
         sparse=True,
         min_w=None,
         max_w=None,
+        n_eigen=12,
         w=0,
         nel_r=1,
-        n_eigval=1,
-        err_max=1e-02,
     ):
 
         """This class is an alternative to build rotors from separated
         sections. Each section has the same number (n) of shaft elements.
 
-        This class will verify the eigenvalues calculation
-        and check its convergence to minimize the numerical errors.
-
         Parameters
         ----------
         leng_data : list
@@ -1913,27 +2059,29 @@ class Rotor(object):
             List with the outer diameters of rotor regions.
         i_d_data : list
             List with the inner diameters of rotor regions.
-        disk_data : list, optional
-            List holding lists of disks datas.
-            Example : disk_data = [[n, material, width, i_d, o_d], [n, ...]]
+        disk_data : dict, optional
+            Dict holding disks datas.
+            Example : disk_data=DiskElement.from_geometry(n=2,
+                                                          material=steel,
+                                                          width=0.07,
+                                                          i_d=0,
+                                                          o_d=0.28
+                                                          )
             ***See 'disk_element.py' docstring for more information***
-        brg_seal_data : list, optional
-            list holding lists of bearings and seals datas.
-            Example : brg_seal_data=[[n, kxx, cxx, kyy=None, kxy=0, kyx=0,
-                                      cyy=None, cxy=0, cyx=0, w=None],
-                                     [n, ...]]
+        brg_seal_data : dict, optional
+            Dict holding lists of bearings and seals datas.
+            Example : brg_seal_data=BearingElement(n=1, kxx=1e6, cxx=0,
+                                                   kyy=1e6, cyy=0, kxy=0,
+                                                   cxy=0, kyx=0, cyx=0)
             ***See 'bearing_seal_element.py' docstring for more information***
         w : float, optional
             Rotor speed.
-        nel_r : int
-            Initial number or elements per shaft region.
+        nel_r : int, optional
+            Number or elements per shaft region.
             Default is 1
-        eigval : int
-            Indicates which eingenvalue convergence to check.
-            default is 1 (1st eigenvalue).
-        err_max : float, optional
-            maximum allowed for eigenvalues calculation.
-            default is 0.01 (or 1%).
+        n_eigen : int, optional
+            Number of eigenvalues calculated by arpack.
+            Default is 12.
 
         Example
         -------
@@ -1945,11 +2093,11 @@ class Rotor(object):
         ...                        DiskElement.from_geometry(n=2, material=steel, width=0.07, i_d=0, o_d=0.35)],
         ...             brg_seal_data=[BearingElement(n=0, kxx=1e6, cxx=0, kyy=1e6, cyy=0, kxy=0, cxy=0, kyx=0, cyx=0),
         ...                            BearingElement(n=3, kxx=1e6, cxx=0, kyy=1e6, cyy=0, kxy=0, cxy=0, kyx=0, cyx=0)],
-        ...             w=0, nel_r=1, n_eigval=1, err_max=1e-07)
+        ...             w=0, nel_r=1)
+        >>> rotor.run()
         >>> rotor.wn[:]
         """
 
-        n_eigen = 2 * n_eigval + 2
         if len(leng_data) != len(o_ds_data) or len(leng_data) != len(i_ds_data):
             raise ValueError("The matrices lenght do not match.")
 
@@ -1999,79 +2147,11 @@ class Rotor(object):
 
             return regions
 
-        el_num = np.array([nel_r * len(leng_data)])
-        eigv_arr = np.array([])
-        error_arr = np.array([0])
-
-        regions0 = rotor_regions(nel_r)
-        rotor0 = Rotor(regions0[0], regions0[1], regions0[2], w=w, n_eigen=n_eigen)
-        rotor0.run()
-
-        eigv_arr = np.append(eigv_arr, rotor0.wn[n_eigval])
-        # this value is up to start the loop while
-        error = 1
-        nel_r *= 2
-
-        while error > err_max:
-
-            regions = rotor_regions(nel_r)
-            rotor = Rotor(regions[0], regions[1], regions[2], w=w, n_eigen=n_eigen)
-            rotor.run()
-
-            eigv_arr = np.append(eigv_arr, rotor.wn[n_eigval])
-            el_num = np.append(el_num, nel_r * len(leng_data))
-
-            error = min(eigv_arr[-1], eigv_arr[-2]) / max(eigv_arr[-1], eigv_arr[-2])
-            error = 1 - error
-            error_arr = np.append(error_arr, 100 * error)
-            nel_r *= 2
-
+        regions = rotor_regions(nel_r)
         shaft_elements = regions[0]
         disk_elements = regions[1]
         bearing_seal_elements = regions[2]
 
-        output_file("convergence.html")
-        source = ColumnDataSource(
-            data=dict(x0=el_num[1:], y0=eigv_arr[1:], x1=el_num[1:], y1=error_arr[1:])
-        )
-
-        TOOLS = "pan,wheel_zoom,box_zoom,hover,reset,save,"
-        TOOLTIPS1 = [("Frquency:", "@y0 Hz"), ("Number of Elements", "@x0")]
-        TOOLTIPS2 = [("Relative Error:", "@y1"), ("Number of Elements", "@x1")]
-        # create a new plot and add a renderer
-        freq_arr = figure(
-            tools=TOOLS,
-            tooltips=TOOLTIPS1,
-            width=500,
-            height=500,
-            title="Frequency Evaluation",
-            x_axis_label="Numer of Elements",
-            y_axis_label="Frequency (Hz)",
-        )
-        freq_arr.line("x0", "y0", source=source, line_width=3, line_color="crimson")
-        freq_arr.circle("x0", "y0", source=source, fill_color="crimson", size=8)
-
-        # create another new plot and add a renderer
-        rel_error = figure(
-            tools=TOOLS,
-            tooltips=TOOLTIPS2,
-            width=500,
-            height=500,
-            title="Relative Error Evaluation",
-            x_axis_label="Number of Rlements",
-            y_axis_label="Relative Rrror",
-        )
-        rel_error.line(
-            "x1", "y1", source=source, line_width=3, line_color="darkslategray"
-        )
-        rel_error.circle("x1", "y1", source=source, fill_color="darkslategray", size=8)
-
-        # put the subplots in a gridplot
-        p = gridplot([[freq_arr, rel_error]])
-
-        # show the results
-        show(p)
-
         return cls(
             shaft_elements,
             disk_elements,
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 7efd0820..6f6f23a0 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -439,7 +439,7 @@ class ShaftElement(Element):
         )
         # bokeh plot - plot the shaft
         bk_ax.quad(top=self.o_d,
-                   bottom=-self.o_d,
+                   bottom=self.i_d,
                    left=position,
                    right=position + self.L,
                    line_color=bokeh_colors[0],
@@ -448,6 +448,15 @@ class ShaftElement(Element):
                    fill_color=bokeh_colors[2],
                    legend="Shaft"
                    )
+        bk_ax.quad(top=-self.i_d,
+                   bottom=-self.o_d,
+                   left=position,
+                   right=position + self.L,
+                   line_color=bokeh_colors[0],
+                   line_width=1,
+                   fill_alpha=0.5,
+                   fill_color=bokeh_colors[2],
+                   )
 
     @classmethod
     def section(


Test patch:
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index 20942387..b20d934e 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -1419,6 +1419,157 @@ def test_freq_response_w_force(rotor4):
     assert_allclose(mag[:4, :4], mag_exp_2_unb)
 
 
+def test_mesh_convergence(rotor3):
+    rotor3.convergence(n_eigval=0, err_max=1e-08)
+
+    assert_allclose(
+            len(rotor3.shaft_elements), 96, atol=0
+    )
+    assert_allclose(
+            rotor3.wn[0], 82.653037335, atol=1e-02
+    )
+    assert_allclose(
+            rotor3.shaft_elements[0].L, 0.015625, atol=1e-06
+    )
+    assert_allclose(
+            rotor3.disk_elements[0].n, 32, atol=0
+    )
+    assert_allclose(
+            rotor3.disk_elements[1].n, 64, atol=0
+    )
+    assert_allclose(
+            rotor3.bearing_seal_elements[0].n, 0, atol=0
+    )
+    assert_allclose(
+            rotor3.bearing_seal_elements[1].n, 96, atol=0
+    )
+    assert rotor3.error_arr[-1] <= 1e-08 * 100
+
+
+def test_static_analysis_rotor3(rotor3):
+    rotor3.static()
+
+    assert_almost_equal(
+        rotor3.disp_y,
+        ([-0.00061784, -0.00108199, -0.00143205, -0.00157464, -0.00147798,
+          -0.00115111, -0.00069521]),
+        decimal=6,
+    )
+    assert_almost_equal(
+        rotor3.Vx,
+        [0, -494.2745, -456.6791, -419.0836, -99.4925, -61.8971, -24.3016,
+         480.9807, 518.5762, 556.1716, 0],
+        decimal=3,
+    )
+    assert_almost_equal(
+        rotor3.Bm,
+        [0, -118.8692, -228.3395, -248.5132, -259.2881, -134.3434, 0],
+        decimal=3,
+    )
+
+
+@pytest.fixture
+def rotor5():
+    #  Rotor without damping with 10 shaft elements 2 disks and 2 bearings
+    i_d = 0
+    o_d = 0.05
+    n = 10
+    L = [0.25 for _ in range(n)]
+
+    shaft_elem = [
+        ShaftElement(
+            l, i_d, o_d, steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
+        )
+        for l in L
+    ]
+
+    disk0 = DiskElement.from_geometry(4, steel, 0.07, 0.05, 0.28)
+    disk1 = DiskElement.from_geometry(6, steel, 0.07, 0.05, 0.35)
+
+    stfx = 1e6
+    stfy = 1e6
+    bearing0 = BearingElement(2, kxx=stfx, kyy=stfy, cxx=0)
+    bearing1 = BearingElement(8, kxx=stfx, kyy=stfy, cxx=0)
+
+    return Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+
+
+def test_static_analysis_rotor5(rotor5):
+    rotor5.static()
+
+    assert_almost_equal(
+        rotor5.disp_y,
+        ([0.00026382, -0.00015021, -0.00056947, -0.00098566, -0.00130592,
+          -0.00143686, -0.00134669, -0.00104446, -0.00063136, -0.00021282,
+          0.0002005]),
+        decimal=6,
+    )
+    assert_almost_equal(
+        rotor5.Vx,
+        [0, 37.5954, 75.1908, -494.2745, -456.6791, -419.08368, -99.4925,
+         -61.8971, -24.3016, 480.9807, 518.5762, 556.1716, -75.1908, -37.5954,
+         0],
+        decimal=3,
+    )
+    assert_almost_equal(
+        rotor5.Bm,
+        [0, 4.6994, 18.7977, -100.0714, -209.5418, -229.7155, -240.4903,
+         -115.5457, 18.7977, 4.6994, 0],
+        decimal=3,
+    )
+
+
+@pytest.fixture
+def rotor6():
+    #  Overhung rotor without damping with 10 shaft elements 
+    #  2 disks and 2 bearings
+    i_d = 0
+    o_d = 0.05
+    n = 10
+    L = [0.25 for _ in range(n)]
+
+    shaft_elem = [
+        ShaftElement(
+            l, i_d, o_d, steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
+        )
+        for l in L
+    ]
+
+    disk0 = DiskElement.from_geometry(5, steel, 0.07, 0.05, 0.28)
+    disk1 = DiskElement.from_geometry(10, steel, 0.07, 0.05, 0.35)
+
+    stfx = 1e6
+    stfy = 1e6
+    bearing0 = BearingElement(2, kxx=stfx, kyy=stfy, cxx=0)
+    bearing1 = BearingElement(8, kxx=stfx, kyy=stfy, cxx=0)
+
+    return Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+
+
+def test_static_analysis_rotor6(rotor6):
+    rotor6.static()
+
+    assert_almost_equal(
+        rotor6.disp_y,
+        ([-2.58459386e-06, -8.83514255e-05, -1.79345202e-04, -2.82280131e-04,
+          -3.83451012e-04, -4.71329601e-04, -5.55753491e-04, -7.08192594e-04,
+          -1.02148266e-03, -1.55824814e-03, -2.22220183e-03]),
+        decimal=6,
+    )
+    assert_almost_equal(
+        rotor6.Vx,
+        [0, 37.5954, 75.1908, -104.1543, -66.5589, -28.9635, 8.6319, 328.2230,
+         365.8184, 403.4139, 441.0093, -580.4733, -542.8778, -505.2824, 0],
+        decimal=3,
+    )
+    assert_almost_equal(
+        rotor6.Bm,
+        [0, 4.6994, 18.7977, -2.5414, -14.4817, -17.0232, 69.7319, 165.8860,
+         271.4389, 131.0200, 0],
+        decimal=3,
+    )
+
+
 def test_save_load():
 
     a = rotor_example()