ENH: Add update state methods to time integrators

Add update state methods to time integrators. Separate classes for different types of time integrators.

Base/TimeIntegrators/imstkBackwardEuler.h

+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+      http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+   =========================================================================*/
+
+#ifndef imstkBackwardEuler_h
+#define imstkBackwardEuler_h
+
+#include "imstkTimeIntegrator.h"
+
+#include <array>
+#include "g3log/g3log.hpp"
+
+namespace imstk
+{
+
+///
+/// \class Backward Euler time integration
+///
+class BackwardEuler : public TimeIntegrator
+{
+
+public:
+    ///
+    /// \brief Constructor
+    ///
+    BackwardEuler(const double dT) : TimeIntegrator(dT), m_type(Type::BackwardEuler){};
+
+    ///
+    /// \brief Destructor
+    ///
+    ~BackwardEuler() = default;
+
+    void updateStateGivenDv(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& dV)
+    {
+        currentState->getQDot() = prevState->getQDot() + dV;
+        currentState->getQ() = prevState->getQ() + m_dT*currentState->getQDot();
+    }
+
+    void updateStateGivenDu(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& dU)
+    {
+        currentState->getQ() = prevState->getQ() + dU;
+        currentState->getQDot() = (currentState->getQ() - prevState->getQ())/m_dT;
+    }
+
+    void updateStateGivenV(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& v)
+    {
+        currentState->getQDot() = v;
+        currentState->getQ() = prevState->getQ() + m_dT*currentState->getQDot();
+    }
+
+    void updateStateGivenU(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& u)
+    {
+        currentState->getQ() = u;
+        currentState->getQDot() = (currentState->getQ() - prevState->getQ()) / m_dT;
+    }
+protected:
+
+    // Coefficients of the time integrator
+    std::array<double, 3> m_alpha = { { 1, 0, 0 } };
+    std::array<double, 3> m_beta = { { 1, -1, 0 } };
+    std::array<double, 3> m_gamma = { { 1, -2, -1 } };
+
+    double m_dT; ///> Delta T
+};
+
+} // imstk
+
+#endif // ifndef imstkBackwardEuler_h

Base/TimeIntegrators/imstkNewmarkBeta.h

+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+      http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+   =========================================================================*/
+
+#ifndef imstkBackwardEuler_h
+#define imstkBackwardEuler_h
+
+#include "imstkTimeIntegrator.h"
+
+#include <array>
+#include "g3log/g3log.hpp"
+
+namespace imstk
+{
+
+///
+/// \class Newmark-beta time integration
+///
+class NewmarkBeta : public TimeIntegrator
+{
+
+public:
+    ///
+    /// \brief Constructor
+    ///
+    NewmarkBeta(const double dT, const double beta = 0.25, const double gamma = 0.5) : TimeIntegrator(dT), m_type(Type::NewmarkBeta), m_gamma(gamma), m_beta(beta){};
+
+    ///
+    /// \brief Destructor
+    ///
+    ~NewmarkBeta() = default;
+
+    void updateStateGivenDv(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& dV)
+    {
+        currentState->getQDot() = prevState->getQDot() + dV;
+        currentState->getQDotDot() = (currentState->getQDot() - prevState->getQDot()) / (m_gamma*m_dT) - (1.0 / m_gamma - 1)*prevState->getQDotDot();
+        currentState->getQ() = prevState->getQ() + m_dT*currentState->getQDot() + 0.5*m_dT*m_dT*((1 - 2 * m_beta)*prevState->getQDotDot() + 2 * m_beta*currentState->getQDotDot());
+    }
+
+    void updateStateGivenDu(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& dU)
+    {
+
+    }
+
+    void updateStateGivenV(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& v)
+    {
+        currentState->getQDot() = v;
+        currentState->getQDotDot() = (currentState->getQDot() - prevState->getQDot()) / (m_gamma*m_dT) - (1.0 / m_gamma - 1)*prevState->getQDotDot();
+        currentState->getQ() = prevState->getQ() + m_dT*currentState->getQDot() + 0.5*m_dT*m_dT*((1 - 2 * m_beta)*prevState->getQDotDot() + 2 * m_beta*currentState->getQDotDot());
+    }
+
+    void updateStateGivenU(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& u)
+    {
+
+    }
+protected:
+    double m_beta;
+    double m_gamma;
+
+    // Coefficients of the time integrator
+    std::array<double, 3> m_alpha = { { 1, 0, 0 } };
+    std::array<double, 3> m_beta = { { 1, -1, 0 } };
+    std::array<double, 3> m_gamma = { { 1, -2, -1 } };
+
+    double m_dT; ///> Delta T
+};
+
+} // imstk
+
+#endif // ifndef imstkBackwardEuler_h

