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>
45 template <Concepts::Material DEV, Concepts::Material VOL>
46 struct DPConstitutiveDriver
50 auto Uhat(
const auto& Estrain)
const {
return vol_.template storedEnergy<strainType>(Estrain); }
52 auto storedEnergy(
const auto& Estrain,
const auto& p) {
53 auto Wdev = dev_.template storedEnergy<strainType>(Estrain);
54 auto uhat = Uhat(Estrain);
55 return Wdev + p * uhat - (1.0 / (2.0 * kappa_)) * Dune::power(p, 2);
58 auto firstDerivatives(
const auto& Estrain,
const auto& p)
const {
59 const auto Sdev = dev_.template stresses<strainType, true>(Estrain);
60 const auto dUdE = vol_.template stresses<strainType, true>(Estrain);
61 const auto S = (Sdev + p * dUdE).eval();
62 return std::make_tuple(dUdE, S);
65 auto secondDerivative(
const auto& Estrain,
const auto& p)
const {
66 const auto Cdev = dev_.template tangentModuli<strainType, true>(Estrain);
67 const auto d2UdE2 = vol_.template tangentModuli<strainType, true>(Estrain);
68 const auto C = (Cdev + p * d2UdE2).eval();
72 double kappa()
const {
return kappa_; }
74 DEV deviatoricMaterial()
const {
return dev_; }
75 VOL volumetricMaterial()
const {
return vol_; }
77 DPConstitutiveDriver(
const DEV& dev,
const VOL& vol,
double kappa)
87 template <
typename STO>
89 auto reboundDEV = dev_.template rebind<STO>();
90 auto reboundVOL = vol_.template rebind<STO>();
91 return DPConstitutiveDriver<decltype(reboundDEV), decltype(reboundVOL)>(reboundDEV, reboundVOL, kappa_);
101template <
typename PreFE,
typename FE,
typename PRE>
102class DisplacementPressure;
108template <Concepts::Material MAT>
114 template <
typename PreFE,
typename FE>
127template <
typename PreFE,
typename FE,
typename PRE>
142 decltype(std::declval<LocalView>().tree().child(Dune::Indices::_0).child(0).finiteElement().localBasis());
144 decltype(std::declval<LocalView>().tree().child(Dune::Indices::_1).finiteElement().localBasis());
146 template <
typename ST>
147 using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
154 template <
template <
typename,
int,
int>
class RT>
157 template <
typename ST =
double>
160 template <
typename ST =
double>
161 using BopType = Eigen::Matrix<ST, strainDim, myDim>;
163 template <
typename ST =
double>
164 using KgType = Eigen::Matrix<ST, myDim, myDim>;
168 template <
typename MAT>
171 template <
typename MAT>
182 : constitutiveDriver_([](const auto& mat) {
186 static_assert(Pre::Material::isHyperelastic,
"DisplacementPressure is only implemented for the hyperelastic case.");
194 const auto& localView = underlying().localView();
195 const auto& element = localView.element();
196 auto& firstChild = localView.tree().child(Dune::Indices::_0);
197 const auto& firstFE = firstChild.child(0).finiteElement();
199 auto& secondChild = localView.tree().child(Dune::Indices::_1);
200 const auto& secondFE = secondChild.finiteElement();
202 geo_ = std::make_shared<const Geometry>(element.geometry());
203 numberOfNodes_ = firstFE.size();
204 order_ = 2 * (firstFE.localBasis().order());
205 localBasisU_ = Dune::CachedLocalBasis(firstFE.localBasis());
206 localBasisP_ = Dune::CachedLocalBasis(secondFE.localBasis());
208 localBasisU_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_), Dune::bindDerivatives(0, 1));
209 localBasisP_.bind(Dune::QuadratureRules<double, myDim>::rule(element.type(), order_), Dune::bindDerivatives(0, 1));
221 template <
typename ScalarType =
double>
225 Ikarus::FEHelper::localSolutionBlockVectorComposite<Traits>(d, underlying().localView(), Dune::Indices::_0, dx);
226 Dune::StandardLocalFunction uFunction(localBasisU_, disp, geo_);
238 template <
typename ScalarType =
double>
242 Ikarus::FEHelper::localSolutionBlockVectorScalar<Traits>(d, underlying().localView(), Dune::Indices::_1, dx);
244 Dune::StandardLocalFunction pFunction(localBasisP_, pressure, geo_);
256 template <
typename ScalarType =
double>
263 [[nodiscard]]
int order()
const {
return order_; }
264 const Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisTypeU>>&
localBasis()
const {
return localBasisU_; }
265 const Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisTypeP>>&
localBasisP()
const {
return localBasisP_; }
267 template <
typename ScalarType =
double>
270 return constitutiveDriver_.template rebind<ScalarType>();
272 return constitutiveDriver_;
282 template <
template <
typename,
int,
int>
class RT>
283 requires(supportsResultType<RT>())
285 return [&](
const Eigen::Vector<double, strainDim>& strainInVoigt,
double p) {
286 if constexpr (isSameResultType<RT, ResultTypes::PK2Stress> or isSameResultType<RT, ResultTypes::PK2StressFull>) {
287 auto [dev, vol] = [&]() {
288 if constexpr (isSameResultType<RT, ResultTypes::PK2StressFull> and PRE::Material::isReduced)
289 return std::make_pair(constitutiveDriver_.deviatoricMaterial().underlying(),
290 constitutiveDriver_.volumetricMaterial().underlying());
292 return std::make_pair(constitutiveDriver_.deviatoricMaterial(), constitutiveDriver_.volumetricMaterial());
296 dev.template stresses<strainType>(enlargeIfReduced<typename PRE::Material>(strainInVoigt)) +
297 p * vol.template stresses<strainType>(enlargeIfReduced<typename PRE::Material>(strainInVoigt))};
311 template <
template <
typename,
int,
int>
class RT>
312 requires(supportsResultType<RT>())
314 Dune::PriorityTag<1>)
const {
315 using namespace Dune::DerivativeDirections;
317 if constexpr (isSameResultType<RT, ResultTypes::PK2Stress> or isSameResultType<RT, ResultTypes::PK2StressFull>) {
321 const auto rFunction = resultFunction<RT>();
322 const auto H = uFunction.evaluateDerivative(local, Dune::wrt(spatialAll), Dune::on(gridElement));
323 const auto p = pFunction.evaluate(local, Dune::on(gridElement)).eval();
324 const auto E = (0.5 * (H.transpose() + H + H.transpose() * H)).eval();
326 return rFunction(
toVoigt(E), p[0]);
332 const auto& underlying()
const {
return static_cast<const FE&
>(*this); }
333 auto& underlying() {
return static_cast<FE&
>(*this); }
334 std::shared_ptr<const Geometry> geo_;
335 Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisTypeU>> localBasisU_;
336 Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisTypeP>> localBasisP_;
338 size_t numberOfNodes_{0};
352 template <
typename ST>
355 return [&](
const FE::template
KgType<ST>& kgIJ,
const int I,
const int J,
const auto& gp) {
356 const auto geo = underlying().localView().element().geometry();
357 const double intElement = geo.integrationElement(gp.position()) * gp.weight();
358 K.template block<FE::myDim, FE::myDim>(I *
FE::myDim, J *
FE::myDim) += kgIJ * intElement;
371 template <
typename ST>
375 [&](
const auto& C,
const BopType<ST>& bopI,
const BopType<ST>& bopJ,
const int I,
const int J,
const auto& gp) {
376 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
377 K.template block<myDim, myDim>(I *
myDim, J *
myDim) += (bopI.transpose() * C * bopJ) * intElement;
392 template <
typename ST>
396 const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
397 force.template segment<myDim>(
myDim * I) += bopI.transpose() * stresses * intElement;
409 template <
typename ST>
412 using namespace Dune::DerivativeDirections;
413 using namespace Dune;
417 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
418 const auto EVoigt = eps.evaluate(gpIndex, on(gridElement));
419 const auto p = pFunction.evaluate(gpIndex, on(gridElement)).eval();
421 auto e = constitutiveDriver<ST>().storedEnergy(EVoigt, p[0]);
422 energy += e * geo_->integrationElement(gp.position()) * gp.weight();
436 template <
typename ScalarType>
438 typename Traits::template MatrixType<> K,
441 static_assert(Dune::AlwaysFalse<ScalarType>::value,
442 "DisplacementPressure element does not support matrix calculations for autodiff scalars");
444 using namespace Dune::DerivativeDirections;
445 using namespace Dune;
448 auto pFE = underlying().localView().tree().child(Dune::Indices::_1);
449 auto plocalBasis = pFE.finiteElement().localBasis();
450 std::vector<Dune::FieldVector<double, 1>> Np;
453 const auto kMFunction = materialStiffnessMatrixFunction<ScalarType>(par, K, dx);
454 const auto kGFunction = geometricStiffnessMatrixFunction<ScalarType>(par, K, dx);
455 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
456 const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
457 const auto p = pFunction.evaluate(gpIndex, on(gridElement)).eval();
458 plocalBasis.evaluateFunction(gp.position(), Np);
459 const auto intElement = geo_->integrationElement(gp.position()) * gp.weight();
461 auto [dUdE, S] = constitutiveDriver<ScalarType>().firstDerivatives(EVoigt, p[0]);
462 auto C = constitutiveDriver<ScalarType>().secondDerivative(EVoigt, p[0]);
465 for (
size_t i = 0; i < numberOfNodes_; ++i) {
466 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
467 for (
size_t j = 0; j < numberOfNodes_; ++j) {
468 const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
469 const auto kgIJ = eps.evaluateDerivative(gpIndex, wrt(coeff(i, j)), along(S), on(gridElement));
470 kMFunction(C, bopI, bopJ, i, j, gp);
471 kGFunction(kgIJ, i, j, gp);
476 for (
size_t i = 0; i < numberOfNodes_; ++i) {
477 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
478 for (
auto j : Dune::range(pFE.size())) {
479 auto jIdx = pFE.localIndex(j);
480 const auto upTerm = (bopI.transpose() * dUdE * Np[j][0]).eval();
481 K.template block<myDim, 1>(i *
myDim, jIdx) += upTerm * intElement;
486 for (
auto i : Dune::range(pFE.size())) {
487 auto iIdx = pFE.localIndex(i);
488 for (
size_t j = 0; j < numberOfNodes_; ++j) {
489 const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
490 const auto upTerm = (bopJ.transpose() * dUdE * Np[i][0]).
transpose().eval();
491 K.template block<1, myDim>(iIdx, j *
myDim) += upTerm * intElement;
496 for (
auto i : Dune::range(pFE.size())) {
497 auto iIdx = pFE.localIndex(i);
498 for (
auto j : Dune::range(pFE.size())) {
499 auto jIdx = pFE.localIndex(j);
500 auto ppTerm = Np[i][0] * (1.0 / constitutiveDriver_.kappa()) * Np[j][0];
501 K(iIdx, jIdx) -= ppTerm * intElement;
507 template <
typename ScalarType>
513 static_assert(Dune::AlwaysFalse<ScalarType>::value,
514 "DisplacementPressure element does not support scalar calculations for autodiff scalars");
518 template <
typename ScalarType>
520 typename Traits::template VectorType<ScalarType> force,
523 static_assert(Dune::AlwaysFalse<ScalarType>::value,
524 "DisplacementPressure element does not support vector calculations for autodiff scalars");
526 using namespace Dune::DerivativeDirections;
527 using namespace Dune;
531 auto pFE = underlying().localView().tree().child(Dune::Indices::_1);
532 auto plocalBasis = pFE.finiteElement().localBasis();
533 std::vector<Dune::FieldVector<double, 1>> Np;
535 const auto fIntFunction = internalForcesFunction<ScalarType>(par, force, dx);
537 for (
const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
538 const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
539 const auto p = pFunction.evaluate(gpIndex, on(gridElement)).eval();
540 plocalBasis.evaluateFunction(gp.position(), Np);
542 auto [dUdE, S] = constitutiveDriver<ScalarType>().firstDerivatives(EVoigt, p[0]);
543 auto Uhat = constitutiveDriver<ScalarType>().Uhat(EVoigt);
546 for (
size_t i = 0; i < numberOfNodes_; ++i) {
547 const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
548 fIntFunction(S, bopI, i, gp);
552 for (
auto i : Dune::range(pFE.size())) {
553 auto iIdx = pFE.localIndex(i);
554 auto pTerm = Np[i][0] * (Uhat - p[0] / constitutiveDriver_.kappa());
555 force(iIdx) += pTerm * geo_->integrationElement(gp.position()) * gp.weight();
567template <Concepts::Material MAT>
576 #error DisplacementPressure depends on dune-localfefunctions, which is not included
Collection of fallback default functions.
Helper for transform between Dune linear algebra types and Eigen.
Helper for the autodiff library.
Header file for various EAS functions.
Definition of several material related enums.
Definition of helper struct and function to decompose a hyperelastic material model.
Header file for types of loads in Ikarus finite element mechanics.
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.
constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept
Converts 2D indices to Voigt notation index.
Definition: tensorutils.hh:182
Definition: assemblermanipulatorbuildingblocks.hh:22
MatrixAffordance
A strongly typed enum class representing the matrix affordance.
Definition: ferequirements.hh:65
auto displacementPressure(const MAT &mat)
A helper function to create a displacement-pressure pre finite element.
Definition: displacementpressure.hh:568
VectorAffordance
A strongly typed enum class representing the vector affordance.
Definition: ferequirements.hh:50
auto transpose(const Eigen::EigenBase< Derived > &A)
ScalarAffordance
A strongly typed enum class representing the scalar affordance.
Definition: ferequirements.hh:39
auto decomposeHyperelasticAndGetMaterialParameters(const MAT &mat)
A helper function to decompose a hyperelastic material model and get the underlying deviatoric materi...
Definition: decomposehyperelastic.hh:47
Definition: utils/dirichletvalues.hh:35
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
std::conditional_t< useFlat, FlatBasis, UntouchedBasis > Basis
Type of the basis version.
Definition: fetraits.hh:39
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
DisplacementPressure class represents a displacement-pressure finite element.
Definition: displacementpressure.hh:129
decltype(std::declval< LocalView >().tree().child(Dune::Indices::_0).child(0).finiteElement().localBasis()) LocalBasisTypeU
Definition: displacementpressure.hh:142
static constexpr auto strainType
Definition: displacementpressure.hh:151
DisplacementPressure(const Pre &pre)
Constructor for the DisplacementPressure class.
Definition: displacementpressure.hh:181
typename Traits::Geometry Geometry
Definition: displacementpressure.hh:136
static constexpr auto stressType
Definition: displacementpressure.hh:152
decltype(auto) constitutiveDriver() const
Definition: displacementpressure.hh:268
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: displacementpressure.hh:313
typename Pre::Material MaterialType
Definition: displacementpressure.hh:166
void calculateVectorImpl(const Requirement &par, VectorAffordance affordance, typename Traits::template VectorType< ScalarType > force, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Definition: displacementpressure.hh:519
int order() const
Definition: displacementpressure.hh:263
std::optional< std::reference_wrapper< const Eigen::VectorX< ST > > > VectorXOptRef
Definition: displacementpressure.hh:147
typename Traits::Basis Basis
Definition: displacementpressure.hh:132
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: displacementpressure.hh:353
typename Traits::Element Element
Definition: displacementpressure.hh:138
Impl::DPConstitutiveDriver< DecomposedDevType< MaterialType >, DecomposedVolType< MaterialType > > DPConstitutiveDriverType
Definition: displacementpressure.hh:175
auto resultFunction() const
Get a lambda function that evaluates the requested result type for a given strain (in Voigt notation)...
Definition: displacementpressure.hh:284
typename Materials::DecomposedMaterialTypes< MaterialType >::DevType DecomposedDevType
Definition: displacementpressure.hh:169
const Geometry & geometry() const
Definition: displacementpressure.hh:261
size_t numberOfNodes() const
Definition: displacementpressure.hh:262
const Dune::CachedLocalBasis< std::remove_cvref_t< LocalBasisTypeU > > & localBasis() const
Definition: displacementpressure.hh:264
typename Traits::GridView GridView
Definition: displacementpressure.hh:137
typename Traits::LocalView LocalView
Definition: displacementpressure.hh:135
void bindImpl()
A helper function to bind the local view to the element.
Definition: displacementpressure.hh:193
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 material part of the stiffness matrix (Ke + Ku) for a given ...
Definition: displacementpressure.hh:372
Eigen::Matrix< ST, myDim, myDim > KgType
Definition: displacementpressure.hh:164
PRE Pre
Definition: displacementpressure.hh:139
typename Materials::DecomposedMaterialTypes< MaterialType >::VolType DecomposedVolType
Definition: displacementpressure.hh:172
static constexpr int strainDim
Definition: displacementpressure.hh:150
typename Traits::FlatBasis FlatBasis
Definition: displacementpressure.hh:133
void calculateMatrixImpl(const Requirement &par, const MatrixAffordance &affordance, typename Traits::template MatrixType<> K, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Calculate the matrix associated with the given Requirement.
Definition: displacementpressure.hh:437
Eigen::Vector< ST, strainDim > StrainType
Definition: displacementpressure.hh:158
auto strainFunction(const Requirement &par, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
The strain function for the given Requirement.
Definition: displacementpressure.hh:257
auto displacementFunction(const Requirement &par, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Get the displacement function for the given Requirement.
Definition: displacementpressure.hh:222
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: displacementpressure.hh:393
decltype(std::declval< LocalView >().tree().child(Dune::Indices::_1).finiteElement().localBasis()) LocalBasisTypeP
Definition: displacementpressure.hh:144
static constexpr int myDim
Definition: displacementpressure.hh:149
auto pressureFunction(const Requirement &par, const VectorXOptRef< ScalarType > &dx=std::nullopt) const
Get the pressure function for the given Requirement.
Definition: displacementpressure.hh:239
Eigen::Matrix< ST, strainDim, myDim > BopType
Definition: displacementpressure.hh:161
auto calculateScalarImpl(const Requirement &par, ScalarAffordance affordance, const VectorXOptRef< ScalarType > &dx=std::nullopt) const -> ScalarType
Definition: displacementpressure.hh:508
const Dune::CachedLocalBasis< std::remove_cvref_t< LocalBasisTypeP > > & localBasisP() const
Definition: displacementpressure.hh:265
auto energyFunction(const Requirement &par, const VectorXOptRef< ST > &dx=std::nullopt) const
Get a lambda function that evaluates the internal energy at a given integration point and its index.
Definition: displacementpressure.hh:410
A PreFE struct for displacement-pressure elements.
Definition: displacementpressure.hh:110
MAT material
Definition: displacementpressure.hh:112
MAT Material
Definition: displacementpressure.hh:111
std::tuple_element_t< 1, TupleType > VolType
Type of the pure volumetric material model.
Definition: decomposehyperelastic.hh:25
std::tuple_element_t< 0, TupleType > DevType
Type of the pure deviatoric material model.
Definition: decomposehyperelastic.hh:24
Definition: utils/dirichletvalues.hh:37
Concept to check if the underlying scalar type is a dual type.
Definition: utils/concepts.hh:625
Implementation of the Hyperelastic material model.