[IN REVIEW] Add basic functionality for dynamic analysis

.clang-format

@@ -94,6 +94,9 @@ IncludeCategories:
   - Regex: '^(<|")autodiff/'
     Priority: 6
     SortPriority: 6
+  - Regex: '^(<|")Spectra/'
+    Priority: 7
+    SortPriority: 7
   - Regex:           '<([A-Za-z0-9.\/-_])+>'
     Priority:        2
     SortPriority: 2

CHANGELOG.md

@@ -56,6 +56,7 @@ SPDX-License-Identifier: LGPL-3.0-or-later
 - Add an About Ikarus page in the documentation ([#291](https://github.com/ikarus-project/ikarus/pull/291))
 - Add new class `Vtk::Writer`, which implements some convenience methods over the existing `dune-vtk` module ([#309](https://github.com/ikarus-project/ikarus/pull/309))
 - Add `VanishingStrain` material (useful for example for plane strain case), also refactor the constructor of `LinearElastic` to take any linear material law ([#317](https://github.com/ikarus-project/ikarus/pull/317))
+- Add a `ModalAnalysis` class to facilitate modal analysis and implement linear mass matrices for all elements. This also includes wrappers for solving the generalized eigenvalue problem with sparse and dense matrices with Eigen and Spectra
 
 ## Release v0.4 (Ganymede)
 

codespellignore

@@ -2,4 +2,6 @@
 # SPDX-License-Identifier: CC0-1.0
 te
 nd
-socio-economic
+socio-economic
+MATA
+MATB

ikarus/assembler/assemblermanipulatorbuildingblocks.hh

@@ -49,6 +49,10 @@
   void bind(F&& f) {
     sfs.emplace_back(std::forward<F>(f));
   }
+
+  /** \brief A helper function to remove all added scalar manipulator functions. */
+  void unbindAllScalarFunctions() { sfs.clear(); }
+
   std::vector<FunctionType> sfs;
 
 protected:
@@ -97,6 +101,10 @@
   void bind(F&& f) {
     vfs.emplace_back(std::forward<F>(f));
   }
+
+  /** \brief A helper function to remove all added vector manipulator functions. */
+  void unbindAllVectorFunctions() { vfs.clear(); }
+
   std::vector<FunctionType> vfs;
 
 protected:
@@ -143,6 +151,7 @@
   using MatrixType = typename Assembler::MatrixType;
   using Interface  = MatrixAssembler<WrappedAssembler, typename Assembler::FEContainer,
                                      typename Assembler::DirichletValuesType, typename Assembler::MatrixType>;
+
   friend Interface;
   using FunctionType =
       std::function<void(const Assembler&, const FERequirement&, MatrixAffordance, DBCOption, MatrixType&)>;
@@ -158,25 +167,44 @@
   void bind(F&& f) {
     mfs.emplace_back(std::forward<F>(f));
   }
+
+  /** \brief A helper function to remove all added matrix manipulator functions. */
+  void unbindAllMatrixFunctions() { mfs.clear(); }
+
   std::vector<FunctionType> mfs;
 
 protected:
   MatrixType& getRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance) {
-    MatrixType& mat = underlying().base_getRawMatrixImpl(feRequirements, affordance);
+    MatrixType& mat = [&]() -> MatrixType& {
+      if constexpr (Assembler::isSparse)
+        return underlying().base_getRawMatrixImpl(feRequirements, affordance, true);
+      else
+        return underlying().base_getRawMatrixImpl(feRequirements, affordance);
+    }();
     for (const auto& mf : mfs)
       mf(underlying().base(), feRequirements, affordance, DBCOption::Raw, mat);
     return mat;
   }
 
   MatrixType& getMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance) {
-    MatrixType& mat = underlying().base_getMatrixImpl(feRequirements, affordance);
+    MatrixType& mat = [&]() -> MatrixType& {
+      if constexpr (Assembler::isSparse)
+        return underlying().base_getMatrixImpl(feRequirements, affordance, true);
+      else
+        return underlying().base_getMatrixImpl(feRequirements, affordance);
+    }();
     for (const auto& mf : mfs)
       mf(underlying().base(), feRequirements, affordance, DBCOption::Full, mat);
     return mat;
   }
 
   MatrixType& getReducedMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance) {
-    MatrixType& mat = underlying().base_getReducedMatrixImpl(feRequirements, affordance);
+    MatrixType& mat = [&]() -> MatrixType& {
+      if constexpr (Assembler::isSparse)
+        return underlying().base_getReducedMatrixImpl(feRequirements, affordance, true);
+      else
+        return underlying().base_getReducedMatrixImpl(feRequirements, affordance);
+    }();
     for (const auto& mf : mfs)
       mf(underlying().base(), feRequirements, affordance, DBCOption::Reduced, mat);
     return mat;

ikarus/assembler/assemblermanipulatorfuser.hh

@@ -84,6 +84,7 @@ public:
   using WrappedAssembler::createFullVector;
 
   using WrappedAssembler::affordanceCollection;
+  using WrappedAssembler::basis;
   using WrappedAssembler::constraintsBelow;
   using WrappedAssembler::dBCOption;
   using WrappedAssembler::estimateOfConnectivity;
@@ -175,6 +176,7 @@ public:
   using WrappedAssembler::createFullVector;
 
   using WrappedAssembler::affordanceCollection;
+  using WrappedAssembler::basis;
   using WrappedAssembler::constraintsBelow;
   using WrappedAssembler::dBCOption;
   using WrappedAssembler::estimateOfConnectivity;
@@ -285,11 +287,13 @@ public:
   using WrappedAssembler::createFullVector;
 
   using WrappedAssembler::affordanceCollection;
+  using WrappedAssembler::basis;
   using WrappedAssembler::constraintsBelow;
   using WrappedAssembler::dBCOption;
   using WrappedAssembler::estimateOfConnectivity;
   using WrappedAssembler::gridView;
   using WrappedAssembler::isConstrained;
+  using WrappedAssembler::isSparse;
   using WrappedAssembler::reducedSize;
   using WrappedAssembler::requirement;
   using WrappedAssembler::size;

ikarus/assembler/interface.hh

@@ -94,6 +94,12 @@ public:
    */
   const auto& dirichletValues() const { return Dune::resolveRef(dirichletValues_); }
 
+  /**
+   * \brief Returns the basis object
+   * \return Reference to the basis
+   */
+  const auto& basis() const { return Dune::resolveRef(dirichletValues_.basis()); }
+
   /**
    * \brief Returns the gridView object.
    * \return Reference to the gridView object.

ikarus/assembler/simpleassemblers.hh

@@ -125,13 +125,18 @@ public:
   using typename Base::VectorType;
   using typename MatrixAssembler<SparseFlatAssembler, FEC, DV, Eigen::SparseMatrix<double>>::MatrixType;
 
+  static constexpr bool isSparse = true;
+
 private:
   void assembleRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance, MatrixType& assemblyMat);
 
 protected:
-  MatrixType& getRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
-  MatrixType& getMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
-  MatrixType& getReducedMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
+  MatrixType& getRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance,
+                               bool resetOccupationPattern = false);
+  MatrixType& getMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance,
+                            bool resetOccupationPattern = false);
+  MatrixType& getReducedMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance,
+                                   bool resetOccupationPattern = false);
 
 private:
   /** Calculates the non-zero entries in the full sparse matrix and passes them to the underlying Eigen sparse matrix.
@@ -205,6 +210,8 @@ public:
   using typename Base::VectorType;
   using typename MatrixAssembler<DenseFlatAssembler, FEC, DV, Eigen::MatrixXd>::MatrixType;
 
+  static constexpr bool isSparse = false;
+
 private:
   void assembleRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance, MatrixType& assemblyMat);
 

ikarus/assembler/simpleassemblers.inl

@@ -120,8 +120,8 @@
 
 template <typename B, typename FEC>
 typename SparseFlatAssembler<B, FEC>::MatrixType& SparseFlatAssembler<B, FEC>::getRawMatrixImpl(
-    const FERequirement& feRequirements, MatrixAffordance affordance) {
-  if (not sparsePreProcessorRaw_) {
+    const FERequirement& feRequirements, MatrixAffordance affordance, bool resetOccupationPattern) {
+  if (not sparsePreProcessorRaw_ or resetOccupationPattern) {
     preProcessSparseMatrix(spMatRaw_);
     sparsePreProcessorRaw_ = true;
   }
@@ -132,8 +132,8 @@
 
 template <typename B, typename FEC>
 typename SparseFlatAssembler<B, FEC>::MatrixType& SparseFlatAssembler<B, FEC>::getMatrixImpl(
-    const FERequirement& feRequirements, MatrixAffordance affordance) {
-  if (not sparsePreProcessor_) {
+    const FERequirement& feRequirements, MatrixAffordance affordance, bool resetOccupationPattern) {
+  if (not sparsePreProcessor_ or resetOccupationPattern) {
     preProcessSparseMatrix(spMat_);
     sparsePreProcessor_ = true;
   }
@@ -152,8 +152,8 @@
 
 template <typename B, typename FEC>
 typename SparseFlatAssembler<B, FEC>::MatrixType& SparseFlatAssembler<B, FEC>::getReducedMatrixImpl(
-    const FERequirement& feRequirements, MatrixAffordance affordance) {
-  if (not sparsePreProcessorReduced_) {
+    const FERequirement& feRequirements, MatrixAffordance affordance, bool resetOccupationPattern) {
+  if (not sparsePreProcessorReduced_ or resetOccupationPattern) {
     preProcessSparseMatrixReduced(spMatReduced_);
     sparsePreProcessorReduced_ = true;
   }

ikarus/finiteelements/ferequirements.hh

@@ -59,7 +59,7 @@ namespace Ikarus {
             geometricstiffness,
             stiffnessdiffBucklingVector,
             microMagneticHessian,
-            mass
+            linearMass
       );
 
   /**
@@ -170,7 +170,7 @@ inline constexpr VectorAffordance forces = VectorAffordance::forces;
 
 inline constexpr MatrixAffordance stiffness                   = MatrixAffordance::stiffness;
 inline constexpr MatrixAffordance stiffnessdiffBucklingVector = MatrixAffordance::stiffnessdiffBucklingVector;
-inline constexpr MatrixAffordance mass                        = MatrixAffordance::mass;
+inline constexpr MatrixAffordance mass                        = MatrixAffordance::linearMass;
 inline constexpr ScalarAffordance potentialEnergy             = ScalarAffordance::mechanicalPotentialEnergy;
 
 auto vectorAffordance(MatrixAffordance affordanceM) {
@@ -203,7 +203,10 @@ auto scalarAffordance(VectorAffordance affordanceV) {
 namespace AffordanceCollections {
   inline constexpr AffordanceCollection elastoStatics(ScalarAffordance::mechanicalPotentialEnergy,
                                                       VectorAffordance::forces, MatrixAffordance::stiffness);
-}
+  inline constexpr Ikarus::AffordanceCollection modalAnalysis(Ikarus::ScalarAffordance::noAffordance,
+                                                              Ikarus::VectorAffordance::noAffordance,
+                                                              Ikarus::MatrixAffordance::linearMass);
+} // namespace AffordanceCollections
 
 namespace Impl {
   template <typename T>

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -162,7 +162,7 @@ protected:
     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

@@ -31,6 +31,7 @@ struct KirchhoffLoveShellPre
 {
   YoungsModulusAndPoissonsRatio material;
   double thickness;
+  double density;
 
   template <typename PreFE, typename FE>
   using Skill = KirchhoffLoveShell<PreFE, FE>;
@@ -99,7 +100,8 @@ public:
    */
   KirchhoffLoveShell(const Pre& pre)
       : mat_{pre.material},
-        thickness_{pre.thickness} {}
+        thickness_{pre.thickness},
+        density_{pre.density} {}
 
 protected:
   /**
@@ -179,6 +181,7 @@ private:
   // DefaultMembraneStrain membraneStrain_;
   YoungsModulusAndPoissonsRatio mat_;
   double thickness_;
+  double density_;
 
   size_t numberOfNodes_{0};
   int order_{};
@@ -238,6 +241,10 @@ protected:
   void calculateMatrixImpl(
       const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<ST> K,
       const std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>& dx = std::nullopt) const {
+    if (affordance == MatrixAffordance::linearMass) {
+      calculateMassImpl<ST>(par, affordance, K);
+      return;
+    }
     if (affordance != MatrixAffordance::stiffness)
       DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));
     using namespace Dune::DerivativeDirections;
@@ -279,6 +286,30 @@ protected:
     K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
   }
 
+  template <typename ScalarType>
+  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();
+      auto& N               = localBasis_.evaluateFunction(gpIndex);
+
+      auto nopI = Eigen::Matrix<double, worldDim, worldDim>::Zero().eval();
+      auto nopJ = Eigen::Matrix<double, worldDim, worldDim>::Zero().eval();
+
+      for (size_t i = 0; i < numberOfNodes_; ++i) {
+        nopI.diagonal().setConstant(N[i]);
+        for (size_t j = 0; j < numberOfNodes_; ++j) {
+          nopJ.diagonal().setConstant(N[j]);
+          M.template block<worldDim, worldDim>(i * worldDim, j * worldDim) +=
+              nopI.transpose() * rhoT * nopJ * intElement;
+        }
+      }
+    }
+  }
+
   template <typename ST>
   void calculateVectorImpl(
       const Requirement& par, const VectorAffordance& affordance, typename Traits::template VectorType<ST> force,
@@ -428,13 +459,14 @@ protected:
 
 /**
  * \brief A struct containing information about the Youngs Modulus,
- * Poisson's ratio and the thickness for the Kirchhoff-Love shell element.
+ * Poisson's ratio, thickness and the density (defaults so 1.0) for the Kirchhoff-Love shell element.
  */
 struct KlArgs
 {
   double youngs_modulus;
   double nu;
   double thickness;
+  double density{1.0};
 };
 
 /**
@@ -443,7 +475,7 @@ struct KlArgs
  * \return A Kirchhoff-Love shell pre finite element.
  */
 auto kirchhoffLoveShell(const KlArgs& args) {
-  KirchhoffLoveShellPre pre({.emodul = args.youngs_modulus, .nu = args.nu}, args.thickness);
+  KirchhoffLoveShellPre pre({.emodul = args.youngs_modulus, .nu = args.nu}, args.thickness, args.density);
 
   return pre;
 }