Skip to content

Commit

Permalink
turn-off displacement-based element while using EAS
Browse files Browse the repository at this point in the history
  • Loading branch information
@tarun-mitruka
tarun-mitruka committed Jan 14, 2025
1 parent 4593681 commit 1729aba
Show file tree
Hide file tree
Showing 7 changed files with 156 additions and 121 deletions.
2 changes: 1 addition & 1 deletion ikarus/finiteelements/feresulttypes.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ public:
}

private:
StoredType value_{StoredType::Zero(rowsAtCompileTime, colsAtCompileTime)};
StoredType value_{};
};

namespace Impl {
Expand Down
2 changes: 1 addition & 1 deletion ikarus/finiteelements/mechanics/eas/eas2d.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ struct E0
explicit E0(const GEO& /*geometry*/) {}

// returns an Eigen zero expression for optimization purposes
static auto calcM(const Dune::FieldVector<double, 2>& /*quadPos*/) { return MType::Zero(); }
static auto calcM(const Dune::FieldVector<double, GEO::mydimension>& /*quadPos*/) { return MType::Zero(); }
};

/**
Expand Down
18 changes: 12 additions & 6 deletions ikarus/finiteelements/mechanics/easvariants.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,18 +55,24 @@ namespace Impl {
*/
template <typename F>
void operator()(F&& f) const {
std::visit(
[&]<typename EAST>(const EAST& easFunction) {
if constexpr (EAST::enhancedStrainSize != 0)
f(easFunction);
},
var_);
std::visit([&]<typename EAST>(const EAST& easFunction) { f(easFunction); }, var_);
}

/**
* \brief A helper function to get the number of EAS parameters.
* \return Number of EAS parameters.
*/
auto numberOfEASParameters() const {
return std::visit([]<typename EAST>(const EAST&) -> int { return EAST::enhancedStrainSize; }, var_);
}

/**
* \brief A helper function to identify if the formulation is purely displacement-based, i.e., number of EAS
* parameters is zero.
* \return A bool indicating if the formulation is displacement-based or not.
*/
bool isDisplacmentBased() const { return numberOfEASParameters() == 0; }

void setEASType(int numberOfEASParameters) {
numberOfEASParameters_ = numberOfEASParameters;
if (geometry_)
Expand Down
78 changes: 37 additions & 41 deletions ikarus/finiteelements/mechanics/enhancedassumedstrains.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,22 +110,21 @@ public:
template <template <typename, int, int> class RT>
auto calculateAtImpl(const Requirement& req, const Dune::FieldVector<double, Traits::mydim>& local,
Dune::PriorityTag<2>) const {
if (isDisplacementBased())
return underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());

auto resultVector = underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
auto strainFunction = underlying().strainFunction(req);
const auto ufunc = underlying().displacementFunction(req);
auto disp = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
auto strainFunction = underlying().strainFunction(req);
const auto ufunc = underlying().displacementFunction(req);
const auto rFunction = underlying().template resultFunction<RT>();
auto disp = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
using RTWrapper = ResultWrapper<RT<typename Traits::ctype, Traits::mydim, Traits::worlddim>, ResultShape::Vector>;
RTWrapper resultWrapper;
auto calculateAtContribution = [&]<typename EAST>(const EAST& easFunction) {
typename EAST::DType D;
calculateDAndLMatrix(easFunction, req, D, L_);
this->alpha_ = (-D.inverse() * L_ * disp).eval();
if constexpr (EAST::enhancedStrainSize != 0) { // compile-time check
typename EAST::DType D;
calculateDAndLMatrix(easFunction, req, D, L_);
this->alpha_ = (-D.inverse() * L_ * disp).eval();
}
const auto enhancedStrain = enhancedStrainFunction(easFunction, strainFunction, local);
resultWrapper = resultVector.asVec() + underlying().stress(enhancedStrain).eval();
resultWrapper = rFunction(enhancedStrain);
};
easVariant_(calculateAtContribution);
return resultWrapper;
Expand Down Expand Up @@ -173,20 +172,20 @@ protected:
const std::remove_reference_t<typename Traits::template VectorType<>>& correction) const {
using ScalarType = Traits::ctype;
easApplicabilityCheck();
if (isDisplacementBased())
return;
const auto& Rtilde = calculateRtilde<ScalarType>(par);
const auto localdxBlock = Ikarus::FEHelper::localSolutionBlockVector<Traits, Eigen::VectorXd, double>(
correction, underlying().localView());
const auto localdx = Dune::viewAsFlatEigenVector(localdxBlock);
decltype(auto) LMat = LMatFunc<ScalarType>();

auto correctAlpha = [&]<typename EAST>(const EAST& easFunction) {
constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize> D;
calculateDAndLMatrix(easFunction, par, D, LMat);
const auto updateAlpha = (D.inverse() * (Rtilde + (LMat * localdx))).eval();
this->alpha_ -= updateAlpha;
if constexpr (EAST::enhancedStrainSize != 0) { // compile-time check
const auto& Rtilde = calculateRtilde<ScalarType>(par);
const auto localdxBlock = Ikarus::FEHelper::localSolutionBlockVector<Traits, Eigen::VectorXd, double>(
correction, underlying().localView());
const auto localdx = Dune::viewAsFlatEigenVector(localdxBlock);
decltype(auto) LMat = LMatFunc<ScalarType>();
constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize> D;
calculateDAndLMatrix(easFunction, par, D, LMat);
const auto updateAlpha = (D.inverse() * (Rtilde + (LMat * localdx))).eval();
this->alpha_ -= updateAlpha;
}
};

easVariant_(correctAlpha);
Expand All @@ -206,24 +205,23 @@ protected:
if (affordance != MatrixAffordance::stiffness)
DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));

