rework review comments pt1

power_grid_model_c/power_grid_model/include/power_grid_model/common/common.hpp

@@ -56,7 +56,6 @@ using DoubleComplex = std::complex<double>;
 using std::numbers::inv_sqrt3;
 using std::numbers::pi;
 using std::numbers::sqrt3;
-using std::numbers::e;
 
 constexpr DoubleComplex a2{-0.5, -sqrt3 / 2.0};
 constexpr DoubleComplex a{-0.5, sqrt3 / 2.0};

power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp

@@ -59,17 +59,6 @@ class MissingCaseForEnumError : public InvalidArguments {
     }
 };
 
-class UnsupportedInputDescriptionAsymLine : public PowerGridError {
-  public:
-    UnsupportedInputDescriptionAsymLine() {
-        append_msg("Invalid or missing parameters supplied for component asym_line. The following input specifications are allowed");
-        append_msg("3 phase x_matrix, 3 phase r_matrix and 3 phase c_matrix");
-        append_msg("3 phase + neutral x_matrix, 3 phase + neutral r_matrix and 3 phase + neutral c_matrix");
-        append_msg("3 phase x_matrix, 3 phase r_matrix and c1, c0");
-        append_msg("3 phase + neutral x_matrix, 3 phase + neutral r_matrix and c1, c0");
-    }
-};
-
 class ConflictVoltage : public PowerGridError {
   public:
     ConflictVoltage(ID id, ID id1, ID id2, double u1, double u2) {

power_grid_model_c/power_grid_model/include/power_grid_model/common/matrix_utils.hpp

@@ -0,0 +1,17 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <[email protected]>
+
+#pragma once
+
+#include "three_phase_tensor.hpp"
+
+namespace power_grid_model {
+
+inline DoubleComplex average_of_diagonal_of_matrix(const ComplexTensor<asymmetric_t> &matrix) {
+    return (matrix(0,0) + matrix(1,1) + matrix(2,2)) / 3.0;
+}
+
+inline DoubleComplex average_of_off_diagonal_of_matrix(const ComplexTensor<asymmetric_t> &matrix) {
+    return (matrix(0,2) + matrix(1,1) + matrix(2,0)) / 3.0;
+}
+
+} // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp

@@ -62,8 +62,7 @@ template <scalar_value T> class Tensor : public Eigen3Tensor<T> {
     // additional constructors
     explicit Tensor(T const& x) { (*this) << x, 0.0, 0.0, 0.0, x, 0.0, 0.0, 0.0, x; }
     explicit Tensor(T const& s, T const& m) { (*this) << s, m, m, m, s, m, m, m, s; }
-    explicit Tensor(T const& x1, T const& x2, T const& x3, T const& x4, T const& x5, T const& x6) { (*this) << x1, x2, x4, x2, x3, x5, x4, x5, x6; }
-    explicit Tensor(T const& x1, T const& x2, T const& x3, T const& x4, T const& x5, T const& x6, T const& x7, T const& x8, T const& x9) { (*this) << x1, x2, x3, x4, x5, x6, x7, x8, x9; }
+    explicit Tensor(T const& s1, T const& s2, T const& s3, T const& m12, T const& m13, T const& m23) { (*this) << s1, m12, m13, m12, s2, m23, m13, m23, s3; }
     explicit Tensor(Vector<T> const& v) { (*this) << v(0), 0.0, 0.0, 0.0, v(1), 0.0, 0.0, 0.0, v(2); }
     // eigen expression
     template <typename OtherDerived> Tensor(Eigen::ArrayBase<OtherDerived> const& other) : Eigen3Tensor<T>{other} {}
@@ -79,7 +78,7 @@ template <scalar_value T> class Tensor4 : public Eigen4Tensor<T> {
     // additional constructors
     explicit Tensor4(T const& x) { (*this) << x, 0.0, 0.0, 0.0, 0.0, x, 0.0, 0.0, 0.0, 0.0, x, 0.0, 0.0, 0.0, 0.0, x; }
     explicit Tensor4(T const& s, T const& m) { (*this) << s, m, m, m, m, s, m, m, m, m, s, m, m, m, m, s; }
-    explicit Tensor4(T const& x1, T const& x2, T const& x3, T const& x4, T const& x5, T const& x6, T const& x7, T const& x8, T const& x9, T const& x10) { (*this) << x1, x2, x4, x7, x2, x3, x5, x8, x4, x5, x6, x9, x7, x8, x9, x10; }
+    explicit Tensor4(T const& s1, T const& s2, T const& s3, T const& s4, const& m12, T const& m13, T const& m14, T const& m23, T const& m24, T const& m34) { (*this) << s1, m12, m13, m14, m12, s2, m23, m24, m13, m23, s3, m34, m14, m24, m34, s4; }
     explicit Tensor4(Vector<T> const& v) { (*this) << v(0), 0.0, 0.0, 0.0, 0.0, v(1), 0.0, 0.0, 0.0, 0.0, v(2), 0.0, 0.0, 0.0, 0.0, v(3); }
     // eigen expression
     template <typename OtherDerived> Tensor4(Eigen::ArrayBase<OtherDerived> const& other) : Eigen4Tensor<T>{other} {}

power_grid_model_c/power_grid_model/include/power_grid_model/component/asym_line.hpp

@@ -11,6 +11,8 @@
 #include "../calculation_parameters.hpp"
 #include "../common/common.hpp"
 #include "../common/three_phase_tensor.hpp"
+#include "../common/matrix_utils.hpp"
+#include "line_utils.hpp"
 
 namespace power_grid_model {
 
@@ -29,12 +31,6 @@ class AsymLine : public Branch {
         ComplexTensor<asymmetric_t> c_matrix = compute_c_matrix_from_input(asym_line_input);
         ComplexTensor<asymmetric_t> z_series = compute_z_series_from_input(asym_line_input);
 
-        ComplexTensor<asymmetric_t> empty_tensor;
-
-        if (c_matrix.isZero() || z_series.isZero()) {
-            throw UnsupportedInputDescriptionAsymLine();
-        }
-
         y_series = inv(z_series);
         y_shunt = 2 * pi * system_frequency * c_matrix * 1.0i;
     }
@@ -52,34 +48,22 @@ class AsymLine : public Branch {
     ComplexTensor<asymmetric_t> y_series;
     ComplexTensor<asymmetric_t> y_shunt;
 
-    ComplexTensor<asymmetric_t> kron_reduction(const ComplexTensor4& matrix_to_reduce) const {
-        ComplexTensor4 Y = matrix_to_reduce;
-        ComplexTensor<asymmetric_t> Y_aa = ComplexTensor<asymmetric_t>(Y(0,0), Y(1,0), Y(1, 1), Y(2,0), Y(2,1), Y(2,2));
-        ComplexValue<asymmetric_t> Y_ab(Y(0,3), Y(1,3), Y(2,3));
-        Eigen::Array<DoubleComplex, 1, 3> Y_ba;
-        Y_ba << Y(3,0), Y(3,1), Y(3,2);
-        DoubleComplex Y_bb_inv = 1.0 / Y(3,3);
-
-        return Y_aa - ((Y_ab * Y_bb_inv).matrix() * Y_ba.matrix()).array();
-    }
-
     ComplexTensor<asymmetric_t> compute_z_series_from_input(const power_grid_model::AsymLineInput& asym_line_input) {
         ComplexTensor<asymmetric_t> z_series_abc;
         if (is_nan(asym_line_input.r_na) && is_nan(asym_line_input.x_na)) {
-            ComplexTensor<asymmetric_t> r_matrix = ComplexTensor<asymmetric_t>(asym_line_input.r_aa, asym_line_input.r_ba, asym_line_input.r_bb, asym_line_input.r_ca, asym_line_input.r_cb, asym_line_input.r_cc);
-            ComplexTensor<asymmetric_t> x_matrix = ComplexTensor<asymmetric_t>(asym_line_input.x_aa, asym_line_input.x_ba, asym_line_input.x_bb, asym_line_input.x_ca, asym_line_input.x_cb, asym_line_input.x_cc);
+            ComplexTensor<asymmetric_t> r_matrix = ComplexTensor<asymmetric_t>(asym_line_input.r_aa, asym_line_input.r_bb, asym_line_input.r_cc, asym_line_input.r_ba, asym_line_input.r_ca, asym_line_input.r_cb);
+            ComplexTensor<asymmetric_t> x_matrix = ComplexTensor<asymmetric_t>(asym_line_input.x_aa, asym_line_input.x_bb, asym_line_input.x_cc, asym_line_input.x_ba, asym_line_input.x_ca, asym_line_input.x_cb);
             z_series_abc = r_matrix + x_matrix * 1.0i;
         } 
         else {
-            ComplexTensor4 r_matrix = ComplexTensor4(asym_line_input.r_aa, asym_line_input.r_ba, asym_line_input.r_bb, asym_line_input.r_ca, asym_line_input.r_cb, asym_line_input.r_cc, asym_line_input.r_na, asym_line_input.r_nb, asym_line_input.r_nc, asym_line_input.r_nn);
-            ComplexTensor4 x_matrix = ComplexTensor4(asym_line_input.x_aa, asym_line_input.x_ba, asym_line_input.x_bb, asym_line_input.x_ca, asym_line_input.x_cb, asym_line_input.x_cc, asym_line_input.x_na, asym_line_input.x_nb, asym_line_input.x_nc, asym_line_input.x_nn);
+            ComplexTensor4 r_matrix = ComplexTensor4(asym_line_input.r_aa, asym_line_input.r_bb, asym_line_input.r_cc, asym_line_input.r_nn, asym_line_input.r_ba, asym_line_input.r_ca, asym_line_input.r_na, asym_line_input.r_cb, asym_line_input.r_nb, asym_line_input.r_nc);
+            ComplexTensor4 x_matrix = ComplexTensor4(asym_line_input.x_aa, asym_line_input.x_bb, asym_line_input.x_cc, asym_line_input.x_nn, asym_line_input.x_ba, asym_line_input.x_ca, asym_line_input.x_na, asym_line_input.x_cb, asym_line_input.x_nb, asym_line_input.x_nc);
 
             ComplexTensor4 y = r_matrix + 1.0i * x_matrix;
             z_series_abc = kron_reduction(y);
         }
-        DoubleComplex a = std::pow(e, 1.0i * (2.0 / 3.0) * pi);
-        ComplexTensor<asymmetric_t> a_matrix = ComplexTensor<asymmetric_t>(1, 1, pow(a, 2), 1, a, pow(a, 2));
-        ComplexTensor<asymmetric_t> a_matrix_inv = (1.0/3.0) * ComplexTensor<asymmetric_t>(1, 1, a, 1, pow(a, 2), a);
+        ComplexTensor<asymmetric_t> a_matrix = get_sym_matrix();
+        ComplexTensor<asymmetric_t> a_matrix_inv = get_sym_matrix_inv();
         ComplexTensor<asymmetric_t> z_series = (a_matrix_inv.matrix() * z_series_abc.matrix() * a_matrix.matrix()).array();
         return z_series;
     }
@@ -94,19 +78,11 @@ class AsymLine : public Branch {
         } 
         else {
             ComplexTensor4 c_matrix_neutral = ComplexTensor4(asym_line_input.c_aa, asym_line_input.c_ba, asym_line_input.c_bb, asym_line_input.c_ca, asym_line_input.c_cb, asym_line_input.c_cc, asym_line_input.c_na, asym_line_input.c_nb, asym_line_input.c_nc, asym_line_input.c_nn);
-            c_matrix = this->kron_reduction(c_matrix_neutral);
+            c_matrix = kron_reduction(c_matrix_neutral);
         }
         return c_matrix;
     }
 
-    DoubleComplex average_of_diagonal_of_matrix(const ComplexTensor<asymmetric_t> &matrix) const {
-        return (matrix(0,0) + matrix(1,1) + matrix(2,2)) / 3.0;
-    }
-
-    DoubleComplex average_of_off_diagonal_of_matrix(const ComplexTensor<asymmetric_t> &matrix) const {
-        return (matrix(0,2) + matrix(1,1) + matrix(2,0)) / 3.0;
-    }
-
     BranchCalcParam<symmetric_t> sym_calc_param() const override {
         DoubleComplex y1_series_ = average_of_diagonal_of_matrix(y_series) - average_of_off_diagonal_of_matrix(y_series);
         DoubleComplex y1_shunt_ = average_of_diagonal_of_matrix(y_shunt) - average_of_off_diagonal_of_matrix(y_shunt);

power_grid_model_c/power_grid_model/include/power_grid_model/component/line_utils.hpp

@@ -0,0 +1,18 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <[email protected]>
+
+#pragma once
+
+#include "../common/three_phase_tensor.hpp"
+
+namespace power_grid_model {
+
+inline ComplexTensor<asymmetric_t> kron_reduction(const ComplexTensor4& matrix_to_reduce) {
+    ComplexTensor4 Y = matrix_to_reduce;
+    ComplexTensor<asymmetric_t> Y_aa = ComplexTensor<asymmetric_t>(Y(0,0), Y(1,0), Y(1, 1), Y(2,0), Y(2,1), Y(2,2));
+    ComplexValue<asymmetric_t> Y_ab(Y(0,3), Y(1,3), Y(2,3));
+    ComplexValue<asymmetric_t> Y_ba(Y(3,0), Y(3,1), Y(3,2));
+    DoubleComplex Y_bb_inv = 1.0 / Y(3,3);
+    return Y_aa - vector_outer_product(Y_ba, Y_ab) * Y_bb_inv;
+}
+
+} // namespace power_grid_model