refactor

Author: rath3t

ikarus/finiteelements/mechanics/kirchhoffloveshell.hh

@@ -36,20 +36,20 @@ namespace Ikarus {
   template <typename Basis_, typename FERequirements_ = FErequirements<>, bool useEigenRef = false>
   class KirchhoffLoveShell : public PowerBasisFE<typename Basis_::FlatBasis> {
   public:
-    using Basis                   = Basis_;
-    using FlatBasis               = typename Basis::FlatBasis;
-    using BasePowerFE             = PowerBasisFE<FlatBasis>;  // Handles globalIndices function
-    using FERequirementType       = FERequirements_;
-    using ResultRequirementsType  = ResultRequirements<FERequirementType>;
-    using LocalView               = typename FlatBasis::LocalView;
-    using Element                 = typename LocalView::Element;
-    using Geometry                = typename Element::Geometry;
-    using GridView                = typename FlatBasis::GridView;
-    using Traits                  = TraitsFromLocalView<LocalView, useEigenRef>;
-    static constexpr int myDim    = Traits::mydim;
-    static constexpr int worlddim = Traits::worlddim;
+    using Basis                             = Basis_;
+    using FlatBasis                         = typename Basis::FlatBasis;
+    using BasePowerFE                       = PowerBasisFE<FlatBasis>;  // Handles globalIndices function
+    using FERequirementType                 = FERequirements_;
+    using ResultRequirementsType            = ResultRequirements<FERequirementType>;
+    using LocalView                         = typename FlatBasis::LocalView;
+    using Element                           = typename LocalView::Element;
+    using Geometry                          = typename Element::Geometry;
+    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 membraneStrainSize = 3;
-    static constexpr int bendingStrainSize = 3;
+    static constexpr int bendingStrainSize  = 3;
 
     /**
      * @brief Constructor for the KirchhoffLoveShell class.
@@ -203,82 +203,81 @@ namespace Ikarus {
 
   protected:
     /**
- * \brief Compute material properties and strains at a given integration point.
- *
- * \param gpPos The position of the integration point.
- * \param gpIndex The index of the integration point.
- * \param geo The geometry object providing position and derivatives.
- * \param uFunction The function representing the displacement field.
- *
- * \tparam gpPos The type of the integration point position.
- * \tparam gpIndex The type of the integration point index.
- * \tparam geo The type of the geometry object.
- * \tparam uFunction The type of the displacement field function.
- *
- * \return A tuple containing the material tangent, membrane strain, curvature variation,
- *         Jacobian matrix, metric tensor, Hessian matrix, second fundamental form,
- *         normalized normal vector, and normal vector at the given integration point.
- */
-auto computeMaterialAndStrains(const Dune::FieldVector<double, 2> &gpPos, int gpIndex, const Geometry &geo,
-                               const auto &uFunction) const {
-    // Deduce the scalar type of the function
-    using ScalarType = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
-
-    using namespace Dune;
-    using namespace Dune::DerivativeDirections;
+     * \brief Compute material properties and strains at a given integration point.
+     *
+     * \param gpPos The position of the integration point.
+     * \param gpIndex The index of the integration point.
+     * \param geo The geometry object providing position and derivatives.
+     * \param uFunction The function representing the displacement field.
+     *
+     * \tparam gpPos The type of the integration point position.
+     * \tparam gpIndex The type of the integration point index.
+     * \tparam geo The type of the geometry object.
+     * \tparam uFunction The type of the displacement field function.
+     *
+     * \return A tuple containing the material tangent, membrane strain, curvature variation,
+     *         Jacobian matrix, metric tensor, Hessian matrix, second fundamental form,
+     *         normalized normal vector, and normal vector at the given integration point.
+     */
+    auto computeMaterialAndStrains(const Dune::FieldVector<double, 2> &gpPos, int gpIndex, const Geometry &geo,
+                                   const auto &uFunction) const {
+      // Deduce the scalar type of the function
+      using ScalarType = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
 
-    // Extract position and derivatives from the geometry
-    const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos);
-    const auto J           = toEigen(Jd);
-    const auto H           = toEigen(Hd);
+      using namespace Dune;
+      using namespace Dune::DerivativeDirections;
 
-    // Compute the metric tensor A
-    const Eigen::Matrix<double, 2, 2> A = J * J.transpose();
+      // Extract position and derivatives from the geometry
+      const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos);
+      const auto J           = toEigen(Jd);
+      const auto H           = toEigen(Hd);
 
-    // Build the metric tensor G in 3D
-    Eigen::Matrix<double, worlddim, worlddim> G;
-    G.setZero();
-    G.template block<2, 2>(0, 0) = A;
-    G(2, 2)             = 1;
+      // Compute the metric tensor A
+      const Eigen::Matrix<double, 2, 2> A = J * J.transpose();
 
-    const Eigen::Matrix<double, worlddim, worlddim> GInv = G.inverse();
+      // Build the metric tensor G in 3D
+      Eigen::Matrix<double, worlddim, worlddim> G;
+      G.setZero();
+      G.template block<2, 2>(0, 0) = A;
+      G(2, 2)                      = 1;
 
-    // Compute the material tangent
-    const auto C = materialTangent(GInv);
+      const Eigen::Matrix<double, worlddim, worlddim> GInv = G.inverse();
 
-    // Compute membrane strain
-    Eigen::Vector3<ScalarType> epsV = membraneStrain.value(gpPos, geo, uFunction);
+      // Compute the material tangent
+      const auto C = materialTangent(GInv);
 
-    const auto &Ndd = localBasis.evaluateSecondDerivatives(gpIndex);
+      // Compute membrane strain
+      Eigen::Vector3<ScalarType> epsV = membraneStrain.value(gpPos, geo, uFunction);
 
-    const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef());
+      const auto &Ndd = localBasis.evaluateSecondDerivatives(gpIndex);
 
-    // Compute the Hessian of the deformed midsurface
-    const auto hessianu = Ndd.transpose().template cast<ScalarType>() * uasMatrix;
-    const Eigen::Matrix3<ScalarType> h = H + hessianu;
+      const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef());
 
