dunefunctionsoperatorassembler.hh

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

#include <type_traits>

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

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

#include <dune/grid/common/mcmgmapper.hh>

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


namespace Dune {
namespace Fufem {


//! Generic global assembler for operators on a gridview
template <class TrialBasis, class AnsatzBasis>
class DuneFunctionsOperatorAssembler
{
  using GridView = typename TrialBasis::GridView;

public:
  //! create assembler for grid
  DuneFunctionsOperatorAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis) :
      trialBasis_(tBasis),
      ansatzBasis_(aBasis)
  {}



  template <class PatternBuilder>
  void assembleBulkPattern(PatternBuilder& patternBuilder) const
  {
    patternBuilder.resize(trialBasis_, ansatzBasis_);

    // create two localViews but use only one if bases are the same
    auto ansatzLocalView = ansatzBasis_.localView();
    auto separateTrialLocalView = trialBasis_.localView();
    auto& trialLocalView = selectTrialLocalView(ansatzLocalView, separateTrialLocalView);

    for (const auto& element : elements(trialBasis_.gridView()))
    {
      // bind the localViews to the element
      bind(ansatzLocalView, trialLocalView, element);

      // Add element stiffness matrix onto the global stiffness matrix
      for (size_t i=0; i<trialLocalView.size(); ++i)
      {
        auto row = trialLocalView.index(i);
        for (size_t j=0; j<ansatzLocalView.size(); ++j)
        {
          auto col = ansatzLocalView.index(j);
          patternBuilder.insertEntry(row,col);
        }
      }
    }
  }



  template <class PatternBuilder>
  void assembleSkeletonPattern(PatternBuilder& patternBuilder) const
  {
    patternBuilder.resize(trialBasis_, ansatzBasis_);

    auto insideTrialLocalView       = trialBasis_.localView();

    auto insideAnsatzLocalView      = ansatzBasis_.localView();

    auto outsideAnsatzLocalView     = ansatzBasis_.localView();

    for (const auto& element : elements(trialBasis_.gridView()))
    {
      insideTrialLocalView.bind(element);

      insideAnsatzLocalView.bind(element);

      /* Coupling on the same element */
      for (size_t i = 0; i < insideTrialLocalView.size(); i++) {
        auto row = insideTrialLocalView.index(i);
        for (size_t j = 0; j < insideAnsatzLocalView.size(); j++) {
          auto col = insideAnsatzLocalView.index(j);
          patternBuilder.insertEntry(row,col);
        }
      }

      for (const auto& is : intersections(trialBasis_.gridView(), element))
      {
        /* Elements on the boundary have been treated above, iterate along the inner edges */
        if (is.neighbor())
        {
          // Current outside element
          const auto& outsideElement = is.outside();

          // Bind the outer parts to the outer element
          outsideAnsatzLocalView.bind(outsideElement);

          // We assume that all basis functions of the inner element couple with all basis functions from the outer one
          for (size_t i=0; i<insideTrialLocalView.size(); ++i)
          {
            auto row = insideTrialLocalView.index(i);
            for (size_t j=0; j<outsideAnsatzLocalView.size(); ++j)
            {
              auto col = outsideAnsatzLocalView.index(j);
              patternBuilder.insertEntry(row,col);
            }
          }
        }
      }
    }
  }



  template <class MatrixBackend, class LocalAssembler>
  void assembleBulkEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const
  {
    // create two localViews but use only one if bases are the same
    auto ansatzLocalView = ansatzBasis_.localView();
    auto separateTrialLocalView = trialBasis_.localView();
    auto& trialLocalView = selectTrialLocalView(ansatzLocalView, separateTrialLocalView);

    using Field = std::decay_t<decltype(matrixBackend(trialLocalView.index(0), ansatzLocalView.index(0)))>;
    using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;

    auto localMatrix = LocalMatrix();

    for (const auto& element : elements(trialBasis_.gridView()))
    {
      // bind the localViews to the element
      bind(ansatzLocalView, trialLocalView, element);

      localMatrix.setSize(trialLocalView.size(), ansatzLocalView.size());
      localAssembler(element, localMatrix, trialLocalView, ansatzLocalView);

      // Add element stiffness matrix onto the global stiffness matrix
      for (size_t localRow=0; localRow<trialLocalView.size(); ++localRow)
      {
        auto row = trialLocalView.index(localRow);
        for (size_t localCol=0; localCol<ansatzLocalView.size(); ++localCol)
        {
          auto col = ansatzLocalView.index(localCol);
          matrixBackend(row,col) += localMatrix[localRow][localCol];
        }
      }
    }
  }


