From 6cc0063aa81f076f9bf2a3cbb9581ad958714ace Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Sun, 1 Dec 2024 18:45:57 +0100 Subject: [PATCH] review part 1 --- ikarus/finiteelements/fehelper.hh | 1 - ikarus/finiteelements/ferequirements.hh | 2 +- .../mechanics/enhancedassumedstrains.hh | 5 +- .../mechanics/kirchhoffloveshell.hh | 3 +- ikarus/finiteelements/mechanics/truss.hh | 26 +++++----- ikarus/python/test/kltest.py | 2 +- ikarus/python/test/modalanalysistest.py | 2 +- .../generalizedeigensolver.hh | 7 ++- ikarus/utils/modalanalysis/lumpingschemes.hh | 2 +- ikarus/utils/modalanalysis/modalanalysis.hh | 4 +- tests/src/dummyproblem.hh | 4 +- tests/src/testeigenvaluesolver.cpp | 51 +++++-------------- tests/src/testmodalanalysis.cpp | 42 +++------------ 13 files changed, 46 insertions(+), 105 deletions(-) diff --git a/ikarus/finiteelements/fehelper.hh b/ikarus/finiteelements/fehelper.hh index 5c9d8a9d7..ad43903ac 100644 --- a/ikarus/finiteelements/fehelper.hh +++ b/ikarus/finiteelements/fehelper.hh @@ -127,5 +127,4 @@ template void globalIndices(const FiniteElement& fe, std::vector& globalIndices) { globalIndicesFromLocalView(fe.localView(), globalIndices); } - } // namespace Ikarus::FEHelper diff --git a/ikarus/finiteelements/ferequirements.hh b/ikarus/finiteelements/ferequirements.hh index a053b98d0..1d6d11409 100644 --- a/ikarus/finiteelements/ferequirements.hh +++ b/ikarus/finiteelements/ferequirements.hh @@ -203,7 +203,7 @@ auto scalarAffordance(VectorAffordance affordanceV) { namespace AffordanceCollections { inline constexpr AffordanceCollection elastoStatics(ScalarAffordance::mechanicalPotentialEnergy, VectorAffordance::forces, MatrixAffordance::stiffness); - constexpr Ikarus::AffordanceCollection dynamics(Ikarus::ScalarAffordance::noAffordance, + constexpr Ikarus::AffordanceCollection modalAnalysis(Ikarus::ScalarAffordance::noAffordance, Ikarus::VectorAffordance::noAffordance, Ikarus::MatrixAffordance::linearMass); } // namespace AffordanceCollections diff --git a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh index ba01a7c0f..bf6461b1e 100644 --- a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh +++ b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh @@ -159,13 +159,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::linearMass) { - return; - } using namespace Dune::DerivativeDirections; using namespace Dune; easApplicabilityCheck(); - if (isDisplacementBased()) + if (isDisplacementBased() or affordance != MatrixAffordance::stiffness) return; decltype(auto) LMat = [this]() -> decltype(auto) { diff --git a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh index 0d8e8673a..c1b74d1e5 100644 --- a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh +++ b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh @@ -290,6 +290,7 @@ protected: 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(); @@ -303,7 +304,7 @@ protected: for (size_t j = 0; j < numberOfNodes_; ++j) { nopJ.diagonal().setConstant(N[j]); M.template block(i * worldDim, j * worldDim) += - nopI.transpose() * density_ * nopJ * intElement; + nopI.transpose() * rhoT * nopJ * intElement; } } } diff --git a/ikarus/finiteelements/mechanics/truss.hh b/ikarus/finiteelements/mechanics/truss.hh index 186891366..5b18bad65 100644 --- a/ikarus/finiteelements/mechanics/truss.hh +++ b/ikarus/finiteelements/mechanics/truss.hh @@ -74,6 +74,7 @@ public: struct KinematicVariables { double L; ///< Length of the reference geometry + Eigen::VectorXd A1; ///< Length of the reference geometry as vector ST l; ///< Length of the deformed geometry ST Elin; ///< Linear strain ST Egl; ///< Green-Lagrange strain @@ -136,8 +137,9 @@ public: const double Lsquared = A1.squaredNorm(); const ST lsquared = (x2 - x1).squaredNorm(); - kin.L = sqrt(Lsquared); - kin.l = sqrt(lsquared); + kin.L = sqrt(Lsquared); + kin.l = sqrt(lsquared); + kin.A1 = A1; // Linear strain kin.Elin = (kin.l - kin.L) / kin.L; @@ -164,12 +166,8 @@ public: * * \return std::pair of length, and T */ - auto computeLengthAndTransformationMatrix() const { - auto& ele = underlying().localView().element(); - const auto X1 = Dune::toEigen(ele.geometry().corner(0)); - const auto X2 = Dune::toEigen(ele.geometry().corner(1)); - const auto A1 = X2 - X1; - const auto L = A1.norm(); + auto computeLengthAndTransformationMatrix(const Requirement& par) const { + const auto [L, A1, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par); const auto A1normed = A1 / L; auto T = Eigen::Matrix, 2, worldDim * 2>::Zero().eval(); @@ -194,8 +192,8 @@ public: requires(supportsResultType()) auto calculateAtImpl(const Requirement& req, [[maybe_unused]] const Dune::FieldVector& local, Dune::PriorityTag<0>) const { - using RTWrapper = ResultWrapper, ResultShape::Vector>; - const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(req); + using RTWrapper = ResultWrapper, ResultShape::Vector>; + const auto [L, _, l, Elin, Egl, dEdu, ddEddu] = computeStrain(req); if constexpr (isSameResultType) { auto N = Eigen::Vector{E_ * A_ * Egl * l / L}; // Axial force in deformed configuration return RTWrapper{N}; @@ -225,14 +223,14 @@ protected: const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K, const std::optional>>& dx = std::nullopt) const { if (affordance == MatrixAffordance::stiffness) { - const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx); + const auto [L, _, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx); K += E_ * A_ * L * (dEdu * dEdu.transpose() + ddEddu * Egl); } else if (affordance == MatrixAffordance::linearMass) { Eigen::Matrix mLoc{ {2, 1}, {1, 2} }; - auto [l, T] = computeLengthAndTransformationMatrix(); + auto [l, T] = computeLengthAndTransformationMatrix(par); mLoc *= (density_ * A_ * l) / 6.0; K += T.transpose() * mLoc * T; } else @@ -243,7 +241,7 @@ protected: auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance, const std::optional>>& dx = std::nullopt) const -> ScalarType { - const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx); + const auto [L, l, _, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx); return 0.5 * E_ * A_ * L * Egl * Egl; } @@ -251,7 +249,7 @@ protected: void calculateVectorImpl( const Requirement& par, VectorAffordance affordance, typename Traits::template VectorType force, const std::optional>>& dx = std::nullopt) const { - const auto [L, l, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx); + const auto [L, l, _, Elin, Egl, dEdu, ddEddu] = computeStrain(par, dx); force += E_ * A_ * Egl * L * dEdu; } }; diff --git a/ikarus/python/test/kltest.py b/ikarus/python/test/kltest.py index f18eb4dd1..d968f42f9 100644 --- a/ikarus/python/test/kltest.py +++ b/ikarus/python/test/kltest.py @@ -69,7 +69,7 @@ def vL(x, lambdaVal): vLoad = iks.finite_elements.volumeLoad3D(vL) klShell = iks.finite_elements.kirchhoffLoveShell( - youngs_modulus=1000, nu=0.0, thickness=thickness, density=density + youngs_modulus=1000, nu=0.0, thickness=thickness ) fes = [] diff --git a/ikarus/python/test/modalanalysistest.py b/ikarus/python/test/modalanalysistest.py index bc70ed8c9..5e75635cb 100644 --- a/ikarus/python/test/modalanalysistest.py +++ b/ikarus/python/test/modalanalysistest.py @@ -32,7 +32,7 @@ def setUp(self): self.fes = [] linMat = iks.materials.LinearElasticity(E=1000, nu=0.2).asPlaneStrain() - linElastic = finite_elements.linearElastic(linMat) + linElastic = finite_elements.linearElastic(linMat, density=100) for e in self.grid.elements: self.fes.append(finite_elements.makeFE(basis, linElastic)) self.fes[-1].bind(e) diff --git a/ikarus/solver/eigenvaluesolver/generalizedeigensolver.hh b/ikarus/solver/eigenvaluesolver/generalizedeigensolver.hh index 6fa66727a..d69f9da96 100644 --- a/ikarus/solver/eigenvaluesolver/generalizedeigensolver.hh +++ b/ikarus/solver/eigenvaluesolver/generalizedeigensolver.hh @@ -3,7 +3,7 @@ /** * \file generalizedeigensolver.hh - * \brief Helper for dune-functions + * \brief Implementation of wrapper classes for solving a generalized eigenvalue problem */ #pragma once @@ -26,7 +26,6 @@ namespace Ikarus { /** * \brief Construct a new make enum object - * */ MAKE_ENUM(EigenValueSolverType, Spectra, Eigen); @@ -191,7 +190,7 @@ struct GeneralizedSymEigenSolver } template - GeneralizedSymEigenSolver(const std::shared_ptr assemblerA, const std::shared_ptr& assemblerB) + GeneralizedSymEigenSolver(const std::shared_ptr& assemblerA, const std::shared_ptr& assemblerB) : GeneralizedSymEigenSolver(assemblerA->matrix(), assemblerB->matrix()) {} /** @@ -305,7 +304,7 @@ struct PartialGeneralizedSymEigenSolver Spectra::SortRule sortRule = Spectra::SortRule::SmallestAlge, ScalarType tolerance = 1e-10, Eigen::Index maxit = 1000) { solver_.init(); - solver_.compute(Spectra::SortRule::SmallestAlge, 1000, 1e-10, Spectra::SortRule::SmallestAlge); + solver_.compute(selection, tolerance, maxit, sortRule); computed_ = solver_.info() == Spectra::CompInfo::Successful; return computed_; diff --git a/ikarus/utils/modalanalysis/lumpingschemes.hh b/ikarus/utils/modalanalysis/lumpingschemes.hh index f889b821f..ad1ecf881 100644 --- a/ikarus/utils/modalanalysis/lumpingschemes.hh +++ b/ikarus/utils/modalanalysis/lumpingschemes.hh @@ -15,7 +15,7 @@ namespace Ikarus::Dynamics::LumpingSchemes { /** - * \brief Implements a row-sum lumping schemes for assemblermanipulator. Works with dense and sprase matrices. + * \brief Implements the operator() for a row-sum lumping scheme that can be passed to the assemblermanipulator to lump the already assembled global mass matrix. */ struct RowSumLumping { diff --git a/ikarus/utils/modalanalysis/modalanalysis.hh b/ikarus/utils/modalanalysis/modalanalysis.hh index 983ccbf2f..81461e790 100644 --- a/ikarus/utils/modalanalysis/modalanalysis.hh +++ b/ikarus/utils/modalanalysis/modalanalysis.hh @@ -66,7 +66,7 @@ struct ModalAnalysis req_.insertGlobalSolution(d_); stiffAssembler_->bind(req_, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced); - massAssembler_->bind(req_, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced); + massAssembler_->bind(req_, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced); lumpedMassAssembler_ = makeAssemblerManipulator(*massAssembler_); } @@ -76,7 +76,7 @@ struct ModalAnalysis * It resets already bound matrix manipulation functions. * * \tparam LumpingScheme The type of the lumping scheme, for example one found at \file - * ikarus/utils/modalanalysis/lumpingschemes.hh.s + * ikarus/utils/modalanalysis/lumpingschemes.hh. * \param ls the instantiated LumpingScheme object (pass either by value or by template definition). */ template diff --git a/tests/src/dummyproblem.hh b/tests/src/dummyproblem.hh index e1926273e..f0ea1e713 100644 --- a/tests/src/dummyproblem.hh +++ b/tests/src/dummyproblem.hh @@ -43,10 +43,8 @@ struct DummyProblem // YASPGrid needs an int, structuresgridfactory an unsigned int explicit DummyProblem( - const std::array, 2>& elementsPerDirection = {10, 10}) + const std::array, 2>& elementsPerDirection = {10, 10}, double Lx = 4.0, double Ly = 4.0) : grid_([&]() { - constexpr double Lx = 4.0; - constexpr double Ly = 4.0; const Dune::FieldVector bbox = {Lx, Ly}; if constexpr (not useYASP) diff --git a/tests/src/testeigenvaluesolver.cpp b/tests/src/testeigenvaluesolver.cpp index 812fcb87d..ff922d797 100644 --- a/tests/src/testeigenvaluesolver.cpp +++ b/tests/src/testeigenvaluesolver.cpp @@ -2,9 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later #include -#include "tests/src/testhelpers.hh" - -#include +#include "dummyproblem.hh" +#include "testhelpers.hh" #include #include @@ -63,39 +62,15 @@ auto testEigenVectors(const SOL1& solver1, const SOL2& solver2, std::shared_ptr< return t; } -auto testRealWorldProblem() { - 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 = {4, 4}; - auto grid = std::make_shared(bbox, elementsPerDirection); - - auto gridView = grid->leafGridView(); - - using namespace Dune::Functions::BasisFactory; - auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>())); - - auto vL = []([[maybe_unused]] auto& globalCoord, auto& lamb) { return Eigen::Vector2d{0, -1}; }; - auto linMat = Ikarus::LinearElasticity(Ikarus::toLamesFirstParameterAndShearModulus({.emodul = 100, .nu = 0.2})); - auto skills_ = Ikarus::skills(Ikarus::linearElastic(Ikarus::planeStress(linMat)), Ikarus::volumeLoad<2>(vL)); - - using FEType = decltype(Ikarus::makeFE(basis, skills_)); - std::vector fes; - for (auto&& element : elements(gridView)) { - fes.emplace_back(Ikarus::makeFE(basis, skills_)); - fes.back().bind(element); - } - - auto req = FEType::Requirement(); - typename FEType::Requirement::SolutionVectorType d; - d.setZero(basis.flat().size()); - req.insertGlobalSolution(d); +auto testEigenSolvers() { + TestSuite t("EigenValueSolver Test"); + using Grid = Dune::YaspGrid<2>; + using GridView = Grid::LeafGridView; - auto dirichletValues = Ikarus::DirichletValues(basis.flat()); - dirichletValues.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }); + DummyProblem testCase({5, 5}); + auto& req = testCase.requirement(); + auto& fes = testCase.finiteElements(); + auto& dirichletValues = testCase.dirichletValues(); auto assM = Ikarus::makeSparseFlatAssembler(fes, dirichletValues); auto assK = Ikarus::makeSparseFlatAssembler(fes, dirichletValues); @@ -103,10 +78,10 @@ auto testRealWorldProblem() { auto assMD = Ikarus::makeDenseFlatAssembler(fes, dirichletValues); auto assKD = Ikarus::makeDenseFlatAssembler(fes, dirichletValues); - assM->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced); + assM->bind(req, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced); assK->bind(req, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced); - assMD->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced); + assMD->bind(req, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced); assKD->bind(req, Ikarus::AffordanceCollections::elastoStatics, Ikarus::DBCOption::Reduced); int nev = 10; // number of requested eigenvalues @@ -192,6 +167,6 @@ int main(int argc, char** argv) { Ikarus::init(argc, argv); TestSuite t; - t.subTest(testRealWorldProblem()); + t.subTest(testEigenSolvers()); return t.exit(); } diff --git a/tests/src/testmodalanalysis.cpp b/tests/src/testmodalanalysis.cpp index 16b66d005..54684f3a5 100644 --- a/tests/src/testmodalanalysis.cpp +++ b/tests/src/testmodalanalysis.cpp @@ -2,10 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later #include +#include "tests/src/dummyproblem.hh" #include "tests/src/testhelpers.hh" -#include - #include #include #include @@ -25,34 +24,14 @@ using Dune::TestSuite; static auto modalAnalysisTest() { - TestSuite t("ModalAnalysis test"); + TestSuite t("ModalAnalysis Test"); using Grid = Dune::YaspGrid<2>; - const double Lx = 4.0; - const double Ly = 4.0; - Dune::FieldVector bbox = {Lx, Ly}; - std::array elementsPerDirection = {5, 5}; - auto grid = std::make_shared(bbox, elementsPerDirection); - - auto gridView = grid->leafGridView(); - - using namespace Dune::Functions::BasisFactory; - auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>())); - - auto linMat = Ikarus::LinearElasticity(Ikarus::toLamesFirstParameterAndShearModulus({.emodul = 100, .nu = 0.2})); - auto skills_ = Ikarus::skills(Ikarus::linearElastic(Ikarus::planeStress(linMat))); - - using FEType = decltype(Ikarus::makeFE(basis, skills_)); - std::vector fes; - for (auto&& element : elements(gridView)) { - fes.emplace_back(Ikarus::makeFE(basis, skills_)); - fes.back().bind(element); - } - - auto dirichletValues = Ikarus::DirichletValues(basis.flat()); - dirichletValues.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }); - - auto mA = Ikarus::Dynamics::ModalAnalysis(std::move(fes), dirichletValues); + DummyProblem testCase({5, 5}); + auto& req = testCase.requirement(); + auto& fes = testCase.finiteElements(); + auto& dirichletValues = testCase.dirichletValues(); + auto mA = Ikarus::Dynamics::ModalAnalysis(std::move(fes), dirichletValues); mA.compute(); auto frequencies = mA.angularFrequencies(); @@ -70,13 +49,8 @@ static auto modalAnalysisTest() { mA.writeEigenModes("eigenforms", 20); - auto req = FEType::Requirement(); - typename FEType::Requirement::SolutionVectorType d; - d.setZero(basis.flat().size()); - req.insertGlobalSolution(d); - auto massAssembler = Ikarus::makeSparseFlatAssembler(fes, dirichletValues); - massAssembler->bind(req, Ikarus::AffordanceCollections::dynamics, Ikarus::DBCOption::Reduced); + massAssembler->bind(req, Ikarus::AffordanceCollections::modalAnalysis, Ikarus::DBCOption::Reduced); auto lumpedAssembler = Ikarus::Dynamics::makeLumpedFlatAssembler(massAssembler); auto M = massAssembler->matrix();