#include "MacroReader.h"
#include "MacroWriter.h"
#include "MacroElement.h"
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
#include <string.h>
#include "FixVec.h"
#include "FixVecConvert.h"
#include "PeriodicMap.h"
#include "ElInfo.h"
#include "Parameters.h"
#include "DOFIterator.h"
#include "SurfaceRegion_ED.h"
#include "ElementRegion_ED.h"
#include "LeafData.h"
#include "VertexVector.h"
#include <map>
#include <iostream>
#include <fstream>

namespace AMDiS {

  MacroInfo* MacroReader::readMacro(const char *filename, 
				    Mesh* mesh,
				    const char *periodicFile,
				    int check)
  {
    FUNCNAME("Mesh::readMacro()");

    TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
   
    MacroInfo *macroInfo = new MacroInfo();
    macroInfo->readAMDiSMacro(filename, mesh);

    std::deque<MacroElement*>::iterator mel = macroInfo->mel.begin();
    int **melVertex = macroInfo->mel_vertex;
    WorldVector<double> *coords = macroInfo->coords;
    DegreeOfFreedom **dof = macroInfo->dof;

    // === read periodic data =================================
    if (periodicFile && (strcmp(periodicFile, "") != 0)) {
      WARNING("periodic boundaries may lead to errors in small meshes if element neighbours not set\n");
    
      FILE *file = fopen(periodicFile, "r");
      TEST_EXIT(file)("can't open file %s\n", periodicFile);

      int n;
      int dim = mesh->getDim();

      int el1, el2;
      int *verticesEl1 = GET_MEMORY(int, dim);
      int *verticesEl2 = GET_MEMORY(int, dim);
      int mode = -1; // 0: drop dofs, 1: associate dofs
      int result;
      BoundaryType boundaryType;

      fscanf(file, "%*s %d", &n);

      fscanf(file, "%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s");

      PeriodicMap periodicMap;
    
      for (int i = 0; i < n; i++) {
	std::map<int, int> vertexMapEl1;
	std::map<int, int> vertexMapEl2;

	result = fscanf(file, "%d", &mode);
	TEST_EXIT(result == 1)("mode?\n");
      
	result = fscanf(file, "%d", &boundaryType);
	TEST_EXIT(result == 1)("boundaryType?\n");
      
	result = fscanf(file, "%d", &el1);
	TEST_EXIT(result == 1)("el1?\n");

	for (int j = 0; j < dim; j++) {
	  result = fscanf(file, "%d", &verticesEl1[j]);
	  TEST_EXIT(result == 1)("vertEl1[%d]\n", j);
	}
	result = fscanf(file, "%d", &el2);
	TEST_EXIT(result == 1)("el2?\n");
	for (int j = 0; j < dim; j++) {
	  result = fscanf(file, "%d", &verticesEl2[j]);
	  TEST_EXIT(result == 1)("vertEl2[%d]\n", j);
	}
	for (int j = 0; j < dim; j++) {
	  if (mode == 0) {
	    periodicMap.setEntry(melVertex[el1][verticesEl1[j]], 
				 melVertex[el2][verticesEl2[j]]);
	  }
	  vertexMapEl1[verticesEl1[j]] = verticesEl2[j];
	  vertexMapEl2[verticesEl2[j]] = verticesEl1[j];
	}

	// calculate sides of periodic vertices
	int sideEl1 = 0, sideEl2 = 0;
	if (dim == 1) {
	  sideEl1 = verticesEl1[0];
	  sideEl2 = verticesEl2[0];
	} else {
	  for (int j = 0; j < dim + 1; j++) {
	    sideEl1 += j;
	    sideEl2 += j;
	  }
	  for (int j = 0; j < dim; j++) {
	    sideEl1 -= verticesEl1[j];
	    sideEl2 -= verticesEl2[j];
	  }
	}
	
	// create periodic info
	DimVec<WorldVector<double> > periodicCoordsEl1(dim - 1, NO_INIT);
	DimVec<WorldVector<double> > periodicCoordsEl2(dim - 1, NO_INIT);

	Element *element1 = const_cast<Element*>((*(mel + el1))->getElement());
	Element *element2 = const_cast<Element*>((*(mel + el2))->getElement());
      
	// for all vertices of this side
	for (int j = 0; j < dim; j++) {
	  periodicCoordsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] = 
	    coords[melVertex[el2][vertexMapEl1[verticesEl1[j]]]];
	  periodicCoordsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
	    coords[melVertex[el1][vertexMapEl2[verticesEl2[j]]]];
	}
      
	// decorate leaf data
	ElementData *ld1 = element1->getElementData();
	ElementData *ld2 = element2->getElementData();

	LeafDataPeriodic *ldp1 = dynamic_cast<LeafDataPeriodic*>(ld1->getElementData(PERIODIC));
	LeafDataPeriodic *ldp2 = dynamic_cast<LeafDataPeriodic*>(ld2->getElementData(PERIODIC));

	if (!ldp1) {
	  ldp1 = new LeafDataPeriodic(ld1);
	  element1->setElementData(ldp1);
	}

	if (!ldp2) {
	  ldp2 = new LeafDataPeriodic(ld2);
	  element2->setElementData(ldp2);
	}

	ldp1->addPeriodicInfo(mode,
			      boundaryType, 
			      sideEl1, 
			      &periodicCoordsEl1);

	ldp2->addPeriodicInfo(mode,
			      boundaryType, 
			      sideEl2, 
			      &periodicCoordsEl2);

	if (mode != 0) {
	  VertexVector *associated = mesh->periodicAssociations[boundaryType];

	  if (!associated) {
	    associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
	    mesh->periodicAssociations[boundaryType] = associated;
	    VertexVector::Iterator it(associated, ALL_DOFS);
	    for (it.reset2(); !it.end(); ++it) {
	      *it = it.getDOFIndex();
	    }
	  }

	  for (int j = 0; j < dim; j++) {
	    (*associated)[melVertex[el1][verticesEl1[j]]] =
	      melVertex[el2][vertexMapEl1[verticesEl1[j]]];
	    (*associated)[melVertex[el2][verticesEl2[j]]] =
	      melVertex[el1][vertexMapEl2[verticesEl2[j]]];
	  }
	}
      }    

      FREE_MEMORY(verticesEl1, int, dim);
      FREE_MEMORY(verticesEl2, int, dim);

      // change periodic vertex dofs
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
	if (periodicMap.getEntry(i) != -1) {
	  mesh->freeDOF(dof[i], VERTEX);
	  dof[i] = dof[periodicMap.getEntry(i)];

	  std::map<BoundaryType, VertexVector*>::iterator assoc;
	  std::map<BoundaryType, VertexVector*>::iterator assocEnd =
	      mesh->periodicAssociations.end();

	  for (assoc = mesh->periodicAssociations.begin(); 
	       assoc != assocEnd; 
	       ++assoc) {

	    DegreeOfFreedom a = (*(assoc->second))[i];
	    if (a != i) {
	      (*(assoc->second))[i] = i;
	      (*(assoc->second))[a] = periodicMap.getEntry(i);
	    }
	  }

	}
      }