  /* This variant of the IntersectionOperatorAssembler that works
   * with dune-functions bases.
   *
   * Currently not supported are constrained bases and lumping.
   *
   * Note that this was written (and tested) having discontinuous Galerkin bases in mind.
   * There may be cases where things do not work as expected for standard (continuous) Galerkin spaces.
   */
  template <class MatrixBackend, class LocalAssembler, class LocalBoundaryAssembler>
  void assembleSkeletonEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler, LocalBoundaryAssembler&& localBoundaryAssembler) const
  {
    Dune::MultipleCodimMultipleGeomTypeMapper<GridView> faceMapper(trialBasis_.gridView(), mcmgElementLayout());

    auto insideTrialLocalView       = trialBasis_.localView();

    auto insideAnsatzLocalView      = ansatzBasis_.localView();

    auto outsideTrialLocalView      = trialBasis_.localView();

    auto outsideAnsatzLocalView     = ansatzBasis_.localView();

    using Field = std::decay_t<decltype(matrixBackend(insideTrialLocalView.index(0), insideAnsatzLocalView.index(0)))>;
    using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;
    using MatrixContainer = Dune::Matrix<LocalMatrix>;
    auto matrixContrainer = MatrixContainer(2,2);

    for (const auto& element : elements(trialBasis_.gridView()))
    {
      insideTrialLocalView.bind(element);

      insideAnsatzLocalView.bind(element);

      // Resize
      matrixContrainer[0][0].setSize(insideTrialLocalView.size(), insideAnsatzLocalView.size());

      for (const auto& is : intersections(trialBasis_.gridView(), element))
      {

        /* Assemble (depending on whether we have a boundary edge or not) */
        if (is.boundary())
        {
          auto& localMatrix = matrixContrainer[0][0];
          localBoundaryAssembler(is, localMatrix, insideTrialLocalView, insideAnsatzLocalView);

          for (size_t i=0; i< insideTrialLocalView.size(); i++) {
            for (size_t j=0; j< insideAnsatzLocalView.size(); j++)
              matrixBackend(insideTrialLocalView.index(i), insideAnsatzLocalView.index(j))+=localMatrix[i][j];
          }
        }
        else if (is.neighbor())
        {
          /* Only handle every intersection once; we choose the case where the inner element has the lower index
           *
           * This should also work with grids where elements can have different geometries. TODO: Actually test this somewhere!*/
          const auto& outsideElement = is.outside();
          if (faceMapper.index(element) > faceMapper.index(outsideElement))
            continue;


          // Bind the outer parts to the outer element
          outsideTrialLocalView.bind(outsideElement);

          outsideAnsatzLocalView.bind(outsideElement);

          matrixContrainer[0][1].setSize(insideTrialLocalView.size(), outsideAnsatzLocalView.size());
          matrixContrainer[1][0].setSize(outsideTrialLocalView.size(), insideAnsatzLocalView.size());
          matrixContrainer[1][1].setSize(outsideTrialLocalView.size(), outsideAnsatzLocalView.size());

          localAssembler(is, matrixContrainer, insideTrialLocalView, insideAnsatzLocalView, outsideTrialLocalView, outsideAnsatzLocalView);

          for (size_t i=0; i < insideTrialLocalView.size(); i++)
          {
            auto row = insideTrialLocalView.index(i);
            // inside x inside
            for (size_t j=0; j < insideTrialLocalView.size(); j++)
            {
                auto col = insideTrialLocalView.index(j);
                matrixBackend(row, col) += matrixContrainer[0][0][i][j];
            }
            // inside x outside
            for (size_t j=0; j < outsideAnsatzLocalView.size(); j++)
            {
              auto col = outsideAnsatzLocalView.index(j);
              matrixBackend(row, col) += matrixContrainer[0][1][i][j];
            }
          }

          for (size_t i=0; i < outsideTrialLocalView.size(); i++)
          {
            auto row = outsideTrialLocalView.index(i);
            // outside x inside
            for (size_t j=0; j < insideAnsatzLocalView.size(); j++)
            {
              auto col = insideAnsatzLocalView.index(j);
              matrixBackend(row, col) += matrixContrainer[1][0][i][j];
            }
            // outside x outside
            for (size_t j=0; j < outsideAnsatzLocalView.size(); j++)
            {
              auto col = outsideAnsatzLocalView.index(j);
              matrixBackend(row, col) += matrixContrainer[1][1][i][j];
            }
          }
        }
      }
    }
  }


  template <class MatrixBackend, class LocalAssembler>
  void assembleBulk(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const
  {
    auto patternBuilder = matrixBackend.patternBuilder();
    assembleBulkPattern(patternBuilder);
    patternBuilder.setupMatrix();
    matrixBackend.assign(0);
    assembleBulkEntries(matrixBackend, std::forward<LocalAssembler>(localAssembler));
  }

private:

//! helper function to select a trialLocalView (possibly a reference to ansatzLocalView if bases are the same)
template<class AnsatzLocalView, class TrialLocalView>
TrialLocalView& selectTrialLocalView(AnsatzLocalView& ansatzLocalView, TrialLocalView& trialLocalView) const
{
   if constexpr (std::is_same<TrialBasis,AnsatzBasis>::value)
     if (&trialBasis_ == &ansatzBasis_)
       return ansatzLocalView;
   return trialLocalView;
}

//! small helper that checks whether two localViews are the same and binds one or both to an element
template<class AnsatzLocalView, class TrialLocalView, class E>
void bind(AnsatzLocalView& ansatzLocalView, TrialLocalView& trialLocalView, const E& e) const
{
  ansatzLocalView.bind(e);
  if constexpr (std::is_same<TrialBasis,AnsatzBasis>::value)
    if (&trialLocalView == &ansatzLocalView)
      return;
  // localViews differ: bind trialLocalView too
  trialLocalView.bind(e);
}

protected:
  const TrialBasis& trialBasis_;
  const AnsatzBasis& ansatzBasis_;
};


/**
 * \brief Create DuneFunctionsOperatorAssembler
 */
template <class TrialBasis, class AnsatzBasis>
auto duneFunctionsOperatorAssembler(const TrialBasis& trialBasis, const AnsatzBasis& ansatzBasis)
{
  return DuneFunctionsOperatorAssembler<TrialBasis, AnsatzBasis>(trialBasis, ansatzBasis);
}



} // namespace Fufem
} // namespace Dune


#endif // DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH