refactor getDisplacementFunction

CHANGELOG.md

@@ -25,6 +25,7 @@ via `pip install pyikarus` ([#152](https://github.com/ikarus-project/ikarus/pull
 - Added getRawMatrix/getRawVector functionality to the assemblers ([#179](https://github.com/ikarus-project/ikarus/pull/179))
 - 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))
+- Refactored load functions and `getDisplacementFunction` in finite elements ([#221](https://github.com/ikarus-project/ikarus/pull/221))
 
 - Improve material library and Python bindings ([#186](https://github.com/ikarus-project/ikarus/pull/176)), default e.g.
   `StVenantKirchhoff` is

docs/CMakeLists.txt

@@ -2,26 +2,26 @@ cmake_minimum_required(VERSION 3.18)
 project(Ikarus-Docs)
 
 find_program(
-        MKDOCS_EXECUTABLE
-        NAMES mkdocs
-        DOC "MkDocs documentation generation tool (http://www.mkdocs.org)" REQUIRED
+  MKDOCS_EXECUTABLE
+  NAMES mkdocs
+  DOC "MkDocs documentation generation tool (http://www.mkdocs.org)" REQUIRED
 )
 
 add_subdirectory(website/doxygen)
 
 add_custom_target(
   site
   COMMAND xvfb-run -a mkdocs build --config-file mkdocs.insiders.yml
-        WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}"
+  WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}"
 )
 
 add_dependencies(site doxygen_ikarus)
 
 add_custom_target(
-        localSite
-        COMMAND ${MKDOCS_EXECUTABLE} serve
-        WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}"
-        DEPENDS mkdocs-macros.py mkdocs.yml
+  localSite
+  COMMAND ${MKDOCS_EXECUTABLE} serve
+  WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}"
+  DEPENDS mkdocs-macros.py mkdocs.yml
 )
 
 add_dependencies(localSite doxygen_ikarus)

ikarus/finiteelements/CMakeLists.txt

@@ -2,7 +2,7 @@
 # SPDX-License-Identifier: LGPL-3.0-or-later
 
 # install headers
-install(FILES ferequirements.hh fetraits.hh feparameter.hh physicshelper.hh
+install(FILES ferequirements.hh fetraits.hh fehelper.hh physicshelper.hh
         DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements
 )
 

ikarus/finiteelements/fehelper.hh

@@ -0,0 +1,34 @@
+// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#pragma once
+
+#include <dune/localfefunctions/impl/standardLocalFunction.hh>
+#include <dune/localfefunctions/manifolds/realTuple.hh>
+
+#include <ikarus/finiteelements/ferequirements.hh>
+
+namespace Ikarus::FEHelper {
+  template <typename Traits, typename FERequirementType, typename DisplacementBasedElement, typename ScalarType>
+  auto getDisplacementFunctionImpl(const DisplacementBasedElement ele, const FERequirementType& par,
+                                   const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) {
+    const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement);
+    Dune::BlockVector<Dune::RealTuple<ScalarType, Traits::worlddim>> disp(ele.dispAtNodes.size());
+
+    if (dx)
+      for (auto i = 0U; i < disp.size(); ++i)
+        for (auto j = 0U; j < Traits::worlddim; ++j)
+          disp[i][j] = dx.value()[i * Traits::worlddim + j]
+                       + d[ele.localView().index(ele.localView().tree().child(j).localIndex(i))[0]];
+    else
+      for (auto i = 0U; i < disp.size(); ++i)
+        for (auto j = 0U; j < Traits::worlddim; ++j)
+          disp[i][j] = d[ele.localView().index(ele.localView().tree().child(j).localIndex(i))[0]];
+
+    auto geo
+        = std::make_shared<const decltype(ele.localView().element().geometry())>(ele.localView().element().geometry());
+
+    Dune::StandardLocalFunction uFunction(ele.localBasis, disp, geo);
+    return std::make_pair(uFunction, disp);
+  }
+}  // namespace Ikarus::FEHelper

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -6,6 +6,7 @@
 #  include <dune/localfefunctions/derivativetransformators.hh>
 #  include <dune/localfefunctions/meta.hh>
 
+#  include <ikarus/finiteelements/fehelper.hh>
 #  include <ikarus/finiteelements/ferequirements.hh>
 #  include <ikarus/utils/eigendunetransformations.hh>
 
@@ -322,12 +323,12 @@ 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;
+      const auto& d             = req.getGlobalSolution(Ikarus::FESolutions::displacement);
+      const auto C              = DisplacementBasedElement::getMaterialTangentFunction(req.getFERequirements());
+      const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes;
 
       Eigen::VectorXd disp(localView().size());
-      for (auto i = 0U; i < numNodes; ++i)
+      for (auto i = 0U; i < numberOfNodes; ++i)
         for (auto k2 = 0U; k2 < Traits::mydim; ++k2)
           disp[i * Traits::mydim + k2] = d[localView().index(localView().tree().child(k2).localIndex(i))[0]];
 
@@ -354,8 +355,8 @@ namespace Ikarus {
     }
 
     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");
 
@@ -411,21 +412,10 @@ 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]];
+      auto strainFunction       = DisplacementBasedElement::getStrainFunction(par, dx);
+      const auto& dispInBlock   = DisplacementBasedElement::getDisplacementFunction(par, dx).second;
+      const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes;
+      const auto disp           = Ikarus::convertDuneBlockToEigenVector(dispInBlock);
 
       using namespace Dune::DerivativeDirections;
 
