review part 1

ikarus/finiteelements/fehelper.hh

@@ -127,5 +127,4 @@ template <typename FiniteElement>
 void globalIndices(const FiniteElement& fe, std::vector<typename FiniteElement::LocalView::MultiIndex>& globalIndices) {
   globalIndicesFromLocalView(fe.localView(), globalIndices);
 }
-
 } // namespace Ikarus::FEHelper

ikarus/finiteelements/ferequirements.hh

@@ -203,7 +203,7 @@ auto scalarAffordance(VectorAffordance affordanceV) {
 namespace AffordanceCollections {
   inline constexpr AffordanceCollection elastoStatics(ScalarAffordance::mechanicalPotentialEnergy,
                                                       VectorAffordance::forces, MatrixAffordance::stiffness);
-  constexpr Ikarus::AffordanceCollection dynamics(Ikarus::ScalarAffordance::noAffordance,
+  constexpr Ikarus::AffordanceCollection modalAnalysis(Ikarus::ScalarAffordance::noAffordance,
                                                   Ikarus::VectorAffordance::noAffordance,
                                                   Ikarus::MatrixAffordance::linearMass);
 } // namespace AffordanceCollections

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -159,13 +159,10 @@ protected:
   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::linearMass) {
-      return;
-    }
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
     easApplicabilityCheck();
-    if (isDisplacementBased())
+    if (isDisplacementBased() or affordance != MatrixAffordance::stiffness)
       return;
 
     decltype(auto) LMat = [this]() -> decltype(auto) {

ikarus/finiteelements/mechanics/kirchhoffloveshell.hh

@@ -290,6 +290,7 @@ protected:
   void calculateMassImpl(const Requirement& par, const MatrixAffordance& affordance,
                          typename Traits::template MatrixType<> M) const {
     const auto geo = underlying().localView().element().geometry();
+    const auto rhoT = thickness_ * density_;
 
     for (const auto& [gpIndex, gp] : localBasis_.viewOverIntegrationPoints()) {
       const auto intElement = geo.integrationElement(gp.position()) * gp.weight();
@@ -303,7 +304,7 @@ protected:
         for (size_t j = 0; j < numberOfNodes_; ++j) {
           nopJ.diagonal().setConstant(N[j]);
           M.template block<worldDim, worldDim>(i * worldDim, j * worldDim) +=
-              nopI.transpose() * density_ * nopJ * intElement;
+              nopI.transpose() * rhoT * nopJ * intElement;
         }
       }
     }

ikarus/finiteelements/mechanics/truss.hh

@@ -74,6 +74,7 @@ public:
   struct KinematicVariables
   {
     double L;                  ///< Length of the reference geometry
+    Eigen::VectorXd A1;        ///< Length of the reference geometry as vector
     ST l;                      ///< Length of the deformed geometry
     ST Elin;                   ///< Linear strain
     ST Egl;                    ///< Green-Lagrange strain
@@ -136,8 +137,9 @@ public:
     const double Lsquared = A1.squaredNorm();
     const ST lsquared     = (x2 - x1).squaredNorm();
 
-    kin.L = sqrt(Lsquared);
-    kin.l = sqrt(lsquared);
+    kin.L  = sqrt(Lsquared);
+    kin.l  = sqrt(lsquared);
+    kin.A1 = A1;
 
     // Linear strain
     kin.Elin = (kin.l - kin.L) / kin.L;
@@ -164,12 +166,8 @@ public:
    *
    * \return std::pair of length, and T
    */
-  auto computeLengthAndTransformationMatrix() const {
-    auto& ele           = underlying().localView().element();
-    const auto X1       = Dune::toEigen(ele.geometry().corner(0));
-    const auto X2       = Dune::toEigen(ele.geometry().corner(1));
-    const auto A1       = X2 - X1;
-    const auto L        = A1.norm();
+  auto computeLengthAndTransformationMatrix(const Requirement& par) const {
+    const auto [L, A1, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par);
     const auto A1normed = A1 / L;
 
     auto T = Eigen::Matrix<typename std::remove_cvref_t<decltype(L)>, 2, worldDim * 2>::Zero().eval();
@@ -194,8 +192,8 @@ public:
   requires(supportsResultType<RT>())
   auto calculateAtImpl(const Requirement& req, [[maybe_unused]] const Dune::FieldVector<double, Traits::mydim>& local,
                        Dune::PriorityTag<0>) const {
-    using RTWrapper                            = ResultWrapper<RT<double, myDim, myDim>, ResultShape::Vector>;
-    const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(req);
+    using RTWrapper                               = ResultWrapper<RT<double, myDim, myDim>, ResultShape::Vector>;
+    const auto [L, _, l, Elin, Egl, dEdu, ddEddu] = computeStrain(req);
     if constexpr (isSameResultType<RT, ResultTypes::cauchyAxialForce>) {
       auto N = Eigen::Vector<double, 1>{E_ * A_ * Egl * l / L}; // Axial force in deformed configuration
       return RTWrapper{N};
@@ -225,14 +223,14 @@ protected:
       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) {
-      const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx);
+      const auto [L, _, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx);
       K += E_ * A_ * L * (dEdu * dEdu.transpose() + ddEddu * Egl);
     } else if (affordance == MatrixAffordance::linearMass) {
       Eigen::Matrix<ScalarType, 2, 2> mLoc{
           {2, 1},
           {1, 2}
       };
-      auto [l, T] = computeLengthAndTransformationMatrix();
+      auto [l, T] = computeLengthAndTransformationMatrix(par);
       mLoc *= (density_ * A_ * l) / 6.0;
       K += T.transpose() * mLoc * T;
     } else
@@ -243,15 +241,15 @@ protected:
   auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
                            const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx =
                                std::nullopt) const -> ScalarType {
-    const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx);
+    const auto [L, l, _, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx);
     return 0.5 * E_ * A_ * L * Egl * Egl;
   }
 
   template <typename ScalarType>
   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 {
-    const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx);
+    const auto [L, l, _, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx);
     force += E_ * A_ * Egl * L * dEdu;
   }
 };

ikarus/python/test/kltest.py

@@ -69,7 +69,7 @@ def vL(x, lambdaVal):
     vLoad = iks.finite_elements.volumeLoad3D(vL)
 
     klShell = iks.finite_elements.kirchhoffLoveShell(
-        youngs_modulus=1000, nu=0.0, thickness=thickness, density=density
+        youngs_modulus=1000, nu=0.0, thickness=thickness
     )
 
     fes = []

ikarus/python/test/modalanalysistest.py

@@ -32,7 +32,7 @@ def setUp(self):
 
         self.fes = []
         linMat = iks.materials.LinearElasticity(E=1000, nu=0.2).asPlaneStrain()
-        linElastic = finite_elements.linearElastic(linMat)
+        linElastic = finite_elements.linearElastic(linMat, density=100)
         for e in self.grid.elements:
             self.fes.append(finite_elements.makeFE(basis, linElastic))
             self.fes[-1].bind(e)

ikarus/solver/eigenvaluesolver/generalizedeigensolver.hh

@@ -3,7 +3,7 @@
 
 /**
  * \file generalizedeigensolver.hh
- * \brief Helper for dune-functions
+ * \brief Implementation of wrapper classes for solving a generalized eigenvalue problem
  */
 
 #pragma once
@@ -26,7 +26,6 @@ namespace Ikarus {
 
 /**
  * \brief Construct a new make enum object
- *
  */
 MAKE_ENUM(EigenValueSolverType, Spectra, Eigen);
 
@@ -191,7 +190,7 @@ struct GeneralizedSymEigenSolver<EigenValueSolverType::Eigen, MT>
   }
 
   template <Concepts::FlatAssembler AssemblerA, Concepts::FlatAssembler AssemblerB>
-  GeneralizedSymEigenSolver(const std::shared_ptr<AssemblerA> assemblerA, const std::shared_ptr<AssemblerB>& assemblerB)
+  GeneralizedSymEigenSolver(const std::shared_ptr<AssemblerA>& assemblerA, const std::shared_ptr<AssemblerB>& assemblerB)
       : GeneralizedSymEigenSolver(assemblerA->matrix(), assemblerB->matrix()) {}
 
   /**
@@ -305,7 +304,7 @@ struct PartialGeneralizedSymEigenSolver
                Spectra::SortRule sortRule = Spectra::SortRule::SmallestAlge, ScalarType tolerance = 1e-10,
                Eigen::Index maxit = 1000) {
     solver_.init();
-    solver_.compute(Spectra::SortRule::SmallestAlge, 1000, 1e-10, Spectra::SortRule::SmallestAlge);
+    solver_.compute(selection, tolerance, maxit, sortRule);
 
     computed_ = solver_.info() == Spectra::CompInfo::Successful;
     return computed_;

ikarus/utils/modalanalysis/lumpingschemes.hh

@@ -15,7 +15,7 @@
 namespace Ikarus::Dynamics::LumpingSchemes {
 
 /**
- * \brief Implements a row-sum lumping schemes for assemblermanipulator. Works with dense and sprase matrices.
+ * \brief Implements the operator() for a row-sum lumping scheme that can be passed to the assemblermanipulator to lump the already assembled global mass matrix.
  */
 struct RowSumLumping
 {

ikarus/utils/modalanalysis/modalanalysis.hh

@@ -66,7 +66,7 @@ struct ModalAnalysis
       req_.insertGlobalSolution(d_);
 
     stiffAssembler_->bind(req_, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced);
-    massAssembler_->bind(req_, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced);
+    massAssembler_->bind(req_, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced);
     lumpedMassAssembler_ = makeAssemblerManipulator(*massAssembler_);
   }
 
@@ -76,7 +76,7 @@ struct ModalAnalysis
    * It resets already bound matrix manipulation functions.
    *
    * \tparam LumpingScheme The type of the lumping scheme, for example one found at \file
-   * ikarus/utils/modalanalysis/lumpingschemes.hh.s
+   * ikarus/utils/modalanalysis/lumpingschemes.hh.
    * \param ls the instantiated LumpingScheme object (pass either by value or by template definition).
    */
   template <typename LumpingScheme>

tests/src/dummyproblem.hh

@@ -43,10 +43,8 @@ struct DummyProblem
 
   // YASPGrid needs an int, structuresgridfactory an unsigned int
   explicit DummyProblem(
-      const std::array<std::conditional_t<useYASP, int, unsigned int>, 2>& elementsPerDirection = {10, 10})
+      const std::array<std::conditional_t<useYASP, int, unsigned int>, 2>& elementsPerDirection = {10, 10}, double Lx = 4.0, double Ly = 4.0)
       : grid_([&]() {
-          constexpr double Lx                     = 4.0;
-          constexpr double Ly                     = 4.0;
           const Dune::FieldVector<double, 2> bbox = {Lx, Ly};
 
           if constexpr (not useYASP)

tests/src/testeigenvaluesolver.cpp

@@ -2,9 +2,8 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 #include <config.h>
 
-#include "tests/src/testhelpers.hh"
-
-#include <vector>
+#include "dummyproblem.hh"
+#include "testhelpers.hh"
 
 #include <dune/common/float_cmp.hh>
 #include <dune/common/test/testsuite.hh>
@@ -63,50 +62,26 @@ auto testEigenVectors(const SOL1& solver1, const SOL2& solver2, std::shared_ptr<
   return t;
 }
 
-auto testRealWorldProblem() {
-  TestSuite t("RealWorldProblem");
-  using Grid = Dune::YaspGrid<2>;
-
-  const double Lx                         = 4.0;
-  const double Ly                         = 4.0;
-  Dune::FieldVector<double, 2> bbox       = {Lx, Ly};
-  std::array<int, 2> elementsPerDirection = {4, 4};
-  auto grid                               = std::make_shared<Grid>(bbox, elementsPerDirection);
-
-  auto gridView = grid->leafGridView();
-
-  using namespace Dune::Functions::BasisFactory;
-  auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>()));
-
-  auto vL      = []([[maybe_unused]] auto& globalCoord, auto& lamb) { return Eigen::Vector2d{0, -1}; };
-  auto linMat  = Ikarus::LinearElasticity(Ikarus::toLamesFirstParameterAndShearModulus({.emodul = 100, .nu = 0.2}));
-  auto skills_ = Ikarus::skills(Ikarus::linearElastic(Ikarus::planeStress(linMat)), Ikarus::volumeLoad<2>(vL));
-
-  using FEType = decltype(Ikarus::makeFE(basis, skills_));
-  std::vector<FEType> fes;
-  for (auto&& element : elements(gridView)) {
-    fes.emplace_back(Ikarus::makeFE(basis, skills_));
-    fes.back().bind(element);
-  }
-
-  auto req = FEType::Requirement();
-  typename FEType::Requirement::SolutionVectorType d;
-  d.setZero(basis.flat().size());
-  req.insertGlobalSolution(d);
+auto testEigenSolvers() {
+  TestSuite t("EigenValueSolver Test");
+  using Grid     = Dune::YaspGrid<2>;
+  using GridView = Grid::LeafGridView;
 
-  auto dirichletValues = Ikarus::DirichletValues(basis.flat());
-  dirichletValues.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; });
+  DummyProblem<Grid, true> testCase({5, 5});
+  auto& req             = testCase.requirement();
+  auto& fes             = testCase.finiteElements();
+  auto& dirichletValues = testCase.dirichletValues();
 
   auto assM = Ikarus::makeSparseFlatAssembler(fes, dirichletValues);
   auto assK = Ikarus::makeSparseFlatAssembler(fes, dirichletValues);
 
   auto assMD = Ikarus::makeDenseFlatAssembler(fes, dirichletValues);
   auto assKD = Ikarus::makeDenseFlatAssembler(fes, dirichletValues);
 
