a tool for extracting variables with boundary conditions. currently only have dirichlet there

active_set.cpp

#include "active_set.h"
#include <set>

template <>
const ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Free>() const {
    return m_active_pair;
}
template <>
const ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Dirichlet>() const {
    return m_dirichlet_pair;
}

template <>
ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Free>() {
    return m_active_pair;
}
template <>
ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Dirichlet>() {
    return m_dirichlet_pair;
}


void ActiveSet::updateActive() {
    std::set<Index> constrained;

    auto&& diri = getIndices<ActiveSet::VariableType::Dirichlet>();
    std::copy(diri.begin(),diri.end(),std::inserter(constrained,constrained.end()));
    int prev = 0;



    auto& active_indices = m_active_pair.indices;
    active_indices.clear();
    active_indices.reserve(active_size());

    for(auto&& d: constrained) {
        if(prev < d) 
            for(; prev < d; ++prev) {
                active_indices.push_back(prev);
            }
        ++prev;
    }
    if(prev != m_total_size) {
        for(; prev < m_total_size; ++prev) {
            active_indices.push_back(prev);
        }
    }


    m_active_pair.updateInverse(m_total_size); 


}
void ActiveSet::IndexPair::updateInverse(size_t total_size) {
    inverse = std::vector<Index>(total_size,-1);
    for(int i = 0; i < indices.size(); ++i) {
        inverse[indices[i]] = i;
    }
}

active_set.h

#ifndef ACTIVE_SET_H
#define ACTIVE_SET_H

#include <vector>
#include <Eigen/Dense>
#include <Eigen/Sparse>


class ActiveSet {
    public:
        enum VariableType { Free, Dirichlet };
        using Index = int;
        using IndexList = std::vector<int>;

        struct IndexPair {
            IndexList indices;
            IndexList inverse;
            size_t size() const { return indices.size(); }
            void updateInverse(size_t total_size);
        };


        size_t active_size() const { return m_active_pair.size(); }
        size_t reduced_size() const { return active_size(); }
        template <VariableType VarType>
            size_t var_size() const { return getIndexPair<VarType>().size();}

        void updateActive();

        ActiveSet(int k): m_total_size(k) {}
        template <VariableType VarType>
            void setVariableType(const IndexList& indices);

        template <VariableType VarType=VariableType::Free>
        const IndexPair& getIndexPair() const;
        template <VariableType VarType=VariableType::Free>
        IndexPair& getIndexPair() ;


        template <VariableType VarType=VariableType::Free>
        const IndexList& getIndices() const { return getIndexPair<VarType>().indices; }
        template <VariableType VarType=VariableType::Free>
        const IndexList& getInverseIndices() const { return getIndexPair<VarType>().inverse; }

        template <VariableType VarType=VariableType::Free>
        Index reduce(Index i) const { return getIndices<VarType>()[i]; }
        template <VariableType VarType=VariableType::Free>
        Index expand(Index i) const { return getInverseIndices<VarType>()[i]; }

        template <VariableType VarType, typename Derived, typename RetType=typename Eigen::internal::remove_all<typename Eigen::internal::eval<Derived>::type>::type>
             RetType reduce(const Eigen::MatrixBase<Derived>& M) const;

        template <VariableType VarType, typename Derived, typename RetType=typename Eigen::internal::remove_all<typename Eigen::internal::eval<Derived>::type>::type>
             RetType expand(const Eigen::MatrixBase<Derived>& M) const;


        template <VariableType VarType, typename T>
            Eigen::SparseMatrix<T> reduce(const Eigen::SparseMatrix<T>& A) const;
        template <ActiveSet::VariableType VarType, ActiveSet::VariableType VarType2, typename T>
            Eigen::SparseMatrix<T> reduceMatrix(const Eigen::SparseMatrix<T>& A) const;
    private:
        Index m_total_size;
        IndexPair m_active_pair;
        IndexPair m_dirichlet_pair;
};

template <>
const ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Free>() const;
template <>
const ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Dirichlet>() const;

template <>
ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Free>() ;
template <>
ActiveSet::IndexPair& ActiveSet::getIndexPair<ActiveSet::VariableType::Dirichlet>() ;




