Add subspaceEntriesToElementEntries functions for mapping localMatrix…

… to elementMatrix

ikarus/finiteelements/fehelper.hh

@@ -127,4 +127,46 @@ template <typename FiniteElement>
 void globalIndices(const FiniteElement& fe, std::vector<typename FiniteElement::LocalView::MultiIndex>& globalIndices) {
   globalIndicesFromLocalView(fe.localView(), globalIndices);
 }
+
+/**
+ * \brief Maps the entries from the localMatrix of the size \f$ n_s \f$ by \f$ n_s \f$ (where \f$ n_s \f$ is the number
+ of
+ * coefficiets of a subspace basis) to the local index in the element matrix.
+ *
+ * \tparam TreePaths type of the indices used to obtain the subspacebasis which corresponds to the localMatrix, can be a
+ * `index_constant`, an integer, or a `Dune::TypeTree::HybridTree` (one or more are accepted corresponding to the
+ * localMatrix)
+ * \param localMatrix the localMatrix containing the entries that are beeing mapped.
+ * \param elementMatrix the elementMatrix where the entries are beeing mapped to.
+ * \param globalBasis the global basis, from where we obtain the subspacebasis.
+ * \param tps the paths to the subspacebasis'.
+ */
+template <typename... TreePaths>
+void subspaceEntriesToElementEntries(const auto& localMatrix, auto& elementMatrix, const auto& globalBasis,
+                                     TreePaths&&... tps) {
+  auto treePathTuple = Dune::TupleVector(std::forward<TreePaths>(tps)...);
+  for (size_t i = 0; i < localMatrix.rows(); ++i) {
+    for (size_t j = 0; j < localMatrix.cols(); ++j) {
+      Dune::Hybrid::forEach(treePathTuple, [&](const auto& treePath) {
+        auto basis = globalBasis.child(treePath);
+        elementMatrix(basis.localIndex(i), basis.localIndex(j)) += localMatrix(i, j);
+      });
+    }
+  }
+}
+
+/**
+ * \brief Forwards the indices for powerbasis to obtain the children to subspaceEntriesToElementEntries()
+ *
+ * \tparam nEnd index from the last child + 1
+ * \tparam nStart optionnally define a starting index
+ * \details see Documentation for subspaceEntriesToElementEntries()
+ */
+template <int nEnd, int nStart = 0>
+void subspaceEntriesToElementEntries(const auto& localMatrix, auto& elementMatrix, const auto& globalBasis) {
+  Dune::Hybrid::forEach(
+      Dune::Hybrid::integralRange(Dune::index_constant<nStart>(), Dune::index_constant<nEnd>()),
+      [&](const auto i) { subspaceEntriesToElementEntries(localMatrix, elementMatrix, globalBasis, i); });
+}
+
 } // namespace Ikarus::FEHelper

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -159,6 +159,9 @@ 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::mass) {
+      return;
+    }
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
     easApplicabilityCheck();

ikarus/finiteelements/mechanics/linearelastic.hh

@@ -224,23 +224,13 @@ protected:
   void calculateMassImpl(const Requirement& par, const MatrixAffordance& affordance,
                          typename Traits::template MatrixType<> M) const {
     double rho = 1.0;
-    using namespace Dune;
-
     for (const auto& [gpIndex, gp] : localBasis_.viewOverIntegrationPoints()) {
       const auto intElement = geo_->integrationElement(gp.position()) * gp.weight();
       auto& N               = localBasis_.evaluateFunction(gpIndex);
       auto Mij              = (rho * N * N.transpose() * intElement).eval();
 
-      for (size_t i = 0; i < numberOfNodes_; ++i) {
-        for (size_t j = i; j < numberOfNodes_; ++j) {
-          Dune::Hybrid::forEach(Dune::Hybrid::integralRange(Dune::index_constant<myDim>()), [&](auto k) {
-            auto& basis = underlying().localView().tree().child(k);
-            M(basis.localIndex(i), basis.localIndex(j)) += Mij(i, j);
-          });
-        }
-      }
+      FEHelper::subspaceEntriesToElementEntries<myDim>(Mij, M, underlying().localView().tree());
     }
-    M.template triangularView<Eigen::StrictlyLower>() = M.transpose();
   }
 
   template <typename ScalarType>

ikarus/finiteelements/mechanics/nonlinearelastic.hh

@@ -242,6 +242,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::mass) {
+      calculateMassImpl<ScalarType>(par, affordance, K);
+      return;
+    }
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
     const auto uFunction = displacementFunction(par, dx);
@@ -264,6 +268,19 @@ protected:
     }
   }
 
+  template <typename ScalarType>
+  void calculateMassImpl(const Requirement& par, const MatrixAffordance& affordance,
+                         typename Traits::template MatrixType<> M) const {
+    double rho = 1.0;
+    for (const auto& [gpIndex, gp] : localBasis_.viewOverIntegrationPoints()) {
+      const auto intElement = geo_->integrationElement(gp.position()) * gp.weight();
+      auto& N               = localBasis_.evaluateFunction(gpIndex);
+      auto Mij              = (rho * N * N.transpose() * intElement).eval();
+
+      FEHelper::subspaceEntriesToElementEntries<myDim>(Mij, M, underlying().localView().tree());
+    }
+  }
+
   template <typename ScalarType>
   auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
                            const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx =

ikarus/utils/dynamics/dynamics.hh

@@ -7,6 +7,8 @@
  */
 
 #pragma once
+#include <dune/vtk/pvdwriter.hh>
+
 #include <Eigen/Core>
 #include <Eigen/SparseCore>
 
