#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
#include "MacroElement.h"
#include "MacroReader.h"
#include "Mesh.h"
#include "Traverse.h"
#include "Parameters.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "MacroWriter.h"
#include "PeriodicMap.h"
#include "Projection.h"
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

namespace AMDiS {

#define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC)

  //**************************************************************************
  //  flags, which information should be present in the elInfo structure     
  //**************************************************************************

  const Flag Mesh::FILL_NOTHING    = 0X00L;
  const Flag Mesh::FILL_COORDS     = 0X01L;
  const Flag Mesh::FILL_BOUND      = 0X02L;
  const Flag Mesh::FILL_NEIGH      = 0X04L;
  const Flag Mesh::FILL_OPP_COORDS = 0X08L;
  const Flag Mesh::FILL_ORIENTATION= 0X10L;
  const Flag Mesh::FILL_DET        = 0X20L;
  const Flag Mesh::FILL_GRD_LAMBDA = 0X40L;
  const Flag Mesh::FILL_ADD_ALL    = 0X80L;


  const Flag Mesh::FILL_ANY_1D= (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
  const Flag Mesh::FILL_ANY_2D= (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
  const Flag Mesh::FILL_ANY_3D= (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L);

  //**************************************************************************
  //  flags for Mesh traversal                                                
  //**************************************************************************

  const Flag Mesh::CALL_EVERY_EL_PREORDER  = 0X0100L;
  const Flag Mesh::CALL_EVERY_EL_INORDER   = 0X0200L;
  const Flag Mesh::CALL_EVERY_EL_POSTORDER = 0X0400L;
  const Flag Mesh::CALL_LEAF_EL            = 0X0800L;
  const Flag Mesh::CALL_LEAF_EL_LEVEL      = 0X1000L;
  const Flag Mesh::CALL_EL_LEVEL           = 0X2000L;
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L ; // used in mg methods 


  // const Flag Mesh::USE_PARAMETRIC          = 0X8000L ; // used in mg methods 

  // ::std::list<Mesh*> Mesh::meshes;
  DOFAdmin* Mesh::compressAdmin = NULL;
  Mesh* Mesh::traversePtr = NULL;
  int Mesh::iadmin = 0;
  ::std::vector<DegreeOfFreedom> Mesh::dof_used;
  const int Mesh::MAX_DOF=100;
  ::std::map<DegreeOfFreedom, DegreeOfFreedom*> Mesh::serializedDOFs;

  struct delmem { 
    DegreeOfFreedom* ptr;
    int              len;
  };


  Mesh::Mesh(const ::std::string& aName, int dimension) 
    : name(aName), 
      dim(dimension), 
      nVertices(0),
      nEdges(0),
      nLeaves(0), 
      nElements(0),
      parametric(NULL), 
      preserveCoarseDOFs(false),
      nDOFEl(0),
      nDOF(dimension, DEFAULT_VALUE, 0),
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
      elementPrototype(NULL),
      elementDataPrototype(NULL),
      elementIndex(-1),
      initialized(false),
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

    FUNCNAME("Mesh::Mesh");

    // set default element prototype
    switch(dim) {
    case 1:
      elementPrototype = NEW Line(this);
      break;
    case 2:
      elementPrototype = NEW Triangle(this);
      break;
    case 3:
      elementPrototype = NEW Tetrahedron(this);
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

    elementIndex=0;
  };

  Mesh::~Mesh()
  {
  };

  void Mesh::addMacroElement(MacroElement* m) {
    macroElements.push_back(m); 
    m->setIndex(macroElements.size());
  };




  int Mesh::traverse(int level, Flag flag, 
		     int (*el_fct)(ElInfo*))
  {
    FUNCNAME("Mesh::traverse()");
    ::std::deque<MacroElement*>::iterator mel;
    ElInfo* elinfo = createNewElInfo();
    Traverse tinfo(this, flag, level, el_fct);
    int sum = 0;
  
    elinfo->setFillFlag(flag);
  
    if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL) || 
	flag.isSet(Mesh::CALL_EL_LEVEL)      || 
	flag.isSet(Mesh::CALL_MG_LEVEL)) {
      TEST(level >= 0)("invalid level: %d\n", level);
    }
  
    for (mel = macroElements.begin(); mel != macroElements.end(); mel++) {
      elinfo->fillMacroInfo(*mel);
      sum += tinfo.recursive(elinfo);
    }

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



  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    int i, j, d, n;
    ::std::vector<DOFAdmin*>::iterator dai;

    localAdmin->setMesh(this);
    n = admin.size();

    dai=::std::find(admin.begin(),admin.end(),localAdmin);
    if (dai!= admin.end()) {
      ERROR("admin %s is already associated to mesh %s\n",
	    localAdmin->getName().c_str(), this->getName().c_str());
    }

    // ===== adding dofs to already existing elements ============================ 

    if (initialized) {
      static bool pnd_1d_0[2] = {true, true};
      static bool pnd_1d_1[1] = {false};
      static bool pnd_2d_0[3] = {true, true, true};
      static bool pnd_2d_1[3] = {true, true, false};
      static bool pnd_2d_2[1] = {false};
      static bool pnd_3d_0[4] = {true, true, true, true};
      static bool pnd_3d_1[6] = {false, true, true, true, true, true};
      static bool pnd_3d_2[4] = {true, true, false, false};
      static bool pnd_3d_3[1] = {false};
      static bool *pnd_1d[2] = {pnd_1d_0, pnd_1d_1};
      static bool *pnd_2d[3] = {pnd_2d_0, pnd_2d_1, pnd_2d_2};
      static bool *pnd_3d[4] = {pnd_3d_0, pnd_3d_1, pnd_3d_2, pnd_3d_3};
      static bool **parentNeedsDOF[4] = {NULL, pnd_1d, pnd_2d, pnd_3d};

     
      ::std::list<struct delmem> delList;
      ::std::map< ::std::set<DegreeOfFreedom>, DegreeOfFreedom*> dofPtrMap;
      const DOFAdmin *vertexAdmin = getVertexAdmin();
      int vertexAdminPreDOFs = vertexAdmin->getNumberOfPreDOFs(VERTEX);

      // finding necessary node number for new admin

      int newNNode=0;
      GeoIndex geoIndex;

      for(d = 0; d < dim+1; d++) {
	geoIndex = INDEX_OF_DIM(d, dim);
      
	if (localAdmin->getNumberOfDOFs(geoIndex)>0||nDOF[geoIndex]>0)
	  newNNode+=getGeo(geoIndex);
      };

      bool extendNodes=(newNNode>nNodeEl);
      int  oldNNodes=nNodeEl;

      nNodeEl=newNNode;

      TraverseStack stack;
      ElInfo *elInfo = NULL;
    
      WARNING("You are using untested code (adding dofs to existing mesh). Please contact\nsoftware administrator if any errors occur in this context.\n");

      elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
      while(elInfo) {
	Element *element = elInfo->getElement();
	DegreeOfFreedom *newDOF, **oldDOF, **dof = 
	  const_cast<DegreeOfFreedom**>(element->getDOF());

	int index = 0;

	if (extendNodes) {
	  oldDOF=dof;
	  element->setDOFPtrs();
	  dof=const_cast<DegreeOfFreedom**>(element->getDOF());
	  int index=0,oldIndex=0;
	  for(d = 0; d < dim+1; d++) {
	    geoIndex = INDEX_OF_DIM(d, dim);
	    if (nDOF[geoIndex]>0) {
	      for(i=0;i<getGeo(geoIndex);++i) 
		dof[index++]=oldDOF[oldIndex++];
	    }
	    else {
	      if (localAdmin->getNumberOfDOFs(geoIndex)>0) 
		index+=getGeo(geoIndex);
	    }
	  }
	
	  FREE_MEMORY(oldDOF, DegreeOfFreedom*, oldNNodes);

	  TEST_EXIT(index==nNodeEl)("ERROR: Number of entered nodes %f != number of nodes %f\n",index,nNodeEl);

	}


	index=0;

	// allocate new memory at elements
	for(d = 0; d < dim+1; d++) {
	  geoIndex = INDEX_OF_DIM(d, dim);
      
	  int numberOfDOFs = localAdmin->getNumberOfDOFs(geoIndex);
	  int numberOfPreDOFs = nDOF[geoIndex];

	  if (numberOfDOFs>0||numberOfPreDOFs>0) {

	    // for all vertices/edges/...
	    for(i = 0; i < getGeo(geoIndex); i++, index++) {
	      ::std::set<DegreeOfFreedom> dofSet;
	      for(j = 0; j < d+1; j++) {
		dofSet.insert(dof[element->getVertexOfPosition(geoIndex, i, j)][vertexAdminPreDOFs]);
	      }
	    
	      if(element->isLeaf() || parentNeedsDOF[dim][d][i]) {
		if(dofPtrMap[dofSet] == NULL) {
		  if(localAdmin->getNumberOfDOFs(geoIndex)) {
		    newDOF = GET_MEMORY(DegreeOfFreedom, numberOfPreDOFs + numberOfDOFs);
		    // copy old dofs to new memory and free old memory
		    if(dof[index]) {
		      for(j = 0; j < numberOfPreDOFs; j++) {
			newDOF[j] = dof[index][j];
		      }
		      //	  FREE_MEMORY(dof[index], DegreeOfFreedom, numberOfPreDOFs);
		      // Do not free memory. The information has to be used to identify the part in other elements.
		      // The memory is only marked for freeing.
		      struct delmem fm;
		      fm.ptr=dof[index];
		      fm.len=numberOfPreDOFs;
		      delList.push_back(fm);
		    }
		    for(j = 0; j < numberOfDOFs; j++) {
		      newDOF[numberOfPreDOFs + j] = localAdmin->getDOFIndex();
		    }
		    dof[index] = newDOF;
		  }
		  dofPtrMap[dofSet] = dof[index];
		} else {
		  dof[index] = dofPtrMap[dofSet];
		}
	      }
	    }
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
  
      // now free the old dof memory:

      ::std::list<struct delmem>::iterator it=delList.begin();
    
      while(it!=delList.end()) {
	FREE_MEMORY((*it).ptr, DegreeOfFreedom, (*it).len);
	it++;
      }

      delList.clear();

    }
    // ============================================================================

    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
    nDOF[VERTEX]  += localAdmin->getNumberOfDOFs(VERTEX);
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

    if(dim > 1) {
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
      nDOF[EDGE]    += localAdmin->getNumberOfDOFs(EDGE);
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

    localAdmin->setNumberOfPreDOFs(CENTER,nDOF[CENTER]);
    nDOF[CENTER]  += localAdmin->getNumberOfDOFs(CENTER);
    nDOFEl += nDOF[CENTER];

    TEST_EXIT(nDOF[VERTEX] > 0)("no vertex dofs\n");

    node[VERTEX]  = 0;
    nNodeEl     = getGeo(VERTEX);

    if(dim > 1) {
      node[EDGE]    = nNodeEl;
      if (nDOF[EDGE] > 0) nNodeEl += getGeo(EDGE);
    }

    if (3==dim){
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
      nDOF[FACE]  += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl     += getGeo(FACE) * nDOF[FACE];
      node[FACE]    = nNodeEl;
      if (nDOF[FACE] > 0) nNodeEl +=  getGeo(FACE);
    }

    node[CENTER]    = nNodeEl;
    if (nDOF[CENTER] > 0) nNodeEl += 1;

    return;
  }


  /****************************************************************************/
  /*  dofCompress: remove holes in dof vectors                                */
  /****************************************************************************/

  void Mesh::dofCompress()
  {
    FUNCNAME("Mesh::dofCompress");
    int            size;
    Flag           fill_flag;

    //   int i;
    //   for(i = 0; i < periodicAssociations[-2]->getSize(); i++) {
    //     MSG("asso %d : dof %d -> dof %d\n", 
    // 	periodicAssociations[-2], 
    // 	i, 
    // 	(*periodicAssociations[-2])[i]);
    //   }

    for (iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++)
      {
	TEST_EXIT((compressAdmin = admin[iadmin]))
	  ("no admin[%d] in mesh\n", iadmin);

	if ((size = compressAdmin->getSize()) < 1) continue;
	if (compressAdmin->getUsedDOFs() < 1)    continue;
	if (compressAdmin->getHoleCount() < 1)    continue;
    
	newDOF.resize(size);

	compressAdmin->compress(newDOF);
    
	if (preserveCoarseDOFs) {
	  fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
	}
	else {
	  fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
	}

	traverse( -1, fill_flag, newDOFFct1);
	traverse( -1, fill_flag, newDOFFct2);

	newDOF.resize(0);
      }

    //   MSG("\n");
    //   for(i = 0; i < periodicAssociations[-2]->getSize(); i++) {
    //     MSG("asso %d : dof %d -> dof %d\n", 
    // 	periodicAssociations[-2], 
    // 	i, 
    // 	(*periodicAssociations[-2])[i]);
    //   }

    return;
  }

  // /****************************************************************************/
  // /*  fill_boundary:                                                          */
  // /*  fills boundary information for the vertices of the macro triangulation  */
  // /*    mel_vertex[k][i] = global index of vertex k of element i              */
  // /****************************************************************************/

  // void Mesh::fillBoundary(int **mel_vertex)
  // {
  //   FUNCNAME("Mesh::fillBoundary");
  //   ::std::deque<MacroElement*>::const_iterator  mel = firstMacroElement();
  //   int         i, j,k;
  //   int       lne = getNumberOfLeaves(), lnv = getNumberOfVertices();
  //   int       ned = getNumberOfEdges();

  //   BoundaryType   *bound = GET_MEMORY(BoundaryType, lnv+(3==dim)?ned:0);

  //   for (i = 0; i < lnv; i++)
  //     bound[i] = INTERIOR;

  //   int facesPlusEdges = 
  //     (*mel)->getElement()->getGeo(FACE) +
  //     (*mel)->getElement()->getGeo(EDGE);

  //   for (i = 0; i < lne; i++) {
  //     for (k = 0; k < getGeo(NEIGH); k++) {
  //       //if ((*(mel+i))->getBoundary(k)) {
  //       if ((*(mel+i))->getBoundary(k) >= DIRICHLET) {
  // 	for (j = 1; j < dim+1; j++)
  // 	  bound[mel_vertex[i][(k+j)%(dim+1)]] = DIRICHLET;
  //       }
  //       else if ((*(mel+i))->getBoundary(k) <= NEUMANN) {
  // 	for (j = 1; j < dim+1; j++)
  // 	  if (bound[mel_vertex[i][(k+j)%(dim+1)]] != DIRICHLET)
  // 	    bound[mel_vertex[i][(k+j)%(dim+1)]] = NEUMANN;
  //       }
  //       //}
  //     }
  //   }

  //   for (i = 0; i < lne; i++)
  //     for (k = 0; k < getGeo(VERTEX); k++) {
  //       (*(mel+i))->setBoundary(facesPlusEdges + k, 
  // 			      bound[mel_vertex[i][k]]);
  //     }
  //   FREE_MEMORY(bound, BoundaryType, lnv+(3==dim)?ned:0);

  //   return;
  // };

  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
    FUNCNAME("Mesh::getDOF");
    const DOFAdmin *localAdmin;
    DegreeOfFreedom *dof;
    int     i, j, n, n0, ndof;

    TEST_EXIT(position >= CENTER && position <= FACE)
      ("unknown position %d\n",position);

    ndof = getNumberOfDOFs(position);
    if (ndof <= 0) return(NULL);

    dof = GET_MEMORY(DegreeOfFreedom, ndof);

    for (i = 0; i < getNumberOfDOFAdmin(); i++)
      {
	localAdmin = &getDOFAdmin(i);
	TEST_EXIT(localAdmin)("no admin[%d]\n", i);

	n  = localAdmin->getNumberOfDOFs(position);
	n0 = localAdmin->getNumberOfPreDOFs(position);

	TEST_EXIT(n+n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);

	for (j = 0; j < n; j++)
	  {
	    dof[n0+j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
	  }
      }
  
    return(dof);
  }


  // const Boundary *Mesh::defaultBoundary(int bound)
  // {
  //   static const Boundary dn[2] = {DIRICHLET,  NEUMANN};
  //   if (bound >= DIRICHLET)
  //     return(dn);
  //   else if (bound <= NEUMANN)
  //     return(dn+1);
  //   else 
  //     return(NULL);
  // }

  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
    FUNCNAME("Mesh::createDOFPtrs");
    int i;
    DegreeOfFreedom **ptrs;

    if (nNodeEl <= 0)
      return(NULL);

    ptrs = GET_MEMORY(DegreeOfFreedom*, nNodeEl);
    for (i = 0; i < nNodeEl; i++)
      ptrs[i] = NULL;

    return(ptrs);
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
    FUNCNAME("Mesh::freeDOFPtrs");

    TEST_EXIT(ptrs)("ptrs=NULL\n");

    if (nNodeEl <= 0)
      return;
  
    FREE_MEMORY(ptrs, DegreeOfFreedom*, nNodeEl);
  }


  const DOFAdmin *Mesh::createDOFAdmin(const ::std::string& lname,DimVec<int> lnDOF)
  {
    FUNCNAME("Mesh::createDOFAdmin");

    DOFAdmin         *localAdmin;
    int              i;

    localAdmin=NEW DOFAdmin(this,lname);

    for (i = 0; i < dim+1; i++)
      localAdmin->setNumberOfDOFs(i,lnDOF[i]);

    addDOFAdmin(localAdmin);

    return(localAdmin);
  }





  // int Mesh::macroType(const ::std::string& filename, const ::std::string& type)
  // {
  //   const char *fn, *t;

  //   if (3==dim) return 0;
  
  //   if (filename.size() <= type.size())
  //     return(false);

  //   fn = filename.data();
  //   while (*fn) fn++;
  //   t = type.data();
  //   while (*t) t++;

  //   while (t != type  &&  *t == *fn) t--;
  
  //   return(t == type);
  // }

  const DOFAdmin* Mesh::getVertexAdmin() const
  {
    int       i;
    const DOFAdmin *localAdmin = NULL;

    for (i = 0; i < static_cast<int>(admin.size()); i++)
      {
	if (admin[i]->getNumberOfDOFs(VERTEX))
	  {
	    if (!localAdmin)  
	      localAdmin = admin[i];
	    else if (admin[i]->getSize() < localAdmin->getSize())
	      localAdmin = admin[i];
	  }
      }
    return(localAdmin);
  }

  void Mesh::freeDOF(DegreeOfFreedom* dof, GeoIndex position)
  {
    FUNCNAME("Mesh::freeDOF");
    DOFAdmin *localAdmin;
    int     i, j, n, n0, ndof;

    TEST_EXIT(position >= CENTER && position <= FACE)
      ("unknown position %d\n",position);

    ndof = nDOF[position];
    if (ndof) 
      {
	if (!dof)
	  {
	    MSG("dof = NULL, but ndof=%d\n", ndof);
	    return;
	  }
      }
    else
      {
	if (dof)
	  {
	    MSG("dof != NULL, but ndof=0\n");
	  }
	return;
      }

    TEST_EXIT(ndof <= MAX_DOF)
      ("ndof too big: ndof=%d, MAX_DOF=%d\n",ndof,MAX_DOF);

    for (i = 0; i < static_cast<int>(admin.size()); i++)
      {
	localAdmin = admin[i];

	n  = localAdmin->getNumberOfDOFs(position);
	n0 = localAdmin->getNumberOfPreDOFs(position);

	TEST_EXIT(n+n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);

	for (j = 0; j < n; j++)
	  {
	    localAdmin->freeDOFIndex(dof[n0+j]);
	  }
      }

    FREE_MEMORY(dof, DegreeOfFreedom, ndof);
    return;  
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
    DELETE el;
  }




  //const Boundary* defaultBoundary(Mesh* m,int i) {return m->defaultBoundary(i);};

  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
    TEST_EXIT(elementPrototype)("no element prototype\n");

    Element *el = parent ? parent->clone() : elementPrototype->clone();

    //el->setIndex(elementIndex++);
    //el->setDOFPtrs();
  
    //   ElementData *elementData = NULL;

    //   if(parent && parent->getElementData()) {
    //     elementData = parent->getElementData()->clone();
    //   } else if(elementDataPrototype) {
    //     elementData = elementDataPrototype->clone();
    //   }
  
    if(!parent && elementDataPrototype) {
      el->setElementData(elementDataPrototype->clone()); 
    } else {
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
    }

    //   if(!elementData) {
    //     WARNING("no element data!\n");
    //   }


    return el;
  }

  ElInfo* Mesh::createNewElInfo()
  {
    switch(dim) {
    case 1:
      return NEW ElInfo1d(this);
      break;
    case 2:
      return NEW ElInfo2d(this);
      break;
    case 3:
      return NEW ElInfo3d(this);
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
    };
  }



  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
			       DimVec<double>&    bary,
			       const MacroElement      *start_mel,
			       const WorldVector<double> *xy0,
			       double            *sp)
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
    ElInfo *mel_info = NULL;
    int  i, k;
    bool inside;

    mel_info = createNewElInfo();

    if (start_mel != NULL)
      mel = start_mel;
    else
      if((mel == NULL)||(mel->getElement()->getMesh() != this))
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
    g_xy  = &xy;
    g_xy0 = xy0;
    g_sp  = sp;

    mel_info->fillMacroInfo(mel);


    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
	mel = mel->getNeighbour(k);
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
    inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
    for (i=0; i<=dim; i++) bary[i] = final_lambda[i];
  
    DELETE mel_info;

    return(inside);
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
				Element           **elp, 
				DimVec<double>&     bary,
				const MacroElement *start_mel,
				const WorldVector<double>  *xy0,
				double             *sp)
  {
    ElInfo *el_info = NULL;
    int val;

    el_info = createNewElInfo();

    val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);

