Hook for custom solve against mass matrix for explicit time integrators

src/gsDynamicSolvers/gsDynamicBase.h

@@ -18,6 +18,7 @@
 
 #include <gsCore/gsLinearAlgebra.h>
 #include <gsIO/gsOptionList.h>
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicMassInverse.h>
 #include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h>
 
 namespace gismo
@@ -53,6 +54,8 @@ class gsDynamicBase
     /// @brief Callback that is evaluated before each timestep.
     typedef std::function<bool(const T /* time */, const T /* dtime */, 
         const gsOptionList & /* solver options */)> StepCallback_t;
+    typedef gsDynamicMassInverse<T> MassInverse_t;
+    typedef typename MassInverse_t::Ptr MassInversePtr;
     virtual ~gsDynamicBase() {};
 
     /// Constructor
@@ -73,7 +76,7 @@ class gsDynamicBase
         m_Tjacobian   = [this](gsVector<T> const & /*x*/, const T /*time*/, gsSparseMatrix<T> & result)     -> bool {return m_stiffness(result);};
         m_Tforce      = [this](                           const T /*time*/, gsVector<T>       & result)     -> bool {return m_force(result);};
         m_Tresidual   = [    ](gsVector<T> const & /*x*/, const T /*time*/, gsVector<T>       & /*result*/) -> bool {GISMO_ERROR("time-dependent residual not available");};
-        _init();
+        _init(false);
     }
 
     /// Constructor
@@ -93,7 +96,7 @@ class gsDynamicBase
         m_Tdamping    = [this](gsVector<T> const & x,     const T /*time*/, gsSparseMatrix<T> & result)     -> bool {return m_damping(x,result);};
         m_Tjacobian   = [this](gsVector<T> const & /*x*/, const T /*time*/, gsSparseMatrix<T> & result)     -> bool {return m_stiffness(result);};
         m_Tresidual   = [    ](gsVector<T> const & /*x*/, const T /*time*/, gsVector<T>       & /*result*/) -> bool {GISMO_ERROR("time-dependent residual not available");};
-        _init();
+        _init(false);
     }
 
     /// Constructor
@@ -114,7 +117,7 @@ class gsDynamicBase
         m_Tjacobian   = [this](gsVector<T> const & x, const T /*time*/, gsSparseMatrix<T> & result)     -> bool {return m_jacobian(x,result);};
         m_Tforce      = [    ](                       const T /*time*/, gsVector<T>       & /*result*/) -> bool {GISMO_ERROR("time-dependent force not available");};
         m_Tresidual   = [this](gsVector<T> const & x, const T /*time*/, gsVector<T>       & result)     -> bool {return m_residual(x,result);};
-        _init();
+        _init(false);
     }
 
     /// Constructor
@@ -133,7 +136,7 @@ class gsDynamicBase
         m_Tmass       = [this](                       const T /*time*/, gsSparseMatrix<T> & result) -> bool {return m_mass(result);};
         m_Tjacobian   = [this](gsVector<T> const & x, const T /*time*/, gsSparseMatrix<T> & result) -> bool {return m_jacobian(x,result);};
         m_Tdamping    = [this](gsVector<T> const & x, const T /*time*/, gsSparseMatrix<T> & result) -> bool {return m_damping(x,result);};
-        _init();
+        _init(false);
     }
 
     /// Constructor
@@ -151,7 +154,7 @@ class gsDynamicBase
     {
         m_Tmass       = [this](                       const T /*time*/, gsSparseMatrix<T> & result) -> bool {return m_mass(result);};
         m_Tdamping    = [this](gsVector<T> const & x, const T /*time*/, gsSparseMatrix<T> & result) -> bool {return m_damping(x,result);};
-        _init();
+        _init(false);
     }
 
     /// Constructor
@@ -167,20 +170,24 @@ 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
         m_numIterations = -1;
         defaultOptions();
 
         m_status = gsStatus::NotStarted;
