dunefunctionsboundaryfunctionalassembler.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSBOUNDARYFUNCTIONALASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSBOUNDARYFUNCTIONALASSEMBLER_HH

#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>

#include <dune/matrix-vector/axpy.hh>

#include "dune/fufem/functionspacebases/functionspacebasis.hh"
#include "dune/fufem/boundarypatch.hh"


namespace Dune {
namespace Fufem {


//! Generic global assembler for functionals on a boundary
template <class TrialBasis>
class DuneFunctionsBoundaryFunctionalAssembler
{
  using GridView = typename TrialBasis::GridView;

public:
  //! create assembler for grid
  DuneFunctionsBoundaryFunctionalAssembler(const TrialBasis& tBasis, const BoundaryPatch<GridView>& boundaryPatch) :
      trialBasis_(tBasis),
      boundaryPatch_(boundaryPatch)
  {}



  template <class VectorBackend, class LocalAssembler>
  void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const
  {
    auto trialLocalView     = trialBasis_.localView();

    using Field = std::decay_t<decltype(vectorBackend(trialLocalView.index(0)))>;
    using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>;

    auto localVector = LocalVector(trialLocalView.maxSize());

    // iterate over boundary
    for (auto it =boundaryPatch_.begin(); it!= boundaryPatch_.end(); ++it)
    {
      const auto& element = it->inside();
      trialLocalView.bind(element);

      localAssembler(it, localVector, trialLocalView);

      for (size_t localRow=0; localRow<trialLocalView.size(); ++localRow)
      {
        auto row = trialLocalView.index(localRow);
        vectorBackend(row) += localVector[localRow];
      }
    }
  }

  template <class VectorBackend, class LocalAssembler>
  void assembleBulk(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const
  {
    vectorBackend.resize(trialBasis_);
    vectorBackend.vector() = 0.0;
    assembleBulkEntries(vectorBackend, std::forward<LocalAssembler>(localAssembler));
  }


protected:
  const TrialBasis& trialBasis_;
  const BoundaryPatch<GridView>& boundaryPatch_;
};


/**
 * \brief Create DuneFunctionsBoundaryFunctionalAssembler
 */
template <class TrialBasis, class BP>
auto duneFunctionsBoundaryFunctionalAssembler(const TrialBasis& trialBasis, const BP& bp)
{
  return DuneFunctionsBoundaryFunctionalAssembler<TrialBasis>(trialBasis, bp);
}



} // namespace Fufem
} // namespace Dune


#endif