Repository: ross-rotordynamics/ross
Pull request #619 (2020-07-07 12:12:36)

Fix Isort deprecation issue
Related to Issue #618.


Patch:
diff --git a/.travis.yml b/.travis.yml
index 8ef5671a..19007d1f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -11,7 +11,7 @@ install:
 - pip install -e .[dev]
 before_script:
 - black --check ross
-- isort --recursive ross --check-only
+- isort ross --check-only
 script:
 - pytest ross --cov=ross
 after_success:
diff --git a/ross/__init__.py b/ross/__init__.py
index bed038f3..cb5f67bf 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -1,6 +1,5 @@
 __version__ = "0.3.3"
-
-import plotly.io as _pio
+from plotly import io as _pio
 
 import ross.plotly_theme
 
diff --git a/ross/api_report.py b/ross/api_report.py
index dc63f8f1..883ffc38 100644
--- a/ross/api_report.py
+++ b/ross/api_report.py
@@ -3,8 +3,8 @@ from copy import copy, deepcopy
 
 import numpy as np
 import pandas as pd
-import plotly.express as px
-import plotly.graph_objects as go
+from plotly import express as px
+from plotly import graph_objects as go
 from plotly.subplots import make_subplots
 from scipy.interpolate import interp1d
 from scipy.signal import argrelextrema
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index c53715aa..072f355c 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -10,9 +10,9 @@ import warnings
 from pathlib import Path
 
 import numpy as np
-import plotly.graph_objects as go
-import scipy.interpolate as interpolate
 import toml
+from plotly import graph_objects as go
+from scipy import interpolate as interpolate
 
 from ross.element import Element
 from ross.fluid_flow import fluid_flow as flow
diff --git a/ross/disk_element.py b/ross/disk_element.py
index 235c56ee..bb5eb025 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -8,8 +8,8 @@ import os
 from pathlib import Path
 
 import numpy as np
-import plotly.graph_objects as go
 import toml
+from plotly import graph_objects as go
 
 from ross.element import Element
 from ross.units import check_units
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index e3bb5265..8371e3cc 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -32,8 +32,6 @@ class FluidFlow:
         Number of points along the Z direction (direction of flow).
     ntheta: int
         Number of points along the direction theta. NOTE: ntheta must be odd.
-    nradius: int
-        Number of points along the direction r.
     length: float
         Length in the Z direction (m).
 
@@ -91,13 +89,11 @@ class FluidFlow:
     dtheta: float
         Range size in the theta direction.
     ntotal: int
-        Number of nodes in the grid., ntheta, nradius, n_interv_z, n_interv_theta,
+        Number of nodes in the grid., ntheta, n_interv_z, n_interv_theta,
     n_interv_z: int
         Number of intervals on Z.
     n_interv_theta: int
       Number of intervals on theta.
-    n_interv_radius: int
-        Number of intervals on r.
     p_mat_analytical : array of shape (nz, ntheta)
         The analytical pressure matrix.
     p_mat_numerical : array of shape (nz, ntheta)
@@ -114,8 +110,6 @@ class FluidFlow:
         The external radius in each position of the grid.
     ri : array of shape (nz, ntheta)
         The internal radius in each position of the grid.
-    z_list : array of shape (1, nz)
-        z along the object. It goes from 0 to lb.
     xre : array of shape (nz, ntheta)
         x value of the external radius.
     xri : array of shape (nz, ntheta)
@@ -124,6 +118,11 @@ class FluidFlow:
         y value of the external radius.
     yri : array of shape (nz, ntheta)
         y value of the internal radius.
+    z_list : array of shape (1, nz)
+        z along the object. It goes from 0 to lb.
+    gama : array of shape (nz, ntheta)
+        Points along the object in the tangential direction.
+        It ranges from 0 to 2 pi, starting at the largest spacing between rotor and stator.
     t : float
         Time.
     xp : float
@@ -156,7 +155,6 @@ class FluidFlow:
     >>> import numpy as np
     >>> nz = 8
     >>> ntheta = 64
-    >>> nradius = 8
     >>> length = 0.01
     >>> omega = 100.*2*np.pi/60
     >>> p_in = 0.
@@ -167,7 +165,7 @@ class FluidFlow:
     >>> density = 860.
     >>> eccentricity = 0.001
     >>> attitude_angle = np.pi