    *elp = el_info->getElement();

    DELETE el_info;

    return(val);
  }



  bool Mesh::findElementAtPointRecursive(ElInfo           *el_info,
					 const DimVec<double>& lambda,
					 int                   outside,
					 ElInfo* final_el_info)
  {
    FUNCNAME("Mesh::findElementAtPointRecursive");
    Element *el = el_info->getElement();
    ElInfo *c_el_info = NULL;
    DimVec<double> c_lambda(dim, NO_INIT);
    int i, inside;
    int     ichild, c_outside;

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
	for (i=0; i<=dim; i++)  final_lambda[i] = lambda[i];
	return(true);
      }
      else 
	{  /* outside */
	  if (g_xy0) 
	    { /* find boundary point of [xy0, xy] */
	      double s;
	      el_info->worldToCoord(*(g_xy0), &c_lambda);
	      s = lambda[outside] / (lambda[outside] - c_lambda[outside]);
	      for (i=0; i<=dim; i++) 
		{
		  final_lambda[i] = s * c_lambda[i] + (1.0-s) * lambda[i];
		}
	      if (g_sp) *(g_sp) = s;
	      if(dim == 3) 
		MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s);

	      return(false);  /* ??? */
	    }
	  else return(false);
	}
    }

    c_el_info = createNewElInfo();

    if(dim == 1) {
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (outside >= 0) {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) ERROR("point outside domain\n");
	} else {
	  c_lambda[0] = lambda[0] - lambda[1];
	  c_lambda[1] = 2.0 * lambda[1];
	}
      } else {
	c_el_info->fillElInfo(1, el_info);
	if (outside >= 0)  {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) ERROR("point outside domain\n");
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

    if(dim == 2) {
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) {
	    ERROR("outside curved boundary child 0\n");
	  }
	} else {
	  c_lambda[0] = lambda[2];
	  c_lambda[1] = lambda[0] - lambda[1];
	  c_lambda[2] = 2.0 * lambda[1];
	}
      } else {
	c_el_info->fillElInfo(1, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) {
	    ERROR("outside curved boundary child 1\n");
	  }
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

    if(dim == 3) {
      if (el->isNewCoordSet()) {
	if (lambda[0] >= lambda[1])
	  ichild = 0;
	else
	  ichild = 1;
	c_el_info->fillElInfo(ichild, el_info);
	c_outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);

	if (c_outside>=0) {  /* test is other child is better... */
	  DimVec<double> c_lambda2(dim, NO_INIT);
	  int c_outside2;
	  ElInfo *c_el_info2 = createNewElInfo();

	  c_el_info2->fillElInfo(1-ichild, el_info);
	  c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
	      ichild, c_outside, c_lambda[0],c_lambda[1],c_lambda[2],c_lambda[3]);
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
	      1-ichild, c_outside2, c_lambda2[0],c_lambda2[1],c_lambda2[2],
	      c_lambda2[3]);

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
	    for (i=0; i<=dim; i++) c_lambda[i] = c_lambda2[i];
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
	  DELETE c_el_info2;
	}
	outside = c_outside;
      } else {  /* no new_coord */
	if (lambda[0] >= lambda[1]) {
	  c_el_info->fillElInfo(0, el_info);
	  c_lambda[0] = lambda[0] - lambda[1];
	  c_lambda[1] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][0][1]];
	  c_lambda[2] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][0][2]];
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][1][1]];
	  c_lambda[2] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][1][2]];
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, 
					 final_el_info);
    DELETE c_el_info;

    return(inside); 
  }


  void Mesh::setDiameter(const WorldVector<double>& w) { diam = w; }

  void Mesh::setDiameter(int i, double w) { diam[i] = w; }


  int Mesh::newDOFFct1(ElInfo* ei) {
    ei->getElement()->newDOFFct1(compressAdmin);
    return 0;
  }

  int Mesh::newDOFFct2(ElInfo* ei) {
    ei->getElement()->newDOFFct2(compressAdmin);
    return 0;
  }

  void Mesh::serialize(::std::ostream &out)
  {
    serializedDOFs.clear();

    // write name
    out << name << ::std::endl;

    // write dim
    out.write(reinterpret_cast<const char*>(&dim), sizeof(int));

    // write nVertices
    out.write(reinterpret_cast<const char*>(&nVertices), sizeof(int));

    // write nEdges
    out.write(reinterpret_cast<const char*>(&nEdges), sizeof(int));

    // write nLeaves
    out.write(reinterpret_cast<const char*>(&nLeaves), sizeof(int));

    // write nElements
    out.write(reinterpret_cast<const char*>(&nElements), sizeof(int));

    // write nFaces
    out.write(reinterpret_cast<const char*>(&nFaces), sizeof(int));

    // write maxEdgeNeigh
    out.write(reinterpret_cast<const char*>(&maxEdgeNeigh), sizeof(int));

    // write diam
    diam.serialize(out);

    // write preserveCoarseDOFs
    out.write(reinterpret_cast<const char*>(&preserveCoarseDOFs), sizeof(bool));

    // write nDOFEl
    out.write(reinterpret_cast<const char*>(&nDOFEl), sizeof(int));

    // write nDOF
    nDOF.serialize(out);

    // write nNodeEl
    out.write(reinterpret_cast<const char*>(&nNodeEl), sizeof(int));

    // write node
    node.serialize(out);

    // write admins
    int i, size = static_cast<int>(admin.size());
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
    for (i = 0; i < size; i++) {
      admin[i]->serialize(out);
    }

    // write macroElements
    size = static_cast<int>(macroElements.size());
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
    for (i = 0; i < size; i++) {
      macroElements[i]->serialize(out);
    }

    // write elementIndex
    out.write(reinterpret_cast<const char*>(&elementIndex), sizeof(int));

    // write initialized
    out.write(reinterpret_cast<const char*>(&initialized), sizeof(bool));

    serializedDOFs.clear();
  }

  void Mesh::deserialize(::std::istream &in)
  {
    serializedDOFs.clear();

    // read name
    in >> name;
    in.get();

    // read dim
    int oldVal = dim;
    in.read(reinterpret_cast<char*>(&dim), sizeof(int));
    TEST_EXIT((oldVal == 0) || (dim == oldVal))("invalid dimension\n");

    // read nVertices
    in.read(reinterpret_cast<char*>(&nVertices), sizeof(int));

    // read nEdges
    in.read(reinterpret_cast<char*>(&nEdges), sizeof(int));

    // read nLeaves
    in.read(reinterpret_cast<char*>(&nLeaves), sizeof(int));

    // read nElements
    in.read(reinterpret_cast<char*>(&nElements), sizeof(int));

    // read nFaces
    in.read(reinterpret_cast<char*>(&nFaces), sizeof(int));

    // read maxEdgeNeigh
    in.read(reinterpret_cast<char*>(&maxEdgeNeigh), sizeof(int));

    // diam
    diam.deserialize(in);

    // read preserveCoarseDOFs
    in.read(reinterpret_cast<char*>(&preserveCoarseDOFs), sizeof(bool));

    // read nDOFEl
    oldVal = nDOFEl;
    in.read(reinterpret_cast<char*>(&nDOFEl), sizeof(int));
    TEST_EXIT((oldVal == 0) || (nDOFEl == oldVal))("invalid nDOFEl\n");

    // read nDOF
    nDOF.deserialize(in);

    // read nNodeEl
    oldVal = nNodeEl;
    in.read(reinterpret_cast<char*>(&nNodeEl), sizeof(int));
    TEST_EXIT((oldVal == 0) || (nNodeEl == oldVal))("invalid nNodeEl\n");

    // read node
    node.deserialize(in);

    // read admins
    int i, size;
    in.read(reinterpret_cast<char*>(&size), sizeof(int));
    admin.resize(size, NULL);
    for (i = 0; i < size; i++) {
      if (!admin[i]) {
	admin[i] = NEW DOFAdmin(this);
      }
      admin[i]->deserialize(in);
    }

    // read macroElements
    in.read(reinterpret_cast<char*>(&size), sizeof(int));

    ::std::vector< ::std::vector<int> > neighbourIndices(size);

    for (i = 0; i < static_cast<int>(macroElements.size()); i++) {
      if (macroElements[i]) {
	DELETE macroElements[i];
      }
    }
    macroElements.resize(size);
    for(i = 0; i < size; i++) {
      macroElements[i] = NEW MacroElement(dim);
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
    }

    // read elementIndex
    in.read(reinterpret_cast<char*>(&elementIndex), sizeof(int));

    // read initialized
    in.read(reinterpret_cast<char*>(&initialized), sizeof(bool));

    // set neighbour pointer in macro elements
    int j, neighs = getGeo(NEIGH);
    for(i = 0; i < static_cast<int>(macroElements.size()); i++) {
      for(j = 0; j < neighs; j++) {
	int index = neighbourIndices[i][j];
	if(index != -1) {
	  macroElements[i]->setNeighbour(j, macroElements[index]);
	} else {
	  macroElements[i]->setNeighbour(j, NULL);
	}
      }
    }

    // set mesh pointer in elements
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
    while(elInfo) {
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
  }

  void Mesh::initialize() 
  {
    ::std::string macroFilename("");
    ::std::string valueFilename("");
    ::std::string periodicFile("");
    int check = 1;

    GET_PARAMETER(0, name + "->macro file name",  &macroFilename);
    GET_PARAMETER(0, name + "->value file name",  &valueFilename);
    GET_PARAMETER(0, name + "->periodic file", &periodicFile);
    GET_PARAMETER(0, name + "->check", "%d", &check);
    GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs);

    if (macroFilename.length()) {
      macroFileInfo_ = MacroReader::readMacro(macroFilename.c_str(), 
					      this,
					      periodicFile == "" ? NULL : periodicFile.c_str(),
					      check);

      // If there is no value file which should be written, we can delete
      // the information of the macro file.
      if (!valueFilename.length()) {
	clearMacroFileInfo();
      }
    }

    initialized = true;
  }

  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
    ::std::map<BoundaryType, VertexVector*>::iterator it;
    ::std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
    for (it = periodicAssociations.begin(); it != end; ++it) {
      if ((*(it->second))[dof1] == dof2)
	return true;
    }
    return false;
  }

  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
    ::std::vector<DegreeOfFreedom> associatedToDOF1;
    int i, size;
    ::std::map<BoundaryType, VertexVector*>::iterator it;
    ::std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
    DegreeOfFreedom dof, assDOF;

    associatedToDOF1.push_back(dof1);
    for(it = periodicAssociations.begin(); it != end; ++it) {
      size = static_cast<int>(associatedToDOF1.size());
      for(i = 0; i < size; i++) {
	dof = associatedToDOF1[i];
	assDOF = (*(it->second))[dof];
	if(assDOF == dof2) {
	  return true;
	} else {
	  if(assDOF != dof) associatedToDOF1.push_back(assDOF);
	}
      }
    }
    return false;
  }

  void Mesh::clearMacroFileInfo()
  {
    macroFileInfo_->clear(getNumberOfEdges(),
			  getNumberOfVertices());
    DELETE macroFileInfo_;
  }
}