add nonlinear version of eas

Author: tarun-mitruka

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -137,27 +137,74 @@ public:
     if (numberOfEASParameters != 0)
       easApplicabilityCheck();
     easVariant_.setEASType(numberOfEASParameters);
+    initializeState();
   }
 
+  /**
+   * \brief Initializes the internal parameter alpha_ based on the number of EAS parameters.
+   */
+  void initializeState() {
+    if (isDisplacementBased())
+      return;
+    alpha_.resize(numberOfEASParameters());
+    alpha_.setZero();
+  }
+
+  /**
+   * \brief Updates the internal parameter alpha_ at the end of an iteration
+   * when NonLinearSolverMessages::SOLUTION_UPDATED is notified by the non-linear solver.
+   *
+   * \param par The Requirement object.
+   */
+  template <typename ScalarType = double>
+  void updateState(const Requirement& par, const Eigen::VectorXd& correction_) {
+    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);
+  }
+
+  const auto& alpha() const { return alpha_; }
+
 protected:
   void bindImpl() {
     assert(underlying().localView().bound());
     easVariant_.bind(underlying().localView().element().geometry());
+    initializeState();
   }
 
-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();
@@ -171,10 +218,30 @@ protected:
         return Eigen::MatrixX<ScalarType>{};
     }();
 
+    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();
     };
@@ -196,45 +263,44 @@ 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 (std::is_same_v<ScalarType, double>)
-          return [this]() -> Eigen::MatrixXd& { return L_; }();
-        else
-          return Eigen::MatrixX<ScalarType>{};
-      }();
-      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); }
@@ -266,6 +332,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/linearelastic.hh

@@ -154,11 +154,15 @@ public:
     return [&]([[maybe_unused]] auto gp) { return materialTangent(); };
   }
 
+  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.
    *

ikarus/finiteelements/mechanics/nonlinearelastic.hh

@@ -160,6 +160,22 @@ public:
     }
   }
 
