13#include <unsupported/Eigen/CXX11/Tensor>
15#include <dune/common/promotiontraits.hh>
31template <
typename Derived,
typename T, auto rank>
33 const std::array<T, rank>& dims) {
34 return Eigen::TensorMap<
const Eigen::TensorFixedSize<
35 const typename Derived::Scalar, Eigen::Sizes<Derived::RowsAtCompileTime, Derived::ColsAtCompileTime>>>(
36 matrix.derived().eval().data(), dims);
47auto dyadic(
const auto& A_ij,
const auto& B_kl) {
48 Eigen::array<Eigen::IndexPair<long>, 0> empty_index_list = {};
49 return A_ij.contract(B_kl, empty_index_list).eval();
59template <
typename ScalarType =
double,
int dim = 3>
61 Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>> idTensor;
62 for (
int i = 0; i < dim; ++i)
63 for (
int j = 0; j < dim; ++j)
64 for (
int k = 0; k < dim; ++k)
65 for (
int l = 0; l < dim; ++l)
66 idTensor(i, j, k, l) = 0.5 * ((i == k) * (j == l) + (i == l) * (j == k));
81template <
typename ScalarType =
double,
int dim = 3>
83 Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>> res;
84 for (
int i = 0; i < dim; ++i)
85 for (
int j = 0; j < dim; ++j)
86 for (
int k = 0; k < dim; ++k)
87 for (
int l = 0; l < dim; ++l)
88 res(i, j, k, l) = 0.5 * (A(i, k) * B(j, l) + A(i, l) * B(j, k));
100template <
typename ScalarType =
double,
int dim = 3>
102 Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>> idTensor;
104 for (
int i = 0; i < dim; ++i)
105 for (
int k = 0; k < dim; ++k)
106 idTensor(i, i, k, k) = 1.0;
121template <
typename AType,
typename BType>
122auto fourthOrderIKJL(
const Eigen::MatrixBase<AType>& A,
const Eigen::MatrixBase<BType>& B) {
123 static_assert(AType::RowsAtCompileTime == BType::RowsAtCompileTime);
124 static_assert(AType::ColsAtCompileTime == BType::ColsAtCompileTime);
125 using ScalarResultType =
typename Dune::PromotionTraits<typename AType::Scalar, typename BType::Scalar>::PromotedType;
126 constexpr int dim = AType::RowsAtCompileTime;
127 Eigen::TensorFixedSize<ScalarResultType, Eigen::Sizes<dim, dim, dim, dim>> res;
128 for (
int i = 0; i < dim; ++i)
129 for (
int j = 0; j < dim; ++j)
130 for (
int k = 0; k < dim; ++k)
131 for (
int l = 0; l < dim; ++l)
132 res(i, j, k, l) = A(i, k) * B(j, l);
144template <
typename ScalarType,
long int dim>
145auto symTwoSlots(
const Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>>& t,
146 const std::array<size_t, 2>& slots) {
147 std::array<size_t, 4> shuffleSlot;
148 std::iota(shuffleSlot.begin(), shuffleSlot.end(), 0);
149 std::swap(shuffleSlot[slots[0]], shuffleSlot[slots[1]]);
150 return (0.5 * (t + t.shuffle(shuffleSlot))).eval();
166constexpr Eigen::Index
toVoigt(Eigen::Index i, Eigen::Index j)
noexcept {
169 if ((i == 1 and j == 2) or (i == 2 and j == 1))
171 if ((i == 0 and j == 2) or (i == 2 and j == 0))
173 if ((i == 0 and j == 1) or (i == 1 and j == 0))
175 assert(i < 3 and j < 3 &&
"For Voigt notation the indices need to be 0,1 or 2.");
176 __builtin_unreachable();
193template <
typename ScalarType =
double>
194Eigen::Matrix<ScalarType, 6, 6>
toVoigt(
const Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>>& ft) {
195 Eigen::Matrix<ScalarType, 6, 6> mat;
196 for (Eigen::Index i = 0; i < 3; ++i)
197 for (Eigen::Index j = 0; j < 3; ++j)
198 for (Eigen::Index k = 0; k < 3; ++k)
199 for (Eigen::Index l = 0; l < 3; ++l)
222template <
typename ST,
int size,
int Options>
223requires(size > 0 and size <= 3)
224auto toVoigt(
const Eigen::Matrix<ST, size, size, Options, size, size>& E,
bool isStrain =
true) {
225 Eigen::Vector<ST, (size * (size + 1)) / 2> EVoigt;
226 EVoigt.template segment<size>(0) = E.diagonal();
228 const ST possibleStrainFactor = isStrain ? 2.0 : 1.0;
229 if constexpr (size == 2)
230 EVoigt(2) = E(0, 1) * possibleStrainFactor;
231 else if constexpr (size == 3) {
232 EVoigt(size) = E(1, 2) * possibleStrainFactor;
233 EVoigt(size + 1) = E(0, 2) * possibleStrainFactor;
234 EVoigt(size + 2) = E(0, 1) * possibleStrainFactor;
254template <
typename ST,
int size>
255requires(size == 1 or size == 3 or size == 6)
256auto fromVoigt(
const Eigen::Vector<ST, size>& EVoigt,
bool isStrain =
true) {
257 constexpr int matrixSize = (-1 +
ct_sqrt(1 + 8 * size)) / 2;
258 Eigen::Matrix<ST, matrixSize, matrixSize> E;
259 E.diagonal() = EVoigt.template head<matrixSize>();
261 const ST possibleStrainFactor = isStrain ? 0.5 : 1.0;
262 if constexpr (matrixSize == 2)
263 E(0, 1) = E(1, 0) = EVoigt(2) * possibleStrainFactor;
264 else if constexpr (matrixSize == 3) {
265 E(2, 1) = E(1, 2) = EVoigt(matrixSize) * possibleStrainFactor;
266 E(2, 0) = E(0, 2) = EVoigt(matrixSize + 1) * possibleStrainFactor;
267 E(1, 0) = E(0, 1) = EVoigt(matrixSize + 2) * possibleStrainFactor;
295 assert(i < 6 &&
"For Voigt notation the indices need to be 0 and 5.");
296 __builtin_unreachable();
311template <
typename ScalarType>
312auto fromVoigt(
const Eigen::Matrix<ScalarType, 6, 6>& CVoigt) {
313 Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>> C;
315 for (
size_t i = 0; i < 6; ++i) {
316 for (
size_t j = 0; j < 6; ++j) {
319 C(firstIndices[0], firstIndices[1], secondIndices[0], secondIndices[1]) = CVoigt(i, j);
335template <
typename Geometry>
337 const auto& referenceElement = Dune::ReferenceElements<double, 2>::general(geometry.type());
338 const auto quadPos0 = referenceElement.position(0, 0);
340 const auto jacobianinvT0 = toEigen(geometry.jacobianInverseTransposed(quadPos0));
341 const auto detJ0 = geometry.integrationElement(quadPos0);
343 auto jaco = (jacobianinvT0).inverse().eval();
344 auto J11 = jaco(0, 0);
345 auto J12 = jaco(0, 1);
346 auto J21 = jaco(1, 0);
347 auto J22 = jaco(1, 1);
351 T0 << J11 * J11, J12 * J12, J11 * J12,
352 J21 * J21, J22 * J22, J21 * J22,
353 2.0 * J11 * J21, 2.0 * J12 * J22, J21 * J12 + J11 * J22;
355 return T0.inverse() * detJ0;
368template <
typename Geometry>
370 const auto& referenceElement = Dune::ReferenceElements<double, 3>::general(geometry.type());
371 const auto quadPos0 = referenceElement.position(0, 0);
373 const auto jacobianinvT0 = toEigen(geometry.jacobianInverseTransposed(quadPos0));
374 const auto detJ0 = geometry.integrationElement(quadPos0);
376 auto jaco = (jacobianinvT0).inverse().eval();
377 auto J11 = jaco(0, 0);
378 auto J12 = jaco(0, 1);
379 auto J13 = jaco(0, 2);
380 auto J21 = jaco(1, 0);
381 auto J22 = jaco(1, 1);
382 auto J23 = jaco(1, 2);
383 auto J31 = jaco(2, 0);
384 auto J32 = jaco(2, 1);
385 auto J33 = jaco(2, 2);
387 Eigen::Matrix<double, 6, 6> T0;
389 T0 << J11 * J11, J12 * J12, J13 * J13, J11 * J12, J11 * J13, J12 * J13,
390 J21 * J21, J22 * J22, J23 * J23, J21 * J22, J21 * J23, J22 * J23,
391 J31 * J31, J32 * J32, J33 * J33, J31 * J32, J31 * J33, J32 * J33,
392 2.0 * J11 * J21, 2.0 * J12 * J22, 2.0 * J13 * J23, J11 * J22 + J21 * J12, J11 * J23 + J21 * J13, J12 * J23 + J22 * J13,
393 2.0 * J11 * J31, 2.0 * J12 * J32, 2.0 * J13 * J33, J11 * J32 + J31 * J12, J11 * J33 + J31 * J13, J12 * J33 + J32 * J13,
394 2.0 * J31 * J21, 2.0 * J32 * J22, 2.0 * J33 * J23, J31 * J22 + J21 * J32, J31 * J23 + J21 * J33, J32 * J23 + J22 * J33;
397 return T0.inverse() * detJ0;
Implementation of math related algorithms.
auto fromVoigt(const Eigen::Vector< ST, size > &EVoigt, bool isStrain=true)
Converts a vector given in Voigt notation to a matrix.
Definition: tensorutils.hh:256
auto symmetricIdentityFourthOrder()
Generates a symmetric identity fourth-order tensor.
Definition: tensorutils.hh:60
auto symmetricFourthOrder(const auto &A, const auto &B)
Generates a symmetric fourth-order tensor based on two second-order input tensors.
Definition: tensorutils.hh:82
constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept
Converts 2D indices to Voigt notation index.
Definition: tensorutils.hh:166
Eigen::Tensor< typename Derived::Scalar, rank > tensorView(const Eigen::EigenBase< Derived > &matrix, const std::array< T, rank > &dims)
View an Eigen matrix as an Eigen Tensor with specified dimensions.
Definition: tensorutils.hh:32
auto fourthOrderIKJL(const Eigen::MatrixBase< AType > &A, const Eigen::MatrixBase< BType > &B)
Computes the IKJL product of two matrices.
Definition: tensorutils.hh:122
auto dyadic(const auto &A_ij, const auto &B_kl)
Computes the dyadic product of two Eigen tensors.
Definition: tensorutils.hh:47
auto identityFourthOrder()
Generates an identity fourth-order tensor.
Definition: tensorutils.hh:101
auto symTwoSlots(const Eigen::TensorFixedSize< ScalarType, Eigen::Sizes< dim, dim, dim, dim > > &t, const std::array< size_t, 2 > &slots)
Creates a symmetric fourth-order tensor in the two specified slots of the input tensor.
Definition: tensorutils.hh:145
Definition: simpleassemblers.hh:22
Eigen::Matrix3d calcTransformationMatrix2D(const Geometry &geometry)
Calculates the 2D transformation matrix.
Definition: tensorutils.hh:336
Eigen::Matrix< double, 6, 6 > calcTransformationMatrix3D(const Geometry &geometry)
Calculates the 3D transformation matrix.
Definition: tensorutils.hh:369
constexpr T ct_sqrt(T x)
Compile-time square root for integer types.
Definition: math.hh:47
Definition: concepts.hh:25