PeriodicBC.hpp

#pragma once
 
#include <map>
#include <optional>
#include <utility>
#include <vector>
 
#include <dune/common/concept.hh>
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
 
#include <amdis/BoundaryCondition.hpp>
#include <amdis/BoundarySubset.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/TreePath.hpp>
 
 
namespace AMDiS
{
  template <class ct, int dow>
  struct FaceTransformation
  {
    using WorldMatrix = FieldMatrix<ct, dow, dow>;
    using WorldVector = FieldVector<ct, dow>;
 
    WorldMatrix matrix_;
 
    WorldVector shift_;
 
    WorldVector evaluate(WorldVector const& x) const
    {
      WorldVector y = shift_;
      matrix_.umv(x, y);
      return y;
    }
 
    WorldVector evaluateInverse(WorldVector const& y) const
    {
      WorldVector ys = y - shift_;
      WorldVector x;
      matrix_.mtv(ys, x);
      return x;
    }
  };
 
 
 
  template <class Basis, class TP>
  class PeriodicBC
  {
    using SubspaceBasis = Dune::Functions::SubspaceBasis<Basis, TP>;
 
  public:
    using Domain = typename SubspaceBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
    using MultiIndex = typename SubspaceBasis::MultiIndex;
    using FaceTrafo = FaceTransformation<typename Domain::field_type, Domain::dimension>;
    using Intersection = typename SubspaceBasis::GridView::Intersection;
 
  public:
    PeriodicBC(SubspaceBasis basis, BoundarySubset<Intersection> boundarySubset,
               FaceTrafo faceTrafo)
      : basis_(std::move(basis))
      , boundarySubset_(std::move(boundarySubset))
      , faceTrafo_(std::move(faceTrafo))
    {}
 
 
    void init();
 
 
    template <class Mat, class Sol, class Rhs>
    void apply(Mat& matrix, Sol& solution, Rhs& rhs);
 
    auto const& associations() const
    {
      return associations_;
    }
 
    auto const& periodicNodes() const
    {
      return periodicNodes_;
    }
 
  protected:
    void initAssociations();
    void initAssociations2();
 
    // Get the coordinates of the DOFs
    template <class Node>
    auto coords(Node const& node, std::vector<std::size_t> const& localIndices) const;
 
  private:
    SubspaceBasis basis_;
    BoundarySubset<Intersection> boundarySubset_;
    FaceTrafo faceTrafo_;
 
    std::vector<bool> periodicNodes_;
    std::map<std::size_t, std::size_t> associations_;
    std::optional<bool> hasConnectedIntersections_;
  };
 
  namespace Impl {
    template <class Basis>
    using FaceTrafo = FaceTransformation<
      typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate::field_type,
      Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate::dimension>;
  }
 
 
  template <class B, class TP>
  auto makePeriodicBC(Dune::Functions::SubspaceBasis<B, TP> basis,
                      BoundarySubset<typename B::GridView::Intersection> boundarySubset,
                      Impl::FaceTrafo<B> trafo)
  {
    using BcType = PeriodicBC<B, TP>;
    return BcType(std::move(basis), std::move(boundarySubset), std::move(trafo));
  }
 
  template <class B, class TP,
    REQUIRES(Concepts::GlobalBasis<B>)>
  auto makePeriodicBC(B const& basis, TP const& tp,
                      BoundarySubset<typename B::GridView::Intersection> boundarySubset,
                      Impl::FaceTrafo<B> trafo)
  {
    return makePeriodicBC(Dune::Functions::subspaceBasis(basis, tp),
                          std::move(boundarySubset), std::move(trafo));
  }
 
  template <class B,
    REQUIRES(Concepts::GlobalBasis<B>)>
  auto makePeriodicBC(B const& basis,
                      BoundarySubset<typename B::GridView::Intersection> boundarySubset,
                      Impl::FaceTrafo<B> trafo)
  {
    return makePeriodicBC(basis, makeTreePath(), std::move(boundarySubset),
                          std::move(trafo));
  }
 
} // end namespace AMDiS
 
#include "PeriodicBC.inc.hpp"