version 0.4.4
assumedstress.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
25
26namespace Ikarus {
27
28template <typename PreFE, typename FE, typename ASType>
29class AssumedStress;
30
35template <typename ASType>
37{
39
40 template <typename PreFE, typename FE>
42};
43
55template <typename PreFE, typename FE, typename ASF>
57 : public std::conditional_t<std::same_as<ASF, PS::LinearStress>,
58 ResultTypeBase<ResultTypes::linearStress, ResultTypes::linearStressFull>,
59 ResultTypeBase<ResultTypes::PK2Stress, ResultTypes::PK2StressFull>>
60
61{
62public:
65 using LocalView = typename Traits::LocalView;
66 using Geometry = typename Traits::Geometry;
67 using GridView = typename Traits::GridView;
70
71 template <typename ST>
72 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
73
74 template <template <typename, int, int> class RT>
76
77 static constexpr int myDim = Traits::mydim;
78 static constexpr int strainDim = myDim * (myDim + 1) / 2;
79 using StrainVector = Eigen::Vector<double, strainDim>;
80 using MaterialMatrix = Eigen::Matrix<double, strainDim, strainDim>;
81
86 explicit AssumedStress(const Pre& pre) {
88 "PS method is only implemented for the total Lagrangian setting.");
89 static_assert(not(FE::strainType == StrainTags::linear) or (std::same_as<ASF, PS::LinearStress>),
90 "If FE::strainType is linear, then the assumed stress must also be linear.");
92 }
93
99 const auto& asVariant() const { return asVariant_; }
100
106 auto numberOfInternalVariables() const { return asVariant_.numberOfInternalVariables(); }
107
121 template <template <typename, int, int> class RT>
122 requires(AssumedStress::template supportsResultType<RT>())
124 Dune::PriorityTag<2>) const {
125 if constexpr (isSameResultType<RT, ResultTypes::linearStress> or isSameResultType<RT, ResultTypes::PK2Stress> or
126 isSameResultType<RT, ResultTypes::linearStressFull> or
127 isSameResultType<RT, ResultTypes::PK2StressFull>) {
128 const auto geo = underlying().localView().element().geometry();
129 const auto ufunc = underlying().displacementFunction(req);
130 auto disp = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
131
132 RTWrapperType<RT> resultWrapper{};
133 auto calculateAtContribution = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
134 Eigen::VectorXd beta;
135 beta.setZero(numberOfInternalVariables());
136 const auto& Rtilde = calculateRtilde<double>(req);
137 typename AssumedStressT::HType H;
138 calculateHAndGMatrix(asFunction, req, H, G_);
139 beta = H.inverse() * (G_ * disp);
140
141 const auto SVoigt = AssumedStressFunction::value(geo, local, asFunction, beta);
142
143 if constexpr ((isSameResultType<RT, ResultTypes::PK2StressFull> or
144 isSameResultType<RT, ResultTypes::linearStressFull>) and
145 requires { underlying().material().underlying(); }) {
146 resultWrapper = enlargeStressAnsatz(SVoigt);
147 } else {
148 resultWrapper = SVoigt;
149 }
150 };
151 asVariant_(calculateAtContribution);
152 return resultWrapper;
153 }
154 DUNE_THROW(Dune::NotImplemented, "The requested result type is not supported");
155 }
156
165 initializeState();
166 }
167
173 const auto& internalVariable() const { return beta_; }
174
175protected:
176 void bindImpl() {
177 assert(underlying().localView().bound());
178 asVariant_.bind(underlying().localView().element().geometry());
179 initializeState();
180 }
181
192 const std::remove_reference_t<typename Traits::template VectorType<>>& correction) {
193 using ScalarType = Traits::ctype;
195 auto correctbeta = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
196 if (correction.size() != par.globalSolution().size())
197 DUNE_THROW(Dune::NotImplemented,
198 "Solution vector and correction vector should be of the same size. Check if DBCOption::Full is "
199 "used. The sizes are " +
200 std::to_string(par.globalSolution().size()) + " and " + std::to_string(correction.size()));
201 const auto localdxBlock = Ikarus::FEHelper::localSolutionBlockVector<Traits, Eigen::VectorXd, double>(
202 correction, underlying().localView());
203 const auto localdx = Dune::viewAsFlatEigenVector(localdxBlock);
204
205 typename AssumedStressT::HType H;
206 calculateHAndGMatrix(asFunction, par, H, G_);
207 const auto& Rtilde = calculateRtilde<ScalarType>(par);
208 this->beta_ += H.inverse() * (Rtilde + (G_ * localdx));
209 };
210
211 asVariant_(correctbeta);
212 calculateMaterialInversion();
213 }
214
215 inline void asApplicabilityCheck() const {
216 const auto& numberOfNodes = underlying().numberOfNodes();
217 if (not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)))
218 DUNE_THROW(Dune::NotImplemented, "AssumedStress is only supported for Q1 and H1 elements" +
219 std::to_string(numberOfNodes) + std::to_string(Traits::mydim));
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 using namespace Dune;
227 using namespace Dune::DerivativeDirections;
228
229 if (affordance != MatrixAffordance::stiffness)
230 DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));
232
233 const auto geo = underlying().localView().element().geometry();
234 const auto& strainFunction = underlying().strainFunction(par, dx);
235 const auto& kGFunction = underlying().template geometricStiffnessMatrixFunction<ScalarType>(par, K, dx);
236 const auto numberOfNodes = underlying().numberOfNodes();
237 const auto& localBasis = underlying().localBasis();
238
239 auto calculateMatrixContribution = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
240 for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
241 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
242
243 for (size_t i = 0; i < numberOfNodes; ++i)
244 for (size_t j = 0; j < numberOfNodes; ++j) {
245 const auto E_dd =
246 strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i, j)), along(SVoigt), on(gridElement));
247 kGFunction(E_dd, i, j, gp);
248 }
249 }
250
251 typename AssumedStressT::HType H;
252 calculateHAndGMatrix(asFunction, par, H, G_, dx);
253
254 K.template triangularView<Eigen::Upper>() += G_.transpose() * H.inverse() * G_;
255 K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
256 };
257 asVariant_(calculateMatrixContribution);
258 }
259
260 template <typename ScalarType>
261 inline ScalarType calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
262 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
263 using namespace Dune;
264 using namespace Dune::DerivativeDirections;
266 ScalarType energy = 0.0;
267
268 // This will not work in the context of autodiff because `AutoDiffFE` class doesn't include static condensation of
269 // the internal variable
270 if constexpr (not Concepts::AutodiffScalar<ScalarType>) {
271 const auto geo = underlying().localView().element().geometry();
272 const auto strainFunction = underlying().strainFunction(par, dx);
273
274 auto calculateScalarContribution = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
275 for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
276 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
277 const auto& Es = strainStates_[gpIndex];
278 const auto EVoigt = strainFunction.evaluate(gp.position(), on(gridElement));
279
280 energy += ((SVoigt.transpose() * EVoigt - ScalarType(0.5) * (SVoigt.transpose() * Es)).eval() *
281 geo.integrationElement(gp.position()) * gp.weight())
282 .eval()[0];
283 }
284 };
285
286 asVariant_(calculateScalarContribution);
287 } else {
288 DUNE_THROW(Dune::NotImplemented,
289 "AssumedStress element does not support scalar calculations for autodiff scalars");
290 }
291 return energy;
292 }
293
294 template <typename ScalarType>
296 typename Traits::template VectorType<ScalarType> force,
297 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
298 using namespace Dune::DerivativeDirections;
299 using namespace Dune;
300
301 if (affordance != VectorAffordance::forces)
302 DUNE_THROW(Dune::NotImplemented, "VectorAffordance not implemented: " + toString(affordance));
304
305 const auto geo = underlying().localView().element().geometry();
306 const auto& uFunction = underlying().displacementFunction(par, dx);
307 const auto& strainFunction = underlying().strainFunction(par, dx);
308 const auto numberOfNodes = underlying().numberOfNodes();
309 const auto& localBasis = underlying().localBasis();
310 const auto fIntFunction = underlying().template internalForcesFunction<ScalarType>(par, force, dx);
311
312 auto calculateForceContribution = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
313 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
314 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
315
316 for (size_t i = 0; i < numberOfNodes; ++i) {
317 const auto E_dI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
318 fIntFunction(SVoigt, E_dI, i, gp);
319 }
320 }
321 typename AssumedStressT::HType H;
322 calculateHAndGMatrix(asFunction, par, H, G_);
323 const auto Rtilde = calculateRtilde(par, dx);
324 force += G_.transpose() * H.inverse() * Rtilde;
325 };
326 asVariant_(calculateForceContribution);
327 }
328
329 template <typename MT, typename BC>
330 void subscribeToImpl(BC& bc) {
331 if constexpr (std::same_as<MT, NonLinearSolverMessages>) {
332 using NLSState = typename BC::State;
333 underlying().subscribe(bc, [&](NonLinearSolverMessages message, const NLSState& state) {
335 this->updateStateImpl(state.domain, state.correction);
336 });
337 }
338 }
339
340private:
342 mutable Eigen::MatrixXd G_;
343 Eigen::VectorXd beta_;
344 std::vector<StrainVector> strainStates_;
345 std::vector<MaterialMatrix> invertedMaterialStates_;
346
347 //> CRTP
348 const auto& underlying() const { return static_cast<const FE&>(*this); }
349 auto& underlying() { return static_cast<FE&>(*this); }
350
354 void initializeState() {
355 beta_.resize(numberOfInternalVariables());
356 beta_.setZero();
357 strainStates_.resize(underlying().localBasis().integrationPointSize(), StrainVector::Zero());
358 invertedMaterialStates_.resize(underlying().localBasis().integrationPointSize(), MaterialMatrix::Zero());
359 calculateMaterialInversion();
360 }
361
362 template <typename ScalarType, int assumedStressSize>
363 void calculateHAndGMatrix(const auto& asFunction, const auto& par,
364 Eigen::Matrix<double, assumedStressSize, assumedStressSize>& H,
365 Eigen::MatrixX<ScalarType>& G, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
366 using namespace Dune::DerivativeDirections;
367 using namespace Dune;
368
369 const auto& uFunction = underlying().displacementFunction(par);
370 const auto& strainFunction = underlying().strainFunction(par, dx);
371 const auto geo = underlying().localView().element().geometry();
372 const auto numberOfNodes = underlying().numberOfNodes();
373 const auto& localBasis = underlying().localBasis();
374 H.setZero();
375 G.setZero(assumedStressSize, underlying().localView().size());
376 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
377 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
378
379 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
380 const auto& D = invertedMaterialStates_[gpIndex];
381
382 const auto S_b = AssumedStressFunction::template firstDerivative(geo, uFunction, localBasis, gpIndex,
383 gp.position(), asFunction, beta_);
384 const auto S_bb = AssumedStressFunction::template secondDerivative(geo, uFunction, localBasis, gpIndex,
385 gp.position(), SVoigt, asFunction, beta_);
386
387 H += (S_b.transpose() * D * S_b + S_bb) * intElement;
388
389 for (size_t i = 0U; i < numberOfNodes; ++i) {
390 const auto E_dI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
391
392 G.template block<assumedStressSize, Traits::worlddim>(0, myDim * i) += S_b.transpose() * E_dI * intElement;
393 }
394 }
395 }
396
397 template <typename ScalarType>
398 Eigen::VectorX<ScalarType> calculateRtilde(const Requirement& par,
399 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
400 using namespace Dune;
401 using namespace Dune::DerivativeDirections;
402
403 const auto geo = underlying().localView().element().geometry();
404 const auto& uFunction = underlying().displacementFunction(par, dx);
405 const auto strainFunction = underlying().strainFunction(par, dx);
406 const auto& localBasis = underlying().localBasis();
407
408 Eigen::VectorX<ScalarType> Rtilde;
409 Rtilde.setZero(numberOfInternalVariables());
410
411 auto calculateRtildeContribution = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
412 for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
413 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
414 const auto S_b = AssumedStressFunction::template firstDerivative(geo, uFunction, localBasis, gpIndex,
415 gp.position(), asFunction, beta_);
416
417 const auto& Es = strainStates_[gpIndex];
418 const auto EVoigt = strainFunction.evaluate(gp.position(), on(gridElement));
419
420 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
421 Rtilde += (S_b.transpose() * (EVoigt - Es)).eval() * intElement;
422 }
423 };
424
425 asVariant_(calculateRtildeContribution);
426 return Rtilde;
427 }
428
429 void calculateMaterialInversion() {
430 auto saveMaterialState = [&]<typename AssumedStressT>(const AssumedStressT& asFunction) {
431 const auto geo = underlying().localView().element().geometry();
432
433 for (const auto& [gpIndex, gp] : underlying().localBasis().viewOverIntegrationPoints()) {
434 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
435 const auto& EsOld = strainStates_[gpIndex];
436 std::tie(invertedMaterialStates_[gpIndex], strainStates_[gpIndex]) =
437 underlying().material().template materialInversion<FE::strainType, true>(SVoigt, EsOld);
438 }
439 };
440 asVariant_(saveMaterialState);
441 }
442
443 auto enlargeStressAnsatz(const StrainVector& SVoigt, size_t qpIndexForIteration = 0) const
444 requires(decltype(underlying().material())::isReduced)
445 {
446 auto [_, Es] = underlying().material().template materialInversion<FE::strainType, true>(
447 SVoigt, strainStates_[qpIndexForIteration]);
448 auto Esfull = enlargeIfReduced<decltype(underlying().material())>(Es).eval();
449 return underlying().material().underlying().template stresses<FE::strainType, true>(Esfull);
450 }
451};
452
458template <typename ASType = PS::LinearStress>
459auto assumedStress(int numberOfInternalVariables) {
460 return AssumedStressPre<ASType>(numberOfInternalVariables);
461}
462
463} // namespace Ikarus
464
465#else
466 #error AssumedStress depends on dune-localfefunctions, which is not included
467#endif
Helper for the autodiff library.
Enums for observer messages.
Definition of the LinearElastic class for finite element mechanics computations.
Definition of the AssumedStress variants.
Header file for various assumed stress functions.
Definition of several material related enums.
decltype(auto) enlargeIfReduced(const Eigen::MatrixBase< Derived > &E)
Enlarges a matrix if it reduced in the context of material laws, i.e., VanishingStress If the materia...
Definition: linearalgebrahelper.hh:617
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
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
auto assumedStress(int numberOfInternalVariables)
A helper function to create an assumed stress pre finite element.
Definition: assumedstress.hh:459
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 Assumed Stress with displacement based elements.
Definition: assumedstress.hh:61
const auto & internalVariable() const
Gets the internal state variable beta for the AssumedStress element.
Definition: assumedstress.hh:173
void calculateMatrixImpl(const Requirement &par, const MatrixAffordance &affordance, typename Traits::template MatrixType<> K, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: assumedstress.hh:223
void subscribeToImpl(BC &bc)
Definition: assumedstress.hh:330
typename Traits::LocalView LocalView
Definition: assumedstress.hh:65
AssumedStress(const Pre &pre)
Constructor for Assunmed Stress elements.
Definition: assumedstress.hh:86
static constexpr int myDim
Definition: assumedstress.hh:77
void setAssumedStressType(int numberOfInternalVariables)
Sets the AssumedStress type for 2D elements.
Definition: assumedstress.hh:162
typename Traits::Geometry Geometry
Definition: assumedstress.hh:66
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 Assumed Stress Method.
Definition: assumedstress.hh:123
void updateStateImpl(const Requirement &par, const std::remove_reference_t< typename Traits::template VectorType<> > &correction)
Updates the internal state variable beta_ at the end of an iteration before the update of the displac...
Definition: assumedstress.hh:191
typename Traits::GridView GridView
Definition: assumedstress.hh:67
void calculateVectorImpl(const Requirement &par, VectorAffordance affordance, typename Traits::template VectorType< ScalarType > force, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: assumedstress.hh:295
Eigen::Matrix< double, strainDim, strainDim > MaterialMatrix
Definition: assumedstress.hh:80
FERequirements< FESolutions::displacement, FEParameter::loadfactor > Requirement
Definition: assumedstress.hh:64
void bindImpl()
Definition: assumedstress.hh:176
ASF AssumedStressFunction
Definition: assumedstress.hh:69
ScalarType calculateScalarImpl(const Requirement &par, ScalarAffordance affordance, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: assumedstress.hh:261
std::optional< std::reference_wrapper< const Eigen::VectorX< ST > > > VectorXOptRef
Definition: assumedstress.hh:72
static constexpr int strainDim
Definition: assumedstress.hh:78
auto numberOfInternalVariables() const
Gets the number of AssumedStress parameters based on the current AssumedStress type.
Definition: assumedstress.hh:106
const auto & asVariant() const
Gets the variant representing the type of Assumed Stress.
Definition: assumedstress.hh:99
Eigen::Vector< double, strainDim > StrainVector
Definition: assumedstress.hh:79
void asApplicabilityCheck() const
Definition: assumedstress.hh:215
A PreFE struct for Assumed Stress.
Definition: assumedstress.hh:37
int numberOfParameters
Definition: assumedstress.hh:38
auto numberOfInternalVariables() const
A helper function to get the number of AssumedStress parameters.
Definition: asvariants.hh:63
void bind(const GEO &geometry)
Definition: asvariants.hh:73
void setAssumedStressType(int numberOfInternalVariables)
Definition: asvariants.hh:68
Definition: utils/dirichletvalues.hh:34
Concept to check if the underlying scalar type is a dual type.
Definition: utils/concepts.hh:625
Concept to check if the underlying strain and stress tag correspond to a total Lagrangian formulation...
Definition: utils/concepts.hh:694
Several concepts.