Repository: ross-rotordynamics/ross
Pull request #628 (2020-07-13 19:50:09)

Add intersection function
This PR adds the intersection function to the utils module.
The Campbell plot and ucs functions will use this to find the intersection point between curves.


Patch:
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index f5d73067..a5545096 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -7,7 +7,6 @@ and 2 element options with 8 or 12 degrees of freedom.
 # fmt: off
 import os
 import warnings
-from pathlib import Path
 
 import numpy as np
 import toml
diff --git a/ross/results.py b/ross/results.py
index 79ceafc1..0f176b9e 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -7,11 +7,11 @@ import copy
 import numpy as np
 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
 from ross.units import Q_
+from ross.utils import intersection
 
 
 class CriticalSpeedResults:
@@ -764,41 +764,27 @@ class CampbellResults:
                 log_dec_i = log_dec[:, i]
 
                 for harm in harmonics:
-                    idx = np.argwhere(np.diff(np.sign(w_i - harm * speed_range)))
-                    idx = idx.flatten()
-                    if len(idx) != 0:
-                        idx = idx[0]
-
-                        interpolated = interpolate.interp1d(
-                            x=[speed_range[idx], speed_range[idx + 1]],
-                            y=[w_i[idx], w_i[idx + 1]],
-                            kind="linear",
-                        )
-                        xnew = np.linspace(
-                            speed_range[idx],
-                            speed_range[idx + 1],
-                            num=30,
-                            endpoint=True,
-                        )
-                        ynew = interpolated(xnew)
-                        idx = np.argwhere(
-                            np.diff(np.sign(ynew - harm * xnew))
-                        ).flatten()
-
-                        fig.add_trace(
-                            go.Scatter(
-                                x=xnew[idx],
-                                y=ynew[idx],
-                                mode="markers",
-                                marker=dict(symbol="x", color="black"),
-                                name="Crit. Speed",
-                                legendgroup="Crit. Speed",
-                                showlegend=False,
-                                hovertemplate=(
-                                    f"Frequency ({frequency_units}): %{{x:.2f}}<br> Critical Speed ({frequency_units}): %{{y:.2f}}"
-                                ),
-                            )
+                    x1 = speed_range
+                    y1 = w_i
+                    x2 = speed_range
+                    y2 = harm * speed_range
+
+                    x, y = intersection(x1, y1, x2, y2)
+
+                    fig.add_trace(
+                        go.Scatter(
+                            x=x,
+                            y=y,
+                            mode="markers",
+                            marker=dict(symbol="x", color="black"),
+                            name="Crit. Speed",
+                            legendgroup="Crit. Speed",
+                            showlegend=False,
+                            hovertemplate=(
+                                f"Frequency ({frequency_units}): %{{x:.2f}}<br>Critical Speed ({frequency_units}): %{{y:.2f}}"
+                            ),
                         )
+                    )
 
                 whirl_mask = whirl_i == whirl_dir
                 if whirl_mask.shape[0] == 0:
@@ -827,7 +813,7 @@ class CampbellResults:
                     x=speed_range,
                     y=h * speed_range,
                     mode="lines",
-                    line=dict(dash="dashdot"),
+                    line=dict(dash="dashdot", color=list(tableau_colors)[j]),
                     name="{}x speed".format(h),
                     hoverinfo="none",
                 )
@@ -1892,7 +1878,7 @@ class ForcedResponseResults:
                 yaxis=dict(title=dict(text=f"Amplitude - X ({displacement_units})")),
                 zaxis=dict(title=dict(text=f"Amplitude - Y ({displacement_units})")),
             ),
-            title=dict(text=f"Deflected Shape<br>" f"Speed = {speed}"),
+            title=dict(text=f"Deflected Shape<br>Speed = {speed}"),
             **kwargs,
         )
 
@@ -2208,7 +2194,7 @@ class StaticResults:
                     name=f"Shaft {count}",
                     showlegend=True,
                     hovertemplate=(
-                        f"Rotor Length ({rotor_length_units}): %{{x:.2f}}<br> Displacement ({deformation_units}): %{{y:.2e}}"
+                        f"Rotor Length ({rotor_length_units}): %{{x:.2f}}<br>Displacement ({deformation_units}): %{{y:.2e}}"
                     ),
                 )
             )
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 1f6251b9..01167326 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -33,7 +33,7 @@ from ross.results import (CampbellResults, ConvergenceResults,
                           StaticResults, SummaryResults, TimeResponseResults)
 from ross.shaft_element import ShaftElement, ShaftElement6DoF
 from ross.units import Q_, check_units
-from ross.utils import convert
+from ross.utils import intersection
 
 # fmt: on
 
@@ -1963,21 +1963,17 @@ class Rotor(object):
         if not np.isnan(bearing0.frequency).all():
             for j in range(rotor_wn.shape[0]):
                 for coeff in ["kxx", "kyy"]:
-                    with warnings.catch_warnings():
-                        warnings.simplefilter("ignore")
-                        frequency_curve = UnivariateSpline(rotor_wn[j], stiffness_log)
-                    coeff_curve = getattr(bearing0, coeff).interpolated
-                    # search diff in several points to have a good precision
-                    search_range = np.linspace(
-                        bearing0.frequency[0], bearing0.frequency[-1], 10000
-                    )
-                    intersection_idx = np.argmin(
-                        abs(frequency_curve(search_range) - coeff_curve(search_range))
-                    )
-                    intersection_points["y"].append(search_range[intersection_idx])
-                    intersection_points["x"].append(
-                        frequency_curve(search_range[intersection_idx])
-                    )
+                    x1 = rotor_wn[j]
+                    y1 = stiffness_log
+                    x2 = bearing0.frequency
+                    y2 = getattr(bearing0, coeff).coefficient
+                    x, y = intersection(x1, y1, x2, y2)
+                    try:
+                        intersection_points["y"].append(float(x))
+                        intersection_points["x"].append(float(y))
+                    except TypeError:
+                        # pass if x/y is empty
+                        pass
 
         return stiffness_log, rotor_wn, bearing0, intersection_points
 
