Fix bad tensor product

src/pysagas/cfd/oblique_prandtl_meyer.py

@@ -8,6 +8,8 @@
 from scipy.optimize import root_scalar
 from pysagas.utilities import add_sens_data
 from pysagas.cfd.solver import FlowSolver, FlowResults, SensitivityResults
+from hypervehicle.utilities import SensitivityStudy, merge_stls
+import time
 
 
 class OPM(FlowSolver):
@@ -60,6 +62,7 @@ def solve(
                 )
 
             # Iterate over all cells
+            F_trace = []
             net_force = Vector(0, 0, 0)
             net_moment = Vector(0, 0, 0)
             bad = 0
@@ -114,6 +117,7 @@ def solve(
                     method = 0
                     M2, p2, T2 = (flow.M, flow.P, flow.T)
 
+
                 total += 1
 
                 # Save results for this cell
@@ -133,6 +137,7 @@ def solve(
 
                 # Calculate force vector
                 F = cell.n * p2 * cell.A
+                F_trace.append(F)
                 net_force += F
                 net_moment += Vector.from_coordinates(np.cross(r.vec, F.vec))
 
@@ -507,7 +512,7 @@ def oblique_rho2_rho1(M1: float, beta: float, gamma: float = 1.4):
         Peter Jacobs
         """
         M1n = M1 * abs(np.sin(beta))
-        rho2_rho1 = (gamma + 1) * M1n**2 / 2 + (gamma - 1) * M1n**2
+        rho2_rho1 = (gamma + 1) * M1n**2 / (2 + (gamma - 1) * M1n**2)
         return rho2_rho1
 
     @staticmethod

src/pysagas/sensitivity/calculator_vec.py

@@ -33,6 +33,7 @@ def sensitivity_calculator_vec(
             np.cross(r.T, sensitivities[p, :, :].T).T
             + np.cross(cells.dcdp[:, p].T, F.T).T
         )
+        breakpoint()
 
     dFdp = np.sum(sensitivities, axis=-1) / (freestream.q * A_ref)
     dMdp = np.sum(moment_sensitivities, axis=-1) / (freestream.q * A_ref * c_ref)