+
+        m_massIsTimeDependent = mass_timedependent;
+        m_massInverse = MassInversePtr(new gsDynamicExplicitMassInverse<T>);
+        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<T>);
+        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<T> & M, gsSparseMatrix<T> & 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<T> & 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<T> & rhs, gsVector<T> & 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<T> eye(M.rows(), M.cols());
-            eye.setIdentity();
-            gsSparseSolver<>::LU solver(M);
-            gsMatrix<T> MinvI = solver.solve(eye);
-            m_massInv = Minv = MinvI.sparseView();
+            gsSparseMatrix<T> 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<T> 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;

src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp

@@ -61,19 +61,19 @@ gsDynamicExplicitEuler<T,_NL>::_step_impl(const T t, const T dt, gsVector<T> & U
   sol.topRows(N) = Uold;
   sol.bottomRows(N) = Vold;
 
-  gsVector<T> F;
-  gsSparseMatrix<T> M, Minv, C, K;
+  gsVector<T> F, massRhs, massSolve;
+  gsSparseMatrix<T> 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<T,_NL>::_step_impl(const T t, const T dt, gsVector<T> & U
   sol.topRows(N) = Uold;
   sol.bottomRows(N) = Vold;
 
-  gsVector<T> R;
-  gsSparseMatrix<T> M, Minv, C, K;
+  gsVector<T> R, massRhs, massSolve;
+  gsSparseMatrix<T> 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<T,_NL>::_stepOutput(const index_t it, const T resnor
   }
 }
 
-} // namespace gismo
+} // namespace gismo

src/gsDynamicSolvers/gsDynamicMassInverse.h

@@ -0,0 +1,94 @@
+/** @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 <gsCore/gsDebug.h>
+#include <gsCore/gsLinearAlgebra.h>
+
+#include <string>
+
+namespace gismo
+{
+
+/// @brief Function object that computes and applies an inverse mass matrix.
+template <class T>
+class gsDynamicMassInverse
+{
+public:
+    typedef memory::shared_ptr<gsDynamicMassInverse> Ptr;
+
+    virtual ~gsDynamicMassInverse() {}
+
+    /// Compute or update the inverse mass representation from \a mass.
+    virtual void compute(const gsSparseMatrix<T> & mass) = 0;
+
+    /// Apply the inverse mass representation, result = M^{-1} rhs.
+    virtual void apply(const gsVector<T> & rhs, gsVector<T> & result) const = 0;
+};
+
+/// @brief Default inverse mass object that explicitly stores M^{-1}.
+template <class T>
+class gsDynamicExplicitMassInverse : public gsDynamicMassInverse<T>
+{
+public:
+    void compute(const gsSparseMatrix<T> & mass)
+    {
+        gsSparseMatrix<T> eye(mass.rows(), mass.cols());
+        eye.setIdentity();
+        typename gsSparseSolver<T>::LU solver(mass);
+        GISMO_ENSURE(solver.info() == gsEigen::Success,
+                     "LU factorization of mass matrix failed.");
+        gsMatrix<T> inverse = solver.solve(eye);
+        GISMO_ENSURE(solver.info() == gsEigen::Success,
+                     "Solving for inverse mass matrix failed.");
+        m_massInv = inverse.sparseView();
+    }
+
+    void apply(const gsVector<T> & rhs, gsVector<T> & result) const
+    {
+        result = m_massInv * rhs;
+    }
+
+private:
+    gsSparseMatrix<T> m_massInv;
+};
+
+/// @brief Inverse mass object backed by a sparse factorization.
+template <class T>
+class gsDynamicSparseMassInverse : public gsDynamicMassInverse<T>
+{
+public:
+    explicit gsDynamicSparseMassInverse(const std::string & solver = "SimplicialLDLT")
+    :
+    m_solver(gsSparseSolver<T>::get(solver))
+    {
+        GISMO_ENSURE(m_solver, "Unknown sparse solver: "<<solver);
+    }
+
+    void compute(const gsSparseMatrix<T> & mass) override
+    {
+        m_solver->compute(mass);
+        GISMO_ENSURE(m_solver->info()==gsEigen::ComputationInfo::Success, "Mass matrix factorization failed.");
+    }
+
+    void apply(const gsVector<T> & rhs, gsVector<T> & result) const override
+    {
+        result = m_solver->solve(rhs);
+    }
+
+private:
+    typename gsSparseSolver<T>::uPtr m_solver;
+};
+
+} // namespace gismo