local_static_condensation.hpp

Go to the documentation of this file.

// Structure to store information for and perform local static condensation
 
#ifndef LOCAL_STATIC_CONDENSATION_HPP
#define LOCAL_STATIC_CONDENSATION_HPP
 
#include <Eigen/Dense>
 
namespace HArDCore2D
{
 
 
  struct LocalStaticCondensation
  {
    LocalStaticCondensation(
                    const Eigen::MatrixXd & Perm,         
                    const std::vector<size_t> & globalDOFs_gl, 
                    const std::vector<size_t> & globalDOFs_sc  
                    )
         : m_Perm(Perm),
           m_globalDOFs_gl(globalDOFs_gl),
           m_globalDOFs_sc(globalDOFs_sc),
           m_dim_gl(globalDOFs_gl.size()),
           m_dim_sc(globalDOFs_sc.size())
           {
            // do nothing
           };
 
 
    std::tuple<Eigen::MatrixXd, Eigen::VectorXd, Eigen::MatrixXd, Eigen::VectorXd>
      compute(const std::pair<Eigen::MatrixXd, Eigen::VectorXd> & lsT)
    {
      // We permute the DOFs that are statically condensed
      Eigen::MatrixXd AT = m_Perm * lsT.first * m_Perm.transpose();
      Eigen::VectorXd bT = m_Perm * lsT.second;
  
      Eigen::MatrixXd AT_gl, AT_sc;
      Eigen::VectorXd bT_gl, bT_sc;
  
      if (m_dim_sc > 0){
        // Extract 4 blocks of AT, bT for static condensation
        Eigen::MatrixXd A11 = AT.topLeftCorner(m_dim_gl, m_dim_gl);  
        Eigen::MatrixXd A12 = AT.topRightCorner(m_dim_gl, m_dim_sc);
        Eigen::MatrixXd A21 = AT.bottomLeftCorner(m_dim_sc, m_dim_gl);
        Eigen::MatrixXd A22 = AT.bottomRightCorner(m_dim_sc, m_dim_sc);
        Eigen::VectorXd b1 = bT.head(m_dim_gl);
        Eigen::VectorXd b2 = bT.tail(m_dim_sc);
        // Create condensed system (AT_reduced, bT_reduced) and SC recovery operator (RT, cT)
        // Eigen::PartialPivLU<Eigen::MatrixXd> A22pivlu(A22);
        Eigen::FullPivLU<Eigen::MatrixXd> A22pivlu(A22);
        if( !A22pivlu.isInvertible() ) {
          std::cout << "[LocalStaticCondensation] Found non invertible local matrix" << std::endl;
          exit(1);
        } // if
        Eigen::MatrixXd A22inv_A21 = A22pivlu.solve(A21);
        Eigen::VectorXd A22inv_b2 = A22pivlu.solve(b2);
        AT_gl = A11 - A12 * A22inv_A21;
        bT_gl = b1 - A12 * A22inv_b2;
        AT_sc = -A22inv_A21;
        bT_sc = A22inv_b2;
      }else{
        // No static condensation
        AT_gl = AT;
        bT_gl = bT;
        AT_sc = Eigen::MatrixXd::Zero(0,0);
        bT_sc = Eigen::VectorXd::Zero(0);
      }
      
      return std::make_tuple(AT_gl, bT_gl, AT_sc, bT_sc);
    };
 
    inline std::vector<size_t> globalDOFs_gl() { return m_globalDOFs_gl;};
    inline size_t dim_gl() { return m_dim_gl;};
    
    inline std::vector<size_t> globalDOFs_sys() { return m_globalDOFs_gl;};
    inline size_t dim_sys() { return m_dim_gl;};
    
    
    inline std::vector<size_t> globalDOFs_sc() { return m_globalDOFs_sc;};
    inline size_t dim_sc() { return m_dim_sc;};
 
    // Members
    Eigen::MatrixXd m_Perm;     
    std::vector<size_t> m_globalDOFs_gl;
    std::vector<size_t> m_globalDOFs_sc;
    size_t m_dim_gl;
    size_t m_dim_sc;
  };
 
 
 
} // end of namespace HArDCore2D
 
#endif