/******************************************************************************
 *
 * 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_REFINEMENT_H
#define EXTENSIONS_REFINEMENT_H

#include "ElementFunction.h"
#include "RefineOperations.h"
#include "MeshFunction_Level.h"

namespace AMDiS { namespace extensions {

/** \brief
 * Base class for Refinement structure to perform local anisotropic refinement
 */
template<typename T, typename T2>
class RefinementLevel
{
public:

  RefinementLevel(const FiniteElemSpace *feSpace_, 
		  MeshRefinementFunction<T,T2>* refineFct_ = NULL) :
    feSpace(feSpace_), 
    refineFct(refineFct_),
    adaptInfo(NULL),
    refineOperation(NULL),
    numRefinements0(15),
    finished(true),
    globalRefined(false)
  {	
    FUNCNAME("RefinementLevel::RefinementLevel()");
    
    mesh = feSpace->getMesh();
    switch (mesh->getDim()) {
    case 1:
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
      break;
    case 2:
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
      break;
    case 3:
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }

    numRefinements = numRefinements0;
    
    refineOperation = new StandardRefineOperation;
  }
  
  virtual ~RefinementLevel() 
  {
//     finalize();
    
    delete coarseningManager;
    delete refinementManager;
    if (refineOperation)
      delete refineOperation;
  }
  
  void finalize()
  {
    if (!finished) {
      #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      Parallel::MeshDistributor::globalMeshDistributor->checkMeshChange();
      #endif
    }
    finished = true;
  }
  
  void setRefineOperation(AdaptInfo* adaptInfo_,
			  StandardRefineOperation* refineOperation_)
  {
    if (refineOperation)
      delete refineOperation;
    
    adaptInfo = adaptInfo_;
    refineOperation = refineOperation_;
  }

  void refine(bool onlyRefine= false, bool meshRepartation = false) 
  {
    FUNCNAME("RefinementLevel::refine()");

    if (!globalRefined && refineFct) {
      MSG("nr of global refinements: %d\n", refineFct->getGlobalSize());
      refinementManager->globalRefine(mesh, refineFct->getGlobalSize());
      globalRefined = true;
    }
    double minH = 0.0, maxH = 1.0;
    int minLevel = 100, maxLevel = 0;
    
    // build mesh for phasefield-function
    bool meshChanged = true;
    Flag markFlag;
    int oldNr = 0, oldOldNr = 0;
    int i = 0;
    while (meshChanged && i < numRefinements) {
      markElements(markFlag);
      meshChanged = refineMesh(markFlag, onlyRefine);
      
      calcMeshSizes(minH, maxH, minLevel, maxLevel, false); 
      int nr = mesh->getNumberOfVertices();
      meshChanged = meshChanged && oldOldNr!=nr && oldNr!=nr;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      int sendValue = static_cast<int>(meshChanged);
      Parallel::mpi::globalAdd(MPI::COMM_WORLD, sendValue);
      meshChanged = static_cast<bool>(sendValue);
#endif
      if (meshChanged) {
	MSG("(local) mesh sizes: [%f, %f], Vs: %d, ELs: %d\n", 
          minH, maxH, nr, mesh->getNumberOfElements());
      }
      i++;
			
      oldOldNr = oldNr;
      oldNr = nr; 
    }
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::MeshDistributor::globalMeshDistributor->checkMeshChange(meshRepartation);
#endif
    calcMeshSizes(minH, maxH, minLevel, maxLevel, true); 
    MSG("Final (global) mesh: [%f, %f], Vs: %d, ELs: %d, Level: [%d, %d]\n",
      minH, maxH, mesh->getNumberOfVertices(), mesh->getNumberOfElements(), minLevel, maxLevel);
  }

  void refine(int numRefinements_, bool onlyRefine= false, bool repartation = false) { 
    numRefinements = numRefinements_; 
    refine(onlyRefine, repartation);  
    numRefinements = numRefinements0; 
  }
  
  
  void initial_refine()
  {
    
  }

  int getNumRefinements(){
      return numRefinements;
  }