diff --git a/ross/utils.py b/ross/utils.py
index 1f77a2a2..cd930f43 100644
--- a/ross/utils.py
+++ b/ross/utils.py
@@ -428,3 +428,112 @@ def convert(name):
 
     s1 = re.sub("(.)([A-Z][a-z]+)", r"\1_\2", name)
     return re.sub("([a-z0-9])([A-Z])", r"\1_\2", s1).lower()
+
+
+def intersection(x1, y1, x2, y2):
+    """
+    Intersection code from https://github.com/sukhbinder/intersection
+    MIT License
+
+    Copyright (c) 2017 Sukhbinder Singh
+
+    Permission is hereby granted, free of charge, to any person obtaining a copy
+    of this software and associated documentation files (the "Software"), to deal
+    in the Software without restriction, including without limitation the rights
+    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+    copies of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be included in all
+    copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+    SOFTWARE.
+
+    INTERSECTIONS Intersections of curves.
+       Computes the (x,y) locations where two curves intersect.  The curves
+       can be broken with NaNs or have vertical segments.
+
+    Usage:
+    x, y = intersection(x1, y1, x2, y2)
+
+    Example:
+    a, b = 1, 2
+    phi = np.linspace(3, 10, 100)
+    x1 = a*phi - b*np.sin(phi)
+    y1 = a - b*np.cos(phi)
+
+    x2=phi
+    y2=np.sin(phi)+2
+    x,y=intersection(x1,y1,x2,y2)
+
+    plt.plot(x1,y1,c='r')
+    plt.plot(x2,y2,c='g')
+    plt.plot(x,y,'*k')
+    plt.show()
+
+    """
+
+    def _rect_inter_inner(x1, x2):
+        n1 = x1.shape[0] - 1
+        n2 = x2.shape[0] - 1
+        X1 = np.c_[x1[:-1], x1[1:]]
+        X2 = np.c_[x2[:-1], x2[1:]]
+        S1 = np.tile(X1.min(axis=1), (n2, 1)).T
+        S2 = np.tile(X2.max(axis=1), (n1, 1))
+        S3 = np.tile(X1.max(axis=1), (n2, 1)).T
+        S4 = np.tile(X2.min(axis=1), (n1, 1))
+        return S1, S2, S3, S4
+
+    def _rectangle_intersection_(x1, y1, x2, y2):
+        S1, S2, S3, S4 = _rect_inter_inner(x1, x2)
+        S5, S6, S7, S8 = _rect_inter_inner(y1, y2)
+
+        C1 = np.less_equal(S1, S2)
+        C2 = np.greater_equal(S3, S4)
+        C3 = np.less_equal(S5, S6)
+        C4 = np.greater_equal(S7, S8)
+
+        ii, jj = np.nonzero(C1 & C2 & C3 & C4)
+        return ii, jj
+
+    x1 = np.asarray(x1)
+    x2 = np.asarray(x2)
+    y1 = np.asarray(y1)
+    y2 = np.asarray(y2)
+
+    ii, jj = _rectangle_intersection_(x1, y1, x2, y2)
+    n = len(ii)
+
+    dxy1 = np.diff(np.c_[x1, y1], axis=0)
+    dxy2 = np.diff(np.c_[x2, y2], axis=0)
+
+    T = np.zeros((4, n))
+    AA = np.zeros((4, 4, n))
+    AA[0:2, 2, :] = -1
+    AA[2:4, 3, :] = -1
+    AA[0::2, 0, :] = dxy1[ii, :].T
+    AA[1::2, 1, :] = dxy2[jj, :].T
+
+    BB = np.zeros((4, n))
+    BB[0, :] = -x1[ii].ravel()
+    BB[1, :] = -x2[jj].ravel()
+    BB[2, :] = -y1[ii].ravel()
+    BB[3, :] = -y2[jj].ravel()
+
+    for i in range(n):
+        try:
+            T[:, i] = np.linalg.solve(AA[:, :, i], BB[:, i])
+        except:
+            T[:, i] = np.Inf
+
+    in_range = (T[0, :] >= 0) & (T[1, :] >= 0) & (T[0, :] <= 1) & (T[1, :] <= 1)
+
+    xy0 = T[2:, in_range]
+    xy0 = xy0.T
+    return xy0[:, 0], xy0[:, 1]


Test patch:
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index 8cda96ca..1c516811 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -1619,9 +1619,9 @@ def test_ucs_calc(rotor8):
     exp_stiffness_range = np.array([1000000.0, 1832980.710832, 3359818.286284])
     exp_rotor_wn = np.array([86.658114, 95.660573, 101.868254])
     exp_intersection_points_x = np.array(
-        [9632232.337582, 9632232.337582, 10333174.417919]
+        [10058123.652648, 10058123.652648, 10363082.398797]
     )
-    exp_intersection_points_y = np.array([107.920792, 107.920792, 411.881188])
+    exp_intersection_points_y = np.array([107.542416, 107.542416, 409.451575])
     stiffness_range, rotor_wn, bearing, intersection_points = rotor8._calc_ucs()
     assert_allclose(stiffness_range[:3], exp_stiffness_range)
     assert_allclose(rotor_wn[0, :3], exp_rotor_wn)