Repository: ross-rotordynamics/ross
Pull request #282 (2019-09-10 19:17:59)

Add matrix visualization
This PR adds a matrix visualization function. This function gives some visualization of a given matrix, displaying values on a heatmap. This can be useful in checking the global matrices’ assembling:

![image](https://user-images.githubusercontent.com/18506378/64643273-71669d80-d3e6-11e9-9239-0e4bb166d0f9.png)

There are also some corrections and formatting of the code.



Patch:
diff --git a/ross/__init__.py b/ross/__init__.py
index fed0633e..67676a1e 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -5,3 +5,4 @@ from .disk_element import *
 from .materials import *
 from .rotor_assembly import *
 from .shaft_element import *
+from .utils import visualize_matrix
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 5d958bb5..5462340d 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -147,7 +147,18 @@ class BearingElement(Element):
     """
 
     def __init__(
-        self, n, kxx, cxx, kyy=None, kxy=0, kyx=0, cyy=None, cxy=0, cyx=0, w=None, tag=None
+        self,
+        n,
+        kxx,
+        cxx,
+        kyy=None,
+        kxy=0,
+        kyx=0,
+        cyy=None,
+        cxy=0,
+        cyx=0,
+        w=None,
+        tag=None,
     ):
 
         args = ["kxx", "kyy", "kxy", "kyx", "cxx", "cyy", "cxy", "cyx"]
@@ -295,7 +306,7 @@ class BearingElement(Element):
         return bearing_elements
 
     def dof_mapping(self):
-        return dict(x0=0, y0=1)
+        return dict(x_0=0, y_0=1)
 
     def M(self):
         M = np.zeros((2, 2))
diff --git a/ross/disk_element.py b/ross/disk_element.py
index c8f05c6f..2cb0b813 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -100,7 +100,7 @@ class DiskElement(Element):
         return disk_elements
 
     def dof_mapping(self):
-        return dict(x0=0, y0=1, alpha0=2, beta0=3)
+        return dict(x_0=0, y_0=1, alpha_0=2, beta_0=3)
 
     def M(self):
         """
@@ -205,10 +205,10 @@ class DiskElement(Element):
         ax.add_patch(mpatches.Polygon(disk_points_l, facecolor=self.color))
 
         ax.add_patch(
-            mpatches.Circle(xy=(zpos, ypos*4), radius=step, color=self.color)
+            mpatches.Circle(xy=(zpos, ypos * 4), radius=step, color=self.color)
         )
         ax.add_patch(
-            mpatches.Circle(xy=(zpos, -ypos*4), radius=step, color=self.color)
+            mpatches.Circle(xy=(zpos, -ypos * 4), radius=step, color=self.color)
         )
 
     def bokeh_patch(self, position, bk_ax):
@@ -306,13 +306,12 @@ class DiskElement(Element):
             ("Polar Moment of Inertia :", "@IP"),
             ("Diametral Moment of Inertia :", "@ID"),
             ("Disk mass :", "@mass"),
-            ("Tag :", "@tag")
+            ("Tag :", "@tag"),
         ]
         hover.mode = "mouse"
 
         return hover
 
-
     @classmethod
     def from_geometry(cls, n, material, width, i_d, o_d, tag=None):
         """A disk element.
@@ -382,9 +381,15 @@ class DiskElement(Element):
         disk : list
             A list of disk objects.
         """
-        parameters = read_table_file(file, 'disk', sheet_name=sheet_name)
+        parameters = read_table_file(file, "disk", sheet_name=sheet_name)
         list_of_disks = []
-        for i in range(0, len(parameters['n'])):
-            list_of_disks.append(cls(n=parameters['n'][i], m=parameters['m'][i],
-                                     Id=parameters['Id'][i], Ip=parameters['Ip'][i]))
+        for i in range(0, len(parameters["n"])):
+            list_of_disks.append(
+                cls(
+                    n=parameters["n"][i],
+                    m=parameters["m"][i],
+                    Id=parameters["Id"][i],
+                    Ip=parameters["Ip"][i],
+                )
+            )
         return list_of_disks