-  assM->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced);
+  assM->bind(req, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced);
   assK->bind(req, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced);
 
-  assMD->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced);
+  assMD->bind(req, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced);
   assKD->bind(req, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced);
 
   int nev = 10; // number of requested eigenvalues
@@ -192,6 +167,6 @@ int main(int argc, char** argv) {
   Ikarus::init(argc, argv);
   TestSuite t;
 
-  t.subTest(testRealWorldProblem());
+  t.subTest(testEigenSolvers());
   return t.exit();
 }

tests/src/testmodalanalysis.cpp

@@ -2,10 +2,9 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 #include <config.h>
 
+#include "tests/src/dummyproblem.hh"
 #include "tests/src/testhelpers.hh"
 
-#include <vector>
-
 #include <dune/common/float_cmp.hh>
 #include <dune/common/test/testsuite.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
@@ -25,34 +24,14 @@
 using Dune::TestSuite;
 
 static auto modalAnalysisTest() {
-  TestSuite t("ModalAnalysis test");
+  TestSuite t("ModalAnalysis Test");
   using Grid = Dune::YaspGrid<2>;
 
-  const double Lx                         = 4.0;
-  const double Ly                         = 4.0;
-  Dune::FieldVector<double, 2> bbox       = {Lx, Ly};
-  std::array<int, 2> elementsPerDirection = {5, 5};
-  auto grid                               = std::make_shared<Grid>(bbox, elementsPerDirection);
-
-  auto gridView = grid->leafGridView();
-
-  using namespace Dune::Functions::BasisFactory;
-  auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>()));
-
-  auto linMat  = Ikarus::LinearElasticity(Ikarus::toLamesFirstParameterAndShearModulus({.emodul = 100, .nu = 0.2}));
-  auto skills_ = Ikarus::skills(Ikarus::linearElastic(Ikarus::planeStress(linMat)));
-
-  using FEType = decltype(Ikarus::makeFE(basis, skills_));
-  std::vector<FEType> fes;
-  for (auto&& element : elements(gridView)) {
-    fes.emplace_back(Ikarus::makeFE(basis, skills_));
-    fes.back().bind(element);
-  }
-
-  auto dirichletValues = Ikarus::DirichletValues(basis.flat());
-  dirichletValues.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; });
-
-  auto mA = Ikarus::Dynamics::ModalAnalysis(std::move(fes), dirichletValues);
+  DummyProblem<Grid, true> testCase({5, 5});
+  auto& req             = testCase.requirement();
+  auto& fes             = testCase.finiteElements();
+  auto& dirichletValues = testCase.dirichletValues();
+  auto mA               = Ikarus::Dynamics::ModalAnalysis(std::move(fes), dirichletValues);
   mA.compute();
 
   auto frequencies = mA.angularFrequencies();
@@ -70,13 +49,8 @@ static auto modalAnalysisTest() {
 
   mA.writeEigenModes("eigenforms", 20);
 
-  auto req = FEType::Requirement();
-  typename FEType::Requirement::SolutionVectorType d;
-  d.setZero(basis.flat().size());
-  req.insertGlobalSolution(d);
-
   auto massAssembler = Ikarus::makeSparseFlatAssembler(fes, dirichletValues);
-  massAssembler->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced);
+  massAssembler->bind(req, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced);
   auto lumpedAssembler = Ikarus::Dynamics::makeLumpedFlatAssembler(massAssembler);
 
   auto M       = massAssembler->matrix();