  void calcMeshSizes(double& minH, double& maxH, int& minLevel, int& maxLevel, bool allReduce = false) 
  {
    FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    minH = 1e15; maxH = 0.0;
    int k = 0;
    minLevel = 100;
    maxLevel = 0;
    while (elInfo) {
      maxLevel = std::max(maxLevel,elInfo->getLevel());
      minLevel = std::min(minLevel,elInfo->getLevel());
      coords = elInfo->getCoords();
      double h = 0.0;
      for (int i = 0; i < coords.getSize(); i++) {
        for (int j = 0; j < coords.getSize(); j++) {
	  if (i != j)
            h = std::max(h,norm(coords[i]-coords[j]));
        }
      }
      minH = std::min(h, minH);
      maxH = std::max(h, maxH);
      elInfo = stack.traverseNext(elInfo);
      k++;
    }
    minLevel += mesh->getMacroElementLevel();
    maxLevel += mesh->getMacroElementLevel();
    
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (allReduce) {
      Parallel::mpi::globalMin(minH);
      Parallel::mpi::globalMax(maxH);
      Parallel::mpi::globalMin(minLevel);
      Parallel::mpi::globalMax(maxLevel);
    }
#endif
  }


  double calcMeshSize(ElInfo *elInfo, bool allReduce = false) 
  {
    FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);
    coords = elInfo->getCoords();
    double h = 0.0;
    for (int i = 0; i < coords.getSize(); i++) {
      for (int j = 0; j < coords.getSize(); j++) {
        if (i != j)
          h = std::max(h,norm(coords[i]-coords[j]));
      }
    }
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (allReduce)
      Parallel::mpi::globalMax(h);
#endif

    return h;
  }


  int calcMark(double refineH, double currentH)
  {
    return (refineH < currentH ? 
        1 : (refineH > currentH * (mesh->getDim() == 1 ? 
          2.0 : (mesh->getDim() == 2 ? 
          std::sqrt(2.0) : 
          std::sqrt(2.0)/2.0 + 0.5)) ? 
        -1 : 
        0));
  }


  virtual int calcMark(int refineLevel, int currentLevel)
  {
    int levelDiff = refineLevel - currentLevel;
    return (levelDiff > 0 ? 1 : (levelDiff < 0 ? -1 : 0));
  }

	
  bool refineMesh(Flag markFlag, bool onlyRefine) 
  {
    finished = false;
    
    int oldSize = mesh->getNumberOfVertices();
    refineOperation->beforeRefine(adaptInfo, markFlag);
    if (markFlag.isSet(1))
      refinementManager->refineMesh(mesh);
    refineOperation->beforeCoarsen(adaptInfo, markFlag);
    if (markFlag.isSet(2) && !onlyRefine)
      coarseningManager->coarsenMesh(mesh);
    refineOperation->afterCoarsen(adaptInfo, markFlag);
    if (markFlag.isSet(1) || markFlag.isSet(2)) {
      int newSize = mesh->getNumberOfVertices();
      if (oldSize != newSize) 
        return true;
    }
    return false;
  }

  virtual void markElements(Flag &markFlag) = 0;

  void setGlobalRefined(bool refined) { globalRefined = refined; }
  RefinementManager* getRefinementManager() { return refinementManager; }
  CoarseningManager* getCoarseningManager() { return coarseningManager; }

protected:

  const FiniteElemSpace *feSpace;
  Mesh* mesh;
  MeshRefinementFunction<T,T2>* refineFct;
  RefinementManager* refinementManager;
  CoarseningManager* coarseningManager;
  
  AdaptInfo* adaptInfo;
  StandardRefineOperation* refineOperation;

  int numRefinements;
  int numRefinements0;
  
  bool finished;
  bool globalRefined;
};



/** \brief
 * Base class for Refinement structure to perform local anisotropic refinement
 */
class RefinementLocal
{
public:

