add nonlinear version of eas

all commits squashed for an easy rebase

Author: tarun-mitruka

CHANGELOG.md

@@ -56,6 +56,18 @@ SPDX-License-Identifier: LGPL-3.0-or-later
 - Add an About Ikarus page in the documentation ([#291](https://github.com/ikarus-project/ikarus/pull/291))
 - Add new class `Vtk::Writer`, which implements some convenience methods over the existing `dune-vtk` module ([#309](https://github.com/ikarus-project/ikarus/pull/309))
 - Add `VanishingStrain` material (useful for example for plane strain case), also refactor the constructor of `LinearElastic` to take any linear material law ([#317](https://github.com/ikarus-project/ikarus/pull/317))
+- Refactor EAS to handle NonLinearElastic ([#325](https://github.com/ikarus-project/ikarus/pull/325))
+    - `updateState()` function is added to `mixin.hh`, which can be used to update the internal state variables of the skills.
+        - `updateStateImpl()` now has to be implemented by every skill.
+        - For EAS, `updateStateImpl()` is used to update the internal variable `alpha` in a nonlinear analysis.
+    - Missing functions like `getStress` and `materialTangentFunction` are added to `LinearElastic` and `NonLinearElastic`, respectively.
+    - A `helperfunctions.hh` file is added that contains the helper functions used by the nonlinear solvers.
+        - This contains the function `updateStates()` that calls the `updateState()` function for every finite element.
+        - `updateStates()` is called by nonlinear solvers during every iteration, whenever the correction vector is updated.
+    - `NonLinearOperator` now has an additional template argument, `Assembler`.
+        - It stores a shared pointer to the underlying assembler.
+        - This helps the nonlinear solvers to have access to the underlying assembler,
+          which in turn provides the finite element container to update the state variables.
 
 ## Release v0.4 (Ganymede)
 

docs/literature.bib

@@ -175,11 +175,29 @@ @book{bonet2008nonlinear
 	doi  = {https://doi.org/10.1017/CBO9780511755446}
 }
 
+@article{bischoffShearDeformableShell1997b,
+  title = {{Shear deformable shell elements for large strains and rotations}},
+  author = {Bischoff, M. and Ramm, E.},
+  year = {1997},
+  journal = {International Journal for Numerical Methods in Engineering},
+  volume = {40},
+  number = {23},
+  pages = {4427--4449},
+  issn = {1097-0207},
+  doi = {10.1002/(SICI)1097-0207(19971215)40:23<4427::AID-NME268>3.0.CO;2-9}
+}
 
-
-
-
-
+@article{klinkelGeometricalNonlinearBrick1997,
+  title = {A Geometrical Non-Linear Brick Element Based on the {{EAS-method}}},
+  author = {Klinkel, S. and Wagner, W.},
+  year = {1997},
+  journal = {International Journal for Numerical Methods in Engineering},
+  volume = {40},
+  number = {24},
+  pages = {4529--4545},
+  issn = {1097-0207},
+  doi = {10.1002/(SICI)1097-0207(19971230)40:24<4529::AID-NME271>3.0.CO;2-I}
+}
 
 @book{trustregion,
 	author = {Conn, Andrew R. and Gould, Nicholas I. M. and Toint, Philippe L.},

ikarus/assembler/assemblermanipulatorfuser.hh

@@ -82,6 +82,7 @@ public:
   using WrappedAssembler::boundToDBCOption;
   using WrappedAssembler::boundToRequirement;
   using WrappedAssembler::createFullVector;
+  using WrappedAssembler::finiteElements;
 
   using WrappedAssembler::affordanceCollection;
   using WrappedAssembler::constraintsBelow;
@@ -173,6 +174,7 @@ public:
   using WrappedAssembler::boundToDBCOption;
   using WrappedAssembler::boundToRequirement;
   using WrappedAssembler::createFullVector;
+  using WrappedAssembler::finiteElements;
 
   using WrappedAssembler::affordanceCollection;
   using WrappedAssembler::constraintsBelow;
@@ -283,6 +285,7 @@ public:
   using WrappedAssembler::boundToDBCOption;
   using WrappedAssembler::boundToRequirement;
   using WrappedAssembler::createFullVector;
+  using WrappedAssembler::finiteElements;
 
   using WrappedAssembler::affordanceCollection;
   using WrappedAssembler::constraintsBelow;

ikarus/finiteelements/ferequirements.hh

@@ -203,7 +203,7 @@ auto scalarAffordance(VectorAffordance affordanceV) {
 namespace AffordanceCollections {
   inline constexpr AffordanceCollection elastoStatics(ScalarAffordance::mechanicalPotentialEnergy,
                                                       VectorAffordance::forces, MatrixAffordance::stiffness);
-}
+} // namespace AffordanceCollections
 
 namespace Impl {
   template <typename T>

ikarus/finiteelements/fetraits.hh

@@ -80,6 +80,11 @@ struct FETraits
   using VectorType =
       std::conditional_t<useEigenRef, Eigen::Ref<Eigen::VectorX<ScalarType>>, Eigen::VectorX<ScalarType>&>;
 
+  /** \brief Type of the vector passed to updateState */
+  template <typename ScalarType = ctype>
+  using VectorTypeConst =
+      std::conditional_t<useEigenRef, const Eigen::Ref<Eigen::VectorX<ScalarType>>, const Eigen::VectorX<ScalarType>&>;
+
   /** \brief Type of the matrix  passed to calculateMatrix */
   template <typename ScalarType = ctype>
   using MatrixType =

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -138,27 +138,73 @@ public:
     if (numberOfEASParameters != 0)
       easApplicabilityCheck();
     easVariant_.setEASType(numberOfEASParameters);
+    initializeState();
   }
 
+  /**
+   * \brief Gets the internal state variable alpha for the EAS element.
+   *
+   * \return Internal state variable (alpha).
+   */
+  const auto& alpha() const { return alpha_; }
+
 protected:
   void bindImpl() {
     assert(underlying().localView().bound());
     easVariant_.bind(underlying().localView().element().geometry());
+    initializeState();
+  }
+
+  /**
+   * \brief Updates the internal state variable alpha_ at the end of an iteration
+   * when NonLinearSolverMessages::SOLUTION_UPDATED is notified by the non-linear solver.
+   * See \cite bischoffShearDeformableShell1997b and \cite klinkelGeometricalNonlinearBrick1997 for implementation
+   * details.
+   *
+   * \param par The Requirement object.
+   * \param correction The correction in displacement (DeltaD) vector passed based on which the internal state variable
+   * alpha is to be updated.
+   */
+  void updateStateImpl(const Requirement& par, typename Traits::template VectorTypeConst<> correction) {
+    using ScalarType = Traits::ctype;
+    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 = [this]() -> decltype(auto) {
+      if constexpr (std::is_same_v<ScalarType, double>)
+        return [this]() -> Eigen::MatrixXd& { return L_; }();
+      else
+        return Eigen::MatrixX<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 decltype(alpha_) updateAlpha = D.inverse() * (Rtilde + (LMat * localdx).eval());
+      this->alpha_ -= updateAlpha;
+    };
+
+    easVariant_(correctAlpha);
   }
 
-public:
-protected:
   inline void easApplicabilityCheck() const {
     const auto& numberOfNodes = underlying().numberOfNodes();
     assert(not(not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)) and
                (not isDisplacementBased())) &&
-           "EAS only supported for Q1 or H1 elements");
+           "EAS is only supported for Q1 or H1 elements");
   }
 
   template <typename ScalarType>
   void calculateMatrixImpl(
       const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K,
       const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+    if (affordance != MatrixAffordance::stiffness)
+      DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
     easApplicabilityCheck();
@@ -172,10 +218,30 @@ protected:
         return [this]() -> Eigen::MatrixXd& { return L_; }();
     }();
 
+    auto strainFunction      = underlying().strainFunction(par, dx);
+    const auto C             = underlying().materialTangentFunction(par);
+    const auto geo           = underlying().localView().element().geometry();
+    const auto numberOfNodes = underlying().numberOfNodes();
+
     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 M            = easFunction.calcM(gp.position());
+        const double intElement = geo.integrationElement(gp.position()) * gp.weight();
+        const auto EVoigt       = strainFunction.evaluate(gpIndex, on(gridElement)).eval();
+        const auto CEval        = C(gpIndex);
+        auto stresses           = (CEval * M * alpha_).eval();
+        for (size_t i = 0; i < numberOfNodes; ++i) {
+          for (size_t j = 0; j < numberOfNodes; ++j) {
+            const auto kgIJ =
+                strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i, j)), along(stresses), on(gridElement));
+            K.template block<Traits::mydim, Traits::mydim>(i * Traits::mydim, j * Traits::mydim) += kgIJ * intElement;
+          }
+        }
+      }
+
       K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
       K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
     };
@@ -197,50 +263,59 @@ protected:
   void calculateVectorImpl(
       const Requirement& par, VectorAffordance affordance, typename Traits::template VectorType<ScalarType> force,
       const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+    if (affordance != VectorAffordance::forces)
+      DUNE_THROW(Dune::NotImplemented, "VectorAffordance not implemented: " + toString(affordance));
     easApplicabilityCheck();
+    if (isDisplacementBased())
+      return;
     using namespace Dune;
     using namespace Dune::DerivativeDirections;
     const auto uFunction      = underlying().displacementFunction(par, dx);
     auto strainFunction       = underlying().strainFunction(par, dx);
     const auto& numberOfNodes = underlying().numberOfNodes();
+    const auto C              = underlying().materialTangentFunction(par);
 
     auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
       typename EAST::DType D;
       calculateDAndLMatrix(easFunction, par, D, L_);
 
-      decltype(auto) LMat = [this]() -> decltype(auto) {
-        if constexpr (Concepts::AutodiffScalar<ScalarType>)
-          return Eigen::MatrixX<ScalarType>{};
-        else
-          return [this]() -> Eigen::MatrixXd& { return L_; }();
-      }();
-      const auto disp  = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef());
-      const auto alpha = (-D.inverse() * L_ * disp).eval();
-      const auto geo   = underlying().localView().element().geometry();
-      auto C           = underlying().materialTangentFunction(par);
-
+      const auto disp    = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef());
+      const auto geo     = underlying().localView().element().geometry();
+      const auto& Rtilde = calculateRtilde(par, dx);
       for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
         const auto M            = easFunction.calcM(gp.position());
         const double intElement = geo.integrationElement(gp.position()) * gp.weight();
         const auto CEval        = C(gpIndex);
-        auto stresses           = (CEval * M * alpha).eval();
+        auto stresses           = (CEval * M * alpha_).eval();
         for (size_t i = 0; i < numberOfNodes; ++i) {
           const auto bopI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
           force.template segment<Traits::worlddim>(Traits::worlddim * i) += bopI.transpose() * stresses * intElement;
         }
       }