diff --git a/ross/element.py b/ross/element.py
index 08343f0f..3cd8bb15 100644
--- a/ross/element.py
+++ b/ross/element.py
@@ -78,14 +78,14 @@ class Element(ABC):
         Example considering a shaft element:
         def dof_mapping(self):
             return dict(
-                x0=0,
-                y0=1,
-                alpha0=2,
-                beta0=3,
-                x1=4,
-                y1=5,
-                alpha1=6,
-                beta1=7,
+                x_0=0,
+                y_0=1,
+                alpha_0=2,
+                beta_0=3,
+                x_1=4,
+                y_1=5,
+                alpha_1=6,
+                beta_1=7,
             )
         """
         pass
@@ -101,11 +101,15 @@ class Element(ABC):
     def dof_global_index(self):
         """Get the global index for a element specific degree of freedom."""
         dof_mapping = self.dof_mapping()
-        dof_tuple = namedtuple("GlobalIndex", dof_mapping)
-
+        global_dof_mapping = {}
         for k, v in dof_mapping.items():
-            dof_mapping[k] = 4 * self.n + v
+            dof_letter, dof_number = k.split("_")
+            global_dof_mapping[dof_letter + "_" + str(int(dof_number) + self.n)] = v
+        dof_tuple = namedtuple("GlobalIndex", global_dof_mapping)
 
-        global_index = dof_tuple(**dof_mapping)
+        for k, v in global_dof_mapping.items():
+            global_dof_mapping[k] = 4 * self.n + v
+
+        global_index = dof_tuple(**global_dof_mapping)
 
         return global_index
diff --git a/ross/results.py b/ross/results.py
index 56ebdff6..82c6eb3b 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -10,14 +10,7 @@ from bokeh.layouts import gridplot, widgetbox
 from bokeh.plotting import figure, show
 from bokeh.transform import linear_cmap
 from bokeh.models.widgets import DataTable, TableColumn, NumberFormatter
-from bokeh.models import (
-        ColumnDataSource,
-        ColorBar,
-        Arrow,
-        NormalHead,
-        Label,
-        HoverTool
-)
+from bokeh.models import ColumnDataSource, ColorBar, Arrow, NormalHead, Label, HoverTool
 
 # set bokeh palette of colors
 bokeh_colors = bp.RdGy[11]
@@ -117,23 +110,27 @@ class CampbellResults:
                         label="Rotor speed",
                     )
 
-                    idx = np.argwhere(np.diff(np.sign(w_i - harm*speed_range_i))).flatten()
+                    idx = np.argwhere(
+                        np.diff(np.sign(w_i - harm * speed_range_i))
+                    ).flatten()
                     if len(idx) != 0:
                         idx = idx[0]
 
                         interpolated = interpolate.interp1d(
-                            x=[speed_range_i[idx], speed_range_i[idx+1]],
-                            y=[w_i[idx], w_i[idx+1]],
-                            kind="linear"
+                            x=[speed_range_i[idx], speed_range_i[idx + 1]],
+                            y=[w_i[idx], w_i[idx + 1]],
+                            kind="linear",
                         )
                         xnew = np.linspace(
                             speed_range_i[idx],
-                            speed_range_i[idx+1],
+                            speed_range_i[idx + 1],
                             num=20,
                             endpoint=True,
                         )
                         ynew = interpolated(xnew)
-                        idx = np.argwhere(np.diff(np.sign(ynew - harm*xnew))).flatten()
+                        idx = np.argwhere(
+                            np.diff(np.sign(ynew - harm * xnew))
+                        ).flatten()
 
                         ax.scatter(xnew[idx], ynew[idx], marker="X", s=30, c="g")
 
@@ -209,10 +206,10 @@ class CampbellResults:
         camp.yaxis.axis_label_text_font_size = "14pt"
 
         color_mapper = linear_cmap(
-                field_name="color",
-                palette=bp.viridis(256),
-                low=min(log_dec_map),
-                high=max(log_dec_map),
+            field_name="color",
+            palette=bp.viridis(256),
+            low=min(log_dec_map),
+            high=max(log_dec_map),
         )
 
         for mark, whirl_dir, legend in zip(
@@ -238,39 +235,41 @@ class CampbellResults:
                         muted_alpha=0.2,
                     )
 
-                    idx = np.argwhere(np.diff(np.sign(w_i - harm*speed_range_i))).flatten()
+                    idx = np.argwhere(
+                        np.diff(np.sign(w_i - harm * speed_range_i))
+                    ).flatten()
                     if len(idx) != 0:
                         idx = idx[0]
 
                         interpolated = interpolate.interp1d(
-                            x=[speed_range_i[idx], speed_range_i[idx+1]],
-                            y=[w_i[idx], w_i[idx+1]],
-                            kind="linear"
+                            x=[speed_range_i[idx], speed_range_i[idx + 1]],
+                            y=[w_i[idx], w_i[idx + 1]],
+                            kind="linear",
                         )
                         xnew = np.linspace(
                             speed_range_i[idx],
-                            speed_range_i[idx+1],
+                            speed_range_i[idx + 1],
                             num=30,
                             endpoint=True,
                         )
                         ynew = interpolated(xnew)
-                        idx = np.argwhere(np.diff(np.sign(ynew - harm*xnew))).flatten()
+                        idx = np.argwhere(
+                            np.diff(np.sign(ynew - harm * xnew))
+                        ).flatten()
 
-                        source = ColumnDataSource(
-                                dict(xnew=xnew[idx], ynew=ynew[idx])
-                        )
+                        source = ColumnDataSource(dict(xnew=xnew[idx], ynew=ynew[idx]))
                         camp.square(
-                                "xnew",
-                                "ynew",
-                                source=source,
-                                size=10,
-                                color=bokeh_colors[9],
-                                name="critspeed"
+                            "xnew",
+                            "ynew",
+                            source=source,
+                            size=10,
+                            color=bokeh_colors[9],
+                            name="critspeed",
                         )
                         hover = HoverTool(names=["critspeed"])
                         hover.tooltips = [
-                                ("Frequency :", "@xnew"),
-                                ("Critical Speed :", "@ynew")
+                            ("Frequency :", "@xnew"),
+                            ("Critical Speed :", "@ynew"),
                         ]
                         hover.mode = "mouse"
 
@@ -1200,7 +1199,7 @@ class ModeShapeResults:
         return fig, ax
 
 
-class StaticResults():
+class StaticResults:
     def __init__(
         self,
         disp_y,
@@ -1244,17 +1243,15 @@ class StaticResults():
         """
         source = ColumnDataSource(
             data=dict(
-                x0=self.nodes_pos,
-                y0=self.disp_y * 1000,
-                y1=[0] * len(self.nodes_pos)
+                x0=self.nodes_pos, y0=self.disp_y * 1000, y1=[0] * len(self.nodes_pos)
             )
         )
 
         TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select,hover"
         TOOLTIPS = [
-            ("Shaft lenght:", "@x0"),
+            ("Shaft lenght:", "@x_0"),
             ("Underformed:", "@y1"),
-            ("Displacement:", "@y0"),
+            ("Displacement:", "@y_0"),
         ]
 
         # create displacement plot
@@ -1271,11 +1268,11 @@ class StaticResults():
         disp_graph.yaxis.axis_label_text_font_size = "14pt"
 
         interpolated = interpolate.interp1d(
-            source.data["x0"], source.data["y0"], kind="cubic"
+            source.data["x_0"], source.data["y_0"], kind="cubic"
         )
         xnew = np.linspace(
-            source.data["x0"][0],
-            source.data["x0"][-1],
+            source.data["x_0"][0],
+            source.data["x_0"][-1],
             num=len(self.nodes_pos) * 20,
             endpoint=True,
         )
@@ -1284,23 +1281,23 @@ class StaticResults():
         auxsource = ColumnDataSource(data=dict(x0=xnew, y0=ynew, y1=[0] * len(xnew)))
 
         disp_graph.line(
-            "x0",
-            "y0",
+            "x_0",
+            "y_0",
             source=auxsource,
             legend="Deformed shaft",
             line_width=3,
             line_color=bokeh_colors[9],
         )
         disp_graph.circle(
-            "x0",
-            "y0",
+            "x_0",
+            "y_0",
             source=source,
             legend="Deformed shaft",
             size=8,
             fill_color=bokeh_colors[9],
         )
         disp_graph.line(
-            "x0",
+            "x_0",
             "y1",
             source=source,
             legend="underformed shaft",
@@ -1308,7 +1305,7 @@ class StaticResults():
             line_color=bokeh_colors[0],
         )
         disp_graph.circle(
-            "x0",
+            "x_0",
             "y1",
             source=source,
             legend="underformed shaft",
@@ -1328,12 +1325,12 @@ class StaticResults():
             title="Free-Body Diagram",
             x_axis_label="shaft lenght",
             x_range=[-0.1 * shaft_end, 1.1 * shaft_end],
-            y_range=[-3*y_start, 3*y_start],
+            y_range=[-3 * y_start, 3 * y_start],
         )
         FBD.yaxis.visible = False
         FBD.xaxis.axis_label_text_font_size = "14pt"
 