      std::map<BoundaryType, VertexVector*>::iterator assoc;
      std::map<BoundaryType, VertexVector*>::iterator assocEnd =
	  mesh->periodicAssociations.end();
      for (assoc = mesh->periodicAssociations.begin(); 
	   assoc != assocEnd; 
	   ++assoc) {

	for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
	  if (i != (*(assoc->second))[i])
	    MSG("association %d: vertex %d -> vertex %d\n", 
		assoc->first, i, (*(assoc->second))[i]);
	}
      }

      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
	if (periodicMap.getEntry(i) != -1) {
	  MSG("identification : vertex %d is now vertex %d\n", i, periodicMap.getEntry(i));
	}
      }
    }

    // =========================================================

    for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
      for (int k = 0; k < mesh->getGeo(VERTEX); k++) {
	(*(mel + i))->setCoord(k, coords[melVertex[i][k]]);

	const_cast<Element*>((*(mel + i))->getElement())->
	  setDOF(k, dof[melVertex[i][k]]);
      }
    }

    if (!macroInfo->neigh_set) {
      TEST_EXIT(!periodicFile)
	("periodic boundary condition => element neighbours must be set\n");
	computeNeighbours(mesh);
    } else {
	/****************************************************************************/
	/* fill MEL oppVertex values when reading neighbour information form file  */
	/****************************************************************************/

      for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
	for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	  MacroElement *neigh = const_cast<MacroElement*>(mel[i]->getNeighbour(k));

	  if (neigh) {
	    int j = 0;
	    for (; j < mesh->getGeo(NEIGH); j++)
	      if (neigh->getNeighbour(j) == *(mel + i))  
		break;
	
	    TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour of neighbour %d\n", 
					       mel[i]->getIndex(), neigh->getIndex());
	    mel[i]->setOppVertex(k, j);
	  } else {
	    mel[i]->setOppVertex(k, -1);
	  }
	}
      }
    }

    if (!macroInfo->bound_set) {
      macroInfo->dirichletBoundary();
    }
  
    if (mesh->getDim() > 1)
      boundaryDOFs(mesh);

    // initial boundary projections
    //if(dim > 1) {
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      MacroElement *macroEl = *(mel+i);
      Projection *projector = macroEl->getProjection(0);
      if (projector && projector->getType() == VOLUME_PROJECTION) {
	for (int j = 0; j <= dim; j++) {
	  projector->project(macroEl->getCoord(j));
	}
      } else {
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
	  projector = macroEl->getProjection(numFaces + j);
	  if (projector) {
	    int vertex0 = Global::getReferenceElement(dim)->getVertexOfEdge(j, 0);
	    int vertex1 = Global::getReferenceElement(dim)->getVertexOfEdge(j, 1);
	    projector->project(macroEl->getCoord(vertex0));
	    projector->project(macroEl->getCoord(vertex1));
	  }
	}
      }
    }
    //}

    macroInfo->fillBoundaryInfo(mesh);

    if (mesh->getNumberOfDOFs(CENTER)) {
      for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
	const_cast<Element*>(mel[i]->getElement())->
	  setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER));
      }
    }

    /****************************************************************************/
    /* domain size                                                              */
    /****************************************************************************/

    WorldVector<double> x_min, x_max;

    for (int j = 0; j < Global::getGeo(WORLD); j++) {
      x_min[j] =  1.E30;
      x_max[j] = -1.E30;
    }

    for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
      for (int j = 0; j < Global::getGeo(WORLD); j++) {
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
      }
    }

    for (int j = 0; j < Global::getGeo(WORLD); j++)
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

      if (mesh->getDim() > 1) {
	char filenew[128];
	strncpy(filenew, filename, 128); 
	filenew[127] = 0;
	strncat(filenew, ".new", 128);   
	filenew[127] = 0;
	macroTest(mesh, filenew);
      }
    }

    return macroInfo;
  }

  /****************************************************************************/
  /*  fill macro info structure and some pointers in mesh ...                 */
  /****************************************************************************/

  void MacroInfo::fill(Mesh *pmesh, int ne, int nv)
  {
    FUNCNAME("MacroInfo::fill()");

    TEST_EXIT(pmesh)("no mesh\n");

    int dim = pmesh->getDim(); 
    mesh = pmesh;

    mesh->setNumberOfElements(ne);
    mesh->setNumberOfLeaves(ne);
    mesh->setNumberOfVertices(nv);

    for (int i = 0; i < ne; i++) {
      MacroElement *newMacro = new MacroElement(mesh->getDim());
      mel.push_back(newMacro);
      mesh->addMacroElement(mel[i]);
    }

    dof = GET_MEMORY(DegreeOfFreedom*, nv);
    coords = GET_MEMORY(WorldVector<double>, nv);
    mel_vertex = GET_MEMORY(int*, ne);

    for (int i = 0; i < ne; i++) {
      mel_vertex[i] = GET_MEMORY(int, mesh->getGeo(VERTEX));
    }

    for (int i = 0; i < nv; i++) {
      dof[i] = mesh->getDOF(VERTEX);
    }

    for (int i = 0; i < ne; i++) {
      mel[i]->element = mesh->createNewElement();
      (mel)[i]->index = i;

      if (dim == 3) {
	(mel)[i]->elType = 0;
      }
    }
    neigh_set = false;
    bound_set = false;
  }

  void MacroInfo::clear(int ne, int nv)
  {
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
      FREE_MEMORY(mel_vertex[i], int, mesh->getGeo(VERTEX));

    FREE_MEMORY(mel_vertex, int*, ne);
    FREE_MEMORY(coords, WorldVector<double>, nv);
    coords = NULL;  
    FREE_MEMORY(dof, DegreeOfFreedom*, nv);
    dof = NULL;

    mesh = NULL;
    neigh_set = false;
  }

  /****************************************************************************/
  /****************************************************************************/
  /*  tool for reading macro triangulations in ALBERT-format                  */
  /****************************************************************************/
  /****************************************************************************/

  /****************************************************************************/
  /*  read_indices()  reads dim+1 indices from  file  into  id[0-dim],        */
  /*    returns true if dim+1 inputs arguments could be read successfully by  */
  /*    fscanf(), else false                                                  */
  /****************************************************************************/

  int  MacroInfo::read_indices(FILE *file, DimVec<int> &id)
  {
    int dim = mesh->getDim();

    for (int i = 0; i <= dim; i++) {
      if (fscanf(file, "%d", &id[i]) != 1)
	return(false);
    }

    return(true);
  }

#define N_KEYS      14
#define N_MIN_KEYS  7
  static const char *keys[N_KEYS] = {
    "DIM",                   //  0 
    "DIM_OF_WORLD",          //  1
    "number of vertices",    //  2
    "number of elements",    //  3
    "vertex coordinates",    //  4
    "element vertices",      //  5
    "element boundaries",    //  6
    "element neighbours",    //  7
    "element type",          //  8
    "projections",           //  9
    "element region",        // 10
    "surface region",        // 11
    "mesh name",             // 12
    "time"                   // 13
  };

  static int get_key_no(const char *key)
  {
    for (int i = 0; i < N_KEYS; i++)
      if (!strcmp(keys[i], key))  
	return(i);

    return(-1);
  }