-    >>> my_fluid_flow = flow.FluidFlow(nz, ntheta, nradius, length,
+    >>> my_fluid_flow = flow.FluidFlow(nz, ntheta, length,
     ...                                omega, p_in, p_out, radius_rotor,
     ...                                radius_stator, viscosity, density,
     ...                                attitude_angle=attitude_angle, eccentricity=eccentricity,
@@ -187,7 +185,6 @@ class FluidFlow:
         self,
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -206,10 +203,8 @@ class FluidFlow:
             sys.exit("Either load or eccentricity must be given.")
         self.nz = nz
         self.ntheta = ntheta
-        self.nradius = nradius
         self.n_interv_z = nz - 1
         self.n_interv_theta = ntheta - 1
-        self.n_interv_radius = nradius - 1
         self.length = length
         self.ltheta = 2.0 * np.pi
         self.dz = length / self.n_interv_z
@@ -269,21 +264,18 @@ class FluidFlow:
         self.yi = self.eccentricity * np.sin(3 * np.pi / 2 + self.attitude_angle)
         self.re = np.zeros([self.nz, self.ntheta])
         self.ri = np.zeros([self.nz, self.ntheta])
-        self.z_list = np.zeros(self.nz)
         self.xre = np.zeros([self.nz, self.ntheta])
         self.xri = np.zeros([self.nz, self.ntheta])
         self.yre = np.zeros([self.nz, self.ntheta])
         self.yri = np.zeros([self.nz, self.ntheta])
+        self.z_list = np.zeros(self.nz)
+        self.gama = np.zeros([self.nz, self.ntheta])
+        self.t = 0
+        self.xp = 0
+        self.yp = 0
         self.p_mat_analytical = np.zeros([self.nz, self.ntheta])
-        self.c1 = np.zeros([self.nz, self.ntheta])
-        self.c2 = np.zeros([self.nz, self.ntheta])
-        self.c0w = np.zeros([self.nz, self.ntheta])
-        self.M = np.zeros([self.ntotal, self.ntotal])
-        self.f = np.zeros([self.ntotal, 1])
-        self.P = np.zeros([self.ntotal, 1])
         self.p_mat_numerical = np.zeros([self.nz, self.ntheta])
-        self.gama = np.zeros([self.nz, self.ntheta])
-        self.calculate_coefficients()
+        self.geometry_description()
         self.analytical_pressure_matrix_available = False
         self.numerical_pressure_matrix_available = False
         self.calculate_pressure_matrix_numerical()
@@ -319,19 +311,19 @@ class FluidFlow:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
                         # fmt: off
-                        self.p_mat_analytical[i][j] = (
+                        self.p_mat_analytical[i, j] = (
                                 ((-3 * self.viscosity * self.omega) / self.radial_clearance ** 2) *
                                 ((i * self.dz - (self.length / 2)) ** 2 - (self.length ** 2) / 4) *
                                 (self.eccentricity_ratio * np.sin(j * self.dtheta)) /
                                 (1 + self.eccentricity_ratio * np.cos(j * self.dtheta)) ** 3)
                         # fmt: on
-                        if self.p_mat_analytical[i][j] < 0:
-                            self.p_mat_analytical[i][j] = 0
+                        if self.p_mat_analytical[i, j] < 0:
+                            self.p_mat_analytical[i, j] = 0
             elif method == 1:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
                         # fmt: off
-                        self.p_mat_analytical[i][j] = (3 * self.viscosity / ((self.radial_clearance ** 2) *
+                        self.p_mat_analytical[i, j] = (3 * self.viscosity / ((self.radial_clearance ** 2) *
                                                                              (1. + self.eccentricity_ratio * np.cos(
                                                                                  j * self.dtheta)) ** 3)) * \
                                                       (-self.eccentricity_ratio * self.omega * np.sin(
@@ -339,18 +331,18 @@ class FluidFlow:
                                                       (((i * self.dz - (self.length / 2)) ** 2) - (
                                                               self.length ** 2) / 4)
                         # fmt: on
-                        if self.p_mat_analytical[i][j] < 0:
-                            self.p_mat_analytical[i][j] = 0
+                        if self.p_mat_analytical[i, j] < 0:
+                            self.p_mat_analytical[i, j] = 0
         elif self.bearing_type == "long_bearing" or force_type == "long":
             if method == 0:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
-                        self.p_mat_analytical[i][j] = (
+                        self.p_mat_analytical[i, j] = (
                             (
                                 6
                                 * self.viscosity
                                 * self.omega
-                                * (self.ri[i][j] / self.radial_clearance) ** 2
+                                * (self.radius_rotor / self.radial_clearance) ** 2
                                 * self.eccentricity_ratio
                                 * np.sin(self.dtheta * j)
                                 * (
@@ -368,8 +360,8 @@ class FluidFlow:
                             )
                             + self.p_in
                         )
-                        if self.p_mat_analytical[i][j] < 0:
-                            self.p_mat_analytical[i][j] = 0
+                        if self.p_mat_analytical[i, j] < 0:
+                            self.p_mat_analytical[i, j] = 0
         elif self.bearing_type == "medium_size":
             raise ValueError(
                 "The pressure matrix for a bearing that is neither short or long can only be calculated "
@@ -379,158 +371,209 @@ class FluidFlow:
         self.analytical_pressure_matrix_available = True
         return self.p_mat_analytical
 
-    def calculate_coefficients(self, direction=None):
-        """This function calculates the constants that form the Poisson equation
-        of the discrete pressure (central differences in the second
-        derivatives). It is executed when the class is instantiated.
+    def geometry_description(self):
+        """This function calculates the geometry description.
+        It is executed when the class is instantiated.
         Examples
         --------
         >>> my_fluid_flow = fluid_flow_example()
-        >>> my_fluid_flow.calculate_coefficients()
-        >>> my_fluid_flow.c0w # doctest: +ELLIPSIS
-        array([[...
+        >>> my_fluid_flow.geometry_description()
         """
         for i in range(0, self.nz):
             zno = i * self.dz
             self.z_list[i] = zno
-            eccentricity_error = False
             for j in range(0, self.ntheta):
                 # fmt: off
-                self.gama[i][j] = j * self.dtheta + np.pi / 2 + self.attitude_angle
-                [radius_external, self.xre[i][j], self.yre[i][j]] = \
-                    external_radius_function(self.gama[i][j], self.radius_stator)
-                [radius_internal, self.xri[i][j], self.yri[i][j]] = \
-                    internal_radius_function(self.gama[i][j], self.attitude_angle, self.radius_rotor,
+                self.gama[i, j] = j * self.dtheta + np.pi / 2 + self.attitude_angle
+                [radius_external, self.xre[i, j], self.yre[i, j]] = \
+                    external_radius_function(self.gama[i, j], self.radius_stator)
+                [radius_internal, self.xri[i, j], self.yri[i, j]] = \
+                    internal_radius_function(self.gama[i, j], self.attitude_angle, self.radius_rotor,
                                              self.eccentricity)
-                self.re[i][j] = radius_external
-                self.ri[i][j] = radius_internal
+                self.re[i, j] = radius_external
+                self.ri[i, j] = radius_internal
+                # fmt: on
+
+    def calculate_coefficients(self, direction=None):
+        """This function calculates the constants that form the Poisson equation
+        of the discrete pressure (central differences in the second
+        derivatives).
+        Parameters
+        ----------
+        direction: str
+            If defined, it calculates the model based on the disturbance in the chosen direction: 'x' or 'y'.
+        Returns
+        --------
+        c1, c2, c0w: matrix of float
+            Constants that form the Poisson equation.
+        >>> my_fluid_flow = fluid_flow_example()
+        >>> my_fluid_flow.calculate_coefficients()# doctest: +ELLIPSIS
+        (array([[...
+        """
+        c1 = np.zeros([self.nz, self.ntheta])
+        c2 = np.zeros([self.nz, self.ntheta])
+        c0w = np.zeros([self.nz, self.ntheta])
+
+        for i in range(0, self.nz):
+            eccentricity_error = False
+            for j in range(0, self.ntheta):
+                # fmt: off
 
                 w = self.omega * self.radius_rotor
 
-                k = (self.re[i][j] ** 2 * (np.log(self.re[i][j]) - 1 / 2) - self.ri[i][j] ** 2 *
-                     (np.log(self.ri[i][j]) - 1 / 2)) / (self.ri[i][j] ** 2 - self.re[i][j] ** 2)
-
-                self.c1[i][j] = (1 / (4 * self.viscosity)) * ((self.re[i][j] ** 2 * np.log(self.re[i][j]) -
-                                                               self.ri[i][j] ** 2 * np.log(self.ri[i][j]) +
-                                                               (self.re[i][j] ** 2 - self.ri[i][j] ** 2) *
-                                                               (k - 1)) - 2 * self.re[i][j] ** 2 * (
-                                                                      (np.log(self.re[i][j]) + k - 1 / 2) * np.log(
-                                                                       self.re[i][j] / self.ri[i][j])))
-
-                self.c2[i][j] = (- self.ri[i][j] ** 2) / (8 * self.viscosity) * \
-                                ((self.re[i][j] ** 2 - self.ri[i][j] ** 2 -
-                                  (self.re[i][j] ** 4 - self.ri[i][j] ** 4) /
-                                  (2 * self.ri[i][j] ** 2)) +
-                                 ((self.re[i][j] ** 2 - self.ri[i][j] ** 2) /
-                                  (self.ri[i][j] ** 2 *
-                                   np.log(self.re[i][j] / self.ri[i][j]))) *
-                                 (self.re[i][j] ** 2 * np.log(self.re[i][j] / self.ri[i][j]) -
-                                  (self.re[i][j] ** 2 - self.ri[i][j] ** 2) / 2))
-
-                self.c0w[i][j] = (- w * self.ri[i][j] *
-                                  (np.log(self.re[i][j] / self.ri[i][j]) *
-                                   (1 + (self.ri[i][j] ** 2) / (self.re[i][j] ** 2 - self.ri[i][j] ** 2)) - 1 / 2))
+                k = (self.re[i, j] ** 2 * (np.log(self.re[i, j]) - 1 / 2) - self.ri[i, j] ** 2 *
+                     (np.log(self.ri[i, j]) - 1 / 2)) / (self.ri[i, j] ** 2 - self.re[i, j] ** 2)
+
+                c1[i, j] = (1 / (4 * self.viscosity)) * ((self.re[i, j] ** 2 * np.log(self.re[i, j]) -
+                                                               self.ri[i, j] ** 2 * np.log(self.ri[i, j]) +
+                                                               (self.re[i, j] ** 2 - self.ri[i, j] ** 2) *
+                                                               (k - 1)) - 2 * self.re[i, j] ** 2 * (
+                                                                      (np.log(self.re[i, j]) + k - 1 / 2) * np.log(
+                                                                       self.re[i, j] / self.ri[i, j])))
+
+                c2[i, j] = (- self.ri[i, j] ** 2) / (8 * self.viscosity) * \
+                                ((self.re[i, j] ** 2 - self.ri[i, j] ** 2 -
+                                  (self.re[i, j] ** 4 - self.ri[i, j] ** 4) /
+                                  (2 * self.ri[i, j] ** 2)) +
+                                 ((self.re[i, j] ** 2 - self.ri[i, j] ** 2) /
+                                  (self.ri[i, j] ** 2 *
+                                   np.log(self.re[i, j] / self.ri[i, j]))) *
+                                 (self.re[i, j] ** 2 * np.log(self.re[i, j] / self.ri[i, j]) -
+                                  (self.re[i, j] ** 2 - self.ri[i, j] ** 2) / 2))
+
+                c0w[i, j] = (- w * self.ri[i, j] *
+                                  (np.log(self.re[i, j] / self.ri[i, j]) *
+                                   (1 + (self.ri[i, j] ** 2) / (self.re[i, j] ** 2 - self.ri[i, j] ** 2)) - 1 / 2))
                 if direction == "x":
-                    a = self.omegap * (self.xp) * np.cos(self.omegap * self.t)
-                    self.c0w[i][j] += self.ri[i][j] * a * np.sin(self.gama[i][j])
+                    a = self.omegap * self.xp * np.cos(self.omegap * self.t)
+                    c0w[i, j] += self.ri[i, j] * a * np.sin(self.gama[i, j])
                 elif direction == "y":
-                    b = self.omegap * (self.yp) * np.cos(self.omegap * self.t)
-                    self.c0w[i][j] -= self.ri[i][j] * b * np.cos(self.gama[i][j])
+                    b = self.omegap * self.yp * np.cos(self.omegap * self.t)
+                    c0w[i, j] -= self.ri[i, j] * b * np.cos(self.gama[i, j])
                 else:
-                    self.c0w[i][j] += 0
+                    c0w[i, j] += 0
                 # fmt: on
                 if not eccentricity_error:
-                    if abs(self.xri[i][j]) > abs(self.xre[i][j]) or abs(
-                        self.yri[i][j]
-                    ) > abs(self.yre[i][j]):
+                    if abs(self.xri[i, j]) > abs(self.xre[i, j]) or abs(
+                        self.yri[i, j]
+                    ) > abs(self.yre[i, j]):
                         eccentricity_error = True
             if eccentricity_error:
                 raise ValueError(
                     "Error: The given parameters create a rotor that is not inside the stator. "
                     "Check parameters and fix accordingly."
                 )
+        return c1, c2, c0w
 
-    def mounting_matrix(self):
+    def mounting_matrix(self, c1, c2, c0w):
         """This function assembles the matrix M and the independent vector f.
+        Parameters
+        ----------
+        c1, c2, c0w: matrix of float
+            Constants that form the Poisson equation.
+        Returns
+        --------
+        M: matrix of float
+            Matrix composed of coefficients that multiply the pressures at each point in the discrete domain.
+        f: array of float
+            Pressure independent terms.
         Examples
         --------
         >>> my_fluid_flow = fluid_flow_example()
-        >>> my_fluid_flow.mounting_matrix()
-        >>> my_fluid_flow.M # doctest: +ELLIPSIS
-        array([[...
+        >>> c1, c2, c0w = my_fluid_flow.calculate_coefficients()
+        >>> my_fluid_flow.mounting_matrix(c1, c2, c0w)# doctest: +ELLIPSIS
+        (array([[...
         """
         # fmt: off
+        M = np.zeros([self.ntotal, self.ntotal])
+        f = np.zeros([self.ntotal, 1])
         count = 0
         for x in range(self.ntheta):
-            self.M[count][count] = 1
-            self.f[count][0] = self.p_in
+            M[count, count] = 1
+            f[count, 0] = self.p_in
             count = count + self.nz - 1
-            self.M[count][count] = 1
-            self.f[count][0] = self.p_out
+            M[count, count] = 1
+            f[count, 0] = self.p_out
             count = count + 1
         count = 0
         for x in range(self.nz - 2):
-            self.M[self.ntotal - self.nz + 1 + count][1 + count] = 1
-            self.M[self.ntotal - self.nz + 1 + count][self.ntotal - self.nz + 1 + count] = -1
+            M[self.ntotal - self.nz + 1 + count, 1 + count] = 1
+            M[self.ntotal - self.nz + 1 + count, self.ntotal - self.nz + 1 + count] = -1
             count = count + 1
         count = 1
         j = 0
         for i in range(1, self.nz - 1):
-            a = (1 / self.dtheta ** 2) * (self.c1[i][self.ntheta - 1])
-            self.M[count][self.ntotal - 2 * self.nz + count] = a
-            b = (1 / self.dz ** 2) * (self.c2[i - 1, j])
-            self.M[count][count - 1] = b
-            c = -((1 / self.dtheta ** 2) * ((self.c1[i][j]) + self.c1[i][self.ntheta - 1])
-                  + (1 / self.dz ** 2) * (self.c2[i][j] + self.c2[i - 1][j]))
-            self.M[count, count] = c
-            d = (1 / self.dz ** 2) * (self.c2[i][j])
-            self.M[count][count + 1] = d
-            e = (1 / self.dtheta ** 2) * (self.c1[i][j])
-            self.M[count][count + self.nz] = e
+            a = (1 / self.dtheta ** 2) * (c1[i, self.ntheta - 1])
+            M[count, self.ntotal - 2 * self.nz + count] = a
+            b = (1 / self.dz ** 2) * (c2[i - 1, j])
+            M[count, count - 1] = b
+            c = -((1 / self.dtheta ** 2) * ((c1[i, j]) + c1[i, self.ntheta - 1])
+                  + (1 / self.dz ** 2) * (c2[i, j] + c2[i - 1, j]))
+            M[count, count] = c
+            d = (1 / self.dz ** 2) * (c2[i, j])
+            M[count, count + 1] = d
+            e = (1 / self.dtheta ** 2) * (c1[i, j])
+            M[count, count + self.nz] = e
             count = count + 1
         count = self.nz + 1
         for j in range(1, self.ntheta - 1):
             for i in range(1, self.nz - 1):
-                a = (1 / self.dtheta ** 2) * (self.c1[i, j - 1])
-                self.M[count][count - self.nz] = a
-                b = (1 / self.dz ** 2) * (self.c2[i - 1][j])
-                self.M[count][count - 1] = b
-                c = -((1 / self.dtheta ** 2) * ((self.c1[i][j]) + self.c1[i][j - 1])
-                      + (1 / self.dz ** 2) * (self.c2[i][j] + self.c2[i - 1][j]))
-                self.M[count, count] = c
-                d = (1 / self.dz ** 2) * (self.c2[i][j])
-                self.M[count][count + 1] = d
-                e = (1 / self.dtheta ** 2) * (self.c1[i][j])
-                self.M[count][count + self.nz] = e
+                a = (1 / self.dtheta ** 2) * (c1[i, j - 1])
+                M[count, count - self.nz] = a
+                b = (1 / self.dz ** 2) * (c2[i - 1, j])
+                M[count, count - 1] = b
+                c = -((1 / self.dtheta ** 2) * ((c1[i, j]) + c1[i, j - 1])
+                      + (1 / self.dz ** 2) * (c2[i, j] + c2[i - 1, j]))
+                M[count, count] = c
+                d = (1 / self.dz ** 2) * (c2[i, j])
+                M[count, count + 1] = d
+                e = (1 / self.dtheta ** 2) * (c1[i, j])
+                M[count, count + self.nz] = e
                 count = count + 1
             count = count + 2
         count = 1
         for j in range(self.ntheta - 1):
             for i in range(1, self.nz - 1):
                 if j == 0:
-                    self.f[count][0] = (self.c0w[i][j] - self.c0w[i][self.ntheta - 1]) / self.dtheta
+                    f[count, 0] = (c0w[i, j] - c0w[i, self.ntheta - 1]) / self.dtheta
                 else:
-                    self.f[count][0] = (self.c0w[i, j] - self.c0w[i, j - 1]) / self.dtheta
+                    f[count, 0] = (c0w[i, j] - c0w[i, j - 1]) / self.dtheta
                 count = count + 1
             count = count + 2
         # fmt: on
+        return M, f
 
-    def resolves_matrix(self):
+    def resolves_matrix(self, M, f):
         """This function resolves the linear system [M]{P} = {f}.
+         Parameters
+        ----------
+        M: matrix of float
+            Matrix composed of coefficients that multiply the pressures at each point in the discrete domain.
+        f: array of float
+            Pressure independent terms.
+        Returns
+        --------
+        P: array of floats
+            Pressure. Unknowns of the finite difference system.
         Examples
         --------
         >>> my_fluid_flow = fluid_flow_example()
-        >>> my_fluid_flow.mounting_matrix()
-        >>> my_fluid_flow.resolves_matrix()
-        >>> my_fluid_flow.P # doctest: +ELLIPSIS
+        >>> c1, c2, c0w = my_fluid_flow.calculate_coefficients()
+        >>> M, f = my_fluid_flow.mounting_matrix(c1, c2, c0w)
+        >>> my_fluid_flow.resolves_matrix(M, f)# doctest: +ELLIPSIS
         array([[...
         """
-        sparse_matrix = sp.sparse.csc_matrix(self.M)
-        self.P = sp.sparse.linalg.spsolve(sparse_matrix, self.f)
-        self.P.shape = (self.P.size, 1)
+        sparse_matrix = sp.sparse.csc_matrix(M)
+        P = sp.sparse.linalg.spsolve(sparse_matrix, f)
+        P.shape = (P.size, 1)
+        return P
 
-    def calculate_pressure_matrix_numerical(self):
+    def calculate_pressure_matrix_numerical(self, direction=None):
         """This function calculates the pressure matrix numerically.
+        Parameters
+        ----------
+        direction: str
+            If defined, it calculates the model based on the disturbance in the chosen direction: 'x' or 'y'.
         Returns
         -------
         p_mat_numerical: matrix of float
@@ -541,15 +584,21 @@ class FluidFlow:
         >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
         array([[...
         """
-        self.mounting_matrix()
-        self.resolves_matrix()
+        if direction == "x":
+            c1, c2, c0w = self.calculate_coefficients(direction="x")
+        elif direction == "y":
+            c1, c2, c0w = self.calculate_coefficients(direction="y")
+        else:
+            c1, c2, c0w = self.calculate_coefficients()
+        M, f = self.mounting_matrix(c1, c2, c0w)
+        P = self.resolves_matrix(M, f)
         for i in range(self.nz):
             for j in range(self.ntheta):
                 k = j * self.nz + i
-                if self.P[k] < 0:
-                    self.p_mat_numerical[i][j] = 0
+                if P[k] < 0:
+                    self.p_mat_numerical[i, j] = 0
                 else:
-                    self.p_mat_numerical[i][j] = self.P[k]
+                    self.p_mat_numerical[i, j] = P[k]
         self.numerical_pressure_matrix_available = True
         return self.p_mat_numerical
 
@@ -575,7 +624,6 @@ def fluid_flow_example():
     my_pressure_matrix = FluidFlow(
         nz=8,
         ntheta=32,
-        nradius=8,
         length=0.04,
         omega=100.0 * 2 * np.pi / 60,
         p_in=0.0,
@@ -611,7 +659,6 @@ def fluid_flow_example2():
     """
     nz = 8
     ntheta = 16
-    nradius = 8
     length = 0.03
     omega = 157.1
     p_in = 0.0
@@ -624,7 +671,6 @@ def fluid_flow_example2():
     return FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index e3c20087..9813904c 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -197,8 +197,8 @@ def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
             * fluid_flow_object.xp
             * np.cos(fluid_flow_object.omegap * fluid_flow_object.t)
         )
-        fluid_flow_object.calculate_coefficients(direction="x")
-        fluid_flow_object.calculate_pressure_matrix_numerical()
+        fluid_flow_object.geometry_description()
+        fluid_flow_object.calculate_pressure_matrix_numerical(direction="x")
         [
             radial_force[i],
             tangential_force[i],
@@ -217,8 +217,8 @@ def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
             * fluid_flow_object.yp
             * np.cos(fluid_flow_object.omegap * fluid_flow_object.t)
         )
-        fluid_flow_object.calculate_coefficients(direction="y")
-        fluid_flow_object.calculate_pressure_matrix_numerical()
+        fluid_flow_object.geometry_description()
+        fluid_flow_object.calculate_pressure_matrix_numerical(direction="y")
         [
             radial_force[i],
             tangential_force[i],
@@ -226,7 +226,7 @@ def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
             force_yy[i],
         ] = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
         move_rotor_center(fluid_flow_object, 0, -delta_y)
-        fluid_flow_object.calculate_coefficients()
+        fluid_flow_object.geometry_description()
         fluid_flow_object.calculate_pressure_matrix_numerical()
 
         X1[i] = [1, dx[i], xdot[i]]
diff --git a/ross/fluid_flow/fluid_flow_graphics.py b/ross/fluid_flow/fluid_flow_graphics.py
index 69d69be7..c19849f7 100644
--- a/ross/fluid_flow/fluid_flow_graphics.py
+++ b/ross/fluid_flow/fluid_flow_graphics.py
@@ -1,5 +1,5 @@
 import numpy as np
-import plotly.graph_objects as go
+from plotly import graph_objects as go
 
 from ross.plotly_theme import tableau_colors
 
@@ -394,11 +394,7 @@ def plot_pressure_theta_cylindrical(
         else:
             p_mat = fluid_flow_object.p_mat_numerical
 
-    r = np.linspace(
-        fluid_flow_object.radius_rotor,
-        fluid_flow_object.radius_stator,
-        fluid_flow_object.nradius,
-    )
+    r = np.linspace(fluid_flow_object.radius_rotor, fluid_flow_object.radius_stator)
     theta = np.linspace(
         0.0, 2.0 * np.pi + fluid_flow_object.dtheta / 2, fluid_flow_object.ntheta
     )
diff --git a/ross/plotly_theme.py b/ross/plotly_theme.py
index b038bf47..d0d3ba6e 100644
--- a/ross/plotly_theme.py
+++ b/ross/plotly_theme.py
@@ -1,5 +1,5 @@
-import plotly.graph_objects as go
-import plotly.io as pio
+from plotly import graph_objects as go
+from plotly import io as pio
 
 # tableau colors
 tableau_colors = {
diff --git a/ross/point_mass.py b/ross/point_mass.py
index e2360c75..8452e3e8 100644
--- a/ross/point_mass.py
+++ b/ross/point_mass.py
@@ -6,8 +6,8 @@ import os
 from pathlib import Path
 
 import numpy as np
-import plotly.graph_objects as go
 import toml
+from plotly import graph_objects as go
 
 from ross.element import Element
 from ross.units import check_units
diff --git a/ross/results.py b/ross/results.py
index 1594dc20..a6d24a31 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -5,10 +5,10 @@ This module returns graphs for each type of analyses in rotor_assembly.py.
 import copy
 
 import numpy as np
-import plotly.graph_objects as go
-import scipy.linalg as la
+from plotly import graph_objects as go
 from plotly.subplots import make_subplots
 from scipy import interpolate
+from scipy import linalg as la
 
 from ross.plotly_theme import tableau_colors
 
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 74adcbcf..a1caedb3 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -8,18 +8,18 @@ from copy import copy, deepcopy
 from itertools import chain, cycle
 from pathlib import Path
 
-import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-import plotly.express as px
-import plotly.graph_objects as go
-import scipy.io as sio
-import scipy.linalg as la
-import scipy.signal as signal
-import scipy.sparse.linalg as las
 import toml
+from matplotlib import pyplot as plt
+from plotly import express as px
+from plotly import graph_objects as go
+from scipy import io as sio
+from scipy import linalg as la
+from scipy import signal as signal
 from scipy.interpolate import UnivariateSpline
 from scipy.optimize import newton
+from scipy.sparse import linalg as las
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
                                        BearingElement6DoF,
@@ -1043,7 +1043,7 @@ class Rotor(object):
         >>> rotor = rotor_example()
         >>> speed_range = rotor._clustering_points(num_modes=12, num_points=5)
         >>> speed_range.shape
-        (60,)
+        (61,)
         """
         critical_speeds = self.run_critical_speed(num_modes=num_modes, rtol=rtol)
         omega = critical_speeds.wd
@@ -1067,6 +1067,7 @@ class Rotor(object):
 
         omega = omega.reshape((len(omega), 1))
         speed_range = np.sort(np.ravel(np.concatenate((omega / a, omega * a))))
+        speed_range = np.insert(speed_range, 0, 0)
 
         return speed_range
 
@@ -1353,7 +1354,7 @@ class Rotor(object):
 
         >>> response = rotor.run_freq_response(cluster_points=True, num_points=5)
         >>> response.speed_range.shape
-        (60,)
+        (61,)
 
         # plot frequency response function:
         >>> fig = response.plot(inp=13, out=13)
@@ -1468,7 +1469,7 @@ class Rotor(object):
         ...     force=force, cluster_points=True, num_modes=12, num_points=5
         ... )
         >>> response.speed_range.shape
-        (60,)
+        (61,)
         """
         if speed_range is None:
             if cluster_points:
@@ -1633,7 +1634,7 @@ class Rotor(object):
         ...     node=3, magnitude=0.01, phase=0.0, cluster_points=True, num_points=5
         ... )
         >>> response2.speed_range.shape
-        (60,)
+        (61,)
 
         plot unbalance response:
         >>> fig = response.plot(dof=13)
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 8ba9f7c4..eb44c9ac 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -8,8 +8,8 @@ import os
 from pathlib import Path
 
 import numpy as np
-import plotly.graph_objects as go
 import toml
+from plotly import graph_objects as go
 
 from ross.element import Element
 from ross.materials import Material, steel
diff --git a/ross/stochastic/st_results.py b/ross/stochastic/st_results.py
index ca1e4e21..c18ae745 100644
--- a/ross/stochastic/st_results.py
+++ b/ross/stochastic/st_results.py
@@ -3,9 +3,9 @@
 This module returns graphs for each type of analyses in st_rotor_assembly.py.
 """
 import numpy as np
-import plotly.express as px
-import plotly.graph_objects as go
-import plotly.io as pio
+from plotly import express as px
+from plotly import graph_objects as go
+from plotly import io as pio
 from plotly.subplots import make_subplots
 
 pio.renderers.default = "browser"
diff --git a/ross/stochastic/st_results_elements.py b/ross/stochastic/st_results_elements.py
index 46d38c82..31eee27c 100644
--- a/ross/stochastic/st_results_elements.py
+++ b/ross/stochastic/st_results_elements.py
@@ -3,8 +3,8 @@
 This modules provides functions to plot the elements statistic data.
 """
 import numpy as np
-import plotly.graph_objects as go
-import plotly.io as pio
+from plotly import graph_objects as go
+from plotly import io as pio
 from plotly.subplots import make_subplots
 from scipy.stats import gaussian_kde
 
diff --git a/ross/utils.py b/ross/utils.py
index bdb986a3..1f77a2a2 100644
--- a/ross/utils.py
+++ b/ross/utils.py
@@ -2,7 +2,7 @@ import re
 
 import numpy as np
 import pandas as pd
-import plotly.graph_objects as go
+from plotly import graph_objects as go
 
 
 class DataNotFoundError(Exception):


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index 2515cd4d..2951e5fb 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -9,6 +9,7 @@ from numpy.testing import assert_allclose
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
                                        BearingElement6DoF,
+                                       MagneticBearingElement,
                                        RollerBearingElement)
 
 # fmt: on
@@ -292,6 +293,41 @@ def test_roller_bearing_element():
     assert_allclose(rollerbearing.G(), G)
 
 
+def test_magnetic_bearing_element():
+    n = 0
+    g0 = 1e-3
+    i0 = 1.0
+    ag = 1e-4
+    nw = 200
+    alpha = 0.392
+    kp_pid = 1.0
+    kd_pid = 1.0
+    k_amp = 1.0
+    k_sense = 1.0
+    tag = "magneticbearing"
+    mbearing = MagneticBearingElement(
+        n=n,
+        g0=g0,
+        i0=i0,
+        ag=ag,
+        nw=nw,
+        alpha=alpha,
+        kp_pid=kp_pid,
+        kd_pid=kd_pid,
+        k_amp=k_amp,
+        k_sense=k_sense,
+    )
+    M = np.array([[0.0, 0.0], [0.0, 0.0]])
+    K = np.array([[-4640.62337718, 0.0], [0.0, -4640.62337718]])
+    C = np.array([[4.64526865, 0.0], [0.0, 4.64526865]])
+    G = np.array([[0.0, 0.0], [0.0, 0.0]])
+
+    assert_allclose(mbearing.M(), M)
+    assert_allclose(mbearing.K(0), K)
+    assert_allclose(mbearing.C(0), C)
+    assert_allclose(mbearing.G(), G)
+
+
 @pytest.fixture
 def bearing_6dof():
     bearing_6dof = BearingElement6DoF(
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 5ea73416..68808ab8 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -2,9 +2,9 @@
 import math
 
 import numpy as np
-import plotly.graph_objects as go
 import pytest
 from numpy.testing import assert_allclose
+from plotly import graph_objects as go
 
 from ross.fluid_flow import fluid_flow as flow
 from ross.fluid_flow.fluid_flow_coefficients import (
@@ -23,7 +23,6 @@ from ross.fluid_flow.fluid_flow_graphics import (
 def fluid_flow_short_eccentricity():
     nz = 8
     ntheta = 32
-    nradius = 8
     omega = 100.0 * 2 * np.pi / 60
     p_in = 0.0
     p_out = 0.0
@@ -36,7 +35,6 @@ def fluid_flow_short_eccentricity():
     return flow.FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -53,7 +51,6 @@ def fluid_flow_short_eccentricity():
 def fluid_flow_short_friswell(set_load=True):
     nz = 8
     ntheta = 32
-    nradius = 8
     length = 0.03
     omega = 157.1
     p_in = 0.0
@@ -68,7 +65,6 @@ def fluid_flow_short_friswell(set_load=True):
         return flow.FluidFlow(
             nz,
             ntheta,
-            nradius,
             length,
             omega,
             p_in,
@@ -84,7 +80,6 @@ def fluid_flow_short_friswell(set_load=True):
         return flow.FluidFlow(
             nz,
             ntheta,
-            nradius,
             length,
             omega,
             p_in,
@@ -185,7 +180,6 @@ def test_damping_matrix():
 def fluid_flow_short_numerical():
     nz = 8
     ntheta = 32
-    nradius = 8
     length = 0.01
     omega = 100.0 * 2 * np.pi / 60
     p_in = 0.0
@@ -199,7 +193,6 @@ def fluid_flow_short_numerical():
     return flow.FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -217,7 +210,6 @@ def fluid_flow_short_numerical():
 def fluid_flow_long_numerical():
     nz = 8
     ntheta = 32
-    nradius = 8
     omega = 100.0 * 2 * np.pi / 60
     p_in = 0.0
     p_out = 0.0
@@ -231,7 +223,6 @@ def fluid_flow_long_numerical():
     return flow.FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
diff --git a/ross/tests/test_stochastic_elements.py b/ross/tests/test_stochastic_elements.py
index 1c89dfcc..373a52e8 100644
--- a/ross/tests/test_stochastic_elements.py
+++ b/ross/tests/test_stochastic_elements.py
@@ -6,8 +6,8 @@ Tests for:
     st_point_mass.py
 """
 import numpy as np
-import plotly.graph_objects as go
 import pytest
+from plotly import graph_objects as go
 
 from ross.stochastic.st_bearing_seal_element import ST_BearingElement
 from ross.stochastic.st_disk_element import ST_DiskElement