Skip to content

Commit

Permalink
refactor and fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@tarun-mitruka
tarun-mitruka committed Jan 9, 2024
1 parent 02fc1d4 commit 9a48bac
Show file tree
Hide file tree
Showing 7 changed files with 107 additions and 116 deletions.
2 changes: 1 addition & 1 deletion docs/website/01_framework/finiteElements.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ calculateScalarImpl(const FERequirementType& par, const Eigen::VectorX<ScalarTyp
calculateVectorImpl(const FERequirementType& par, const Eigen::VectorX<ScalarType>& dx, Eigen::VectorX<ScalarType>& force)
```
Inside `getStrainFunction(const FERequirementType& par, const Eigen::VectorX<ScalarType>& dx)` is used to get the desired strain measure.
Inside `strainFunction(const FERequirementType& par, const Eigen::VectorX<ScalarType>& 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.
Expand Down
47 changes: 29 additions & 18 deletions ikarus/finiteelements/fehelper.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,22 +3,29 @@

#pragma once

#include <dune/localfefunctions/cachedlocalBasis/cachedlocalBasis.hh>
#include <dune/localfefunctions/impl/standardLocalFunction.hh>
#include <dune/localfefunctions/linearAlgebraHelper.hh>
#include <dune/localfefunctions/manifolds/realTuple.hh>

#include <ikarus/finiteelements/ferequirements.hh>
#include <ikarus/finiteelements/physicshelper.hh>
#include <ikarus/utils/eigendunetransformations.hh>

namespace Ikarus::FEHelper {
/**
* \brief Views a dune fieldvector as an Eigen::Vector as Map, no copies take place!
* @brief Gets the local displacement Dune block vector
*
* @tparam LocalView Type of the local view.
* @tparam FERequirementType Type of the FE requirements.
* @tparam ScalarType Scalar type for the local displacement vector.
*
* @param localView Local view of the element.
* @param par The FERequirementType object.
* @param dx Optional displacement vector.
*
* @return A Dune block vector representing the displacement quantities at each node.
* */
template <typename LocalView, typename FERequirementType, typename ScalarType>
auto getLocalDisplacementBlockVector(const LocalView& localView, const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
auto localDisplacementBlockVector(const LocalView& localView, const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
using Traits = TraitsFromLocalView<LocalView>;
constexpr int worldDim = Traits::worlddim;
const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement);
Expand All @@ -36,18 +43,22 @@ namespace Ikarus::FEHelper {
return disp;
}

/**
* @brief Gets the local displacement flat Eigen vector
*
* @tparam LocalView Type of the local view.
* @tparam FERequirementType Type of the FE requirements.
* @tparam ScalarType Scalar type for the local displacement vector.
*
* @param localView Local view of the element.
* @param par The FERequirementType object.
* @param dx Optional displacement vector.
*
* @return The local displacement vector viewed as flat Eigen vector.
* */
template <typename LocalView, typename FERequirementType, typename ScalarType = double>
auto getLocalDisplacementVector(const LocalView& localView, const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
return Dune::viewAsFlatEigenVector(getLocalDisplacementBlockVector(localView, par, dx));
}

template <typename DispFunctionType, typename LocalView, typename FERequirementType, typename ScalarType>
void updateDisplacementCoEffs(const DispFunctionType& uFunction, const LocalView& localView,
const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
auto disp = getLocalDisplacementBlockVector(localView, par, dx);
auto coEffs = uFunction.coefficientsRef();
coEffs = disp;
inline auto localDisplacementFlatVector(const LocalView& localView, const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
return Dune::viewAsFlatEigenVector(localDisplacementBlockVector(localView, par, dx));
}
} // namespace Ikarus::FEHelper
14 changes: 7 additions & 7 deletions ikarus/finiteelements/mechanics/enhancedassumedstrains.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,9 +469,9 @@ namespace Ikarus {
if (isDisplacementBased()) return;

const auto& par = req.getFERequirements();
const auto C = DisplacementBasedElement::getMaterialTangentFunction(req.getFERequirements());
const auto C = DisplacementBasedElement::materialTangentFunction(req.getFERequirements());
const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes;
const auto disp = Ikarus::FEHelper::getLocalDisplacementVector(this->localView(), par);
const auto disp = Ikarus::FEHelper::localDisplacementFlatVector(this->localView(), par);

std::visit(
[&]<typename EAST>(const EAST& easFunction) {
Expand Down Expand Up @@ -558,13 +558,13 @@ namespace Ikarus {
DisplacementBasedElement::calculateVectorImpl(par, force, dx);
if (isDisplacementBased()) return;
using namespace Dune;
auto strainFunction = DisplacementBasedElement::getStrainFunction(par, dx);
auto strainFunction = DisplacementBasedElement::strainFunction(par, dx);
const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes;
const auto disp = Ikarus::FEHelper::getLocalDisplacementVector(this->localView(), par, dx);
const auto disp = Ikarus::FEHelper::localDisplacementFlatVector(this->localView(), par, dx);

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
Expand Down Expand Up @@ -605,8 +605,8 @@ namespace Ikarus {
using namespace Dune;
using namespace Dune::DerivativeDirections;

auto strainFunction = DisplacementBasedElement::getStrainFunction(par);
const auto C = DisplacementBasedElement::getMaterialTangentFunction(par);
auto strainFunction = DisplacementBasedElement::strainFunction(par);
const auto C = DisplacementBasedElement::materialTangentFunction(par);
const auto geo = localView().element().geometry();
const auto numberOfNodes = DisplacementBasedElement::numberOfNodes;
DMat.setZero();
Expand Down
26 changes: 10 additions & 16 deletions ikarus/finiteelements/mechanics/kirchhoffloveshell.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,16 +103,14 @@ namespace Ikarus {
neumannBoundary{p_neumannBoundary},
emod_{emod},
nu_{nu},
thickness_{thickness},
uFunction(localBasis, dispAtNodes, geo_) {
thickness_{thickness} {
this->localView().bind(element);
auto& first_child = this->localView().tree().child(0);
const auto& fe = first_child.finiteElement();
geo_ = std::make_shared<const Geometry>(this->localView().element().geometry());
numberOfNodes = fe.size();
dispAtNodes.resize(numberOfNodes);
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));
Expand Down Expand Up @@ -143,10 +141,11 @@ namespace Ikarus {
* @return A pair containing the displacement function and nodal displacements.
*/
template <typename ScalarType = double>
inline void getDisplacementFunction(const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx
= std::nullopt) const {
Ikarus::FEHelper::updateDisplacementCoEffs(uFunction, this->localView(), par, dx);
auto displacementFunction(const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
auto disp = Ikarus::FEHelper::localDisplacementBlockVector(this->localView(), par, dx);
Dune::StandardLocalFunction uFunction(localBasis, disp, geo_);
return uFunction;
}

/**
Expand Down Expand Up @@ -180,16 +179,12 @@ namespace Ikarus {
std::function<Eigen::Vector<double, worldDim>(const Eigen::Vector<double, worldDim>&, const double&)>
neumannBoundaryLoad;
const BoundaryPatch<GridView>* neumannBoundary;
mutable Dune::BlockVector<Dune::RealTuple<double, worldDim>> dispAtNodes;
double emod_;
double nu_;
double thickness_;
size_t numberOfNodes{0};
int order{};
std::shared_ptr<const Geometry> geo_;
Dune::StandardLocalFunction<std::remove_cvref_t<LocalBasisType>,
Dune::BlockVector<Dune::RealTuple<double, worldDim>>, Geometry>
uFunction;

protected:
/**
Expand All @@ -208,11 +203,10 @@ namespace Ikarus {
= std::nullopt) const -> ScalarType {
using namespace Dune::DerivativeDirections;
using namespace Dune;
getDisplacementFunction(par, dx);
const auto uNodes = Ikarus::FEHelper::getLocalDisplacementBlockVector(this->localView(), par, dx);
const auto uFunction = displacementFunction(par, dx);
const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor);
ScalarType energy = 0.0;
const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uNodes);
const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef());

for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
const auto [X, Jd, Hd] = geo_->impl().zeroFirstAndSecondDerivativeOfPosition(gp.position());
Expand Down
66 changes: 29 additions & 37 deletions ikarus/finiteelements/mechanics/linearelastic.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,19 +74,14 @@ namespace Ikarus {
LinearElastic(const Basis& globalBasis, const typename LocalView::Element& element, double emod, double nu,
VolumeLoad p_volumeLoad = {}, const BoundaryPatch<GridView>* p_neumannBoundary = nullptr,
NeumannBoundaryLoad p_neumannBoundaryLoad = {})
: BaseDisp(globalBasis.flat(), element),
neumannBoundary{p_neumannBoundary},
emod_{emod},
nu_{nu},
uFunction(localBasis, dispAtNodes, geo_) {
: BaseDisp(globalBasis.flat(), element), neumannBoundary{p_neumannBoundary}, emod_{emod}, nu_{nu} {
this->localView().bind(element);
auto& first_child = this->localView().tree().child(0);
const auto& fe = first_child.finiteElement();
geo_ = std::make_shared<const Geometry>(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 * (this->localView().tree().child(0).finiteElement().localBasis().order());
localBasis = Dune::CachedLocalBasis(this->localView().tree().child(0).finiteElement().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));
Expand All @@ -102,7 +97,7 @@ namespace Ikarus {
neumannBoundaryLoad = p_neumannBoundaryLoad;

assert(((not p_neumannBoundary and not neumannBoundaryLoad) or (p_neumannBoundary and neumannBoundaryLoad))
& &"If you pass a Neumann boundary you should also pass the function for the Neumann load!");
&& "If you pass a Neumann boundary you should also pass the function for the Neumann load!");
}
/**
* @brief Gets the displacement function for the given FERequirementType and optional displacement vector.
Expand All @@ -113,10 +108,11 @@ namespace Ikarus {
* @return The displacement function.
*/
template <typename ScalarType = double>
inline void getDisplacementFunction(const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx
= std::nullopt) const {
Ikarus::FEHelper::updateDisplacementCoEffs(uFunction, this->localView(), par, dx);
auto displacementFunction(const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
auto disp = Ikarus::FEHelper::localDisplacementBlockVector(this->localView(), par, dx);
Dune::StandardLocalFunction uFunction(localBasis, disp, geo_);
return uFunction;
}
/**
* @brief Gets the strain function for the given FERequirementType and optional displacement vector.
Expand All @@ -127,18 +123,17 @@ namespace Ikarus {
* @return The strain function.
*/
template <class ScalarType = double>
auto getStrainFunction(const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
getDisplacementFunction(par, dx);
return Dune::linearStrains(uFunction);
auto strainFunction(const FERequirementType& par,
const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
return Dune::linearStrains(displacementFunction(par, dx));
}

/**
* @brief Gets the material tangent matrix for the linear elastic material.
*
* @return The material tangent matrix.
*/
auto getMaterialTangent() const {
auto materialTangent() const {
if constexpr (myDim == 2)
return planeStressLinearElasticMaterialTangent(emod_, nu_);
else if constexpr (myDim == 3)
Expand All @@ -151,8 +146,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(); };
}

/**
Expand Down Expand Up @@ -180,11 +175,11 @@ 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 C = materialTangent();
for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
for (size_t i = 0; i < numberOfNodes; ++i) {
Expand All @@ -210,8 +205,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();
Expand All @@ -227,27 +222,23 @@ namespace Ikarus {
std::function<Eigen::Vector<double, worldDim>(const Eigen::Vector<double, worldDim>&, const double&)>
neumannBoundaryLoad;
const BoundaryPatch<GridView>* neumannBoundary;
mutable Dune::BlockVector<Dune::RealTuple<double, worldDim>> dispAtNodes;
double emod_;
double nu_;
size_t numberOfNodes{0};
int order{};
std::shared_ptr<const Geometry> geo_;
Dune::StandardLocalFunction<std::remove_cvref_t<LocalBasisType>,
Dune::BlockVector<Dune::RealTuple<double, worldDim>>, Geometry>
uFunction;

protected:
template <typename ScalarType>
auto calculateScalarImpl(const FERequirementType& par, const std::optional<const Eigen::VectorX<ScalarType>>& dx
= std::nullopt) const -> ScalarType {
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();

ScalarType energy = 0.0;
for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
Expand Down Expand Up @@ -295,12 +286,13 @@ namespace Ikarus {
template <typename ScalarType>
void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType<ScalarType> force,
const std::optional<const Eigen::VectorX<ScalarType>>& 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 C = materialTangent();

// Internal forces
for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
Expand All @@ -314,7 +306,7 @@ namespace Ikarus {

// External forces volume forces over the domain
if (volumeLoad) {
getDisplacementFunction(par, dx);
displacementFunction(par, dx);
for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(toEigen(geo_->global(gp.position())), lambda);
for (size_t i = 0; i < numberOfNodes; ++i) {
Expand All @@ -327,7 +319,7 @@ namespace Ikarus {

// External forces, boundary forces, i.e., at the Neumann boundary
if (not neumannBoundary) return;
getDisplacementFunction(par, dx);
displacementFunction(par, dx);
auto element = this->localView().element();
for (auto&& intersection : intersections(neumannBoundary->gridView(), element)) {
if (not neumannBoundary->contains(intersection)) continue;
Expand Down
Loading

0 comments on commit 9a48bac

Please sign in to comment.