-        FBD.line("x0", "y1", source=source, line_width=5, line_color=bokeh_colors[0])
+        FBD.line("x_0", "y1", source=source, line_width=5, line_color=bokeh_colors[0])
 
         # FBD - plot arrows indicating shaft weight distribution
         text = str("%.1f" % sh_weight)
@@ -1345,8 +1342,8 @@ class StaticResults():
         )
 
         ini = self.nodes_pos[0]
-        fin = self.nodes_pos[-1]        
-        arrows_list = np.arange(ini, 1.01*fin, fin/5.)
+        fin = self.nodes_pos[-1]
+        arrows_list = np.arange(ini, 1.01 * fin, fin / 5.0)
         for node in arrows_list:
             FBD.add_layout(
                 Arrow(
@@ -1391,7 +1388,7 @@ class StaticResults():
                         line_color=bokeh_colors[0],
                     ),
                     x_start=self.nodes_pos[node],
-                    y_start=-2*y_start,
+                    y_start=-2 * y_start,
                     x_end=self.nodes_pos[node],
                     y_end=0,
                 )
@@ -1399,8 +1396,8 @@ class StaticResults():
             FBD.add_layout(
                 Label(
                     x=self.nodes_pos[node],
-                    y=-2*y_start,
-                    angle=np.pi/2,
+                    y=-2 * y_start,
+                    angle=np.pi / 2,
                     text="Fb = " + text + "N",
                     text_font_style="bold",
                     text_font_size="10pt",
@@ -1425,7 +1422,7 @@ class StaticResults():
                             line_color=bokeh_colors[0],
                         ),
                         x_start=self.nodes_pos[node],
-                        y_start=2*y_start,
+                        y_start=2 * y_start,
                         x_end=self.nodes_pos[node],
                         y_end=0,
                     )
