/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2019 Adrian Kummerländer
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

#ifndef REFINEMENT_GRID_2D_HH
#define REFINEMENT_GRID_2D_HH

#include "grid2D.h"
#include "coupler2D.hh"

#include "utilities/vectorHelpers.h"
#include "communication/heuristicLoadBalancer.h"

namespace olb {


template <typename T, template<typename> class DESCRIPTOR>
Grid2D<T,DESCRIPTOR>::Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, int resolution, T tau, int re):
  _domainF(std::move(domainF)),
  _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
               resolution, // resolution: number of voxels per charPhysL
               tau,        // latticeRelaxationTime: relaxation time, has to be greater than 0.5!
               T{1},       // charPhysLength: reference length of simulation geometry
               T{1},       // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
               T{1./re},   // physViscosity: physical kinematic viscosity in __m^2 / s__
               T{1},       // physDensity: physical density in __kg / m^3__
               T{1})),
  _cuboids(new CuboidGeometry2D<T>(
             *_domainF,
             _converter->getConversionFactorLength(),
#ifdef PARALLEL_MODE_MPI
             singleton::mpi().getSize()
#else
             1
#endif
           )),
  _balancer(new HeuristicLoadBalancer<T>(
              *_cuboids)),
  _geometry(new SuperGeometry2D<T>(
              *_cuboids,
              *_balancer,
              2)),
  _lattice(new SuperLattice2D<T,DESCRIPTOR>(
             *_geometry))
{
  _converter->print();
}

template <typename T, template<typename> class DESCRIPTOR>
UnitConverter<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getConverter()
{
  return *_converter;
}

template <typename T, template<typename> class DESCRIPTOR>
CuboidGeometry2D<T>& Grid2D<T,DESCRIPTOR>::getCuboidGeometry()
{
  return *_cuboids;
}

template <typename T, template<typename> class DESCRIPTOR>
LoadBalancer<T>& Grid2D<T,DESCRIPTOR>::getLoadBalancer()
{
  return *_balancer;
}

template <typename T, template<typename> class DESCRIPTOR>
SuperGeometry2D<T>& Grid2D<T,DESCRIPTOR>::getSuperGeometry()
{
  return *_geometry;
}

template <typename T, template<typename> class DESCRIPTOR>
SuperLattice2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getSuperLattice()
{
  return *_lattice;
}

template <typename T, template<typename> class DESCRIPTOR>
void Grid2D<T,DESCRIPTOR>::collideAndStream()
{
  for ( auto& fineCoupler : _fineCouplers ) {
    fineCoupler->store();
  }

  this->getSuperLattice().collideAndStream();

  for ( auto& fineGrid : _fineGrids ) {
    fineGrid->collideAndStream();
  }

  for ( auto& fineCoupler : _fineCouplers ) {
    fineCoupler->interpolate();
    fineCoupler->couple();
  }

  for ( auto& fineGrid : _fineGrids ) {
    fineGrid->collideAndStream();
  }

  for ( auto& fineCoupler : _fineCouplers ) {
    fineCoupler->store();
    fineCoupler->couple();
  }

  for ( auto& coarseCoupler : _coarseCouplers ) {
    coarseCoupler->couple();
  }
}

template <typename T, template<typename> class DESCRIPTOR>
FineCoupler2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::addFineCoupling(
  Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend)
{
  _fineCouplers.emplace_back(
    new FineCoupler2D<T,DESCRIPTOR>(
      *this, fineGrid, origin, extend));
  return *_fineCouplers.back();
}

template <typename T, template<typename> class DESCRIPTOR>
CoarseCoupler2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::addCoarseCoupling(
  Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend)
{
  _coarseCouplers.emplace_back(
    new CoarseCoupler2D<T,DESCRIPTOR>(
      *this, fineGrid, origin, extend));
  return *_coarseCouplers.back();
}

