12#if HAVE_DUNE_LOCALFEFUNCTIONS
13 #include <dune/fufem/boundarypatch.hh>
14 #include <dune/geometry/quadraturerules.hh>
15 #include <dune/geometry/type.hh>
16 #include <dune/localfefunctions/cachedlocalBasis/cachedlocalBasis.hh>
17 #include <dune/localfefunctions/expressions/greenLagrangeStrains.hh>
18 #include <dune/localfefunctions/impl/standardLocalFunction.hh>
19 #include <dune/localfefunctions/manifolds/realTuple.hh>
33template <
typename PreFE,
typename FE,
typename PRE>
34class NonLinearElastic;
40template <
typename MAT>
46 template <
typename PreFE,
typename FE>
59template <
typename PreFE,
typename FE,
typename PRE>
74 using LocalBasisType =
decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
84 : mat_{pre.material} {}
91 const auto& localView = underlying().localView();
92 const auto& element = localView.element();
93 auto& firstChild = localView.tree().child(0);
94 const auto& fe = firstChild.finiteElement();
95 geo_ = std::make_shared<const Geometry>(element.geometry());
96 numberOfNodes_ = fe.size();
97 order_ = 2 * (fe.localBasis().order());
98 localBasis_ = Dune::CachedLocalBasis(fe.localBasis());
99 if constexpr (
requires { element.impl().getQuadratureRule(order_); })
100 if (element.impl().isTrimmed())
101 localBasis_.bind(element.impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1));
103 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_),
104 Dune::bindDerivatives(0, 1));
106 localBasis_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_), Dune::bindDerivatives(0, 1));
118 template <
typename ScalarType =
double>
121 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
123 auto disp = Ikarus::FEHelper::localSolutionBlockVector<Traits>(d, underlying().localView(), dx);
124 Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_);
136 template <
typename ScalarType =
double>
139 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
152 template <
typename ScalarType,
int strainDim,
bool voigt = true>
154 if constexpr (std::is_same_v<ScalarType, double>)
155 return mat_.template tangentModuli<strainType, voigt>(strain);
157 decltype(
auto) matAD = mat_.template rebind<ScalarType>();
158 return matAD.template tangentModuli<strainType, voigt>(strain);
170 template <
typename ScalarType,
int strainDim>
172 if constexpr (std::is_same_v<ScalarType, double>)
173 return mat_.template storedEnergy<strainType>(strain);
175 decltype(
auto) matAD = mat_.template rebind<ScalarType>();
176 return matAD.template storedEnergy<strainType>(strain);
189 template <
typename ScalarType,
int strainDim,
bool voigt = true>
190 auto getStress(
const Eigen::Vector<ScalarType, strainDim>& strain)
const {
191 if constexpr (std::is_same_v<ScalarType, double>)
192 return mat_.template stresses<strainType, voigt>(strain);
194 decltype(
auto) matAD = mat_.template rebind<ScalarType>();
195 return matAD.template stresses<strainType, voigt>(strain);
201 [[nodiscard]]
int order()
const {
return order_; }
206 template <
template <
typename,
int,
int>
class RT>
207 static consteval bool canProvideResultType() {
208 return isSameResultType<RT, ResultTypes::PK2Stress>;
221 template <
template <
typename,
int,
int>
class RT>
222 requires(canProvideResultType<RT>())
224 Dune::PriorityTag<1>)
const {
225 using namespace Dune::DerivativeDirections;
228 if constexpr (isSameResultType<RT, ResultTypes::PK2Stress>) {
230 const auto H = uFunction.evaluateDerivative(local, Dune::wrt(spatialAll), Dune::on(gridElement));
231 const auto E = (0.5 * (H.transpose() + H + H.transpose() * H)).eval();
233 return RTWrapper{mat_.template stresses<StrainTags::greenLagrangian>(
toVoigt(E))};
239 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
240 auto& underlying() {
return static_cast<FE&
>(*this); }
241 std::shared_ptr<const Geometry> geo_;
242 Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>> localBasis_;
244 size_t numberOfNodes_{0};
255 template <
typename ScalarType>
258 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
259 using namespace Dune::DerivativeDirections;
260 using namespace Dune;
263 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
264 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
265 const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
266 const auto u = (uFunction.evaluate(gpIndex, on(gridElement))).eval();
270 for (
size_t i = 0; i < numberOfNodes_; ++i) {
271 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
272 for (
size_t j = 0; j < numberOfNodes_; ++j) {
273 const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
274 const auto kgIJ = eps.evaluateDerivative(gpIndex, wrt(coeff(i, j)), along(stresses), on(gridElement));
275 K.template block<myDim, myDim>(i *
myDim, j *
myDim) += (bopI.transpose() * C * bopJ + kgIJ) * intElement;
281 template <
typename ScalarType>
283 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx =
284 std::nullopt)
const -> ScalarType {
285 using namespace Dune::DerivativeDirections;
286 using namespace Dune;
290 ScalarType energy = 0.0;
292 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
293 const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
295 energy += internalEnergy * geo_->integrationElement(gp.position()) * gp.weight();
301 template <
typename ScalarType>
303 const FERequirementType& par,
typename Traits::template VectorType<ScalarType> force,
304 const std::optional<std::reference_wrapper<
const Eigen::VectorX<ScalarType>>>& dx = std::nullopt)
const {
305 using namespace Dune::DerivativeDirections;
306 using namespace Dune;
310 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
311 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
312 const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
314 for (
size_t i = 0; i < numberOfNodes_; ++i) {
315 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
316 force.template segment<myDim>(
myDim * i) += bopI.transpose() * stresses * intElement;
328template <
typename MAT>
338 #error NonLinearElastic depends on dune-localfefunctions, which is not included
Helper for transform between Dune linear algebra types and Eigen.
Collection of fallback default functions.
Helper for the autodiff library.
Material property functions and conversion utilities.
Contains the FE class, which is used as a base class for all finite elements. It provides information...
Definition of the LinearElastic class for finite element mechanics computations.
Definition of several material related enums.
Header file for types of loads in Ikarus finite element mechanics.
constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept
Converts 2D indices to Voigt notation index.
Definition: tensorutils.hh:166
Definition: simpleassemblers.hh:22
auto nonLinearElastic(const MAT &mat)
A helper function to create a non-linear elastic pre finite element.
Definition: nonlinearelastic.hh:329
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
typename Basis::GridView GridView
Type of the grid view.
Definition: fetraits.hh:49
std::conditional_t< useFlat, FlatBasis, UntouchedBasis > Basis
Type of the basis version.
Definition: fetraits.hh:43
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
NonLinearElastic class represents a non-linear elastic finite element.
Definition: nonlinearelastic.hh:61
auto getInternalEnergy(const Eigen::Vector< ScalarType, strainDim > &strain) const
Get the internal energy for the given strain.
Definition: nonlinearelastic.hh:171
typename Traits::Basis Basis
Definition: nonlinearelastic.hh:64
typename Traits::Geometry Geometry
Definition: nonlinearelastic.hh:68
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: nonlinearelastic.hh:302
typename Traits::FlatBasis FlatBasis
Definition: nonlinearelastic.hh:65
PRE::Material Material
Definition: nonlinearelastic.hh:71
void calculateMatrixImpl(const FERequirementType &par, typename Traits::template MatrixType<> K, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Calculate the matrix associated with the given FERequirementType.
Definition: nonlinearelastic.hh:256
auto calculateScalarImpl(const FERequirementType &par, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const -> ScalarType
Definition: nonlinearelastic.hh:282
typename Traits::GridView GridView
Definition: nonlinearelastic.hh:69
void bindImpl()
A helper function to bind the local view to the element.
Definition: nonlinearelastic.hh:90
const Geometry & geometry() const
Definition: nonlinearelastic.hh:199
auto strainFunction(const FERequirementType &par, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
The strain function for the given FERequirementType.
Definition: nonlinearelastic.hh:137
static constexpr int myDim
Definition: nonlinearelastic.hh:76
NonLinearElastic(const Pre &pre)
Constructor for the NonLinearElastic class.
Definition: nonlinearelastic.hh:83
typename Traits::Element Element
Definition: nonlinearelastic.hh:70
PRE Pre
Definition: nonlinearelastic.hh:72
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: nonlinearelastic.hh:223
auto displacementFunction(const FERequirementType &par, const std::optional< std::reference_wrapper< const Eigen::VectorX< ScalarType > > > &dx=std::nullopt) const
Get the displacement function for the given FERequirementType.
Definition: nonlinearelastic.hh:119
std::tuple< decltype(makeRT< ResultTypes::PK2Stress >())> SupportedResultTypes
Definition: nonlinearelastic.hh:203
auto materialTangent(const Eigen::Vector< ScalarType, strainDim > &strain) const
Get the material tangent for the given strain.
Definition: nonlinearelastic.hh:153
typename Traits::FERequirementType FERequirementType
Definition: nonlinearelastic.hh:66
static constexpr auto strainType
Definition: nonlinearelastic.hh:77
typename Traits::LocalView LocalView
Definition: nonlinearelastic.hh:67
int order() const
Definition: nonlinearelastic.hh:201
decltype(std::declval< LocalView >().tree().child(0).finiteElement().localBasis()) LocalBasisType
Definition: nonlinearelastic.hh:74
auto getStress(const Eigen::Vector< ScalarType, strainDim > &strain) const
Get the stress for the given strain.
Definition: nonlinearelastic.hh:190
size_t numberOfNodes() const
Definition: nonlinearelastic.hh:200
A PreFE struct for non-linear elastic elements.
Definition: nonlinearelastic.hh:42
MAT Material
Definition: nonlinearelastic.hh:43
MAT material
Definition: nonlinearelastic.hh:44
Definition: utils/dirichletvalues.hh:30