/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_MESH_FUNCTION_LEVEL_H
#define EXTENSIONS_MESH_FUNCTION_LEVEL_H

using namespace AMDiS;


/** \brief
 * Base-Implementation of a MeshRefinementFunction for refinement of phase-field
 * interfaces. operator() must be overwritten in subclasses. Class provides
 * basic parameters for the refinement like refinement-level values on interface
 * or in the interesting domain and defines the refinement range for the phase-
 * field variable, i.e. interface is defined in [0.05, 0.95]
 **/
template<typename T>
class PhaseFieldRefinementBase : public MeshRefinementFunction<T, int>
{
public:

  PhaseFieldRefinementBase(Mesh *mesh) : 
    MeshRefinementFunction<T, int>(mesh),
    lInner(10),
    lOuter(10),
    lInterface(14),
    minPhase(0.05),
    maxPhase(0.95),
    minOuterPhase(0.001),
    maxOuterPhase(0.999)
  {
      Parameters::get("mesh->refinement->level in inner domain",lInner);
      Parameters::get("mesh->refinement->level in outer domain",lOuter);
      Parameters::get("mesh->refinement->level on interface",lInterface);

      lInterface-= mesh->getMacroElementLevel();
      lInner-= mesh->getMacroElementLevel();
      lOuter-= mesh->getMacroElementLevel();

      int local_globalSize = 10;
      Parameters::get("mesh->refinement->initial level", local_globalSize);
      MeshRefinementFunction<T, int>::globalSize = local_globalSize;

      Parameters::get("mesh->refinement->min interface value",minPhase);
      Parameters::get("mesh->refinement->max interface value",maxPhase);
      Parameters::get("mesh->refinement->min outer interface value",minOuterPhase);
      Parameters::get("mesh->refinement->max outer interface value",maxOuterPhase);
  }

  double meshSize() { return levelToH(MeshRefinementFunction<T, int>::globalSize); }

protected:

  int lInner;
  int lOuter;
  int lInterface;
  
  double minPhase;
  double maxPhase;
  double minOuterPhase;
  double maxOuterPhase;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 **/
class PhaseFieldRefinement : public PhaseFieldRefinementBase<double>
{
public:

    PhaseFieldRefinement(Mesh* mesh_) : PhaseFieldRefinementBase<double>(mesh_) {};
    
    int operator()(const double& phase) const {
      int result= lOuter;
      if (minPhase < phase && phase < maxPhase)
        result= lInterface;         // auf dem Interface
      else if (phase > minOuterPhase && phase <= minPhase)
        result= std::max(lOuter, static_cast<int>(floor((lOuter+lInterface)/2.0)));
      else if (phase < maxOuterPhase && phase >= maxPhase)
        result= std::max(lInner, static_cast<int>(floor((lInner+lInterface)/2.0)));
      else if (phase > (minPhase+maxPhase)/2.0)
        result= lInner;
      return result;
    }
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 * Additionally to the phase-field value a list of refinement points can be given
 * where the mesh should be refined up to the a given level, e.g. in corners of 
 * the macro-mesh.
 **/
class PhaseFieldCoordsRefinement : public PhaseFieldRefinementBase< std::pair<WorldVector<double>, double> >
{
public:
    PhaseFieldCoordsRefinement(Mesh* mesh_, std::vector<WorldVector<double> > points_, double radius_) :
      PhaseFieldRefinementBase< std::pair<WorldVector<double>, double> >(mesh_),
      lPoints(14),
      points(points_),
      radius(radius_)
    {
        Parameters::get("mesh->refinement->level on points",lPoints);
        lPoints-= mesh->getMacroElementLevel();
    }
    
    int operator()(const std::pair<WorldVector<double>, double>& data) const {
      double phase = data.second;
      int result= lOuter;
      if (minPhase < phase && phase < maxPhase)
        result = lInterface;         // auf dem Interface
      else if ((phase > minOuterPhase && phase <= minPhase) || (phase < maxOuterPhase && phase >= maxPhase))
        result = std::max(lInner, static_cast<int>(floor((lOuter+lInterface)/2.0)));
      else if (phase > 0.5)
        result = lInner;

      WorldVector<double> x = data.first;
      double minDist = 1.e15;
      for (unsigned i = 0; i < points.size(); ++i) {
        minDist = std::min(minDist, norm(points[i]-x));
      }
      double lambda = std::max(0.0, 1.0 - minDist/radius);
      result = std::max(result, static_cast<int>(floor(lOuter+lambda*(lPoints-lOuter))));

      return result;
    }
    
private:

    int lPoints;

    std::vector<WorldVector<double> > points;

    double radius;
};


/** \brief
 * Implementation of MeshRefinementFunction. Provides the method operator()
 * that defines a list of refinement points where the mesh should be refined 
 * up to the a given level, e.g. in corners of the macro-mesh.
 **/
class CoordsRefinement : public MeshRefinementFunction< WorldVector<double>, int >
{
public:
    CoordsRefinement(Mesh* mesh_, std::vector<WorldVector<double> > points_, double radius_) 
    :
    MeshRefinementFunction< WorldVector<double>, int >(mesh_),
    points(points_),
    radius(radius_)
    {
        lInner= 15; lPoints = 20;
        Parameters::get("mesh->refinement->level on points",lPoints);
        Parameters::get("mesh->refinement->level in inner domain",lInner);
        lPoints-= mesh->getMacroElementLevel();
        lInner-= mesh->getMacroElementLevel();

        globalSize= lInner;
    }
    