  RefinementLocal(Mesh *mesh_) :
    mesh(mesh_), 
    adaptInfo(NULL),
    refineOperation(NULL),
    numRefinements0(15),
    finished(true),
    globalRefined(false)
  {	
    FUNCNAME("RefinementLevel::RefinementLevel()");
    
    switch (mesh->getDim()) {
    case 1:
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
      break;
    case 2:
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
      break;
    case 3:
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }

    numRefinements = numRefinements0;
    
    refineOperation = new StandardRefineOperation;
  }
  
  virtual ~RefinementLocal() 
  {
//     finalize();
    
    delete coarseningManager;
    delete refinementManager;
    if (refineOperation)
      delete refineOperation;
  }
  
  void finalize()
  {
    if (!finished) {
      #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      Parallel::MeshDistributor::globalMeshDistributor->checkMeshChange();
      #endif
    }
    finished = true;
  }
  
  void setRefineOperation(AdaptInfo* adaptInfo_,
			  StandardRefineOperation* refineOperation_)
  {
    if (refineOperation)
      delete refineOperation;
    
    adaptInfo = adaptInfo_;
    refineOperation = refineOperation_;
  }

  void setOnlyRefine(bool onlyRefine_)
  {
    onlyRefine = onlyRefine_;
  }
  
  template<typename MeshRefFunc>
  void refine(MeshRefFunc* refineFct) 
  {
    FUNCNAME("RefinementLevel::refine()");

    refineFct->init(mesh);
    if (!globalRefined) {
      MSG("nr of global refinements: %d\n", refineFct->getGlobalSize());
      refinementManager->globalRefine(mesh, refineFct->getGlobalSize());
      globalRefined = true;
    }
    double minH = 0.0, maxH = 1.0;
    int minLevel = 100, maxLevel = 0;
    
    // build mesh for phasefield-function
    bool meshChanged = true;
    Flag markFlag;
    int oldNr = 0, oldOldNr = 0;
    int i = 0;
    while (meshChanged && i < numRefinements) {
      markElements(refineFct, markFlag);
      meshChanged = refineMesh(markFlag, onlyRefine);
      
      calcMeshSizes(minH, maxH, minLevel, maxLevel, false); 
      int nr = mesh->getNumberOfVertices();
      meshChanged = meshChanged && oldOldNr!=nr && oldNr!=nr;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      int sendValue = static_cast<int>(meshChanged);
      Parallel::mpi::globalAdd(MPI::COMM_WORLD, sendValue);
      meshChanged = static_cast<bool>(sendValue);
#endif
      if (meshChanged) {
	MSG("(local) mesh sizes: [%f, %f], Vs: %d, ELs: %d\n", 
          minH, maxH, nr, mesh->getNumberOfElements());
      }
      i++;
			
      oldOldNr = oldNr;
      oldNr = nr; 
    }
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::MeshDistributor::globalMeshDistributor->checkMeshChange();
#endif
    calcMeshSizes(minH, maxH, minLevel, maxLevel, true); 
    MSG("Final (global) mesh: [%f, %f], Vs: %d, ELs: %d, Level: [%d, %d]\n",
      minH, maxH, mesh->getNumberOfVertices(), mesh->getNumberOfElements(), minLevel, maxLevel);
  }

  template<typename MeshRefFunc>
  void refine(int numRefinements_, MeshRefFunc* refineFct, bool onlyRefine_= false) { 
    numRefinements = numRefinements_; 
    setOnlyRefine(onlyRefine_);
    refine(refineFct);  
    numRefinements = numRefinements0; 
  }
  
  template<typename MeshRefFunc>
  void markElements(MeshRefFunc* refineFct, Flag &markFlag) 
  {    
    bool elMarkRefine = false, elMarkCoarsen = false;

    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
    while(elInfo) {
      refineFct->init(elInfo);
      int refineLevel = refineFct->getLevel();
      int oldLevel = elInfo->getLevel();
      elInfo->getElement()->setMark( calcMark(refineLevel, oldLevel) );

      elMarkRefine |= elInfo->getElement()->getMark() == 1;
      elMarkCoarsen |= elInfo->getElement()->getMark() == -1;
      elInfo = stack.traverseNext(elInfo);
    }
      
    markFlag= 0;
    if(elMarkRefine) 
	    markFlag= 1;
    if(elMarkCoarsen) 
	    markFlag|= 2;
  }

