version 0.4.2
enhancedassumedstrains.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2024 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
21
22namespace Ikarus {
23
24template <typename PreFE, typename FE>
25class EnhancedAssumedStrains;
26
31{
33
34 template <typename PreFE, typename FE>
36};
37
48template <typename PreFE, typename FE>
50{
51public:
55 using LocalView = typename Traits::LocalView;
56 using Geometry = typename Traits::Geometry;
57 using GridView = typename Traits::GridView;
59
60 using SupportedResultTypes = std::tuple<decltype(makeRT<ResultTypes::linearStress>())>;
61
66 explicit EnhancedAssumedStrains(const Pre& pre) { this->setEASType(pre.numberOfParameters); }
67
73 bool isDisplacementBased() const { return easVariant_.isDisplacmentBased(); }
74
80 const auto& easVariant() const { return easVariant_; }
81
87 auto numberOfEASParameters() const { return easVariant_.numberOfEASParameters(); }
88
103 template <template <typename, int, int> class RT>
105 Dune::PriorityTag<2>) const {
107 return underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
108
109 using namespace Dune::Indices;
110 using namespace Dune::DerivativeDirections;
111 auto resultVector = underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
112 if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
114 RTWrapper resultWrapper;
115
116 auto calculateAtContribution = [&]<typename EAST>(const EAST& easFunction) {
117 typename EAST::DType D;
118 calculateDAndLMatrix(easFunction, req, D, L_);
119 const auto ufunc = underlying().displacementFunction(req);
120 const auto disp = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
121 const auto C = underlying().materialTangentFunction(req);
122 const auto alpha = (-D.inverse() * L_ * disp).eval();
123 const auto M = easFunction.calcM(local);
124 const auto CEval = C(local);
125 auto easStress = (CEval * M * alpha).eval();
126 resultWrapper = resultVector.asVec() + easStress;
127 };
128 easVariant_(calculateAtContribution);
129 return resultWrapper;
130 }
131 }
132
139 if (numberOfEASParameters != 0)
141 easVariant_.setEASType(numberOfEASParameters);
142 }
143
144protected:
145 void bindImpl() {
146 assert(underlying().localView().bound());
147 easVariant_.bind(underlying().localView().element().geometry());
148 }
149
150public:
151protected:
152 inline void easApplicabilityCheck() const {
153 const auto& numberOfNodes = underlying().numberOfNodes();
154 assert(not(not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)) and
155 (not isDisplacementBased())) &&
156 "EAS only supported for Q1 or H1 elements");
157 }
158
159 template <typename ScalarType>
161 const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K,
162 const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
163 using namespace Dune::DerivativeDirections;
164 using namespace Dune;
167 return;
168
169 decltype(auto) LMat = [this]() -> decltype(auto) {
170 if constexpr (std::is_same_v<ScalarType, double>)
171 return [this]() -> Eigen::MatrixXd& { return L_; }();
172 else
173 return Eigen::MatrixX<ScalarType>{};
174 }();
175
176 auto calculateMatrixContribution = [&]<typename EAST>(const EAST& easFunction) {
177 typename EAST::DType D;
178 calculateDAndLMatrix(easFunction, par, D, LMat, dx);
179
180 K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
181 K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
182 };
183 easVariant_(calculateMatrixContribution);
184 }
185
186 template <typename ScalarType>
187 inline ScalarType calculateScalarImpl(
188 const Requirement& par, ScalarAffordance affordance,
189 const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
192 return 0.0;
193 DUNE_THROW(Dune::NotImplemented,
194 "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
195 }
196
197 template <typename ScalarType>
199 const Requirement& par, VectorAffordance affordance, typename Traits::template VectorType<ScalarType> force,
200 const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
202 using namespace Dune;
203 using namespace Dune::DerivativeDirections;
204 const auto uFunction = underlying().displacementFunction(par, dx);
205 auto strainFunction = underlying().strainFunction(par, dx);
206 const auto& numberOfNodes = underlying().numberOfNodes();
207
208 auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
209 typename EAST::DType D;
210 calculateDAndLMatrix(easFunction, par, D, L_);
211
212 decltype(auto) LMat = [this]() -> decltype(auto) {
213 if constexpr (std::is_same_v<ScalarType, double>)
214 return [this]() -> Eigen::MatrixXd& { return L_; }();
215 else
216 return Eigen::MatrixX<ScalarType>{};
217 }();
218 const auto disp = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef());
219 const auto alpha = (-D.inverse() * L_ * disp).eval();
220 const auto geo = underlying().localView().element().geometry();
221 auto C = underlying().materialTangentFunction(par);
222
223 for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
224 const auto M = easFunction.calcM(gp.position());
225 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
226 const auto CEval = C(gpIndex);
227 auto stresses = (CEval * M * alpha).eval();
228 for (size_t i = 0; i < numberOfNodes; ++i) {
229 const auto bopI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
230 force.template segment<Traits::worlddim>(Traits::worlddim * i) += bopI.transpose() * stresses * intElement;
231 }
232 }
233 };
234 easVariant_(calculateForceContribution);
235 }
236
237private:
238 EAS::Impl::EASVariant<Geometry> easVariant_;
239 mutable Eigen::MatrixXd L_;
240
241 //> CRTP
242 const auto& underlying() const { return static_cast<const FE&>(*this); }
243 auto& underlying() { return static_cast<FE&>(*this); }
244
245 template <typename ScalarType, int enhancedStrainSize>
246 void calculateDAndLMatrix(
247 const auto& easFunction, const auto& par, Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
248 Eigen::MatrixX<ScalarType>& LMat,
249 const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
250 using namespace Dune;
251 using namespace Dune::DerivativeDirections;
252
253 auto strainFunction = underlying().strainFunction(par, dx);
254 const auto C = underlying().materialTangentFunction(par);
255 const auto geo = underlying().localView().element().geometry();
256 const auto numberOfNodes = underlying().numberOfNodes();
257 DMat.setZero();
258 LMat.setZero(enhancedStrainSize, underlying().localView().size());
259 for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
260 const auto M = easFunction.calcM(gp.position());
261 const auto CEval = C(gpIndex);
262 const double detJTimesW = geo.integrationElement(gp.position()) * gp.weight();
263 DMat += M.transpose() * CEval * M * detJTimesW;
264 for (size_t i = 0U; i < numberOfNodes; ++i) {
265 const size_t I = Traits::worlddim * i;
266 const auto Bi = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
267 LMat.template block<enhancedStrainSize, Traits::worlddim>(0, I) += M.transpose() * CEval * Bi * detJTimesW;
268 }
269 }
270 }
271};
272
278auto eas(int numberOfEASParameters = 0) {
279 EnhancedAssumedStrainsPre pre(numberOfEASParameters);
280
281 return pre;
282}
283
284} // namespace Ikarus
285
286#else
287 #error EnhancedAssumedStrains depends on dune-localfefunctions, which is not included
288#endif
Definition of the LinearElastic class for finite element mechanics computations.
Definition of the EAS variants.
Implementation of the finite element CRTP mixin class.
auto viewAsFlatEigenVector(Dune::BlockVector< ValueType > &blockedVector)
View Dune::BlockVector as an Eigen::Vector.
Definition: linearalgebrahelper.hh:57
Definition: dirichletbcenforcement.hh:6
MatrixAffordance
A strongly typed enum class representing the matrix affordance.
Definition: ferequirements.hh:63
auto eas(int numberOfEASParameters=0)
A helper function to create an enhanced assumed strain pre finite element.
Definition: enhancedassumedstrains.hh:278
VectorAffordance
A strongly typed enum class representing the vector affordance.
Definition: ferequirements.hh:48
ScalarAffordance
A strongly typed enum class representing the scalar affordance.
Definition: ferequirements.hh:37
Definition: utils/dirichletvalues.hh:28
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:252
Container that is used for FE Results. It gives access to the stored value, but can also be used to a...
Definition: feresulttypes.hh:155
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
typename Basis::GridView GridView
Type of the grid view.
Definition: fetraits.hh:45
static constexpr int worlddim
Dimension of the world space.
Definition: fetraits.hh:60
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:50
typename Traits::GridView GridView
Definition: enhancedassumedstrains.hh:57
void calculateMatrixImpl(const Requirement &par, const MatrixAffordance &affordance, typename Traits::template MatrixType<> K, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:160
typename Traits::LocalView LocalView
Definition: enhancedassumedstrains.hh:55
EnhancedAssumedStrains(const Pre &pre)
Constructor for Enhanced Assumed Strains elements.
Definition: enhancedassumedstrains.hh:66
const auto & easVariant() const
Gets the variant representing the type of Enhanced Assumed Strains (EAS).
Definition: enhancedassumedstrains.hh:80
void easApplicabilityCheck() const
Definition: enhancedassumedstrains.hh:152
void calculateVectorImpl(const Requirement &par, VectorAffordance affordance, typename Traits::template VectorType< ScalarType > force, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:198
void bindImpl()
Definition: enhancedassumedstrains.hh:145
bool isDisplacementBased() const
Checks if the element is displacement-based and the EAS is turned off.
Definition: enhancedassumedstrains.hh:73
std::tuple< decltype(makeRT< ResultTypes::linearStress >())> SupportedResultTypes
Definition: enhancedassumedstrains.hh:60
void setEASType(int numberOfEASParameters)
Sets the EAS type for 2D elements.
Definition: enhancedassumedstrains.hh:138
auto numberOfEASParameters() const
Gets the number of EAS parameters based on the current EAS type.
Definition: enhancedassumedstrains.hh:87
ScalarType calculateScalarImpl(const Requirement &par, ScalarAffordance affordance, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Definition: enhancedassumedstrains.hh:187
typename Traits::Geometry Geometry
Definition: enhancedassumedstrains.hh:56
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:104
A PreFE struct for Enhanced Assumed Strains.
Definition: enhancedassumedstrains.hh:31
int numberOfParameters
Definition: enhancedassumedstrains.hh:32
Definition: utils/dirichletvalues.hh:30