From 73d1aeda37c64acd653a5d6063af61468f54dee1 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Wed, 21 Aug 2019 17:40:57 +0200
Subject: [PATCH] [MaL] Extract fourth order rotation to KelvinVec.
Previously only used in the linear elastic orthotropic
model, but the transformation itself is more general
and is being used in other (not yet available) code
parts.
---
.../SolidModels/LinearElasticOrthotropic.cpp | 9 +++-
.../SolidModels/LinearElasticOrthotropic.h | 41 -------------------
MathLib/KelvinVector.cpp | 32 +++++++++++++++
MathLib/KelvinVector.h | 7 ++++
4 files changed, 47 insertions(+), 42 deletions(-)
diff --git a/MaterialLib/SolidModels/LinearElasticOrthotropic.cpp b/MaterialLib/SolidModels/LinearElasticOrthotropic.cpp
index d4b0dd3076e..f15f3e45125 100644
--- a/MaterialLib/SolidModels/LinearElasticOrthotropic.cpp
+++ b/MaterialLib/SolidModels/LinearElasticOrthotropic.cpp
@@ -68,7 +68,14 @@ LinearElasticOrthotropic<DisplacementDim>::getElasticTensor(
// clang-format on
KelvinMatrixType<3> const C_ortho = S_ortho.inverse();
- KelvinMatrixType<3> const Q = fourthOrderRotationMatrix(x);
+ auto const Q = [this, &x]() -> KelvinMatrixType<3> {
+ if (!_local_coordinate_system)
+ {
+ return MathLib::KelvinVector::KelvinMatrixType<3>::Identity();
+ }
+ return MathLib::KelvinVector::fourthOrderRotationMatrix(
+ _local_coordinate_system->transformation<3>(x));
+ }();
// Rotate the forth-order tenser in Kelvin mapping with Q*C_ortho*Q^T and
// return the top left corner block of size 4x4 for two-dimensional case or
diff --git a/MaterialLib/SolidModels/LinearElasticOrthotropic.h b/MaterialLib/SolidModels/LinearElasticOrthotropic.h
index 524ce694f4b..13c5b8140b7 100644
--- a/MaterialLib/SolidModels/LinearElasticOrthotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticOrthotropic.h
@@ -172,47 +172,6 @@ public:
MaterialProperties getMaterialProperties() const { return _mp; }
-private:
- /// Rotation matrix of a forth order tensor in 3D.
- MathLib::KelvinVector::KelvinMatrixType<3> fourthOrderRotationMatrix(
- ParameterLib::SpatialPosition const& x) const
- {
- if (!_local_coordinate_system)
- {
- return MathLib::KelvinVector::KelvinMatrixType<3>::Identity();
- }
-
- // 1-based index access for convenience.
- auto Q = [R = _local_coordinate_system->transformation<3>(x)](
- int const i, int const j) { return R(i - 1, j - 1); };
-
- MathLib::KelvinVector::KelvinMatrixType<3> R;
- R << Q(1, 1) * Q(1, 1), Q(1, 2) * Q(1, 2), Q(1, 3) * Q(1, 3),
- std::sqrt(2) * Q(1, 1) * Q(1, 2), std::sqrt(2) * Q(1, 2) * Q(1, 3),
- std::sqrt(2) * Q(1, 1) * Q(1, 3), Q(2, 1) * Q(2, 1),
- Q(2, 2) * Q(2, 2), Q(2, 3) * Q(2, 3),
- std::sqrt(2) * Q(2, 1) * Q(2, 2), std::sqrt(2) * Q(2, 2) * Q(2, 3),
- std::sqrt(2) * Q(2, 1) * Q(2, 3), Q(3, 1) * Q(3, 1),
- Q(3, 2) * Q(3, 2), Q(3, 3) * Q(3, 3),
- std::sqrt(2) * Q(3, 1) * Q(3, 2), std::sqrt(2) * Q(3, 2) * Q(3, 3),
- std::sqrt(2) * Q(3, 1) * Q(3, 3), std::sqrt(2) * Q(1, 1) * Q(2, 1),
- std::sqrt(2) * Q(1, 2) * Q(2, 2), std::sqrt(2) * Q(1, 3) * Q(2, 3),
- Q(1, 1) * Q(2, 2) + Q(1, 2) * Q(2, 1),
- Q(1, 2) * Q(2, 3) + Q(1, 3) * Q(2, 2),
- Q(1, 1) * Q(2, 3) + Q(1, 3) * Q(2, 1),
- std::sqrt(2) * Q(2, 1) * Q(3, 1), std::sqrt(2) * Q(2, 2) * Q(3, 2),
- std::sqrt(2) * Q(2, 3) * Q(3, 3),
- Q(2, 1) * Q(3, 2) + Q(2, 2) * Q(3, 1),
- Q(2, 2) * Q(3, 3) + Q(2, 3) * Q(3, 2),
- Q(2, 1) * Q(3, 3) + Q(2, 3) * Q(3, 1),
- std::sqrt(2) * Q(1, 1) * Q(3, 1), std::sqrt(2) * Q(1, 2) * Q(3, 2),
- std::sqrt(2) * Q(1, 3) * Q(3, 3),
- Q(1, 1) * Q(3, 2) + Q(1, 2) * Q(3, 1),
- Q(1, 2) * Q(3, 3) + Q(1, 3) * Q(3, 2),
- Q(1, 1) * Q(3, 3) + Q(1, 3) * Q(3, 1);
- return R;
- }
-
protected:
MaterialProperties _mp;
boost::optional<ParameterLib::CoordinateSystem> const&
diff --git a/MathLib/KelvinVector.cpp b/MathLib/KelvinVector.cpp
index 575b0334403..eaed17384d4 100644
--- a/MathLib/KelvinVector.cpp
+++ b/MathLib/KelvinVector.cpp
@@ -178,5 +178,37 @@ kelvinVectorToSymmetricTensor(Eigen::Matrix<double,
"or 6, but a vector of size %d was given.",
v.size());
}
+
+KelvinMatrixType<3> fourthOrderRotationMatrix(
+ Eigen::Matrix3d const& transformation)
+{
+ // 1-based index access for convenience.
+ auto Q = [&](int const i, int const j) {
+ return transformation(i - 1, j - 1);
+ };
+
+ MathLib::KelvinVector::KelvinMatrixType<3> R;
+ R << Q(1, 1) * Q(1, 1), Q(1, 2) * Q(1, 2), Q(1, 3) * Q(1, 3),
+ std::sqrt(2) * Q(1, 1) * Q(1, 2), std::sqrt(2) * Q(1, 2) * Q(1, 3),
+ std::sqrt(2) * Q(1, 1) * Q(1, 3), Q(2, 1) * Q(2, 1), Q(2, 2) * Q(2, 2),
+ Q(2, 3) * Q(2, 3), std::sqrt(2) * Q(2, 1) * Q(2, 2),
+ std::sqrt(2) * Q(2, 2) * Q(2, 3), std::sqrt(2) * Q(2, 1) * Q(2, 3),
+ Q(3, 1) * Q(3, 1), Q(3, 2) * Q(3, 2), Q(3, 3) * Q(3, 3),
+ std::sqrt(2) * Q(3, 1) * Q(3, 2), std::sqrt(2) * Q(3, 2) * Q(3, 3),
+ std::sqrt(2) * Q(3, 1) * Q(3, 3), std::sqrt(2) * Q(1, 1) * Q(2, 1),
+ std::sqrt(2) * Q(1, 2) * Q(2, 2), std::sqrt(2) * Q(1, 3) * Q(2, 3),
+ Q(1, 1) * Q(2, 2) + Q(1, 2) * Q(2, 1),
+ Q(1, 2) * Q(2, 3) + Q(1, 3) * Q(2, 2),
+ Q(1, 1) * Q(2, 3) + Q(1, 3) * Q(2, 1), std::sqrt(2) * Q(2, 1) * Q(3, 1),
+ std::sqrt(2) * Q(2, 2) * Q(3, 2), std::sqrt(2) * Q(2, 3) * Q(3, 3),
+ Q(2, 1) * Q(3, 2) + Q(2, 2) * Q(3, 1),
+ Q(2, 2) * Q(3, 3) + Q(2, 3) * Q(3, 2),
+ Q(2, 1) * Q(3, 3) + Q(2, 3) * Q(3, 1), std::sqrt(2) * Q(1, 1) * Q(3, 1),
+ std::sqrt(2) * Q(1, 2) * Q(3, 2), std::sqrt(2) * Q(1, 3) * Q(3, 3),
+ Q(1, 1) * Q(3, 2) + Q(1, 2) * Q(3, 1),
+ Q(1, 2) * Q(3, 3) + Q(1, 3) * Q(3, 2),
+ Q(1, 1) * Q(3, 3) + Q(1, 3) * Q(3, 1);
+ return R;
+}
} // namespace KelvinVector
} // namespace MathLib
diff --git a/MathLib/KelvinVector.h b/MathLib/KelvinVector.h
index e65b1ba1868..fc83db6d785 100644
--- a/MathLib/KelvinVector.h
+++ b/MathLib/KelvinVector.h
@@ -205,6 +205,13 @@ symmetricTensorToKelvinVector(Eigen::MatrixBase<Derived> const& v)
}
return result;
}
+
+/// Rotation tensor for Kelvin mapped vectors and tensors. It is meant to be
+/// used for rotation of stress/strain tensors epsilon:Q and tangent stiffness
+/// tensors Q*C*Q^t.
+KelvinMatrixType<3> fourthOrderRotationMatrix(
+ Eigen::Matrix3d const& transformation);
+
} // namespace KelvinVector
} // namespace MathLib
--
GitLab