@@ -40,13 +42,15 @@ namespace Impl {
 
 template <typename ScalarType = double>
 inline auto lumpingSparse() {
-  return
-      [](const auto&, const auto&, auto, auto, Eigen::SparseMatrix<ScalarType>& mat) { return Impl::lumpingSparseImpl(mat); };
+  return [](const auto&, const auto&, auto, auto, Eigen::SparseMatrix<ScalarType>& mat) {
+    return Impl::lumpingSparseImpl(mat);
+  };
 }
 
 template <typename ScalarType = double>
 inline auto lumpingDense() {
-  return [](const auto&, const auto&, auto, auto, Eigen::MatrixX<ScalarType>& mat) { return Impl::lumpingDenseImpl(mat); };
+  return
+      [](const auto&, const auto&, auto, auto, Eigen::MatrixX<ScalarType>& mat) { return Impl::lumpingDenseImpl(mat); };
 }
 
 template <Concepts::FlatAssembler AS>
@@ -76,4 +80,24 @@ void writeEigenformsToVTK(const Eigensolver& solver, std::shared_ptr<Assembler>
   writer.write(filename);
 }
 
+template <Concepts::EigenValueSolver Eigensolver, Concepts::FlatAssembler Assembler>
+void writeEigenformsToPVD(const Eigensolver& solver, std::shared_ptr<Assembler> assembler, const std::string& filename,
+                          std::optional<Eigen::Index> nev_ = std::nullopt) {
+  auto nev          = nev_.value_or(solver.nev());
+  auto eigenvectors = solver.eigenvectors();
+  auto basis        = assembler->basis();
+
+  auto writer    = std::make_shared<decltype(Ikarus::Vtk::Writer(assembler))>(Ikarus::Vtk::Writer(assembler));
+  auto pvdWriter = Dune::Vtk::PvdWriter(writer);
+
+  Eigen::VectorXd evG(assembler->size());
+  writer->addInterpolation(evG, basis, "EF");
+
+  for (auto i : Dune::range(nev)) {
+    evG = assembler->createFullVector(eigenvectors.col(i));
+    pvdWriter.writeTimestep(i, filename, "data", false);
+  }
+  pvdWriter.write(filename);
+}
+
 } // namespace Ikarus::Dynamics

tests/src/testdynamics.cpp

@@ -41,7 +41,7 @@ static auto dynamicsTest() {
   const double Lx                         = 4.0;
   const double Ly                         = 4.0;
   Dune::FieldVector<double, 2> bbox       = {Lx, Ly};
-  std::array<int, 2> elementsPerDirection = {6, 6};
+  std::array<int, 2> elementsPerDirection = {4, 4};
   auto grid                               = std::make_shared<Grid>(bbox, elementsPerDirection);
 
   auto gridView = grid->leafGridView();
@@ -91,7 +91,7 @@ static auto dynamicsTest() {
   auto mat1 = assMLumped->matrix();
   auto mat2 = assMLumped->matrix();
 
-  int nev = 10; // number of requested eigenvalues
+  int nev = 50; // number of requested eigenvalues
   using Ikarus::Dynamics::EigenSolverTypeTag;
 
   auto solver  = Ikarus::Dynamics::PartialSparseGeneralSymEigenSolver(assK, assM, nev);
@@ -110,6 +110,8 @@ static auto dynamicsTest() {
   std::cout << "Angular frequencies found (n=" << nev << "):\n" << eigenvalues2.array().sqrt() << std::endl;
 
   Ikarus::Dynamics::writeEigenformsToVTK(solver2, assM, "eigenforms");
+  Ikarus::Dynamics::writeEigenformsToPVD(solver, assM, "eigenforms_full");
+  // Ikarus::Dynamics::writeEigenformsToPVD(solver2, assM, "eigenforms_lumped");
 
   auto solver3 = Ikarus::Dynamics::makeGeneralSymEigenSolver<EigenSolverTypeTag::Eigen>(assKD, assMD);
   solver3.compute();
@@ -127,7 +129,8 @@ static auto dynamicsTest() {
   t.check(isApproxSame(eigenvalues5, eigenvalues4, 1e-10));
   t.check(isApproxSame(eigenvalues4.head(nev).eval(), eigenvalues, 1e-10));
 
-  Ikarus::Dynamics::writeEigenformsToVTK(solver4, assM, "eigenforms_all");
+  Ikarus::Dynamics::writeEigenformsToVTK(solver4, assM, "eigenforms_vtk");
+  Ikarus::Dynamics::writeEigenformsToPVD(solver4, assM, "eigenforms_pvd");
 
   static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver)>);
   static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver2)>);
@@ -138,10 +141,14 @@ static auto dynamicsTest() {
   return t;
 }
 
+
+
+
 int main(int argc, char** argv) {
   Ikarus::init(argc, argv);
   TestSuite t;
 
+
   t.subTest(dynamicsTest());
   return t.exit();
 }

tests/src/testvtkwriter.cpp

@@ -50,11 +50,11 @@ auto vtkWriterTest() {
 
   DummyProblem<Grid> testCase{};
 
-  auto& gridView       = testCase.gridView();
+  auto& gridView        = testCase.gridView();
   auto& sparseAssembler = testCase.sparseAssembler();
-  auto& req            = testCase.requirement();
-  auto& basis          = testCase.basis();
-  auto D_Glob         = req.globalSolution();
+  auto& req             = testCase.requirement();
+  auto& basis           = testCase.basis();
+  auto D_Glob           = req.globalSolution();
 
   // Tests
   Dune::Vtk::DiscontinuousDataCollector dc{gridView};