Refactor getDisplacementFunction (#223)

ikarus-project · Jan 17, 2024 · cca5ed0 · cca5ed0
1 parent 17bb21e
commit cca5ed0
Show file tree

Hide file tree

Showing 10 changed files with 222 additions and 248 deletions.
diff --git 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
@@ -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.

diff --git 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
@@ -0,0 +1,42 @@
+// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#pragma once
+
+#include <dune/localfefunctions/manifolds/realTuple.hh>
+
+#include <ikarus/finiteelements/ferequirements.hh>
+#include <ikarus/finiteelements/physicshelper.hh>
+
+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 <typename FERequirementType, typename LocalView, typename ScalarType>
+  auto localSolutionBlockVector(const typename FERequirementType::SolutionVectorTypeRaw& x, const LocalView& localView,
+                                const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
+    using Traits           = TraitsFromLocalView<LocalView>;
+    constexpr int worldDim = Traits::worlddim;
+    const auto& fe         = localView.tree().child(0).finiteElement();
+    Dune::BlockVector<Dune::RealTuple<ScalarType, worldDim>> 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
@@ -10,8 +10,10 @@
 #pragma once
 #if HAVE_DUNE_LOCALFEFUNCTIONS
 #  include <dune/localfefunctions/derivativetransformators.hh>
+#  include <dune/localfefunctions/linearAlgebraHelper.hh>
 #  include <dune/localfefunctions/meta.hh>
 
+#  include <ikarus/finiteelements/fehelper.hh>
 #  include <ikarus/finiteelements/ferequirements.hh>
 #  include <ikarus/utils/eigendunetransformations.hh>
 
@@ -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(
           [&]<typename EAST>(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<ScalarType> 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>(Traits::mydim * i) += bopI.transpose() * stresses * intElement;
+                  force.template segment<Traits::worlddim>(Traits::worlddim * i)
+                      += bopI.transpose() * stresses * intElement;
                 }
               }
             }
@@ -619,18 +608,18 @@ 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()) {
         const auto M      = easFunction.calcM(gp.position());
         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<enhancedStrainSize, Traits::worlddim>(0, I)