DirichletBC.inc.hpp

#pragma once

#include <set>
#include <type_traits>

#include <dune/functions/functionspacebases/subentitydofs.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>

#include <amdis/LinearAlgebra.hpp>
#include <amdis/functions/EntitySet.hpp>
#include <amdis/interpolators/SimpleInterpolator.hpp>
#include <amdis/typetree/TreePath.hpp>

namespace AMDiS {

template <class B, class Row, class Col, class V>
void DirichletBC<B, Row, Col, V>::
init()
{
  if (row_ == col_) {
    auto rowBasis = Dune::Functions::subspaceBasis(*basis_, col_);

    std::set<typename B::MultiIndex> rowSet;

    auto rowLocalView = rowBasis.localView();
    auto rowDOFs = Dune::Functions::subEntityDOFs(rowBasis);
    auto const& gridView = basis_->gridView();
    for (auto const& element : entitySet(*basis_)) {
      if (element.hasBoundaryIntersections()) {
        rowLocalView.bind(element);
        for(auto const& intersection: intersections(gridView, element)) {
          if (intersection.boundary() && boundarySubset_(intersection)) {
            for(auto localIndex: rowDOFs.bind(rowLocalView,intersection))
              rowSet.insert(rowLocalView.index(localIndex));
          }
        }
      }
    }

    rowIndices_.clear(); rowIndices_.reserve(rowSet.size());
    rowIndices_.insert(rowIndices_.begin(), rowSet.begin(), rowSet.end());
    colIndices_ = rowIndices_;
  } else {
    auto rowBasis = Dune::Functions::subspaceBasis(*basis_, row_);
    auto colBasis = Dune::Functions::subspaceBasis(*basis_, col_);

    std::set<typename B::MultiIndex> rowSet, colSet;

    auto rowLocalView = rowBasis.localView();
    auto colLocalView = colBasis.localView();
    auto rowDOFs = Dune::Functions::subEntityDOFs(rowBasis);
    auto colDOFs = Dune::Functions::subEntityDOFs(colBasis);
    auto const& gridView = basis_->gridView();
    for (auto const& element : entitySet(*basis_)) {
      if (element.hasBoundaryIntersections()) {
        rowLocalView.bind(element);
        colLocalView.bind(element);
        for(auto const& intersection: intersections(gridView, element)) {
          if (intersection.boundary() && boundarySubset_(intersection)) {
            for(auto localIndex: rowDOFs.bind(rowLocalView,intersection))
              rowSet.insert(rowLocalView.index(localIndex));
            for(auto localIndex: colDOFs.bind(colLocalView,intersection))
              colSet.insert(colLocalView.index(localIndex));
          }
        }
      }
    }

    rowIndices_.clear(); rowIndices_.reserve(rowSet.size());
    rowIndices_.insert(rowIndices_.begin(), rowSet.begin(), rowSet.end());

    colIndices_.clear(); colIndices_.reserve(colSet.size());
    colIndices_.insert(colIndices_.begin(), colSet.begin(), colSet.end());
  }
}


template <class B, class Row, class Col, class V>
  template <class Mat, class Sol, class Rhs>
void DirichletBC<B, Row, Col, V>::
apply(Mat& matrix, Sol& solution, Rhs& rhs, bool symmetric)
{
  std::vector<typename Sol::value_type> solutionValues;
  solutionValues.reserve(colIndices_.size());
  {
    // create a copy the solution, because the valueGridFct_ might involve the solution vector.
    Sol boundarySolution{solution};
    valueOf(boundarySolution, col_).interpolate_noalias(valueGridFct_, tag::assign{});
    boundarySolution.gather(colIndices_,solutionValues);
  }

  // copy boundary solution dirichlet data to solution vector
  solution.scatter(colIndices_,solutionValues);
  solution.finish();

  if (symmetric)
    matrix.zeroRowsColumns(rowIndices_, colIndices_, true, solution, rhs);
  else
    matrix.zeroRows(rowIndices_, colIndices_, true);

  // copy boundary solution dirichlet data to rhs vector
  rhs.scatter(rowIndices_,solutionValues);
  rhs.finish();
}

} // end namespace AMDiS