+  /**
+   * \brief Gets the material tangent function for the given Requirement.
+   *
+   * \param par The Requirement object.
+   * \return The material tangent function.
+   */
+  auto materialTangentFunction([[maybe_unused]] const Requirement& par) const {
+    return [&]([[maybe_unused]] auto gpIndex) {
+      using namespace Dune::DerivativeDirections;
+      using namespace Dune;
+      const auto eps    = strainFunction(par);
+      const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
+      return materialTangent(EVoigt);
+    };
+  }
+
   /**
    * \brief Get the internal energy for the given strain.
    *

tests/src/CMakeLists.txt

@@ -27,6 +27,7 @@ set(programSourceFiles
     testlinearelasticity.cpp
     testmanifolds.cpp
     testmaterial.cpp
+    testnonlineareas.cpp
     testnonlinearelasticityneohooke.cpp
     testnonlinearelasticitysvk.cpp
     testnonlinearoperator.cpp

tests/src/testnonlineareas.cpp

@@ -0,0 +1,135 @@
+// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#include <config.h>
+
+#include "testcommon.hh"
+
+#include <dune/common/test/testsuite.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <ikarus/finiteelements/autodifffe.hh>
+#include <ikarus/finiteelements/fefactory.hh>
+#include <ikarus/finiteelements/mechanics/enhancedassumedstrains.hh>
+#include <ikarus/finiteelements/mechanics/materials.hh>
+#include <ikarus/finiteelements/mechanics/nonlinearelastic.hh>
+#include <ikarus/utils/basis.hh>
+#include <ikarus/utils/init.hh>
+
+using namespace Ikarus;
+using namespace Dune::Indices;
+using Dune::TestSuite;
+
+template <int gridDim, typename GV, typename MAT, typename BH, typename VectorType>
+auto nonlinearEASAutoDiffTest(const GV& gridView, const MAT& mat, const BH& basis, VectorType& d,
+                              const int numberOfEASParameters) {
+  TestSuite t("Nonlinear EAS AutoDiff with gridDim = " + std::to_string(gridDim) +
+              " and numberOfEASParameters = " + std::to_string(numberOfEASParameters));
+  double lambda = 7.3;
+  auto sk       = skills(nonLinearElastic(mat), eas(numberOfEASParameters));
+  auto fe       = Ikarus::makeFE(basis, sk);
+  using FE      = decltype(fe);
+  auto req      = typename FE::Requirement();
+  req.insertGlobalSolution(d).insertParameter(lambda);
+
+  auto affordance = AffordanceCollections::elastoStatics;
+
+  for (auto element : elements(gridView)) {
+    auto localView = basis.flat().localView();
+    localView.bind(element);
+    auto nDOF        = localView.size();
+    const double tol = 1e-10;
+
+    fe.bind(element);
+    fe.updateState(req, d); // here d = correction vector (DeltaD)
+
+    const std::string feClassName = Dune::className(fe);
+
+    using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe), true>;
+    AutoDiffBasedFE feAutoDiff(fe);
+
+    Eigen::MatrixXd K, KAutoDiff;
+    K.setZero(nDOF, nDOF);
+    KAutoDiff.setZero(nDOF, nDOF);
+
+    calculateMatrix(fe, req, affordance.matrixAffordance(), K);
+    calculateMatrix(feAutoDiff, req, affordance.matrixAffordance(), KAutoDiff);
+
+    t.check(fe.numberOfEASParameters() == feAutoDiff.realFE().numberOfEASParameters())
+        << "Number of EAS parameters for FE(" << fe.numberOfEASParameters()
+        << ") and number of EAS parameters for AutodiffFE(" << feAutoDiff.realFE().numberOfEASParameters()
+        << ") are not equal";
+
+    t.check(fe.alpha().isApprox(feAutoDiff.alpha(), tol), "Mismatch in alpha.")
+        << "alpha is\n"
+        << fe.alpha().transpose() << "\n alphaAutoDiff is\n"
+        << feAutoDiff.alpha();
+
+    t.check(K.isApprox(KAutoDiff, tol),
+            "Mismatch between the stiffness matrices obtained from explicit implementation and the one based on "
+            "automatic differentiation.")
+        << "K is \n"
+        << K << "\nKAutoDiff is \n"
+        << KAutoDiff << "\nThe difference is\n"
+        << (K - KAutoDiff);
+
+    if (numberOfEASParameters != 0) {
+      t.checkThrow<Dune::NotImplemented>(
+          [&]() {
+            Eigen::VectorXd RAutoDiff;
+            RAutoDiff.setZero(nDOF);
+            calculateVector(feAutoDiff, req, affordance.vectorAffordance(), RAutoDiff);
+          },
+          "calculateVector with AutoDiff should throw a Dune::NotImplemented");
+
+      t.checkThrow<Dune::NotImplemented>(
+          [&]() { auto energyAutoDiff = calculateScalar(feAutoDiff, req, affordance.scalarAffordance()); },
+          "calculateScalar with AutoDiff should throw a Dune::NotImplemented");
+    }
+  }
+  return t;
+}
+
+template <int gridDim, typename TestSuitType>
+void easAutoDiffTest(const TestSuitType& t) {
+  std::array<int, (gridDim == 2) ? 4 : 3> easParameters;
+  if constexpr (gridDim == 2)
+    easParameters = {0, 4, 5, 7};
+  else
+    easParameters = {0, 9, 21};
+
+  auto grid = createUGGridFromCorners<gridDim>(CornerDistortionFlag::randomlyDistorted);
+  grid->globalRefine(2);
+  auto gridView     = grid->leafGridView();
+  auto matParameter = toLamesFirstParameterAndShearModulus({.emodul = 100.0, .nu = 0.3});
+  StVenantKirchhoff matSVK(matParameter);
+  auto reducedMat = planeStress(matSVK);
+
+  using namespace Dune::Functions::BasisFactory;
+  auto basis = Ikarus::makeBasis(gridView, power<gridDim>(lagrange<1>()));
+  Eigen::VectorXd d;
+
+  auto autoDiffTestFunctor = [&](int numberOfEASParameters) {
+    if constexpr (gridDim == 2)
+      nonlinearEASAutoDiffTest<gridDim>(gridView, reducedMat, basis, d, numberOfEASParameters);
+    else
+      nonlinearEASAutoDiffTest<gridDim>(gridView, matSVK, basis, d, numberOfEASParameters);
+  };
+
+  for (const int numberOfEASParameters : easParameters) {
+    d.setZero(basis.flat().size());
+    autoDiffTestFunctor(numberOfEASParameters);
+    d.setRandom(basis.flat().size());
+    autoDiffTestFunctor(numberOfEASParameters);
+  }
+}
+
+int main(int argc, char** argv) {
+  Ikarus::init(argc, argv);
+  Dune::TestSuite t("Nonlinear EAS Test");
+  easAutoDiffTest<2>(t);
+  easAutoDiffTest<3>(t);
+  return t.exit();
+}