#include <ctype.h>

  static const char *read_key(const char *line)
  {
    static char key[100];
    char *k = key;

    while (isspace(*line)) 
	line++;
    while ((*k++ = *line++) != ':');
    *--k = '\0';
  
    return(const_cast<const char *>( key));
  }

  /****************************************************************************/
  /*  read_albert_macro():                                                    */
  /*    read macro triangulation from ascii file in ALBERT format             */
  /*    fills macro_info structure                                            */
  /*    called by read_macro(), fills missing information                     */
  /****************************************************************************/


  void MacroInfo::readAMDiSMacro(const char *filename, Mesh* mesh)
  {
    FUNCNAME("MacroInfo::readAMDiSMacro()");

    FILE *file;
    int dim;
    int dow, nv, ne, j, k;
    double dbl;
    char name[128], line[256];
    int line_no, n_keys, i_key, sort_key[N_KEYS], nv_key, ne_key;
    int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0};
    const char *key;
    DimVec<int> *ind = NULL;

    TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
    TEST_EXIT(strlen(filename) < static_cast<unsigned int>(127))
      ("can only handle filenames up to 127 characters\n");

    file = fopen(filename, "r");
    TEST_EXIT(file)("cannot open file %s\n", filename);
    strncpy(name, filename, 127);

    /****************************************************************************/
    /*  looking for all keys in the macro file ...                              */
    /****************************************************************************/

    line_no = n_keys = 0;
    while (fgets(line, 255, file)) {
	line_no++;
	if (!strchr(line, ':'))  continue;
	key = read_key(line);
	i_key = get_key_no(key);
	TEST_EXIT(i_key >= 0)
	  ("macro file %s must not contain key %s on line %d\n",
	   name, key, line_no);
	TEST_EXIT(!key_def[i_key])
	  ("key %s defined second time on line %d in file %s\n");

	sort_key[n_keys++] = i_key;
	key_def[i_key] = true;
    }

    fclose(file);
    for (i_key = 0; i_key < N_MIN_KEYS; i_key++) {
	for (j = 0; j < n_keys; j++)
	    if (sort_key[j] == i_key)  break;
	TEST_EXIT(j < n_keys)("You do not have specified data for %s in %s\n",
			      keys[i_key], name);

	for (j = 0; j < n_keys; j++)
	  if (sort_key[j] == 2)  break;
	nv_key = j;
	for (j = 0; j < n_keys; j++)
	  if (sort_key[j] == 3)  break;
	ne_key = j;
    
	switch(i_key) {
	    case 0:
	    case 1:
		TEST_EXIT(sort_key[i_key] < 2)
		    ("You have to specify DIM or mesh->getGeo(WORLD) before all other data\n");
		break;
	    case 4: 
		TEST_EXIT(nv_key < i_key)
		    ("Before reading data for %s, you have to specify the %s in file\n",
		     keys[4], keys[2], name);
		break;
	    case 5: 
		TEST_EXIT(nv_key < i_key  &&  ne_key < i_key)
		    ("Before reading data for %s, you have to specify the %s and %s in file %s\n",
		     keys[5], keys[3], keys[2], name);
	    case 6:
	    case 7:
	    case 8:
		TEST_EXIT(ne_key < i_key)
		    ("Before reading data for %s, you have to specify the %s in file %s\n",
		     keys[i_key], keys[3], name);
	}
    }

    for (i_key = 0; i_key < N_KEYS; i_key++)
      key_def[i_key] = false;

    /****************************************************************************/
    /*  and now, reading data ...                                               */
    /****************************************************************************/
	
    file = fopen(name, "r");
    TEST_EXIT(file)("cannot open file %s\n",name);

    int result;

    for (i_key = 0; i_key < n_keys; i_key++) {

	switch(sort_key[i_key]) {
	  case 0:
	    result = fscanf(file, "%*s %d", &dim);
	    TEST_EXIT(result == 1)
	      ("cannot read DIM correctly in file %s\n", name);

	    ind = new DimVec<int>(dim, NO_INIT);

	    key_def[0] = true;
	    break;
	  case 1:
	    result = fscanf(file, "%*s %d", &dow);
	    TEST_EXIT(result == 1)
	      ("cannot read Global::getGeo(WORLD) correctly in file %s\n", name);
	    TEST_EXIT(dow == Global::getGeo(WORLD))
	      ("dimension of world = %d != Global::getGeo(WORLD) = %d\n", 
	       dow, Global::getGeo(WORLD));

	    key_def[1] = true;
	    break;
	  case 2:
	    result = fscanf(file, "%*s %*s %*s %d", &nv);
	    TEST_EXIT(result == 1)
	      ("cannot read number of vertices correctly in file %s\n", name);
	    TEST_EXIT(nv > 0)
	      ("number of vertices = %d must be bigger than 0\n", nv);

	    key_def[2] = true;
	    if (key_def[3])
	      fill(mesh, ne, nv);
	    break;
	  case 3:
	    result = fscanf(file, "%*s %*s %*s %d", &ne);
	    TEST_EXIT(result == 1)
	      ("cannot read number of elements correctly in file %s\n", name);
	    TEST_EXIT(ne > 0)
	      ("number of elements = %d must be bigger than 0\n", ne);

	    key_def[3] = true;
	    if (key_def[2])
	      fill(mesh, ne, nv);
	    break;
	  case 4:
	    fscanf(file, "%*s %*s");
	    for (int i = 0; i < nv; i++) {
	      for (j = 0; j <Global::getGeo(WORLD) ; j++) {
		result = fscanf(file, "%lf", &dbl);
		TEST_EXIT(result == 1)
		  ("error while reading coordinates, check file %s\n", name);
		coords[i][j] = dbl;
	      }
	    }
	    key_def[4] = true;
	    break;
	  case 5:
	    fscanf(file, "%*s %*s");
	    /****************************************************************************/
	    /* global number of vertices on a single element                            */
	    /****************************************************************************/

	    for (int i = 0; i < ne; i++) {
	      result = read_indices(file, *ind);
	      TEST_EXIT(result)
		("cannot read vertex indices of element %d in file %s\n",  i, name);

	      for (k = 0; k < mesh->getGeo(VERTEX); k++)
		mel_vertex[i][k] = (*ind)[k];
	    }

	    key_def[5] = true;
	    break;
	  case 6:
	    fscanf(file, "%*s %*s");
	    /****************************************************************************/
	    /* MEL boundary pointers                                                    */
	    /****************************************************************************/
	    for (int i = 0; i < ne; i++) {  
                // boundary information of ith element 

	      result = read_indices(file, *ind);
	      TEST_EXIT(result)
		("cannot read boundary type of element %d in file %s\n", i, name);

	      // fill boundary of macro-element
	      MacroReader::fillMelBoundary(mesh, 
					   mel[i], 
					   VecConv<int,NEIGH,PARTS>::convertVec((*ind), mesh));
	    }

	    this->fillBoundaryInfo(mesh);
                   
	    bound_set = true;
	    key_def[6] = true;
	    break;
	  case 7:
	    fscanf(file, "%*s %*s");
	    /****************************************************************************/
	    /* fill MEL neighbour pointers:                                             */
	    /*   if they are specified in the file: read them from file,                */
	    /*   else init them by a call of fill_mel_neighbour()                       */
	    /****************************************************************************/
	    neigh_set = true;
	    for (int i = 0; i < ne; i++) {
                //  neighbour information about ith element

		if (read_indices(file, *ind))
		  MacroReader::fillMelNeigh(mel[i], mel, 
					    VecConv<int,NEIGH,PARTS>::convertVec((*ind), 
										 mesh));
		else {
		  neigh_set = false; /* setting of neighbours fails :-( */
		  break;
		}
	      }

	    key_def[7] = true;
	    break;
	  case 8:
	    fscanf(file, "%*s %*s");
	    /****************************************************************************/
	    /* MEL elType                                                               */
	    /****************************************************************************/

	    if (dim == 2 || dim == 1)
	      ERROR("there is no element type in 2d and 2d; ignoring data for elType\n");

	    for (int i = 0; i < ne; i++) {
	      result = fscanf(file, "%d", &j);
	      TEST_EXIT(result == 1)
		("cannot read elType of element %d in file %s\n",
		 i, name);
	      if (dim == 3) {
		(mel)[i]->elType = j;
	      }
	    }

	    key_def[8] = true;
	    break;
	  case 9:
	    {
	      fscanf(file, "%*s");

	      int numFaces = mesh->getGeo(FACE);
	      int numEdgesAtBoundary = 0;

	      for (k = 1; k < dim; k++) {
		numEdgesAtBoundary += k;
	      }

	      for (int i = 0; i < ne; i++) {
		result = read_indices(file, *ind);
		TEST_EXIT(result)
		  ("cannot read boundary projector of element %d in file %s\n", i, name);
	
		Projection *projector = Projection::getProjection((*ind)[0]);

		if(projector && projector->getType() == VOLUME_PROJECTION) {
		  mel[i]->setProjection(0, projector);
		} else { // boundary projection
		  for(j = 0; j < mesh->getGeo(NEIGH); j++) {
		    projector = Projection::getProjection((*ind)[j]);
		    if(projector) {
		      mel[i]->setProjection(j, projector);
		      if(dim > 2) {
			for(k = 0; k < numEdgesAtBoundary; k++) {
			  int edgeNr = Global::getReferenceElement(dim)->getEdgeOfFace(j, k);
			  mel[i]->setProjection(numFaces + edgeNr, projector);
			}
		      }
		    }
		  }
		}
	      }
	    }
	    key_def[9] = true;
	    break;
	  case 10:
	    fscanf(file, "%*s %*s");
	    /****************************************************************************/
	    /* MEL regions                                                              */
	    /****************************************************************************/

	    for (int i = 0; i < ne; i++) {
	      result = fscanf(file, "%d", &j);
	      TEST_EXIT(result == 1)
		("cannot read region of element %d in file %s\n", i, name);
	      if (j >= 0) {
		Element *el = mel[i]->getElement();
		ElementRegion_ED *elementRegion = 
		  new ElementRegion_ED(el->getElementData());
		elementRegion->setRegion(j);
		el->setElementData(elementRegion);
	      }
	    }
	    key_def[10] = true;
	    break;
	  case 11:
	    fscanf(file, "%*s %*s");
	    for (int i = 0; i < ne; i++) {
	      result = read_indices(file, *ind);
	      TEST_EXIT(result)
		("cannot read surface regions of element %d in file %s\n", i, name);

	      Element *el = mel[i]->getElement();

	      for(j = 0; j < mesh->getGeo(NEIGH); j++) {
		if((*ind)[j] >= 0) {
		  SurfaceRegion_ED *surfaceRegion = 
		    new SurfaceRegion_ED(el->getElementData());
		  surfaceRegion->setSide(j);
		  surfaceRegion->setRegion((*ind)[j]);
		  el->setElementData(surfaceRegion);
		}
	      }
	    }
	    key_def[11] = true;
	    break;
	  case 12:
	    fscanf(file, "%*s %*s %*s");
	    break;
	  case 13:
	    fscanf(file, "%*s %*s");
	    break;
	  }
      }

    if (ind) {
      delete ind;
    }

    fclose(file);
  }


  int macro_type(const char *filename, const char *type)
  {
    const char *fn, *t;
  
    if (strlen(filename) <= strlen(type))
      return(false);
  
    fn = filename;
    while (*fn) fn++;
    t = type;
    while (*t) t++;
  
    while (t != type  &&  *t == *fn) t--;
  
    return(t == type);
  }


  /****************************************************************************/
  /*  sets the boundary of all edges/faces with no neigbour to a straight     */
  /*  line/face with Dirichlet boundary type                                  */
  /****************************************************************************/

  void MacroInfo::dirichletBoundary()
  {
    for (int i = 0; i < static_cast<int>( mel.size()); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	if (mel[i]->neighbour[k])
	  mel[i]->boundary[k] = INTERIOR;
	else
	  mel[i]->boundary[k] = DIRICHLET;
      }
    }
  }


  void MacroInfo::fillBoundaryInfo(Mesh *mesh)
  {
    int i,j,k, nv = mesh->getNumberOfVertices();

    std::deque<MacroElement*>::iterator melIt;

    BoundaryType *bound = GET_MEMORY(BoundaryType, nv);

    int dim = mesh->getDim();

    switch(dim) {
    case 1:
      break;
    case 2:
      for (i = 0; i < nv; i++)
	bound[i] = INTERIOR;

      for (i=0, melIt = mesh->firstMacroElement(); 
	   melIt != mesh->endOfMacroElements(); 
	   ++melIt, ++i) 
	{
	  for (j = 0; j < mesh->getGeo(NEIGH); j++) {
	    if ((*melIt)->getBoundary(j) != INTERIOR) {
	      if ((*melIt)->getBoundary(j) >= DIRICHLET) {
		int j1 = mel_vertex[i][(j+1)%3];
		int j2 = mel_vertex[i][(j+2)%3];

		bound[j1] = 
		  max(bound[j1], (*melIt)->getBoundary(j));
		bound[j2] = 
		  max(bound[j2], (*melIt)->getBoundary(j));
	      } 
	      else if ((*melIt)->getBoundary(j) <= NEUMANN) {
		int j1 = mel_vertex[i][(j+1)%3];
		int j2 = mel_vertex[i][(j+2)%3];

		if (bound[j1] != INTERIOR)
		  bound[j1] = 
		    max(bound[j1], (*melIt)->getBoundary(j));
		else
		  bound[j1] = (*melIt)->getBoundary(j);

		if (bound[j2] != INTERIOR)
		  bound[j2] = 
		    max(bound[j2], (*melIt)->getBoundary(j));
		else
		  bound[j2] = (*melIt)->getBoundary(j);
	      }
	    }
	  }
	}

      for (i=0, melIt = mesh->firstMacroElement(); 
	   melIt != mesh->endOfMacroElements(); 
	   ++melIt, i++) 
	{
	  for (j = 0; j < mesh->getGeo(VERTEX); j++)
	    (*melIt)->setBoundary(3 + j, bound[mel_vertex[i][j]]);
	}
      break;
    case 3:
      for (i = 0; i < nv; i++)
	bound[i] = INTERIOR;

      for (i=0, melIt = mesh->firstMacroElement(); 
	   melIt != mesh->endOfMacroElements(); 
	   ++melIt, i++) 
	{
	  for (j = 0; j < mesh->getGeo(NEIGH); j++) {
	    for (k = 1; k < 4; k++)
	      bound[mel_vertex[i][(j+k)%4]] =
		((*melIt)->getBoundary(j) != INTERIOR) ?
		newBound((*melIt)->getBoundary(j),
			 bound[mel_vertex[i][(j+k)%4]]) :
		//(*melIt)->getBoundary(j)->
		//newVal(bound[data->mel_vertex[i][(j+k)%4]]) :
		bound[mel_vertex[i][(j+k)%4]];
	  }
	}

      for (i = 0, melIt = mesh->firstMacroElement(); 
	   melIt != mesh->endOfMacroElements(); 
	   ++melIt, i++) 
	{
	  for (j = 0; j < mesh->getGeo(VERTEX); j++)
	    (*melIt)->setBoundary(10 + j, bound[mel_vertex[i][j]]);
	}
      break;
    default: ERROR_EXIT("invalid dim\n");
    }

    FREE_MEMORY(bound, BoundaryType, nv);
  }

  void MacroReader::computeNeighbours(Mesh *mesh)
  {
    FUNCNAME("MacroReader::computeNeighbours()");

    int dim = mesh->getDim();
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);

    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
	mesh->getMacroElement(i)->setNeighbour(k, NULL);
      }
    }

    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	if (mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
	  mesh->getMacroElement(i)->setNeighbour(k, NULL);
	  mesh->getMacroElement(i)->setOppVertex(k, -1);
	  continue;
	}

	if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
	  if (dim == 1) {
	    dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						  getElement()->getDOF(k));
	  } else {
	    for (int l = 0; l < dim; l++)
	      dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						    getElement()->
						    getDOF((k + l + 1) % (dim + 1)));
	  }
	  
	  int j = 0;
	  for (j = i + 1; j < mesh->getNumberOfLeaves(); j++) {
	    int m = mesh->getMacroElement(j)->getElement()->oppVertex(dof);
	    if (m != -1) {
	      mesh->getMacroElement(i)->setNeighbour(k, mesh->getMacroElement(j));
	      mesh->getMacroElement(j)->setNeighbour(m, mesh->getMacroElement(i));
	      mesh->getMacroElement(i)->setOppVertex(k, m);
	      mesh->getMacroElement(j)->setOppVertex(m, k);
	      break;
	    }
	  }

	  if (j >= mesh->getNumberOfLeaves()) {
	    std::cout << "----------- ERROR ------------" << std::endl;
	    std::cout << "Cannot find neighbour " << k << " of element " << i << std::endl;
	    std::cout << "  dim = " << dim << std::endl;
	    std::cout << "  coords of element = ";
	    for (int l = 0; l <= dim; l++) {
	      std::cout << mesh->getMacroElement(i)->getCoord(l);
	      if (l < dim) {
		std::cout << " / ";
	      }
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
	    for (int l = 0; l < dim; l++) {
	      std::cout << *(dof[l]) << " ";
	    }
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
	  }    
	}
      }
    }
  }


  /****************************************************************************/
  /*  boundaryDOFs:                                                           */
  /*  adds dof's at the edges of a given macro triangulation and calculates   */
  /*  the number of edges                                                     */
  /****************************************************************************/

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
    std::deque<MacroElement*>::iterator mel;
    const MacroElement* neigh;
    DegreeOfFreedom *dof;

    mesh->setNumberOfEdges(0);
    mesh->setNumberOfFaces(0);

    int dim = mesh->getDim();

    switch(dim) {
    case 2:
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);

	if (ed) {
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
	  for (it = periodicInfos.begin(); it != end; ++it) {
	    if (it->type != 0) {
	      periodic[it->elementSide] = true;
	    }
	  }
	}

	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	  if (!(*mel)->getNeighbour(i) || 
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

	    if (mesh->getNumberOfDOFs(EDGE)) {
	      dof = el->setDOF(lnode + i, mesh->getDOF(EDGE));
      
	      if ((*mel)->getNeighbour(i)) {
		Element *neigh = const_cast<Element*>((*mel)->getNeighbour(i)->getElement());

		if (periodic[i]) {
		  neigh->setDOF(lnode + (*mel)->getOppVertex(i), mesh->getDOF(EDGE));
		} else {
		  neigh->setDOF(lnode + (*mel)->getOppVertex(i), dof);
		}
	      }
	    }
	  }  
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
      for (int i = 0; i < lne; i++) {

	// check for periodic boundary
	Element *el = const_cast<Element*>((*(mel+i))->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
      
	if(ed) {
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
	  for(it = periodicInfos.begin(); it != end; ++it) {
	    if(it->type != 0) {
	      periodic[it->elementSide] = true;
	    }
	  }
	}

	for (k = 0; k < mesh->getGeo(EDGE); k++) {      
	  /*********************************************************************/
	  /* check for not counted edges                                       */
	  /*********************************************************************/
	  n_neigh = 1;

	  if (newEdge(mesh, (*(mel+i)), k, &n_neigh/*, periodicEdge*/)) {
	    mesh->incrementNumberOfEdges(1);
	    max_n_neigh = max(max_n_neigh, n_neigh);
	  }
	}
      
	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
	  neigh = (*(mel+i))->getNeighbour(k);
	  /*********************************************************************/
	  /* face is counted and dof is added by the element with bigger index */
	  /*********************************************************************/
	  if (neigh  &&  (neigh->getIndex() > (*(mel+i))->getIndex()))  continue;
	
	  mesh->incrementNumberOfFaces(1);
	
	  if (mesh->getNumberOfDOFs(FACE)) {
	    TEST_EXIT(!(*(mel+i))->getElement()->getDOF(lnode+k))
	      ("dof %d on element %d already set\n", 
	       lnode+k, (*(mel+i))->getIndex());
	  
	    const_cast<Element*>((*(mel+i))->getElement())->setDOF(lnode+k, 
								   mesh->getDOF(FACE));

	    if (neigh) {
	      ov = (*(mel+i))->getOppVertex(k);
	      TEST_EXIT(!neigh->getElement()->getDOF(lnode+ov))
		("dof %d on neighbour %d already set\n", 
		 lnode+ov, neigh->getIndex());
	    
	      Element *neighEl = 
		const_cast<Element*>((*(mel+i))->getNeighbour(k)->getElement());

	      if (periodic[k]) {
		neighEl->setDOF(lnode+ov, mesh->getDOF(FACE));
	      } else {
		neighEl->setDOF(lnode+ov, const_cast<int*>((*(mel+i))->getElement()->
							   getDOF(lnode+k)));
	      }
	    }
	  }
	}
      }
      break;
    default: ERROR_EXIT("invalid dim\n");
    }
    
    if (3 == dim) {
      mesh->setMaxEdgeNeigh(std::max(8, 2*max_n_neigh));
    } else {
      mesh->setMaxEdgeNeigh(dim-1);    
    }
  }

  /* 
     testet mesh auf auftretende Zyklen
  
     wenn Zyklus auftritt:
     ordnet Eintraege in MacroElement-Struktur um, so dass kein Zyklus auftritt
     erzeugt neue Macro-Datei nameneu mit umgeordnetem Netz 
     (wenn nameneu=NULL wird keine MAcro-Datei erzeugt)
  */      

  void MacroReader::macroTest(Mesh *mesh, const char *nameneu)
  {
    FUNCNAME("MacroReader::macroTest()");
   
    int i = macrotest(mesh);

    if (i >= 0) {
      ERROR("There is a cycle beginning in macro element %d\n", i);
      ERROR("Entries in MacroElement structures get reordered\n");
      umb(NULL, mesh, umbVkantMacro);

      if (nameneu) {
	ERROR_EXIT("mesh->feSpace\n");
      }
    }
  }
  
  /****************************************************************************/
  /*  macro_test():                              Author: Thomas Kastl (1998)  */
  /****************************************************************************/
  /*
    testet mesh auf auftretende Zyklen
  
    wenn mesh zyklenfrei -> -1
    sonst ->  globaler Index des Macroelementes bei dem ein Zyklus beginnt 
  */

  int MacroReader::macrotest(Mesh *mesh)
  {
    FUNCNAME("MacroReader::macrotest()");

    int *test;
    int *zykl;
    std::deque<MacroElement*>::const_iterator macro,mac;
    int flg;
    std::deque<MacroElement*>::const_iterator macrolfd;
    int zykstart;
    int dim = mesh->getDim();

    test = GET_MEMORY(int, mesh->getNumberOfMacros());
    zykl = GET_MEMORY(int, mesh->getNumberOfMacros());
 
    for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
      test[i] = 0;
    }

    zykstart = -1;
    macrolfd = mesh->firstMacroElement();

    while (macrolfd != mesh->endOfMacroElements()) {
      if (test[(*macrolfd)->getIndex()] == 1) {
	macrolfd++;
      } else {
	for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
	  zykl[i] = 0;
	}
    
	macro = macrolfd;
	flg = 2;
	do {
	  if (zykl[(*macro)->getIndex()] == 1) {
	    flg = 0;
	    zykstart = (*macro)->getIndex();
	  } else {
	    zykl[(*macro)->getIndex()] = 1;
      
	    if (test[(*macro)->getIndex()] == 1) {
	      flg = 1;
	    } else if ((*macro)->getNeighbour(dim) == NULL) {
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
	    }
	    else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim)) {
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
	      test[(*macro)->getNeighbour(dim)->getIndex()] = 1;
	    } else {
	      for (mac = mesh->firstMacroElement();
		   (*mac)!=(*macro)->getNeighbour(dim);
		   mac++);
	      macro = mac;
	    } 
	  }	  
	} while(flg == 2);
 
	if (flg == 1) {
	  macrolfd++;
	} else  { 
	  macrolfd=mesh->endOfMacroElements();
	}
      }
    }
  
    FREE_MEMORY(zykl, int, mesh->getNumberOfMacros());
    FREE_MEMORY(test, int, mesh->getNumberOfMacros());
 
    return zykstart;
  }

  //   waehlt geeignete Verfeinerungskanten, so dass kein Zyklus auftritt (recumb)

  //   ele     Integer-Vektor der Dimension Anzahl der Macro-Elemente
  //           zur Speicherung der neuen Verfeinerungskanten
  //           (wird nur benoetigt, wenn umbvk=umb_vkant_macrodat) 
  
  //   umbvk   Fkt. zur Neuordnung der Verfeinerungskanten
  //           = umb_vkant_macro :
  //               Eintraege in MacroElement-Struktur und Eintraege in macro->el
  //               werden tatsaechlich umgeordnet
  //               -> ALBERT-Routine write_macro kann zum erzeugen einer
  //                  neuen Macro-Datei angewendet werden 
  //           = umb_vkant_macrodat :
  //               Eintraege werden nicht veraendert, es werden nur die lokalen
  //               Indices der Kanten, die zu Verfeinerungskanten werden im
  //               Integer-Vektor ele abgespeichert
  //               -> print_Macrodat zur Erzeugung einer zyklenfreien neuen
  //                  Macro-Datei kann angewendet werden

  void MacroReader::umb(int *ele, Mesh *mesh,
			void (*umbvk)(Mesh*,MacroElement*,int,int*))
  {
    FUNCNAME("MacroReader::umb");

    int *test;
    int i;
  
    test=GET_MEMORY(int, mesh->getNumberOfMacros());
  
    for (i=0; i < static_cast<int>(mesh->getNumberOfMacros()); i++)
      test[i]=0;

    recumb(mesh, (*mesh->firstMacroElement()), NULL,test,0,0,ele,umbvk);

    FREE_MEMORY(test, int, mesh->getNumberOfMacros());
  }

  bool MacroReader::newEdge(Mesh *mesh, MacroElement *mel,
			    int mel_edge_no, int *n_neigh)
  {
    FUNCNAME("MacroElement::newEdge"); 
    MacroElement    *nei;
    const DegreeOfFreedom *dof[2];
    DegreeOfFreedom *edge_dof = NULL;
    int             j, k, opp_v, mel_index, node=0;
    BoundaryType    lbound = INTERIOR;
    Projection *lproject = NULL;
    const int       max_no_list_el = 100;
    BoundaryType *list_bound[100];
    Projection **list_project[100];
    Element *el = const_cast<Element*>(mel->getElement());
    int edge_no = mel_edge_no;

    static int  next_el[6][2] = {{3,2},{1,3},{1,2},{0,3},{0,2},{0,1}};

    int vertices = mesh->getGeo(VERTEX);

    mel_index = mel->getIndex();

    list_bound[0] = &(mel->boundary[mesh->getGeo(FACE)+edge_no]);
    list_project[0] = &(mel->projection[mesh->getGeo(FACE)+edge_no]);

    if (mesh->getNumberOfDOFs(EDGE)) {
      node = mesh->getNode(EDGE);
      if (el->getDOF(node+edge_no)) {
	/****************************************************************************/
	/*  edge was counted by another macro element and dof was added on the      */
	/*  complete patch                                                          */
	/****************************************************************************/
	return false;
      } else {
	edge_dof = el->setDOF(node+edge_no,mesh->getDOF(EDGE));
      }
    }

    for (j = 0; j < 2; j++) {
      dof[j] = el->getDOF(el->getVertexOfEdge(edge_no, j));
    }


    /****************************************************************************/
    /*  first look for all neighbours in one direction until a boundary is      */
    /*  reached :-( or we are back to mel :-)                                   */
    /*  if the index of a neighbour is smaller than the element index, the edge */
    /*  is counted by this neighbour, return 0.                                 */
    /*  If we are back to element, return 1, to count the edge                  */
    /****************************************************************************/

    nei = mel->getNeighbour(next_el[edge_no][0]);
    opp_v = mel->getOppVertex(next_el[edge_no][0]);


    if(mel->getBoundary(next_el[edge_no][0])) {
      lbound = newBound(mel->getBoundary(next_el[edge_no][0]), lbound);
      lproject = mel->getProjection(next_el[edge_no][0]);
    }

    while (nei  &&  nei != mel) {
      for (j = 0; j < vertices; j++)
	if (nei->getElement()->getDOF(j) == dof[0])  break;
      for (k = 0; k < vertices; k++)
	if (nei->getElement()->getDOF(k) == dof[1])  break;

      // check for periodic boundary
      if(j == 4 || k == 4) {
	nei = NULL;
	break;
      }

      if (mesh->getNumberOfDOFs(EDGE))
	TEST_EXIT(nei->index > mel_index)
	  ("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);

      if (!mesh->getNumberOfDOFs(EDGE) &&  nei->getIndex() < mel_index)  return false;

    
      edge_no = Tetrahedron::edgeOfDOFs[j][k];

      TEST_EXIT(*n_neigh < max_no_list_el)
	("too many neigbours for local list\n");

      list_bound[(*n_neigh)] = 
	&(nei->boundary[mesh->getGeo(FACE)+edge_no]);

      list_project[(*n_neigh)++] = 
	&(nei->projection[mesh->getGeo(FACE)+edge_no]);

      if (mesh->getNumberOfDOFs(EDGE))
	{
	  //       if(periodic) {
	  // 	nei->element->setDOF(node+edge_no, mesh->getDOF(EDGE));
	  //       } else {
	  nei->element->setDOF(node+edge_no,edge_dof);
	  //       }
	}

      if (next_el[edge_no][0] != opp_v)
	{
	  if(nei->getBoundary(next_el[edge_no][0])) {
	    lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
	    Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
	    if(!lproject)
	      lproject = neiProject;
	    else {
	      if(neiProject && (lproject->getID() < neiProject->getID())) {
		lproject = neiProject;
	      }
	    }
	  }
	  opp_v = nei->getOppVertex(next_el[edge_no][0]);
	  nei = nei->getNeighbour(next_el[edge_no][0]);
	}
      else
	{
	  if(nei->getBoundary(next_el[edge_no][1])) {
	    lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
	    Projection *neiProject = nei->getProjection(next_el[edge_no][1]);	
	    if(!lproject)
	      lproject = neiProject;
	    else {
	      if(neiProject && (lproject->getID() < neiProject->getID())) {
		lproject = neiProject;
	      }
	    }
	  }
	  opp_v = nei->getOppVertex(next_el[edge_no][1]);
	  nei = nei->getNeighbour(next_el[edge_no][1]);
	}
    }
    if (!nei)
      {
	/****************************************************************************/
	/*  while looping around the edge the domain's boundary was reached. Now,   */
	/*  loop in the other direction to the domain's boundary		    */
	/****************************************************************************/
	edge_no = mel_edge_no;

	nei = mel->getNeighbour(next_el[edge_no][1]);
	opp_v = mel->getOppVertex(next_el[edge_no][1]);
	if(mel->getBoundary(next_el[edge_no][1])) {
	  lbound = newBound(mel->getBoundary(next_el[edge_no][1]), lbound); 
	  Projection *neiProject =  mel->getProjection(next_el[edge_no][1]);
	  if(!lproject)
	    lproject = neiProject;
	  else {
	    if(neiProject && (lproject->getID() < neiProject->getID())) {
	      lproject = neiProject;
	    }
	  }
	}
	while (nei)
	  {
	    for (j = 0; j < vertices; j++)
	      if (nei->getElement()->getDOF(j) == dof[0])  break;
	    for (k = 0; k < vertices; k++)
	      if (nei->getElement()->getDOF(k) == dof[1])  break;

	    // check for periodic boundary
	    if(j == 4 || k == 4) {
	      return false;
	    }

	    if (mesh->getNumberOfDOFs(EDGE))
	      TEST_EXIT(nei->getIndex() > mel_index)
		("neighbour index %d < element index %d\n", nei->getIndex(),
		 mel_index);

	    if (nei->getIndex() < mel_index)  
	      return false;

	    edge_no = Tetrahedron::edgeOfDOFs[j][k];

	    TEST_EXIT(*n_neigh < max_no_list_el)
	      ("too many neigbours for local list\n");

	    list_bound[(*n_neigh)] = 
	      &(nei->boundary[mesh->getGeo(FACE)+edge_no]);

	    list_project[(*n_neigh)++] = 
	      &(nei->projection[mesh->getGeo(FACE)+edge_no]);

	    if (mesh->getNumberOfDOFs(EDGE)) {
	      TEST_EXIT(!nei->getElement()->getDOF(node+edge_no))
		("dof %d on element %d is already set, but not on element %d\n",
		 node + edge_no, nei->getIndex(), mel_index);
	
	      // 	if(periodic) {
	      // 	  nei->element->setDOF(node+edge_no, mesh->getDOF(EDGE));
	      // 	} else {
	      nei->element->setDOF(node+edge_no,edge_dof);
	      // 	}
	    }

	    if (next_el[edge_no][0] != opp_v)
	      {
		if(nei->getBoundary(next_el[edge_no][0])) {
		  lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
		  Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
		  if(!lproject)
		    lproject = neiProject;
		  else {
		    if(neiProject &&( lproject->getID() < neiProject->getID())) {
		      lproject = neiProject;
		    }
		  }
		}

		opp_v = nei->getOppVertex(next_el[edge_no][0]);
		nei = nei->getNeighbour(next_el[edge_no][0]);
	      }
	    else {
	      if(nei->getBoundary(next_el[edge_no][1])) {
		lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound); 
		Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
		if(!lproject)
		  lproject = neiProject;
		else {
		  if(neiProject && (lproject->getID() < neiProject->getID())) {
		    lproject = neiProject;
		  }
		}
	      }

	      opp_v = nei->getOppVertex(next_el[edge_no][1]);
	      nei = nei->getNeighbour(next_el[edge_no][1]);
	    }
	  }
      }
  
    for (j = 0; j < *n_neigh; j++) {
      *(list_bound[j]) = lbound;
      *(list_project[j]) = lproject;
    }
  
    return true;
  }

  void MacroReader::fillMelBoundary(Mesh *mesh, MacroElement *mel, 
				    FixVec<BoundaryType ,NEIGH> ind)
  {
    for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
      mel->boundary[i] = ind[i];    
    }
  }


  void MacroReader::fillMelNeigh(MacroElement *mel,
				 std::deque<MacroElement*>& macro_elements, 
				 FixVec<int,NEIGH> ind)
  {
    int dim = mel->element->getMesh()->getDim();

    for (int k = 0; k < Global::getGeo(NEIGH, dim); k++) {
      if (ind[k] >= 0) 
	mel->neighbour[k] = macro_elements[ind[k]];
      else
	mel->neighbour[k] = NULL;
    }
  }


  //   ordnet Eintraege in Macro-Element macro bzgl. Verfeinerungskante ka um
  //   (coord, bound, boundary, neigh, oppVertex)

  //   ordnet Eintraege in macro->el bzgl. Verfeinerungskante ka um
  //   (Element-DOF's (macro->el->dof) in Ecken und auf Kanten,
  //    wenn NEIGH_IN_EL macro->el->neigh, macro->el->oppVertex)
  //   (wird fuer ALBERT-Routine write_macro benoetigt)

  //   ele wird nicht benoetigt (es kann NULL uebergeben werden)    


  void MacroReader::umbVkantMacro(Mesh *mesh, MacroElement* me, int ka, int *)
  {
    MacroElement* macr=new MacroElement(mesh->getDim());
    int i;
    int n0;
    DegreeOfFreedom *d[7];
  
    int vertices = mesh->getGeo(VERTEX);
    int facesPlusEdges = mesh->getGeo(EDGE) + mesh->getGeo(FACE);

    if (ka == 2);
    else { 
      for (i=0; i < 3; i++) {
	macr->coord[i]=me->coord[i];
	macr->setBoundary(facesPlusEdges + i, me->getBoundary(facesPlusEdges + i));
	macr->setBoundary(i, me->getBoundary(i));
	macr->setNeighbour(i, me->getNeighbour(i));
	macr->setOppVertex(i,me->getOppVertex(i));
      }    
  
      for (i=0; i < 7; i++) {
	d[i]=const_cast<int*>(me->getElement()->getDOF(i));
      }

      if (ka == 1) { 
	me->coord[0] = macr->coord[2];
	me->coord[1] = macr->coord[0];
	me->coord[2] = macr->coord[1];

	me->setBoundary(facesPlusEdges + 0,macr->getBoundary(facesPlusEdges + 2));
	me->setBoundary(facesPlusEdges + 1,macr->getBoundary(facesPlusEdges + 0));
	me->setBoundary(facesPlusEdges + 2,macr->getBoundary(facesPlusEdges + 1));

	me->setBoundary(0, macr->getBoundary(2));
	me->setBoundary(1, macr->getBoundary(0));
	me->setBoundary(2, macr->getBoundary(1));

	me->setNeighbour(0,const_cast<MacroElement*>(macr->getNeighbour(2)));
	me->setNeighbour(1,const_cast<MacroElement*>(macr->getNeighbour(0)));
	me->setNeighbour(2,const_cast<MacroElement*>(macr->getNeighbour(1)));

	me->setOppVertex(0,macr->getOppVertex(2));
	me->setOppVertex(1,macr->getOppVertex(0));
	me->setOppVertex(2,macr->getOppVertex(1));


	if (mesh->getNumberOfDOFs(VERTEX))                /* Ecken */
	  {
	    n0=mesh->getNode(VERTEX);              
        
	    const_cast<Element*>(me->getElement())->setDOF(n0,d[n0+2]);     
	    const_cast<Element*>(me->getElement())->setDOF(n0+1,d[n0]);  
	    const_cast<Element*>(me->getElement())->setDOF(n0+2,d[n0+1]);   
	  }
 
	if (mesh->getNumberOfDOFs(EDGE))                  /* Kanten */
	  { 
	    n0=mesh->getNode(EDGE);    
       
	    const_cast<Element*>(me->getElement())->setDOF(n0,d[n0+2]);  
	    const_cast<Element*>(me->getElement())->setDOF(n0+1,d[n0]);  
	    const_cast<Element*>(me->getElement())->setDOF(n0+2,d[n0+1]);
	  } 

      } else {
	me->coord[0] = macr->coord[1];
	me->coord[1] = macr->coord[2];
	me->coord[2] = macr->coord[0];

	me->setBoundary(facesPlusEdges + 0,macr->getBoundary(facesPlusEdges + 1));
	me->setBoundary(facesPlusEdges + 1,macr->getBoundary(facesPlusEdges + 2));
	me->setBoundary(facesPlusEdges + 2,macr->getBoundary(facesPlusEdges + 0));

	me->setBoundary(0, macr->getBoundary(1));
	me->setBoundary(1, macr->getBoundary(2));
	me->setBoundary(2, macr->getBoundary(0));

	me->setNeighbour(0,const_cast<MacroElement*>(macr->getNeighbour(1)));
	me->setNeighbour(1,const_cast<MacroElement*>(macr->getNeighbour(2)));
	me->setNeighbour(2,const_cast<MacroElement*>(macr->getNeighbour(0)));

	me->setOppVertex(0,macr->getOppVertex(1));
	me->setOppVertex(1,macr->getOppVertex(2));
	me->setOppVertex(2,macr->getOppVertex(0));
    
	if (mesh->getNumberOfDOFs(VERTEX))                /* Ecken */
	  {
	    n0=mesh->getNode(VERTEX);              
        
	    const_cast<Element*>(me->getElement())->setDOF(n0,d[n0+1]);     
	    const_cast<Element*>(me->getElement())->setDOF(n0+1,d[n0+2]);  
	    const_cast<Element*>(me->getElement())->setDOF(n0+2,d[n0]);   
	  }
 
	if (mesh->getNumberOfDOFs(EDGE))                  /* Kanten */
	  { 
	    n0=mesh->getNode(EDGE);    
       
	    const_cast<Element*>(me->getElement())->setDOF(n0,d[n0+1]);  
	    const_cast<Element*>(me->getElement())->setDOF(n0+1,d[n0+2]);  
	    const_cast<Element*>(me->getElement())->setDOF(n0+2,d[n0]);
	  } 
      }
  

      for (i=0; i < vertices; i++)   /* oppVertex der Nachbarn umsetzen*/  
	{
	  if (me->getNeighbour(i))
	    {
	      const_cast<MacroElement*>(me->getNeighbour(i))->setOppVertex(me->getOppVertex(i),i);
	    }
	}
    }
    delete macr;
  }


  //   durchlaeuft rek. Macro-Triangulierung mit Start auf Macro-Element macro und
  //   waehlt geignete Verfeinerungskanten, so dass kein Zyklus auftritt
  
  //   (Umbennenung der lokalen Indices der Macro-Elemente
  //    oder Speichern der neuen Verfeinerungskanten mit Fkt. umbvk)
  
  //   waehlt als neue Verfeinerungskante die laengste Kante eines Elementes
  //   fuehrt "fiktiv verlaengerte" Kanten ein, 
  //   wenn ein Element 2 laengste Kanten besitzt  
 
  //   macroalt  Zeiger auf "Rekursions-Vorgaenger"-Macroelement
  //   lg        Laenge der Verfeinerungskante von "Vorgaenger"
  //   ka = 1, wenn Verfeinerungskante von "Vorgaenger" fiktiv verlaengert
  //   ka = 0, sonst
  //   test      Vektor von Flags, die angeben, ob Macro-Element schon getestet (=1)
  //   ele       Integer-Vektor der Dimension Anzahl der Macro-Elemente
  //             zur Speicherung der neuen Verfeinerungskanten
  //             (wird nur benoetigt, wenn umbvk=umb_vkant_macrodat) 

  void MacroReader::recumb(Mesh *mesh,  
			   MacroElement *mel, MacroElement *macroalt,
			   int *test, double lg, int ka, int *ele, 
			   void (*umbvk)(Mesh *mesh, MacroElement*, int k, int *el))
  {
    int i;
    double l[3];
    int v[3];
    int k = 0;
    MacroElement *n;

    int vertices = mesh->getGeo(VERTEX);

    if (!mel || test[mel->getIndex()]==1) {
      return;
    } else {
      test[mel->getIndex()]=1; 
      
      laengstekante(mel->coord,l,v);
      
      if (v[1] == mesh->getGeo(VERTEX))              /*nur eine laengste Kante*/
	{
	  umbvk(mesh,mel,v[0],ele);      
	  
	  for (i=0; i < vertices; i++) {
	    recumb(mesh,
		   mel->neighbour[i],
		   mel,
		   test,
		   l[0],
		   0,
		   ele,
		   umbvk);
	  }
	  return;
	}
      else
	{ 
	  if (ka == 1)
	    { 
	      if (fabs((l[0]-lg)/lg) < 0.0000001)
		{  
		  for (i=0; i < vertices; i++)
		    {
		      if (mel->neighbour[i] == macroalt)
			{
			  k=i;
			}
		    }
		  
		  umbvk(mesh,mel,k,ele);
		  
		  for (i=0; i < vertices; i++) {
		    recumb(mesh,
			   mel->neighbour[i],
			   mel,
			   test,
			   l[0],
			   0,
			   ele,
			   umbvk);
		  }
		  return;
		}
	    }
	  
	  n=const_cast<MacroElement*>(mel->getNeighbour(v[0])); 
	  /*Nachbar an fiktiv verlaengerter Kante*/
	  umbvk(mesh,mel,v[0],ele);   
	  
	  recumb(mesh, n, mel,test,l[0],1,ele,umbvk);
	  for (i=0; i < vertices; i++) { 
	    recumb(mesh,
		   mel->neighbour[i],
		   mel,
		   test,
		   l[0],
		   0,
		   ele,
		   umbvk);
	  }
	  return;
	}
    }             
  }
  
  //   berechnet aus Koord. der Eckpunkte eines Elementes 
  //   die Laenge der laengsten Kante

  //   l[0] = Laenge der laengsten Kante
  //   v[0] = lokaler Index der laengsten Kante
  //   v[1] = anderer lokaler Index, wenn ex. 2 laengste Kanten
  //        = 3, wenn ex. nur eine laengste Kante
  //   v[2] = dritter lokaler Index, wenn ex. 3 laengste Kanten
  //        = 3, wenn ex. nur eine oder zwei laengste Kanten      

  void MacroReader::laengstekante(FixVec<WorldVector<double>,VERTEX> coord, double *l, int *v)
  {
    int dim = coord.getSize() - 1;

    int i;
    int k;
    double lg;
    double eps;
    double kz; 

    int vertices = Global::getGeo(VERTEX,dim);

    l[0]=absteukl(coord[1],coord[2]);
    l[1]=absteukl(coord[0],coord[2]);
    l[2]=absteukl(coord[0],coord[1]);

    lg=l[0];
    kz=l[0];
    for (i=1; i < vertices; i++)
      {
	if (l[i] > lg) lg=l[i];
	if (l[i] < kz) kz=l[i];
      }

    eps=std::min(0.000001,kz/10000);
    k=0;
    for (i=0; i < vertices; i++)
      { 
	if (fabs(l[i]-lg) < eps)
	  {
	    v[k]=i;
	    k++;
	  } 
      }
    for (i=k; i < vertices; i++)
      {
	v[i]=Global::getGeo(VERTEX,dim);
      }

    l[0]=lg;
  }

  void MacroReader::checkMesh(Mesh *mesh)
  {
    FUNCNAME("MacroReader::checkMesh()");
    
    int       i, nused, nfree;
    DOFAdmin  *localAdmin=mesh->admin[0];
    Flag     fill_flag;
    int       error_detected = 0;

    MSG("checking mesh\n");

    fill_flag = Mesh::CALL_EVERY_EL_INORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

    Mesh::traversePtr = mesh;
    error_detected += mesh->traverse(-1, fill_flag|Mesh::FILL_ADD_ALL, basicCheckFct);

    if (mesh->preserveCoarseDOFs)
      fill_flag = Mesh::CALL_EVERY_EL_INORDER;
    else 
      fill_flag = Mesh::CALL_LEAF_EL;
  
    fill_flag |= Mesh::FILL_NEIGH;

    for (int iadmin = 0; iadmin < static_cast<int>(mesh->admin.size()); iadmin++) {
      localAdmin = mesh->admin[iadmin];
   
      if (localAdmin->getSize() > 0) {
	if (static_cast<int>(mesh->dof_used.size()) < localAdmin->getSize()) {
	  mesh->dof_used.resize(localAdmin->getSize() + 1000);
	}
	for (i = 0; i < static_cast<int>(mesh->dof_used.size()); i++) 
	  mesh->dof_used[i] = 0;
	
	nused = nfree = 0;

	TraverseStack stack;
	ElInfo *elInfo = stack.traverseFirst(mesh, -1, fill_flag | Mesh::FILL_ADD_ALL);
	while (elInfo) {
	  basicDOFCheckFct(elInfo, mesh, iadmin);
	  elInfo = stack.traverseNext(elInfo);
	}
	
	DOFIteratorBase it(localAdmin, USED_DOFS);
	for (it.reset(); !it.end(); ++it) {
	  nused++;
	  if (!mesh->dof_used[it.getDOFIndex()])  {
	    error_detected++;
	    MSG("dof[%d] not used??\n",it.getDOFIndex());
	  }
	}
	
	DOFIteratorBase freeIt(localAdmin, FREE_DOFS);
	for (freeIt.reset(); !freeIt.end(); ++freeIt) {
	  nfree++;
	  if(mesh->dof_used[freeIt.getDOFIndex()]) {
	    error_detected++;
	    MSG("dof[%d] used??\n",freeIt.getDOFIndex());
	  }
	}
	
	TEST(nused + nfree == localAdmin->getSize())
	  ("nused = %d, nfree = %d, admin.size = %d ????\n",
	   nused, nfree, localAdmin->getSize());
	TEST(nused == localAdmin->getUsedDOFs())
	  ("nused = %d, admin.used_count = %d ?????\n",
	   nused, localAdmin->getUsedDOFs());
      }
    }

    if (!error_detected) {
      MSG("checking done; no error detected\n");
    } else {
      MSG("checking done; %d error%s detected\n", error_detected,
	  error_detected == 1 ? "" : "s");

      Mesh::traversePtr = mesh;
      mesh->traverse(-1, Mesh::CALL_EVERY_EL_INORDER, basicNodeFct);
      WAIT_REALLY;
    }
  }

  int MacroReader::basicCheckFct(ElInfo* elinfo)
  {
    FUNCNAME("MacroReader::basicCheckFct()");

    int j, k, opp_v;
    Element *el = elinfo->getElement();
    const Element *neig;
    int error_detected=0;

    Mesh* mesh = Mesh::traversePtr;

    int dim = mesh->getDim();

    elinfo->testFlag(Mesh::FILL_NEIGH);

    for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
      if ((neig = elinfo->getNeighbour(i))) {
	if(elinfo->getBoundary(i) > 0) { // < 0 => periodic boundary
	  if (!error_detected)
	    MSG("error detected!!!\n");
	  error_detected++;
	  MSG("interior (*boundary)[%d] non NULL on element = %d\n", 
	      i, el->getIndex());
	}
	
	opp_v = elinfo->getOppVertex(i);
	if(opp_v < 0  ||  opp_v >= mesh->getGeo(NEIGH)) {
	  if (!error_detected)
	    MSG("error detected!!!\n");
	  error_detected++;
	  MSG("opp_v = %d\n", opp_v);
	}
	
	if(elinfo->getBoundary(i) > 0) {  // < 0 => periodic boundary
	  if(dim == 1) {
	    if(el->getDOF(i) != neig->getDOF(opp_v)) {
	      if (!error_detected)
		MSG("error detected!!!\n");
	      error_detected++;
	      MSG("neighbour error\n");
	    }
	  } else {
	    for (j = 1; j < mesh->getGeo(VERTEX); j++) {
	      for (k = 1; k < mesh->getGeo(VERTEX); k++)
		if (el->getDOF((i+j) % mesh->getGeo(VERTEX)) == 
		    neig->getDOF((opp_v+k) % mesh->getGeo(VERTEX))) 
		  break;
	      
	      if (k >= mesh->getGeo(VERTEX)) {
		if (!error_detected)
		  MSG("error detected!!!\n");
		error_detected++;
		MSG("dof %d of el %d at face %d isn't dof of neigh %d at face %d\n",
		    el->getDOF((i+j) % 3,0), el->getIndex(), i, neig->getIndex(), 
		    opp_v);
	      }	   
	    }
	  }
	}
      } else {
	if (elinfo->getBoundary(i) == INTERIOR) {
	  if (!error_detected)
	    MSG("error detected!!!\n");
	  error_detected++;
	  MSG("(*boundary)[%d] on domains boundary is NULL on element = %d\n",
	      i, el->getIndex());
	}
      }
    }

    return error_detected;
  }

  void MacroReader::basicDOFCheckFct(ElInfo* elinfo, Mesh *mesh, int iadmin)
  {
    FUNCNAME("MacroReader::basicDOFCheckFct()");
    
    Element* el = elinfo->getElement();
    const DOFAdmin& admin = mesh->getDOFAdmin(iadmin);
    const Element *neig;
    const DegreeOfFreedom *dof;

    if (0 == mesh->dof_used.size()) 
      return;

    int ndof = admin.getNumberOfDOFs(VERTEX);
    if (ndof) {
      int j0 = admin.getNumberOfPreDOFs(VERTEX);
      TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(VERTEX))
	("admin.getNumberOfPreDOFs(VERTEX) %d + nDOF %d > mesh->nDOF %d\n",
	 j0, ndof, mesh->getNumberOfDOFs(VERTEX));
      int i0 = mesh->getNode(VERTEX);
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	if ((dof = el->getDOF(i0+i)) == NULL) {
	  ERROR("no vertex dof %d on element %d\n", i, el->getIndex());
	} else {
	  for (int j = 0; j < ndof; j++) {
	    int jdof = dof[j0 + j];
	    TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
	      ("vertex dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
	    mesh->dof_used[jdof]++;
	  }
	}
      }
      /* neighbour vertex dofs have been checked in check_fct() */
    }

  
    if (mesh->getDim() > 1) {
      ndof = admin.getNumberOfDOFs(EDGE);
      if (ndof) {
	int j0 = admin.getNumberOfPreDOFs(EDGE);

	TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(EDGE))
	  ("admin.getNumberOfPreDOFs(EDGE) %d + nDOF %d > mesh->nDOF %d\n",
	   j0, ndof, mesh->getNumberOfDOFs(EDGE));

	int i0 = mesh->getNode(EDGE);
	
	for (int i = 0; i < mesh->getGeo(EDGE); i++) {
	  dof = el->getDOF(i0 + i);

	  if (dof == NULL) {
	    ERROR("no edge dof %d on element %d\n", i, el->getIndex());
	  } else {
	    for (int j = 0; j < ndof; j++) {
	      int jdof = dof[j0 + j];
	      TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
		("edge dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
	      mesh->dof_used[jdof]++;
	    }
	  }

	  if (el->getFirstChild() == NULL) {
	    if (mesh->getDim() == 2) {
	      neig = elinfo->getNeighbour(i);
	      if (neig) {
		int ov = elinfo->getOppVertex(i);

		TEST(neig->getDOF(i0 + ov) == dof)
		  ("el %d edge %d dof %8X: wrong dof %8X in neighbour %d edge %d\n",
		   el->getIndex(), i, dof, neig->getDOF(i0 + ov), 
		   neig->getIndex(), ov);
	      }
	    } else { // dim == 3
	      for (int in = 0; in < mesh->getGeo(NEIGH); in++) {
		if ((in != el->getVertexOfEdge(i,0)) && 
		    (in != el->getVertexOfEdge(i,1)) &&
		    (neig = elinfo->getNeighbour(in))) {
		  int found = 0;
		  for (int k = 0; k < mesh->getGeo(EDGE); k++) {
		    if (neig->getDOF(i0 + k) == dof) found++;
		  }
		  TEST(found==1)("el %d edge %d dof found=%d in neighbour %d\n",
				 el->getIndex(), i, found, neig->getIndex());
		}
	      }
	    }
	  }
	}
      }
    }

    if (mesh->getDim() == 3) {
      ndof = admin.getNumberOfDOFs(FACE);
      if (ndof) {
	int j0 = admin.getNumberOfPreDOFs(FACE);
	TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(FACE))
	  ("admin->n0_dof[FACE] %d + nDOF %d > mesh->nDOF %d\n",
	   j0, ndof, mesh->getNumberOfDOFs(FACE));
	int i0 = mesh->getNode(FACE);
	for (int i = 0; i < mesh->getGeo(FACE); i++) {
	  TEST(dof = el->getDOF(i0 + i))("no face dof %d ???\n", i);
	  for (int j = 0; j < ndof; j++) {
	    int jdof = dof[j0 + j];
	    TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
	      ("face dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
	    mesh->dof_used[jdof]++;
	  }

	  if (el->getChild(0) == NULL) {
	    if ((neig = elinfo->getNeighbour(i))) {
	      int ov = elinfo->getOppVertex(i);

	      TEST(neig->getDOF(i0 + ov) == dof)
		("el %d face %d dof %8X: wrong dof %8X in neighbour %d face %d\n",
		 el->getIndex(), i, dof, neig->getDOF(i0 + ov), neig->getIndex(),
		 ov);
	    }
	  }
	}
      }
    }

    ndof = admin.getNumberOfDOFs(CENTER);
    if (ndof) {
      int i0 = mesh->getNode(CENTER);
      TEST(dof = el->getDOF(i0))("no center dof???\n");
      int j0 = admin.getNumberOfPreDOFs(CENTER);
      TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(CENTER))
      ("admin.getNumberOfPreDOFs(CENTER) %d + nDOF %d > mesh->nDOF %d\n",
       j0, ndof, mesh->getNumberOfDOFs(CENTER));
      for (int j = 0; j < ndof; j++) {
	int jdof = dof[j0 + j];
	TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
	  ("center dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
	mesh->dof_used[jdof]++;
      }
    }
  }

  int MacroReader::basicNodeFct(ElInfo* elinfo)
  {
    FUNCNAME("MacroReader::basicNodeFct()");

    Element *el = elinfo->getElement();
    Mesh *mesh = Mesh::traversePtr;

    MSG("el %4d: ", el->getIndex());
    for (int i = 0; i < mesh->getGeo(VERTEX); i++)
      Msg::print("%4d%s", el->getDOF(i,0), 
		 i < mesh->getGeo(VERTEX)-1 ? ", " : "\n");

    return 0;
  }
}