version 0.4.4
linearelastic.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
14 #include <optional>
15 #include <type_traits>
16
17 #include <dune/geometry/quadraturerules.hh>
18 #include <dune/localfefunctions/expressions/linearStrainsExpr.hh>
19 #include <dune/localfefunctions/impl/standardLocalFunction.hh>
20
28
29namespace Ikarus {
30
31template <typename PreFE, typename FE, typename PRE>
32class LinearElastic;
33
38template <Concepts::GeometricallyLinearMaterial MAT>
40{
41 using Material = MAT;
43
44 template <typename PreFE, typename FE>
46};
47
56template <typename PreFE, typename FE, typename PRE>
57class LinearElastic : public ResultTypeBase<ResultTypes::linearStress, ResultTypes::linearStressFull>
58{
59public:
62 using FlatBasis = typename Traits::FlatBasis;
64 using LocalView = typename Traits::LocalView;
65 using Geometry = typename Traits::Geometry;
66 using GridView = typename Traits::GridView;
67 using Element = typename Traits::Element;
68 using Material = PRE::Material;
69 using Pre = PRE;
70
71 using LocalBasisType = decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
72
73 template <typename ST>
74 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
75
76 static constexpr int myDim = Traits::mydim;
77 static constexpr int strainDim = myDim * (myDim + 1) / 2;
78 static constexpr auto strainType = StrainTags::linear;
79 static constexpr auto stressType = StressTags::linear;
80
81 template <template <typename, int, int> class RT>
83
84 template <typename ST = double>
85 using StrainType = Eigen::Vector<ST, strainDim>;
86
87 template <typename ST = double>
88 using BopType = Eigen::Matrix<ST, strainDim, myDim>;
89
90 template <typename ST = double>
91 using KgType = Eigen::Matrix<ST, myDim, myDim>;
92
97 explicit LinearElastic(const Pre& pre)
98 : mat_{pre.material} {}
99
100protected:
104 void bindImpl() {
105 const auto& localView = underlying().localView();
106 const auto& element = localView.element();
107 auto& firstChild = localView.tree().child(0);
108 const auto& fe = firstChild.finiteElement();
109 geo_ = std::make_shared<const Geometry>(element.geometry());
110 numberOfNodes_ = fe.size();
111 order_ = 2 * (fe.localBasis().order());
112 localBasis_ = Dune::CachedLocalBasis(fe.localBasis());
113 if constexpr (requires { element.impl().getQuadratureRule(order_); })
114 if (element.impl().isTrimmed())
115 localBasis_.bind(element.impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1));
116 else
117 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_),
118 Dune::bindDerivatives(0, 1));
119 else
120 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_), Dune::bindDerivatives(0, 1));
121 }
122
123public:
132 template <typename ScalarType = double>
133 auto displacementFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
134 const auto& d = par.globalSolution();
135 auto disp = Ikarus::FEHelper::localSolutionBlockVector<Traits>(d, underlying().localView(), dx);
136 Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_);
137 return uFunction;
138 }
148 template <class ScalarType = double>
149 auto strainFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
150 return Dune::linearStrains(displacementFunction(par, dx));
151 }
152
161 template <typename ScalarType, int strainDim>
162 auto internalEnergy(const Eigen::Vector<ScalarType, strainDim>& strain) const {
163 const auto C = materialTangent<ScalarType, strainDim>(strain);
164 return (0.5 * strain.dot(C * strain));
165 }
166
176 template <typename ScalarType, int strainDim, bool voigt = true>
177 auto stress(const Eigen::Vector<ScalarType, strainDim>& strain) const {
178 const auto C = materialTangent<ScalarType, strainDim, voigt>(strain);
179 return (C * strain).eval();
180 }
181
191 template <typename ScalarType, int strainDim, bool voigt = true>
193 const Eigen::Vector<ScalarType, strainDim>& strain = Eigen::Vector<ScalarType, strainDim>::Zero()) const {
194 // Since that material is independent of the strains, a zero strain is passed here
195 return material<ScalarType>().template tangentModuli<strainType, voigt>(strain);
196 }
197
198 const Geometry& geometry() const { return *geo_; }
199 [[nodiscard]] size_t numberOfNodes() const { return numberOfNodes_; }
200 [[nodiscard]] int order() const { return order_; }
201 const Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>>& localBasis() const { return localBasis_; }
202
203 template <typename ScalarType = double>
204 decltype(auto) material() const {
206 return mat_.template rebind<ScalarType>();
207 else
208 return mat_;
209 }
210
211public:
217 template <template <typename, int, int> class RT>
218 requires(LinearElastic::template supportsResultType<RT>())
219 auto resultFunction() const {
220 return [&](const Eigen::Vector<double, strainDim>& strainInVoigt) {
221 if constexpr (isSameResultType<RT, ResultTypes::linearStress> or
222 isSameResultType<RT, ResultTypes::linearStressFull>) {
223 decltype(auto) mat = [&]() {
224 if constexpr (isSameResultType<RT, ResultTypes::linearStressFull> and requires { mat_.underlying(); })
225 return mat_.underlying();
226 else
227 return mat_;
228 }();
229 return RTWrapperType<RT>{mat.template stresses<StrainTags::linear>(enlargeIfReduced<Material>(strainInVoigt))};
230 }
231 };
232 }
233
242 template <template <typename, int, int> class RT>
243 requires(LinearElastic::template supportsResultType<RT>())
245 Dune::PriorityTag<1>) const {
246 if constexpr (FE::isMixed())
247 return RTWrapperType<RT>{};
248 if constexpr (isSameResultType<RT, ResultTypes::linearStress> or
249 isSameResultType<RT, ResultTypes::linearStressFull>) {
250 const auto rFunction = resultFunction<RT>();
251 const auto eps = strainFunction(req);
252 auto epsVoigt = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
253 return rFunction(epsVoigt);
254 }
255 }
256
257private:
258 //> CRTP
259 const auto& underlying() const { return static_cast<const FE&>(*this); }
260 auto& underlying() { return static_cast<FE&>(*this); }
261
262 std::shared_ptr<const Geometry> geo_;
263 Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>> localBasis_;
264 Material mat_;
265 size_t numberOfNodes_{0};
266 int order_{};
267
268public:
280 template <typename ST>
281 auto geometricStiffnessMatrixFunction(const Requirement& par, typename Traits::template MatrixType<ST>& K,
282 const VectorXOptRef<ST>& dx = std::nullopt) const {
283 return [&](const FE::template KgType<ST>& kgIJ, const int I, const int J, const auto& gp) {
284 const auto geo = underlying().localView().element().geometry();
285 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
286 K.template block<FE::myDim, FE::myDim>(I * FE::myDim, J * FE::myDim) += kgIJ * intElement;
287 };
288 }
299 template <typename ST>
300 auto materialStiffnessMatrixFunction(const Requirement& par, typename Traits::template MatrixType<ST>& K,
301 const VectorXOptRef<ST>& dx = std::nullopt) const {
302 return [&](const StrainType<ST>& strain, const BopType<ST>& bopI, const BopType<ST>& bopJ, const int I, const int J,
303 const auto& gp) {
304 const auto C = materialTangent(strain);
305 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
306 K.template block<myDim, myDim>(I * myDim, J * myDim) += (bopI.transpose() * C * bopJ) * intElement;
307 };
308 }
309
321 template <typename ST>
322 auto internalForcesFunction(const Requirement& par, typename Traits::template VectorType<ST>& force,
323 const VectorXOptRef<ST>& dx = std::nullopt) const {
324 return [&](const StrainType<ST>& stresses, const BopType<ST>& bopI, const int I, const auto& gp) {
325 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
326 force.template segment<myDim>(myDim * I) += bopI.transpose() * stresses * intElement;
327 };
328 }
329
338 template <typename ScalarType>
339 auto energyFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
340 return [&]() -> ScalarType {
341 using namespace Dune::DerivativeDirections;
342 using namespace Dune;
343 ScalarType energy = 0.0;
344 const auto eps = strainFunction(par, dx);
345 for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
346 const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
347 energy += internalEnergy(epsVoigt) * geo_->integrationElement(gp.position()) * gp.weight();
348 }
349 return energy;
350 };
351 }
352
353protected:
354 template <typename ScalarType>
355 void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
356 typename Traits::template MatrixType<> K,
357 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
358 if constexpr (FE::isMixed()) {
359 return;
360 }
361 using namespace Dune::DerivativeDirections;
362 using namespace Dune;
363 const auto eps = strainFunction(par, dx);
364 const auto kFunction = materialStiffnessMatrixFunction<ScalarType>(par, K, dx);
365
366 for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
367 const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
368 for (size_t i = 0; i < numberOfNodes_; ++i) {
369 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
370 for (size_t j = 0; j < numberOfNodes_; ++j) {
371 const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
372 kFunction(epsVoigt, bopI, bopJ, i, j, gp);
373 }
374 }
375 }
376 }
377
378 template <typename ScalarType>
380 const VectorXOptRef<ScalarType>& dx = std::nullopt) const -> ScalarType {
381 if constexpr (FE::isMixed())
382 return ScalarType{0.0};
383 return energyFunction(par, dx)();
384 }
385
386 template <typename ScalarType>
388 typename Traits::template VectorType<ScalarType> force,
389 const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
390 if constexpr (FE::isMixed())
391 return;
392 using namespace Dune::DerivativeDirections;
393 using namespace Dune;
394 const auto eps = strainFunction(par, dx);
395 const auto fIntFunction = internalForcesFunction<ScalarType>(par, force, dx);
396
397 // Internal forces
398 for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
399 const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
400 const auto stresses = stress(epsVoigt);
401 for (size_t i = 0; i < numberOfNodes_; ++i) {
402 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
403 fIntFunction(stresses, bopI, i, gp);
404 }
405 }
406 }
407};
408
415template <Concepts::GeometricallyLinearMaterial MAT>
416auto linearElastic(const MAT& mat) {
417 LinearElasticPre<MAT> pre(mat);
418
419 return pre;
420}
421} // namespace Ikarus
422
423#else
424 #error LinearElastic depends on dune-localfefunctions, which is not included
425#endif
Definition of the LinearElastic class for finite element mechanics computations.
Header file for various EAS functions.
Definition of several material related enums.
Definitions of ResultTypes used for finite element results.
Material property functions and conversion utilities.
Definition: assemblermanipulatorbuildingblocks.hh:22
MatrixAffordance
A strongly typed enum class representing the matrix affordance.
Definition: ferequirements.hh:65
VectorAffordance
A strongly typed enum class representing the vector affordance.
Definition: ferequirements.hh:50
ScalarAffordance
A strongly typed enum class representing the scalar affordance.
Definition: ferequirements.hh:39
auto linearElastic(const MAT &mat)
A helper function to create a linear elastic pre finite element.
Definition: linearelastic.hh:416
Definition: utils/dirichletvalues.hh:32
FE class is a base class for all finite elements.
Definition: febase.hh:79
static constexpr int myDim
Definition: febase.hh:95
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
Base class for element definitions that provides common functionality for ResultTypes.
Definition: feresulttypes.hh:277
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
BH BasisHandler
Type of the basis of the finite element.
Definition: fetraits.hh:27
typename Basis::GridView GridView
Type of the grid view.
Definition: fetraits.hh:45
typename BasisHandler::FlatBasis FlatBasis
Type of the flat basis.
Definition: fetraits.hh:33
typename LocalView::Element Element
Type of the grid element.
Definition: fetraits.hh:48
static constexpr int mydim
Dimension of the geometry.
Definition: fetraits.hh:63
LinearElastic class represents a linear elastic finite element.
Definition: linearelastic.hh:58
void bindImpl()
A helper function to bind the local view to the element.
Definition: linearelastic.hh:104
auto internalForcesFunction(const Requirement &par, typename Traits::template VectorType< ST > &force, const VectorXOptRef< ST > &dx=std::nullopt) const
Get a lambda function that evaluates the internal force vector for a given strain,...
Definition: linearelastic.hh:322
auto energyFunction(const Requirement &par, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Get a lambda function that evaluates the internal energy at a given integration point and its index.
Definition: linearelastic.hh:339
static constexpr int myDim
Definition: linearelastic.hh:76
Eigen::Vector< ST, strainDim > StrainType
Definition: linearelastic.hh:85
Eigen::Matrix< ST, strainDim, myDim > BopType
Definition: linearelastic.hh:88
void calculateMatrixImpl(const Requirement &par, const MatrixAffordance &affordance, typename Traits::template MatrixType<> K, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: linearelastic.hh:355
Eigen::Matrix< ST, myDim, myDim > KgType
Definition: linearelastic.hh:91
typename Traits::Element Element
Definition: linearelastic.hh:67
auto resultFunction() const
Get a lambda function that evaluates the requested result type for a given strain (in Voigt notation)...
Definition: linearelastic.hh:219
static constexpr auto strainType
Definition: linearelastic.hh:78
PRE Pre
Definition: linearelastic.hh:69
static constexpr auto stressType
Definition: linearelastic.hh:79
auto calculateScalarImpl(const Requirement &par, ScalarAffordance affordance, const VectorXOptRef< ScalarType > &dx=std::nullopt) const -> ScalarType
Definition: linearelastic.hh:379
size_t numberOfNodes() const
Definition: linearelastic.hh:199
int order() const
Definition: linearelastic.hh:200
auto materialStiffnessMatrixFunction(const Requirement &par, typename Traits::template MatrixType< ST > &K, const VectorXOptRef< ST > &dx=std::nullopt) const
Get a lambda function that evaluates the elastic stiffness matrix for a given strain,...
Definition: linearelastic.hh:300
auto materialTangent(const Eigen::Vector< ScalarType, strainDim > &strain=Eigen::Vector< ScalarType, strainDim >::Zero()) const
Get the linear elastic material tangent for the given strain for the given Requirement.
Definition: linearelastic.hh:192
auto calculateAtImpl(const Requirement &req, const Dune::FieldVector< double, Traits::mydim > &local, Dune::PriorityTag< 1 >) const
Calculates a requested result at a specific local position.
Definition: linearelastic.hh:244
auto geometricStiffnessMatrixFunction(const Requirement &par, typename Traits::template MatrixType< ST > &K, const VectorXOptRef< ST > &dx=std::nullopt) const
Get a lambda function that evaluates the geometric part of the stiffness matrix (Kg) for a given inte...
Definition: linearelastic.hh:281
typename Traits::GridView GridView
Definition: linearelastic.hh:66
typename Traits::BasisHandler BasisHandler
Definition: linearelastic.hh:61
typename Traits::LocalView LocalView
Definition: linearelastic.hh:64
std::optional< std::reference_wrapper< const Eigen::VectorX< ST > > > VectorXOptRef
Definition: linearelastic.hh:74
typename Traits::FlatBasis FlatBasis
Definition: linearelastic.hh:62
auto displacementFunction(const Requirement &par, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Gets the displacement function for the given Requirement and optional displacement vector.
Definition: linearelastic.hh:133
LinearElastic(const Pre &pre)
Constructor for the LinearElastic class.
Definition: linearelastic.hh:97
PRE::Material Material
Definition: linearelastic.hh:68
static constexpr int strainDim
Definition: linearelastic.hh:77
decltype(std::declval< LocalView >().tree().child(0).finiteElement().localBasis()) LocalBasisType
Definition: linearelastic.hh:71
decltype(auto) material() const
Definition: linearelastic.hh:204
typename Traits::Geometry Geometry
Definition: linearelastic.hh:65
auto internalEnergy(const Eigen::Vector< ScalarType, strainDim > &strain) const
Get the internal energy for the given strain.
Definition: linearelastic.hh:162
const Dune::CachedLocalBasis< std::remove_cvref_t< LocalBasisType > > & localBasis() const
Definition: linearelastic.hh:201
auto strainFunction(const Requirement &par, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Gets the strain function for the given Requirement and optional di splacement vector.
Definition: linearelastic.hh:149
const Geometry & geometry() const
Definition: linearelastic.hh:198
void calculateVectorImpl(const Requirement &par, VectorAffordance affordance, typename Traits::template VectorType< ScalarType > force, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: linearelastic.hh:387
auto stress(const Eigen::Vector< ScalarType, strainDim > &strain) const
Get the stress for the given strain.
Definition: linearelastic.hh:177
A PreFE struct for linear elastic elements.
Definition: linearelastic.hh:40
MAT Material
Definition: linearelastic.hh:41
MAT material
Definition: linearelastic.hh:42
static constexpr bool isMixed()
Definition: mixin.hh:110
Definition: utils/dirichletvalues.hh:34
Concept to check if the underlying scalar type is a dual type.
Definition: utils/concepts.hh:625
Header file for material models in Ikarus finite element mechanics.