@@ -447,9 +437,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;
                 }
               }
             }
@@ -469,18 +460,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::getStrainFunction(par);
+      const auto C             = DisplacementBasedElement::getMaterialTangentFunction(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)

ikarus/finiteelements/mechanics/kirchhoffloveshell.hh

@@ -15,6 +15,7 @@
 #include <autodiff/forward/dual/eigen.hpp>
 
 #include <ikarus/finiteelements/febases/powerbasisfe.hh>
+#include <ikarus/finiteelements/fehelper.hh>
 #include <ikarus/finiteelements/ferequirements.hh>
 #include <ikarus/finiteelements/fetraits.hh>
 #include <ikarus/finiteelements/mechanics/materials.hh>
@@ -53,7 +54,7 @@ namespace Ikarus {
     using GridView                = typename FlatBasis::GridView;
     using Traits                  = TraitsFromLocalView<LocalView, useEigenRef>;
     static constexpr int myDim    = Traits::mydim;
-    static constexpr int worlddim = Traits::worlddim;
+    static constexpr int worldDim = Traits::worlddim;
 
     template <typename VolumeLoad = LoadDefault, typename NeumannBoundaryLoad = LoadDefault>
     KirchhoffLoveShell(const Basis& globalBasis, const typename LocalView::Element& element, double emod, double nu,
@@ -69,7 +70,7 @@ namespace Ikarus {
       auto& first_child = this->localView().tree().child(0);
       const auto& fe    = first_child.finiteElement();
       numberOfNodes     = fe.size();
-      dispAtNodes.resize(fe.size());
+      dispAtNodes.resize(numberOfNodes);
       order      = 2 * (fe.localBasis().order());
       localBasis = Dune::CachedLocalBasis(fe.localBasis());
       if constexpr (requires { this->localView().element().impl().getQuadratureRule(order); })
@@ -93,23 +94,7 @@ namespace Ikarus {
     template <typename ScalarType = double>
     auto getDisplacementFunction(const FERequirementType& par,
                                  const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
-      const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement);
-
-      Dune::BlockVector<Dune::RealTuple<ScalarType, Traits::worlddim>> disp(dispAtNodes.size());
-      if (dx)
-        for (auto i = 0U; i < disp.size(); ++i)
-          for (auto k2 = 0U; k2 < worlddim; ++k2)
-            disp[i][k2] = dx.value()[i * worlddim + k2]
-                          + d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]];
-      else
-        for (auto i = 0U; i < disp.size(); ++i)
-          for (auto k2 = 0U; k2 < worlddim; ++k2)
-            disp[i][k2] = d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]];
-
-      auto geo = std::make_shared<const typename GridView::GridView::template Codim<0>::Entity::Geometry>(
-          this->localView().element().geometry());
-      Dune::StandardLocalFunction uFunction(localBasis, disp, geo);
-      return std::make_pair(uFunction, disp);
+      return Ikarus::FEHelper::getDisplacementFunctionImpl<Traits>(*this, par, dx);
     }
 
     inline double calculateScalar(const FERequirementType& par) const { return calculateScalarImpl<double>(par); }
@@ -123,14 +108,11 @@ namespace Ikarus {
     Dune::CachedLocalBasis<
         std::remove_cvref_t<decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis())>>
         localBasis;