-    // Compute the gradient and Jacobian matrices
-    const Eigen::Matrix<ScalarType, worlddim, myDim> gradu =
-        toEigen(uFunction.evaluateDerivative(gpIndex, Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement)));
-    const Eigen::Matrix<ScalarType, myDim, worlddim> j = J + gradu.transpose();
+      // Compute the Hessian of the deformed midsurface
+      const auto hessianu                = Ndd.transpose().template cast<ScalarType>() * uasMatrix;
+      const Eigen::Matrix3<ScalarType> h = H + hessianu;
 
-    // Compute the normal vector
-    const Eigen::Vector3<ScalarType> a3N = (j.row(0).cross(j.row(1))).normalized();
+      // Compute the gradient and Jacobian matrices
+      const Eigen::Matrix<ScalarType, worlddim, myDim> gradu = toEigen(uFunction.evaluateDerivative(
+          gpIndex, Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement)));
+      const Eigen::Matrix<ScalarType, myDim, worlddim> j     = J + gradu.transpose();
 
-    // Compute the normalized normal vector a3
-    const Eigen::Vector3<ScalarType> a3 = a3N.normalized();
+      // Compute the normal vector
+      const Eigen::Vector3<ScalarType> a3N = (j.row(0).cross(j.row(1))).normalized();
 
-    // Compute the vector bV
-    Eigen::Vector<ScalarType, worlddim> bV = h * a3;
-    bV(2) *= 2;  // Voigt notation requires the factor of 2 here
+      // Compute the normalized normal vector a3
+      const Eigen::Vector3<ScalarType> a3 = a3N.normalized();
 
-    // Compute BV and kappaV
-    const auto BV     = toVoigt(toEigen(geo.impl().secondFundamentalForm(gpPos)));
-    const auto kappaV = (BV - bV).eval();
+      // Compute the vector bV
+      Eigen::Vector<ScalarType, worlddim> bV = h * a3;
+      bV(2) *= 2;  // Voigt notation requires the factor of 2 here
 
-    return std::make_tuple(C, epsV, kappaV, j, J, h, H, a3N, a3);
-}
+      // Compute BV and kappaV
+      const auto BV     = toVoigt(toEigen(geo.impl().secondFundamentalForm(gpPos)));
+      const auto kappaV = (BV - bV).eval();
 
+      return std::make_tuple(C, epsV, kappaV, j, J, h, H, a3N, a3);
+    }
 
     template <typename ScalarType>
     void calculateMatrixImpl(const FERequirementType &par, typename Traits::template MatrixType<ScalarType> K,
@@ -296,7 +295,7 @@ auto computeMaterialAndStrains(const Dune::FieldVector<double, 2> &gpPos, int gp
         const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3]
             = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction);
         const Eigen::Vector<ScalarType, membraneStrainSize> membraneForces = thickness_ * C * epsV;
