#include "RefinementManager.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "DOFAdmin.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FixVec.h"
#include "RCNeighbourList.h"
#include "ProblemStatBase.h"
#include "DOFIndexed.h"
#include "Projection.h"
#include "LeafData.h"
#include "VertexVector.h"

namespace AMDiS {

  ElInfo* RefinementManager2d::refineFunction(ElInfo* el_info)
  {
    FUNCNAME("RefinementManager::refineFunction()");

    int n_neigh;
    bool bound = false;

    DegreeOfFreedom *edge[2];

    RCNeighbourList* refineList = NEW RCNeighbourList(2);

    if (el_info->getElement()->getMark() <= 0)  
      return(el_info);    /*   element may not be refined   */

    refineList->setElement(0, el_info->getElement());
    n_neigh = 1;

    if (el_info->getProjection(0) && 
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION) {
      newCoords = true;
    }
  
    /****************************************************************************/
    /*  give the refinement edge the right orientation                          */
    /****************************************************************************/

    if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0)) {
      edge[0] = const_cast<int*>( el_info->getElement()->getDOF(0));
      edge[1] = const_cast<int*>( el_info->getElement()->getDOF(1));
    } else {
      edge[1] = const_cast<int*>( el_info->getElement()->getDOF(0));
      edge[0] = const_cast<int*>( el_info->getElement()->getDOF(1));
    }

    /****************************************************************************/
    /*  get the refinement patch                                                */
    /****************************************************************************/
    if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) {
      /****************************************************************************/
      /*  domain's boundary was reached while looping around the refinement edge  */
      /****************************************************************************/
      getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
      bound = true;
    }
    
    // ==========================================================================
    // === check for periodic boundary ==========================================
    // ==========================================================================

    DegreeOfFreedom *next_edge[2];
    DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
    DegreeOfFreedom *last_edge[2] = {NULL, NULL};
    int n_neigh_periodic;

    DegreeOfFreedom newDOF = -1;
    DegreeOfFreedom lastNewDOF = -1;
    DegreeOfFreedom firstNewDOF = -1;

    RCNeighbourList *periodicList;
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();

    while (edge[0] != NULL) {
      periodicList = refineList->periodicSplit(edge, 
					       next_edge,
					       &n_neigh,
					       &n_neigh_periodic);

      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");

      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);

      if (firstNewDOF == -1) {
	firstNewDOF = newDOF;
      }
      
      if (lastNewDOF != -1) {
	for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	  if (it->second) {
	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
		 (*(it->second))[edge[1][0]] == last_edge[1][0]) || 
		((*(it->second))[edge[0][0]] == last_edge[1][0] &&
		 (*(it->second))[edge[1][0]] == last_edge[0][0]))
	      {
		(*(it->second))[lastNewDOF] = newDOF;
		(*(it->second))[newDOF] = lastNewDOF;
	      } 
	  }
	}
      }
      lastNewDOF = newDOF;

      last_edge[0] = edge[0];
      last_edge[1] = edge[1];

      edge[0] = next_edge[0];
      edge[1] = next_edge[1];
    }

    if (lastNewDOF != firstNewDOF) {
      for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	if (it->second) {
	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
	       (*(it->second))[first_edge[1][0]] == last_edge[1][0]) || 
	      ((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
	       (*(it->second))[first_edge[1][0]] == last_edge[0][0]))
	    {
	      (*(it->second))[lastNewDOF] = firstNewDOF;
	      (*(it->second))[firstNewDOF] = lastNewDOF;
	    }
	}
      }
    }
  

    /****************************************************************************/
    /*  and now refine the patch                                                */
    /****************************************************************************/

    DELETE refineList;

    return(el_info);
  }


  int RefinementManager2d::newCoordsFct(ElInfo *el_info)
  {
    Element *el = el_info->getElement();
    int dow = Global::getGeo(WORLD);

    Projection *projector = el_info->getProjection(0);

    if (!projector || projector->getType() != VOLUME_PROJECTION) {
      projector = el_info->getProjection(2);
    }

    if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
      WorldVector<double> *new_coord = NEW WorldVector<double>;
      
      for (int j = 0; j < dow; j++)
	(*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j])*0.5;
      
      projector->project(*new_coord);

      el->setNewCoord(new_coord);
    }
    
    return 0;
  }

  void RefinementManager2d::setNewCoords()
  {
    Flag fillFlag = Mesh::CALL_EVERY_EL_PREORDER|
      Mesh::FILL_BOUND|
      Mesh::FILL_COORDS;
  
    mesh->traverse(-1, fillFlag, newCoordsFct);
  }


  DegreeOfFreedom RefinementManager2d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
    FUNCNAME("RefinementManager2d::refinePatch()");
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};
    Triangle        *el = dynamic_cast<Triangle*>(const_cast<Element*>( refineList->getElement(0)));
    Triangle        *neigh = dynamic_cast<Triangle*>(const_cast<Element*>( refineList->getElement(1)));

    /****************************************************************************/
    /*  there is one new vertex in the refinement edge                          */
    /****************************************************************************/

    dof[0] = mesh->getDOF(VERTEX);

    mesh->incrementNumberOfVertices(1);
    mesh->incrementNumberOfEdges(1);

    if (mesh->getNumberOfDOFs(EDGE)) {
      /****************************************************************************/
      /*  there are two additional dofs in the refinement edge                    */
      /****************************************************************************/
      dof[1] = mesh->getDOF(EDGE);
      dof[2] = mesh->getDOF(EDGE);
    }

    /****************************************************************************/
    /*  first refine the element                                                */
    /****************************************************************************/

    bisectTriangle(el, dof);

    if (neigh) {
      DegreeOfFreedom *tmp = dof[1];
      /****************************************************************************/
      /*  there is a neighbour; refine it also, but first exchange dof[1] and     */
      /*  dof[2]; thus, dof[1] is always added on child[0]!                       */
      /****************************************************************************/
      dof[1] = dof[2];
      dof[2] = tmp;

      bisectTriangle(neigh, dof);
    } else {
      newCoords = true;
    }

    /****************************************************************************/
    /*  if there are functions to interpolate data to the finer grid, do so     */
    /****************************************************************************/
  
    int iadmin;
    int nrAdmin = mesh->getNumberOfDOFAdmin();
    for(iadmin = 0; iadmin < nrAdmin; iadmin++) {
      std::list<DOFIndexedBase*>::iterator it;
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
      for(it = admin->beginDOFIndexed(); it != end; it++) {
	(*it)->refineInterpol(*refineList, n_neigh);
      }
    }


    if (!mesh->queryCoarseDOFs()) {
      /****************************************************************************/
      /*  if there should be no dof information on interior leaf elements remove  */
      /*  dofs from edges and the centers of parents                              */
      /****************************************************************************/
      if (mesh->getNumberOfDOFs(EDGE)) {
	
	int node = mesh->getNode(EDGE);
	
	/****************************************************************************/
	/*  the only DOF that can be freed is that in the refinement edge; all other*/
	/*  DOFs are handed on the children                                         */
	/****************************************************************************/
	
	mesh->freeDOF(const_cast<int*>( el->getDOF(node+2)), EDGE);
      }
      if (mesh->getNumberOfDOFs(EDGE)  ||  mesh->getNumberOfDOFs(CENTER)) {
	refineList->removeDOFParents(n_neigh);
      }
    }

    return dof[0][0];
  }


  void RefinementManager2d::bisectTriangle(Triangle *el, 
					   DegreeOfFreedom* newDOFs[3])
  {
    FUNCNAME("RefinementManager2d::bisectTriangle()");
    TEST_EXIT_DBG(mesh)("no mesh!\n");
  
    Triangle *child[2];

    child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
  
    signed char newMark = max(0, el->getMark()-1);
    child[0]->setMark(newMark);
    child[1]->setMark(newMark);
    el->setMark(0);

    /****************************************************************************/
    /*  transfer hidden data from parent to children                            */
    /****************************************************************************/

    el->refineElementData(child[0], child[1]); // call of subclass-method
    el->setFirstChild(child[0]);
    el->setSecondChild(child[1]);

    if (newMark > 0) 
      doMoreRecursiveRefine = true;

    /****************************************************************************/
    /*  vertex 2 is the newest vertex                                           */
    /****************************************************************************/

    child[0]->setDOF(2, newDOFs[0]);
    child[1]->setDOF(2, newDOFs[0]);

    /****************************************************************************/
    /*  the other vertices are handed on from the parent                        */
    /****************************************************************************/

    for (int i_child = 0; i_child < 2; i_child++) {
      child[i_child]->setDOF(i_child, const_cast<int*>( el->getDOF(2)));
      child[i_child]->setDOF(1-i_child, const_cast<int*>( el->getDOF(i_child)));
    }

    /****************************************************************************/
    /*  there is one more leaf element, two hierachical elements and one more   */
    /*  edge                                        			    */
    /****************************************************************************/

    mesh->incrementNumberOfEdges(1);
    mesh->incrementNumberOfLeaves(1);
    mesh->incrementNumberOfElements(2);

    if (mesh->getNumberOfDOFs(EDGE)) {
      /****************************************************************************/
      /*  there are dof's in the midpoint of the edges			    */
      /****************************************************************************/
      DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);
      
      child[0]->setDOF(4, newEdgeDOFs);
      child[1]->setDOF(3, newEdgeDOFs);
      
      /****************************************************************************/
      /*  dofs handed on by the parent                                            */
      /****************************************************************************/
      child[0]->setDOF(5, const_cast<int*>( el->getDOF(4)));
      child[1]->setDOF(5, const_cast<int*>( el->getDOF(3)));
      
      /****************************************************************************/
      /*  dofs in the refinement edge                                             */
      /****************************************************************************/
      child[0]->setDOF(3, newDOFs[1]);
      child[1]->setDOF(4, newDOFs[2]);
    }
    
    if (mesh->getNumberOfDOFs(CENTER)) {
      int node = mesh->getNode(CENTER);
      
      /****************************************************************************/
      /* there are dofs at the barycenter of the triangles                        */
      /****************************************************************************/
      child[0]->setDOF(node, mesh->getDOF(CENTER));
      child[1]->setDOF(node, mesh->getDOF(CENTER));
    }
  }

  int RefinementManager2d::getRefinePatch(ElInfo **el_info, 
					  DegreeOfFreedom *edge[2],
					  int dir, RCNeighbourList* refineList, 
					  int *n_neigh)
  {
    FUNCNAME("RefinementManager2d::getRefinePatch()");

    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( (*el_info)->getElement()));
    int opp_vertex = 0;

    if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
      /****************************************************************************/
      /*  neighbour is not compatible devisible; refine neighbour first, store the*/
      /*  opp_vertex to traverse back to el                                       */
      /****************************************************************************/
      opp_vertex = (*el_info)->getOppVertex(2);
      
      ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
      neigh_info = refineFunction(neigh_info);
      
      /****************************************************************************/
      /*  now go back to the original element and refine the patch                */
      /****************************************************************************/
      *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
      TEST_EXIT_DBG(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
    }
  
    if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
      TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }

    return(0);
  }

}