-    std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
-                                                          const double&)>
-        volumeLoad;
-    std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
-                                                          const double&)>
+    std::function<Eigen::Vector<double, worldDim>(const Eigen::Vector<double, worldDim>&, const double&)> volumeLoad;
+    std::function<Eigen::Vector<double, worldDim>(const Eigen::Vector<double, worldDim>&, const double&)>
         neumannBoundaryLoad;
     const BoundaryPatch<GridView>* neumannBoundary;
-    mutable Dune::BlockVector<Dune::RealTuple<double, Traits::dimension>> dispAtNodes;
+    mutable Dune::BlockVector<Dune::RealTuple<double, worldDim>> dispAtNodes;
     double emod_;
     double nu_;
     double thickness_;

ikarus/finiteelements/mechanics/linearelastic.hh

@@ -19,6 +19,7 @@
 #  include <autodiff/forward/dual/eigen.hpp>
 
 #  include <ikarus/finiteelements/febases/powerbasisfe.hh>
+#  include <ikarus/finiteelements/fehelper.hh>
 #  include <ikarus/finiteelements/ferequirements.hh>
 #  include <ikarus/finiteelements/fetraits.hh>
 #  include <ikarus/finiteelements/mechanics/materials.hh>
@@ -28,7 +29,6 @@
 #  include <ikarus/utils/linearalgebrahelper.hh>
 
 namespace Ikarus {
-
   template <typename Basis_, typename FERequirements_ = FErequirements<>, bool useEigenRef = false>
   class LinearElastic : public PowerBasisFE<typename Basis_::FlatBasis> {
   public:
@@ -43,7 +43,8 @@ namespace Ikarus {
 
     using Traits = TraitsFromLocalView<LocalView, useEigenRef>;
 
-    static constexpr int myDim = Traits::mydim;
+    static constexpr int myDim    = Traits::mydim;
+    static constexpr int worldDim = Traits::worlddim;
 
     template <typename VolumeLoad = LoadDefault, typename NeumannBoundaryLoad = LoadDefault>
     LinearElastic(const Basis& globalBasis, const typename LocalView::Element& element, double emod, double nu,
@@ -74,34 +75,16 @@ namespace Ikarus {
              && "If you pass a Neumann boundary you should also pass the function for the Neumann load!");
     }
 
-  public:
     template <typename ScalarType>
     auto getDisplacementFunction(const FERequirementType& par,
                                  const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
-      const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement);
-
-      Dune::BlockVector<Dune::RealTuple<ScalarType, Traits::dimension>> disp(dispAtNodes.size());
-      if (dx) {
-        for (auto i = 0U; i < disp.size(); ++i)
-          for (auto k2 = 0U; k2 < myDim; ++k2)
-            disp[i][k2] = dx.value()[i * myDim + k2]
-                          + d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]];
-      } else
-        for (auto i = 0U; i < disp.size(); ++i)
-          for (auto k2 = 0U; k2 < myDim; ++k2)
-            disp[i][k2] = d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]];
-
-      auto geo = std::make_shared<const typename GridView::GridView::template Codim<0>::Entity::Geometry>(
-          this->localView().element().geometry());
-      Dune::StandardLocalFunction uFunction(localBasis, disp, geo);
-
-      return uFunction;
+      return Ikarus::FEHelper::getDisplacementFunctionImpl<Traits>(*this, par, dx);
     }
 
     template <class ScalarType = double>
     auto getStrainFunction(const FERequirementType& par,
                            const std::optional<const Eigen::VectorX<ScalarType>>& dx = std::nullopt) const {
-      return linearStrains(getDisplacementFunction(par, dx));
+      return linearStrains(getDisplacementFunction(par, dx).first);
     }
 
     auto getMaterialTangent() const {
@@ -162,14 +145,11 @@ namespace Ikarus {
     Dune::CachedLocalBasis<
         std::remove_cvref_t<decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis())>>
         localBasis;
-    std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
-                                                          const double&)>
-        volumeLoad;
-    std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
-                                                          const double&)>
+    std::function<Eigen::Vector<double, worldDim>(const Eigen::Vector<double, worldDim>&, const double&)> volumeLoad;
+    std::function<Eigen::Vector<double, worldDim>(const Eigen::Vector<double, worldDim>&, const double&)>
         neumannBoundaryLoad;
     const BoundaryPatch<GridView>* neumannBoundary;
-    mutable Dune::BlockVector<Dune::RealTuple<double, Traits::dimension>> dispAtNodes;
+    mutable Dune::BlockVector<Dune::RealTuple<double, worldDim>> dispAtNodes;
     double emod_;
     double nu_;
     size_t numberOfNodes{0};