Repository: ross-rotordynamics/ross (python)

Issue:
Include bearing housing dynamics
It is important to be able to include in the modelling the stiffness and inertia of the bearing housing, external supports, squeeze filme dampers or any other kind of support dynamics other than the actual oil film.

We usually have mass and stiffness parameters of these supports. So a easy way to include this in the model is create a point mass element and connect it to the rotor via the oil film sitffness and to the ground via the support stiffness.

It would be best if the software allowed the user to create multiple point mass elements and connect them to any other two nodes, that might be a rotor node, another point mass node or the ground.
This way we can model the full system including, for example, the squeeze filme damper, the bearing housing and baseplate, when needed.

Pull request description:
Add rotor link
The purpose of this PR is to add the possibility of simulating bearing housing and supports, as requested on #245.

It is still a working in progress.

- [x] Create class `PointMass`;
- [x] Add an additional argument `n_link` to bearing;
- [x] Change rotor assembly to check for links when assembling global matrices;
- [x] Add patches for drawing.

Closes #245.

Patch:
diff --git a/ross/__init__.py b/ross/__init__.py
index 67676a1e..b9c92489 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -6,3 +6,4 @@ from .materials import *
 from .rotor_assembly import *
 from .shaft_element import *
 from .utils import visualize_matrix
+from .point_mass import *
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 2c7be88f..60cf2d9a 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -7,6 +7,7 @@ import matplotlib.lines as mlines
 import numpy as np
 import scipy.interpolate as interpolate
 
+from collections import namedtuple
 from ross.utils import read_table_file
 from ross.element import Element
 from ross.fluid_flow import fluid_flow as flow
@@ -125,9 +126,13 @@ class BearingElement(Element):
         (defaults to 0)
     w: array, optional
         Array with the speeds (rad/s).
-     tag : str, optional
-         A tag to name the element
-         Default is None
+    tag: str, optional
+        A tag to name the element
+        Default is None.
+    n_link: int, optional
+        Node to which the bearing will connect. If None the bearing is
+        connected to ground.
+        Default is None.
     Examples
     --------
     >>> # A bearing element located in the first rotor node, with these
@@ -159,6 +164,7 @@ class BearingElement(Element):
         cyx=0,
         w=None,
         tag=None,
+        n_link=None,
     ):
 
         args = ["kxx", "kyy", "kxy", "kyx", "cxx", "cyy", "cxy", "cyx"]
@@ -199,6 +205,7 @@ class BearingElement(Element):
             setattr(self, k, v)
 
         self.n = n
+        self.n_link = n_link
         self.n_l = n
         self.n_r = n
 
@@ -259,6 +266,7 @@ class BearingElement(Element):
             "cxy",
             "cyx",
             "w",
+            "n",
         ]
         if isinstance(other, self.__class__):
             return all(
@@ -270,6 +278,9 @@ class BearingElement(Element):
             )
         return False
 
+    def __hash__(self):
+        return hash(self.tag)
+
     def save(self, file_name):
         """Saves a bearing element in a toml format. It works as an auxiliary function of
         the save function in the Rotor class.
@@ -387,6 +398,19 @@ class BearingElement(Element):
         """
         return dict(x_0=0, y_0=1)
 
+    def dof_global_index(self):
+        """Get the global index for a element specific degree of freedom."""
+        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}"] = 4 * self.n_link
+            global_index[f"y_{self.n_link}"] = 4 * self.n_link + 1
+            dof_tuple = namedtuple("GlobalIndex", global_index)
+            global_index = dof_tuple(**global_index)
+
+        return global_index
+
     def M(self):
         """Returns the mass matrix.
         Parameters
@@ -403,7 +427,7 @@ class BearingElement(Element):
         array([[0., 0.],
                [0., 0.]])
         """
-        M = np.zeros((2, 2))
+        M = np.zeros_like(self.K(0))
 
         return M
 
@@ -432,6 +456,12 @@ class BearingElement(Element):
 
         K = np.array([[kxx, kxy], [kyx, kyy]])
 
