version 0.4.4
enhancedassumedstrains.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@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 ESF>
55 : public std::conditional_t<std::same_as<ESF, 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;
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(
84 not(FE::strainType == StrainTags::linear) or (std::same_as<EnhancedStrainFunction, EAS::LinearStrain>),
85 "If FE::strainType is linear, then enhancedStrain must also be linear.");
87 }
88
94 bool isDisplacementBased() const { return easVariant_.isDisplacmentBased(); }
95
101 const auto& easVariant() const { return easVariant_; }
102
108 auto numberOfInternalVariables() const { return easVariant_.numberOfInternalVariables(); }
109
124 template <template <typename, int, int> class RT>
125 requires(EnhancedAssumedStrains::template supportsResultType<RT>())
127 Dune::PriorityTag<2>) const {
128 if constexpr (isSameResultType<RT, ResultTypes::linearStress> or isSameResultType<RT, ResultTypes::PK2Stress> or
129 isSameResultType<RT, ResultTypes::linearStressFull> or
130 isSameResultType<RT, ResultTypes::PK2StressFull>) {
131 const auto ufunc = underlying().displacementFunction(req);
132 const auto rFunction = underlying().template resultFunction<RT>();
133 auto disp = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
134 const auto geo = underlying().localView().element().geometry();
135
136 RTWrapperType<RT> resultWrapper{};
137 auto calculateAtContribution = [&]<typename EAST>(const EAST& easFunction) {
138 Eigen::VectorXd alpha;
139 alpha.setZero(numberOfInternalVariables());
140 if constexpr (EAST::enhancedStrainSize != 0) {
141 typename EAST::DType D;
142 calculateDAndLMatrix(easFunction, req, D, L_);
143 alpha = -D.inverse() * L_ * disp;
144 }
145 const auto enhancedStrain = EnhancedStrainFunction::value(geo, ufunc, local, easFunction, alpha);
146 resultWrapper = rFunction(enhancedStrain);
147 };
148 easVariant_(calculateAtContribution);
149 return resultWrapper;
150 }
151 DUNE_THROW(Dune::NotImplemented, "The requested result type is not supported");
152 }
153
163 initializeState();
164 }
165
171 const auto& internalVariable() const { return alpha_; }
172
173protected:
174 void bindImpl() {
175 assert(underlying().localView().bound());
176 easVariant_.bind(underlying().localView().element().geometry());
177 initializeState();
178 }
179
190 const std::remove_reference_t<typename Traits::template VectorType<>>& correction) {
191 using ScalarType = Traits::ctype;
193 auto correctAlpha = [&]<typename EAST>(const EAST& easFunction) {
194 if constexpr (EAST::enhancedStrainSize != 0) {
195 if (correction.size() != par.globalSolution().size())
196 DUNE_THROW(Dune::NotImplemented,
197 "Solution vector and correction vector should be of the same size. Check if DBCOption::Full is "
198 "used. The sizes are " +
199 std::to_string(par.globalSolution().size()) + " and " + std::to_string(correction.size()));
200 const auto& Rtilde = calculateRtilde<ScalarType>(par);
201 const auto localdxBlock = Ikarus::FEHelper::localSolutionBlockVector<Traits, Eigen::VectorXd, double>(
202 correction, underlying().localView());
203 const auto localdx = Dune::viewAsFlatEigenVector(localdxBlock);
204 constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
205 Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize> D;
206 calculateDAndLMatrix(easFunction, par, D, L_);
207 this->alpha_ -= D.inverse() * (Rtilde + (L_ * localdx));
208 }
209 };
210
211 easVariant_(correctAlpha);
212 }
213
214 inline void easApplicabilityCheck() const {
215 const auto& numberOfNodes = underlying().numberOfNodes();
216 if (not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 9 and Traits::mydim == 2) or
217 (numberOfNodes == 8 and Traits::mydim == 3)) and
218 (not isDisplacementBased()))
219 DUNE_THROW(Dune::NotImplemented, "EAS is only supported for Q1, Q2 and H1 elements");
220 }
221
222 template <typename ScalarType>
223 void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
224 typename Traits::template MatrixType<> K,
225 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
226 if (affordance != MatrixAffordance::stiffness)
227 DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));
229
230 const auto geo = underlying().localView().element().geometry();
231 const auto& uFunction = underlying().displacementFunction(par, dx);
232 const auto& kMFunction = underlying().template materialStiffnessMatrixFunction<ScalarType>(par, K, dx);
233 const auto& kGFunction = underlying().template geometricStiffnessMatrixFunction<ScalarType>(par, K, dx);
234 const auto numberOfNodes = underlying().numberOfNodes();
235 const auto& localBasis = underlying().localBasis();
236
237 auto calculateMatrixContribution = [&]<typename EAST>(const EAST& easFunction) {
238 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
239 const auto enhancedStrain = EnhancedStrainFunction::value(geo, uFunction, gp.position(), easFunction, alpha_);
240 const auto stresses = underlying().stress(enhancedStrain);
241 for (size_t i = 0; i < numberOfNodes; ++i) {
242 const auto E_dI = EnhancedStrainFunction::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
243 gp.position(), easFunction, alpha_, i);
244 for (size_t j = 0; j < numberOfNodes; ++j) {
245 const auto E_dJ = EnhancedStrainFunction::template firstDerivative<0>(
246 geo, uFunction, localBasis, gpIndex, gp.position(), easFunction, alpha_, j);
247 const auto E_dd = EnhancedStrainFunction::template secondDerivative<0>(
248 geo, uFunction, localBasis, gpIndex, gp.position(), stresses, easFunction, alpha_, i, j);
249 kMFunction(enhancedStrain, E_dI, E_dJ, i, j, gp);
250 kGFunction(E_dd, i, j, gp);
251 }
252 }
253 }
254
255 if constexpr (EAST::enhancedStrainSize != 0) {
256 typename EAST::DType D;
257 calculateDAndLMatrix(easFunction, par, D, L_);
258 K.template triangularView<Eigen::Upper>() -= L_.transpose() * D.inverse() * L_;
259 K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
260 }
261 };
262 easVariant_(calculateMatrixContribution);
263 }
264
265 template <typename ScalarType>
266 inline ScalarType calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
267 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
270 return underlying().template energyFunction<ScalarType>(par, dx)();
271
272 DUNE_THROW(Dune::NotImplemented,
273 "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
274 }
275
276 template <typename ScalarType>
278 typename Traits::template VectorType<ScalarType> force,
279 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
280 if (affordance != VectorAffordance::forces)
281 DUNE_THROW(Dune::NotImplemented, "VectorAffordance not implemented: " + toString(affordance));
283
284 const auto geo = underlying().localView().element().geometry();
285 const auto& uFunction = underlying().displacementFunction(par, dx);
286 const auto numberOfNodes = underlying().numberOfNodes();
287 const auto& localBasis = underlying().localBasis();
288 const auto fIntFunction = underlying().template internalForcesFunction<ScalarType>(par, force, dx);
289
290 auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
291 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
292 const auto enhancedStrain = EnhancedStrainFunction::value(geo, uFunction, gp.position(), easFunction, alpha_);
293 const auto stresses = underlying().stress(enhancedStrain);
294 for (size_t i = 0; i < numberOfNodes; ++i) {
295 const auto E_dI = EnhancedStrainFunction::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
296 gp.position(), easFunction, alpha_, i);
297 fIntFunction(stresses, E_dI, i, gp);
298 }
299 }
300 if constexpr (EAST::enhancedStrainSize != 0) {
301 // Here, L_ and D doesn't see the ScalarType, while Rtilde sees it. This is needed because only then when we use
302 // automatic differentiation w.r.t to the displacements, we get the correct static-condensed stiffness matrix.
303 // If L_ and D sees the ScalarType, then their derivatives w.r.t. the displacements is also performed which is
304 // then not equivalent to the static-condensed stiffness matrix in general.
305 typename EAST::DType D;
306 calculateDAndLMatrix(easFunction, par, D, L_);
307 const auto& Rtilde = calculateRtilde(par, dx);
308 force -= L_.transpose() * D.inverse() * Rtilde;
309 }
310 };
311 easVariant_(calculateForceContribution);
312 }
313
314 template <typename MT, typename BC>
315 void subscribeToImpl(BC& bc) {
316 if constexpr (std::same_as<MT, NonLinearSolverMessages>) {
317 using NLSState = typename BC::State;
318 underlying().subscribe(bc, [&](NonLinearSolverMessages message, const NLSState& state) {
320 this->updateStateImpl(state.domain, state.correction);
321 });
322 }
323 }
324
325private:
327 mutable Eigen::MatrixXd L_;
328 Eigen::VectorXd alpha_;
329
330 //> CRTP
331 const auto& underlying() const { return static_cast<const FE&>(*this); }
332 auto& underlying() { return static_cast<FE&>(*this); }
333
337 void initializeState() {
338 alpha_.resize(numberOfInternalVariables());
339 alpha_.setZero();
340 }
341
342 template <int enhancedStrainSize>
343 void calculateDAndLMatrix(const auto& easFunction, const auto& par,
344 Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
345 Eigen::MatrixX<double>& LMat) const {
346 using namespace Dune;
347 using namespace Dune::DerivativeDirections;
348
349 const auto& uFunction = underlying().displacementFunction(par);
350 const auto geo = underlying().localView().element().geometry();
351 const auto numberOfNodes = underlying().numberOfNodes();
352 const auto& localBasis = underlying().localBasis();
353 DMat.setZero();
354 LMat.setZero(enhancedStrainSize, underlying().localView().size());
355 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
356 const auto enhancedStrain = EnhancedStrainFunction::value(geo, uFunction, gp.position(), easFunction, alpha_);
357 const auto C = underlying().materialTangent(enhancedStrain);
358 const auto stresses = underlying().stress(enhancedStrain);
359 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
360 const auto E_a = EnhancedStrainFunction::template firstDerivative<1>(geo, uFunction, localBasis, gpIndex,
361 gp.position(), easFunction, alpha_);
362 const auto E_aa = EnhancedStrainFunction::template secondDerivative<1>(
363 geo, uFunction, localBasis, gpIndex, gp.position(), stresses, easFunction, alpha_);
364 DMat += (E_a.transpose() * C * E_a + E_aa) * intElement;
365 for (size_t i = 0U; i < numberOfNodes; ++i) {
366 const auto E_dI = EnhancedStrainFunction::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
367 gp.position(), easFunction, alpha_, i);
368 const auto E_ad = EnhancedStrainFunction::template secondDerivative<2>(
369 geo, uFunction, localBasis, gpIndex, gp.position(), stresses, easFunction, alpha_, i);
370 LMat.template block<enhancedStrainSize, myDim>(0, myDim * i) +=
371 (E_a.transpose() * C * E_dI + E_ad) * intElement;
372 }
373 }
374 }
375
376 template <typename ScalarType>
377 Eigen::VectorX<ScalarType> calculateRtilde(const Requirement& par,
378 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
379 const auto geo = underlying().localView().element().geometry();
380 const auto& uFunction = underlying().displacementFunction(par, dx);
381 const auto& localBasis = underlying().localBasis();
382 Eigen::VectorX<ScalarType> Rtilde;
383 Rtilde.setZero(numberOfInternalVariables());
384
385 auto calculateRtildeContribution = [&]<typename EAST>(const EAST& easFunction) {
386 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
387 const auto enhancedStrain = EnhancedStrainFunction::value(geo, uFunction, gp.position(), easFunction, alpha_);
388 const auto E_a = EnhancedStrainFunction::template firstDerivative<1>(geo, uFunction, localBasis, gpIndex,
389 gp.position(), easFunction, alpha_);
390 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
391 auto stresses = underlying().stress(enhancedStrain);
392 Rtilde += (E_a.transpose() * stresses).eval() * intElement;
393 }
394 };
395
396 easVariant_(calculateRtildeContribution);
397 return Rtilde;
398 }
399};
400
407template <typename ES = EAS::LinearStrain>
408auto eas(int numberOfInternalVariables = 0) {
409 EnhancedAssumedStrainsPre<ES> pre(numberOfInternalVariables);
410
411 return pre;
412}
413
414} // namespace Ikarus
415
416#else
417 #error EnhancedAssumedStrains depends on dune-localfefunctions, which is not included
418#endif
Enums for observer messages.
Definition of the LinearElastic class for finite element mechanics computations.
Definition of the EAS variants.
Definition of several material related enums.
auto viewAsFlatEigenVector(Dune::BlockVector< ValueType > &blockedVector)
View Dune::BlockVector as an Eigen::Vector.
Definition: linearalgebrahelper.hh:56
Definition: assemblermanipulatorbuildingblocks.hh:22
MatrixAffordance
A strongly typed enum class representing the matrix affordance.
Definition: ferequirements.hh:65
auto eas(int numberOfInternalVariables=0)
A helper function to create an enhanced assumed strain pre finite element.
Definition: enhancedassumedstrains.hh:408
NonLinearSolverMessages
Enum class defining non-linear solver-related messages.
Definition: broadcastermessages.hh:22
VectorAffordance
A strongly typed enum class representing the vector affordance.
Definition: ferequirements.hh:50
constexpr std::string toString(DBCOption _e)
Definition: dirichletbcenforcement.hh:8
ScalarAffordance
A strongly typed enum class representing the scalar affordance.
Definition: ferequirements.hh:39
Definition: utils/dirichletvalues.hh:32
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:224
const SolutionVectorType & globalSolution() const
Get the global solution vector.
Definition: ferequirements.hh:314
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
std::optional< std::reference_wrapper< const Eigen::VectorX< ST > > > VectorXOptRef
Definition: enhancedassumedstrains.hh:71
FERequirements< FESolutions::displacement, FEParameter::loadfactor > Requirement
Definition: enhancedassumedstrains.hh:61
void calculateMatrixImpl(const Requirement &par, const MatrixAffordance &affordance, typename Traits::template MatrixType<> K, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:223
void easApplicabilityCheck() const
Definition: enhancedassumedstrains.hh:214
void calculateVectorImpl(const Requirement &par, VectorAffordance affordance, typename Traits::template VectorType< ScalarType > force, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:277
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:126
void bindImpl()
Definition: enhancedassumedstrains.hh:174
ScalarType calculateScalarImpl(const Requirement &par, ScalarAffordance affordance, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:266
ESF EnhancedStrainFunction
Definition: enhancedassumedstrains.hh:66
bool isDisplacementBased() const
Checks if the element is displacement-based and the EAS is turned off.
Definition: enhancedassumedstrains.hh:94
void setEASType(int numberOfInternalVariables)
Sets the EAS type for 2D elements.
Definition: enhancedassumedstrains.hh:159
EnhancedAssumedStrains(const Pre &pre)
Constructor for Enhanced Assumed Strains elements.
Definition: enhancedassumedstrains.hh:80
typename Traits::Geometry Geometry
Definition: enhancedassumedstrains.hh:63
const auto & internalVariable() const
Gets the internal state variable alpha for the EAS element.
Definition: enhancedassumedstrains.hh:171
static constexpr int myDim
Definition: enhancedassumedstrains.hh:68
void subscribeToImpl(BC &bc)
Definition: enhancedassumedstrains.hh:315
const auto & easVariant() const
Gets the variant representing the type of Enhanced Assumed Strains (EAS).
Definition: enhancedassumedstrains.hh:101
typename Traits::GridView GridView
Definition: enhancedassumedstrains.hh:64
typename Traits::LocalView LocalView
Definition: enhancedassumedstrains.hh:62
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:189
auto numberOfInternalVariables() const
Gets the number of EAS parameters based on the current EAS type.
Definition: enhancedassumedstrains.hh:108
A PreFE struct for Enhanced Assumed Strains.
Definition: enhancedassumedstrains.hh:35
int numberOfParameters
Definition: enhancedassumedstrains.hh:36
void bind(const GEO &geometry)
Definition: easvariants.hh:95
auto numberOfInternalVariables() const
A helper function to get the number of EAS parameters.
Definition: easvariants.hh:79
bool isDisplacmentBased() const
A helper function to identify if the formulation is purely displacement-based, i.e....
Definition: easvariants.hh:88
void setEASType(int numberOfInternalVariables)
Definition: easvariants.hh:90
Definition: utils/dirichletvalues.hh:34
Concept to check if the underlying strain and stress tag correspond to a total Lagrangian formulation...
Definition: utils/concepts.hh:694
Several concepts.