#pragma once

#include <dune/functions/functionspacebases/interpolate.hh>

#include "Log.hpp"

namespace AMDiS
{
  template <class WorldVector>
    template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
  void DirichletBC<WorldVector>::init(bool apply, 
				      RowFeSpace const& rowFeSpace, 
				      ColFeSpace const& colFeSpace, 
				      DOFMatrix* matrix, 
				      DOFVector* rhs, 
				      DOFVector* solution)
  {
    using Dune::Functions::interpolate;
    
    if (!initialized) {
      interpolate(rowFeSpace, dirichletNodes, predicate);
      initialized = true;
    }
  }
    
    
  template <class WorldVector>
    template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
  void DirichletBC<WorldVector>::finish(bool apply, 
					RowFeSpace const& rowFeSpace, 
					ColFeSpace const& colFeSpace, 
					DOFMatrix* matrix, 
					DOFVector* rhs, 
					DOFVector* solution)
  {
    using Dune::Functions::interpolate;
    
    AMDIS_TEST_EXIT( initialized, "Boundary condition not initialized!" );
    
    // loop over the matrix rows
    for (size_t i = 0; i < matrix->N(); ++i) {
	if (dirichletNodes[i]) {
	    auto cIt = (*matrix)[i].begin();
	    auto cEndIt = (*matrix)[i].end();
	    // loop over nonzero matrix entries in current row
	    for (; cIt != cEndIt; ++cIt)
		*cIt = (apply && i == cIt.index()) ? 1.0 : 0.0;
	}
    }

    if (apply) {
	interpolate(rowFeSpace, *rhs, values, dirichletNodes);
	interpolate(colFeSpace, *solution, values, dirichletNodes);
    }
  }

} // end namespace AMDiS