+      force -= L_.transpose() * D.inverse() * Rtilde;
     };
     easVariant_(calculateForceContribution);
   }
 
 private:
   EAS::Impl::EASVariant<Geometry> easVariant_;
   mutable Eigen::MatrixXd L_;
+  Eigen::VectorXd alpha_;
 
   //> CRTP
   const auto& underlying() const { return static_cast<const FE&>(*this); }
   auto& underlying() { return static_cast<FE&>(*this); }
 
+  /**
+   * \brief Initializes the internal state variable alpha_ based on the number of EAS parameters.
+   */
+  void initializeState() {
+    if (isDisplacementBased())
+      return;
+    alpha_.resize(numberOfEASParameters());
+    alpha_.setZero();
+  }
+
   template <typename ScalarType, int enhancedStrainSize>
   void calculateDAndLMatrix(
       const auto& easFunction, const auto& par, Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
@@ -267,6 +342,34 @@ private:
       }
     }
   }
+
+  template <typename ScalarType>
+  Eigen::VectorX<ScalarType> calculateRtilde(
+      const Requirement& par,
+      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+    using namespace Dune;
+    using namespace Dune::DerivativeDirections;
+    const auto geo      = underlying().localView().element().geometry();
+    auto strainFunction = underlying().strainFunction(par, dx);
+    const auto C        = underlying().materialTangentFunction(par);
+    Eigen::VectorX<ScalarType> Rtilde;
+    Rtilde.setZero(numberOfEASParameters());
+
+    auto calculateRtildeContribution = [&]<typename EAST>(const EAST& easFunction) {
+      for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
+        const auto M            = easFunction.calcM(gp.position());
+        const double intElement = geo.integrationElement(gp.position()) * gp.weight();
+        const auto EVoigt       = (strainFunction.evaluate(gpIndex, on(gridElement))).eval();
+        const auto CEval        = C(gpIndex);
+        auto stresses           = ((CEval * M * alpha_).template cast<ScalarType>()).eval();
+        stresses += underlying().getStress(EVoigt).eval();
+        Rtilde += (M.transpose() * stresses).eval() * intElement;
+      }
+    };
+
+    easVariant_(calculateRtildeContribution);
+    return Rtilde;
+  }
 };
 
 /**

ikarus/finiteelements/mechanics/kirchhoffloveshell.hh

@@ -234,6 +234,8 @@ protected:
     return kin;
   }
 
+  void updateStateImpl(const Requirement& /* par */, typename Traits::template VectorTypeConst<> /* correction */) {}
+
   template <typename ST>
   void calculateMatrixImpl(
       const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<ST> K,

ikarus/finiteelements/mechanics/linearelastic.hh

@@ -154,11 +154,24 @@ public:
     return [&]([[maybe_unused]] auto gp) { return materialTangent(); };
   }
 
