#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 "DOFVector.h"
#include "PeriodicBC.h"
#include "VertexVector.h"

namespace AMDiS {

  void RefinementManager3d::bisectTetrahedron(RCNeighbourList* ref_list, 
					      int index,
					      DegreeOfFreedom* dof[3],
					      DegreeOfFreedom *edge[2])
  {
    Tetrahedron *el = dynamic_cast<Tetrahedron*>(const_cast<Element*>( ref_list->getElement(index))), *child[2];
    int          i, node;
    int          el_type = ref_list->getType(index);

    child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
  
    int mark = max(0, el->getMark()-1);
    child[0]->setMark(mark);
    child[1]->setMark(mark);
    el->setMark(0);

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

    el->refineElementData(child[0], child[1], el_type);

    el->setFirstChild(child[0]);
    el->setSecondChild(child[1]);

    if (child[0]->getMark() > 0)  doMoreRecursiveRefine = true;

    int n_vertices = mesh->getGeo(VERTEX);
    child[0]->setDOF(n_vertices-1, dof[0]);
    child[1]->setDOF(n_vertices-1, dof[0]);

    for (i = 0; i < n_vertices-1; i++)
      {
	child[0]->
	  setDOF(i, const_cast<int*>( el->getDOF(Tetrahedron::childVertex[el_type][0][i])));
	child[1]->
	  setDOF(i, const_cast<int*>( el->getDOF(Tetrahedron::childVertex[el_type][1][i])));
      }
    /****************************************************************************/
    /*  there is one more leaf element and two more hierachical elements        */
    /****************************************************************************/

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

    /****************************************************************************/
    /* first set those dof pointers for higher order without neighbour          */
    /* information                                                              */
    /****************************************************************************/

    if (mesh->getNumberOfDOFs(EDGE))
      {
	node = mesh->getNode(EDGE);

	/****************************************************************************/
	/*  set pointers to those dof's that are handed on from the parant          */
	/****************************************************************************/

	child[0]->
	  setDOF(node, 
		 const_cast<int*>( el->getDOF(node+Tetrahedron::childEdge[el_type][0][0])));
	child[1]->
	  setDOF(node, 
		 const_cast<int*>( el->getDOF(node+Tetrahedron::childEdge[el_type][1][0])));
	child[0]->
	  setDOF(node+1, 
		 const_cast<int*>( el->getDOF(node+Tetrahedron::childEdge[el_type][0][1])));
	child[1]->
	  setDOF(node+1, 
		 const_cast<int*>( el->getDOF(node+Tetrahedron::childEdge[el_type][1][1])));
	child[0]->
	  setDOF(node+3, 
		 const_cast<int*>( el->getDOF(node+Tetrahedron::childEdge[el_type][0][3])));
	child[1]->
	  setDOF(node+3, 
		 const_cast<int*>( el->getDOF(node+Tetrahedron::childEdge[el_type][1][3])));

	/****************************************************************************/
	/*  adjust pointers to the dof's in the refinement edge                     */
	/****************************************************************************/

	if (el->getDOF(0) == edge[0])
	  {
	    child[0]->setDOF(node+2, dof[1]);
	    child[1]->setDOF(node+2, dof[2]);
	  }
	else
	  {
	    child[0]->setDOF(node+2, dof[2]);
	    child[1]->setDOF(node+2, dof[1]);
	  }
      }

    if (mesh->getNumberOfDOFs(FACE))
      {
	node = mesh->getNode(FACE);    

	/****************************************************************************/
	/*  set pointers to those dof's that are handed on from the parant          */
	/****************************************************************************/

	child[0]->setDOF(node+3, const_cast<int*>( el->getDOF(node+1)));
	child[1]->setDOF(node+3, const_cast<int*>( el->getDOF(node+0)));

	/****************************************************************************/
	/*  get new dof for the common face of child0 and child1                    */
	/****************************************************************************/

	DegreeOfFreedom *newDOF = mesh->getDOF(FACE);
	child[0]->setDOF(node, static_cast<int*>( newDOF));
	child[1]->setDOF(node, static_cast<int*>( newDOF));
      }
  
    if (mesh->getNumberOfDOFs(CENTER))
      {
	node = mesh->getNode(CENTER);
	child[0]->setDOF(node, const_cast<int*>( mesh->getDOF(CENTER)));
	child[1]->setDOF(node, const_cast<int*>( mesh->getDOF(CENTER)));
      }

    if (mesh->getNumberOfDOFs(EDGE)  || mesh->getNumberOfDOFs(FACE))
      fillPatchConnectivity(ref_list, index);

    //   MSG("%d -> %d %d\n", 
    //       el->getIndex(),
    //       child[0]->getIndex(),
    //       child[1]->getIndex());
  }