@@ -1433,8 +1430,8 @@ class StaticResults():
                 FBD.add_layout(
                     Label(
                         x=self.nodes_pos[node],
-                        y=2*y_start,
-                        angle=np.pi/2,
+                        y=2 * y_start,
+                        angle=np.pi / 2,
                         text="Fd = " + text + "N",
                         text_font_style="bold",
                         text_font_size="10pt",
@@ -1496,7 +1493,9 @@ class StaticResults():
             interpolated_BM = interpolate.interp1d(
                 self.nodes_pos[i : i + 3], self.Bm[i : i + 3], kind="quadratic"
             )
-            xnew_BM = np.linspace(self.nodes_pos[i], self.nodes_pos[i + 2], num=42, endpoint=True)
+            xnew_BM = np.linspace(
+                self.nodes_pos[i], self.nodes_pos[i + 2], num=42, endpoint=True
+            )
 
             ynew_BM = interpolated_BM(xnew_BM)
             auxsource_BM = ColumnDataSource(data=dict(x=xnew_BM, y=ynew_BM))
@@ -1545,8 +1544,8 @@ class ConvergenceResults:
         )
 
         TOOLS = "pan,wheel_zoom,box_zoom,hover,reset,save,"
-        TOOLTIPS1 = [("Frequency:", "@y0"), ("Number of Elements", "@x0")]
-        TOOLTIPS2 = [("Relative Error:", "@y1"), ("Number of Elements", "@x0")]
+        TOOLTIPS1 = [("Frequency:", "@y_0"), ("Number of Elements", "@x_0")]
+        TOOLTIPS2 = [("Relative Error:", "@y1"), ("Number of Elements", "@x_0")]
 
         # create a new plot and add a renderer
         freq_arr = figure(
@@ -1561,8 +1560,8 @@ class ConvergenceResults:
         freq_arr.xaxis.axis_label_text_font_size = "14pt"
         freq_arr.yaxis.axis_label_text_font_size = "14pt"
 
-        freq_arr.line("x0", "y0", source=source, line_width=3, line_color="crimson")
-        freq_arr.circle("x0", "y0", source=source, size=8, fill_color="crimson")
+        freq_arr.line("x_0", "y_0", source=source, line_width=3, line_color="crimson")
+        freq_arr.circle("x_0", "y_0", source=source, size=8, fill_color="crimson")
 
         # create another new plot and add a renderer
         rel_error = figure(
@@ -1578,9 +1577,9 @@ class ConvergenceResults:
         rel_error.yaxis.axis_label_text_font_size = "14pt"
 
         rel_error.line(
-            "x0", "y1", source=source, line_width=3, line_color="darkslategray"
+            "x_0", "y1", source=source, line_width=3, line_color="darkslategray"
         )
-        rel_error.circle("x0", "y1", source=source, fill_color="darkslategray", size=8)
+        rel_error.circle("x_0", "y1", source=source, fill_color="darkslategray", size=8)
 
         # put the subplots in a gridplot
         plot = gridplot([[freq_arr, rel_error]])
@@ -1643,7 +1642,8 @@ class TimeResponseResults:
             tools="pan, box_zoom, wheel_zoom, reset, save",
             width=1200,
             height=900,
-            title="Response for node %s and degree of freedom %s" % (self.dof // 4, obs_dof),
+            title="Response for node %s and degree of freedom %s"
+            % (self.dof // 4, obs_dof),
             x_axis_label="Time (s)",
             y_axis_label="Amplitude (%s)" % amp,
         )
@@ -1652,10 +1652,7 @@ class TimeResponseResults:
 
         # bokeh plot - plot shaft centerline
         bk_ax.line(
-                self.t,
-                self.yout[:, self.dof],
-                line_width=3,
-                line_color=bokeh_colors[0]
+            self.t, self.yout[:, self.dof], line_width=3, line_color=bokeh_colors[0]
         )
 
         show(bk_ax)
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 9cee7b73..56a5a239 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -243,36 +243,44 @@ class ShaftElement(Element):
         shaft : list
             A list of shaft objects.
         """
-        parameters = read_table_file(file, 'shaft', sheet_name=sheet_name, sheet_type=sheet_type)
+        parameters = read_table_file(
+            file, "shaft", sheet_name=sheet_name, sheet_type=sheet_type
+        )
         list_of_shafts = []
         if sheet_type == "Model":
-            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],
-                                          material=parameters[parameters['material'][i]],
-                                          n=parameters['n'][i],
-                                          axial_force=parameters['axial_force'][i],
-                                          torque=parameters['torque'][i],
-                                          shear_effects=parameters['shear_effects'][i],
-                                          rotary_inertia=parameters['rotary_inertia'][i],
-                                          gyroscopic=parameters['gyroscopic'][i],
-                                          shear_method_calc=parameters['shear_method_calc'][i]
-                                          ))
+            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],
+                        material=parameters[parameters["material"][i]],
+                        n=parameters["n"][i],
+                        axial_force=parameters["axial_force"][i],
+                        torque=parameters["torque"][i],
+                        shear_effects=parameters["shear_effects"][i],
+                        rotary_inertia=parameters["rotary_inertia"][i],
+                        gyroscopic=parameters["gyroscopic"][i],
+                        shear_method_calc=parameters["shear_method_calc"][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],
-                                          material=parameters['material'][i],
-                                          n=parameters['n'][i],
-                                          axial_force=parameters['axial_force'][i],
-                                          torque=parameters['torque'][i],
-                                          shear_effects=parameters['shear_effects'][i],
-                                          rotary_inertia=parameters['rotary_inertia'][i],
-                                          gyroscopic=parameters['gyroscopic'][i],
-                                          shear_method_calc=parameters['shear_method_calc'][i]
-                                          ))
+            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],
+                        material=parameters["material"][i],
+                        n=parameters["n"][i],
+                        axial_force=parameters["axial_force"][i],
+                        torque=parameters["torque"][i],
+                        shear_effects=parameters["shear_effects"][i],
+                        rotary_inertia=parameters["rotary_inertia"][i],
+                        gyroscopic=parameters["gyroscopic"][i],
+                        shear_method_calc=parameters["shear_method_calc"][i],
+                    )
+                )
         return list_of_shafts
 
     @property
@@ -287,7 +295,9 @@ class ShaftElement(Element):
             self.n_r = value + 1
 
     def dof_mapping(self):
-        return dict(x0=0, y0=1, alpha0=2, beta0=3, x1=4, y1=5, alpha1=6, beta1=7)
+        return dict(
+            x_0=0, y_0=1, alpha_0=2, beta_0=3, x_1=4, y_1=5, alpha_1=6, beta_1=7
+        )
 
     def __repr__(self):
         return (
@@ -841,7 +851,12 @@ class ShaftTaperedElement(Element):
         a1 = 2 * np.pi * (roj * delta_ro - rij * delta_ri) / A_l
         a2 = np.pi * (roj ** 3 * delta_ro - rij ** 3 * delta_ri) / Ie_l
         b1 = np.pi * (delta_ro ** 2 - delta_ri ** 2) / A_l
-        b2 = 3 * np.pi * (roj ** 2 * delta_ro ** 2 - rij ** 2 * delta_ri ** 2) / (2 * Ie_l)
+        b2 = (
+            3
+            * np.pi
+            * (roj ** 2 * delta_ro ** 2 - rij ** 2 * delta_ri ** 2)
+            / (2 * Ie_l)
+        )
         gama = np.pi * (roj * delta_ro ** 3 - rij * delta_ri ** 3) / Ie_l
         delta = np.pi * (delta_ro ** 4 - delta_ri ** 4) / (4 * Ie_l)
 
@@ -953,7 +968,9 @@ class ShaftTaperedElement(Element):
             self.n_r = value + 1
 
     def dof_mapping(self):
-        return dict(x0=0, y0=1, alpha0=2, beta0=3, x1=4, y1=5, alpha1=6, beta1=7)
+        return dict(
+            x_0=0, y_0=1, alpha_0=2, beta_0=3, x_1=4, y_1=5, alpha_1=6, beta_1=7
+        )
 
     def __repr__(self):
         return (
diff --git a/ross/utils.py b/ross/utils.py
index f29ca979..f617e2c9 100644
--- a/ross/utils.py
+++ b/ross/utils.py
@@ -1,13 +1,19 @@
+import numpy as np
 import pandas as pd
 import xlrd
 import sys
 from ross.materials import Material
+from bokeh.io import show, output_notebook
+from bokeh.plotting import figure, ColumnDataSource
+from bokeh.palettes import Viridis256
+from bokeh.models import LinearColorMapper, LogColorMapper
 
 
 class DataNotFoundError(Exception):
     """
     An exception indicating that the data could not be found in the file.
     """
+
     pass
 
 
@@ -42,56 +48,66 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
     # Assign specific values to variables
     parameter_columns = {}
     optional_parameter_columns = {}
-    header_key_word = ''
+    header_key_word = ""
     default_dictionary = {}
     parameters = {}
-    if element == 'bearing':
-        header_key_word = 'kxx'
-        parameter_columns['kxx'] = ['kxx']
-        parameter_columns['cxx'] = ['cxx']
-        optional_parameter_columns['kyy'] = ['kyy']
-        optional_parameter_columns['kxy'] = ['kxy']
-        optional_parameter_columns['kyx'] = ['kyx']
-        optional_parameter_columns['cyy'] = ['cyy']
-        optional_parameter_columns['cxy'] = ['cxy']
-        optional_parameter_columns['cyx'] = ['cyx']
-        optional_parameter_columns['w'] = ['w', 'speed']
-        default_dictionary['kyy'] = None
-        default_dictionary['kxy'] = 0
-        default_dictionary['kyx'] = 0
-        default_dictionary['cyy'] = None
-        default_dictionary['cxy'] = 0
-        default_dictionary['cyx'] = 0
-        default_dictionary['w'] = None
-    elif element == 'shaft':
-        if sheet_type == 'Model':
-            header_key_word = 'od_left'
+    if element == "bearing":
+        header_key_word = "kxx"
+        parameter_columns["kxx"] = ["kxx"]
+        parameter_columns["cxx"] = ["cxx"]
+        optional_parameter_columns["kyy"] = ["kyy"]
+        optional_parameter_columns["kxy"] = ["kxy"]
+        optional_parameter_columns["kyx"] = ["kyx"]
+        optional_parameter_columns["cyy"] = ["cyy"]
+        optional_parameter_columns["cxy"] = ["cxy"]
+        optional_parameter_columns["cyx"] = ["cyx"]
+        optional_parameter_columns["w"] = ["w", "speed"]
+        default_dictionary["kyy"] = None
+        default_dictionary["kxy"] = 0
+        default_dictionary["kyx"] = 0
+        default_dictionary["cyy"] = None
+        default_dictionary["cxy"] = 0
+        default_dictionary["cyx"] = 0
+        default_dictionary["w"] = None
+    elif element == "shaft":
+        if sheet_type == "Model":
+            header_key_word = "od_left"
         else:
-            header_key_word = 'material'
-        parameter_columns['L'] = ['length']
-        parameter_columns['i_d'] = ['i_d', 'id', 'id_left']
-        parameter_columns['o_d'] = ['o_d', 'od', 'od_left']
-        parameter_columns['material'] = ['material', 'matnum']
-        optional_parameter_columns['n'] = ['n', 'elemnum']
-        optional_parameter_columns['axial_force'] = ['axial_force', 'axial force', 'axial']
-        optional_parameter_columns['torque'] = ['torque']
-        optional_parameter_columns['shear_effects'] = ['shear_effects', 'shear effects']
-        optional_parameter_columns['rotary_inertia'] = ['rotary_inertia', 'rotary inertia']
-        optional_parameter_columns['gyroscopic'] = ['gyroscopic']
-        optional_parameter_columns['shear_method_calc'] = ['shear_method_calc', 'shear method calc']
-        default_dictionary['n'] = None
-        default_dictionary['axial_force'] = 0
-        default_dictionary['torque'] = 0
-        default_dictionary['shear_effects'] = True
-        default_dictionary['rotary_inertia'] = True
-        default_dictionary['gyroscopic'] = True
-        default_dictionary['shear_method_calc'] = 'cowper'
-    elif element == 'disk':
-        header_key_word = 'ip'
-        parameter_columns['n'] = ['unnamed: 0', 'n']
-        parameter_columns['m'] = ['m', 'mass']
-        parameter_columns['Id'] = ['it', 'id']
-        parameter_columns['Ip'] = ['ip']
+            header_key_word = "material"
+        parameter_columns["L"] = ["length"]
+        parameter_columns["i_d"] = ["i_d", "id", "id_left"]
+        parameter_columns["o_d"] = ["o_d", "od", "od_left"]
+        parameter_columns["material"] = ["material", "matnum"]
+        optional_parameter_columns["n"] = ["n", "elemnum"]
+        optional_parameter_columns["axial_force"] = [
+            "axial_force",
+            "axial force",
+            "axial",
+        ]
+        optional_parameter_columns["torque"] = ["torque"]
+        optional_parameter_columns["shear_effects"] = ["shear_effects", "shear effects"]
+        optional_parameter_columns["rotary_inertia"] = [
+            "rotary_inertia",
+            "rotary inertia",
+        ]
+        optional_parameter_columns["gyroscopic"] = ["gyroscopic"]
+        optional_parameter_columns["shear_method_calc"] = [
+            "shear_method_calc",
+            "shear method calc",
+        ]
+        default_dictionary["n"] = None
+        default_dictionary["axial_force"] = 0
+        default_dictionary["torque"] = 0
+        default_dictionary["shear_effects"] = True
+        default_dictionary["rotary_inertia"] = True
+        default_dictionary["gyroscopic"] = True
+        default_dictionary["shear_method_calc"] = "cowper"
+    elif element == "disk":
+        header_key_word = "ip"
+        parameter_columns["n"] = ["unnamed: 0", "n"]
+        parameter_columns["m"] = ["m", "mass"]
+        parameter_columns["Id"] = ["it", "id"]
+        parameter_columns["Ip"] = ["ip"]
 
     # Find table header and define if conversion is needed
     header_index = -1
@@ -105,25 +121,27 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
                     if row[i].lower() == header_key_word:
                         header_index = index
                         header_found = True
-                if 'inches' in row[i].lower() or 'lbm' in row[i].lower():
+                if "inches" in row[i].lower() or "lbm" in row[i].lower():
                     convert_to_metric = True
-                if 'rpm' in row[i].lower():
+                if "rpm" in row[i].lower():
                     convert_to_rad_per_sec = True
                 if header_found and convert_to_metric and convert_to_rad_per_sec:
                     break
         if header_found and convert_to_metric:
             break
     if not header_found:
-        raise ValueError("Could not find the header. Make sure the table has a header "
-                         "containing the names of the columns. In the case of a " + element + ", "
-                         "there should be a column named " + header_key_word + ".")
+        raise ValueError(
+            "Could not find the header. Make sure the table has a header "
+            "containing the names of the columns. In the case of a " + element + ", "
+            "there should be a column named " + header_key_word + "."
+        )
 
     # Get specific data from the file
     new_materials = {}
-    if element == 'shaft' and sheet_type == 'Model':
+    if element == "shaft" and sheet_type == "Model":
         material_header_index = -1
         material_header_found = False
-        material_header_key_word = 'matno'
+        material_header_key_word = "matno"
         for index, row in df.iterrows():
             for i in range(0, row.size):
                 if isinstance(row[i], str):
@@ -134,20 +152,24 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
             if material_header_found:
                 break
         if not material_header_found:
-            raise ValueError("Could not find the header for the materials. Make sure the table has a header "
-                             "with the parameters for the materials that will be used. There should be a column "
-                             "named " + material_header_key_word + ".")
-        df_material = pd.read_excel(file, header=material_header_index, sheet_name=sheet_name)
+            raise ValueError(
+                "Could not find the header for the materials. Make sure the table has a header "
+                "with the parameters for the materials that will be used. There should be a column "
+                "named " + material_header_key_word + "."
+            )
+        df_material = pd.read_excel(
+            file, header=material_header_index, sheet_name=sheet_name
+        )
         material_name = []
         material_rho = []
         material_e = []
         material_g_s = []
         for index, row in df_material.iterrows():
-            if not pd.isna(row['matno']):
-                material_name.append(int(row['matno']))
-                material_rho.append(row['rhoa'])
-                material_e.append(row['ea'])
-                material_g_s.append(row['ga'])
+            if not pd.isna(row["matno"]):
+                material_name.append(int(row["matno"]))
+                material_rho.append(row["rhoa"])
+                material_e.append(row["ea"])
+                material_g_s.append(row["ga"])
             else:
                 break
         if convert_to_metric:
@@ -157,12 +179,12 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
                 material_g_s[i] = material_g_s[i] * 6894.757
         for i in range(0, len(material_name)):
             new_material = Material(
-                name='shaft_mat_' + str(material_name[i]),
+                name="shaft_mat_" + str(material_name[i]),
                 rho=material_rho[i],
                 E=material_e[i],
                 G_s=material_g_s[i],
             )
-            new_materials['shaft_mat_' + str(material_name[i])] = new_material
+            new_materials["shaft_mat_" + str(material_name[i])] = new_material
 
     df = pd.read_excel(file, header=header_index, sheet_name=sheet_name)
     df.columns = df.columns.str.lower()
@@ -172,15 +194,22 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
     last_data_row_found = False
     indexes_to_drop = []
     for index, row in df.iterrows():
-        if not first_data_row_found \
-                and (isinstance(row[header_key_word], int) or isinstance(row[header_key_word], float)) \
-                and not pd.isna(row[header_key_word]):
+        if (
+            not first_data_row_found
+            and (
+                isinstance(row[header_key_word], int)
+                or isinstance(row[header_key_word], float)
+            )
+            and not pd.isna(row[header_key_word])
+        ):
             first_data_row_found = True
         elif not first_data_row_found:
             indexes_to_drop.append(index)
         elif first_data_row_found and not last_data_row_found:
-            if (isinstance(row[header_key_word], int) or isinstance(row[header_key_word], float))\
-                    and not pd.isna(row[header_key_word]):
+            if (
+                isinstance(row[header_key_word], int)
+                or isinstance(row[header_key_word], float)
+            ) and not pd.isna(row[header_key_word]):
                 continue
             else:
                 last_data_row_found = True
@@ -188,14 +217,16 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
         elif last_data_row_found:
             indexes_to_drop.append(index)
     if not first_data_row_found:
-        raise DataNotFoundError("Could not find the data. Make sure you have at least one row containing "
-                                "data below the header.")
+        raise DataNotFoundError(
+            "Could not find the data. Make sure you have at least one row containing "
+            "data below the header."
+        )
     if len(indexes_to_drop) > 0:
         df = df.drop(indexes_to_drop)
 
     # Build parameters list
-    if element == 'bearing':
-        parameters['n'] = n
+    if element == "bearing":
+        parameters["n"] = n
     for key, value in parameter_columns.items():
         for name in value:
             try:
@@ -203,7 +234,9 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
                 break
             except KeyError:
                 if name == value[-1]:
-                    raise ValueError("Could not find a column with one of these names: " + str(value))
+                    raise ValueError(
+                        "Could not find a column with one of these names: " + str(value)
+                    )
                 continue
     for key, value in optional_parameter_columns.items():
         for name in value:
@@ -212,49 +245,165 @@ def read_table_file(file, element, sheet_name=0, n=0, sheet_type="Model"):
                 break
             except KeyError:
                 if name == value[-1]:
-                    parameters[key] = [default_dictionary[key]]*df.shape[0]
+                    parameters[key] = [default_dictionary[key]] * df.shape[0]
                 else:
                     continue
-    if element == 'shaft':
-        new_n = parameters['n']
+    if element == "shaft":
+        new_n = parameters["n"]
         for i in range(0, df.shape[0]):
             new_n[i] -= 1
-        parameters['n'] = new_n
-        if sheet_type == 'Model':
-            new_material = parameters['material']
+        parameters["n"] = new_n
+        if sheet_type == "Model":
+            new_material = parameters["material"]
             for i in range(0, df.shape[0]):
-                new_material[i] = 'shaft_mat_' + str(int(new_material[i]))
-            parameters['material'] = new_material
+                new_material[i] = "shaft_mat_" + str(int(new_material[i]))
+            parameters["material"] = new_material
     if convert_to_metric:
         for i in range(0, df.shape[0]):
-            if element == 'bearing':
-                parameters['kxx'][i] = parameters['kxx'][i] * 175.1268369864
-                parameters['cxx'][i] = parameters['cxx'][i] * 175.1268369864
-                parameters['kyy'][i] = parameters['kyy'][i] * 175.1268369864
-                parameters['kxy'][i] = parameters['kxy'][i] * 175.1268369864
-                parameters['kyx'][i] = parameters['kyx'][i] * 175.1268369864
-                parameters['cyy'][i] = parameters['cyy'][i] * 175.1268369864
-                parameters['cxy'][i] = parameters['cxy'][i] * 175.1268369864
-                parameters['cyx'][i] = parameters['cyx'][i] * 175.1268369864
-            if element == 'shaft':
-                parameters['L'][i] = parameters['L'][i] * 0.0254
-                parameters['i_d'][i] = parameters['i_d'][i] * 0.0254
-                parameters['o_d'][i] = parameters['o_d'][i] * 0.0254
-                parameters['axial_force'][i] = parameters['axial_force'][i] * 4.448221615255
-            elif element == 'disk':
-                parameters['m'][i] = parameters['m'][i] * 0.45359237
-                parameters['Id'][i] = parameters['Id'][i] * 0.0002926397
-                parameters['Ip'][i] = parameters['Ip'][i] * 0.0002926397
+            if element == "bearing":
+                parameters["kxx"][i] = parameters["kxx"][i] * 175.126_836_986_4
+                parameters["cxx"][i] = parameters["cxx"][i] * 175.126_836_986_4
+                parameters["kyy"][i] = parameters["kyy"][i] * 175.126_836_986_4
+                parameters["kxy"][i] = parameters["kxy"][i] * 175.126_836_986_4
+                parameters["kyx"][i] = parameters["kyx"][i] * 175.126_836_986_4
+                parameters["cyy"][i] = parameters["cyy"][i] * 175.126_836_986_4
+                parameters["cxy"][i] = parameters["cxy"][i] * 175.126_836_986_4
+                parameters["cyx"][i] = parameters["cyx"][i] * 175.126_836_986_4
+            if element == "shaft":
+                parameters["L"][i] = parameters["L"][i] * 0.0254
+                parameters["i_d"][i] = parameters["i_d"][i] * 0.0254
+                parameters["o_d"][i] = parameters["o_d"][i] * 0.0254
+                parameters["axial_force"][i] = (
+                    parameters["axial_force"][i] * 4.448_221_615_255
+                )
+            elif element == "disk":
+                parameters["m"][i] = parameters["m"][i] * 0.453_592_37
+                parameters["Id"][i] = parameters["Id"][i] * 0.000_292_639_7
+                parameters["Ip"][i] = parameters["Ip"][i] * 0.000_292_639_7
     if convert_to_rad_per_sec:
         for i in range(0, df.shape[0]):
-            if element == 'bearing':
-                parameters['w'][i] = parameters['w'][i] * 0.1047197551197
+            if element == "bearing":
+                parameters["w"][i] = parameters["w"][i] * 0.104_719_755_119_7
     parameters.update(new_materials)
     return parameters
 
 
+def visualize_matrix(rotor, matrix=None, frequency=None):
+    """Visualize global matrix.
+
+    This function gives some visualization of a given matrix, displaying
+    values on a heatmap.
+
+    Parameters
+    ----------
+    rotor: rs.Rotor
+    
+    matrix: str
+        String for the desired matrix.
+    frequency: float, optional
+        Excitation frequency. Defaults to rotor speed.
+    
+    Returns
+    -------
+    fig: bokeh.figure
+    """
+    if frequency is None:
+        frequency = rotor.w
+
+    A = np.zeros((rotor.ndof, rotor.ndof))
+    # E will store element's names and contributions to the global matrix
+    E = np.zeros((rotor.ndof, rotor.ndof), dtype=np.object)
+
+    M, N = E.shape
+    for i in range(M):
+        for j in range(N):
+            E[i, j] = []
+
+    for elm in rotor.elements:
+        g_dofs = elm.dof_global_index()
+        l_dofs = elm.dof_local_index()
+        try:
+            elm_matrix = getattr(elm, matrix)(frequency)
+        except TypeError:
+            elm_matrix = getattr(elm, matrix)()
+
+        A[np.ix_(g_dofs, g_dofs)] += elm_matrix
+
+        for l0, g0 in zip(l_dofs, g_dofs):
+            for l1, g1 in zip(l_dofs, g_dofs):
+                if elm_matrix[l0, l1] != 0:
+                    E[g0, g1].append(
+                        "\n"
+                        + elm.__class__.__name__
+                        + f"(n={elm.n})"
+                        + f": {elm_matrix[l0, l1]:.2E}"
+                    )
+
+    # list for dofs -> ['x0', 'y0', 'alpha0', 'beta0', 'x1'...]
+    dof_list = [0 for i in range(rotor.ndof)]
+
+    for elm in rotor.elements:
+        for k, v in elm.dof_global_index()._asdict().items():
+            dof_list[v] = k
+
+    data = {"row": [], "col": [], "value": [], "pos_value": [], "elements": []}
+    for i, dof_row in enumerate(dof_list):
+        for j, dof_col in enumerate(dof_list):
+            data["row"].append(i)
+            data["col"].append(j)
+            data["value"].append(A[i, j])
+            data["pos_value"].append(abs(A[i, j]))
+            data["elements"].append(E[i, j])
+
+    TOOLTIPS = """
+        <div style="width:150px;">
+            <span style="font-size: 10px;">Value: @value</span>
+            <br>
+            <span style="font-size: 10px;">Elements: </span>
+            <br>
+            <span style="font-size: 10px;">@elements</span>
+        </div>
+    """
+    fig = figure(
+        tools="hover",
+        tooltips=TOOLTIPS,
+        x_range=(-0.5, len(A) - 0.5),
+        y_range=(len(A) - 0.5, -0.5),
+        x_axis_location="above",
+    )
+    fig.plot_width = 500
+    fig.plot_height = 500
+    fig.grid.grid_line_color = None
+    fig.axis.axis_line_color = None
+    fig.axis.major_tick_line_color = None
+    fig.axis.major_label_text_font_size = "5pt"
+    fig.axis.major_label_standoff = 0
+    # p.xaxis.major_label_orientation = np.pi/3
+    fig.xaxis.ticker = [i for i in range(len(dof_list))]
+    fig.yaxis.ticker = [i for i in range(len(dof_list))]
 
+    dof_string_list = []
+    sub = str.maketrans("0123456789", "₀₁₂₃₄₅₆₇₈₉")
+    for d in dof_list:
+        d = d.replace("alpha", "α")
+        d = d.replace("beta", "β")
+        d = d.replace("_", "")
+        d = d.translate(sub)
+        dof_string_list.append(d)
 
+    fig.xaxis.major_label_overrides = {k: v for k, v in enumerate(dof_string_list)}
+    fig.yaxis.major_label_overrides = {k: v for k, v in enumerate(dof_string_list)}
 
+    mapper = LogColorMapper(palette=Viridis256, low=0, high=A.max())
 
+    fig.rect(
+        "row",
+        "col",
+        0.95,
+        0.95,
+        source=data,
+        fill_color={"field": "value", "transform": mapper},
+        line_color=None,
+    )
 
+    return fig


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index 2390fbb2..cf234de0 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -80,11 +80,11 @@ def bearing1():
 
 def test_index(bearing1):
     assert bearing1.dof_local_index()[0] == 0
-    assert bearing1.dof_local_index().x0 == 0
+    assert bearing1.dof_local_index().x_0 == 0
     assert bearing1.dof_local_index()[1] == 1
-    assert bearing1.dof_local_index().y0 == 1
-    assert bearing1.dof_global_index().x0 == 16
-    assert bearing1.dof_global_index().y0 == 17
+    assert bearing1.dof_local_index().y_0 == 1
+    assert bearing1.dof_global_index().x_4 == 16
+    assert bearing1.dof_global_index().y_4 == 17
     assert bearing1.dof_global_index()[-1] == 17
 
 
diff --git a/ross/tests/test_disk_element.py b/ross/tests/test_disk_element.py
index d84c31cb..aa6238e6 100644
--- a/ross/tests/test_disk_element.py
+++ b/ross/tests/test_disk_element.py
@@ -13,14 +13,14 @@ def disk():
 
 
 def test_index(disk):
-    assert disk.dof_local_index().x0 == 0
-    assert disk.dof_local_index().y0 == 1
-    assert disk.dof_local_index().alpha0 == 2
-    assert disk.dof_local_index().beta0 == 3
-    assert disk.dof_global_index().x0 == 0
-    assert disk.dof_global_index().y0 == 1
-    assert disk.dof_global_index().alpha0 == 2
-    assert disk.dof_global_index().beta0 == 3
+    assert disk.dof_local_index().x_0 == 0
+    assert disk.dof_local_index().y_0 == 1
+    assert disk.dof_local_index().alpha_0 == 2
+    assert disk.dof_local_index().beta_0 == 3
+    assert disk.dof_global_index().x_0 == 0
+    assert disk.dof_global_index().y_0 == 1
+    assert disk.dof_global_index().alpha_0 == 2
+    assert disk.dof_global_index().beta_0 == 3
 
 
 def test_mass_matrix_disk(disk):
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index 76183a51..7d5bb96a 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -22,14 +22,14 @@ def eb():
 
 
 def test_index(eb):
-    assert eb.dof_local_index().x0 == 0
-    assert eb.dof_local_index().y0 == 1
-    assert eb.dof_local_index().alpha0 == 2
-    assert eb.dof_local_index().beta0 == 3
-    assert eb.dof_global_index().x0 == 12
-    assert eb.dof_global_index().y0 == 13
-    assert eb.dof_global_index().alpha0 == 14
-    assert eb.dof_global_index().beta0 == 15
+    assert eb.dof_local_index().x_0 == 0
+    assert eb.dof_local_index().y_0 == 1
+    assert eb.dof_local_index().alpha_0 == 2
+    assert eb.dof_local_index().beta_0 == 3
+    assert eb.dof_global_index().x_3 == 12
+    assert eb.dof_global_index().y_3 == 13
+    assert eb.dof_global_index().alpha_3 == 14
+    assert eb.dof_global_index().beta_3 == 15
 
 
 def test_parameters_eb(eb):
@@ -137,7 +137,9 @@ def test_from_table():
         os.path.dirname(os.path.realpath(__file__)) + "/data/shaft_us.xls",
         os.path.dirname(os.path.realpath(__file__)) + "/data/shaft_si.xls",
     ]:
-        shaft = ShaftElement.from_table(shaft_file, sheet_type="Model", sheet_name="Model")
+        shaft = ShaftElement.from_table(
+            shaft_file, sheet_type="Model", sheet_name="Model"
+        )
         el0 = shaft[0]
         assert el0.n == 0
         assert_allclose(el0.i_d, 0.1409954)
@@ -160,25 +162,42 @@ def test_from_table():
 def tap_tim():
     #  Timoshenko element
     L = 0.4
-    i_d_l = 0.
-    i_d_r = 0.
+    i_d_l = 0.0
+    i_d_r = 0.0
     o_d_l = 0.25
     o_d_r = 0.10
 
     return ShaftTaperedElement(
-        L, i_d_l, i_d_r, o_d_l, o_d_r, steel, shear_effects=True, rotary_inertia=True, n=3
+        L,
+        i_d_l,
+        i_d_r,
+        o_d_l,
+        o_d_r,
+        steel,
+        shear_effects=True,
+        rotary_inertia=True,
+        n=3,
     )
 
 
 def test_tapered_index(tap_tim):
-    assert tap_tim.dof_local_index().x0 == 0
-    assert tap_tim.dof_local_index().y0 == 1
-    assert tap_tim.dof_local_index().alpha0 == 2
-    assert tap_tim.dof_local_index().beta0 == 3
-    assert tap_tim.dof_global_index().x0 == 12
-    assert tap_tim.dof_global_index().y0 == 13
-    assert tap_tim.dof_global_index().alpha0 == 14
-    assert tap_tim.dof_global_index().beta0 == 15
+    assert tap_tim.dof_local_index().x_0 == 0
+    assert tap_tim.dof_local_index().y_0 == 1
+    assert tap_tim.dof_local_index().alpha_0 == 2
+    assert tap_tim.dof_local_index().beta_0 == 3
+    assert tap_tim.dof_local_index().x_1 == 4
+    assert tap_tim.dof_local_index().y_1 == 5
+    assert tap_tim.dof_local_index().alpha_1 == 6
+    assert tap_tim.dof_local_index().beta_1 == 7
+    print(tap_tim.dof_global_index())
+    assert tap_tim.dof_global_index().x_3 == 12
+    assert tap_tim.dof_global_index().y_3 == 13
+    assert tap_tim.dof_global_index().alpha_3 == 14
+    assert tap_tim.dof_global_index().beta_3 == 15
+    assert tap_tim.dof_global_index().x_4 == 16
+    assert tap_tim.dof_global_index().y_4 == 17
+    assert tap_tim.dof_global_index().alpha_4 == 18
+    assert tap_tim.dof_global_index().beta_4 == 19
 
 
 def test_parameters_tap_tim(tap_tim):
@@ -194,7 +213,7 @@ def test_parameters_tap_tim(tap_tim):
     assert tap_tim.material.rho == 7810
     assert_almost_equal(tap_tim.material.Poisson, 0.29926108)
     assert_almost_equal(tap_tim.A, 0.024052818754)
-    assert_almost_equal(tap_tim.Ie * 1e5,  4.60385984)
+    assert_almost_equal(tap_tim.Ie * 1e5, 4.60385984)
     assert_almost_equal(tap_tim.phi, 0.4208816002)
 
 
@@ -247,8 +266,8 @@ def test_gyroscopic_matrix_tap_tim(tap_tim):
 def tap2():
     #  Timoshenko element
     L = 0.4
-    i_d_l = 0.
-    i_d_r = 0.
+    i_d_l = 0.0
+    i_d_r = 0.0
     o_d_l = 0.25
     o_d_r = 0.25
 
@@ -263,9 +282,7 @@ def tim2():
     le_ = 0.4
     i_d_ = 0
     o_d_ = 0.25
-    return ShaftElement(
-            le_, i_d_, o_d_, steel, rotary_inertia=True, shear_effects=True
-    )
+    return ShaftElement(le_, i_d_, o_d_, steel, rotary_inertia=True, shear_effects=True)
 
 
 def test_match_mass_matrix(tap2, tim2):