diff --git a/ikarus/finiteelements/fehelper.hh b/ikarus/finiteelements/fehelper.hh index ad43903ac..8893cc9c1 100644 --- a/ikarus/finiteelements/fehelper.hh +++ b/ikarus/finiteelements/fehelper.hh @@ -127,4 +127,46 @@ template void globalIndices(const FiniteElement& fe, std::vector& 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 +void subspaceEntriesToElementEntries(const auto& localMatrix, auto& elementMatrix, const auto& globalBasis, + TreePaths&&... tps) { + auto treePathTuple = Dune::TupleVector(std::forward(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 +void subspaceEntriesToElementEntries(const auto& localMatrix, auto& elementMatrix, const auto& globalBasis) { + Dune::Hybrid::forEach( + Dune::Hybrid::integralRange(Dune::index_constant(), Dune::index_constant()), + [&](const auto i) { subspaceEntriesToElementEntries(localMatrix, elementMatrix, globalBasis, i); }); +} + } // namespace Ikarus::FEHelper diff --git a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh index a9ecbe525..0ebfc9867 100644 --- a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh +++ b/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>>& dx = std::nullopt) const { + if (affordance == MatrixAffordance::mass) { + return; + } using namespace Dune::DerivativeDirections; using namespace Dune; easApplicabilityCheck(); diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh index 01c1bee5b..40bf93969 100644 --- a/ikarus/finiteelements/mechanics/linearelastic.hh +++ b/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()), [&](auto k) { - auto& basis = underlying().localView().tree().child(k); - M(basis.localIndex(i), basis.localIndex(j)) += Mij(i, j); - }); - } - } + FEHelper::subspaceEntriesToElementEntries(Mij, M, underlying().localView().tree()); } - M.template triangularView() = M.transpose(); } template diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh index 8a5697919..fd40312ba 100644 --- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh +++ b/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>>& dx = std::nullopt) const { + if (affordance == MatrixAffordance::mass) { + calculateMassImpl(par, affordance, K); + return; + } using namespace Dune::DerivativeDirections; using namespace Dune; const auto uFunction = displacementFunction(par, dx); @@ -264,6 +268,19 @@ protected: } } + template + 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(Mij, M, underlying().localView().tree()); + } + } + template auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance, const std::optional>>& dx = diff --git a/ikarus/utils/dynamics/dynamics.hh b/ikarus/utils/dynamics/dynamics.hh index a81095d53..7c9c6ed69 100644 --- a/ikarus/utils/dynamics/dynamics.hh +++ b/ikarus/utils/dynamics/dynamics.hh @@ -7,6 +7,8 @@ */ #pragma once +#include + #include #include @@ -40,13 +42,15 @@ namespace Impl { template inline auto lumpingSparse() { - return - [](const auto&, const auto&, auto, auto, Eigen::SparseMatrix& mat) { return Impl::lumpingSparseImpl(mat); }; + return [](const auto&, const auto&, auto, auto, Eigen::SparseMatrix& mat) { + return Impl::lumpingSparseImpl(mat); + }; } template inline auto lumpingDense() { - return [](const auto&, const auto&, auto, auto, Eigen::MatrixX& mat) { return Impl::lumpingDenseImpl(mat); }; + return + [](const auto&, const auto&, auto, auto, Eigen::MatrixX& mat) { return Impl::lumpingDenseImpl(mat); }; } template @@ -76,4 +80,24 @@ void writeEigenformsToVTK(const Eigensolver& solver, std::shared_ptr writer.write(filename); } +template +void writeEigenformsToPVD(const Eigensolver& solver, std::shared_ptr assembler, const std::string& filename, + std::optional nev_ = std::nullopt) { + auto nev = nev_.value_or(solver.nev()); + auto eigenvectors = solver.eigenvectors(); + auto basis = assembler->basis(); + + auto writer = std::make_shared(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 diff --git a/tests/src/testdynamics.cpp b/tests/src/testdynamics.cpp index 84a65199c..7b8a7e44a 100644 --- a/tests/src/testdynamics.cpp +++ b/tests/src/testdynamics.cpp @@ -41,7 +41,7 @@ static auto dynamicsTest() { const double Lx = 4.0; const double Ly = 4.0; Dune::FieldVector bbox = {Lx, Ly}; - std::array elementsPerDirection = {6, 6}; + std::array elementsPerDirection = {4, 4}; auto grid = std::make_shared(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(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); static_assert(Ikarus::Concepts::EigenValueSolver); @@ -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(); } diff --git a/tests/src/testvtkwriter.cpp b/tests/src/testvtkwriter.cpp index 1fa3d6958..53b515022 100644 --- a/tests/src/testvtkwriter.cpp +++ b/tests/src/testvtkwriter.cpp @@ -50,11 +50,11 @@ auto vtkWriterTest() { DummyProblem 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};