diff --git a/src/gsDynamicSolvers/gsDynamicBase.h b/src/gsDynamicSolvers/gsDynamicBase.h index 00630a2..4a3e7c0 100644 --- a/src/gsDynamicSolvers/gsDynamicBase.h +++ b/src/gsDynamicSolvers/gsDynamicBase.h @@ -18,6 +18,7 @@ #include #include +#include #include namespace gismo @@ -53,6 +54,8 @@ class gsDynamicBase /// @brief Callback that is evaluated before each timestep. typedef std::function StepCallback_t; + typedef gsDynamicMassInverse MassInverse_t; + typedef typename MassInverse_t::Ptr MassInversePtr; virtual ~gsDynamicBase() {}; /// Constructor @@ -73,7 +76,7 @@ class gsDynamicBase m_Tjacobian = [this](gsVector const & /*x*/, const T /*time*/, gsSparseMatrix & result) -> bool {return m_stiffness(result);}; m_Tforce = [this]( const T /*time*/, gsVector & result) -> bool {return m_force(result);}; m_Tresidual = [ ](gsVector const & /*x*/, const T /*time*/, gsVector & /*result*/) -> bool {GISMO_ERROR("time-dependent residual not available");}; - _init(); + _init(false); } /// Constructor @@ -93,7 +96,7 @@ class gsDynamicBase m_Tdamping = [this](gsVector const & x, const T /*time*/, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; m_Tjacobian = [this](gsVector const & /*x*/, const T /*time*/, gsSparseMatrix & result) -> bool {return m_stiffness(result);}; m_Tresidual = [ ](gsVector const & /*x*/, const T /*time*/, gsVector & /*result*/) -> bool {GISMO_ERROR("time-dependent residual not available");}; - _init(); + _init(false); } /// Constructor @@ -114,7 +117,7 @@ class gsDynamicBase m_Tjacobian = [this](gsVector const & x, const T /*time*/, gsSparseMatrix & result) -> bool {return m_jacobian(x,result);}; m_Tforce = [ ]( const T /*time*/, gsVector & /*result*/) -> bool {GISMO_ERROR("time-dependent force not available");}; m_Tresidual = [this](gsVector const & x, const T /*time*/, gsVector & result) -> bool {return m_residual(x,result);}; - _init(); + _init(false); } /// Constructor @@ -133,7 +136,7 @@ class gsDynamicBase m_Tmass = [this]( const T /*time*/, gsSparseMatrix & result) -> bool {return m_mass(result);}; m_Tjacobian = [this](gsVector const & x, const T /*time*/, gsSparseMatrix & result) -> bool {return m_jacobian(x,result);}; m_Tdamping = [this](gsVector const & x, const T /*time*/, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; - _init(); + _init(false); } /// Constructor @@ -151,7 +154,7 @@ class gsDynamicBase { m_Tmass = [this]( const T /*time*/, gsSparseMatrix & result) -> bool {return m_mass(result);}; m_Tdamping = [this](gsVector const & x, const T /*time*/, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; - _init(); + _init(false); } /// Constructor @@ -167,13 +170,13 @@ class gsDynamicBase m_Tjacobian(TJacobian), m_Tresidual(TResidual) { - _init(); + _init(true); } - gsDynamicBase() {_init();} + gsDynamicBase() {_init(true);} protected: - void _init() + void _init(const bool mass_timedependent) { m_time = 0; // initialize variables @@ -181,6 +184,10 @@ class gsDynamicBase defaultOptions(); m_status = gsStatus::NotStarted; + + m_massIsTimeDependent = mass_timedependent; + m_massInverse = MassInversePtr(new gsDynamicExplicitMassInverse); + m_massInverseComputed = false; } // General functions @@ -263,6 +270,21 @@ class gsDynamicBase m_preStepCallback = StepCallback_t(); } + /// @brief Set a custom inverse mass function object. + virtual void setMassInverse(const MassInversePtr & massInverse) + { + GISMO_ENSURE(massInverse, "Cannot set an empty inverse mass function object."); + m_massInverse = massInverse; + m_massInverseComputed = false; + } + + /// @brief Restore the default explicit inverse mass function object. + virtual void clearMassInverse() + { + m_massInverse = MassInversePtr(new gsDynamicExplicitMassInverse); + m_massInverseComputed = false; + } + /// Access the options virtual gsOptionList & options() {return m_options;}; @@ -302,19 +324,34 @@ class gsDynamicBase throw 2; } - /// Compute the mass matrix - virtual void _computeMassInverse(const gsSparseMatrix & M, gsSparseMatrix & Minv) const + /// Return true if the inverse mass representation needs a fresh mass matrix. + virtual bool _massInverseNeedsUpdate() const + { + return m_massIsTimeDependent || !m_massInverseComputed; + } + + /// Compute or update the inverse mass representation. + virtual void _computeMassInverse(const gsSparseMatrix & M) const + { + if (this->_massInverseNeedsUpdate()) + { + GISMO_ENSURE(M.rows()!=0 && M.cols()!=0, "Cannot compute inverse mass from an empty mass matrix."); + m_massInverse->compute(M); + m_massInverseComputed = true; + } + } + + /// Apply the inverse mass matrix without necessarily forming it explicitly. + virtual void _applyMassInverse(const T time, const gsVector & rhs, gsVector & result) const { - if ((m_mass==nullptr) || (m_massInv.rows()==0 || m_massInv.cols()==0)) // then mass is time-dependent or the mass inverse is not stored, compute it + if (this->_massInverseNeedsUpdate()) { - gsSparseMatrix eye(M.rows(), M.cols()); - eye.setIdentity(); - gsSparseSolver<>::LU solver(M); - gsMatrix MinvI = solver.solve(eye); - m_massInv = Minv = MinvI.sparseView(); + gsSparseMatrix M; + this->_computeMass(time,M); + this->_computeMassInverse(M); } - else - Minv = m_massInv; + result.setZero(rhs.size()); + m_massInverse->apply(rhs,result); } /// Compute the damping matrix @@ -342,8 +379,10 @@ class gsDynamicBase index_t m_numDofs; Mass_t m_mass; - mutable gsSparseMatrix m_massInv; TMass_t m_Tmass; + bool m_massIsTimeDependent; + mutable MassInversePtr m_massInverse; + mutable bool m_massInverseComputed; Damping_t m_damping; TDamping_t m_Tdamping; diff --git a/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp b/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp index ad9aa78..ad3c0af 100644 --- a/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp +++ b/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp @@ -61,19 +61,19 @@ gsDynamicExplicitEuler::_step_impl(const T t, const T dt, gsVector & U sol.topRows(N) = Uold; sol.bottomRows(N) = Vold; - gsVector F; - gsSparseMatrix M, Minv, C, K; + gsVector F, massRhs, massSolve; + gsSparseMatrix C, K; // Computed at t=t0 - this->_computeMass(t,M); - this->_computeMassInverse(M,Minv); this->_computeForce(t,F); this->_computeDamping(U,t,C); this->_computeJacobian(U,t,K); this->_initOutput(); sol.topRows(N) += dt * Vold; - sol.bottomRows(N) += dt * Minv * (F - K * Uold - C * Vold); + massRhs = F - K * Uold - C * Vold; + this->_applyMassInverse(t, massRhs, massSolve); + sol.bottomRows(N) += dt * massSolve; this->_stepOutput(0,sol.norm(),0.); gsDebugVar(sol.transpose()); @@ -101,18 +101,18 @@ gsDynamicExplicitEuler::_step_impl(const T t, const T dt, gsVector & U sol.topRows(N) = Uold; sol.bottomRows(N) = Vold; - gsVector R; - gsSparseMatrix M, Minv, C, K; + gsVector R, massRhs, massSolve; + gsSparseMatrix C, K; // Computed at t=t0 - this->_computeMass(t,M); - this->_computeMassInverse(M,Minv); this->_computeDamping(Uold,t,C); this->_computeResidual(Uold,t,R); this->_initOutput(); sol.topRows(N) += dt * Vold; - sol.bottomRows(N) += dt * Minv * ( - R - C * Vold); + massRhs = - R - C * Vold; + this->_applyMassInverse(t, massRhs, massSolve); + sol.bottomRows(N) += dt * massSolve; this->_stepOutput(0,sol.norm(),0.); U = sol.topRows(N); @@ -158,4 +158,4 @@ void gsDynamicExplicitEuler::_stepOutput(const index_t it, const T resnor } } -} // namespace gismo \ No newline at end of file +} // namespace gismo diff --git a/src/gsDynamicSolvers/gsDynamicMassInverse.h b/src/gsDynamicSolvers/gsDynamicMassInverse.h new file mode 100644 index 0000000..17aa183 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicMassInverse.h @@ -0,0 +1,95 @@ +/** @file gsDynamicMassInverse.h + + @brief Inverse mass function objects for dynamic solvers + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): M.M. Hodzelmans +*/ + +#pragma once + +#include +#include + +#include + +namespace gismo +{ + +/// @brief Function object that computes and applies an inverse mass matrix. +template +class gsDynamicMassInverse +{ +public: + typedef memory::shared_ptr Ptr; + + virtual ~gsDynamicMassInverse() {} + + /// Compute or update the inverse mass representation from \a mass. + virtual void compute(const gsSparseMatrix & mass) = 0; + + /// Apply the inverse mass representation, result = M^{-1} rhs. + virtual void apply(const gsVector & rhs, gsVector & result) const = 0; +}; + +/// @brief Default inverse mass object that explicitly stores M^{-1}. +template +class gsDynamicExplicitMassInverse : public gsDynamicMassInverse +{ +public: + void compute(const gsSparseMatrix & mass) + { + gsSparseMatrix eye(mass.rows(), mass.cols()); + eye.setIdentity(); + typename gsSparseSolver::LU solver; + solver.compute(mass); + GISMO_ENSURE(solver.info() == gsEigen::Success, + "LU factorization of mass matrix failed."); + gsMatrix inverse = solver.solve(eye); + GISMO_ENSURE(solver.info() == gsEigen::Success, + "Solving for inverse mass matrix failed."); + m_massInv = inverse.sparseView(); + } + + void apply(const gsVector & rhs, gsVector & result) const + { + result = m_massInv * rhs; + } + +private: + gsSparseMatrix m_massInv; +}; + +/// @brief Inverse mass object backed by a sparse factorization. +template +class gsDynamicSparseMassInverse : public gsDynamicMassInverse +{ +public: + explicit gsDynamicSparseMassInverse(const std::string & solver = "SimplicialLDLT") + : + m_solver(gsSparseSolver::get(solver)) + { + GISMO_ENSURE(m_solver, "Unknown sparse solver: "< & mass) override + { + m_solver->compute(mass); + GISMO_ENSURE(m_solver->info()==gsEigen::ComputationInfo::Success, "Mass matrix factorization failed."); + } + + void apply(const gsVector & rhs, gsVector & result) const override + { + result = m_solver->solve(rhs); + } + +private: + typename gsSparseSolver::uPtr m_solver; +}; + +} // namespace gismo diff --git a/src/gsDynamicSolvers/gsDynamicRK4.hpp b/src/gsDynamicSolvers/gsDynamicRK4.hpp index d6b691a..fc1c228 100644 --- a/src/gsDynamicSolvers/gsDynamicRK4.hpp +++ b/src/gsDynamicSolvers/gsDynamicRK4.hpp @@ -31,12 +31,11 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector sol.topRows(N) = Uold; sol.bottomRows(N) = Vold; - gsVector F, R; //R is the residual vector, and the new _computeForce give the inline factor - gsSparseMatrix M, Minv, C, K; + //R is the residual vector, and the new _computeForce give the inline factor + gsVector F, R, massRhs, massSolve; + gsSparseMatrix C, K; // Computed at t=t0 - this->_computeMass(t,M); - this->_computeMassInverse(M,Minv); // this->_computeForce(t,F); this->_computeDamping(U,t,C); //C is damping this->_computeJacobian(U,t,K); @@ -50,7 +49,9 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector _computeForce(t, F); R = F - K * Uold; k1.topRows(N) = Vold; - k1.bottomRows(N) = Minv * (R - C * Vold); + massRhs = R - C * Vold; + this->_applyMassInverse(t, massRhs, massSolve); + k1.bottomRows(N) = massSolve; //Step2 (calculate k2) Utmp = Uold + dt/2. * k1.topRows(N); @@ -58,7 +59,9 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector _computeForce(t + dt/2.,F); R = F - K * Utmp; k2.topRows(N) = Vtmp; - k2.bottomRows(N) = Minv * ( R - C * Vtmp); + massRhs = R - C * Vtmp; + this->_applyMassInverse(t + dt/2., massRhs, massSolve); + k2.bottomRows(N) = massSolve; //Step3 (calculate k3) Utmp = Uold + dt/2. * k2.topRows(N); @@ -66,7 +69,9 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector _computeForce(t + dt/2., F); R = F - K * Utmp; k3.topRows(N) = Vtmp; - k3.bottomRows(N) = Minv * ( R - C * Vtmp); + massRhs = R - C * Vtmp; + this->_applyMassInverse(t + dt/2., massRhs, massSolve); + k3.bottomRows(N) = massSolve; //Step4 (calculate k4) Utmp = Uold + dt * k3.topRows(N); @@ -74,7 +79,9 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector _computeForce(t + dt, F); R = F - K * Utmp; k4.topRows(N) = Vtmp; - k4.bottomRows(N) = Minv * ( R - C * Vtmp); + massRhs = R - C * Vtmp; + this->_applyMassInverse(t + dt, massRhs, massSolve); + k4.bottomRows(N) = massSolve; sol += 1./6 * dt * (k1 + 2.*k2 + 2.*k3 + k4); @@ -104,12 +111,11 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector sol.topRows(N) = Uold; sol.bottomRows(N) = Vold; - gsVector F, R; //R is the residual vector, and the new _computeForce give the inline factor - gsSparseMatrix M, Minv, C, K; + //R is the residual vector, and the new _computeForce give the inline factor + gsVector F, R, massRhs, massSolve; + gsSparseMatrix C, K; // Computed at t=t0 - this->_computeMass(t,M); - this->_computeMassInverse(M,Minv); // this->_computeForce(t,F); this->_computeDamping(U,t,C); //C is damping this->_computeJacobian(U,t,K); @@ -122,28 +128,36 @@ gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector //Step1 (calculate k1) _computeResidual(Uold, t, R); k1.topRows(N) = Vold; - k1.bottomRows(N) = Minv * (R - C * Vold); + massRhs = R - C * Vold; + this->_applyMassInverse(t, massRhs, massSolve); + k1.bottomRows(N) = massSolve; //Step2 (calculate k2) Utmp = Uold + dt/2. * k1.topRows(N); Vtmp = Vold + dt/2. * k1.bottomRows(N); _computeResidual(Utmp,t + dt/2.,R); k2.topRows(N) = Vtmp; - k2.bottomRows(N) = Minv * ( R - C * Vtmp); + massRhs = R - C * Vtmp; + this->_applyMassInverse(t + dt/2., massRhs, massSolve); + k2.bottomRows(N) = massSolve; //Step3 (calculate k3) Utmp = Uold + dt/2. * k2.topRows(N); Vtmp = Vold + dt/2. * k2.bottomRows(N); _computeResidual(Utmp,t + dt/2.,R); k3.topRows(N) = Vtmp; - k3.bottomRows(N) = Minv * ( R - C * Vtmp); + massRhs = R - C * Vtmp; + this->_applyMassInverse(t + dt/2., massRhs, massSolve); + k3.bottomRows(N) = massSolve; //Step4 (calculate k4) Utmp = Uold + dt * k3.topRows(N); Vtmp = Vold + dt * k3.bottomRows(N); _computeResidual(Utmp,t + dt,R); k4.topRows(N) = Vtmp; - k4.bottomRows(N) = Minv * ( R - C * Vtmp); + massRhs = R - C * Vtmp; + this->_applyMassInverse(t + dt, massRhs, massSolve); + k4.bottomRows(N) = massSolve; sol += 1./6 * dt * (k1 + 2.*k2 + 2.*k3 + k4);