Constraints.hpp

#pragma once

#include <amdis/Output.hpp>
#include <amdis/linearalgebra/Constraints.hpp>
#include <amdis/linearalgebra/istl/MatrixBackend.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp>

namespace AMDiS
{
  template <class T, class C>
  struct Constraints<ISTLBCRSMatrix<T,C>>
  {
    using Matrix = ISTLBCRSMatrix<T,C>;
    using Vector = ISTLBlockVector<T>;

    template <class BitVector>
    static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
    {
      // loop over the matrix rows
      for (std::size_t i = 0; i < mat.matrix().N(); ++i) {
        if (nodes[i]) {
          auto cIt = mat.matrix()[i].begin();
          auto cEndIt = mat.matrix()[i].end();
          // loop over nonzero matrix entries in current row
          for (; cIt != cEndIt; ++cIt)
            *cIt = (setDiagonal && i == cIt.index() ? T(1) : T(0));
        }
      }

      // copy solution dirichlet data to rhs vector
      for (std::size_t 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");
    }
  };

} // end namespace AMDiS