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>
28template <
typename PreFE,
typename FE>
38 template <
typename PreFE,
typename FE>
50template <
typename PreFE,
typename FE>
65 using LocalBasisType =
decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
72 : mat_{pre.material} {}
79 const auto& localView = underlying().localView();
80 const auto& element = localView.element();
81 auto& firstChild = localView.tree().child(0);
82 const auto& fe = firstChild.finiteElement();
83 geo_ = std::make_shared<const Geometry>(element.geometry());
84 numberOfNodes_ = fe.size();
85 order_ = 2 * (fe.localBasis().order());
86 localBasis_ = Dune::CachedLocalBasis(fe.localBasis());
87 if constexpr (
requires { element.impl().getQuadratureRule(order_); })
88 if (element.impl().isTrimmed())
89 localBasis_.bind(element.impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1));
91 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_),
92 Dune::bindDerivatives(0, 1));
94 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_), Dune::bindDerivatives(0, 1));
106 template <
typename ScalarType =
double>
109 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
111 auto disp = Ikarus::FEHelper::localSolutionBlockVector<Traits>(d, underlying().localView(), dx);
112 Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_);
123 template <
class ScalarType =
double>
126 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
136 if constexpr (
myDim == 2)
138 else if constexpr (
myDim == 3)
154 [[nodiscard]]
int order()
const {
return order_; }
158 template <
template <
typename,
int,
int>
class RT>
160 return isSameResultType<RT, ResultTypes::linearStress>;
174 template <
template <
typename,
int,
int>
class RT>
175 requires(canProvideResultType<RT>())
177 Dune::PriorityTag<1>)
const {
179 if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
182 auto epsVoigt = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
184 return RTWrapper{(C * epsVoigt).eval()};
190 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
191 auto& underlying() {
return static_cast<FE&
>(*this); }
193 std::shared_ptr<const Geometry> geo_;
194 Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>> localBasis_;
195 YoungsModulusAndPoissonsRatio mat_;
196 size_t numberOfNodes_{0};
200 template <
typename ScalarType>
203 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
205 using namespace Dune::DerivativeDirections;
206 using namespace Dune;
209 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
210 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
211 for (
size_t i = 0; i < numberOfNodes_; ++i) {
212 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
213 for (
size_t j = 0; j < numberOfNodes_; ++j) {
214 const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
215 K.template block<myDim, myDim>(i *
myDim, j *
myDim) += bopI.transpose() * C * bopJ * intElement;
221 template <
typename ScalarType>
223 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx =
224 std::nullopt)
const -> ScalarType {
228 using namespace Dune::DerivativeDirections;
229 using namespace Dune;
233 ScalarType energy = 0.0;
234 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
235 const auto EVoigt = eps.evaluate(gpIndex, on(gridElement));
236 energy += (0.5 * EVoigt.dot(C * EVoigt)) * geo_->integrationElement(gp.position()) * gp.weight();
241 template <
typename ScalarType>
243 const FERequirementType& par,
typename Traits::template VectorType<ScalarType> force,
244 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
246 using namespace Dune::DerivativeDirections;
247 using namespace Dune;
252 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
253 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
254 auto stresses = (C * eps.evaluate(gpIndex, on(gridElement))).eval();
255 for (
size_t i = 0; i < numberOfNodes_; ++i) {
256 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
257 force.template segment<myDim>(
myDim * i) += bopI.transpose() * stresses * intElement;
276 #error LinearElastic depends on dune-localfefunctions, which is not included
Material property functions and conversion utilities.
Definition of the LinearElastic class for finite element mechanics computations.
Header file for material models in Ikarus finite element mechanics.
Definition: simpleassemblers.hh:22
Eigen::Matrix3d planeStressLinearElasticMaterialTangent(double E, double nu)
Computes the plane stress linear elastic material tangent matrix.
Definition: physicshelper.hh:23
Eigen::Matrix< double, 6, 6 > linearElasticMaterialTangent3D(double E, double nu)
Computes the 3D linear elastic material tangent matrix.
Definition: physicshelper.hh:40
auto linearElastic(const YoungsModulusAndPoissonsRatio &mat)
A helper function to create a linear elastic pre finite element.
Definition: linearelastic.hh:268
Definition: utils/dirichletvalues.hh:28
FE class is a base class for all finite elements.
Definition: febase.hh:81
FETraits< BH, FER, useEigenRef, useFlat > Traits
Definition: febase.hh:40
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:26
FER FERequirementType
Type of the requirements for the finite element.
Definition: fetraits.hh:31
typename Basis::LocalView LocalView
Type of the local view.
Definition: fetraits.hh:46
BH BasisHandler
Type of the basis of the finite element.
Definition: fetraits.hh:28
typename Basis::GridView GridView
Type of the grid view.
Definition: fetraits.hh:49
typename BasisHandler::FlatBasis FlatBasis
Type of the flat basis.
Definition: fetraits.hh:37
static constexpr int mydim
Dimension of the geometry.
Definition: fetraits.hh:67
typename LocalView::Element Element
Type of the grid element.
Definition: fetraits.hh:52
typename Element::Geometry Geometry
Type of the element geometry.
Definition: fetraits.hh:55
LinearElastic class represents a linear elastic finite element.
Definition: linearelastic.hh:52
typename Traits::FlatBasis FlatBasis
Definition: linearelastic.hh:56
typename Traits::Geometry Geometry
Definition: linearelastic.hh:59
typename Traits::Element Element
Definition: linearelastic.hh:61
auto strainFunction(const FERequirementType &par, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Gets the strain function for the given FERequirementType and optional displacement vector.
Definition: linearelastic.hh:124
auto materialTangent() const
Gets the material tangent matrix for the linear elastic material.
Definition: linearelastic.hh:135
int order() const
Definition: linearelastic.hh:154
void calculateVectorImpl(const FERequirementType &par, typename Traits::template VectorType< ScalarType > force, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Definition: linearelastic.hh:242
static consteval bool canProvideResultType()
Definition: linearelastic.hh:159
typename Traits::LocalView LocalView
Definition: linearelastic.hh:58
auto materialTangentFunction(const FERequirementType &par) const
Gets the material tangent function for the given FERequirementType.
Definition: linearelastic.hh:148
typename Traits::GridView GridView
Definition: linearelastic.hh:60
LinearElastic(const Pre &pre)
Constructor for the LinearElastic class.
Definition: linearelastic.hh:71
typename Traits::BasisHandler BasisHandler
Definition: linearelastic.hh:55
const Geometry & geometry() const
Definition: linearelastic.hh:152
decltype(std::declval< LocalView >().tree().child(0).finiteElement().localBasis()) LocalBasisType
Definition: linearelastic.hh:65
auto calculateScalarImpl(const FERequirementType &par, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const -> ScalarType
Definition: linearelastic.hh:222
typename Traits::FERequirementType FERequirementType
Definition: linearelastic.hh:57
std::tuple< decltype(makeRT< ResultTypes::linearStress >())> SupportedResultTypes
Definition: linearelastic.hh:156
void bindImpl()
A helper function to bind the local view to the element.
Definition: linearelastic.hh:78
static constexpr int myDim
Definition: linearelastic.hh:64
size_t numberOfNodes() const
Definition: linearelastic.hh:153
auto calculateAtImpl(const FERequirementType &req, const Dune::FieldVector< double, Traits::mydim > &local, Dune::PriorityTag< 1 >) const
Calculates a requested result at a specific local position.
Definition: linearelastic.hh:176
auto displacementFunction(const FERequirementType &par, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Gets the displacement function for the given FERequirementType and optional displacement vector.
Definition: linearelastic.hh:107
void calculateMatrixImpl(const FERequirementType &par, typename Traits::template MatrixType<> K, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Definition: linearelastic.hh:201
A PreFE struct for linear elastic elements.
Definition: linearelastic.hh:35
YoungsModulusAndPoissonsRatio material
Definition: linearelastic.hh:36
Structure representing Young's modulus and shear modulus.
Definition: physicshelper.hh:53
double emodul
Definition: physicshelper.hh:54
double nu
Definition: physicshelper.hh:55
Definition: utils/dirichletvalues.hh:30