version 0.4.1
enhancedassumedstrains.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
10#pragma once
11
12#if HAVE_DUNE_LOCALFEFUNCTIONS
13 #include <dune/localfefunctions/derivativetransformators.hh>
14 #include <dune/localfefunctions/linearAlgebraHelper.hh>
15 #include <dune/localfefunctions/meta.hh>
16
23
24namespace Ikarus {
25
26template <typename PreFE, typename FE, typename ES>
27class EnhancedAssumedStrains;
28
33template <typename ES>
35{
37
38 template <typename PreFE, typename FE>
40};
41
53template <typename PreFE, typename FE, typename ES>
55 : public std::conditional_t<std::same_as<ES, EAS::LinearStrain>,
56 ResultTypeBase<ResultTypes::linearStress, ResultTypes::linearStressFull>,
57 ResultTypeBase<ResultTypes::PK2Stress, ResultTypes::PK2StressFull>>
58{
59public:
62 using LocalView = typename Traits::LocalView;
63 using Geometry = typename Traits::Geometry;
64 using GridView = typename Traits::GridView;
66 using EnhancedStrain = ES;
67
68 static constexpr int myDim = Traits::mydim;
69
70 template <typename ST>
71 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
72
73 template <template <typename, int, int> class RT>
75
80 explicit EnhancedAssumedStrains(const Pre& pre) {
82 "EAS method is only implemented for the total Lagrangian setting.");
83 static_assert(not(FE::strainType == StrainTags::linear) or (std::same_as<EnhancedStrain, EAS::LinearStrain>),
84 "If FE::strainType is linear, then enhancedStrain must also be linear.");
86 }
87
93 bool isDisplacementBased() const { return easVariant_.isDisplacmentBased(); }
94
100 const auto& easVariant() const { return easVariant_; }
101
107 auto numberOfEASParameters() const { return easVariant_.numberOfEASParameters(); }
108
123 template <template <typename, int, int> class RT>
124 requires(EnhancedAssumedStrains::template supportsResultType<RT>())
126 Dune::PriorityTag<2>) const {
127 if constexpr (isSameResultType<RT, ResultTypes::linearStress> or isSameResultType<RT, ResultTypes::PK2Stress> or
128 isSameResultType<RT, ResultTypes::linearStressFull> or
129 isSameResultType<RT, ResultTypes::PK2StressFull>) {
130 const auto ufunc = underlying().displacementFunction(req);
131 const auto rFunction = underlying().template resultFunction<RT>();
132 auto disp = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
133 const auto geo = underlying().localView().element().geometry();
134
135 RTWrapperType<RT> resultWrapper{};
136 auto calculateAtContribution = [&]<typename EAST>(const EAST& easFunction) {
137 Eigen::VectorXd alpha;
138 alpha.setZero(numberOfEASParameters());
139 if constexpr (EAST::enhancedStrainSize != 0) {
140 typename EAST::DType D;
141 calculateDAndLMatrix(easFunction, req, D, L_);
142 alpha = -D.inverse() * L_ * disp;
143 }
144 const auto enhancedStrain = EnhancedStrain::value(geo, ufunc, local, easFunction, alpha);
145 resultWrapper = rFunction(enhancedStrain);
146 };
147 easVariant_(calculateAtContribution);
148 return resultWrapper;
149 }
150 DUNE_THROW(Dune::NotImplemented, "The requested result type is not supported");
151 }
152
159 if (numberOfEASParameters != 0)
162 initializeState();
163 }
164
170 const auto& alpha() const { return alpha_; }
171
172protected:
173 void bindImpl() {
174 assert(underlying().localView().bound());
175 easVariant_.bind(underlying().localView().element().geometry());
176 initializeState();
177 }
178
189 const std::remove_reference_t<typename Traits::template VectorType<>>& correction) {
190 using ScalarType = Traits::ctype;
192 auto correctAlpha = [&]<typename EAST>(const EAST& easFunction) {
193 if constexpr (EAST::enhancedStrainSize != 0) {
194 const auto& Rtilde = calculateRtilde<ScalarType>(par);
195 const auto localdxBlock = Ikarus::FEHelper::localSolutionBlockVector<Traits, Eigen::VectorXd, double>(
196 correction, underlying().localView());
197 const auto localdx = Dune::viewAsFlatEigenVector(localdxBlock);
198 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
199 Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize> D;
200 calculateDAndLMatrix(easFunction, par, D, L_);
201 this->alpha_ -= D.inverse() * (Rtilde + (L_ * localdx));
202 }
203 };
204
205 easVariant_(correctAlpha);
206 }
207
208 inline void easApplicabilityCheck() const {
209 const auto& numberOfNodes = underlying().numberOfNodes();
210 if (not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 9 and Traits::mydim == 2) or
211 (numberOfNodes == 8 and Traits::mydim == 3)) and
212 (not isDisplacementBased()))
213 DUNE_THROW(Dune::NotImplemented, "EAS is only supported for Q1, Q2 and H1 elements");
214 }
215
216 template <typename ScalarType>
217 void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
218 typename Traits::template MatrixType<> K,
219 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
220 if (affordance != MatrixAffordance::stiffness)
221 DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));
223
224 const auto geo = underlying().localView().element().geometry();
225 const auto& uFunction = underlying().displacementFunction(par, dx);
226 const auto& kFunction = underlying().template stiffnessMatrixFunction<ScalarType>(par, K, dx);
227 const auto numberOfNodes = underlying().numberOfNodes();
228 const auto& localBasis = underlying().localBasis();
229
230 auto calculateMatrixContribution = [&]<typename EAST>(const EAST& easFunction) {
231 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
232 const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
233 const auto stresses = underlying().stress(enhancedStrain);
234 for (size_t i = 0; i < numberOfNodes; ++i) {
235 const auto E_dI = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
236 gp.position(), easFunction, alpha_, i);
237 for (size_t j = 0; j < numberOfNodes; ++j) {
238 const auto E_dJ = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
239 gp.position(), easFunction, alpha_, j);
240 const auto E_dd = EnhancedStrain::template secondDerivative<0>(
241 geo, uFunction, localBasis, gpIndex, gp.position(), stresses, easFunction, alpha_, i, j);
242 kFunction(enhancedStrain, E_dI, E_dJ, E_dd, i, j, gp);
243 }
244 }
245 }
246
247 if constexpr (EAST::enhancedStrainSize != 0) {
248 typename EAST::DType D;
249 calculateDAndLMatrix(easFunction, par, D, L_);
250 K.template triangularView<Eigen::Upper>() -= L_.transpose() * D.inverse() * L_;
251 K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
252 }
253 };
254 easVariant_(calculateMatrixContribution);
255 }
256
257 template <typename ScalarType>
258 inline ScalarType calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
259 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
262 return underlying().template energyFunction<ScalarType>(par, dx)();
263
264 DUNE_THROW(Dune::NotImplemented,
265 "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
266 }
267
268 template <typename ScalarType>
270 typename Traits::template VectorType<ScalarType> force,
271 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
272 if (affordance != VectorAffordance::forces)
273 DUNE_THROW(Dune::NotImplemented, "VectorAffordance not implemented: " + toString(affordance));
275
276 const auto geo = underlying().localView().element().geometry();
277 const auto& uFunction = underlying().displacementFunction(par, dx);
278 const auto numberOfNodes = underlying().numberOfNodes();
279 const auto& localBasis = underlying().localBasis();
280 const auto fIntFunction = underlying().template internalForcesFunction<ScalarType>(par, force, dx);
281
282 auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
283 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
284 const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
285 for (size_t i = 0; i < numberOfNodes; ++i) {
286 const auto E_dI = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
287 gp.position(), easFunction, alpha_, i);
288 fIntFunction(enhancedStrain, E_dI, i, gp);
289 }
290 }
291 if constexpr (EAST::enhancedStrainSize != 0) {
292 // Here, L_ and D doesn't see the ScalarType, while Rtilde sees it. This is needed because only then when we use
293 // automatic differentiation w.r.t to the displacements, we get the correct static-condensed stiffness matrix.
294 // If L_ and D sees the ScalarType, then their derivatives w.r.t. the displacements is also performed which is
295 // then not equivalent to the static-condensed stiffness matrix in general.
296 typename EAST::DType D;
297 calculateDAndLMatrix(easFunction, par, D, L_);
298 const auto& Rtilde = calculateRtilde(par, dx);
299 force -= L_.transpose() * D.inverse() * Rtilde;
300 }
301 };
302 easVariant_(calculateForceContribution);
303 }
304
305 template <typename MT, typename BC>
306 void subscribeToImpl(BC& bc) {
307 if constexpr (std::same_as<MT, NonLinearSolverMessages>) {
308 using NLSState = typename BC::State;
309 underlying().subscribe(bc, [&](NonLinearSolverMessages message, const NLSState& state) {
311 this->updateStateImpl(state.domain, state.correction);
312 });
313 }
314 }
315
316private:
318 mutable Eigen::MatrixXd L_;
319 Eigen::VectorXd alpha_;
320
321 //> CRTP
322 const auto& underlying() const { return static_cast<const FE&>(*this); }
323 auto& underlying() { return static_cast<FE&>(*this); }
324
328 void initializeState() {
329 alpha_.resize(numberOfEASParameters());
330 alpha_.setZero();
331 }
332
333 template <int enhancedStrainSize>
334 void calculateDAndLMatrix(const auto& easFunction, const auto& par,
335 Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
336 Eigen::MatrixX<double>& LMat) const {
337 using namespace Dune;
338 using namespace Dune::DerivativeDirections;
339
340 const auto& uFunction = underlying().displacementFunction(par);
341 const auto geo = underlying().localView().element().geometry();
342 const auto numberOfNodes = underlying().numberOfNodes();
343 const auto& localBasis = underlying().localBasis();
344 DMat.setZero();
345 LMat.setZero(enhancedStrainSize, underlying().localView().size());
346 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
347 const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
348 const auto C = underlying().materialTangent(enhancedStrain);
349 const auto stresses = underlying().stress(enhancedStrain);
350 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
351 const auto E_a = EnhancedStrain::template firstDerivative<1>(geo, uFunction, localBasis, gpIndex, gp.position(),
352 easFunction, alpha_);
353 const auto E_aa = EnhancedStrain::template secondDerivative<1>(geo, uFunction, localBasis, gpIndex, gp.position(),
354 stresses, easFunction, alpha_);
355 DMat += (E_a.transpose() * C * E_a + E_aa) * intElement;
356 for (size_t i = 0U; i < numberOfNodes; ++i) {
357 const auto E_dI = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
358 gp.position(), easFunction, alpha_, i);
359 const auto E_ad = EnhancedStrain::template secondDerivative<2>(geo, uFunction, localBasis, gpIndex,
360 gp.position(), stresses, easFunction, alpha_, i);
361 LMat.template block<enhancedStrainSize, myDim>(0, myDim * i) +=
362 (E_a.transpose() * C * E_dI + E_ad) * intElement;
363 }
364 }
365 }
366
367 template <typename ScalarType>
368 Eigen::VectorX<ScalarType> calculateRtilde(const Requirement& par,
369 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
370 const auto geo = underlying().localView().element().geometry();
371 const auto& uFunction = underlying().displacementFunction(par, dx);
372 const auto& localBasis = underlying().localBasis();
373 Eigen::VectorX<ScalarType> Rtilde;
374 Rtilde.setZero(numberOfEASParameters());
375
376 auto calculateRtildeContribution = [&]<typename EAST>(const EAST& easFunction) {
377 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
378 const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
379 const auto E_a = EnhancedStrain::template firstDerivative<1>(geo, uFunction, localBasis, gpIndex, gp.position(),
380 easFunction, alpha_);
381 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
382 auto stresses = underlying().stress(enhancedStrain);
383 Rtilde += (E_a.transpose() * stresses).eval() * intElement;
384 }
385 };
386
387 easVariant_(calculateRtildeContribution);
388 return Rtilde;
389 }
390};
391
398template <typename ES = EAS::LinearStrain>
399auto eas(int numberOfEASParameters = 0) {
400 EnhancedAssumedStrainsPre<ES> pre(numberOfEASParameters);
401
402 return pre;
403}
404
405} // namespace Ikarus
406
407#else
408 #error EnhancedAssumedStrains depends on dune-localfefunctions, which is not included
409#endif
Enums for observer messages.
Definition of the LinearElastic class for finite element mechanics computations.
Definition of several material related enums.
Definition of the EAS variants.
NonLinearSolverMessages
Enum class defining non-linear solver-related messages.
Definition: broadcastermessages.hh:22
auto viewAsFlatEigenVector(Dune::BlockVector< ValueType > &blockedVector)
View Dune::BlockVector as an Eigen::Vector.
Definition: linearalgebrahelper.hh:57
Definition: assemblermanipulatorbuildingblocks.hh:22
MatrixAffordance
A strongly typed enum class representing the matrix affordance.
Definition: ferequirements.hh:64
auto eas(int numberOfEASParameters=0)
A helper function to create an enhanced assumed strain pre finite element.
Definition: enhancedassumedstrains.hh:399
VectorAffordance
A strongly typed enum class representing the vector affordance.
Definition: ferequirements.hh:49
constexpr std::string toString(DBCOption _e)
Definition: dirichletbcenforcement.hh:8
ScalarAffordance
A strongly typed enum class representing the scalar affordance.
Definition: ferequirements.hh:38
Definition: utils/dirichletvalues.hh:30
FE class is a base class for all finite elements.
Definition: febase.hh:79
FETraits< BH, useEigenRef, useFlat > Traits
Definition: febase.hh:38
Class representing the requirements for finite element calculations.
Definition: ferequirements.hh:223
Container that is used for FE Results. It gives access to the stored value, but can also be used to a...
Definition: feresulttypes.hh:164
Traits for handling finite elements.
Definition: fetraits.hh:25
typename Basis::LocalView LocalView
Type of the local view.
Definition: fetraits.hh:42
typename Element::Geometry Geometry
Type of the element geometry.
Definition: fetraits.hh:51
double ctype
Type used for coordinates.
Definition: fetraits.hh:57
typename Basis::GridView GridView
Type of the grid view.
Definition: fetraits.hh:45
static constexpr int mydim
Dimension of the geometry.
Definition: fetraits.hh:63
Wrapper class for using Enhanced Assumed Strains (EAS) with displacement based elements.
Definition: enhancedassumedstrains.hh:58
typename Traits::Geometry Geometry
Definition: enhancedassumedstrains.hh:63
typename Traits::LocalView LocalView
Definition: enhancedassumedstrains.hh:62
FERequirements< FESolutions::displacement, FEParameter::loadfactor > Requirement
Definition: enhancedassumedstrains.hh:61
void updateStateImpl(const Requirement &par, const std::remove_reference_t< typename Traits::template VectorType<> > &correction)
Updates the internal state variable alpha_ at the end of an iteration before the update of the displa...
Definition: enhancedassumedstrains.hh:188
void calculateMatrixImpl(const Requirement &par, const MatrixAffordance &affordance, typename Traits::template MatrixType<> K, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:217
void setEASType(int numberOfEASParameters)
Sets the EAS type for 2D elements.
Definition: enhancedassumedstrains.hh:158
bool isDisplacementBased() const
Checks if the element is displacement-based and the EAS is turned off.
Definition: enhancedassumedstrains.hh:93
void bindImpl()
Definition: enhancedassumedstrains.hh:173
const auto & easVariant() const
Gets the variant representing the type of Enhanced Assumed Strains (EAS).
Definition: enhancedassumedstrains.hh:100
void calculateVectorImpl(const Requirement &par, VectorAffordance affordance, typename Traits::template VectorType< ScalarType > force, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:269
const auto & alpha() const
Gets the internal state variable alpha for the EAS element.
Definition: enhancedassumedstrains.hh:170
ScalarType calculateScalarImpl(const Requirement &par, ScalarAffordance affordance, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:258
void subscribeToImpl(BC &bc)
Definition: enhancedassumedstrains.hh:306
std::optional< std::reference_wrapper< const Eigen::VectorX< ST > > > VectorXOptRef
Definition: enhancedassumedstrains.hh:71
auto calculateAtImpl(const Requirement &req, const Dune::FieldVector< double, Traits::mydim > &local, Dune::PriorityTag< 2 >) const
Calculates a requested result at a specific local position using the Enhanced Assumed Strains (EAS) m...
Definition: enhancedassumedstrains.hh:125
typename Traits::GridView GridView
Definition: enhancedassumedstrains.hh:64
ES EnhancedStrain
Definition: enhancedassumedstrains.hh:66
static constexpr int myDim
Definition: enhancedassumedstrains.hh:68
auto numberOfEASParameters() const
Gets the number of EAS parameters based on the current EAS type.
Definition: enhancedassumedstrains.hh:107
void easApplicabilityCheck() const
Definition: enhancedassumedstrains.hh:208
EnhancedAssumedStrains(const Pre &pre)
Constructor for Enhanced Assumed Strains elements.
Definition: enhancedassumedstrains.hh:80
A PreFE struct for Enhanced Assumed Strains.
Definition: enhancedassumedstrains.hh:35
int numberOfParameters
Definition: enhancedassumedstrains.hh:36
void setEASType(int numberOfEASParameters)
Definition: easvariants.hh:73
auto numberOfEASParameters() const
A helper function to get the number of EAS parameters.
Definition: easvariants.hh:62
bool isDisplacmentBased() const
A helper function to identify if the formulation is purely displacement-based, i.e....
Definition: easvariants.hh:71
void bind(const GEO &geometry)
Definition: easvariants.hh:78
Definition: utils/dirichletvalues.hh:32
Concept to check if the underlying strain and stress tag correspond to a total Lagrangian formulation...
Definition: utils/concepts.hh:650
Several concepts.