From 936b375d6ceef468bc81c692398ab428667905df Mon Sep 17 00:00:00 2001 From: Maarten Hodzelmans Date: Fri, 26 Jun 2026 15:22:35 +0200 Subject: [PATCH 1/6] Add hook for custom solve against mass matrix. Externalize existing explicit inversion. --- src/gsDynamicSolvers/gsDynamicBase.h | 76 ++++++++++++---- .../gsDynamicExplicitEuler.hpp | 22 ++--- src/gsDynamicSolvers/gsDynamicMassInverse.h | 86 +++++++++++++++++++ src/gsDynamicSolvers/gsDynamicRK4.hpp | 46 ++++++---- 4 files changed, 184 insertions(+), 46 deletions(-) create mode 100644 src/gsDynamicSolvers/gsDynamicMassInverse.h diff --git a/src/gsDynamicSolvers/gsDynamicBase.h b/src/gsDynamicSolvers/gsDynamicBase.h index 00630a2..e1e41fd 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; + + mass_is_timedependent = 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,33 @@ 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 mass_is_timedependent || !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; + m_massInverse->apply(rhs,result); } /// Compute the damping matrix @@ -342,8 +378,10 @@ class gsDynamicBase index_t m_numDofs; Mass_t m_mass; - mutable gsSparseMatrix m_massInv; TMass_t m_Tmass; + bool mass_is_timedependent; + 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..d105b05 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicMassInverse.h @@ -0,0 +1,86 @@ +/** @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 + +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(mass); + gsMatrix inverse = solver.solve(eye); + 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)) + {} + + void compute(const gsSparseMatrix & mass) + { + m_solver->compute(mass); + } + + void apply(const gsVector & rhs, gsVector & result) const + { + 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); From 654904be10b27e2b9725e9ac9ff5202fbff896e5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 26 Jun 2026 13:55:58 +0000 Subject: [PATCH 2/6] Rename mass_is_timedependent to m_massIsTimeDependent --- src/gsDynamicSolvers/gsDynamicBase.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gsDynamicSolvers/gsDynamicBase.h b/src/gsDynamicSolvers/gsDynamicBase.h index e1e41fd..a5b415a 100644 --- a/src/gsDynamicSolvers/gsDynamicBase.h +++ b/src/gsDynamicSolvers/gsDynamicBase.h @@ -185,7 +185,7 @@ class gsDynamicBase m_status = gsStatus::NotStarted; - mass_is_timedependent = mass_timedependent; + m_massIsTimeDependent = mass_timedependent; m_massInverse = MassInversePtr(new gsDynamicExplicitMassInverse); m_massInverseComputed = false; } @@ -327,7 +327,7 @@ class gsDynamicBase /// Return true if the inverse mass representation needs a fresh mass matrix. virtual bool _massInverseNeedsUpdate() const { - return mass_is_timedependent || !m_massInverseComputed; + return m_massIsTimeDependent || !m_massInverseComputed; } /// Compute or update the inverse mass representation. @@ -379,7 +379,7 @@ class gsDynamicBase Mass_t m_mass; TMass_t m_Tmass; - bool mass_is_timedependent; + bool m_massIsTimeDependent; mutable MassInversePtr m_massInverse; mutable bool m_massInverseComputed; From f999b7814f973d3622ea60eaa4b1596ee1a94f25 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 26 Jun 2026 13:58:08 +0000 Subject: [PATCH 3/6] Initialize result vector before calling apply in _applyMassInverse --- src/gsDynamicSolvers/gsDynamicBase.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/gsDynamicSolvers/gsDynamicBase.h b/src/gsDynamicSolvers/gsDynamicBase.h index a5b415a..4a3e7c0 100644 --- a/src/gsDynamicSolvers/gsDynamicBase.h +++ b/src/gsDynamicSolvers/gsDynamicBase.h @@ -350,6 +350,7 @@ class gsDynamicBase this->_computeMass(time,M); this->_computeMassInverse(M); } + result.setZero(rhs.size()); m_massInverse->apply(rhs,result); } From 05fd719f0761a9738754f43068293f9d6d45833a Mon Sep 17 00:00:00 2001 From: Maarten Hodzelmans Date: Fri, 26 Jun 2026 15:58:46 +0200 Subject: [PATCH 4/6] Check solver success/failure Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/gsDynamicSolvers/gsDynamicMassInverse.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/gsDynamicSolvers/gsDynamicMassInverse.h b/src/gsDynamicSolvers/gsDynamicMassInverse.h index d105b05..d33ff21 100644 --- a/src/gsDynamicSolvers/gsDynamicMassInverse.h +++ b/src/gsDynamicSolvers/gsDynamicMassInverse.h @@ -67,14 +67,17 @@ class gsDynamicSparseMassInverse : public gsDynamicMassInverse explicit gsDynamicSparseMassInverse(const std::string & solver = "SimplicialLDLT") : m_solver(gsSparseSolver::get(solver)) - {} + { + GISMO_ENSURE(m_solver, "Unknown sparse solver: "< & mass) + void compute(const gsSparseMatrix & 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 + void apply(const gsVector & rhs, gsVector & result) const override { result = m_solver->solve(rhs); } From e8efe32cb2f668384e5f8135bc47baef30780721 Mon Sep 17 00:00:00 2001 From: Maarten Hodzelmans Date: Fri, 26 Jun 2026 16:03:28 +0200 Subject: [PATCH 5/6] Check solver success in explicit inversion --- src/gsDynamicSolvers/gsDynamicMassInverse.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/gsDynamicSolvers/gsDynamicMassInverse.h b/src/gsDynamicSolvers/gsDynamicMassInverse.h index d33ff21..4fab0b0 100644 --- a/src/gsDynamicSolvers/gsDynamicMassInverse.h +++ b/src/gsDynamicSolvers/gsDynamicMassInverse.h @@ -13,6 +13,7 @@ #pragma once +#include #include #include @@ -46,7 +47,11 @@ class gsDynamicExplicitMassInverse : public gsDynamicMassInverse gsSparseMatrix eye(mass.rows(), mass.cols()); eye.setIdentity(); typename gsSparseSolver::LU solver(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(); } From 000fa90b6b6decc54684570ceb38ac16896449fa Mon Sep 17 00:00:00 2001 From: Maarten Hodzelmans Date: Fri, 26 Jun 2026 16:11:17 +0200 Subject: [PATCH 6/6] Copilot suggestion, compute LU by call rather than construction Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/gsDynamicSolvers/gsDynamicMassInverse.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gsDynamicSolvers/gsDynamicMassInverse.h b/src/gsDynamicSolvers/gsDynamicMassInverse.h index 4fab0b0..17aa183 100644 --- a/src/gsDynamicSolvers/gsDynamicMassInverse.h +++ b/src/gsDynamicSolvers/gsDynamicMassInverse.h @@ -46,7 +46,8 @@ class gsDynamicExplicitMassInverse : public gsDynamicMassInverse { gsSparseMatrix eye(mass.rows(), mass.cols()); eye.setIdentity(); - typename gsSparseSolver::LU solver(mass); + typename gsSparseSolver::LU solver; + solver.compute(mass); GISMO_ENSURE(solver.info() == gsEigen::Success, "LU factorization of mass matrix failed."); gsMatrix inverse = solver.solve(eye);