version 0.4.8
displacementgradienttransposed.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2026 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
12#pragma once
13
14#include <dune/common/fvector.hh>
15#include <dune/localfefunctions/eigenDuneTransformations.hh>
16#include <dune/localfefunctions/expressions/greenLagrangeStrains.hh>
17
18#include <Eigen/Core>
19
20#include <ikarus/utils/tensorutils.hh>
21
22namespace Ikarus::EAS {
30struct DisplacementGradientTransposed
31{
46 template <typename GEO, typename EAST>
47 static auto computeDisplacementGradient(const GEO& geo, const auto& uFunction,
48 const Dune::FieldVector<double, GEO::mydimension>& gpPos,
49 const EAST& easFunction, const auto& alpha) {
50 using ST = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
51 constexpr int myDim = GEO::mydimension;
52 using MatrixType = Eigen::Matrix<ST, myDim, myDim>;
53 const MatrixType Hc = compatibleDisplacementGradient<GEO>(uFunction, gpPos);
54 const MatrixType Htilde = enhancedDisplacementGradient(geo, uFunction, gpPos, easFunction, alpha);
55 const MatrixType H = Hc + Htilde;
56 return H;
57 }
58
73 template <typename GEO, typename EAST>
74 static auto value(const GEO& geo, const auto& uFunction, const Dune::FieldVector<double, GEO::mydimension>& gpPos,
75 const EAST& easFunction, const auto& alpha) {
76 using ST = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
77 constexpr int myDim = GEO::mydimension;
78 const Eigen::Matrix<ST, myDim, myDim> H = computeDisplacementGradient(geo, uFunction, gpPos, easFunction, alpha);
79 const Eigen::Matrix<ST, myDim, myDim> E = 0.5 * (H + H.transpose() + H.transpose() * H);
80 return toVoigt(E, true).eval();
81 }
82
103 template <int wrtCoeff, typename GEO, typename EAST>
104 static auto firstDerivative(const GEO& geo, const auto& uFunction, const auto& localBasis, const auto& gpIndex,
105 const Dune::FieldVector<double, GEO::mydimension>& gpPos, const EAST& easFunction,
106 const auto& alpha, const int node = sNaN) {
107 using namespace Dune::DerivativeDirections;
108 using namespace Dune;
109 using ST = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
110 constexpr int myDim = GEO::mydimension;
111 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
112 static_assert(myDim == 2 or myDim == 3,
113 "An enhancement of displacement gradient (transposed) can only be performed for the 2D or 3D case.");
114 constexpr int strainSize = myDim * (myDim + 1) / 2;
115 using MatrixType = Eigen::Matrix<ST, myDim, myDim>;
116 const auto centerPos = centerPosition(geo);
117
118 const MatrixType F = computeDeformationGradient(geo, uFunction, gpPos, easFunction, alpha)
119 .transpose()
120 .eval(); // transposed such that the rows are u_{,1} and u_{,2}
121 const typename EAST::HType Harray = easFunction(gpPos);
122 const auto g0 = F.row(0);
123 const auto g1 = F.row(1);
124
125 if constexpr (wrtCoeff == 0) { // E,d
126 const auto gradUdI = uFunction.evaluateDerivative(gpPos, wrt(spatialAll, coeff(node)), on(gridElement));
127 const auto gradUdI0 = uFunction.evaluateDerivative(centerPos, wrt(spatialAll, coeff(node)), on(gridElement));
128 Eigen::Vector<ST, myDim> dN, dN0, dNtilde;
129 dNtilde.setZero();
130 for (const auto i : Dune::range(myDim)) {
131 dN[i] = getDiagonalEntry(gradUdI[i], 0);
132 dN0[i] = getDiagonalEntry(gradUdI0[i], 0);
133 }
134
135 for (const auto i : Dune::range(myDim))
136 for (const auto t : Dune::range(enhancedStrainSize))
137 for (const auto p : Dune::range(myDim))
138 dNtilde[i] += Harray[t](i, p) * alpha[t] * dN0[p];
139
140 dN += dNtilde;
141 Eigen::Matrix<ST, strainSize, myDim> bopI;
142 bopI.row(0) = dN[0] * g0;
143 bopI.row(1) = dN[1] * g1;
144 if constexpr (myDim == 2)
145 bopI.row(2) = dN[1] * g0 + dN[0] * g1;
146 else {
147 const auto g2 = F.row(2);
148 bopI.row(2) = dN[2] * g2;
149 bopI.row(3) = dN[2] * g1 + dN[1] * g2;
150 bopI.row(4) = dN[2] * g0 + dN[0] * g2;
151 bopI.row(5) = dN[1] * g0 + dN[0] * g1;
152 }
153 return bopI.eval();
154 } else if constexpr (wrtCoeff == 1) { // E,a
155 const MatrixType Fc0 = MatrixType::Identity() + compatibleDisplacementGradient<GEO>(uFunction, centerPos);
156 Eigen::Matrix<ST, strainSize, enhancedStrainSize> M;
157 M.setZero(); // zeros returned if Harray is empty
158 for (const auto q : Dune::range(enhancedStrainSize)) {
159 const MatrixType Htilde = Fc0 * Harray[q].transpose();
160 M(0, q) = Htilde.col(0).dot(g0);
161 M(1, q) = Htilde.col(1).dot(g1);
162 if constexpr (myDim == 2)
163 M(2, q) = Htilde.col(1).dot(g0) + Htilde.col(0).dot(g1);
164 else {
165 const auto g2 = F.row(2);
166 M(2, q) = Htilde.col(2).dot(g2);
167 M(3, q) = Htilde.col(2).dot(g1) + Htilde.col(1).dot(g2);
168 M(4, q) = Htilde.col(2).dot(g0) + Htilde.col(0).dot(g2);
169 M(5, q) = Htilde.col(1).dot(g0) + Htilde.col(0).dot(g1);
170 }
171 }
172 return M.eval();
173 } else
174 static_assert(Dune::AlwaysFalse<GEO>::value,
175 "firstDerivative can only be called with wrtCoeff as 0 and 1 indicating first derivative of the "
176 "Green-Lagrange strain w.r.t d and "
177 "alpha, respectively.");
178 }
179
204 template <int wrtCoeff, typename GEO, typename ST, typename EAST>
205 static auto secondDerivative(const GEO& geo, const auto& uFunction, const auto& localBasis, const auto& gpIndex,
206 const Dune::FieldVector<double, GEO::mydimension>& gpPos,
207 const Eigen::Vector<ST, GEO::mydimension*(GEO::mydimension + 1) / 2>& S,
208 const EAST& easFunction, const auto& alpha, const int I = sNaN, const int J = sNaN) {
209 constexpr int myDim = GEO::mydimension;
210 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
211 using namespace Dune::DerivativeDirections;
212 using namespace Dune;
213 using MatrixType = Eigen::Matrix<ST, myDim, myDim>;
214 const auto centerPos = centerPosition(geo);
215 static_assert(myDim == 2 or myDim == 3,
216 "An enhancement of displacement gradient (transposed) can only be performed for the 2D or 3D case.");
217 if constexpr (wrtCoeff == 0) { // E,dd
218 const typename EAST::HType Harray = easFunction(gpPos);
219 const auto gradUdI = uFunction.evaluateDerivative(gpPos, wrt(spatialAll, coeff(I)), on(gridElement));
220 const auto gradUdI0 = uFunction.evaluateDerivative(centerPos, wrt(spatialAll, coeff(I)), on(gridElement));
221 const auto gradUdJ = uFunction.evaluateDerivative(gpPos, wrt(spatialAll, coeff(J)), on(gridElement));
222 const auto gradUdJ0 = uFunction.evaluateDerivative(centerPos, wrt(spatialAll, coeff(J)), on(gridElement));
223
224 Eigen::Vector<ST, myDim> dNI, dNI0;
225 Eigen::Vector<ST, myDim> dNJ, dNJ0;
226 for (const auto i : Dune::range(myDim)) {
227 dNI[i] = getDiagonalEntry(gradUdI[i], 0);
228 dNI0[i] = getDiagonalEntry(gradUdI0[i], 0);
229 dNJ[i] = getDiagonalEntry(gradUdJ[i], 0);
230 dNJ0[i] = getDiagonalEntry(gradUdJ0[i], 0);
231 }
232
233 Eigen::Vector<ST, myDim> dNItilde = Eigen::Vector<ST, myDim>::Zero();
234 Eigen::Vector<ST, myDim> dNJtilde = Eigen::Vector<ST, myDim>::Zero();
235
236 for (const auto t : Dune::range(enhancedStrainSize)) {
237 dNItilde += Harray[t] * (alpha[t] * dNI0);
238 dNJtilde += Harray[t] * (alpha[t] * dNJ0);
239 }
240
241 dNI += dNItilde;
242 dNJ += dNJtilde;
243 ST val{0.0};
244 if constexpr (myDim == 2)
245 val = S[0] * dNI[0] * dNJ[0] + S[1] * dNI[1] * dNJ[1] + S[2] * (dNI[0] * dNJ[1] + dNI[1] * dNJ[0]);
246 else
247 val = S[0] * dNI[0] * dNJ[0] + S[1] * dNI[1] * dNJ[1] + S[2] * dNI[2] * dNJ[2] +
248 S[3] * (dNI[2] * dNJ[1] + dNI[1] * dNJ[2]) + S[4] * (dNI[0] * dNJ[2] + dNI[2] * dNJ[0]) +
249 S[5] * (dNI[0] * dNJ[1] + dNI[1] * dNJ[0]);
250
251 const MatrixType kg = createScaledIdentityMatrix<ST, myDim, myDim>(val);
252 return kg;
253 } else if constexpr (wrtCoeff == 1) { // E,aa
254 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
255 const typename EAST::HType Harray = easFunction(gpPos);
256 auto kg = Eigen::Matrix<ST, enhancedStrainSize, enhancedStrainSize>::Zero().eval();
257 const MatrixType Fc0 = MatrixType::Identity() + compatibleDisplacementGradient<GEO>(uFunction, centerPos);
258
259 for (const auto i : Dune::range(myDim)) {
260 for (const auto j : Dune::range(myDim)) {
261 const auto sij = 0.5 * S[toVoigt<myDim>(i, j)];
262 for (const auto p : Dune::range(enhancedStrainSize)) {
263 const auto Hp = Fc0 * Harray[p].transpose();
264 for (const auto q : Dune::range(enhancedStrainSize)) {
265 const auto Hq = Fc0 * Harray[q].transpose();
266 kg(p, q) += sij * (Hp.col(i).transpose() * Hq.col(j) + Hp.col(j).transpose() * Hq.col(i)).sum();
267 }
268 }
269 }
270 }
271 return kg;
272 } else if constexpr (wrtCoeff == 2) { // E,ad
273 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
274 const typename EAST::HType Harray = easFunction(gpPos);
275 auto kg = Eigen::Matrix<ST, enhancedStrainSize, myDim>::Zero().eval();
276 const auto gradUdI = uFunction.evaluateDerivative(gpIndex, wrt(spatialAll, coeff(I)), on(gridElement));
277 const auto gradUdI0 = uFunction.evaluateDerivative(centerPos, wrt(spatialAll, coeff(I)), on(gridElement));
278 const MatrixType Fc0 = MatrixType::Identity() + compatibleDisplacementGradient<GEO>(uFunction, centerPos);
279 const MatrixType F = computeDeformationGradient(geo, uFunction, gpPos, easFunction, alpha);
280
281 Eigen::Vector<ST, myDim> dN0;
282 for (const auto i : Dune::range(myDim))
283 dN0[i] = getDiagonalEntry(gradUdI0[i], 0);
284
285 Eigen::Vector<ST, myDim> dNtilde = Eigen::Vector<ST, myDim>::Zero();
286 for (const auto t : Dune::range(enhancedStrainSize))
287 dNtilde += Harray[t] * (alpha[t] * dN0);
288
289 std::array<MatrixType, myDim> dNIdT, dNIdT0;
290 for (const auto i : Dune::range(myDim)) {
291 dNIdT[i].setZero();
292 dNIdT0[i].setZero();
293 for (const auto j : Dune::range(myDim)) {
294 dNIdT[i](i, j) = getDiagonalEntry(gradUdI[j], 0) + dNtilde[j];
295 dNIdT0[i](i, j) = getDiagonalEntry(gradUdI0[j], 0);
296 }
297 }
298
299 std::array<std::array<MatrixType, enhancedStrainSize>, myDim> dNXHtilde;
300 for (auto& row : dNXHtilde)
301 std::ranges::fill(row, MatrixType::Zero());
302 for (const auto p : Dune::range(enhancedStrainSize)) {
303 const auto Htilde = Harray[p].transpose();
304 for (const auto q : Dune::range(myDim))
305 dNXHtilde[q][p] += dNIdT0[q] * Htilde;
306 }
307
308 for (const auto i : Dune::range(myDim)) {
309 for (const auto j : Dune::range(myDim)) {
310 const auto sij = 0.5 * S[toVoigt<myDim>(i, j)];
311 for (const auto p : Dune::range(enhancedStrainSize)) {
312 const auto Hp = Fc0 * Harray[p].transpose();
313 for (const auto q : Dune::range(myDim)) {
314 kg(p, q) +=
315 sij * (Hp.col(i).transpose() * dNIdT[q].col(j) + Hp.col(j).transpose() * dNIdT[q].col(i) +
316 dNXHtilde[q][p].col(i).transpose() * F.col(j) + dNXHtilde[q][p].col(j).transpose() * F.col(i))
317 .sum();
318 }
319 }
320 }
321 }
322 return kg;
323 } else
324 static_assert(
325 Dune::AlwaysFalse<GEO>::value,
326 "secondDerivative can only be called with wrtCoeff as 0, 1 and 2 indicating second derivative of the "
327 "Green-Lagrange strain w.r.t d and "
328 "alpha and the mixed derivative, respectively.");
329 }
330
332 static constexpr auto name() { return std::string("Displacement Gradient (Transposed)"); }
333
334private:
335 static constexpr int sNaN = std::numeric_limits<int>::signaling_NaN();
336
337 template <typename GEO>
338 static Dune::FieldVector<double, GEO::mydimension> centerPosition(const GEO& geo) {
339 const auto& referenceElement = Dune::ReferenceElements<double, GEO::mydimension>::general(geo.type());
340 const auto centerPos = referenceElement.position(0, 0);
341 return centerPos;
342 }
343
344 template <typename GEO>
345 static auto compatibleDisplacementGradient(const auto& uFunction,
346 const Dune::FieldVector<double, GEO::mydimension>& gpPos) {
347 using namespace Dune::DerivativeDirections;
348 using namespace Dune;
349 return toEigen(uFunction.evaluateDerivative(gpPos, wrt(spatialAll), on(gridElement))).eval();
350 }
351
352 template <typename GEO, typename EAST>
353 static auto enhancedDisplacementGradient(const GEO& geo, const auto& uFunction,
354 const Dune::FieldVector<double, GEO::mydimension>& gpPos,
355 const EAST& easFunction, const auto& alpha) {
356 using ST = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
357 constexpr int myDim = GEO::mydimension;
358 using MatrixType = Eigen::Matrix<ST, myDim, myDim>;
359 const auto centerPos = centerPosition(geo);
360 const MatrixType Hc0 = compatibleDisplacementGradient<GEO>(uFunction, centerPos);
361 const MatrixType Fc0 = (MatrixType::Identity() + Hc0).eval();
362 const typename EAST::HType Harray = easFunction(gpPos);
363 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
364 typename EAST::AnsatzType Htilde;
365 Htilde.setZero(); // zeros returned if Harray is empty
366 for (const auto p : Dune::range(enhancedStrainSize))
367 Htilde += Harray[p] * alpha[p];
368 const MatrixType Henhanced = (Fc0 * Htilde.transpose()).eval();
369 return Henhanced;
370 }
371
372 template <typename GEO, typename EAST>
373 static auto computeDeformationGradient(const GEO& geo, const auto& uFunction,
374 const Dune::FieldVector<double, GEO::mydimension>& gpPos,
375 const EAST& easFunction, const auto& alpha) {
376 using ST = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
377 constexpr int myDim = GEO::mydimension;
378 using MatrixType = Eigen::Matrix<ST, myDim, myDim>;
379 const MatrixType H = computeDisplacementGradient(geo, uFunction, gpPos, easFunction, alpha);
380 const MatrixType F = (MatrixType::Identity() + H).eval();
381 return F;
382 }
383};
384
385} // namespace Ikarus::EAS
tensorutils.hh
Helper for the Eigen::Tensor types.
Ikarus::toVoigt
constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept
Converts 2D indices to Voigt notation index.
Definition: tensorutils.hh:182
Ikarus::EAS
Definition: easfunctions/displacementgradient.hh:21
Dune
Definition: utils/dirichletvalues.hh:36
Ikarus::EAS::DisplacementGradientTransposed
A struct computing the value, first and second derivatives of the Green-Lagrange strain tensor (in Vo...
Definition: displacementgradienttransposed.hh:31
Ikarus::EAS::DisplacementGradientTransposed::firstDerivative
static auto firstDerivative(const GEO &geo, const auto &uFunction, const auto &localBasis, const auto &gpIndex, const Dune::FieldVector< double, GEO::mydimension > &gpPos, const EAST &easFunction, const auto &alpha, const int node=sNaN)
Compute the first derivative of the Green-Lagrange strain w.r.t d or alpha for a given node and integ...
Definition: displacementgradienttransposed.hh:104
Ikarus::EAS::DisplacementGradientTransposed::secondDerivative
static auto secondDerivative(const GEO &geo, const auto &uFunction, const auto &localBasis, const auto &gpIndex, const Dune::FieldVector< double, GEO::mydimension > &gpPos, const Eigen::Vector< ST, GEO::mydimension *(GEO::mydimension+1)/2 > &S, const EAST &easFunction, const auto &alpha, const int I=sNaN, const int J=sNaN)
Compute the second derivative of the Green-Lagrange strain w.r.t d or alpha for a given node and inte...
Definition: displacementgradienttransposed.hh:205
Ikarus::EAS::DisplacementGradientTransposed::computeDisplacementGradient
static auto computeDisplacementGradient(const GEO &geo, const auto &uFunction, const Dune::FieldVector< double, GEO::mydimension > &gpPos, const EAST &easFunction, const auto &alpha)
Compute the displacement gradient at a given integration point.
Definition: displacementgradienttransposed.hh:47
Ikarus::EAS::DisplacementGradientTransposed::value
static auto value(const GEO &geo, const auto &uFunction, const Dune::FieldVector< double, GEO::mydimension > &gpPos, const EAST &easFunction, const auto &alpha)
Compute the strain vector at a given integration point.
Definition: displacementgradienttransposed.hh:74
Ikarus::EAS::DisplacementGradientTransposed::name
static constexpr auto name()
The name of the strain measure enhanced w.r.t EAS method.
Definition: displacementgradienttransposed.hh:332
Dune::FieldVector
Definition: utils/dirichletvalues.hh:38