updated WEC.from_impedance to reflect correct dimension order

tests/test_integration.py

@@ -120,10 +120,9 @@ def wec_from_floatingbody(f1, nfreq, fb, pto):
 @pytest.fixture(scope='module')
 def wec_from_impedance(bem, pto, fb):
     """Simple WEC: 1 DOF, no constraints."""
-    bemc = bem.copy().transpose(
-        "radiating_dof", "influenced_dof", "omega", "wave_direction")
+    bemc = bem.copy()
     omega = bemc['omega'].values
-    w = np.expand_dims(omega, [0, 1])
+    w = np.expand_dims(omega, [1,2])
     A = bemc['added_mass'].values
     B = bemc['radiation_damping'].values
     fb.center_of_mass = [0, 0, 0]

wecopttool/core.py

@@ -503,7 +503,7 @@ def from_impedance(
             Frequency vector [:math:`Hz`] not including the zero frequency,
             :python:`freqs = [f1, 2*f1, ..., nfreq*f1]`.
         impedance
-            Complex impedance of size :python:`(ndof, ndof, nfreq)`.
+            Complex impedance of size :python:`(nfreq, ndof, ndof)`.
         exc_coeff
             Complex excitation transfer function of size
             :python:`(ndof, nfreq)`.
@@ -537,17 +537,18 @@ def from_impedance(
         # impedance matrix shape
         shape = impedance.shape
         ndim = impedance.ndim
-        if (ndim!=3) or (shape[0]!=shape[1]) or (shape[2]!=nfreq):
+        if (ndim!=3) or (shape[1]!=shape[2]) or (shape[0]!=nfreq):
             raise ValueError(
-                "'impedance' must have shape '(ndof, ndof, nfreq)'.")
+                "'impedance' must have shape '(nfreq, ndof, ndof)'.")
 
         impedance = check_impedance(impedance, min_damping)
 
         # impedance force
         omega = freqs * 2*np.pi
-        transfer_func = impedance * (1j*omega)
+        omega0 = np.expand_dims(omega, [1,2])
+        transfer_func = impedance * (1j*omega0)
         transfer_func0 = np.expand_dims(hydrostatic_stiffness, 2)
-        transfer_func = np.concatenate([transfer_func0, transfer_func], 2)
+        transfer_func = np.concatenate([transfer_func0, transfer_func], axis=0)
         transfer_func = -1 * transfer_func  # RHS of equation: ma = Σf
         force_impedance = force_from_rao_transfer_function(transfer_func)
 
@@ -564,7 +565,7 @@ def from_impedance(
 
         # wec
         wec = WEC(f1, nfreq, forces, constraints,
-                  inertia_in_forces=True, ndof=shape[0])
+                  inertia_in_forces=True, ndof=shape[1])
         return wec
 
     def residual(self, x_wec: ndarray, x_opt: ndarray, waves: Dataset,
@@ -1481,7 +1482,7 @@ def mimo_transfer_mat(
     For example, it can be an impedance matrix or an RAO transfer
     matrix.
     The input complex impedance matrix has shape
-    :python`(ndof, ndof, nfreq)`.
+    :python`(nfreq, ndof, ndof)`.
 
     Returns the 2D real matrix that transform the state representation
     of the input variable variable to the state representation of the
@@ -1502,20 +1503,20 @@ def mimo_transfer_mat(
     zero_freq
         Whether the first frequency should be zero.
     """
-    ndof = transfer_mat.shape[0]
-    assert transfer_mat.shape[1] == ndof
+    ndof = transfer_mat.shape[1]
+    assert transfer_mat.shape[2] == ndof
     elem = [[None]*ndof for _ in range(ndof)]
     def block(re, im): return np.array([[re, -im], [im, re]])
     for idof in range(ndof):
         for jdof in range(ndof):
             if zero_freq:
-                Zp0 = transfer_mat[idof, jdof, 0]
+                Zp0 = transfer_mat[0, idof, jdof]
                 assert np.all(np.isreal(Zp0))
                 Zp0 = np.real(Zp0)
-                Zp = transfer_mat[idof, jdof, 1:]
+                Zp = transfer_mat[1:, idof, jdof]
             else:
                 Zp0 = [0.0]
-                Zp = transfer_mat[idof, jdof, :]
+                Zp = transfer_mat[:, idof, jdof]
             re = np.real(Zp)
             im = np.imag(Zp)
             blocks = [block(ire, iim) for (ire, iim) in zip(re[:-1], im[:-1])]
@@ -1962,7 +1963,8 @@ def force_from_impedance(
     --------
     force_from_rao_transfer_function,
     """
-    return force_from_rao_transfer_function(impedance*(1j*omega), False)
+    omega0 = np.expand_dims(omega, [1,2])
+    return force_from_rao_transfer_function(impedance*(1j*omega0), False)
 
 
 def force_from_waves(force_coeff: ArrayLike,
@@ -1997,7 +1999,7 @@ def inertia(
     inertia_matrix
         Inertia matrix.
     """
-    omega = np.reshape(frequency(f1, nfreq, False)*2*np.pi, [1,1,-1])
+    omega = np.expand_dims(frequency(f1, nfreq, False)*2*np.pi, [1,2])
     inertia_matrix = np.expand_dims(inertia_matrix, -1)
     rao_transfer_function = -1*omega**2*inertia_matrix + 0j
     inertia_fun = force_from_rao_transfer_function(
@@ -2018,8 +2020,6 @@ def standard_forces(hydro_data: Dataset) -> TForceDict:
     hydro_data
         Linear hydrodynamic data.
     """
-    hydro_data = hydro_data.transpose(
-         "omega", "wave_direction", "radiating_dof", "influenced_dof")
 
     # intrinsic impedance
     w = hydro_data['omega']