diff --git a/CHANGELOG.md b/CHANGELOG.md index 7873e3813..8f5625ff6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ via `pip install pyikarus` ([#152](https://github.com/ikarus-project/ikarus/pull - Add Clang 16 support ([#186](https://github.com/ikarus-project/ikarus/pull/176)) - Added a default adaptive step sizing possibility and refactored loggers ([#193](https://github.com/ikarus-project/ikarus/pull/193)) - Added a wrapper to fix Dirichlet BCs for Lagrange Nodes ([#222](https://github.com/ikarus-project/ikarus/pull/222)) +- Refactored `getDisplacementFunction` in finite elements ([#223](https://github.com/ikarus-project/ikarus/pull/223)) - Improve material library and Python bindings ([#186](https://github.com/ikarus-project/ikarus/pull/176)), default e.g. `StVenantKirchhoff` is diff --git a/docs/website/01_framework/finiteElements.md b/docs/website/01_framework/finiteElements.md index 1d7c505ef..0fb5949e1 100644 --- a/docs/website/01_framework/finiteElements.md +++ b/docs/website/01_framework/finiteElements.md @@ -110,7 +110,7 @@ calculateScalarImpl(const FERequirementType& par, const Eigen::VectorX& dx, Eigen::VectorX& force) ``` -Inside `getStrainFunction(const FERequirementType& par, const Eigen::VectorX& dx)` is used to get the desired strain measure. +Inside `strainFunction(const FERequirementType& par, const Eigen::VectorX& dx)` is used to get the desired strain measure. It can be used to toggle between the geometrically linear and non-linear cases. `LinearStrains` are used for the geometrically linear case, while `GreenLagrangianStrains` are used for the non-linear case. These strain measures are defined as expressions in `dune-localfefunctions`. Refer to [Expressions](localFunctions.md#expressions) for more details. diff --git a/ikarus/finiteelements/CMakeLists.txt b/ikarus/finiteelements/CMakeLists.txt index b0e7daa7c..d9b7fb7a8 100644 --- a/ikarus/finiteelements/CMakeLists.txt +++ b/ikarus/finiteelements/CMakeLists.txt @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-3.0-or-later # install headers -install(FILES ferequirements.hh fetraits.hh physicshelper.hh +install(FILES ferequirements.hh fetraits.hh fehelper.hh physicshelper.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements ) diff --git a/ikarus/finiteelements/fehelper.hh b/ikarus/finiteelements/fehelper.hh new file mode 100644 index 000000000..f035cd466 --- /dev/null +++ b/ikarus/finiteelements/fehelper.hh @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once + +#include + +#include +#include + +namespace Ikarus::FEHelper { + /** + * @brief Gets the local solution Dune block vector + * + * @tparam FERequirementType Type of the FE requirements. + * @tparam LocalView Type of the local view. + * @tparam ScalarType Scalar type for the local solution vector. + * + * @param x The global solution vector. + * @param localView Local view of the element. + * @param dx Optional global solution vector. + * + * @return A Dune block vector representing the solution quantities at each node. + * */ + template + auto localSolutionBlockVector(const typename FERequirementType::SolutionVectorTypeRaw& x, const LocalView& localView, + const std::optional>& dx = std::nullopt) { + using Traits = TraitsFromLocalView; + constexpr int worldDim = Traits::worlddim; + const auto& fe = localView.tree().child(0).finiteElement(); + Dune::BlockVector> localX(fe.size()); + if (dx) + for (auto i = 0U; i < localX.size(); ++i) + for (auto j = 0U; j < worldDim; ++j) + localX[i][j] = dx.value()[i * worldDim + j] + x[localView.index(localView.tree().child(j).localIndex(i))[0]]; + else + for (auto i = 0U; i < localX.size(); ++i) + for (auto j = 0U; j < worldDim; ++j) + localX[i][j] = x[localView.index(localView.tree().child(j).localIndex(i))[0]]; + return localX; + } +} // namespace Ikarus::FEHelper diff --git a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh index 6a29c4767..c9fe9ee4f 100644 --- a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh +++ b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh @@ -10,8 +10,10 @@ #pragma once #if HAVE_DUNE_LOCALFEFUNCTIONS # include +# include # include +# include # include # include @@ -467,14 +469,11 @@ namespace Ikarus { if (isDisplacementBased()) return; - const auto& d = req.getGlobalSolution(Ikarus::FESolutions::displacement); - const auto C = DisplacementBasedElement::getMaterialTangentFunction(req.getFERequirements()); - const auto& numNodes = DisplacementBasedElement::numberOfNodes; - - Eigen::VectorXd disp(localView().size()); - for (auto i = 0U; i < numNodes; ++i) - for (auto k2 = 0U; k2 < Traits::mydim; ++k2) - disp[i * Traits::mydim + k2] = d[localView().index(localView().tree().child(k2).localIndex(i))[0]]; + const auto& par = req.getFERequirements(); + const auto C = DisplacementBasedElement::materialTangentFunction(req.getFERequirements()); + const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes; + auto uFunction = DisplacementBasedElement::displacementFunction(par); + const auto disp = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef()); std::visit( [&](const EAST& easFunction) { @@ -504,8 +503,8 @@ namespace Ikarus { * @param numberOfEASParameters_ The number of EAS parameters */ void setEASType(int numberOfEASParameters_) { - const auto& numNodes = DisplacementBasedElement::numberOfNodes; - if (not((numNodes == 4 and Traits::mydim == 2) or (numNodes == 8 and Traits::mydim == 3)) + const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes; + if (not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)) and (not isDisplacementBased())) DUNE_THROW(Dune::NotImplemented, "EAS only supported for Q1 or H1 elements"); @@ -561,25 +560,14 @@ namespace Ikarus { DisplacementBasedElement::calculateVectorImpl(par, force, dx); if (isDisplacementBased()) return; using namespace Dune; - const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); - auto strainFunction = DisplacementBasedElement::getStrainFunction(par, dx); - Eigen::VectorX disp(localView().size()); - const auto& numNodes = DisplacementBasedElement::numberOfNodes; - - // FIXME this should not be needed in the future strainFunction should be able to hand out this vector - if (dx) - for (auto i = 0U; i < numNodes; ++i) - for (auto k2 = 0U; k2 < Traits::mydim; ++k2) - disp[i * Traits::mydim + k2] = dx.value()[i * Traits::mydim + k2] - + d[localView().index(localView().tree().child(k2).localIndex(i))[0]]; - else - for (auto i = 0U; i < numNodes; ++i) - for (auto k2 = 0U; k2 < Traits::mydim; ++k2) - disp[i * Traits::mydim + k2] = d[localView().index(localView().tree().child(k2).localIndex(i))[0]]; + const auto uFunction = DisplacementBasedElement::displacementFunction(par, dx); + auto strainFunction = DisplacementBasedElement::strainFunction(par, dx); + const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes; + const auto disp = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef()); using namespace Dune::DerivativeDirections; - auto C = DisplacementBasedElement::getMaterialTangentFunction(par); + auto C = DisplacementBasedElement::materialTangentFunction(par); const auto geo = localView().element().geometry(); // Internal forces from enhanced strains @@ -597,9 +585,10 @@ namespace Ikarus { const double intElement = geo.integrationElement(gp.position()) * gp.weight(); const auto CEval = C(gpIndex); auto stresses = (CEval * M * alpha).eval(); - for (size_t i = 0; i < numNodes; ++i) { + for (size_t i = 0; i < numberOfNodes; ++i) { const auto bopI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement)); - force.template segment(Traits::mydim * i) += bopI.transpose() * stresses * intElement; + force.template segment(Traits::worlddim * i) + += bopI.transpose() * stresses * intElement; } } } @@ -619,10 +608,10 @@ namespace Ikarus { using namespace Dune; using namespace Dune::DerivativeDirections; - auto strainFunction = DisplacementBasedElement::getStrainFunction(par); - const auto C = DisplacementBasedElement::getMaterialTangentFunction(par); - const auto geo = localView().element().geometry(); - const auto numNodes = DisplacementBasedElement::numberOfNodes; + auto strainFunction = DisplacementBasedElement::strainFunction(par); + const auto C = DisplacementBasedElement::materialTangentFunction(par); + const auto geo = localView().element().geometry(); + const auto numberOfNodes = DisplacementBasedElement::numberOfNodes; DMat.setZero(); LMat.setZero(enhancedStrainSize, localView().size()); for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) { @@ -630,7 +619,7 @@ namespace Ikarus { const auto CEval = C(gpIndex); const double detJ = geo.integrationElement(gp.position()); DMat += M.transpose() * CEval * M * detJ * gp.weight(); - for (size_t i = 0U; i < numNodes; ++i) { + for (size_t i = 0U; i < numberOfNodes; ++i) { const size_t I = Traits::worlddim * i; const auto Bi = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement)); LMat.template block(0, I) diff --git a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh index d14cb040b..56d324529 100644 --- a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh +++ b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -36,18 +37,20 @@ namespace Ikarus { template , bool useEigenRef = false> class KirchhoffLoveShell : public PowerBasisFE { public: - using Basis = Basis_; - using FlatBasis = typename Basis::FlatBasis; - using BasePowerFE = PowerBasisFE; // Handles globalIndices function - using FERequirementType = FERequirements_; - using ResultRequirementsType = ResultRequirements; - using LocalView = typename FlatBasis::LocalView; - using Element = typename LocalView::Element; - using Geometry = typename Element::Geometry; - using GridView = typename FlatBasis::GridView; - using Traits = TraitsFromLocalView; - static constexpr int myDim = Traits::mydim; - static constexpr int worlddim = Traits::worlddim; + using Basis = Basis_; + using FlatBasis = typename Basis::FlatBasis; + using BasePowerFE = PowerBasisFE; // Handles globalIndices function + using FERequirementType = FERequirements_; + using ResultRequirementsType = ResultRequirements; + using LocalView = typename FlatBasis::LocalView; + using Element = typename LocalView::Element; + using Geometry = typename Element::Geometry; + using GridView = typename FlatBasis::GridView; + using Traits = TraitsFromLocalView; + static constexpr int myDim = Traits::mydim; + static constexpr int worldDim = Traits::worlddim; + using LocalBasisType = decltype(std::declval().tree().child(0).finiteElement().localBasis()); + static constexpr int membraneStrainSize = 3; static constexpr int bendingStrainSize = 3; @@ -102,10 +105,10 @@ namespace Ikarus { this->localView().bind(element); auto& first_child = this->localView().tree().child(0); const auto& fe = first_child.finiteElement(); + geo_ = std::make_shared(this->localView().element().geometry()); numberOfNodes = fe.size(); - dispAtNodes.resize(fe.size()); - order = 2 * (fe.localBasis().order()); - localBasis = Dune::CachedLocalBasis(fe.localBasis()); + order = 2 * (fe.localBasis().order()); + localBasis = Dune::CachedLocalBasis(fe.localBasis()); if constexpr (requires { this->localView().element().impl().getQuadratureRule(order); }) if (this->localView().element().impl().isTrimmed()) localBasis.bind(this->localView().element().impl().getQuadratureRule(order), Dune::bindDerivatives(0, 1, 2)); @@ -135,26 +138,12 @@ namespace Ikarus { * @param dx Optional additional displacement vector. * @return The displacement function. */ - template - auto getDisplacementFunction(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const { - const auto& d = par.getGlobalSolution(FESolutions::displacement); - auto geo = std::make_shared::Entity::Geometry>( - this->localView().element().geometry()); - Dune::BlockVector> disp(dispAtNodes.size()); - - if (dx) - for (auto i = 0U; i < disp.size(); ++i) - for (auto k2 = 0U; k2 < worlddim; ++k2) - disp[i][k2] = dx.value()[i * worlddim + k2] - + d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - else - for (auto i = 0U; i < disp.size(); ++i) - for (auto k2 = 0U; k2 < worlddim; ++k2) - disp[i][k2] = d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - - Dune::StandardLocalFunction uFunction(localBasis, disp, geo); - + template + auto displacementFunction(const FERequirementType& par, + const std::optional>& dx = std::nullopt) const { + const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); + auto disp = Ikarus::FEHelper::localSolutionBlockVector(d, this->localView(), dx); + Dune::StandardLocalFunction uFunction(localBasis, disp, geo_); return uFunction; } @@ -205,17 +194,13 @@ namespace Ikarus { DUNE_THROW(Dune::NotImplemented, "No results are implemented"); } - Dune::CachedLocalBasis< - std::remove_cvref_t().tree().child(0).finiteElement().localBasis())>> - localBasis; - std::function(const Dune::FieldVector&, - const double&)> + std::shared_ptr geo_; + Dune::CachedLocalBasis> localBasis; + std::function(const Dune::FieldVector&, const double&)> volumeLoad; - std::function(const Dune::FieldVector&, - const double&)> + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; - mutable Dune::BlockVector> dispAtNodes; DefaultMembraneStrain membraneStrain; double emod_; double nu_; @@ -250,7 +235,7 @@ namespace Ikarus { KinematicVariables kin; using namespace Dune; using namespace Dune::DerivativeDirections; - const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos); + const auto [X, Jd, Hd] = geo_->impl().zeroFirstAndSecondDerivativeOfPosition(gpPos); kin.J = toEigen(Jd); kin.H = toEigen(Hd); const Eigen::Matrix A = kin.J * kin.J.transpose(); @@ -274,7 +259,7 @@ namespace Ikarus { kin.a3 = kin.a3N.normalized(); Eigen::Vector bV = kin.h * kin.a3; bV(2) *= 2; // Voigt notation requires the two here - const auto BV = toVoigt(toEigen(geo.impl().secondFundamentalForm(gpPos))); + const auto BV = toVoigt(toEigen(geo_->impl().secondFundamentalForm(gpPos))); kin.kappaV = BV - bV; return kin; } @@ -285,35 +270,34 @@ namespace Ikarus { K.setZero(); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto uFunction = getDisplacementFunction(par, dx); + const auto uFunction = displacementFunction(par, dx); const auto& lambda = par.getParameter(FEParameter::loadfactor); - const auto geo = this->localView().element().geometry(); for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { - const auto intElement = geo.integrationElement(gp.position()) * gp.weight(); + const auto intElement = geo_->integrationElement(gp.position()) * gp.weight(); const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] - = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + = computeMaterialAndStrains(gp.position(), gpIndex, *geo_, uFunction); const Eigen::Vector membraneForces = thickness_ * C * epsV; const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; const auto& Nd = localBasis.evaluateJacobian(gpIndex); const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); for (size_t i = 0; i < numberOfNodes; ++i) { - Eigen::Matrix bopIMembrane - = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, i); + Eigen::Matrix bopIMembrane + = membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, i); - Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); + Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); for (size_t j = i; j < numberOfNodes; ++j) { - auto KBlock = K.template block(worlddim * i, worlddim * j); - Eigen::Matrix bopJMembrane - = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, j); - Eigen::Matrix bopJBending = bopBending(jE, h, Nd, Ndd, j, a3N, a3); + auto KBlock = K.template block(worldDim * i, worldDim * j); + Eigen::Matrix bopJMembrane + = membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, j); + Eigen::Matrix bopJBending = bopBending(jE, h, Nd, Ndd, j, a3N, a3); KBlock += thickness_ * bopIMembrane.transpose() * C * bopJMembrane * intElement; KBlock += Dune::power(thickness_, 3) / 12.0 * bopIBending.transpose() * C * bopJBending * intElement; - Eigen::Matrix kgMembraneIJ - = membraneStrain.secondDerivative(gp.position(), Nd, geo, uFunction, localBasis, membraneForces, i, j); - Eigen::Matrix kgBendingIJ + Eigen::Matrix kgMembraneIJ = membraneStrain.secondDerivative( + gp.position(), Nd, *geo_, uFunction, localBasis, membraneForces, i, j); + Eigen::Matrix kgBendingIJ = kgBending(jE, h, Nd, Ndd, a3N, a3, moments, i, j); KBlock += kgMembraneIJ * intElement; KBlock += kgBendingIJ * intElement; @@ -329,14 +313,13 @@ namespace Ikarus { force.setZero(); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto uFunction = getDisplacementFunction(par, dx); + const auto uFunction = displacementFunction(par, dx); const auto& lambda = par.getParameter(FEParameter::loadfactor); - const auto geo = this->localView().element().geometry(); // Internal forces for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] - = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + = computeMaterialAndStrains(gp.position(), gpIndex, *geo_, uFunction); const Eigen::Vector membraneForces = thickness_ * C * epsV; const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; @@ -344,24 +327,24 @@ namespace Ikarus { const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); for (size_t i = 0; i < numberOfNodes; ++i) { Eigen::Matrix bopIMembrane - = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, i); + = membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, i); Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); force.template segment<3>(3 * i) - += bopIMembrane.transpose() * membraneForces * geo.integrationElement(gp.position()) * gp.weight(); + += bopIMembrane.transpose() * membraneForces * geo_->integrationElement(gp.position()) * gp.weight(); force.template segment<3>(3 * i) - += bopIBending.transpose() * moments * geo.integrationElement(gp.position()) * gp.weight(); + += bopIBending.transpose() * moments * geo_->integrationElement(gp.position()) * gp.weight(); } } // External forces volume forces over the domain if (volumeLoad) { - const auto u = getDisplacementFunction(par, dx); + const auto u = displacementFunction(par, dx); for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { - Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); + Eigen::Vector fext = volumeLoad(geo_->global(gp.position()), lambda); for (size_t i = 0; i < numberOfNodes; ++i) { const auto udCi = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); - force.template segment(worlddim * i) - -= udCi * fext * geo.integrationElement(gp.position()) * gp.weight(); + force.template segment(worldDim * i) + -= udCi * fext * geo_->integrationElement(gp.position()) * gp.weight(); } } } @@ -385,7 +368,7 @@ namespace Ikarus { // Value of the Neumann data at the current position auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); - force.template segment(worlddim * i) -= udCi * neumannValue * intElement; + force.template segment(worldDim * i) -= udCi * neumannValue * intElement; } } } @@ -396,25 +379,24 @@ namespace Ikarus { = std::nullopt) const -> ScalarType { using namespace Dune::DerivativeDirections; using namespace Dune; - const auto uFunction = getDisplacementFunction(par, dx); - const auto& lambda = par.getParameter(FEParameter::loadfactor); - const auto geo = this->localView().element().geometry(); + const auto uFunction = displacementFunction(par, dx); + const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); ScalarType energy = 0.0; for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto [C, epsV, kappaV, j, J, h, H, a3N, a3] - = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + = computeMaterialAndStrains(gp.position(), gpIndex, *geo_, uFunction); const ScalarType membraneEnergy = 0.5 * thickness_ * epsV.dot(C * epsV); const ScalarType bendingEnergy = 0.5 * Dune::power(thickness_, 3) / 12.0 * kappaV.dot(C * kappaV); - energy += (membraneEnergy + bendingEnergy) * geo.integrationElement(gp.position()) * gp.weight(); + energy += (membraneEnergy + bendingEnergy) * geo_->integrationElement(gp.position()) * gp.weight(); } if (volumeLoad) { for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { - const auto u = uFunction.evaluate(gpIndex); - const Eigen::Vector fExt = volumeLoad(geo.global(gp.position()), lambda); - energy -= u.dot(fExt) * geo.integrationElement(gp.position()) * gp.weight(); + const auto u = uFunction.evaluate(gpIndex); + const Eigen::Vector fExt = volumeLoad(geo_->global(gp.position()), lambda); + energy -= u.dot(fExt) * geo_->integrationElement(gp.position()) * gp.weight(); } } diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh index 567ce62c3..865773c66 100644 --- a/ikarus/finiteelements/mechanics/linearelastic.hh +++ b/ikarus/finiteelements/mechanics/linearelastic.hh @@ -10,6 +10,7 @@ #pragma once #if HAVE_DUNE_LOCALFEFUNCTIONS + # include # include # include @@ -21,6 +22,7 @@ # include # include +# include # include # include # include @@ -41,18 +43,19 @@ namespace Ikarus { template , bool useEigenRef = false> class LinearElastic : public PowerBasisFE { public: - using Basis = Basis_; - using FlatBasis = typename Basis::FlatBasis; - using BaseDisp = PowerBasisFE; // Handles globalIndices function - using FERequirementType = FERequirements_; - using ResultRequirementsType = ResultRequirements; - using LocalView = typename FlatBasis::LocalView; - using Element = typename LocalView::Element; - using GridView = typename FlatBasis::GridView; - - using Traits = TraitsFromLocalView; - - static constexpr int myDim = Traits::mydim; + using Basis = Basis_; + using FlatBasis = typename Basis::FlatBasis; + using BaseDisp = PowerBasisFE; // Handles globalIndices function + using FERequirementType = FERequirements_; + using ResultRequirementsType = ResultRequirements; + using LocalView = typename FlatBasis::LocalView; + using Element = typename LocalView::Element; + using Traits = TraitsFromLocalView; + using Geometry = typename Element::Geometry; + using GridView = typename FlatBasis::GridView; + static constexpr int myDim = Traits::mydim; + static constexpr int worldDim = Traits::worlddim; + using LocalBasisType = decltype(std::declval().tree().child(0).finiteElement().localBasis()); /** * @brief Constructor for the LinearElastic class. @@ -75,10 +78,10 @@ namespace Ikarus { this->localView().bind(element); auto& first_child = this->localView().tree().child(0); const auto& fe = first_child.finiteElement(); + geo_ = std::make_shared(this->localView().element().geometry()); numberOfNodes = fe.size(); - dispAtNodes.resize(numberOfNodes); - order = 2 * (this->localView().tree().child(0).finiteElement().localBasis().order()); - localBasis = Dune::CachedLocalBasis(this->localView().tree().child(0).finiteElement().localBasis()); + order = 2 * (fe.localBasis().order()); + localBasis = Dune::CachedLocalBasis(fe.localBasis()); if constexpr (requires { this->localView().element().impl().getQuadratureRule(order); }) if (this->localView().element().impl().isTrimmed()) localBasis.bind(this->localView().element().impl().getQuadratureRule(order), Dune::bindDerivatives(0, 1)); @@ -104,26 +107,12 @@ namespace Ikarus { * @param dx Optional displacement vector. * @return The displacement function. */ - template - auto getDisplacementFunction(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const { + template + auto displacementFunction(const FERequirementType& par, + const std::optional>& dx = std::nullopt) const { const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); - - Dune::BlockVector> disp(dispAtNodes.size()); - if (dx) { - for (auto i = 0U; i < disp.size(); ++i) - for (auto k2 = 0U; k2 < myDim; ++k2) - disp[i][k2] = dx.value()[i * myDim + k2] - + d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - } else - for (auto i = 0U; i < disp.size(); ++i) - for (auto k2 = 0U; k2 < myDim; ++k2) - disp[i][k2] = d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - - auto geo = std::make_shared::Entity::Geometry>( - this->localView().element().geometry()); - Dune::StandardLocalFunction uFunction(localBasis, disp, geo); - + auto disp = Ikarus::FEHelper::localSolutionBlockVector(d, this->localView(), dx); + Dune::StandardLocalFunction uFunction(localBasis, disp, geo_); return uFunction; } /** @@ -135,9 +124,9 @@ namespace Ikarus { * @return The strain function. */ template - auto getStrainFunction(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const { - return linearStrains(getDisplacementFunction(par, dx)); + auto strainFunction(const FERequirementType& par, + const std::optional>& dx = std::nullopt) const { + return Dune::linearStrains(displacementFunction(par, dx)); } /** @@ -145,7 +134,7 @@ namespace Ikarus { * * @return The material tangent matrix. */ - auto getMaterialTangent() const { + auto materialTangent() const { if constexpr (myDim == 2) return planeStressLinearElasticMaterialTangent(emod_, nu_); else if constexpr (myDim == 3) @@ -158,8 +147,8 @@ namespace Ikarus { * @param par The FERequirementType object. * @return The material tangent function. */ - auto getMaterialTangentFunction([[maybe_unused]] const FERequirementType& par) const { - return [&]([[maybe_unused]] auto gp) { return getMaterialTangent(); }; + auto materialTangentFunction([[maybe_unused]] const FERequirementType& par) const { + return [&]([[maybe_unused]] auto gp) { return materialTangent(); }; } /** @@ -187,14 +176,13 @@ namespace Ikarus { * @param K Matrix to store the calculated stiffness. */ void calculateMatrix(const FERequirementType& par, typename Traits::template MatrixType<> K) const { - const auto eps = getStrainFunction(par); + const auto eps = strainFunction(par); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto C = getMaterialTangent(); - const auto geo = this->localView().element().geometry(); + const auto C = materialTangent(); for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { - const double intElement = geo.integrationElement(gp.position()) * gp.weight(); + const double intElement = geo_->integrationElement(gp.position()) * gp.weight(); for (size_t i = 0; i < numberOfNodes; ++i) { const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement)); for (size_t j = 0; j < numberOfNodes; ++j) { @@ -218,8 +206,8 @@ namespace Ikarus { using namespace Dune::DerivativeDirections; using namespace Dune; - const auto eps = getStrainFunction(req.getFERequirements()); - const auto C = getMaterialTangent(); + const auto eps = strainFunction(req.getFERequirements()); + const auto C = materialTangent(); auto epsVoigt = eps.evaluate(local, on(gridElement)); auto linearStress = (C * epsVoigt).eval(); @@ -230,17 +218,13 @@ namespace Ikarus { DUNE_THROW(Dune::NotImplemented, "The requested result type is NOT implemented."); } - Dune::CachedLocalBasis< - std::remove_cvref_t().tree().child(0).finiteElement().localBasis())>> - localBasis; - std::function(const Dune::FieldVector&, - const double&)> + std::shared_ptr geo_; + Dune::CachedLocalBasis> localBasis; + std::function(const Dune::FieldVector&, const double&)> volumeLoad; - std::function(const Dune::FieldVector&, - const double&)> + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; - mutable Dune::BlockVector> dispAtNodes; double emod_; double nu_; size_t numberOfNodes{0}; @@ -250,28 +234,27 @@ namespace Ikarus { template auto calculateScalarImpl(const FERequirementType& par, const std::optional>& dx = std::nullopt) const -> ScalarType { - const auto u = getDisplacementFunction(par, dx); - const auto eps = getStrainFunction(par, dx); - const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); + const auto uFunction = displacementFunction(par, dx); + const auto eps = strainFunction(par, dx); + const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto C = getMaterialTangent(); + const auto C = materialTangent(); - const auto geo = this->localView().element().geometry(); ScalarType energy = 0.0; for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { const auto EVoigt = eps.evaluate(gpIndex, on(gridElement)); - energy += (0.5 * EVoigt.dot(C * EVoigt)) * geo.integrationElement(gp.position()) * gp.weight(); + energy += (0.5 * EVoigt.dot(C * EVoigt)) * geo_->integrationElement(gp.position()) * gp.weight(); } // External forces volume forces over the domain if (volumeLoad) { for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { - const auto uVal = u.evaluate(gpIndex); - Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); - energy -= uVal.dot(fext) * geo.integrationElement(gp.position()) * gp.weight(); + const auto uVal = uFunction.evaluate(gpIndex); + Eigen::Vector fext = volumeLoad(geo_->global(gp.position()), lambda); + energy -= uVal.dot(fext) * geo_->integrationElement(gp.position()) * gp.weight(); } } @@ -291,7 +274,7 @@ namespace Ikarus { const double integrationElement = intersection.geometry().integrationElement(curQuad.position()); // The value of the local function - const auto uVal = u.evaluate(quadPos); + const auto uVal = uFunction.evaluate(quadPos); // Value of the Neumann data at the current position auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); @@ -305,17 +288,17 @@ namespace Ikarus { template void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, const std::optional>& dx = std::nullopt) const { - const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); - const auto eps = getStrainFunction(par, dx); + const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); + const auto uFunction = displacementFunction(par, dx); + const auto eps = strainFunction(par, dx); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto C = getMaterialTangent(); - const auto geo = this->localView().element().geometry(); + const auto C = materialTangent(); // Internal forces for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { - const double intElement = geo.integrationElement(gp.position()) * gp.weight(); + const double intElement = geo_->integrationElement(gp.position()) * gp.weight(); auto stresses = (C * eps.evaluate(gpIndex, on(gridElement))).eval(); for (size_t i = 0; i < numberOfNodes; ++i) { const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement)); @@ -325,21 +308,18 @@ namespace Ikarus { // External forces volume forces over the domain if (volumeLoad) { - const auto u = getDisplacementFunction(par, dx); - for (const auto& [gpIndex, gp] : u.viewOverIntegrationPoints()) { - Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); + for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { + Eigen::Vector fext = volumeLoad(geo_->global(gp.position()), lambda); for (size_t i = 0; i < numberOfNodes; ++i) { - const auto udCi = u.evaluateDerivative(gpIndex, wrt(coeff(i))); + const auto udCi = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); force.template segment(myDim * i) - -= udCi * fext * geo.integrationElement(gp.position()) * gp.weight(); + -= udCi * fext * geo_->integrationElement(gp.position()) * gp.weight(); } } } - // External forces, boundary forces, i.e. at the Neumann boundary + // External forces, boundary forces, i.e., at the Neumann boundary if (not neumannBoundary) return; - - const auto u = getDisplacementFunction(par, dx); auto element = this->localView().element(); for (auto&& intersection : intersections(neumannBoundary->gridView(), element)) { if (not neumannBoundary->contains(intersection)) continue; @@ -352,9 +332,9 @@ namespace Ikarus { const double integrationElement = intersection.geometry().integrationElement(curQuad.position()); - // The value of the local function wrt the i-th coef + // The value of the local function wrt the i-th coeff for (size_t i = 0; i < numberOfNodes; ++i) { - const auto udCi = u.evaluateDerivative(quadPos, wrt(coeff(i))); + const auto udCi = uFunction.evaluateDerivative(quadPos, wrt(coeff(i))); // Value of the Neumann data at the current position auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh index 49b42abc1..d5a90b5ac 100644 --- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh +++ b/ikarus/finiteelements/mechanics/nonlinearelastic.hh @@ -19,6 +19,7 @@ # include # include +# include # include # include # include @@ -54,7 +55,9 @@ namespace Ikarus { using GridView = typename FlatBasis::GridView; using Traits = TraitsFromLocalView; static constexpr int myDim = Traits::mydim; + static constexpr int worldDim = Traits::worlddim; static constexpr auto strainType = StrainTags::greenLagrangian; + using LocalBasisType = decltype(std::declval().tree().child(0).finiteElement().localBasis()); /** * @brief Constructor for the NonLinearElastic class. @@ -76,10 +79,10 @@ namespace Ikarus { this->localView().bind(element); auto& first_child = this->localView().tree().child(0); const auto& fe = first_child.finiteElement(); + geo_ = std::make_shared(this->localView().element().geometry()); numberOfNodes = fe.size(); - dispAtNodes.resize(numberOfNodes); - order = 2 * (this->localView().tree().child(0).finiteElement().localBasis().order()); - localBasis = Dune::CachedLocalBasis(this->localView().tree().child(0).finiteElement().localBasis()); + order = 2 * (fe.localBasis().order()); + localBasis = Dune::CachedLocalBasis(fe.localBasis()); if constexpr (requires { this->localView().element().impl().getQuadratureRule(order); }) if (this->localView().element().impl().isTrimmed()) localBasis.bind(this->localView().element().impl().getQuadratureRule(order), Dune::bindDerivatives(0, 1)); @@ -110,22 +113,8 @@ namespace Ikarus { auto displacementFunction(const FERequirementType& par, const std::optional>& dx = std::nullopt) const { const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); - - Dune::BlockVector> disp(dispAtNodes.size()); - // If optional displacement vector is provided, apply it - if (dx) - for (auto i = 0U; i < disp.size(); ++i) - for (auto k2 = 0U; k2 < myDim; ++k2) - disp[i][k2] = dx.value()[i * myDim + k2] - + d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - else - for (auto i = 0U; i < disp.size(); ++i) - for (auto k2 = 0U; k2 < myDim; ++k2) - disp[i][k2] = d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - - auto geo = std::make_shared::Entity::Geometry>( - this->localView().element().geometry()); - Dune::StandardLocalFunction uFunction(localBasis, disp, geo); + auto disp = Ikarus::FEHelper::localSolutionBlockVector(d, this->localView(), dx); + Dune::StandardLocalFunction uFunction(localBasis, disp, geo_); return uFunction; } @@ -138,9 +127,9 @@ namespace Ikarus { * @return The strain function calculated using greenLagrangeStrains. */ template - auto strainFunction(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const { - return greenLagrangeStrains(displacementFunction(par, dx)); + inline auto strainFunction(const FERequirementType& par, + const std::optional>& dx = std::nullopt) const { + return Dune::greenLagrangeStrains(displacementFunction(par, dx)); } /** @@ -153,7 +142,7 @@ namespace Ikarus { * @return The material tangent calculated using the material's tangentModuli function. */ template - auto getMaterialTangent(const Eigen::Vector& strain) const { + auto materialTangent(const Eigen::Vector& strain) const { if constexpr (std::is_same_v) return mat.template tangentModuli(strain); else { @@ -230,12 +219,10 @@ namespace Ikarus { using namespace Dune::DerivativeDirections; using namespace Dune; const auto eps = strainFunction(par); - const auto geo = this->localView().element().geometry(); - for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { - const double intElement = geo.integrationElement(gp.position()) * gp.weight(); + const double intElement = geo_->integrationElement(gp.position()) * gp.weight(); const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval(); - const auto C = getMaterialTangent(EVoigt); + const auto C = materialTangent(EVoigt); const auto stresses = getStress(EVoigt); for (size_t i = 0; i < numberOfNodes; ++i) { const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement)); @@ -272,18 +259,14 @@ namespace Ikarus { DUNE_THROW(Dune::NotImplemented, "The requested result type is NOT implemented."); } - Dune::CachedLocalBasis< - std::remove_cvref_t().tree().child(0).finiteElement().localBasis())>> - localBasis; - std::function(const Dune::FieldVector&, - const double&)> + std::shared_ptr geo_; + Dune::CachedLocalBasis> localBasis; + std::function(const Dune::FieldVector&, const double&)> volumeLoad; - std::function(const Dune::FieldVector&, - const double&)> + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; Material mat; - mutable Dune::BlockVector> dispAtNodes; size_t numberOfNodes{0}; int order{}; @@ -296,20 +279,19 @@ namespace Ikarus { const auto uFunction = displacementFunction(par, dx); const auto eps = strainFunction(par, dx); const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); - const auto geo = this->localView().element().geometry(); ScalarType energy = 0.0; for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval(); const auto internalEnergy = getInternalEnergy(EVoigt); - energy += internalEnergy * geo.integrationElement(gp.position()) * gp.weight(); + energy += internalEnergy * geo_->integrationElement(gp.position()) * gp.weight(); } if (volumeLoad) { for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto u = uFunction.evaluate(gpIndex); - const Eigen::Vector fExt = volumeLoad(geo.global(gp.position()), lambda); - energy -= u.dot(fExt) * geo.integrationElement(gp.position()) * gp.weight(); + const Eigen::Vector fExt = volumeLoad(geo_->global(gp.position()), lambda); + energy -= u.dot(fExt) * geo_->integrationElement(gp.position()) * gp.weight(); } } @@ -345,13 +327,13 @@ namespace Ikarus { const std::optional>& dx = std::nullopt) const { using namespace Dune::DerivativeDirections; using namespace Dune; - const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); - const auto eps = strainFunction(par, dx); - const auto geo = this->localView().element().geometry(); + const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); + const auto uFunction = displacementFunction(par, dx); + const auto eps = strainFunction(par, dx); // Internal forces for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { - const double intElement = geo.integrationElement(gp.position()) * gp.weight(); + const double intElement = geo_->integrationElement(gp.position()) * gp.weight(); const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval(); const auto stresses = getStress(EVoigt); for (size_t i = 0; i < numberOfNodes; ++i) { @@ -362,12 +344,11 @@ namespace Ikarus { // External forces volume forces over the domain if (volumeLoad) { - const auto u = displacementFunction(par, dx); - for (const auto& [gpIndex, gp] : u.viewOverIntegrationPoints()) { - const double intElement = geo.integrationElement(gp.position()) * gp.weight(); - const Eigen::Vector fExt = volumeLoad(geo.global(gp.position()), lambda); + for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { + const double intElement = geo_->integrationElement(gp.position()) * gp.weight(); + const Eigen::Vector fExt = volumeLoad(geo_->global(gp.position()), lambda); for (size_t i = 0; i < numberOfNodes; ++i) { - const auto udCi = u.evaluateDerivative(gpIndex, wrt(coeff(i))); + const auto udCi = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); force.template segment(myDim * i) -= udCi * fExt * intElement; } } @@ -375,8 +356,6 @@ namespace Ikarus { // External forces, boundary forces, i.e., at the Neumann boundary if (not neumannBoundary and not neumannBoundaryLoad) return; - - const auto u = displacementFunction(par, dx); const auto& element = this->localView().element(); for (auto&& intersection : intersections(neumannBoundary->gridView(), element)) { if (not neumannBoundary->contains(intersection)) continue; @@ -391,7 +370,7 @@ namespace Ikarus { // The value of the local function wrt the i-th coefficient for (size_t i = 0; i < numberOfNodes; ++i) { - const auto udCi = u.evaluateDerivative(quadPos, wrt(coeff(i))); + const auto udCi = uFunction.evaluateDerivative(quadPos, wrt(coeff(i))); // Value of the Neumann data at the current position const auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); diff --git a/ikarus/python/finiteelements/registerelement.hh b/ikarus/python/finiteelements/registerelement.hh index c9e6cf033..4dad09755 100644 --- a/ikarus/python/finiteelements/registerelement.hh +++ b/ikarus/python/finiteelements/registerelement.hh @@ -119,8 +119,8 @@ namespace Ikarus::Python { }, pybind11::arg("FERequirements"), pybind11::arg("elementMatrix").noconvert()); - if constexpr (requires { std::declval().getMaterialTangent(); }) - cls.def("getMaterialTangent", [](FE& self) { return self.getMaterialTangent(); }); + if constexpr (requires { std::declval().materialTangent(); }) + cls.def("materialTangent", [](FE& self) { return self.materialTangent(); }); cls.def( "resultAt", diff --git a/tests/src/testklshell.cpp b/tests/src/testklshell.cpp index f77384eeb..d73ba4f3f 100644 --- a/tests/src/testklshell.cpp +++ b/tests/src/testklshell.cpp @@ -147,9 +147,9 @@ static auto NonLinearKLShellLoadControlTR() { const auto maxDisp = std::ranges::max(d); std::cout << std::setprecision(16) << maxDisp << std::endl; - t.check(Dune::FloatCmp::eq(0.2087574597947082, maxDisp, 1e-6)) - << std::setprecision(16) << "The maximum displacement is " << maxDisp << "but it should be " << 0.2087574597947082 - << ". The difference is " << 0.2087574597947082 - maxDisp; + t.check(Dune::FloatCmp::eq(0.2087577577946809, maxDisp, 1e-6)) + << std::setprecision(16) << "The maximum displacement is " << maxDisp << "but it should be " << 0.2087577577946809 + << ". The difference is " << 0.2087577577946809 - maxDisp; return t; } @@ -208,4 +208,5 @@ int main(int argc, char** argv) { TestSuite t("Kirchhoff-Love"); t.subTest(singleElementTest()); t.subTest(NonLinearKLShellLoadControlTR()); + return t.exit(); }