+        if self.n_link is not None:
+            # fmt: off
+            K = np.vstack((np.hstack([K, -K]),
+                           np.hstack([-K, K])))
+            # fmt: on
+
         return K
 
     def C(self, w):
@@ -459,6 +489,12 @@ class BearingElement(Element):
 
         C = np.array([[cxx, cxy], [cyx, cyy]])
 
+        if self.n_link is not None:
+            # fmt: off
+            C = np.vstack((np.hstack([C, -C]),
+                           np.hstack([-C, C])))
+            # fmt: on
+
         return C
 
     def G(self):
@@ -477,7 +513,7 @@ class BearingElement(Element):
         array([[0., 0.],
                [0., 0.]])
         """
-        G = np.zeros((2, 2))
+        G = np.zeros_like(self.K(0))
 
         return G
 
@@ -649,19 +685,20 @@ class BearingElement(Element):
         bk_ax.line(x=x_top, y=yu_top, legend="Bearing", **kwargs)
 
         # plot ground
-        zl_g = [zs0 - step, zs1 + step]
-        yl_g = [yl_top[0], yl_top[0]]
-        yu_g = [-y for y in yl_g]
-        bk_ax.line(x=zl_g, y=yl_g, **kwargs)
-        bk_ax.line(x=zl_g, y=yu_g, **kwargs)
-
-        step2 = (zl_g[1] - zl_g[0]) / n
-        for i in range(n + 1):
-            zl_g2 = [(zs0 - step) + step2 * (i), (zs0 - step) + step2 * (i + 1)]
-            yl_g2 = [yl_g[0], 1.1 * yl_g[0]]
-            yu_g2 = [-y for y in yl_g2]
-            bk_ax.line(x=zl_g2, y=yl_g2, **kwargs)
-            bk_ax.line(x=zl_g2, y=yu_g2, **kwargs)
+        if self.n_link is None:
+            zl_g = [zs0 - step, zs1 + step]
+            yl_g = [yl_top[0], yl_top[0]]
+            yu_g = [-y for y in yl_g]
+            bk_ax.line(x=zl_g, y=yl_g, **kwargs)
+            bk_ax.line(x=zl_g, y=yu_g, **kwargs)
+
+            step2 = (zl_g[1] - zl_g[0]) / n
+            for i in range(n + 1):
+                zl_g2 = [(zs0 - step) + step2 * (i), (zs0 - step) + step2 * (i + 1)]
+                yl_g2 = [yl_g[0], 1.1 * yl_g[0]]
+                yu_g2 = [-y for y in yl_g2]
+                bk_ax.line(x=zl_g2, y=yl_g2, **kwargs)
+                bk_ax.line(x=zl_g2, y=yu_g2, **kwargs)
 
         # plot spring
         z_spring = np.array([zs0, zs0, zs0, zs0])
@@ -957,6 +994,7 @@ class SealElement(BearingElement):
     >>> seal.C(w) # doctest: +ELLIPSIS
     array([[[200., 200., ...
     """
+
     def __init__(
         self,
         n,
@@ -1119,6 +1157,3 @@ def seal_example():
     w = np.linspace(0, 200, 11)
     seal = SealElement(n=0, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, w=w)
     return seal
-
-
-
diff --git a/ross/disk_element.py b/ross/disk_element.py
index 2cb0b813..91a7bec9 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -77,6 +77,9 @@ class DiskElement(Element):
             f"n={self.n}, tag={self.tag!r})"
         )
 