Check warning on line 206 in ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

View check run for this annotation

Codecov / codecov/patch

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh#L206 

Added line #L206 was not covered by tests
easApplicabilityCheck();
if (isDisplacementBased())
return;

auto strainFunction = underlying().strainFunction(par, dx);
const auto kFunction = underlying().template stiffnessMatrixFunction<ScalarType>(par, K, dx);
decltype(auto) LMat = LMatFunc<ScalarType>();

auto calculateMatrixContribution = [&]<typename EAST>(const EAST& easFunction) {
typename EAST::DType D;
calculateDAndLMatrix(easFunction, par, D, LMat, dx);

for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
const auto enhancedStrain = enhancedStrainFunction(easFunction, strainFunction, gp.position());
kFunction(enhancedStrain, gpIndex, gp);
}

K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
if constexpr (EAST::enhancedStrainSize != 0) { // compile-time check
typename EAST::DType D;
decltype(auto) LMat = LMatFunc<ScalarType>();
calculateDAndLMatrix(easFunction, par, D, LMat, dx);
K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
}
};
easVariant_(calculateMatrixContribution);
}
Expand All @@ -233,7 +231,8 @@ protected:
const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
easApplicabilityCheck();
if (isDisplacementBased())
return 0.0;
return underlying().template energyFunction<ScalarType>(par, dx)();

DUNE_THROW(Dune::NotImplemented,
"EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
}
Expand All @@ -245,21 +244,20 @@ protected:
if (affordance != VectorAffordance::forces)
DUNE_THROW(Dune::NotImplemented, "VectorAffordance not implemented: " + toString(affordance));

Check warning on line 245 in ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

View check run for this annotation

Codecov / codecov/patch

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh#L245 

Added line #L245 was not covered by tests
easApplicabilityCheck();
if (isDisplacementBased())
return;
auto strainFunction = underlying().strainFunction(par, dx);
const auto fIntFunction = underlying().template internalForcesFunction<ScalarType>(par, force, dx);

auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
typename EAST::DType D;
calculateDAndLMatrix(easFunction, par, D, L_);

const auto& Rtilde = calculateRtilde(par, dx);
for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
const auto enhancedStrain = enhancedStrainFunction(easFunction, strainFunction, gp.position());
fIntFunction(enhancedStrain, gpIndex, gp);
}
force -= L_.transpose() * D.inverse() * Rtilde;
if constexpr (EAST::enhancedStrainSize != 0) { // compile-time check
typename EAST::DType D;
calculateDAndLMatrix(easFunction, par, D, L_);
const auto& Rtilde = calculateRtilde(par, dx);
force -= L_.transpose() * D.inverse() * Rtilde;
}
};
easVariant_(calculateForceContribution);
}
Expand All @@ -277,8 +275,6 @@ private:
* \brief Initializes the internal state variable alpha_ based on the number of EAS parameters.
*/
void initializeState() {
if (isDisplacementBased())
return;
alpha_.resize(numberOfEASParameters());
alpha_.setZero();
}
Expand Down
85 changes: 51 additions & 34 deletions ikarus/finiteelements/mechanics/linearelastic.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,28 +187,40 @@ public:
[[nodiscard]] size_t numberOfNodes() const { return numberOfNodes_; }
[[nodiscard]] int order() const { return order_; }

