12#if HAVE_DUNE_LOCALFEFUNCTIONS
15 #include <type_traits>
17 #include <dune/geometry/quadraturerules.hh>
18 #include <dune/localfefunctions/expressions/linearStrainsExpr.hh>
19 #include <dune/localfefunctions/impl/standardLocalFunction.hh>
31template <
typename PreFE,
typename FE,
typename PRE>
38template <Concepts::GeometricallyLinearMaterial MAT>
44 template <
typename PreFE,
typename FE>
56template <
typename PreFE,
typename FE,
typename PRE>
71 using LocalBasisType =
decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
73 template <
typename ST>
74 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
81 template <
template <
typename,
int,
int>
class RT>
84 template <
typename ST =
double>
87 template <
typename ST =
double>
88 using BopType = Eigen::Matrix<ST, strainDim, myDim>;
90 template <
typename ST =
double>
91 using KgType = Eigen::Matrix<ST, myDim, myDim>;
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));
117 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_),
118 Dune::bindDerivatives(0, 1));
120 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_), Dune::bindDerivatives(0, 1));
132 template <
typename ScalarType =
double>
135 auto disp = Ikarus::FEHelper::localSolutionBlockVector<Traits>(d, underlying().localView(), dx);
136 Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_);
148 template <
class ScalarType =
double>
161 template <
typename ScalarType,
int strainDim>
163 const auto C = materialTangent<ScalarType, strainDim>(strain);
164 return (0.5 * strain.dot(C * strain));
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();
191 template <
typename ScalarType,
int strainDim,
bool voigt = true>
193 const Eigen::Vector<ScalarType, strainDim>& strain = Eigen::Vector<ScalarType, strainDim>::Zero())
const {
195 return material<ScalarType>().template tangentModuli<strainType, voigt>(strain);
200 [[nodiscard]]
int order()
const {
return order_; }
201 const Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>>&
localBasis()
const {
return localBasis_; }
203 template <
typename ScalarType =
double>
206 return mat_.template rebind<ScalarType>();
217 template <
template <
typename,
int,
int>
class RT>
218 requires(LinearElastic::template supportsResultType<RT>())
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();
229 return RTWrapperType<RT>{mat.template stresses<StrainTags::linear>(enlargeIfReduced<Material>(strainInVoigt))};
242 template <
template <
typename,
int,
int>
class RT>
243 requires(LinearElastic::template supportsResultType<RT>())
245 Dune::PriorityTag<1>)
const {
248 if constexpr (isSameResultType<RT, ResultTypes::linearStress> or
249 isSameResultType<RT, ResultTypes::linearStressFull>) {
250 const auto rFunction = resultFunction<RT>();
252 auto epsVoigt = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
253 return rFunction(epsVoigt);
259 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
260 auto& underlying() {
return static_cast<FE&
>(*this); }
262 std::shared_ptr<const Geometry> geo_;
263 Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>> localBasis_;
265 size_t numberOfNodes_{0};
280 template <
typename ST>
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;
299 template <
typename ST>
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;
321 template <
typename ST>
325 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
326 force.template segment<myDim>(
myDim * I) += bopI.transpose() * stresses * intElement;
338 template <
typename ScalarType>
340 return [&]() -> ScalarType {
341 using namespace Dune::DerivativeDirections;
342 using namespace Dune;
343 ScalarType energy = 0.0;
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();
354 template <
typename ScalarType>
356 typename Traits::template MatrixType<> K,
361 using namespace Dune::DerivativeDirections;
362 using namespace Dune;
364 const auto kFunction = materialStiffnessMatrixFunction<ScalarType>(par, K, dx);
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);
378 template <
typename ScalarType>
382 return ScalarType{0.0};
386 template <
typename ScalarType>
388 typename Traits::template VectorType<ScalarType> force,
392 using namespace Dune::DerivativeDirections;
393 using namespace Dune;
395 const auto fIntFunction = internalForcesFunction<ScalarType>(par, force, dx);
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);
415template <Concepts::GeometricallyLinearMaterial MAT>
424 #error LinearElastic depends on dune-localfefunctions, which is not included
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.