fix test

ikarus-project · Nov 27, 2024 · 9f92ea6 · 9f92ea6
1 parent 4e546d3
commit 9f92ea6
Show file tree

Hide file tree

Showing 4 changed files with 46 additions and 62 deletions.
diff --git a/ikarus/solver/eigenvaluesolver/generaleigensolver.hh b/ikarus/solver/eigenvaluesolver/generaleigensolver.hh
@@ -116,7 +116,6 @@ struct GeneralSymEigenSolver<EigenSolverTypeTag::Spectra, MT>
     return eigenvectors_;
   }
 
-
   Eigen::Index nev() const { return nev_; }
 
 private:
@@ -148,13 +147,11 @@ struct GeneralSymEigenSolver<EigenSolverTypeTag::Eigen, MT>
   template <typename MATA, typename MATB>
   requires(Concepts::EigenMatrix<std::remove_cvref_t<MATA>>)
   GeneralSymEigenSolver(MATA&& A, MATB&& B)
-      : matA_(A),
-        matB_(B),
+      : matA_(std::forward<MATA>(A)),
+        matB_(std::forward<MATB>(B)),
         solver_(A.size()) {
     if (A.cols() != B.cols())
       DUNE_THROW(Dune::IOError, "GeneralSymEigenSolver: The passed matrices should have the same size");
-    std::cout << A << std::endl;
-    std::cout << B << std::endl;
   }
 
   template <Concepts::FlatAssembler AssemblerA, Concepts::FlatAssembler AssemblerB>
@@ -199,7 +196,6 @@ struct GeneralSymEigenSolver<EigenSolverTypeTag::Eigen, MT>
     return solver_.eigenvectors();
   }
 
-
 private:
   MatrixType matA_;
   MatrixType matB_;
@@ -294,7 +290,6 @@ struct PartialGeneralSymEigenSolver
     return solver_.eigenvectors(_nev.value_or(nev_));
   }
 
-
   Eigen::Index nev() const { return nev_; }
 
 private:

diff --git a/ikarus/utils/dynamics/modalanalysis.hh b/ikarus/utils/dynamics/modalanalysis.hh
@@ -130,12 +130,8 @@ struct ModalAnalysis
 
   auto nev() const { return solver_->nev(); }
 
-  auto& stiffnessAssembler() const {
-    return stiffAssembler_;
-  }
-  auto& massAssembler() const {
-    return lumpedMassAssembler_;
-  }
+  auto& stiffnessAssembler() const { return stiffAssembler_; }
+  auto& massAssembler() const { return lumpedMassAssembler_; }
 
 private:
   FEContainer fes_;