+    def __hash__(self):
+        return hash(self.tag)
+
     def save(self, file_name):
         data = self.load_data(file_name)
         data["DiskElement"][str(self.n)] = {
diff --git a/ross/element.py b/ross/element.py
index 3cd8bb15..4a4a82c6 100644
--- a/ross/element.py
+++ b/ross/element.py
@@ -40,7 +40,7 @@ class Element(ABC):
 
     @abstractmethod
     def G(self):
-        """Giroscopic matrix."""
+        """Gyroscopic matrix."""
         pass
 
     def summary(self):
@@ -94,9 +94,9 @@ class Element(ABC):
         """Get the local index for a element specific degree of freedom."""
         dof_mapping = self.dof_mapping()
         dof_tuple = namedtuple("LocalIndex", dof_mapping)
-        global_index = dof_tuple(**dof_mapping)
+        local_index = dof_tuple(**dof_mapping)
 
-        return global_index
+        return local_index
 
     def dof_global_index(self):
         """Get the global index for a element specific degree of freedom."""
diff --git a/ross/point_mass.py b/ross/point_mass.py
new file mode 100644
index 00000000..44a57fa4
--- /dev/null
+++ b/ross/point_mass.py
@@ -0,0 +1,84 @@
+"""Point mass module.
+
+This module defines the PointMass class which will be used to link elements.
+"""
+import numpy as np
+from ross.element import Element
+
+__all__ = ["PointMass"]
+
+
+class PointMass(Element):
+    """A point mass element.
+
+    This class will create a point mass element.
+    This element can be used to link other elements in the analysis.
+    The mass provided to the element can be different on the x and y direction
+    (e.g. different support inertia for x and y directions).
+
+    Parameters
+    ----------
+    n: int
+        Node which the bearing will be located in.
+    m: float, optional
+        Mass for the element.
+    mx: float, optional
+        Mass for the element on the x direction.
+    my: float, optional
+        Mass for the element on the y direction.
+
+    Examples
+    --------
+    >>> p0 = PointMass(n=0, m=2)
+    >>> p0.M()
+    array([[2, 0],
+           [0, 2]])
+    >>> p1 = PointMass(n=0, mx=2, my=3)
+    >>> p1.M()
+    array([[2, 0],
+           [0, 3]])
+    """
+
+    def __init__(self, n=None, m=None, mx=None, my=None, tag=None):
+        self.n = n
+        self.m = m
+
+        if mx is None and my is None:
+            mx = m
+            my = m
+
+        self.mx = mx
+        self.my = my
+
+        if tag is None:
+            tag = self.__class__.__name__ + " " + str(self.n)
+        self.tag = tag
+
+    def __hash__(self):
+        return hash(self.tag)
+
+    def M(self):
+        """Mass matrix."""
+        mx = self.mx
+        my = self.my
+        # fmt: off
+        M = np.array([[mx, 0],
+                      [0, my]])
+        # fmt: on
+
+        return M
+
+    def C(self):
+        """Damping coefficients matrix."""
+        return np.zeros((2, 2))
+
+    def K(self):
+        """Stiffness coefficients matrix."""
+        return np.zeros((2, 2))
+
+    def G(self):
+        """Gyroscopic matrix."""
+        return np.zeros((2, 2))
+
+    def dof_mapping(self):
+        return dict(x_0=0, y_0=1)
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 454b031d..67507760 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -4,6 +4,7 @@ import warnings
 from collections import Iterable
 from copy import copy, deepcopy
 from pathlib import Path
+from itertools import cycle
 
 import bokeh.palettes as bp
 import matplotlib as mpl
@@ -74,6 +75,8 @@ class Rotor(object):
         List with the disk elements
     bearing_seal_elements : list
         List with the bearing elements
+    point_mass_elements: list
+        List with the point mass elements
     w : float, optional
         Rotor speed. Defaults to 0.
     sparse : bool, optional
@@ -131,6 +134,7 @@ class Rotor(object):
         shaft_elements,
         disk_elements=None,
         bearing_seal_elements=None,
+        point_mass_elements=None,
         w=0,
         sparse=True,
         n_eigen=12,
@@ -179,11 +183,15 @@ class Rotor(object):
         for i, sh in enumerate(shaft_elements):
             if sh.n is None:
                 sh.n = i
+            if sh.tag is None:
+                sh.tag = sh.__class__.__name__ + " " + str(i) + f" i_d={sh.i_d}"
 
         if disk_elements is None:
             disk_elements = []
         if bearing_seal_elements is None:
             bearing_seal_elements = []
+        if point_mass_elements is None:
+            point_mass_elements = []
 
         for i, disk in enumerate(disk_elements):
             if disk.tag is None:
@@ -198,10 +206,16 @@ class Rotor(object):
         self.shaft_elements = shaft_elements
         self.bearing_seal_elements = bearing_seal_elements
         self.disk_elements = disk_elements
+        self.point_mass_elements = point_mass_elements
         self.elements = [
             el
             for el in flatten(
-                [self.shaft_elements, self.disk_elements, self.bearing_seal_elements]
+                [
+                    self.shaft_elements,
+                    self.disk_elements,
+                    self.bearing_seal_elements,
+                    self.point_mass_elements,
+                ]
             )
         ]
 
@@ -211,9 +225,11 @@ class Rotor(object):
         columns = [
             "type",
             "n",
+            "n_link",
             "L",
             "node_pos",
             "node_pos_r",
+            "y_pos",
             "i_d",
             "o_d",
             "i_d_r",
@@ -228,11 +244,20 @@ class Rotor(object):
         df_shaft = pd.DataFrame([el.summary() for el in self.shaft_elements])
         df_disks = pd.DataFrame([el.summary() for el in self.disk_elements])
         df_bearings = pd.DataFrame(
-                [el.summary() for el in self.bearing_seal_elements if (
-                        el.__class__.__name__ == 'BearingElement')])
+            [
+                el.summary()
+                for el in self.bearing_seal_elements
+                if (el.__class__.__name__ == "BearingElement")
+            ]
+        )
         df_seals = pd.DataFrame(
-                [el.summary() for el in self.bearing_seal_elements if (
-                        el.__class__.__name__ == 'SealElement')])
+            [
+                el.summary()
+                for el in self.bearing_seal_elements
+                if (el.__class__.__name__ == "SealElement")
+            ]
+        )
+        df_point_mass = pd.DataFrame([el.summary() for el in self.point_mass_elements])
 
         nodes_pos_l = np.zeros(len(df_shaft.n_l))
         nodes_pos_r = np.zeros(len(df_shaft.n_l))