    int operator()(const WorldVector<double>& x) const {
      double minDist = 1.e15;
      for (unsigned i = 0; i < points.size(); ++i) {
        minDist = std::min(minDist, norm(points[i]-x));
      }
      int result = lInner;
      double lambda = std::max(0.0, 1.0 - minDist/radius);
      result = std::max(result, static_cast<int>(floor(lInner+lambda*(lPoints-lInner))));

      return result;
    }
    double meshSize() { return levelToH(lInner); }
    
private:
    int lInner, lPoints;
    std::vector<WorldVector<double> > points;
    double radius;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the level of refinement depending on a list of phase-fields, e.g.
 * a diffuse-domain representation and a Cahn-Hilliard field.
 **/
class PhaseFieldChRefinement : public PhaseFieldRefinementBase< std::vector<double> >
{
public:
  PhaseFieldChRefinement(Mesh* mesh_) 
  :
  PhaseFieldRefinementBase< std::vector<double> >(mesh_),
  lInterfaceDomain(14)
  {
    Parameters::get("mesh->refinement->level on domain interface",lInterfaceDomain);
    lInterfaceDomain-= mesh->getMacroElementLevel();
    
    Parameters::get("mesh->refinement->initial level", globalSize);
  }
  
  int operator()(const std::vector<double>& phases) const 
  {
    int result= lOuter;
    for (unsigned i = 0; i < phases.size(); ++i) {
      int localResult= lOuter;
      if (minPhase < phases[i] && phases[i] < maxPhase
        && (i < 1 || phases[i-1] > 0.5))
        localResult= (i==0 ? lInterfaceDomain : lInterface);         // auf dem Interface
      else if (phases[i] > 0.5
        && (i < 1 || phases[i-1] > 0.5))
        localResult= lInner;            // im Innern des Gebietes
      if (((phases[i] > minOuterPhase && phases[i] <= minPhase) || (phases[i] < maxOuterPhase && phases[i] >= maxPhase))
        && (i < 1 || phases[i-1] > 0.5))
        localResult= static_cast<int>(floor((lInner + lInterface) / 2.0));
      result= std::max(result, localResult);
    }
    return result;
  }
    
private:

  int lInterfaceDomain;
};


/** \brief
 * Base-Implementation of a MeshRefinementFunction for refinement of phase-field
 * interfaces. operator() must be overwritten in subclasses. Class provides
 * basic parameters for the refinement like refinement-level values on interface
 * or in the interesting domain and defines the refinement range for the phase-
 * field variable, i.e. interface is defined in [0.05, 0.95]
 **/
template<typename T>
class SignedDistRefinementBase : public MeshRefinementFunction<T, int>
{
public:

  SignedDistRefinementBase(Mesh *mesh) :
    MeshRefinementFunction<T, int>(mesh),
    lInner(10),
    lOuter(10),
    lInterface(14),
    interfaceWidth(0.2),
    fadeOutWidth(0.0),
    signInInnerDomain(-1.0)
  {
      Parameters::get("mesh->refinement->level in inner domain",lInner);
      Parameters::get("mesh->refinement->level in outer domain",lOuter);
      Parameters::get("mesh->refinement->level on interface",lInterface);

      lInterface-= mesh->getMacroElementLevel();
      lInner-= mesh->getMacroElementLevel();
      lOuter-= mesh->getMacroElementLevel();

      int local_globalSize = 10;
      Parameters::get("mesh->refinement->initial level", local_globalSize);
      MeshRefinementFunction<T, int>::globalSize = local_globalSize;

      Parameters::get("mesh->refinement->interface width",interfaceWidth);
      Parameters::get("mesh->refinement->fade out width",fadeOutWidth);
      Parameters::get("mesh->refinement->sign in inner domain",signInInnerDomain);
  }

  double meshSize() { return levelToH(MeshRefinementFunction<T, int>::globalSize); }

protected:

  int lInner;
  int lOuter;
  int lInterface;

  double interfaceWidth;
  double fadeOutWidth;
  double signInInnerDomain;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 **/
class SignedDistRefinement : public SignedDistRefinementBase<double>
{
public:

    SignedDistRefinement(Mesh* mesh_) : SignedDistRefinementBase<double>(mesh_) {};

    int operator()(const double& dist) const {
      int result= lOuter;
      if (abs(dist) < interfaceWidth)
        result= lInterface;         // auf dem Interface
      else if (abs(dist) < fadeOutWidth+interfaceWidth && signInInnerDomain*dist < 0.0) {
	double lambda = abs(dist)/fadeOutWidth - interfaceWidth/fadeOutWidth;
        result= static_cast<int>(floor(lambda*lInner+(1.0-lambda)*lInterface));
      } else if (abs(dist) < fadeOutWidth+interfaceWidth && signInInnerDomain*dist > 0.0) {
	double lambda = abs(dist)/fadeOutWidth - interfaceWidth/fadeOutWidth;
        result= static_cast<int>(floor(lambda*lOuter+(1.0-lambda)*lInterface));
      } else if (signInInnerDomain*dist < 0.0)
        result= lInner;
      return result;
    }
};


class ESIndicator : public MeshRefinementFunction< std::vector<double>, int >
{
public:
  ESIndicator(Mesh* mesh_) 
  :
  MeshRefinementFunction< std::vector<double>, int >(mesh_)
  {
    Parameters::get("mesh->refinement->initial level", globalSize);
  }
  
  int operator()(const std::vector<double>& means) const 
  {
    double tol = means[0];
    return (abs(means[1]-means[2]) > tol ? 0 : -1);
  }
    
  double indicator(const std::vector<double>& means) const 
  {
    return abs(means[1]-means[2]);
  }
};

#endif // EXTENSIONS_MESH_FUNCTION_LEVEL_H