12#if HAVE_DUNE_LOCALFEFUNCTIONS
13 #include <dune/localfefunctions/derivativetransformators.hh>
14 #include <dune/localfefunctions/linearAlgebraHelper.hh>
15 #include <dune/localfefunctions/meta.hh>
24template <
typename PreFE,
typename FE>
25class EnhancedAssumedStrains;
34 template <
typename PreFE,
typename FE>
48template <
typename PreFE,
typename FE>
103 template <
template <
typename,
int,
int>
class RT>
105 Dune::PriorityTag<2>)
const {
107 return underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
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;
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);
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;
128 easVariant_(calculateAtContribution);
129 return resultWrapper;
146 assert(underlying().localView().bound());
147 easVariant_.bind(underlying().localView().element().geometry());
153 const auto& numberOfNodes = underlying().numberOfNodes();
156 "EAS only supported for Q1 or H1 elements");
159 template <
typename ScalarType>
162 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
163 using namespace Dune::DerivativeDirections;
164 using namespace Dune;
169 decltype(
auto) LMat = [
this]() ->
decltype(
auto) {
170 if constexpr (std::is_same_v<ScalarType, double>)
171 return [
this]() -> Eigen::MatrixXd& {
return L_; }();
173 return Eigen::MatrixX<ScalarType>{};
176 auto calculateMatrixContribution = [&]<
typename EAST>(
const EAST& easFunction) {
177 typename EAST::DType D;
178 calculateDAndLMatrix(easFunction, par, D, LMat, dx);
180 K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
181 K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
183 easVariant_(calculateMatrixContribution);
186 template <
typename ScalarType>
189 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
193 DUNE_THROW(Dune::NotImplemented,
194 "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
197 template <
typename ScalarType>
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();
208 auto calculateForceContribution = [&]<
typename EAST>(
const EAST& easFunction) {
209 typename EAST::DType D;
210 calculateDAndLMatrix(easFunction, par, D, L_);
212 decltype(
auto) LMat = [
this]() ->
decltype(
auto) {
213 if constexpr (std::is_same_v<ScalarType, double>)
214 return [
this]() -> Eigen::MatrixXd& {
return L_; }();
216 return Eigen::MatrixX<ScalarType>{};
219 const auto alpha = (-D.inverse() * L_ * disp).eval();
220 const auto geo = underlying().localView().element().geometry();
221 auto C = underlying().materialTangentFunction(par);
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;
234 easVariant_(calculateForceContribution);
238 EAS::Impl::EASVariant<Geometry> easVariant_;
239 mutable Eigen::MatrixXd L_;
242 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
243 auto& underlying() {
return static_cast<FE&
>(*this); }
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;
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();
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) {
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;
278auto eas(
int numberOfEASParameters = 0) {
287 #error EnhancedAssumedStrains depends on dune-localfefunctions, which is not included
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