@@ -251,21 +276,26 @@ class Rotor(object):
         df_shaft["nodes_pos_l"] = nodes_pos_l
         df_shaft["nodes_pos_r"] = nodes_pos_r
 
-        df = pd.concat([df_shaft, df_disks, df_bearings, df_seals], sort=True)
+        df = pd.concat(
+            [df_shaft, df_disks, df_bearings, df_point_mass, df_seals], sort=True
+        )
         df = df.sort_values(by="n_l")
         df = df.reset_index(drop=True)
 
         self.df_disks = df_disks
         self.df_bearings = df_bearings
         self.df_shaft = df_shaft
+        self.df_point_mass = df_point_mass
         self.df_seals = df_seals
 
         # check consistence for disks and bearings location
-        if df.n_l.max() > df_shaft.n_r.max():
+        if len(df_point_mass) > 0:
+            max_loc_point_mass = df_point_mass.n.max()
+        else:
+            max_loc_point_mass = 0
+        max_location = max(df_shaft.n_r.max(), max_loc_point_mass)
+        if df.n_l.max() > max_location:
             raise ValueError("Trying to set disk or bearing outside shaft")
-
-        self.df = df
-
         # nodes axial position and diameter
         nodes_pos = list(df_shaft.groupby("n_l")["nodes_pos_l"].max())
         nodes_pos.append(df_shaft["nodes_pos_r"].iloc[-1])
@@ -292,7 +322,7 @@ class Rotor(object):
         self.m_shaft = np.sum([sh_el.m for sh_el in self.shaft_elements])
         self.m = self.m_disks + self.m_shaft
 
-        # values for evalues and evectors will be calculated by self._calc_system
+        # values for evalues and evectors will be calculated by self.run_modal
         self.evalues = None
         self.evectors = None
         self.wn = None
@@ -302,14 +332,77 @@ class Rotor(object):
         self._v0 = None  # used to call eigs
 
         # number of dofs
