PeriodicBC.inc.hpp

#pragma once

#include <algorithm>
#include <array>
#include <cmath>
#include <limits>

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

#include <amdis/LinearAlgebra.hpp>
#include <amdis/Output.hpp>
#include <amdis/functions/FunctionFromCallable.hpp>
#include <amdis/typetree/Traversal.hpp>

namespace AMDiS {

template <class Mat, class Sol, class Rhs, class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>::
init()
{
  if (!bool(hasConnectedIntersections_)) {
    hasConnectedIntersections_ = true;
    const auto& gridView = basis_.gridView().grid().levelGridView(0);
    for (auto const& element : elements(gridView)) {
      for (const auto& intersection : intersections(gridView, element)) {
        if (!boundarySubset_(intersection))
          continue;

        if (!intersection.neighbor()) {
          hasConnectedIntersections_ = false;
          break;
        }
      }
    }
  }

  if (*hasConnectedIntersections_)
    initAssociations();
  else
    initAssociations2();
}


template <class Mat, class Sol, class Rhs, class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>::
initAssociations()
{
  using std::sqrt;
  static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
  periodicNodes_.assign(basis_.dimension(), false);

  auto localView = basis_.localView();
  auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
  const auto& gridView = basis_.gridView();
  for (auto const& element : elements(gridView)) {
    if (!element.hasBoundaryIntersections())
      continue;

    localView.bind(element);
    for (const auto& intersection : intersections(gridView, element)) {
      if (!boundarySubset_(intersection))
        continue;

      test_exit(intersection.neighbor(),
        "Neighbors of periodic intersection not assigned. Use a periodic Grid instead!");

      // collect DOFs of inside element on intersection
      localView.bind(intersection.inside());
      seDOFs.bind(localView, intersection.indexInInside(), 1);
      std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
      std::vector<MultiIndex> insideGlobalDOFs(insideDOFs.size());
      std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
        [&](std::size_t localIndex) { return localView.index(localIndex); });
      auto insideCoords = coords(localView.tree(), insideDOFs);

      // collect DOFs of ouside element on intersection
      localView.bind(intersection.outside());
      seDOFs.bind(localView, intersection.indexInOutside(), 1);
      std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
      auto outsideCoords = coords(localView.tree(), outsideDOFs);

      // compare mapped coords of DOFs
      assert(insideDOFs.size() == outsideDOFs.size());
      for (std::size_t i = 0; i < insideCoords.size(); ++i) {
        auto x = faceTrafo_.evaluate(insideCoords[i]);
        for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
          auto const& y = outsideCoords[j];
          if (distance(x,y) < tol) {
            periodicNodes_[insideGlobalDOFs[i]] = true;
            associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
          }
        }
      }
    }
  }
}

namespace Impl
{
  template <class D>
  class CompareTol
  {
    using ctype = typename D::field_type;

  public:
    CompareTol(ctype tol)
      : tol_(tol)
    {}

    bool operator()(D const& lhs, D const& rhs) const
    {
      using std::abs;
      for (int i = 0; i < D::dimension; i++) {
        if (abs(lhs[i] - rhs[i]) < tol_)
          continue;
        return lhs[i] < rhs[i];
      }
      return false;
    }

  private:
    ctype tol_;
  };
}


template <class Mat, class Sol, class Rhs, class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>::
initAssociations2()
{
  using std::sqrt;
  static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
  periodicNodes_.assign(basis_.dimension(), false);

  // Dune::ReservedVector<std::pair<EntitySeed,int>,2>
  using EntitySeed = typename SubspaceBasis::GridView::Grid::template Codim<0>::EntitySeed;
  using CoordAssoc = std::map<Domain, std::vector<std::pair<EntitySeed,int>>, Impl::CompareTol<Domain>>;
  CoordAssoc coordAssoc(Impl::CompareTol<Domain>{tol});

  // generate periodic element pairs
  const auto& gridView = basis_.gridView();
  for (auto const& element : elements(gridView)) {
    for (const auto& intersection : intersections(gridView, element)) {
      if (!boundarySubset_(intersection))
        continue;

      auto x = intersection.geometry().center();
      auto seed = intersection.inside().seed();
      coordAssoc[x].push_back({seed,intersection.indexInInside()});
      coordAssoc[faceTrafo_.evaluate(x)].push_back({seed,intersection.indexInInside()});
    }
  }

  auto localView = basis_.localView();
  auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
  for (auto const& assoc : coordAssoc) {
    auto const& seeds = assoc.second;
    if (seeds.size() != 2)
      continue;

    // collect DOFs of inside element on intersection
    auto el1 = gridView.grid().entity(seeds[0].first);
    localView.bind(el1);
    seDOFs.bind(localView, seeds[0].second, 1);
    std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
    std::vector<MultiIndex> insideGlobalDOFs(insideDOFs.size());
    std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
      [&](std::size_t localIndex) { return localView.index(localIndex); });
    auto insideCoords = coords(localView.tree(), insideDOFs);

    // collect DOFs of outside element on intersection
    auto el2 = gridView.grid().entity(seeds[1].first);
    localView.bind(el2);
    seDOFs.bind(localView, seeds[1].second, 1);
    std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
    auto outsideCoords = coords(localView.tree(), outsideDOFs);

    // compare mapped coords of DOFs
    assert(insideDOFs.size() == outsideDOFs.size());
    for (std::size_t i = 0; i < insideCoords.size(); ++i) {
      auto x = faceTrafo_.evaluate(insideCoords[i]);
      for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
        auto const& y = outsideCoords[j];
        if (distance(x,y) < tol) {
          periodicNodes_[insideGlobalDOFs[i]] = true;
          associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
        }
      }
    }
  }
}


template <class Mat, class Sol, class Rhs, class Basis, class TP>
  template <class Node>
auto PeriodicBC<Mat, Sol, Rhs, Basis, TP>::
coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
{
  std::vector<Domain> dofCoords(localIndices.size());
  for_each_leaf_node(tree, [&](auto const& node, auto&&)
  {
    std::size_t size = node.finiteElement().size();
    auto geometry = node.element().geometry();

    auto const& localInterpol = node.finiteElement().localInterpolation();
    using FiniteElement = TYPEOF(node.finiteElement());
    using DomainType = typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
    using RangeType = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;

    std::array<std::vector<typename RangeType::field_type>, Domain::dimension> coeffs;
    for (int d = 0; d < Domain::dimension; ++d) {
      auto evalCoord = [&](DomainType const& local) -> RangeType { return geometry.global(local)[d]; };
      auto evalCoordFct = functionFromCallable<RangeType(DomainType)>(evalCoord);
      localInterpol.interpolate(evalCoordFct, coeffs[d]);
    }

    for (std::size_t j = 0; j < localIndices.size(); ++j) {
      Domain x;
      for (std::size_t i = 0; i < size; ++i) {
        if (node.localIndex(i) == localIndices[j]) {
          for (int d = 0; d < Domain::dimension; ++d)
            x[d] = coeffs[d][i];
          break;
        }
      }
      dofCoords[j] = std::move(x);
    }
  });

  return dofCoords;
}


template <class Mat, class Sol, class Rhs, class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>::
apply(Mat& matrix, Sol& solution, Rhs& rhs)
{
  periodicBC(matrix, solution, rhs, periodicNodes_, associations_);
}

} // end namespace AMDiS