diff --git a/tests/src/testdynamics.cpp b/tests/src/testdynamics.cpp
@@ -72,7 +72,7 @@ static auto dynamicsTest() {
   mA.registerLumpingScheme<Ikarus::Dynamics::LumpingSchemes::RowSumLumping>();
   mA.compute();
 
-  auto frequenciesLumped = mA.angularFrequencies();
+  auto frequenciesLumped  = mA.angularFrequencies();
   auto frequenciesLumped2 = mA.frequencies(Ikarus::Dynamics::ModalAnalysisResultType::angularFrequency);
 
   t.check(isApproxSame(frequenciesLumped, frequenciesLumped2, 1e-16));

diff --git a/tests/src/testeigenvaluesolver.cpp b/tests/src/testeigenvaluesolver.cpp
@@ -34,14 +34,10 @@
 using Dune::TestSuite;
 
 template <Ikarus::Concepts::EigenValueSolver SOL1, Ikarus::Concepts::EigenValueSolver SOL2>
-auto testSolversAgainstEachOther(const SOL1& solver1, const SOL2& solver2,
-                                 std::optional<Eigen::Index> nev_ = std::nullopt) {
-  TestSuite t("Test " + Dune::className<SOL1>() + ", " + Dune::className<SOL2>());
-  auto eigenvalues1  = solver1.eigenvalues();
-  auto eigenvectors1 = solver1.eigenvectors();
-
-  auto eigenvalues2  = solver2.eigenvalues();
-  auto eigenvectors2 = solver2.eigenvectors();
+auto testEigenValues(const SOL1& solver1, const SOL2& solver2, std::optional<Eigen::Index> nev_ = std::nullopt) {
+  TestSuite t("Test Eigenvalues" + Dune::className<SOL1>() + ", " + Dune::className<SOL2>());
+  auto eigenvalues1 = solver1.eigenvalues();
+  auto eigenvalues2 = solver2.eigenvalues();
 
   if (not nev_.has_value()) {
     t.check(isApproxSame(eigenvalues1, eigenvalues2, 1e-10)) << testLocation();
@@ -50,44 +46,38 @@ auto testSolversAgainstEachOther(const SOL1& solver1, const SOL2& solver2,
         << testLocation() << "\n"
         << eigenvalues1.head(*nev_).eval() << "\n\n"
         << eigenvalues2.head(*nev_).eval();
-
-    // t.check(isApproxSame(eigenvectors1.leftCols(*nev_).eval(), eigenvectors2.leftCols(*nev_).eval(), 1e-10))
-    //     // << testLocation() << "\n"
-    //     << "\n"
-    //     << eigenvectors1.leftCols(*nev_).eval() << "\n\n"
-    //     << eigenvectors2.leftCols(*nev_).eval();
   }
+
+  return t;
 }
 
+// as1 and as2 correspond to solver1, which is tested against solver2
 template <Ikarus::Concepts::EigenValueSolver SOL1, Ikarus::Concepts::EigenValueSolver SOL2,
           Ikarus::Concepts::FlatAssembler AS1, Ikarus::Concepts::FlatAssembler AS2>
-auto testSparseAgainstDense(const SOL1& solver1, const SOL2& solver2, std::shared_ptr<AS1> as1,
-                            std::shared_ptr<AS2> as2) {
-  // We assume as1 and as2 are sparse assemblers
-  Eigen::MatrixXd K = as1->matrix().toDense();
-  Eigen::MatrixXd M = as2->matrix().toDense();
+auto testEigenVectors(const SOL1& solver1, const SOL2& solver2, std::shared_ptr<AS1> as1, std::shared_ptr<AS2> as2) {
+  TestSuite t("Test Eigenvalues " + Dune::className<SOL1>() + ", " + Dune::className<SOL2>());
 
-  auto MK = M.inverse() * K;
+  Eigen::MatrixXd K = as1->matrix();
+  Eigen::MatrixXd M = as2->matrix();
 
   Eigen::MatrixXd evecs = solver1.eigenvectors();
   Eigen::VectorXd evals = solver2.eigenvalues();
 
-  Eigen::MatrixXd T = evecs.transpose() * MK * evecs;
+  Eigen::VectorXd T = (evecs.transpose() * K * evecs).diagonal();
 
-  std::cout << T.diagonal() << std::endl;
-  std::cout << evals << std::endl;
+  t.check(isApproxSame(T, evals, 1e-10)) << testLocation() << "\n" << T << "\n\n" << evals;
 
-  // return t;
+  return t;
 }
 
 auto testRealWorldProblem() {
-  TestSuite t("EigenvalueTest");
+  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 = {3, 3};
+  std::array<int, 2> elementsPerDirection = {4, 4};
   auto grid                               = std::make_shared<Grid>(bbox, elementsPerDirection);
 
   auto gridView = grid->leafGridView();
@@ -126,36 +116,39 @@ auto testRealWorldProblem() {
   assMD->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced);
   assKD->bind(req, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced);
 
-  int nev = 4; // number of requested eigenvalues
+  int nev = 10; // number of requested eigenvalues
   using Ikarus::EigenSolverTypeTag;
 
-  // auto solver = Ikarus::PartialGeneralSymEigenSolver(assK, assM, nev);
-  // t.checkThrow([&]() { solver.eigenvalues(); }) << testLocation();
-  // bool success = solver.compute();
-  // t.check(success) << testLocation();
+  auto partialSolver = Ikarus::PartialGeneralSymEigenSolver(assK, assM, nev);
+  t.checkThrow([&]() { partialSolver.eigenvalues(); }) << testLocation();
+  bool success = partialSolver.compute();
+  t.check(success) << testLocation();
 
-  // auto solver3 = Ikarus::makeGeneralSymEigenSolver<EigenSolverTypeTag::Eigen>(assKD, assMD);
-  // t.check(solver3.compute()) << testLocation();
+  auto solver1 = Ikarus::makeGeneralSymEigenSolver<EigenSolverTypeTag::Eigen>(assKD, assMD);
+  t.check(solver1.compute()) << testLocation();
 
-  auto solver4 = Ikarus::makeGeneralSymEigenSolver<EigenSolverTypeTag::Spectra>(assK, assM);
-  t.check(solver4.compute()) << testLocation();
+  auto solver2 = Ikarus::makeGeneralSymEigenSolver<EigenSolverTypeTag::Spectra>(assK, assM);
+  t.check(solver2.compute()) << testLocation();
 
-  auto solver5 = Ikarus::makeGeneralSymEigenSolver<EigenSolverTypeTag::Spectra>(assKD, assMD);
-  t.check(solver5.compute()) << testLocation();
+  auto solver3 = Ikarus::makeGeneralSymEigenSolver<EigenSolverTypeTag::Spectra>(assKD, assMD);
+  t.check(solver3.compute()) << testLocation();
 
-  testSparseAgainstDense(solver4, solver4, assK, assM);
+  t.subTest(testEigenVectors(solver2, solver3, assK, assM));
+  t.subTest(testEigenVectors(solver3, solver2, assKD, assMD));
+  t.subTest(testEigenVectors(solver1, solver2, assKD, assMD));
+  t.subTest(testEigenVectors(solver2, solver1, assK, assM));
 
-  // t.subTest(testSolversAgainstEachOther(solver, solver3, nev));
-  // t.subTest(testSolversAgainstEachOther(solver, solver4, nev));
-  // t.subTest(testSolversAgainstEachOther(solver, solver5, nev));
-  // t.subTest(testSolversAgainstEachOther(solver3, solver4));
-  // t.subTest(testSolversAgainstEachOther(solver3, solver5));
-  // t.subTest(testSolversAgainstEachOther(solver4, solver5));
+  t.subTest(testEigenValues(partialSolver, solver1, nev));
+  t.subTest(testEigenValues(partialSolver, solver2, nev));
+  t.subTest(testEigenValues(partialSolver, solver3, nev));
+  t.subTest(testEigenValues(solver1, solver2));
+  t.subTest(testEigenValues(solver1, solver3));
+  t.subTest(testEigenValues(solver2, solver3));
 
-  // static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver)>);
-  // static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver3)>);
-  // static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver4)>);
-  // static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver5)>);
+  static_assert(Ikarus::Concepts::EigenValueSolver<decltype(partialSolver)>);
+  static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver1)>);
+  static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver2)>);
+  static_assert(Ikarus::Concepts::EigenValueSolver<decltype(solver3)>);
 
   return t;
 }