-        self.ndof = 4 * max([el.n for el in shaft_elements]) + 8
+        self.ndof = (
+            4 * max([el.n for el in shaft_elements])
+            + 8
+            + 2 * len([el for el in point_mass_elements])
+        )
 
         #  values for static analysis will be calculated by def static
         self.Vx = None
         self.Bm = None
         self.disp_y = None
 
-        #  diameter at node position
+        # define positions for disks
+        for disk in disk_elements:
+            z_pos = nodes_pos[disk.n]
+            y_pos = nodes_o_d[disk.n]
+            df.loc[df.tag == disk.tag, "nodes_pos_l"] = z_pos
+            df.loc[df.tag == disk.tag, "nodes_pos_r"] = z_pos
+            df.loc[df.tag == disk.tag, "y_pos"] = y_pos
+
+        # define positions for bearings
+        # check if there are bearings without location
+        bearings_no_zloc = {
+            b
+            for b in bearing_seal_elements
+            if pd.isna(df.loc[df.tag == b.tag, "nodes_pos_l"]).all()
+        }
+        # cycle while there are bearings without a z location
+        for b in cycle(bearing_seal_elements):
+            if bearings_no_zloc:
+                if b in bearings_no_zloc:
+                    # first check if b.n is on list, if not, check for n_link
+                    node_l = df.loc[(df.n_l == b.n) & (df.tag != b.tag), "nodes_pos_l"]
+                    node_r = df.loc[(df.n_r == b.n) & (df.tag != b.tag), "nodes_pos_r"]
+                    if len(node_l) == 0 and len(node_r) == 0:
+                        node_l = df.loc[
+                            (df.n_link == b.n) & (df.tag != b.tag), "nodes_pos_l"
+                        ]
+                        node_r = node_l
+                    if len(node_l):
+                        df.loc[df.tag == b.tag, "nodes_pos_l"] = node_l.values[0]
+                        df.loc[df.tag == b.tag, "nodes_pos_r"] = node_l.values[0]
+                        bearings_no_zloc.discard(b)
+                    elif len(node_r):
+                        df.loc[df.tag == b.tag, "nodes_pos_l"] = node_r.values[0]
+                        df.loc[df.tag == b.tag, "nodes_pos_r"] = node_r.values[0]
+                        bearings_no_zloc.discard(b)
+            else:
+                break
+
+        dfb = df[df.type == "BearingElement"]
+        z_positions = [pos for pos in dfb["nodes_pos_l"]]
+        for z_pos in z_positions:
+            dfb_z_pos = dfb[dfb.nodes_pos_l == z_pos]
+            dfb_z_pos = dfb_z_pos.sort_values(by="n_l")
+            y_pos = nodes_o_d[int(dfb_z_pos.iloc[0]["n_l"])] / 2
+            mean_od = np.mean(nodes_o_d)
+            for t in dfb_z_pos.tag:
+                df.loc[df.tag == t, "y_pos"] = y_pos
+                y_pos += mean_od
+
+        # define position for point mass elements
+        dfb = df[df.type == "BearingElement"]
+        for p in point_mass_elements:
+            z_pos = dfb[dfb.n_l == p.n]["nodes_pos_l"].values[0]
+            y_pos = dfb[dfb.n_l == p.n]["y_pos"].values[0]
+            df.loc[df.tag == p.tag, "nodes_pos_l"] = z_pos
+            df.loc[df.tag == p.tag, "nodes_pos_r"] = z_pos
+            df.loc[df.tag == p.tag, "y_pos"] = y_pos
+
+        self.df = df
+
         self.run_modal()
 
     def __eq__(self, other):
@@ -476,9 +569,7 @@ class Rotor(object):
 
         for elm in self.elements:
             dofs = elm.dof_global_index()
-            n0 = dofs[0]
-            n1 = dofs[-1] + 1  # +1 to include this dof in the slice
-            M0[n0:n1, n0:n1] += elm.M()
+            M0[np.ix_(dofs, dofs)] += elm.M()
 
         return M0
 