  int getNumRefinements() {
      return numRefinements;
  }

  void calcMeshSizes(double& minH, double& maxH, int& minLevel, int& maxLevel, bool allReduce = false) 
  {
    FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    minH = 1e15; maxH = 0.0;
    int k = 0;
    minLevel = 100;
    maxLevel = 0;
    while (elInfo) {
      maxLevel = std::max(maxLevel,elInfo->getLevel());
      minLevel = std::min(minLevel,elInfo->getLevel());
      coords = elInfo->getCoords();
      double h = 0.0;
      for (int i = 0; i < coords.getSize(); i++) {
        for (int j = 0; j < coords.getSize(); j++) {
	  if (i != j)
            h = std::max(h,norm(coords[i]-coords[j]));
        }
      }
      minH = std::min(h, minH);
      maxH = std::max(h, maxH);
      elInfo = stack.traverseNext(elInfo);
      k++;
    }
    minLevel += mesh->getMacroElementLevel();
    maxLevel += mesh->getMacroElementLevel();
    
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (allReduce) {
      Parallel::mpi::globalMin(minH);
      Parallel::mpi::globalMax(maxH);
      Parallel::mpi::globalMin(minLevel);
      Parallel::mpi::globalMax(maxLevel);
    }
#endif
  }


  double calcMeshSize(ElInfo *elInfo, bool allReduce = false) 
  {
    FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);
    coords = elInfo->getCoords();
    double h = 0.0;
    for (int i = 0; i < coords.getSize(); i++) {
      for (int j = 0; j < coords.getSize(); j++) {
        if (i != j)
          h = std::max(h,norm(coords[i]-coords[j]));
      }
    }
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (allReduce)
      Parallel::mpi::globalMax(h);
#endif

    return h;
  }


  int calcMark(double refineH, double currentH)
  {
    return (refineH < currentH ? 
        1 : (refineH > currentH * (mesh->getDim() == 1 ? 
          2.0 : (mesh->getDim() == 2 ? 
          std::sqrt(2.0) : 
          std::sqrt(2.0)/2.0 + 0.5)) ? 
        -1 : 
        0));
  }


  virtual int calcMark(int refineLevel, int currentLevel)
  {
    int levelDiff = refineLevel - currentLevel;
    return (levelDiff > 0 ? 1 : (levelDiff < 0 ? -1 : 0));
  }

	
  bool refineMesh(Flag markFlag, bool onlyRefine) 
  {
    finished = false;
    
    int oldSize = mesh->getNumberOfVertices();
    refineOperation->beforeRefine(adaptInfo, markFlag);
    if (markFlag.isSet(1))
      refinementManager->refineMesh(mesh);
    refineOperation->beforeCoarsen(adaptInfo, markFlag);
    if (markFlag.isSet(2) && !onlyRefine)
      coarseningManager->coarsenMesh(mesh);
    refineOperation->afterCoarsen(adaptInfo, markFlag);
    if (markFlag.isSet(1) || markFlag.isSet(2)) {
      int newSize = mesh->getNumberOfVertices();
      if (oldSize != newSize) 
        return true;
    }
    return false;
  }

  void setGlobalRefined(bool refined) { globalRefined = refined; }
  RefinementManager* getRefinementManager() { return refinementManager; }
  CoarseningManager* getCoarseningManager() { return coarseningManager; }

protected:

  Mesh* mesh;
  RefinementManager* refinementManager;
  CoarseningManager* coarseningManager;
  
  AdaptInfo* adaptInfo;
  StandardRefineOperation* refineOperation;

  int numRefinements;
  int numRefinements0;
  
  bool finished;
  bool globalRefined;
  bool onlyRefine;
};

} }

using namespace AMDiS::extensions;

#include "Refinement_Level.h"
// #include "Refinement_MeshSize.h"

#endif // EXTENSIONS_REFINEMENT_H