+  /**
+   * \brief Get the stress for the given strain.
+   *
+   * \tparam ScalarType The scalar type for the material and strain.
+   * \tparam strainDim The dimension of the strain vector.
+   * \tparam voigt A boolean indicating whether to use the Voigt notation for stress.
+   * \param strain The strain vector.
+   * \return The stress vector calculated using the material's stresses function.
+   */
+  template <typename ScalarType, int strainDim, bool voigt = true>
+  auto getStress(const Eigen::Vector<ScalarType, strainDim>& strain) const {
+    return (materialTangent() * strain).eval();
+  }
+
   const Geometry& geometry() const { return *geo_; }
   [[nodiscard]] size_t numberOfNodes() const { return numberOfNodes_; }
   [[nodiscard]] int order() const { return order_; }
 
-public:
   /**
    * \brief Calculates a requested result at a specific local position.
    *
@@ -195,6 +208,8 @@ private:
   int order_{};
 
 protected:
+  void updateStateImpl(const Requirement& /* par */, typename Traits::template VectorTypeConst<> /* correction */) {}
+
   template <typename ScalarType>
   void calculateMatrixImpl(
       const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K,

ikarus/finiteelements/mechanics/loads/traction.hh

@@ -143,6 +143,8 @@ protected:
       const Requirement&, MatrixAffordance, typename Traits::template MatrixType<>,
       const std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>& = std::nullopt) const {}
 
+  void updateStateImpl(const Requirement& /* par */, typename Traits::template VectorTypeConst<> /* correction */) {}
+
 private:
   std::function<Eigen::Vector<double, worldDim>(const Dune::FieldVector<double, worldDim>&, const double&)>
       neumannBoundaryLoad_;

ikarus/finiteelements/mechanics/loads/volume.hh

@@ -111,6 +111,8 @@ protected:
       const Requirement& par, MatrixAffordance, typename Traits::template MatrixType<> K,
       const std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>& dx = std::nullopt) const {}
 
+  void updateStateImpl(const Requirement& /* par */, typename Traits::template VectorTypeConst<> /* correction */) {}
+
 private:
   std::function<Eigen::Vector<double, worldDim>(const Dune::FieldVector<double, worldDim>&, const double&)> volumeLoad_;
   //> CRTP

ikarus/finiteelements/mechanics/membranestrains.hh

@@ -7,7 +7,6 @@
  */
 
 #pragma once
-#include <ranges>
 
 #include <dune/common/fvector.hh>