  void RefinementManager3d::fillPatchConnectivity(RCNeighbourList* ref_list, 
						  int index)
  {
    FUNCNAME("RefinementManager3d::fillPatchConnectivity");
    Element *el = ref_list->getElement(index), *neigh;
    int     dir, el_type = ref_list->getType(index), n_type = 0;
    int     adjc, i, j, i_neigh, j_neigh;
    int     node0, node1, opp_v = 0;

    for (dir = 0; dir < 2; dir++)
      {
	neigh = ref_list->getNeighbourElement(index, dir);
	if (neigh)
	  {
	    n_type = ref_list->getType( ref_list->getNeighbourNr(index, dir) );
	    opp_v = ref_list->getOppVertex(index, dir);
	  }

	if (!neigh  ||  neigh->isLeaf())
	  {

	    /****************************************************************************/
	    /*  get new dof's in the midedge of the face of el and for the two midpoints*/
	    /*  of the sub-faces. If face is an interior face those pointers have to be */
	    /*  adjusted by the neighbour element also (see below)                      */
	    /****************************************************************************/

	    if (mesh->getNumberOfDOFs(EDGE))
	      {
		node0 = node1 = mesh->getNode(EDGE);
		node0 += Tetrahedron::nChildEdge[el_type][0][dir];
		node1 += Tetrahedron::nChildEdge[el_type][1][dir];
		DegreeOfFreedom *newDOF = mesh->getDOF(EDGE);
		(const_cast<Element*>( el->getFirstChild()))->setDOF(node0, newDOF);
		(const_cast<Element*>( el->getSecondChild()))->setDOF(node1, newDOF);
	      }
	    if (mesh->getNumberOfDOFs(FACE))
	      {
		node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
		(const_cast<Element*>( el->getFirstChild()))->setDOF(node0, mesh->getDOF(FACE));
		node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
		(const_cast<Element*>( el->getSecondChild()))->setDOF(node1, mesh->getDOF(FACE));
	      }
	  }
	else     /*   if (!neigh  ||  !neigh->child[0])                         */
	  {
	    /****************************************************************************/
	    /*  interior face and neighbour has been refined, look for position at the  */
	    /*  refinement edge                                                         */
	    /****************************************************************************/
      
	    if (el->getDOF(0) == neigh->getDOF(0))
	      {
		/****************************************************************************/
		/* same position at refinement edge                                         */
		/****************************************************************************/
		adjc = 0;
	      }
	    else
	      {
		/****************************************************************************/
		/* different position at refinement edge                                    */
		/****************************************************************************/
		adjc = 1;
	      }

	    for (i = 0; i < 2; i++)
	      {
		j = Tetrahedron::adjacentChild[adjc][i];

		i_neigh = Tetrahedron::nChildFace[el_type][i][dir];
		j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v-2];

		/****************************************************************************/
		/*  adjust dof pointer in the edge in the common face of el and neigh and   */
		/*  the dof pointer in the sub-face child_i-child_j (allocated by neigh!)   */
		/****************************************************************************/

		if (mesh->getNumberOfDOFs(EDGE))
		  {
		    node0 = 
		      mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir];
		    node1 = 
		      mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v-2];

		    TEST_EXIT(neigh->getChild(j)->getDOF(node1))
		      ("no dof on neighbour %d at node %d\n", 
		       neigh->getChild(j)->getIndex(), node1);

		    (const_cast<Element*>( el->getChild(i)))->
		      setDOF(node0, const_cast<int*>( neigh->getChild(j)->getDOF(node1)));
		  }
		if (mesh->getNumberOfDOFs(FACE))
		  {
		    node0 = mesh->getNode(FACE) + i_neigh;
		    node1 = mesh->getNode(FACE) + j_neigh;

		    TEST_EXIT(neigh->getChild(j)->getDOF(node1))
		      ("no dof on neighbour %d at node %d\n",
		       neigh->getChild(j)->getIndex(), node1);

		    (const_cast<Element*>( el->getChild(i)))->
		      setDOF(node0, const_cast<int*>( neigh->getChild(j)->getDOF(node1)));
		  }

	      }  /*   for (i = 0; i < 2; i++)                                       */
	  }    /*   else of   if (!neigh  ||  !neigh->child[0])                   */
      }      /*   for (dir = 0; dir < 2; dir++)                                 */
  }


  int RefinementManager3d::newCoordsFct(ElInfo *el_info)
  {
    Element *el = el_info->getElement();
    int i, j, n_neigh;
    DegreeOfFreedom *edge[2];
    ElInfo *elinfo =  el_info;

    int dow = Global::getGeo(WORLD);

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

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

    if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
      WorldVector<double> *new_coord = NEW WorldVector<double>;

      for (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);
      /****************************************************************************/
      /*  now, information should be passed on to patch neighbours...             */
      /*  get the refinement patch                                                */
      /****************************************************************************/
      refList->setElement(0, el);
      refList->setElType(0, dynamic_cast<ElInfo3d*>(el_info)->getType());
      n_neigh = 1;
      for (i = 0; i < 2; i++)
	edge[i] = const_cast<int*>( el_info->getElement()->getDOF(i));
      if (getRefinePatch(&elinfo, edge, 0, refList, &n_neigh))
	{
	  /****************************************************************************/
	  /*  domain's boundary was reached while looping around the refinement edge  */
	  /****************************************************************************/
	  getRefinePatch(&elinfo, edge, 1, refList, &n_neigh);
	}
      for (i = 1; i < n_neigh; i++)            /* start with 1, as list[0]=el */
	{
	  TEST(!refList->getElement(i)->isNewCoordSet())
	    ("non-nil new_coord in el %d ref_list[%d] el %d (n_neigh=%d)\n",
	     el->getIndex(), i, refList->getElement(i)->getIndex(), n_neigh);

	  refList->getElement(i)->setNewCoord(el->getNewCoord());
	}
    }

    return 0;
  }

  void RefinementManager3d::setNewCoords()
  {
    ElInfo *el_info;

    Flag fillFlag = Mesh::CALL_EVERY_EL_PREORDER|
      Mesh::FILL_BOUND|
      Mesh::FILL_COORDS;

    if (refList)
      DELETE refList;
 
    refList = NEW RCNeighbourList(mesh->getMaxEdgeNeigh());

    fillFlag |= Mesh::FILL_NEIGH;
    el_info = stack->traverseFirst(mesh, -1, fillFlag);
    while (el_info) {
      newCoordsFct(el_info);
      el_info = stack->traverseNext(el_info);
    }
  }


  DegreeOfFreedom RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
    int    i;
    Tetrahedron *el = dynamic_cast<Tetrahedron*>(const_cast<Element*>( refineList->getElement(0)));
    /* first element in the list */
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};

    /****************************************************************************/
    /*  get new dof's in the refinement edge                                    */
    /****************************************************************************/

    dof[0] = mesh->getDOF(VERTEX);
    mesh->incrementNumberOfVertices(1);
  
    if (mesh->getNumberOfDOFs(EDGE))
      {
	dof[1] = mesh->getDOF(EDGE);
	dof[2] = mesh->getDOF(EDGE);
      }

    for (i = 0; i < n_neigh; i++) {
      bisectTetrahedron(refineList, i, dof, edge);
    }

    /****************************************************************************/
    /*  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, faces and the centers of parents                       */
	/****************************************************************************/
	if (mesh->getNumberOfDOFs(EDGE))
	  {
	    /****************************************************************************/
	    /*  remove dof of the midpoint of the common refinement edge                */
	    /****************************************************************************/

	    el = dynamic_cast<Tetrahedron*>(const_cast<Element*>( refineList->getElement(0)));
	    mesh->freeDOF(const_cast<int*>( el->getDOF(mesh->getNode(EDGE))), EDGE);
	  }

	if (mesh->getNumberOfDOFs(EDGE)  ||  
	    mesh->getNumberOfDOFs(FACE)  ||  
	    mesh->getNumberOfDOFs(CENTER))
	  {
	    for (i = 0; i < n_neigh; i++)
	      refineList->removeDOFParent(i);
	  }
      }


    /****************************************************************************/
    /*  update the number of edges and faces; depends whether refinement edge   */
    /*  is a boundary edge or not                                               */
    /****************************************************************************/

    if (bound)
      {
	mesh->incrementNumberOfEdges(n_neigh + 2);
	mesh->incrementNumberOfFaces(2*n_neigh + 1);
	newCoords = true;
      }
    else
      {
	mesh->incrementNumberOfEdges(n_neigh + 1);
	mesh->incrementNumberOfFaces(2*n_neigh);
      }

    return dof[0][0];
  }


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

    Tetrahedron *el, *neigh;
    int          i, j, k, opp_v = 0;
    int          edge_no, neigh_el_type;

    el = dynamic_cast<Tetrahedron*>(const_cast<Element*>( (*el_info)->getElement()));
  
    if ((*el_info)->getNeighbour(3-dir) == NULL)  {return(1);}
  
    opp_v = (*el_info)->getOppVertex(3-dir);
    neigh_info = stack->traverseNeighbour3d((*el_info), 3-dir);
    neigh_el_type = dynamic_cast<ElInfo3d*>(neigh_info)->getType();

    neigh = dynamic_cast<Tetrahedron*>(const_cast<Element*>( neigh_info->getElement()));
  
    int vertices = mesh->getGeo(VERTEX);

    while (neigh != el) {
      for (j = 0; j < vertices; j++)
	if (neigh->getDOF(j) == edge[0])  break;
      for (k = 0; k < vertices; k++)
	if (neigh->getDOF(k) == edge[1])  break;


      if(j > 3 || k > 3) {
	for (j = 0; j < vertices; j++)
	  if (mesh->associated(neigh->getDOF(j, 0), edge[0][0]))  break;
	for (k = 0; k < vertices; k++)
	  if (mesh->associated(neigh->getDOF(k, 0), edge[1][0]))  break;
	    
	if(j > 3 || k > 3) {
	  for (j = 0; j < vertices; j++)
	    if (mesh->indirectlyAssociated(neigh->getDOF(j, 0), edge[0][0]))  break;
	  for (k = 0; k < vertices; k++)
	    if (mesh->indirectlyAssociated(neigh->getDOF(k, 0), edge[1][0]))  break;
	    
	  TEST_EXIT(j < vertices && k < vertices)
	    ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
	     edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0),
	     neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0));

	}
      }

      //       LeafDataPeriodic *pd = 
      // 	dynamic_cast<LeafDataPeriodic*>(el->getElementData(PERIODIC));

      //       if(pd) {
      // 	::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
      // 	::std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = 
      // 	  pd->getInfoList().end();
	  
      // 	for(it = pd->getInfoList().begin(); it != end; ++it) {
      // 	  if(it->periodicMode != 0) {
      // 	    PeriodicBC *periodicCondition = mesh->getPeriodicBCMap()[it->type];
      // 	    if(periodicCondition) {
      // 	      VertexVector *associated = mesh->getPeriodicAssociations()[it->type];
      // 	      for (j = 0; j < vertices; j++)
      // 		if (neigh->getDOF(j, 0) == (*associated)[edge[0][0]])  break;
      // 	      for (k = 0; k < vertices; k++)
      // 		if (neigh->getDOF(k, 0) == (*associated)[edge[1][0]])  break;
	      
      // 	      TEST_EXIT(j < vertices && k < vertices)
      // 		("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
      // 		 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0),
      // 		 neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0));

      // 	    } else {
      // 	      ERROR_EXIT("???\n");
      // 	    }
      // 	  }
      // 	}
      //       }

      TEST_EXIT(j <  mesh->getGeo(VERTEX) &&
		k < mesh->getGeo(VERTEX))
	("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
	 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0),
	 neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0));

      if ((edge_no = Tetrahedron::edgeOfDOFs[j][k]))
	{
	  /****************************************************************************/
	  /*  neigh has not a compatible refinement edge                              */
	  /****************************************************************************/
	  neigh->setMark(max(neigh->getMark(), 1));
	  neigh_info = refineFunction(neigh_info);

	  /****************************************************************************/
	  /*  now, go to a child at the edge and determine the opposite vertex for    */
	  /*  this child; continue the looping around the edge with this element      */
	  /****************************************************************************/

	  neigh_info = stack->traverseNext(neigh_info);
	  neigh_el_type = dynamic_cast<ElInfo3d*>(neigh_info)->getType();

	  switch (edge_no)
	    {
	    case 1: 
	      opp_v = opp_v == 1 ? 3 : 2;
	      break;
	    case 2: 
	      opp_v = opp_v == 2 ? 1 : 3;
	      break;
	    case 3: 
	      neigh_info = stack->traverseNext(neigh_info);
	      neigh_el_type = dynamic_cast<ElInfo3d*>(neigh_info)->getType();
	      if (neigh_el_type != 1)
		opp_v = opp_v == 0 ? 3 : 2;
	      else
		opp_v = opp_v == 0 ? 3 : 1;
	      break;
	    case 4:
	      neigh_info = stack->traverseNext(neigh_info);
	      neigh_el_type = dynamic_cast<ElInfo3d*>(neigh_info)->getType();
	      if (neigh_el_type != 1)
		opp_v = opp_v == 0 ? 3 : 1;
	      else
		opp_v = opp_v == 0 ? 3 : 2;
	      break;
	    case 5:
	      if (neigh_el_type != 1)
		{
		  neigh_info = stack->traverseNext(neigh_info);
		  neigh_el_type = (dynamic_cast<ElInfo3d*>(neigh_info))->getType();
		}
	      opp_v = 3;
	      break;
	    }
	  neigh = dynamic_cast<Tetrahedron*>(const_cast<Element*>( neigh_info->getElement()));
	}
      else
	{
	  /****************************************************************************/
	  /*  neigh is compatible devisible; put neigh to the list of patch elements; */
	  /*  go to next neighbour                                                    */
	  /****************************************************************************/
	  TEST_EXIT(*n_neigh < mesh->getMaxEdgeNeigh())
	    ("too many neighbours %d in refine patch\n", *n_neigh);

      
	  refineList->setElement(*n_neigh, neigh);
	  refineList->setElType(*n_neigh, neigh_el_type);

	  /****************************************************************************/
	  /*  we have to go back to the starting element via opp_v values             */
	  /*  correct information is produce by get_neigh_on_patch()                  */
	  /****************************************************************************/
	  refineList->setOppVertex(*n_neigh, 0, opp_v); 

	  ++*n_neigh;

	  if (opp_v != 3)
	    i = 3;
	  else
	    i = 2;

	  if (neigh_info->getNeighbour(i))
	    {
	      opp_v = neigh_info->getOppVertex(i);
	      neigh_info = stack->traverseNeighbour3d(neigh_info, i);
	      neigh_el_type = (dynamic_cast<ElInfo3d*>(neigh_info))->getType();
	      neigh = dynamic_cast<Tetrahedron*>(const_cast<Element*>( neigh_info->getElement()));
	    }
	  else
	    break;
	}
    }

    if (neigh == el)
      {
	(*el_info) = neigh_info;
	return(0);
      }

    /****************************************************************************/
    /*  the domain's boundary is reached; loop back to the starting el          */
    /****************************************************************************/

    i = *n_neigh-1;
    opp_v = refineList->getOppVertex(i, 0);
    do
      {
	TEST_EXIT(neigh_info->getNeighbour(opp_v)  &&  i > 0)
	  ("while looping back domains boundary was reached or i == 0\n");
	opp_v = refineList->getOppVertex(i--, 0);
	neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
      } while(neigh_info->getElement() != el);
    (*el_info) = neigh_info;

    return(1);
  }


  ElInfo* RefinementManager3d::refineFunction(ElInfo* el_info)
  {
    FUNCNAME("RefinementManager3d::refineFunction()");
    int              n_neigh, bound = false;
    DegreeOfFreedom  *edge[2];
    RCNeighbourList  *ref_list;

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

    /****************************************************************************/
    /*  get memory for a list of all elements at the refinement edge            */
    /****************************************************************************/
    ref_list = NEW RCNeighbourList(mesh->getMaxEdgeNeigh());
    ref_list->setElType(0, (dynamic_cast<ElInfo3d*>(el_info))->getType());
    ref_list->setElement(0, el_info->getElement());
    n_neigh = 1;

    /****************************************************************************/
    /*  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                                                */
    /****************************************************************************/
    //   MSG("index %d\n", el_info->getElement()->getIndex());
    if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh))
      {
	/****************************************************************************/
	/*  domain's boundary was reached while looping around the refinement edge  */
	/****************************************************************************/
	getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh);
	bound = true;
      }

    /****************************************************************************/
    /*  fill neighbour information inside the patch in the refinement list      */
    /****************************************************************************/
    ref_list->getNeighOnPatch(n_neigh, bound);

    // ==========================================================================
    // === 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();

    //   for(int i = 0; i < n_neigh; i++) {
    //     MSG("%d ", ref_list->getElement(i)->getIndex());
    //   }
    //   MSG("\n");

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

      TEST_EXIT(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                                                */
    /****************************************************************************/
    //refinePatch(edge, ref_list, n_neigh, bound);

    stack->update();

    DELETE ref_list;

    return(el_info);
  }

}