template <typename T, template<typename> class DESCRIPTOR>
Vector<T,2> Grid2D<T,DESCRIPTOR>::alignOriginToGrid(Vector<T,2> physR) const
{
  Vector<int,3> latticeR{};
  _cuboids->getLatticeR(physR, latticeR);
  return _cuboids->getPhysR(latticeR.toStdVector());
}

template <typename T, template<typename> class DESCRIPTOR>
Vector<T,2> Grid2D<T,DESCRIPTOR>::alignExtendToGrid(Vector<T,2> extend) const
{
  const T deltaX = _converter->getPhysDeltaX();
  return util::floor(extend / deltaX) * deltaX;
}

template <typename T, template<typename> class DESCRIPTOR>
RefiningGrid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(
  Vector<T,2> wantedOrigin, Vector<T,2> wantedExtend, bool addCouplers)
{
  if (addCouplers) {
    auto& fineGrid = refine(wantedOrigin, wantedExtend, false);

    const Vector<T,2> origin = fineGrid.getOrigin();
    const Vector<T,2> extend = fineGrid.getExtend();

    const Vector<T,2> extendX = {extend[0],0};
    const Vector<T,2> extendY = {0,extend[1]};

    addFineCoupling(fineGrid, origin,           extendY);
    addFineCoupling(fineGrid, origin + extendX, extendY);
    addFineCoupling(fineGrid, origin + extendY, extendX);
    addFineCoupling(fineGrid, origin,           extendX);

    const T coarseDeltaX = getConverter().getPhysDeltaX();
    const Vector<T,2> innerOrigin  = origin + coarseDeltaX;
    const Vector<T,2> innerExtendX = extendX - Vector<T,2> {2*coarseDeltaX,0};
    const Vector<T,2> innerExtendY = extendY - Vector<T,2> {0,2*coarseDeltaX};

    addCoarseCoupling(fineGrid, innerOrigin,                innerExtendY);
    addCoarseCoupling(fineGrid, innerOrigin + innerExtendX, innerExtendY);
    addCoarseCoupling(fineGrid, innerOrigin + innerExtendY, innerExtendX);
    addCoarseCoupling(fineGrid, innerOrigin,                innerExtendX);

    return fineGrid;
  }
  else {
    const Vector<T,2> origin = alignOriginToGrid(wantedOrigin);
    const Vector<T,2> extend = alignExtendToGrid(wantedExtend);

    _fineGrids.emplace_back(
      new RefiningGrid2D<T,DESCRIPTOR>(*this, origin, extend));
    return *_fineGrids.back();
  }
}


template <typename T, template<typename> class DESCRIPTOR>
RefiningGrid2D<T,DESCRIPTOR>::RefiningGrid2D(
  Grid2D<T,DESCRIPTOR>& parentGrid, Vector<T,2> origin, Vector<T,2> extend):
  Grid2D<T,DESCRIPTOR>(
    std::unique_ptr<IndicatorF2D<T>>(new IndicatorCuboid2D<T>(extend, origin)),
    2*parentGrid.getConverter().getResolution(),
    2*parentGrid.getConverter().getLatticeRelaxationTime() - 0.5,
    parentGrid.getConverter().getReynoldsNumber()),
  _origin(origin),
  _extend(extend),
  _parentGrid(parentGrid)
{ }

template <typename T, template<typename> class DESCRIPTOR>
Vector<T,2> RefiningGrid2D<T,DESCRIPTOR>::getOrigin() const
{
  return _origin;
}

template <typename T, template<typename> class DESCRIPTOR>
Vector<T,2> RefiningGrid2D<T,DESCRIPTOR>::getExtend() const
{
  return _extend;
}

template <typename T, template<typename> class DESCRIPTOR>
std::unique_ptr<IndicatorF2D<T>> RefiningGrid2D<T,DESCRIPTOR>::getRefinedOverlap() const
{
  const T coarseDeltaX = _parentGrid.getConverter().getPhysDeltaX();

  return std::unique_ptr<IndicatorF2D<T>>(
           new IndicatorCuboid2D<T>(_extend - 4*coarseDeltaX, _origin + 2*coarseDeltaX));
}

}

#endif