From 9f92ea6f96faeaee63ca9d9fdde75889071097c9 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Wed, 27 Nov 2024 19:34:04 +0100 Subject: [PATCH] fix test --- .../eigenvaluesolver/generaleigensolver.hh | 9 +- ikarus/utils/dynamics/modalanalysis.hh | 8 +- tests/src/testdynamics.cpp | 2 +- tests/src/testeigenvaluesolver.cpp | 89 +++++++++---------- 4 files changed, 46 insertions(+), 62 deletions(-) diff --git a/ikarus/solver/eigenvaluesolver/generaleigensolver.hh b/ikarus/solver/eigenvaluesolver/generaleigensolver.hh index 347805f90..2f1d2dff2 100644 --- a/ikarus/solver/eigenvaluesolver/generaleigensolver.hh +++ b/ikarus/solver/eigenvaluesolver/generaleigensolver.hh @@ -116,7 +116,6 @@ struct GeneralSymEigenSolver return eigenvectors_; } - Eigen::Index nev() const { return nev_; } private: @@ -148,13 +147,11 @@ struct GeneralSymEigenSolver template requires(Concepts::EigenMatrix>) GeneralSymEigenSolver(MATA&& A, MATB&& B) - : matA_(A), - matB_(B), + : matA_(std::forward(A)), + matB_(std::forward(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 @@ -199,7 +196,6 @@ struct GeneralSymEigenSolver 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 index addf160d0..4e4a694b5 100644 --- 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 index c9368d837..c68a04afa 100644 --- a/tests/src/testdynamics.cpp +++ b/tests/src/testdynamics.cpp @@ -72,7 +72,7 @@ static auto dynamicsTest() { mA.registerLumpingScheme(); 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 index 91d8ccfca..c14a7a86b 100644 --- a/tests/src/testeigenvaluesolver.cpp +++ b/tests/src/testeigenvaluesolver.cpp @@ -34,14 +34,10 @@ using Dune::TestSuite; template -auto testSolversAgainstEachOther(const SOL1& solver1, const SOL2& solver2, - std::optional nev_ = std::nullopt) { - TestSuite t("Test " + Dune::className() + ", " + Dune::className()); - 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 nev_ = std::nullopt) { + TestSuite t("Test Eigenvalues" + Dune::className() + ", " + Dune::className()); + 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 -auto testSparseAgainstDense(const SOL1& solver1, const SOL2& solver2, std::shared_ptr as1, - std::shared_ptr 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, std::shared_ptr as2) { + TestSuite t("Test Eigenvalues " + Dune::className() + ", " + Dune::className()); - 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 bbox = {Lx, Ly}; - std::array elementsPerDirection = {3, 3}; + std::array elementsPerDirection = {4, 4}; auto grid = std::make_shared(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(assKD, assMD); - // t.check(solver3.compute()) << testLocation(); + auto solver1 = Ikarus::makeGeneralSymEigenSolver(assKD, assMD); + t.check(solver1.compute()) << testLocation(); - auto solver4 = Ikarus::makeGeneralSymEigenSolver(assK, assM); - t.check(solver4.compute()) << testLocation(); + auto solver2 = Ikarus::makeGeneralSymEigenSolver(assK, assM); + t.check(solver2.compute()) << testLocation(); - auto solver5 = Ikarus::makeGeneralSymEigenSolver(assKD, assMD); - t.check(solver5.compute()) << testLocation(); + auto solver3 = Ikarus::makeGeneralSymEigenSolver(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); - // static_assert(Ikarus::Concepts::EigenValueSolver); - // static_assert(Ikarus::Concepts::EigenValueSolver); - // static_assert(Ikarus::Concepts::EigenValueSolver); + static_assert(Ikarus::Concepts::EigenValueSolver); + static_assert(Ikarus::Concepts::EigenValueSolver); + static_assert(Ikarus::Concepts::EigenValueSolver); + static_assert(Ikarus::Concepts::EigenValueSolver); return t; }