imstkReducedStVKBodyModel.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.

=========================================================================*/

#pragma once

#include <memory>

#include "imstkDynamicalModel.h"
#include "imstkVectorizedState.h"
#include "imstkNonLinearSystem.h"
#include "imstkLogger.h"

namespace vega
{
class VolumetricMesh;
class ReducedStVKForceModel;
class StVKReducedInternalForces;
class ModalMatrix;
} // namespace vega

namespace imstk
{
class InternalForceModel;
class SolverBase;
class TaskNode;
class TimeIntegrator;
class VegaMeshIO;

struct ReducedStVKConfig
{
    std::string m_cubicPolynomialFilename;
    std::string m_modesFileName;
    // int         r;

    double m_dampingMassCoefficient      = 0.1;
    double m_dampingStiffnessCoefficient = 0.01;
    double m_dampingLaplacianCoefficient = 0.0;
    double m_deformationCompliance       = 1.0;
    double m_gravity = 9.81;

    int m_numberOfThreads = 4;
};

class ReducedStVK : public DynamicalModel<FeDeformBodyState>
{
using kinematicState = FeDeformBodyState;
// using Matrixd = Eigen::MatrixXd;
public:
    ReducedStVK();

    virtual ~ReducedStVK() override;

public:
    void configure(const std::string& configFileName);
    void configure(std::shared_ptr<ReducedStVKConfig> config = std::make_shared<ReducedStVKConfig>());

    bool initialize() override;

    void setForceModelConfiguration(std::shared_ptr<ReducedStVKConfig> config);
    std::shared_ptr<ReducedStVKConfig> getForceModelConfiguration() const;

    void setInternalForceModel(std::shared_ptr<vega::StVKReducedInternalForces> fm);
    std::shared_ptr<imstk::InternalForceModel> getInternalForceModel() const;

    void setTimeIntegrator(std::shared_ptr<TimeIntegrator> timeIntegrator);
    std::shared_ptr<TimeIntegrator> getTimeIntegrator() const;

    void loadInitialStates();

    bool initializeForceModel();

    bool initializeMassMatrix();

    bool initializeDampingMatrix();

    bool initializeTangentStiffness();

    bool initializeGravityForce();

    bool initializeExplicitExternalForces();

    static void initializeEigenMatrixFromStdVector(const std::vector<double>& v,
                                                   Matrixd&                   eigenMatrix);

    void computeImplicitSystemRHS(kinematicState&       prevState,
                                  kinematicState&       newState,
                                  const StateUpdateType updateType);

    void computeSemiImplicitSystemRHS(kinematicState&       stateAtT,
                                      kinematicState&       newState,
                                      const StateUpdateType updateType);

    void computeImplicitSystemLHS(const kinematicState& prevState,
                                  kinematicState&       newState,
                                  const StateUpdateType updateType);

    void updateDampingMatrix();

    void updateMassMatrix();

    void applyBoundaryConditions(Matrixd& M, const bool withCompliance = false) const;
    void applyBoundaryConditions(Vectord& x) const;

    void updatePhysicsGeometry() override;

    void updateBodyStates(const Vectord& solution, const StateUpdateType updateType);
    void updateBodyIntermediateStates(const Vectord& solution, const StateUpdateType updateType);

    void updateBodyPreviousStates();

    NonLinearSystem<Matrixd>::VectorFunctionType getFunction();

    NonLinearSystem<Matrixd>::UpdateFunctionType          getUpdateFunction();
    NonLinearSystem<Matrixd>::UpdatePrevStateFunctionType getUpdatePrevStateFunction();

    NonLinearSystem<Matrixd>::MatrixFunctionType getFunctionGradient();

    Vectord& getContactForce();

    Vectord& getUnknownVec()
    {
        return m_qSol;
    }

    void setUpdateType(const StateUpdateType& updateType)
    {
        m_updateType = updateType;
    }

    const StateUpdateType& getUpdateType() const
    {
        return m_updateType;
    }

    std::vector<std::size_t>& getFixNodeIds()
    {
        return m_fixedNodeIds;
    }

    virtual void setTimeStep(const double timeStep);

    virtual double getTimeStep() const;

    void setFixedSizeTimeStepping();

    void enableFixedBC()
    {
        LOG(WARNING) << "Boundary conditions are not allowed to change";
        m_implementFixedBC = true;
    };
    void disableFixedBC()
    {
        // LOG(FATAL, "Boundary conditions are not allowed to change");
        m_implementFixedBC = false;
    };
    bool isFixedBCImplemented() const
    {
        return m_implementFixedBC;
    };

    void prolongate(const Vectord& uReduced, Vectord& u) const;

    void prolongate(kinematicState& uReduced, kinematicState& u) const;

    void project(const Vectord& u, Vectord& uReduced) const;
    void readModalMatrix(const std::string& fname);

    std::shared_ptr<TaskNode> getSolveNode() const { return m_solveNode; }

    std::shared_ptr<SolverBase> getSolver() const { return m_solver; }
    void setSolver(std::shared_ptr<SolverBase> solver) { this->m_solver = solver; }

protected:
    void initGraphEdges(std::shared_ptr<TaskNode> source, std::shared_ptr<TaskNode> sink) override;

protected:
    std::shared_ptr<SolverBase> m_solver = nullptr;
    std::shared_ptr<vega::StVKReducedInternalForces> m_internalForceModel; 
    std::shared_ptr<vega::ReducedStVKForceModel>     m_forceModel;
    std::shared_ptr<TimeIntegrator> m_timeIntegrator;                      
    std::shared_ptr<NonLinearSystem<Matrixd>> m_nonLinearSystem;           
    std::shared_ptr<vega::ModalMatrix> m_modalMatrix;

    std::shared_ptr<ReducedStVKConfig> m_config;

    Matrixd m_M;    
    Matrixd m_C;    
    Matrixd m_K;    
    Matrixd m_Keff; 

    // full-space vectors
    Vectord m_Fcontact;          
    Vectord m_Fgravity;          
    Vectord m_FexplicitExternal; 
    Vectord m_qSol;              

    // reduced-space vectors
    Vectord m_Feff;                     
    Vectord m_Finternal;                
    Vectord m_qSolReduced;              
    Vectord m_FcontactReduced;
    Vectord m_FgravityReduced;          
    Vectord m_FexplicitExternalReduced; 
    size_t  m_numDOFReduced;

    std::shared_ptr<vega::VolumetricMesh> m_vegaPhysicsMesh;       
    std::vector<double> m_massMatrix;                              
    std::vector<double> m_stiffnessMatrix;                         
    std::vector<double> m_dampingMatrix;                           

    std::vector<std::size_t> m_fixedNodeIds;                       

    StateUpdateType m_updateType = StateUpdateType::DeltaVelocity; 

    bool m_damped = false;                                         

    // If this is true, the tangent stiffness and force vector will be modified to
    // accommodate (the rows and columns will be nullified) the fixed boundary conditions
    bool m_implementFixedBC = false;

    std::shared_ptr<kinematicState> m_initialStateReduced;
    std::shared_ptr<kinematicState> m_previousStateReduced;
    std::shared_ptr<kinematicState> m_currentStateReduced;

private:
    std::shared_ptr<TaskNode> m_solveNode = nullptr;
};
}  // namespace imstk