-        const Eigen::Vector<ScalarType, bendingStrainSize> moments        = Dune::power(thickness_, 3) / 12.0 * C * kappaV;
+        const Eigen::Vector<ScalarType, bendingStrainSize> moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV;
 
         const auto &Nd  = localBasis.evaluateJacobian(gpIndex);
         const auto &Ndd = localBasis.evaluateSecondDerivatives(gpIndex);
@@ -306,18 +305,17 @@ auto computeMaterialAndStrains(const Dune::FieldVector<double, 2> &gpPos, int gp
 
           Eigen::Matrix<ScalarType, bendingStrainSize, worlddim> bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3);
           for (size_t j = i; j < numberOfNodes; ++j) {
-            auto KBlock = K.template block<worlddim,worlddim>(worlddim * i, worlddim * j);
+            auto KBlock = K.template block<worlddim, worlddim>(worlddim * i, worlddim * j);
             Eigen::Matrix<ScalarType, membraneStrainSize, worlddim> bopJMembrane
                 = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, j);
             Eigen::Matrix<ScalarType, bendingStrainSize, worlddim> bopJBending = bopBending(jE, h, Nd, Ndd, j, a3N, a3);
-            KBlock
-                += thickness_ * bopIMembrane.transpose() * C * bopJMembrane * intElement;
-            KBlock
-                += Dune::power(thickness_, 3) / 12.0 * bopIBending.transpose() * C * bopJBending * intElement;
+            KBlock += thickness_ * bopIMembrane.transpose() * C * bopJMembrane * intElement;
+            KBlock += Dune::power(thickness_, 3) / 12.0 * bopIBending.transpose() * C * bopJBending * intElement;
 
-            Eigen::Matrix<ScalarType, worlddim,worlddim> kgMembraneIJ
+            Eigen::Matrix<ScalarType, worlddim, worlddim> kgMembraneIJ
                 = membraneStrain.secondDerivative(gp.position(), Nd, geo, uFunction, localBasis, membraneForces, i, j);
-            Eigen::Matrix<ScalarType, worlddim, worlddim> kgBendingIJ = kgBending(jE, h, Nd, Ndd, a3N, a3, moments, i, j);
+            Eigen::Matrix<ScalarType, worlddim, worlddim> kgBendingIJ
+                = kgBending(jE, h, Nd, Ndd, a3N, a3, moments, i, j);
             KBlock += kgMembraneIJ * intElement;
             KBlock += kgBendingIJ * intElement;
           }

ikarus/finiteelements/mechanics/membranestrains.hh

@@ -11,23 +11,21 @@
 
 #include <dune/common/fvector.hh>
 
