12#if HAVE_DUNE_LOCALFEFUNCTIONS
13 #include <dune/localfefunctions/derivativetransformators.hh>
14 #include <dune/localfefunctions/linearAlgebraHelper.hh>
15 #include <dune/localfefunctions/meta.hh>
28template <
typename PreFE,
typename FE,
typename ASType>
35template <
typename ASType>
40 template <
typename PreFE,
typename FE>
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>>
71 template <
typename ST>
72 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
74 template <
template <
typename,
int,
int>
class RT>
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.");
99 const auto&
asVariant()
const {
return asVariant_; }
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);
133 auto calculateAtContribution = [&]<
typename AssumedStressT>(
const AssumedStressT& asFunction) {
134 Eigen::VectorXd beta;
136 const auto& Rtilde = calculateRtilde<double>(req);
137 typename AssumedStressT::HType H;
138 calculateHAndGMatrix(asFunction, req, H, G_);
139 beta = H.inverse() * (G_ * disp);
141 const auto SVoigt = AssumedStressFunction::value(geo, local, asFunction, beta);
143 if constexpr ((isSameResultType<RT, ResultTypes::PK2StressFull> or
144 isSameResultType<RT, ResultTypes::linearStressFull>) and
145 requires { underlying().material().underlying(); }) {
146 resultWrapper = enlargeStressAnsatz(SVoigt);
148 resultWrapper = SVoigt;
151 asVariant_(calculateAtContribution);
152 return resultWrapper;
154 DUNE_THROW(Dune::NotImplemented,
"The requested result type is not supported");
177 assert(underlying().localView().bound());
178 asVariant_.
bind(underlying().localView().element().geometry());
192 const std::remove_reference_t<
typename Traits::template VectorType<>>& correction) {
195 auto correctbeta = [&]<
typename AssumedStressT>(
const AssumedStressT& asFunction) {
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());
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));
211 asVariant_(correctbeta);
212 calculateMaterialInversion();
216 const auto& numberOfNodes = underlying().numberOfNodes();
218 DUNE_THROW(Dune::NotImplemented,
"AssumedStress is only supported for Q1 and H1 elements" +
219 std::to_string(numberOfNodes) + std::to_string(
Traits::mydim));
222 template <
typename ScalarType>
224 typename Traits::template MatrixType<> K,
226 using namespace Dune;
227 using namespace Dune::DerivativeDirections;
230 DUNE_THROW(Dune::NotImplemented,
"MatrixAffordance not implemented: " +
toString(affordance));
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();
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_);
243 for (
size_t i = 0; i < numberOfNodes; ++i)
244 for (
size_t j = 0; j < numberOfNodes; ++j) {
246 strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i, j)), along(SVoigt), on(gridElement));
247 kGFunction(E_dd, i, j, gp);
251 typename AssumedStressT::HType H;
252 calculateHAndGMatrix(asFunction, par, H, G_, dx);
254 K.template triangularView<Eigen::Upper>() += G_.transpose() * H.inverse() * G_;
255 K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
257 asVariant_(calculateMatrixContribution);
260 template <
typename ScalarType>
263 using namespace Dune;
264 using namespace Dune::DerivativeDirections;
266 ScalarType energy = 0.0;
271 const auto geo = underlying().localView().element().geometry();
272 const auto strainFunction = underlying().strainFunction(par, dx);
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));
280 energy += ((SVoigt.transpose() * EVoigt - ScalarType(0.5) * (SVoigt.transpose() * Es)).eval() *
281 geo.integrationElement(gp.position()) * gp.weight())
286 asVariant_(calculateScalarContribution);
288 DUNE_THROW(Dune::NotImplemented,
289 "AssumedStress element does not support scalar calculations for autodiff scalars");
294 template <
typename ScalarType>
296 typename Traits::template VectorType<ScalarType> force,
298 using namespace Dune::DerivativeDirections;
299 using namespace Dune;
302 DUNE_THROW(Dune::NotImplemented,
"VectorAffordance not implemented: " +
toString(affordance));
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);
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_);
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);
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;
326 asVariant_(calculateForceContribution);
329 template <
typename MT,
typename BC>
331 if constexpr (std::same_as<MT, NonLinearSolverMessages>) {
332 using NLSState =
typename BC::State;
342 mutable Eigen::MatrixXd G_;
343 Eigen::VectorXd beta_;
344 std::vector<StrainVector> strainStates_;
345 std::vector<MaterialMatrix> invertedMaterialStates_;
348 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
349 auto& underlying() {
return static_cast<FE&
>(*this); }
354 void initializeState() {
357 strainStates_.resize(underlying().localBasis().integrationPointSize(), StrainVector::Zero());
358 invertedMaterialStates_.resize(underlying().localBasis().integrationPointSize(), MaterialMatrix::Zero());
359 calculateMaterialInversion();
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;
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();
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();
379 const auto SVoigt = AssumedStressFunction::value(geo, gp.position(), asFunction, beta_);
380 const auto& D = invertedMaterialStates_[gpIndex];
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_);
387 H += (S_b.transpose() * D * S_b + S_bb) * intElement;
389 for (
size_t i = 0U; i < numberOfNodes; ++i) {
390 const auto E_dI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
392 G.template block<assumedStressSize, Traits::worlddim>(0,
myDim * i) += S_b.transpose() * E_dI * intElement;
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;
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();
408 Eigen::VectorX<ScalarType> Rtilde;
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_);
417 const auto& Es = strainStates_[gpIndex];
418 const auto EVoigt = strainFunction.evaluate(gp.position(), on(gridElement));
420 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
421 Rtilde += (S_b.transpose() * (EVoigt - Es)).eval() * intElement;
425 asVariant_(calculateRtildeContribution);
429 void calculateMaterialInversion() {
430 auto saveMaterialState = [&]<
typename AssumedStressT>(
const AssumedStressT& asFunction) {
431 const auto geo = underlying().localView().element().geometry();
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);
440 asVariant_(saveMaterialState);
443 auto enlargeStressAnsatz(
const StrainVector& SVoigt,
size_t qpIndexForIteration = 0) const
444 requires(decltype(underlying().material())::isReduced)
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);
458template <
typename ASType = PS::LinearStress>
466 #error AssumedStress depends on dune-localfefunctions, which is not included
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