Traverse.cc

#include "Traverse.h"
#include "Element.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "ElInfo.h"
#include "Mesh.h"
#include "Flag.h"
#include "MacroElement.h"
#include "FixVec.h"
#include "Parametric.h"
#include "OpenMP.h"

namespace AMDiS {

  ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level,
				       Flag fill_flag)
  {
    FUNCNAME("TraverseStack::traverseFirst()");

    traverse_mesh = mesh;
    traverse_level = level;
    traverse_fill_flag = fill_flag; 

    if (stack_size < 1)
      enlargeTraverseStack();

    Flag FILL_ANY = Mesh::getFillAnyFlag(mesh->getDim());

    for (int i = 0; i < stack_size; i++)
      elinfo_stack[i]->setFillFlag(fill_flag & FILL_ANY);


    elinfo_stack[0]->setMesh(mesh);
    elinfo_stack[1]->setMesh(mesh);

    if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
      TEST_EXIT_DBG(level >= 0)("invalid level: %d\n", level);
    }

    traverse_mel = NULL;
    stack_used = 0;

    return(traverseNext(NULL));  
  }

  ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old)
  {
    FUNCNAME("TraverseStack::traverseNext()");

    ElInfo *elinfo = NULL;
    ElInfo::traverseId[omp_get_thread_num()] = id;
    Parametric *parametric = traverse_mesh->getParametric();

    if (stack_used) {
      if (parametric) 
	elinfo_old = parametric->removeParametricInfo(elinfo_old);

      TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
    }
    else {
      TEST_EXIT_DBG(elinfo_old == NULL)("invalid old elinfo != nil\n");
    }

    if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL)) {
      elinfo = traverseLeafElement();
    } else {
      if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL))
	elinfo = traverseLeafElementLevel();
      else if (traverse_fill_flag.isSet(Mesh::CALL_EL_LEVEL))
	elinfo = traverseElementLevel();
      else if (traverse_fill_flag.isSet(Mesh::CALL_MG_LEVEL))
	elinfo = traverseMultiGridLevel();
      else
	if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER))
	  elinfo = traverseEveryElementPreorder();
	else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_INORDER))
	  elinfo = traverseEveryElementInorder();
	else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_POSTORDER))
	  elinfo = traverseEveryElementPostorder();
	else
	  ERROR_EXIT("invalid traverse_flag\n");
    }

    if (elinfo) {
      if (parametric) {
	elinfo = parametric->addParametricInfo(elinfo);
      }
      elinfo->fillDetGrdLambda();
    }

    return(elinfo);
  }



  /****************************************************************************/
  /* common (static) traversal routines for 2d and 3d                         */
  /* to be #included in traverse_r.c                                          */
  /****************************************************************************/

  /* traverse hierarchical mesh,  call el_fct() for each/some element */

  /****************************************************************************/
  /*   recursive_traverse:                                		    */
  /*   -------------------                                		    */
  /*   recursively traverse the mesh hierarchy tree       		    */
  /*   call the routine traverse_info->el_fct(el_info) with    		    */
  /*    - all tree leaves                                 		    */
  /*    - all tree leaves at a special level              		    */
  /*    - all tree elements in pre-/in-/post-order        		    */
  /*   depending on the traverse_info->level variable     		    */
  /****************************************************************************/
  int Traverse::recursive(ElInfoStack *elInfoStack)
  {
    FUNCNAME("Traverse::recursive()");

    ElInfo *elinfo = elInfoStack->getCurrentElement();

    Element *el = elinfo->getElement();
    int mg_level, sum = 0;
    Parametric *parametric = mesh->getParametric();
    ElInfo::traverseId[omp_get_thread_num()] = id;

    if (flag.isSet(Mesh::CALL_LEAF_EL)) {
      if (el->getFirstChild()) {
	ElInfo* elinfo_new = elInfoStack->getNextElement();
	elinfo_new->fillElInfo(0, elinfo);
	sum += recursive(elInfoStack);
	elinfo_new->fillElInfo(1, elinfo);
	sum += recursive(elInfoStack);
	elInfoStack->getBackElement();
      } else {
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      }
      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
    }

    if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
      if (el->getFirstChild()) {
	if ((elinfo->getLevel() >= level)) {
	  return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
	}
	ElInfo* elinfo_new = elInfoStack->getNextElement();
	elinfo_new->fillElInfo(0, elinfo);
	sum += recursive(elInfoStack);
	elinfo->fillElInfo(1, elinfo);
	sum += recursive(elInfoStack);
	elInfoStack->getBackElement();
      } else {
	if ((elinfo->getLevel() == level) && (el_fct!=NULL)) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      }
      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
    }

    if (flag.isSet(Mesh::CALL_EL_LEVEL)) {
      if (elinfo->getLevel() == level) {
	if (NULL != el_fct) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo); 
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      } else {
	if (elinfo->getLevel() > level){
	  return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
	} else if (el->getFirstChild()) {
	  ElInfo* elinfo_new = elInfoStack->getNextElement();
	  elinfo_new->fillElInfo(0, elinfo);
	  sum += recursive(elInfoStack);
	  elinfo_new->fillElInfo(1, elinfo);
	  sum += recursive(elInfoStack);
	  elInfoStack->getBackElement();
	}
      }

      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
    }

    if (flag.isSet(Mesh::CALL_MG_LEVEL)) {

      mg_level = (elinfo->getLevel() + mesh->getDim() - 1) / mesh->getDim();

      if (mg_level > level)  {
	return 0;
      }

      if (!(el->getFirstChild())) {
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
	return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
      }
    
      if ((mg_level == level) && ((elinfo->getLevel() % mesh->getDim()) == 0)) {
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo); 
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
	return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
      }

      ElInfo* elinfo_new = elInfoStack->getNextElement();
 
      elinfo_new->fillElInfo(0, elinfo);
      sum += recursive(elInfoStack);

      elinfo_new->fillElInfo(1, elinfo);
      sum += recursive(elInfoStack);

      elInfoStack->getBackElement();

      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;   
    }

    if (flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) {
      if (el_fct != NULL) {
	if (parametric) 
	  elinfo = parametric->addParametricInfo(elinfo);
	elinfo->fillDetGrdLambda();
	sum += el_fct(elinfo);
	if (parametric) 
	  elinfo = parametric->removeParametricInfo(elinfo);
      }
    }

    if (el->getFirstChild()) {
      ElInfo* elinfo_new = elInfoStack->getNextElement();
   
      elinfo_new->fillElInfo(0, elinfo);
      sum += recursive(elInfoStack);

      if (flag.isSet(Mesh::CALL_EVERY_EL_INORDER)) 
	if (el_fct!=NULL) { 
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      elinfo_new->fillElInfo(1, elinfo);
      sum += recursive(elInfoStack);

      elInfoStack->getBackElement();
    } else {
      if (flag.isSet(Mesh::CALL_EVERY_EL_INORDER)) 
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
    }

    if (flag.isSet(Mesh::CALL_EVERY_EL_POSTORDER)) 
      if (el_fct != NULL) {
	if (parametric) 
	  elinfo = parametric->addParametricInfo(elinfo);
	elinfo->fillDetGrdLambda();
	sum += el_fct(elinfo);
	if (parametric) 
	  elinfo = parametric->removeParametricInfo(elinfo);
      }

    return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
  }

  void TraverseStack::enlargeTraverseStack()
  {
    FUNCNAME("TraverseStack::enlargeTraverseStack()");

    int new_stack_size = stack_size + 10;

    elinfo_stack.resize(new_stack_size, NULL);
 
    // create new elinfos
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(elinfo_stack[i] == NULL)("???\n");
      elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }

    if (stack_size > 0) {
      for (int i = stack_size; i < new_stack_size; i++) {
	elinfo_stack[i]->setFillFlag(elinfo_stack[0]->getFillFlag());
      }
    }

    info_stack.resize(new_stack_size);
    save_elinfo_stack.resize(new_stack_size, NULL);

    // create new elinfos
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(save_elinfo_stack[i] == NULL)("???\n");
      save_elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }
    save_info_stack.resize(new_stack_size);

    stack_size = new_stack_size;    
  }

  
  ElInfo* TraverseStack::traverseLeafElement()
  {
    FUNCNAME("TraverseStack::traverseLeafElement()");

    Element *el = NULL;

    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();

      while (((*currentMacro)->getIndex() % maxThreads_ != myThreadId_) &&
      	     (currentMacro != traverse_mesh->endOfMacroElements())) {
      	currentMacro++;
      }

      if (currentMacro == traverse_mesh->endOfMacroElements()) {
	return NULL;
      }

      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

      el = elinfo_stack[stack_used]->getElement();
      if ((el == NULL) || (el->getFirstChild() == NULL)) {
	return (elinfo_stack[stack_used]);
      }
    } else {
      el = elinfo_stack[stack_used]->getElement();
      
      /* go up in tree until we can go down again */
      while ((stack_used > 0) && 
	     ((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
      
      /* goto next macro element */
      if (stack_used < 1) {
	do {	
	  currentMacro++;
	} while ((currentMacro != traverse_mesh->endOfMacroElements()) && 
		 ((*currentMacro)->getIndex() % maxThreads_ != myThreadId_));

	if (currentMacro == traverse_mesh->endOfMacroElements()) {
	  return NULL;
	}
	traverse_mel = *currentMacro;
	
	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
	info_stack[stack_used] = 0;
	
	el = elinfo_stack[stack_used]->getElement();
	if ((el == NULL) || (el->getFirstChild() == NULL)) {
	  return (elinfo_stack[stack_used]);
	}
      }
    }
   
    /* go down tree until leaf */
    while (el->getFirstChild()) {
      if (stack_used >= stack_size - 1) 
	enlargeTraverseStack();
      
      int i = info_stack[stack_used];
      el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
      info_stack[stack_used]++;
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
      stack_used++;

      TEST_EXIT_DBG(stack_used < stack_size)
	("stack_size=%d too small, level=(%d,%d)\n",
	 stack_size, elinfo_stack[stack_used]->getLevel());
      
      info_stack[stack_used] = 0;
    }

    return elinfo_stack[stack_used];
  }

  ElInfo* TraverseStack::traverseLeafElementLevel()
  {
    FUNCNAME("TraverseStack::traverseLeafElementLevel()");
    ERROR_EXIT("not yet implemented\n");

    return NULL;
  }

  ElInfo* TraverseStack::traverseElementLevel()
  {
    FUNCNAME("TraverseStack::traverseElementLevel()");

    ElInfo *elInfo;
    do {
      elInfo = traverseEveryElementPreorder();
    } while (elInfo != NULL && elInfo->getLevel() != traverse_level);

    return elInfo;
  }

  ElInfo* TraverseStack::traverseMultiGridLevel()
  {
    FUNCNAME("TraverseStack::traverseMultiGridLevel()");

    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      traverse_mel = *currentMacro;
      if (traverse_mel == NULL)  return NULL;
      
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
      
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
	return(elinfo_stack[stack_used]);
    }
  
    Element *el = elinfo_stack[stack_used]->getElement();

    /* go up in tree until we can go down again */
    while ((stack_used > 0) && 
	   ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL))) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return(NULL);

      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
	return(elinfo_stack[stack_used]);
    }


    /* go down tree */

    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
    info_stack[stack_used]++;

    elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
    stack_used++;

    TEST_EXIT_DBG(stack_used < stack_size)
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
    if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	(elinfo_stack[stack_used]->getLevel() < traverse_level && 
	 elinfo_stack[stack_used]->getElement()->isLeaf()))
      return(elinfo_stack[stack_used]);

    return traverseMultiGridLevel();
  }

  ElInfo* TraverseStack::traverseEveryElementPreorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementPreorder()");

    Element *el;

    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      traverse_mel = *currentMacro;
      if (traverse_mel == NULL)  
	return NULL;
      
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

      return(elinfo_stack[stack_used]);
    }
  
    el = elinfo_stack[stack_used]->getElement();

    /* go up in tree until we can go down again */
    while((stack_used > 0) && 
	  ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL)))
      {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return(NULL);
      traverse_mel = *currentMacro;

      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

      return(elinfo_stack[stack_used]);
    }


    /* go down tree */

    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
    info_stack[stack_used]++;

    elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
    stack_used++;

    TEST_EXIT_DBG(stack_used < stack_size)
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
    return(elinfo_stack[stack_used]);
  }

  ElInfo* TraverseStack::traverseEveryElementInorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementInorder");
    ERROR_EXIT("not yet implemented\n");
    return NULL;
  }

  ElInfo* TraverseStack::traverseEveryElementPostorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementPostorder");

    Element *el;
    int i;

    if (stack_used == 0)   /* first call */
      {
	currentMacro = traverse_mesh->firstMacroElement();
	if(currentMacro == traverse_mesh->endOfMacroElements()) return NULL;
	traverse_mel = *currentMacro;

	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
	info_stack[stack_used] = 0;

	//return(elinfo_stack[stack_used]);
      }
    else { /* don't go up on first call */
  
      el = elinfo_stack[stack_used]->getElement();

      /* go up in tree until we can go down again */          /* postorder!!! */
      while((stack_used > 0) && 
	    ((info_stack[stack_used] >= 3) || (el->getFirstChild()==NULL)))
	{
	  stack_used--;
	  el = elinfo_stack[stack_used]->getElement();
	}


      /* goto next macro element */
      if (stack_used < 1) {
	currentMacro++;
	if(currentMacro == traverse_mesh->endOfMacroElements()) return NULL;
	traverse_mel = *currentMacro;

	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
	info_stack[stack_used] = 0;

	/*    return(elinfo_stack+stack_used); */
      }
    }
    /* go down tree */

    while(elinfo_stack[stack_used]->getElement()->getFirstChild()
	  && (info_stack[stack_used] < 2)) {
      if (stack_used >= stack_size-1) 
	enlargeTraverseStack();

      i = info_stack[stack_used];
      info_stack[stack_used]++;
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
      stack_used++;
      info_stack[stack_used] = 0;
    }
  
    info_stack[stack_used]++;      /* postorder!!! */

    return(elinfo_stack[stack_used]);
  }

  ElInfo* TraverseStack::traverseNeighbour(ElInfo* elinfo_old, int neighbour)
  {
    int dim = elinfo_old->getMesh()->getDim();
    switch(dim) {
    case 1:
      ERROR_EXIT("invalid dim\n");
      break;
    case 2:
      return traverseNeighbour2d(elinfo_old, neighbour);
      break;
    case 3:
      return traverseNeighbour3d(elinfo_old, neighbour);
      break;
    default: ERROR_EXIT("invalid dim\n");
    }
    return NULL;
  }

  ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
  {
    FUNCNAME("TraverseStackLLtraverseNeighbour3d()");

    Element *el2;
    ElInfo *old_elinfo, *elinfo, *elinfo2;
    const DegreeOfFreedom *dof;
    int i, nb, opp_vertex, stack2_used, typ = 0;
    int sav_index, sav_neighbour = neighbour;

    static int coarse_nb[3][3][4] = {{{-2,-2,-2,-2}, {-1,2,3,1}, {-1,3,2,0}},
				     {{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}},
				     {{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}}};

    TEST_EXIT_DBG(stack_used > 0)("no current element\n");

    Parametric *parametric = traverse_mesh->getParametric();

    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);

    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])
      ("invalid old elinfo\n");

    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
      ("FILL_NEIGH not set");
    Element *el = elinfo_stack[stack_used]->getElement();
    sav_index = el->getIndex();

    /* first, goto to leaf level, if necessary... */
    if ((traverse_fill_flag & Mesh::CALL_LEAF_EL).isAnySet()) {
      if ((el->getChild(0)) && (neighbour < 2)) {
	if (stack_used >= stack_size - 1)
	  enlargeTraverseStack();
	i = 1 - neighbour;
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
	info_stack[stack_used] = i + 1;
	stack_used++;
	info_stack[stack_used] = 0;
	neighbour = 3;
      }
    }


    /* save information about current element and its position in the tree */
    save_traverse_mel = traverse_mel;
    save_stack_used = stack_used;

    nb = neighbour;

    while (stack_used > 1) { /* go up in tree until we can go down again */
      stack_used--;
      typ = dynamic_cast<ElInfo3d*>( elinfo_stack[stack_used])->getType();
      nb = coarse_nb[typ][info_stack[stack_used]][nb];
      if (nb == -1) 
	break;
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }

    /* save hierarchy information about current element */

    for (i = stack_used; i <= save_stack_used; i++) {
      save_info_stack[i] = info_stack[i];
      *(save_elinfo_stack[i]) = *(elinfo_stack[i]);
    }

    old_elinfo = save_elinfo_stack[save_stack_used];
    opp_vertex = old_elinfo->getOppVertex(neighbour);

    if (nb >= 0) {                           /* go to macro element neighbour */

      i = traverse_mel->getOppVertex(nb);
      traverse_mel = traverse_mel->getNeighbour(nb);
      if (traverse_mel == NULL)  return(NULL);
    
      if ((nb < 2) && (save_stack_used > 1)) {
	stack2_used = 2;                /* go down one level in OLD hierarchy */
      }
      else {
	stack2_used = 1;
      }
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
      nb = i;
    } else {                                                /* goto other child */
      stack2_used = stack_used + 1;
      if (save_stack_used > stack2_used) {
	stack2_used++;                  /* go down one level in OLD hierarchy */
      }
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

      if (stack_used >= stack_size-1)
	enlargeTraverseStack();
      i = 2 - info_stack[stack_used];
      info_stack[stack_used] = i+1;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
      stack_used++;
      info_stack[stack_used] = 0;
      nb = 0;

    }


    elinfo = elinfo_stack[stack_used];
    el = elinfo->getElement();

    while(el->getChild(0)) {
 
      if (nb < 2) {                         /* go down one level in hierarchy */
	if (stack_used >= stack_size-1)
	  enlargeTraverseStack();
	info_stack[stack_used] = 2-nb;
	elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
	stack_used++;
	info_stack[stack_used] = 0;
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
	nb = 3;
      }

      if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
	TEST_EXIT_DBG(el->getChild(0))("invalid NEW refinement?\n");

	if (el->getDOF(0) == el2->getDOF(0))
	  i = save_info_stack[stack2_used] - 1;
	else if (el->getDOF(1) == el2->getDOF(0))
	  i = 2 - save_info_stack[stack2_used];
	else {
	  if (traverse_mesh->associated(el->getDOF(0, 0), el2->getDOF(0, 0)))
	    i = save_info_stack[stack2_used] - 1;
	  else if (traverse_mesh->associated(el->getDOF(1, 0), el2->getDOF(0, 0)))
	    i = 2 - save_info_stack[stack2_used];
	  else {
	    ERROR_EXIT("no common refinement edge\n");
	  }
	}

	if (el->getChild(0)  &&  
	    (el->getChild(i)->getDOF(1) == el->getDOF(nb) ||
	     traverse_mesh->associated(el->getChild(i)->getDOF(1, 0), el->getDOF(nb, 0)))) 
	  {   /*  el->child[0]  Kuni  22.08.96  */
	    nb = 1;
	  }
	else {
	  nb = 2;
	}

	info_stack[stack_used] = i+1;
	if (stack_used >= stack_size-1)
	  enlargeTraverseStack();
	elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
	stack_used++;
	info_stack[stack_used] = 0;

	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();

	stack2_used++;
	elinfo2 = save_elinfo_stack[stack2_used];
	el2 = elinfo2->getElement();
	if (save_stack_used > stack2_used) {
	  dof = el2->getDOF(1);
	  if (dof != el->getDOF(1) && dof != el->getDOF(2) &&
	      !traverse_mesh->associated(dof[0], el->getDOF(1, 0)) &&
	      !traverse_mesh->associated(dof[0], el->getDOF(2, 0))) {
	  
	    stack2_used++;               /* go down one level in OLD hierarchy */
	    elinfo2 = save_elinfo_stack[stack2_used];
	    el2 = elinfo2->getElement();
	  } 
	}

      }
      else {   /* now we're done... */
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
	break;
      }
    }


    if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
      MSG(" looking for neighbour %d of element %d at %p\n",
	  neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
      MSG(" originally: neighbour %d of element %d at %p\n",
	  sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
      MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
	  elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
	  opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
	("didn't succeed !?!?!?\n");
    }


    if ((traverse_fill_flag & Mesh::CALL_EVERY_EL_POSTORDER).isAnySet())
      info_stack[stack_used] = 3;
    else if ((traverse_fill_flag & Mesh::CALL_EVERY_EL_INORDER).isAnySet())
      info_stack[stack_used] = 1;  /* ??? */

    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);
      elinfo->fillDetGrdLambda();
    }
    return(elinfo);
  }

  ElInfo* TraverseStack::traverseNeighbour2d(ElInfo* elinfo_old, int neighbour)
  {
    FUNCNAME("TraverseStack::traverseNeighbour2d()");

    Element *el, *el2, *sav_el;
    ElInfo *old_elinfo, *elinfo, *elinfo2;
    int i, nb, opp_vertex, stack2_used;

    static int coarse_nb[3][3] = {{-2,-2,-2}, {2,-1,1}, {-1,2,0}};
    /* father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j] */

    int sav_index, sav_neighbour = neighbour;

    TEST_EXIT_DBG(stack_used > 0)("no current element");
    TEST_EXIT_DBG(traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL))("invalid traverse_fill_flag");

    Parametric *parametric = traverse_mesh->getParametric();

    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);

    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");

    elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
    el = elinfo_stack[stack_used]->getElement();
    sav_index = el->getIndex();
    sav_el = el;

    /* first, goto to leaf level, if necessary... */
    if (!(el->isLeaf()) && (neighbour < 2)) {
      if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
      i = 1 - neighbour;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
      info_stack[stack_used] = i + 1;
      stack_used++;
      neighbour = 2;
    }

    /* save information about current element and its position in the tree */
    save_traverse_mel = traverse_mel;
    save_stack_used = stack_used;
    for (i = 0; i<=stack_used; i++)
      save_info_stack[i] = info_stack[i];
    for (i = 0; i<=stack_used; i++)
      (*(save_elinfo_stack[i])) = (*(elinfo_stack[i]));
    old_elinfo = save_elinfo_stack[stack_used];
    opp_vertex = old_elinfo->getOppVertex(neighbour);


    /****************************************************************************/
    /* First phase: go up in tree until we can go down again.                   */
    /*                                                                          */
    /* During this first phase, nb is the neighbour index which points from an  */
    /* element of the OLD hierarchy branch to the NEW branch                    */
    /****************************************************************************/

    nb = neighbour;

    while (stack_used > 1) {
      stack_used--;
      nb = coarse_nb[info_stack[stack_used]][nb];
      if (nb == -1) 
	break;
      
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }

    /****************************************************************************/
    /* Now, goto neighbouring element at the local hierarchy entry              */
    /* This is either a macro element neighbour or the other child of parent.   */
    /* initialize nb for second phase (see below)                               */
    /****************************************************************************/

    if (nb >= 0) {                        /* go to macro element neighbour */

      if ((nb < 2) && (save_stack_used > 1)) {
	stack2_used = 2;           /* go down one level in OLD hierarchy */
      } else {
	stack2_used = 1;
      }
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

      i = traverse_mel->getOppVertex(nb);
      traverse_mel = traverse_mel->getNeighbour(nb);
      if (traverse_mel == NULL)
	return NULL;
      nb = i;
    
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
      info_stack[stack_used] = 0;
    
    } else {                                               /* goto other child */

      stack2_used = stack_used + 1;
      if (save_stack_used > stack2_used) {
	stack2_used++;               /* go down one level in OLD hierarchy */
      }
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();


      if (stack_used >= stack_size - 1) {
	enlargeTraverseStack();
      }
      i = 2 - info_stack[stack_used];
      info_stack[stack_used] = i + 1;
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
      stack_used++;
      nb = 1 - i;
    }

    /****************************************************************************/
    /* Second phase: go down in a new hierarchy branch until leaf level.        */
    /* Now, nb is the neighbour index which points from an element of the       */
    /* NEW hierarchy branch to the OLD branch.                                  */
    /****************************************************************************/

    elinfo = elinfo_stack[stack_used];
    el = elinfo->getElement();

    while (el->getFirstChild()) {

      if (nb < 2) {   /* go down one level in hierarchy */
	if (stack_used >= stack_size - 1)
	  enlargeTraverseStack();
	elinfo_stack[stack_used + 1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
	info_stack[stack_used] = 2 - nb;
	stack_used++;
	nb = 2;
      }

      if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
	TEST_EXIT_DBG(el->getFirstChild())("invalid NEW refinement?");

	/* use child i, neighbour of el2->child[nb-1] */
	i = 2 - save_info_stack[stack2_used];
	TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
	info_stack[stack_used] = i + 1;
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
	stack_used++;
	nb = i;

	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();

	stack2_used++;
	if (save_stack_used > stack2_used) {
	  stack2_used++;                /* go down one level in OLD hierarchy */
	}
	elinfo2 = save_elinfo_stack[stack2_used];
	el2 = elinfo2->getElement();

      } else {   /* now we're done... */
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
      }
    }


    if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
      MSG(" looking for neighbour %d of element %d at %8X\n",
	  neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
      MSG(" originally: neighbour %d of element %d at %8X\n",
	  sav_neighbour, sav_index, sav_el);
      MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
	  elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
	  elinfo->getNeighbour(opp_vertex)?elinfo->getNeighbour(opp_vertex)->getIndex():-1);
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
	("didn't succeed !?!?!?");
    }

    if (elinfo->getElement()->getFirstChild()) {
      MSG(" looking for neighbour %d of element %d at %8X\n",
	  neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
      MSG(" originally: neighbour %d of element %d at %8X\n",
	  sav_neighbour, sav_index, sav_el);
      MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
	  elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
	  elinfo->getNeighbour(opp_vertex)->getIndex());
      MSG("got no leaf element\n");
      WAIT_REALLY;
    }

    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);

      elinfo->fillDetGrdLambda();
    }

    return elinfo;
  }

  void TraverseStack::update()
  {
    FUNCNAME("TraverseStack::update()");

    TEST_EXIT_DBG(traverse_mesh->getDim() == 3)("update only in 3d\n");

    for (int i = stack_used; i > 0; i--) {
      dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
    }
  }

  void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo, 
				      DimMat<double> *coords)
  {
    int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo->getLevel();

    upperElInfo->getRefSimplexCoords(coords);
    
    for (int i = 1; i <= levelDif; i++) {
      upperElInfo->getSubElementCoords(coords,
				       elinfo_stack[stack_used - levelDif + i]->getIChild());
    }
  }
}