-
 #include <Eigen/Core>
 namespace Ikarus {
 
   struct DefaultMembraneStrain {
-
     /**
- * \brief Compute the strain vector at a given integration point.
- *
- * \param gpPos The position of the integration point.
- * \param geo The geometry object providing the transposed Jacobian.
- * \param uFunction The function representing the displacement field.
- *
- * \tparam Geometry The type of the geometry object.
- *
- * \return The strain vector at the given integration point.
- */
+     * \brief Compute the strain vector at a given integration point.
+     *
+     * \param gpPos The position of the integration point.
+     * \param geo The geometry object providing the transposed Jacobian.
+     * \param uFunction The function representing the displacement field.
+     *
+     * \tparam Geometry The type of the geometry object.
+     *
+     * \return The strain vector at the given integration point.
+     */
     template <typename Geometry>
     auto value(const Dune::FieldVector<double, 2> &gpPos, const Geometry &geo, const auto &uFunction) const
         -> Eigen::Vector3<typename std::remove_cvref_t<decltype(uFunction)>::ctype> {
@@ -47,21 +45,21 @@ namespace Ikarus {
     }
 
     /**
- * \brief Compute the strain-displacement matrix for a given node and integration point.
- *
- * \param gpPos The position of the integration point.
- * \param jcur The Jacobian matrix.
- * \param dNAtGp The derivative of the shape functions at the integration point.
- * \param geo The geometry object of the finite element.
- * \param uFunction The function representing the displacement field.
- * \param localBasis The local basis object.
- * \param node The FE node index.
- *
- * \tparam Geometry The type of the geometry object.
- * \tparam ScalarType The scalar type for the matrix elements.
- *
- * \return The strain-displacement matrix for the given node and integration point.
- */
+     * \brief Compute the strain-displacement matrix for a given node and integration point.
+     *
+     * \param gpPos The position of the integration point.
+     * \param jcur The Jacobian matrix.
+     * \param dNAtGp The derivative of the shape functions at the integration point.
+     * \param geo The geometry object of the finite element.
+     * \param uFunction The function representing the displacement field.
+     * \param localBasis The local basis object.
+     * \param node The FE node index.
+     *
+     * \tparam Geometry The type of the geometry object.
+     * \tparam ScalarType The scalar type for the matrix elements.
+     *
+     * \return The strain-displacement matrix for the given node and integration point.
+     */
     template <typename Geometry, typename ScalarType>
     auto derivative(const Dune::FieldVector<double, 2> &gpPos, const Eigen::Matrix<ScalarType, 2, 3> &jcur,
                     const auto &dNAtGp, const Geometry &geo, const auto &uFunction, const auto &localBasis,
@@ -75,23 +73,24 @@ namespace Ikarus {
     }
 
     /**
- * \brief Compute the second derivative of the membrane strain for a given node pair and integration point.
- * \details This function computes the geometric tangent stiffness for a given node pair at a given integration point.
- *
- * \param gpPos The position of the integration point.
- * \param dNAtGp The derivative of the shape functions at the integration point.
- * \param geo The geometry object.
- * \param uFunction The function representing the displacement field.
- * \param localBasis The local basis object.
- * \param S The strain vector.
- * \param I The index of the first node.
- * \param J The index of the second node.
- *
- * \tparam Geometry The type of the geometry object.
- * \tparam ScalarType The scalar type for the matrix elements.
- *
- * \return The second derivative of the membrane strain.
- */
+     * \brief Compute the second derivative of the membrane strain for a given node pair and integration point.
+     * \details This function computes the geometric tangent stiffness for a given node pair at a given integration
+     * point.
+     *
+     * \param gpPos The position of the integration point.
+     * \param dNAtGp The derivative of the shape functions at the integration point.
+     * \param geo The geometry object.
+     * \param uFunction The function representing the displacement field.
+     * \param localBasis The local basis object.
+     * \param S The strain vector.
+     * \param I The index of the first node.
+     * \param J The index of the second node.
+     *
+     * \tparam Geometry The type of the geometry object.
+     * \tparam ScalarType The scalar type for the matrix elements.
+     *
+     * \return The second derivative of the membrane strain.
+     */
     template <typename Geometry, typename ScalarType>
     auto secondDerivative(const Dune::FieldVector<double, 2> &gpPos, const auto &dNAtGp, const Geometry &geo,
                           const auto &uFunction, const auto &localBasis, const Eigen::Vector3<ScalarType> &S, int I,

ikarus/finiteelements/mechanics/nonlinearelastic.hh

@@ -333,8 +333,7 @@ namespace Ikarus {
           const auto u = uFunction.evaluate(quadPos);
 
           // Value of the Neumann data at the current position
-          const auto neumannValue
-              = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
+          const auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
           energy -= neumannValue.dot(u) * curQuad.weight() * intElement;
         }
       }
@@ -395,8 +394,7 @@ namespace Ikarus {
             const auto udCi = u.evaluateDerivative(quadPos, wrt(coeff(i)));
 
             // Value of the Neumann data at the current position
-            const auto neumannValue
-                = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
+            const auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
             force.template segment<myDim>(myDim * i) -= udCi * neumannValue * curQuad.weight() * intElement;
           }
         }

tests/src/checkfebyautodiff.hh

@@ -1,25 +1,25 @@
 // SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
 // SPDX-License-Identifier: LGPL-3.0-or-later
 
-
 #pragma once
 #include <dune/common/test/testsuite.hh>
+#include <ikarus/utils/basis.hh>
+#include <ikarus/finiteelements/febases/autodifffe.hh>
+#include <Eigen/Core>
 
-template <template<typename ,typename > class FE, typename GridView,typename PreBasis, typename... ElementArgsType>
-auto checkFEByAutoDiff( const GridView& gridView, const PreBasis& pb,const ElementArgsType&... eleArgs) {
-
-  using namespace Dune::Functions::BasisFactory;
+template <template <typename...> class FE, typename GridView, typename PreBasis, typename... ElementArgsType>
+auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const ElementArgsType&... eleArgs) {
   auto basis       = Ikarus::makeBasis(gridView, pb);
   auto element     = gridView.template begin<0>();
   auto nDOF        = basis.flat().size();
   const double tol = 1e-10;
 
-  FE fe(basis,*element,eleArgs...);
+  FE fe(basis, *element, eleArgs...);
 
   const std::string feClassName = Dune::className(fe);
   Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation for gridDim = "
-            + feClassName);
-  using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe)>;
+                    + feClassName);
+  using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe),typename decltype(fe)::FERequirementType,false,true>;
   AutoDiffBasedFE feAutoDiff(fe);
 
   Eigen::VectorXd d;