From dcc5e3eea28dc5ef5ee81241a382cb64c933108a Mon Sep 17 00:00:00 2001
From: Manuel Weberndorfer <manuel.weberndorfer@id.ethz.ch>
Date: Tue, 30 Mar 2021 14:10:30 +0200
Subject: [PATCH] use Tensor to define 3D stiffness matrix
---
.../ae108/elements/TimoshenkoBeamElement.h | 66 ++++++++-----------
1 file changed, 29 insertions(+), 37 deletions(-)
diff --git a/elements/src/include/ae108/elements/TimoshenkoBeamElement.h b/elements/src/include/ae108/elements/TimoshenkoBeamElement.h
index b5855739..3a494161 100644
--- a/elements/src/include/ae108/elements/TimoshenkoBeamElement.h
+++ b/elements/src/include/ae108/elements/TimoshenkoBeamElement.h
@@ -79,7 +79,6 @@ template <>
inline Eigen::Matrix<double, 12, 12, Eigen::RowMajor>
stiffness_matrix<double, 3>(const Properties<double, 3> &properties,
const double length) {
-
const auto L = length;
const auto A = properties.area;
const auto E = properties.young_modulus;
@@ -90,45 +89,39 @@ stiffness_matrix<double, 3>(const Properties<double, 3> &properties,
const auto J_x = properties.polar_moment_x;
const auto k_z = properties.shear_correction_factor_z;
- const double Phi_y = 12 * E * I_z * k_y / A / G / L / L;
- const double Phi_z = 12 * E * I_y * k_z / A / G / L / L;
-
- auto K = Eigen::Matrix<double, 12, 12, Eigen::RowMajor>::Zero().eval();
+ const double phi_y = 12 * E * I_z * k_y / A / G / L / L;
+ const double phi_z = 12 * E * I_y * k_z / A / G / L / L;
const auto X = A * E / L;
- const auto Y1 = 12 * E * I_z / (1 + Phi_y) / L / L / L;
- const auto Y2 = 6 * E * I_z / (1 + Phi_y) / L / L;
- const auto Y3 = (4 + Phi_y) * E * I_z / (1 + Phi_y) / L;
- const auto Y4 = (2 - Phi_y) * E * I_z / (1 + Phi_y) / L;
- const auto Z1 = 12 * E * I_y / (1 + Phi_z) / L / L / L;
- const auto Z2 = 6 * E * I_y / (1 + Phi_z) / L / L;
- const auto Z3 = (4 + Phi_z) * E * I_y / (1 + Phi_z) / L;
- const auto Z4 = (2 - Phi_z) * E * I_y / (1 + Phi_z) / L;
+ const auto Y1 = 12 * E * I_z / (1 + phi_y) / L / L / L;
+ const auto Y2 = 6 * E * I_z / (1 + phi_y) / L / L;
+ const auto Y3 = (4 + phi_y) * E * I_z / (1 + phi_y) / L;
+ const auto Y4 = (2 - phi_y) * E * I_z / (1 + phi_y) / L;
+ const auto Z1 = 12 * E * I_y / (1 + phi_z) / L / L / L;
+ const auto Z2 = 6 * E * I_y / (1 + phi_z) / L / L;
+ const auto Z3 = (4 + phi_z) * E * I_y / (1 + phi_z) / L;
+ const auto Z4 = (2 - phi_z) * E * I_y / (1 + phi_z) / L;
const auto S = G * J_x / L;
- K(0, 0) = K(6, 6) = X;
- K(1, 1) = K(7, 7) = Y1;
- K(5, 5) = K(11, 11) = Y3;
- K(2, 2) = K(8, 8) = Z1;
- K(3, 3) = K(9, 9) = S;
- K(4, 4) = K(10, 10) = Z3;
-
- K(6, 0) = -X;
- K(7, 1) = -Y1;
- K(5, 1) = Y2;
- K(11, 1) = Y2;
- K(7, 5) = -Y2;
- K(11, 7) = -Y2;
- K(11, 5) = Y4;
- K(8, 2) = -Z1;
- K(4, 2) = -Z2;
- K(10, 2) = -Z2;
- K(8, 4) = Z2;
- K(10, 8) = Z2;
- K(10, 4) = Z4;
- K(9, 3) = -S;
-
- return properties.weight * K.template selfadjointView<Eigen::Upper>();
+ const auto _ = 0.;
+
+ const tensor::Tensor<double, 12, 12> matrix = {{
+ {{X, 0., 0., 0., 0., 0., -X, 0., 0., 0., 0., 0.}},
+ {{_, Y1, 0., 0., 0., Y2, 0., -Y1, 0., 0., 0., Y2}},
+ {{_, _, Z1, 0., -Z2, 0., 0., 0., -Z1, 0., -Z2, 0.}},
+ {{_, _, _, S, 0., 0., 0., 0., 0., -S, 0., 0.}},
+ {{_, _, _, _, Z3, 0., 0., 0., Z2, 0., Z4, 0.}},
+ {{_, _, _, _, _, Y3, 0. - Y2, 0., 0., 0., Y4}},
+ {{_, _, _, _, _, _, X, 0., 0., 0., 0., 0.}},
+ {{_, _, _, _, _, _, _, Y1, 0., 0., 0., -Y2}},
+ {{_, _, _, _, _, _, _, _, Z1, 0., Z2, 0.}},
+ {{_, _, _, _, _, _, _, _, _, S, 0., 0.}},
+ {{_, _, _, _, _, _, _, _, _, _, Z3, 0.}},
+ {{_, _, _, _, _, _, _, _, _, _, _, Y3}},
+ }};
+
+ return properties.weight * tensor::as_matrix_of_rows(&matrix)
+ .template selfadjointView<Eigen::Upper>();
}
// refer to Cook et. al (2002), "Concepts and applications of Finite Element
@@ -137,7 +130,6 @@ template <>
inline Eigen::Matrix<double, 6, 6, Eigen::RowMajor>
stiffness_matrix<double, 2>(const Properties<double, 2> &properties,
const double length) {
-
const auto L = length;
const auto A = properties.area;
const auto E = properties.young_modulus;
--
GitLab