/**
* \brief Get a lambda function that evaluates the requested result type for a given strain (in Voigt notation).
* \tparam RT The type representing the requested result.
* \return A lambda function that evaluates the requested result type for a given strain (in Voigt notation).
*/
template <template <typename, int, int> class RT>
requires(supportsResultType<RT>())
auto resultFunction() const {
return [&]<int strainDim>(const Eigen::Vector<double, strainDim>& strainInVoigt) {
using RTWrapper = ResultWrapper<RT<typename Traits::ctype, myDim, Traits::worlddim>, ResultShape::Vector>;
if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
return RTWrapper{stress(strainInVoigt)};
}
};
}

/**
* \brief Calculates a requested result at a specific local position.
*
* \param req The Requirement object holding the global solution.
* \param local Local position vector.
* \tparam RT The requested result type
* \return calculated result
* \return Calculated result.
*
* \tparam RT The type representing the requested result.
*/
template <template <typename, int, int> class RT>
requires(supportsResultType<RT>())
auto calculateAtImpl(const Requirement& req, const Dune::FieldVector<double, Traits::mydim>& local,
Dune::PriorityTag<1>) const {
using RTWrapper = ResultWrapper<RT<typename Traits::ctype, myDim, Traits::worlddim>, ResultShape::Vector>;
if (usesEASSkill())
return RTWrapper{};
if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
const auto eps = strainFunction(req);
auto epsVoigt = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
return RTWrapper{stress(epsVoigt)};
}
if constexpr (hasEAS)
return;
const auto rFunction = resultFunction<RT>();
const auto eps = strainFunction(req);
auto epsVoigt = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
return rFunction(epsVoigt);
}

private:
Expand All @@ -230,13 +242,6 @@ private:
return mat_;
}

bool usesEASSkill() const {
if constexpr (hasEAS)
return underlying().numberOfEASParameters() != 0;
else
return false;
}

public:
/**
* \brief Get a lambda function that evaluates the stiffness matrix for a given strain, integration point and its
Expand Down Expand Up @@ -294,12 +299,35 @@ public:
};
}

/**
* \brief Get a lambda function that evaluates the internal energy at a given integration point and its index.
*
* \tparam ScalarType The scalar type for the material and strain.
* \param par The Requirement object.
* \param dx Optional displacement vector.
* \return A lambda function that returns the intenral energy at a given integration point and its index.
*/
template <typename ScalarType>
auto energyFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
return [&]() -> ScalarType {
using namespace Dune::DerivativeDirections;
using namespace Dune;
ScalarType energy = 0.0;
const auto eps = strainFunction(par, dx);
for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
energy += internalEnergy(epsVoigt) * geo_->integrationElement(gp.position()) * gp.weight();
}
return energy;
};
}

protected:
template <typename ScalarType>
void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
typename Traits::template MatrixType<> K,
const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
if (usesEASSkill())
if constexpr (hasEAS)
return;
const auto eps = strainFunction(par, dx);
const auto kFunction = stiffnessMatrixFunction<ScalarType>(par, K, dx);
Expand All @@ -315,27 +343,16 @@ protected:
template <typename ScalarType>
auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
const VectorXOptRef<ScalarType>& dx = std::nullopt) const -> ScalarType {
ScalarType energy = 0.0;
if (usesEASSkill())
return energy;
const auto uFunction = displacementFunction(par, dx);
const auto eps = strainFunction(par, dx);
const auto& lambda = par.parameter();
using namespace Dune::DerivativeDirections;
using namespace Dune;

for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
energy += internalEnergy(epsVoigt) * geo_->integrationElement(gp.position()) * gp.weight();
}
return energy;
if constexpr (hasEAS)
return ScalarType{0.0};
return energyFunction(par, dx)();
}

template <typename ScalarType>
void calculateVectorImpl(const Requirement& par, VectorAffordance affordance,
typename Traits::template VectorType<ScalarType> force,
const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
if (usesEASSkill())
if constexpr (hasEAS)
return;
const auto eps = strainFunction(par, dx);
const auto fIntFunction = internalForcesFunction<ScalarType>(par, force, dx);
Expand Down
Loading

0 comments on commit 1729aba

Please sign in to comment.