Base/TimeIntegrators/imstkTimeIntegrator.cpp

-/*=========================================================================
-
-Library: iMSTK
-
-Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
-& Imaging in Medicine, Rensselaer Polytechnic Institute.
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-http://www.apache.org/licenses/LICENSE-2.0.txt
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-
-=========================================================================*/
-
-#include "imstkTimeIntegrator.h"
-
-#include <g3log/g3log.hpp>
-
-namespace imstk
-{
-
-TimeIntegrator::TimeIntegrator(const TimeIntegrator::Type type, double dT) : m_type(type), m_dT(dT){}
-
-void
-TimeIntegrator::setType(const TimeIntegrator::Type type)
-{
-    m_type = type;
-    this->setCoefficients(type);
-}
-
-const TimeIntegrator::Type&
-TimeIntegrator::getType() const
-{
-    return m_type;
-}
-
-void
-TimeIntegrator::setCoefficients(const TimeIntegrator::Type type)
-{
-    switch (type)
-    {
-    case TimeIntegrator::Type::BackwardEuler:
-        m_alpha = { { 1, 0, 0 } };
-        m_beta = { { 1, -1, 0 } };
-        m_gamma = { { 1, -2, -1 } };
-        break;
-
-    case TimeIntegrator::Type::ForwardEuler:
-    case TimeIntegrator::Type::NewmarkBeta:
-    case TimeIntegrator::Type::CentralDifference:
-        LOG(WARNING) << "TimeIntegrator::setCoefficients error: type of the time integrator not supported.";
-        break;
-
-    default:
-        LOG(WARNING) << "TimeIntegrator::setCoefficients error: type of the time integrator not identified!";
-        break;
-    }
-}
-
-void
-TimeIntegrator::setTimestepSize(const double dT)
-{
-    m_dT = dT;
-}
-
-}

Base/TimeIntegrators/imstkTimeIntegrator.h

@@ -25,6 +25,10 @@
 #include <array>
 #include "g3log/g3log.hpp"
 
+#include "Eigen/sparse"
+
+#include "imstkProblemState.h"
+
 namespace imstk
 {
 
@@ -38,12 +42,14 @@ namespace imstk
 ///
 class TimeIntegrator
 {
+public:
     enum class Type
     {
         ForwardEuler,
         BackwardEuler,
         NewmarkBeta,
         CentralDifference,
+        NoTimeStepper,
         none
     };
 
@@ -51,7 +57,7 @@ public:
     ///
     /// \brief Constructor
     ///
-    TimeIntegrator(const TimeIntegrator::Type type, double dT);
+    TimeIntegrator(double dT) : m_type(Type::none), m_dT(dT){};
 
     ///
     /// \brief Destructor
@@ -59,33 +65,26 @@ public:
     ~TimeIntegrator() = default;
 
     ///
-    /// \brief Set/Get type of the time integrator
+    /// \brief Return the type of the time integrator
     ///
-    void setType(const TimeIntegrator::Type type);
-    const Type& getType() const;
+    TimeIntegrator::Type getType() const { return m_type; };
 
     ///
-    /// \brief Set coefficients for a given time integrator type
+    /// \brief Set/Get the time step size
     ///
-    void setCoefficients(const TimeIntegrator::Type type);
+    void setTimestepSize(const double dT){ m_dT = dT; }
+    double getTimestepSize() const { return m_dT; }
 
     ///
-    /// \brief Set the time step size
+    /// \brief Update states given the updates in different forms
     ///
-    void setTimestepSize(const double dT);
-    double getTimestepSize() const
-    {
-        return m_dT;
-    }
+    virtual void updateStateGivenDv(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& dV) = 0;
+    virtual void updateStateGivenDu(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& dU) = 0;
+    virtual void updateStateGivenV(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& v) = 0;
+    virtual void updateStateGivenU(std::shared_ptr<ProblemState> prevState, std::shared_ptr<ProblemState> currentState, Vectord& u) = 0;
 
 protected:
-    TimeIntegrator::Type m_type; ///> Type of the time integrator
-
-    // Coefficients of the time integrator
-    std::array<double,3> m_alpha;
-    std::array<double,3> m_gamma;
-    std::array<double,3> m_beta;
-
+    Type m_type; ///> Type of the time integrator
     double m_dT; ///> Delta T
 };