template <ActiveSet::VariableType VarType>
void ActiveSet::setVariableType(const IndexList& indices) {
    auto&& pair = getIndexPair<VarType>();
    pair.indices = indices;
    pair.updateInverse(m_total_size);
}

template <ActiveSet::VariableType VarType, typename Derived, typename RetType>
RetType ActiveSet::expand(const Eigen::MatrixBase<Derived>& M) const {
    RetType R = RetType::Zero(m_total_size,M.cols());

    auto&& pair = getIndexPair<VarType>();
    assert(M.rows() == pair.size() );
    for(int i = 0; i < pair.size(); ++i) {

        R.row(pair.indices[i]) = M.row(i);
    }

    return R;
}

template <ActiveSet::VariableType VarType, typename Derived, typename RetType>
RetType ActiveSet::reduce(const Eigen::MatrixBase<Derived>& M) const {
    RetType R = RetType::Zero(active_size(),M.cols());

    assert(M.rows() == m_total_size);

    auto&& pair = getIndexPair<VarType>();
    for(int i = 0; i < pair.size(); ++i) {

        R.row(i) = M.row(pair.indices[i]);
    }

    return R;
}






template <ActiveSet::VariableType VarType, typename T>
Eigen::SparseMatrix<T> ActiveSet::reduce(const Eigen::SparseMatrix<T>& M) const {
    return reduceMatrix<VarType,VarType,T>(M);
}

template <ActiveSet::VariableType VarType, ActiveSet::VariableType VarType2, typename T>
Eigen::SparseMatrix<T> ActiveSet::reduceMatrix(const Eigen::SparseMatrix<T>& M) const {
    auto&& pair = getIndexPair<VarType>();
    auto&& pair2 = getIndexPair<VarType2>();


    assert(M.rows() == m_total_size);
    assert(M.cols() == m_total_size);

    Eigen::SparseMatrix<T> A(pair.size(), pair2.size());
    std::vector<Eigen::Triplet<T>> trips;
    for(int i = 0; i < pair2.size(); ++i) {

        int k = pair2.indices[i];

        for(typename Eigen::SparseMatrix<T>::InnerIterator it(M,k); it; ++it) {
            if(pair.inverse[it.row()] != -1) {
                int r = pair.inverse[it.row()];
                int c = pair2.inverse[it.col()];
                trips.emplace_back(r,c,it.value());
            }
        }
    }

    A.setFromTriplets(trips.begin(),trips.end());
    return A;
}


#endif//ACTIVE_SET_H

test.cpp

#include "active_set.h"
#include <iostream>

int main(int argc, char* argv[]) {
    std::vector<int> bad;
    int full_size = 20;
    std::vector<int> dir = {0,3,4,10};
  ActiveSet coords(20);

  coords.setVariableType<ActiveSet::VariableType::Dirichlet>(dir);
  coords.updateActive();
  for(auto&& i: coords.getIndices()) {
      std::cout <<i << " ";
  }
  std::cout << std::endl;


  Eigen::SparseMatrix<double> M(full_size,full_size);
    std::vector<Eigen::Triplet<double>> trips;
    for(int i = 0; i < full_size; ++i) {
        trips.emplace_back(i,i,i+1);
    }
    for(int i = 0; i < full_size-1; ++i) {
        trips.emplace_back(i,i+1,-i-1);
        trips.emplace_back(i+1,i,-i-1);
    }

    M.setFromTriplets(trips.begin(),trips.end());
    std::cout << M << std::endl;

    std::cout << "Reduced: " << std::endl;

    std::cout << coords.reduce<ActiveSet::VariableType::Dirichlet>(M) << std::endl;
    std::cout << coords.reduce<ActiveSet::VariableType::Free>(M) << std::endl;
    std::cout << coords.reduceMatrix<ActiveSet::VariableType::Dirichlet,ActiveSet::VariableType::Free>(M) << std::endl;
    std::cout << coords.reduceMatrix<ActiveSet::VariableType::Free,ActiveSet::VariableType::Dirichlet>(M) << std::endl;


}