@@ -503,19 +594,17 @@ class Rotor(object):
                [ 0., -6.,  1.,  0.],
                [ 6.,  0.,  0.,  1.]])
         """
-        K0 = np.zeros((self.ndof, self.ndof))
-
         if frequency is None:
             frequency = self.w
 
+        K0 = np.zeros((self.ndof, self.ndof))
+
         for elm in self.elements:
             dofs = elm.dof_global_index()
-            n0 = dofs[0]
-            n1 = dofs[-1] + 1  # +1 to include this dof in the slice
             try:
-                K0[n0:n1, n0:n1] += elm.K(frequency)
+                K0[np.ix_(dofs, dofs)] += elm.K(frequency)
             except TypeError:
-                K0[n0:n1, n0:n1] += elm.K()
+                K0[np.ix_(dofs, dofs)] += elm.K()
 
         return K0
 
@@ -540,19 +629,18 @@ class Rotor(object):
                [0., 0., 0., 0.],
                [0., 0., 0., 0.]])
         """
-        C0 = np.zeros((self.ndof, self.ndof))
-
         if frequency is None:
             frequency = self.w
 
+        C0 = np.zeros((self.ndof, self.ndof))
+
         for elm in self.elements:
             dofs = elm.dof_global_index()
-            n0 = dofs[0]
-            n1 = dofs[-1] + 1  # +1 to include this dof in the slice
+
             try:
-                C0[n0:n1, n0:n1] += elm.C(frequency)
+                C0[np.ix_(dofs, dofs)] += elm.C(frequency)
             except TypeError:
-                C0[n0:n1, n0:n1] += elm.C()
+                C0[np.ix_(dofs, dofs)] += elm.C()
 
         return C0
 
@@ -576,9 +664,7 @@ class Rotor(object):
 
         for elm in self.elements:
             dofs = elm.dof_global_index()
-            n0 = dofs[0]
-            n1 = dofs[-1] + 1  # +1 to include this dof in the slice
-            G0[n0:n1, n0:n1] += elm.G()
+            G0[np.ix_(dofs, dofs)] += elm.G()
 
         return G0
 
@@ -1234,7 +1320,7 @@ class Rotor(object):
             position = (self.nodes_pos[disk.n], self.nodes_o_d[disk.n] / 2)
             disk.patch(position, ax)
 
-        mean_od = np.mean(self.nodes_o_d)
+        mean_od = np.mean(self.nodes_o_d) / 2
         # plot bearings
         for bearing in self.bearing_seal_elements:
             position = (self.nodes_pos[bearing.n], self.nodes_o_d[bearing.n] / 2)
@@ -1337,7 +1423,9 @@ class Rotor(object):
         mean_od = np.mean(self.nodes_o_d) / 2
         # plot bearings
         for bearing in self.bearing_seal_elements:
-            position = (self.nodes_pos[bearing.n], self.nodes_o_d[bearing.n] / 2)
+            z_pos = self.df[self.df.tag == bearing.tag]["nodes_pos_l"].values[0]
+            y_pos = self.df[self.df.tag == bearing.tag]["y_pos"].values[0]
+            position = (z_pos, y_pos)
             bearing.bokeh_patch(position, mean_od, bk_ax)
 
         show(bk_ax)
