Constraints.hpp

#pragma once

#include <list>

#include <Eigen/SparseCore>

#include <amdis/common/Index.hpp>
#include <amdis/linearalgebra/Constraints.hpp>
#include <amdis/linearalgebra/eigen/MatrixBackend.hpp>
#include <amdis/linearalgebra/eigen/VectorBackend.hpp>

namespace AMDiS
{
  template <class T, int O>
  struct Constraints<EigenSparseMatrix<T,O>>
  {
    using Matrix = EigenSparseMatrix<T,O>;
    using Vector = EigenVector<T>;

    template <class BitVector>
    static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
    {
      clearDirichletRow(mat.matrix(), nodes, setDiagonal);

      // copy solution dirichlet data to rhs vector
      for (typename Vector::size_type i = 0; i < sol.vector().size(); ++i) {
        if (nodes[i])
          rhs.vector()[i] = sol.vector()[i];
      }
    }

    template <class BitVector, class Associations>
    static void periodicBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& left, Associations const& left2right,
      bool setDiagonal = true)
    {
      error_exit("Not implemented");
    }

  protected:
    template <class BitVector>
    static void clearDirichletRow(Eigen::SparseMatrix<T, Eigen::ColMajor>& mat, BitVector const& nodes, bool setDiagonal)
    {
      using Mat = Eigen::SparseMatrix<T, Eigen::ColMajor>;
      for (typename Mat::Index c = 0; c < mat.outerSize(); ++c) {
        for (typename Mat::InnerIterator it(mat, c); it; ++it) {
          if (nodes[it.row()]) {
            it.valueRef() = (setDiagonal && it.row() == it.col() ? T(1) : T(0));
            break;
          }
        }
      }
    }

    template <class BitVector>
    static void clearDirichletRow(Eigen::SparseMatrix<T, Eigen::RowMajor>& mat, BitVector const& nodes, bool setDiagonal)
    {
      using Mat = Eigen::SparseMatrix<T, Eigen::RowMajor>;
      for (typename Mat::Index r = 0; r < mat.outerSize(); ++r) {
        if (nodes[r]) {
          for (typename Mat::InnerIterator it(mat, r); it; ++it) {
            it.valueRef() = (setDiagonal && it.row() == it.col() ? T(1) : T(0));
          }
        }
      }
    }
  };

} // end namespace AMDiS