12#if HAVE_DUNE_LOCALFEFUNCTIONS
13 #include <dune/localfefunctions/derivativetransformators.hh>
14 #include <dune/localfefunctions/linearAlgebraHelper.hh>
15 #include <dune/localfefunctions/meta.hh>
26template <
typename PreFE,
typename FE,
typename ES>
27class EnhancedAssumedStrains;
38 template <
typename PreFE,
typename FE>
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>>
70 template <
typename ST>
71 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
73 template <
template <
typename,
int,
int>
class RT>
82 "EAS method is only implemented for the total Lagrangian setting.");
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.");
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>();
134 const auto geo = underlying().localView().element().geometry();
137 auto calculateAtContribution = [&]<
typename EAST>(
const EAST& easFunction) {
138 Eigen::VectorXd alpha;
140 if constexpr (EAST::enhancedStrainSize != 0) {
141 typename EAST::DType D;
142 calculateDAndLMatrix(easFunction, req, D, L_);
143 alpha = -D.inverse() * L_ * disp;
145 const auto enhancedStrain = EnhancedStrainFunction::value(geo, ufunc, local, easFunction, alpha);
146 resultWrapper = rFunction(enhancedStrain);
148 easVariant_(calculateAtContribution);
149 return resultWrapper;
151 DUNE_THROW(Dune::NotImplemented,
"The requested result type is not supported");
175 assert(underlying().localView().bound());
176 easVariant_.
bind(underlying().localView().element().geometry());
190 const std::remove_reference_t<
typename Traits::template VectorType<>>& correction) {
193 auto correctAlpha = [&]<
typename EAST>(
const EAST& easFunction) {
194 if constexpr (EAST::enhancedStrainSize != 0) {
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());
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));
211 easVariant_(correctAlpha);
215 const auto& numberOfNodes = underlying().numberOfNodes();
219 DUNE_THROW(Dune::NotImplemented,
"EAS is only supported for Q1, Q2 and H1 elements");
222 template <
typename ScalarType>
224 typename Traits::template MatrixType<> K,
227 DUNE_THROW(Dune::NotImplemented,
"MatrixAffordance not implemented: " +
toString(affordance));
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();
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);
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();
262 easVariant_(calculateMatrixContribution);
265 template <
typename ScalarType>
270 return underlying().template energyFunction<ScalarType>(par, dx)();
272 DUNE_THROW(Dune::NotImplemented,
273 "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
276 template <
typename ScalarType>
278 typename Traits::template VectorType<ScalarType> force,
281 DUNE_THROW(Dune::NotImplemented,
"VectorAffordance not implemented: " +
toString(affordance));
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);
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);
300 if constexpr (EAST::enhancedStrainSize != 0) {
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;
311 easVariant_(calculateForceContribution);
314 template <
typename MT,
typename BC>
316 if constexpr (std::same_as<MT, NonLinearSolverMessages>) {
317 using NLSState =
typename BC::State;
327 mutable Eigen::MatrixXd L_;
328 Eigen::VectorXd alpha_;
331 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
332 auto& underlying() {
return static_cast<FE&
>(*this); }
337 void initializeState() {
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;
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();
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;
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;
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;
396 easVariant_(calculateRtildeContribution);
407template <
typename ES = EAS::LinearStrain>
408auto eas(
int numberOfInternalVariables = 0) {
417 #error EnhancedAssumedStrains depends on dune-localfefunctions, which is not included
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