@@ -1377,12 +1465,10 @@ class Rotor(object):
         """
         if plot_type == "matplotlib":
             return self._plot_rotor_matplotlib(
-                    nodes=nodes, check_sld=False, *args, **kwargs
+                nodes=nodes, check_sld=False, *args, **kwargs
             )
         elif plot_type == "bokeh":
-            return self._plot_rotor_bokeh(
-                    nodes=nodes, check_sld=False, *args, **kwargs
-            )
+            return self._plot_rotor_bokeh(nodes=nodes, check_sld=False, *args, **kwargs)
         else:
             raise ValueError(f"{plot_type} is not a valid plot type.")
 
@@ -1421,12 +1507,10 @@ class Rotor(object):
 
         if plot_type == "matplotlib":
             return self._plot_rotor_matplotlib(
-                    nodes=nodes, check_sld=True, *args, **kwargs
+                nodes=nodes, check_sld=True, *args, **kwargs
             )
         elif plot_type == "bokeh":
-            return self._plot_rotor_bokeh(
-                    nodes=nodes, check_sld=True, *args, **kwargs
-            )
+            return self._plot_rotor_bokeh(nodes=nodes, check_sld=True, *args, **kwargs)
         else:
             raise ValueError(f"{plot_type} is not a valid plot type.")
 
@@ -1921,7 +2005,11 @@ class Rotor(object):
                     -disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0]
                 )
                 BrgForceToReturn.append(
-                    np.around(-disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0], decimals=1)
+                    np.around(
+                        -disp_y[node]
+                        * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0],
+                        decimals=1,
+                    )
                 )
 
             # Disk Forces
@@ -1982,9 +2070,7 @@ class Rotor(object):
                 Bm = np.append(Bm, Bm[i] + Mx[i])
             self.Bm = Bm
 
-            sh_weight = np.around(
-                    sum(self.df_shaft["m"].values) * 9.8065, decimals=1
-            )
+            sh_weight = np.around(sum(self.df_shaft["m"].values) * 9.8065, decimals=1)
 
             force_data = {
                 "Static displacement vector": disp_y,
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 56a5a239..cb0c8c13 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -106,6 +106,7 @@ class ShaftElement(Element):
         rotary_inertia=True,
         gyroscopic=True,
         shear_method_calc="cowper",
+        tag=None,
     ):
 
         if type(material) is str:
@@ -126,6 +127,7 @@ class ShaftElement(Element):
             self.n_r = n + 1
 
         self.shear_method_calc = shear_method_calc
+        self.tag = tag
 
         self.L = float(L)
 
@@ -191,6 +193,9 @@ class ShaftElement(Element):
         else:
             return False
 
+    def __hash__(self):
+        return hash(self.tag)
+
     def save(self, file_name):
         data = self.load_data(file_name)
         data["ShaftElement"][str(self.n)] = {
@@ -796,6 +801,7 @@ class ShaftTaperedElement(Element):
         rotary_inertia=True,
         gyroscopic=True,
         shear_method_calc="cowper",
+        tag=None,
     ):
 
         if type(material) is str:
@@ -815,6 +821,8 @@ class ShaftTaperedElement(Element):
         if n is not None:
             self.n_r = n + 1
 
+        self.tag = tag
+
         self.shear_method_calc = shear_method_calc
 
         self.L = float(L)
@@ -924,6 +932,9 @@ class ShaftTaperedElement(Element):
         else:
             return False
 
+    def __hash__(self):
+        return hash(self.tag)
+
     def save(self, file_name):
         data = self.load_data(file_name)
         data["ShaftTaperedElement"][str(self.n)] = {
diff --git a/ross/utils.py b/ross/utils.py
index f617e2c9..553c31ec 100644
--- a/ross/utils.py
+++ b/ross/utils.py
@@ -1,12 +1,10 @@
 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.models import LogColorMapper
 from bokeh.palettes import Viridis256
-from bokeh.models import LinearColorMapper, LogColorMapper
+from bokeh.plotting import figure
+
+from ross.materials import Material
 
 
 class DataNotFoundError(Exception):
@@ -376,9 +374,8 @@ def visualize_matrix(rotor, matrix=None, frequency=None):
     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_text_font_size = "7pt"
     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))]
 
@@ -402,7 +399,7 @@ def visualize_matrix(rotor, matrix=None, frequency=None):
         0.95,
         0.95,
         source=data,
-        fill_color={"field": "value", "transform": mapper},
+        fill_color={"field": "pos_value", "transform": mapper},
         line_color=None,
     )
 


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index cf234de0..d1f9c46c 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -203,3 +203,28 @@ def test_from_table():
     bearing = BearingElement.from_table(0, bearing_file)
     assert bearing.n == 0
     assert math.isclose(bearing.w[2], 523.6, abs_tol=0.1)
+
+
+def test_bearing_link_global_index():
+    b0 = BearingElement(n=0, n_link=3, kxx=1, cxx=1)
+    idx = b0.dof_global_index()
+    assert idx.x_0 == 0
+    assert idx.y_0 == 1
+    print(idx)
+    assert idx.x_3 == 12
+    assert idx.y_3 == 13
+
+
+def test_bearing_link_matrices():
+    b0 = BearingElement(n=0, n_link=3, kxx=1, cxx=1)
+    # fmt: off
+    M = np.array(
+        [[1, 0, -1, 0],
+         [0, 1, 0, -1],
+         [-1, 0, 1, 0],
+         [0, -1, 0, 1]]
+    )
+    # fmt: on
+
+    assert_allclose(b0.K(0), M)
+    assert_allclose(b0.C(0), M)
diff --git a/ross/tests/test_point_mass.py b/ross/tests/test_point_mass.py
new file mode 100644
index 00000000..270c1050
--- /dev/null
+++ b/ross/tests/test_point_mass.py
@@ -0,0 +1,17 @@
+import numpy as np
+from ross.point_mass import PointMass
+from numpy.testing import assert_allclose
+
+
+def test_point_mass():
+    # fmt: off
+    m0 = np.array([[10.0, 0],
+                   [0, 10.0]])
+    # fmt: on
+
+    p = PointMass(n=0, m=10.0)
+
+    assert_allclose(m0, p.M())
+    assert_allclose(np.zeros((2, 2)), p.K())
+    assert_allclose(np.zeros((2, 2)), p.C())
+    assert_allclose(np.zeros((2, 2)), p.G())
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index fee2a570..99cd3856 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -10,6 +10,7 @@ from ross.materials import steel
 from ross.rotor_assembly import *
 from ross.rotor_assembly import MAC_modes
 from ross.shaft_element import *
+from ross.point_mass import PointMass
 
 test_dir = os.path.dirname(__file__)
 
@@ -717,7 +718,7 @@ def test_static_analysis_rotor3(rotor3):
                 -9.18114192e-04,
                 -8.08560219e-04,
                 -4.68788888e-04,
-                -5.56171636e-12
+                -5.56171636e-12,
             ]
         ),
         decimal=6,
@@ -789,7 +790,7 @@ def test_static_analysis_rotor5(rotor5):
                 -4.23416398e-04,
                 -6.31362481e-12,
                 4.28859620e-04,
-                8.52492302e-04
+                8.52492302e-04,
             ]
         ),
         decimal=6,
@@ -878,7 +879,7 @@ def test_static_analysis_rotor6(rotor6):
                 1.72933811e-04,
                 -1.02148266e-11,
                 -3.96409257e-04,
-                -9.20006704e-04
+                -9.20006704e-04,
             ]
         ),
         decimal=6,
@@ -1223,3 +1224,43 @@ def test_save_load():
     Rotor.remove("teste00000000000000001")
 
     assert a == b
+
+
+def test_rotor_link():
+    i_d = 0
+    o_d = 0.05
+    n = 6
+    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(
+        n=2, material=steel, width=0.07, i_d=0.05, o_d=0.28
+    )
+    disk1 = DiskElement.from_geometry(
+        n=4, material=steel, width=0.07, i_d=0.05, o_d=0.28
+    )
+
+    stfx = 1e6
+    stfy = 0.8e6
+    bearing0 = BearingElement(0, n_link=7, kxx=stfx, kyy=stfy, cxx=0)
+    support = BearingElement(7, kxx=stfx, kyy=stfy, cxx=0)
+    bearing1 = BearingElement(6, kxx=stfx, kyy=stfy, cxx=0)
+
+    point_mass = PointMass(7, m=1.0)
+
+    rotor = Rotor(
+        shaft_elem, [disk0, disk1], [bearing0, support, bearing1], [point_mass]
+    )
+    rotor.w = 1000 / 9.5492965964254
+
+    wn_expected = np.array(
+        [82.43809938, 87.5639844, 239.36021417, 261.15187937, 646.59858922, 688.428424]
+    )
+
+    assert_allclose(rotor.wn, wn_expected)