From 6640e0b5d606026d9df57d8c2c3bb71a453f4496 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 15 Apr 2008 13:58:25 +0000
Subject: [PATCH] * VTK Output possible for Lagrange elements of second order

---
 AMDiS/src/BasisFunction.h   |   4 +-
 AMDiS/src/DOFAdmin.cc       |  13 +-
 AMDiS/src/DataCollector.cc  | 238 ++++++++++---------
 AMDiS/src/DataCollector.h   |  14 +-
 AMDiS/src/ElInfo.cc         |   4 +-
 AMDiS/src/ElInfo.h          |   6 +-
 AMDiS/src/ElementInfo.h     |  10 +-
 AMDiS/src/FileWriter.cc     |  10 +-
 AMDiS/src/FiniteElemSpace.h |  24 +-
 AMDiS/src/Global.cc         |  13 +-
 AMDiS/src/Global.h          |  14 +-
 AMDiS/src/MacroWriter.cc    | 452 +-----------------------------------
 AMDiS/src/MacroWriter.h     | 238 ++++++++-----------
 AMDiS/src/ValueWriter.cc    | 214 ++---------------
 AMDiS/src/ValueWriter.h     | 142 +++++------
 AMDiS/src/VertexInfo.h      |  51 ++--
 AMDiS/src/VtkWriter.cc      | 238 +++++++++++++------
 AMDiS/src/VtkWriter.h       |  90 ++++++-
 18 files changed, 675 insertions(+), 1100 deletions(-)

diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index 0cd8a50a..b9393a11 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -96,14 +96,14 @@ namespace AMDiS {
      * compares two BasisFunction objects.
      */
     virtual bool operator==(const BasisFunction& a) const {
-      return a.getName()==name;
+      return a.getName() == name;
     };
 
     /** \brief
      * Returns !(*this == b)
      */
     inline bool operator!=(const BasisFunction& b) const {
-      return !operator==(b);
+      return !(operator == (b));
     };
 
     /** \brief
diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc
index 79290284..6783cab2 100755
--- a/AMDiS/src/DOFAdmin.cc
+++ b/AMDiS/src/DOFAdmin.cc
@@ -303,18 +303,21 @@ namespace AMDiS {
   }
 
   const int DOFAdmin::getNumberOfPreDOFs(int i) const { 
-    TEST_EXIT((0<=i)&&(4>i))("");
+    TEST_EXIT((0 <= i) && (4 > i))("");
+
     return nr0DOF[i]; 
   }
 
   void DOFAdmin::setNumberOfDOFs(int i,int v) { 
-    TEST_EXIT((0<=i)&&(4>i))("");
-    nrDOF[i]=v; 
+    TEST_EXIT((0 <= i) && (4 > i))("");
+
+    nrDOF[i] = v; 
   }
 
   void DOFAdmin::setNumberOfPreDOFs(int i, int v) { 
-    TEST_EXIT((0<=i)&&(4>i))(""); 
-    nr0DOF[i]=v; 
+    TEST_EXIT((0 <= i) && (4 > i))(""); 
+
+    nr0DOF[i] = v; 
   }
 
   DOFAdmin::~DOFAdmin()
diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index 2f377464..eb5299a3 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -13,7 +13,7 @@ namespace AMDiS {
 			       bool (*writeElem)(ElInfo*))
     : writeElem_(writeElem)
   {
-    FUNCNAME("DataCollector::DataCollector");
+    FUNCNAME("DataCollector::DataCollector()");
 
     TEST_EXIT(feSpace)("no feSpace\n");
     TEST_EXIT(values)("no value Vector\n");
@@ -27,6 +27,7 @@ namespace AMDiS {
     // === create vertex info vector ===
     vertexInfos_ = NEW DOFVector< ::std::list<VertexInfo> >(feSpace, "vertex infos");
     interpPointInd_ = NEW DOFVector<int>(feSpace, "interpolation point indices");
+    interpPointCoords_ = NEW DOFVector< ::std::list<WorldVector<double> > >(feSpace, "interpolation point coordinates");
     dofCoords_ = NEW DOFVector< ::std::list<WorldVector<double> > >(feSpace, "dof coords");
 
     dim_ = mesh_->getDim();
@@ -51,12 +52,13 @@ namespace AMDiS {
   {
     DELETE vertexInfos_;
     DELETE interpPointInd_;
+    DELETE interpPointCoords_;
     DELETE dofCoords_;
   }
 
   int DataCollector::startCollectingElementData()
   {
-    FUNCNAME("DataCollector::startCollectingElementData");
+    FUNCNAME("DataCollector::startCollectingElementData()");
 
     Flag flag = traverseFlag_;
     flag |= 
@@ -71,16 +73,16 @@ namespace AMDiS {
 
     // Traverse elements to create continuous element indices
     ElInfo *elInfo = stack.traverseFirst(mesh_, level_, flag);
-    while(elInfo) {
-      if(!writeElem_ || writeElem_(elInfo))
+    while (elInfo) {
+      if (!writeElem_ || writeElem_(elInfo))
 	outputIndices_[elInfo->getElement()->getIndex()] = nElements_++;
       elInfo = stack.traverseNext(elInfo);
     }
 
     // Traverse elements to create element information
     elInfo = stack.traverseFirst(mesh_, level_, flag);
-    while(elInfo) {
-      if(!writeElem_ || writeElem_(elInfo))
+    while (elInfo) {
+      if (!writeElem_ || writeElem_(elInfo))
 	addElementData(elInfo);
       elInfo = stack.traverseNext(elInfo);
     }
@@ -92,10 +94,10 @@ namespace AMDiS {
 
   int DataCollector::startCollectingValueData()
   {
-    FUNCNAME("DataCollector::startCollectingValueData");
+    FUNCNAME("DataCollector::startCollectingValueData()");
 
     DOFVector<int>::Iterator intPointIt(interpPointInd_, USED_DOFS);
-    for(intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
+    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
       (*intPointIt) = -1;
     }
 
@@ -109,7 +111,7 @@ namespace AMDiS {
 					 level_, 
 					 traverseFlag_ | Mesh::FILL_COORDS);
     while(elInfo) {
-      if(!writeElem_ || writeElem_(elInfo))
+      if (!writeElem_ || writeElem_(elInfo))
 	addValueData(elInfo);
       elInfo = stack.traverseNext(elInfo);
     }
@@ -117,16 +119,16 @@ namespace AMDiS {
     // Remove all interpolation marks and, instead, set to each
     // interpolation point its continous index starting from 0.
     int i = 0;
-    for(intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
-      if(*intPointIt == -3) {
+    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
+      if (*intPointIt == -3) {
 	*intPointIt = i++;
       }
     }
 
     // Traverse elements to create interpolation values.
     elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
-    while(elInfo) {
-      if(!writeElem_ || writeElem_(elInfo))
+    while (elInfo) {
+      if (!writeElem_ || writeElem_(elInfo))
 	addInterpData(elInfo);
       elInfo = stack.traverseNext(elInfo);
     }
@@ -138,22 +140,21 @@ namespace AMDiS {
 
   int DataCollector::startCollectingPeriodicData()
   {    
-    FUNCNAME("DataCollector::startCollectingPeriodicData");
+    FUNCNAME("DataCollector::startCollectingPeriodicData()");
 
     periodicConnections_.clear();
     
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
-    while(elInfo) {
-      if(!writeElem_ || writeElem_(elInfo)) {
+    while (elInfo) {
+      if (!writeElem_ || writeElem_(elInfo)) {
 	LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
 	  (elInfo->getElement()->
 	   getElementData()->
 	   getElementData(PERIODIC));
 	
-	if(ldp) {
-	  nConnections_ += 
-	    dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
+	if (ldp) {
+	  nConnections_ += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
 	}
 
 	periodicConnections_.push_back(DimVec<bool>(dim_, DEFAULT_VALUE, false));
@@ -173,8 +174,8 @@ namespace AMDiS {
       Mesh::FILL_BOUND;
 
     elInfo = stack.traverseFirst(mesh_, level_, flag);
-    while(elInfo) {
-      if(!writeElem_ || writeElem_(elInfo))
+    while (elInfo) {
+      if (!writeElem_ || writeElem_(elInfo))
 	addPeriodicData(elInfo);
       elInfo = stack.traverseNext(elInfo);
     }   
@@ -186,9 +187,8 @@ namespace AMDiS {
 
   int DataCollector::addElementData(ElInfo* elInfo)
   {
-    FUNCNAME("DataCollector::addElementData");
+    FUNCNAME("DataCollector::addElementData()");
 
-    int i;
     const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
     DegreeOfFreedom vertexDOF;
     WorldVector<double> vertexCoords;
@@ -197,10 +197,9 @@ namespace AMDiS {
     ElementInfo elementInfo(dim_);
     
     // read element region
-    ElementData *ed = 
-      elInfo->getElement()->getElementData(ELEMENT_REGION);
+    ElementData *ed = elInfo->getElement()->getElementData(ELEMENT_REGION);
     
-    if(ed) {
+    if (ed) {
       elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
     } else {
       elementInfo.elementRegion = -1;
@@ -211,14 +210,14 @@ namespace AMDiS {
     
     elementInfo.surfaceRegions.set(-1);
     
-    while(ed) {
+    while (ed) {
       SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
       elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
       ed = ed->getDecorated(SURFACE_REGION);
     }
 
     // for all vertices
-    for(i = 0; i < dim_ + 1; i++) {
+    for (int i = 0; i < dim_ + 1; i++) {
       // get coords of this vertex
       vertexCoords = elInfo->getCoord(i);
       
@@ -232,7 +231,7 @@ namespace AMDiS {
 	       vertexCoords);
       
       // coords not yet in list?
-      if(it == (*vertexInfos_)[vertexDOF].end()) {
+      if (it == (*vertexInfos_)[vertexDOF].end()) {
 	// create new vertex info and increment number of vertices
 	VertexInfo newVertexInfo = {vertexCoords, nVertices_};
 
@@ -255,7 +254,7 @@ namespace AMDiS {
 	-1;
     }
     
-    if(dim_ == 3) {
+    if (dim_ == 3) {
       elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());
     }
 
@@ -267,16 +266,14 @@ namespace AMDiS {
 
   int DataCollector::addValueData(ElInfo *elInfo)
   {
-    FUNCNAME("DataCollector::addValueData");
+    FUNCNAME("DataCollector::addValueData()");
 
-    int i, j, k;
-    int node_offset, dof_offset, num_dofs, node, dof_index;
     const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
 
     /* vertex dofs */
-    dof_offset  = localAdmin_->getNumberOfPreDOFs(VERTEX);
-    
-    for (i = 0; i < mesh_->getGeo(VERTEX); i++) {
+    int dof_offset = localAdmin_->getNumberOfPreDOFs(VERTEX);
+
+    for (int i = 0; i < mesh_->getGeo(VERTEX); i++) {
       (*interpPointInd_)[dof[i][dof_offset]] = -2; // mark as vertex
       
       // get coords of this vertex
@@ -289,53 +286,76 @@ namespace AMDiS {
 	       vertexCoords);
       
       // coords not yet in list?
-      if(it == (*dofCoords_)[dof[i][dof_offset]].end()) {
+      if (it == (*dofCoords_)[dof[i][dof_offset]].end()) {
 	// add new coords to list
 	(*dofCoords_)[dof[i][dof_offset]].push_back(vertexCoords);
       }
     }
 
-    for(i = 1; i <= dim_; i++) {
-      num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
-      node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
-      dof_offset  = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
-      for (j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
-	node = node_offset + j;
-	for(k = 0; k < num_dofs; k++) {
-	  dof_index = dof_offset + k;
+    int nInterpPoints = 0;
+    const BasisFunction *basisFcts = feSpace_->getBasisFcts();
+
+    for (int i = 1; i <= dim_; i++) {
+      int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
+      int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
+      dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
+
+      for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
+	int node = node_offset + j;
 
+	for (int k = 0; k < num_dofs; k++) {
+	  int dof_index = dof_offset + k;
+
+
+	  WorldVector<double> interpolCoords;
+	  elInfo->coordToWorld((*basisFcts->getCoords(mesh_->getGeo(VERTEX) + nInterpPoints)), 
+			       &interpolCoords);
+
+	  nInterpPoints++;
+	  
 	  if ((*interpPointInd_)[dof[node][dof_index]] == -1) {
-	    // mark as interp. point 
+	    // mark as interpolation point
 	    (*interpPointInd_)[dof[node][dof_index]] = -3; 
+
+	    // search for interpolation point coordinates, and insert them to the 
+	    // dof-entry, if not contained in the list
+	    ::std::list<WorldVector<double> >::iterator it =
+		find((*interpPointCoords_)[dof[node][dof_index]].begin(),
+		     (*interpPointCoords_)[dof[node][dof_index]].end(),
+		     interpolCoords);
+	    if (it == (*interpPointCoords_)[dof[node][dof_index]].end()) {
+	      (*interpPointCoords_)[dof[node][dof_index]].push_back(interpolCoords); 
+	    }
+
+
 	    nInterpPoints_++;
 	  }
 	}
       }
     }
-   
+  
     return(0);
   }
 
   int DataCollector::addInterpData(ElInfo *elInfo)
   {
-    FUNCNAME("DataCollector::addInterpData");
+    FUNCNAME("DataCollector::addInterpData()");
 
-    int i, j, k;
-    int node_offset, dof_offset, num_dofs, node, dof_index;
     const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
 
-    ::std::list<int> elemInterpPoints;
-
+    ::std::vector<int> elemInterpPoints;
     elemInterpPoints.clear();
 
-    for(i = 1; i <= dim_; i++) {
-      num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
-      node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
-      dof_offset  = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
-      for (j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
-	node = node_offset + j;
-	for(k = 0; k < num_dofs; k++) {
-	  dof_index = dof_offset + k;
+    for (int i = 1; i <= dim_; i++) {
+      int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
+      int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
+      int dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
+
+      for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
+	int node = node_offset + j;
+
+	for (int k = 0; k < num_dofs; k++) {
+	  int dof_index = dof_offset + k;
 	  
 	  elemInterpPoints.push_back((*interpPointInd_)[dof[node][dof_index]]);
 	}
@@ -355,57 +375,54 @@ namespace AMDiS {
        getElementData()->
        getElementData(PERIODIC));
     
-    if(ldp) {
+    if (ldp) {
       ::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
       
-      for(it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
-	  it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
-	  ++it)
-	{
-	  int outputIndex = outputIndices_[elInfo->getElement()->getIndex()];
-	  int neighIndex  = outputIndices_[elInfo->
-					   getNeighbour(it->elementSide)->
-					   getIndex()];
+      for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
+	   it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
+	   ++it) {
+
+	int outputIndex = outputIndices_[elInfo->getElement()->getIndex()];
+	int neighIndex  = outputIndices_[elInfo->
+					 getNeighbour(it->elementSide)->
+					 getIndex()];
+	
+	if (!periodicConnections_[outputIndex][it->elementSide]) {
+	  PeriodicInfo periodicInfo;
 	  
-	  if(!periodicConnections_[outputIndex][it->elementSide]) {
-	    PeriodicInfo periodicInfo;
-
-	    periodicInfo.mode = it->periodicMode;
-	    periodicInfo.type = it->type;
-	    periodicInfo.outputIndex = outputIndex;
-	    periodicInfo.neighIndex = neighIndex;
-	    periodicInfo.vertexMap.clear();
-
-	    periodicConnections_[outputIndex][it->elementSide] = true;
-	    periodicConnections_
-	      [neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
-
-	    int index1, index2, dof1, dof2, i, j;
-
-	    for(i = 0; i < dim_; i++) {
-	      index1 = 
-		elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
-							  it->elementSide,
-							  i);
-	      dof1 = elInfo->getElement()->getDOF(index1, nPreDofs_);
-      
-	      for(j = 0; j < dim_; j++) {
-		index2 = 
-		  elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
-							    elInfo->getOppVertex(it->elementSide),
-							    j);
-		dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, nPreDofs_);
-
-		if((dof1 == dof2) || (mesh_->associated(dof1, dof2))) {
-		  periodicInfo.vertexMap[index1] = index2;
-		  break;
-		} 
-	      }
+	  periodicInfo.mode = it->periodicMode;
+	  periodicInfo.type = it->type;
+	  periodicInfo.outputIndex = outputIndex;
+	  periodicInfo.neighIndex = neighIndex;
+	  periodicInfo.vertexMap.clear();
+	  
+	  periodicConnections_[outputIndex][it->elementSide] = true;
+	  periodicConnections_
+	    [neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
+	    
+  
+	  for (int i = 0; i < dim_; i++) {
+	    int index1 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
+								   it->elementSide,
+								   i);
+	    int dof1 = elInfo->getElement()->getDOF(index1, nPreDofs_);
+	    
+	    for (int j = 0; j < dim_; j++) {
+	      int index2 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
+								     elInfo->getOppVertex(it->elementSide),
+								     j);
+	      int dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, nPreDofs_);
+	      
+	      if ((dof1 == dof2) || (mesh_->associated(dof1, dof2))) {
+		periodicInfo.vertexMap[index1] = index2;
+		break;
+	      } 
 	    }
-
-	    periodicInfos_.push_back(periodicInfo);
 	  }
+	  
+	  periodicInfos_.push_back(periodicInfo);
 	}
+      }
     }
       
     return(0);
@@ -511,7 +528,16 @@ namespace AMDiS {
     return interpPointInd_;
   }
 
-  ::std::list< ::std::list<int> >* DataCollector::getInterpPoints()
+  DOFVector< ::std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
+  {
+    if (!valueDataCollected_) {
+      startCollectingValueData();
+    }
+
+    return interpPointCoords_;
+  }
+
+  ::std::vector< ::std::vector<int> >* DataCollector::getInterpPoints()
   {
     if (!valueDataCollected_) {
       startCollectingValueData();
diff --git a/AMDiS/src/DataCollector.h b/AMDiS/src/DataCollector.h
index a531bad3..a7a23fb1 100644
--- a/AMDiS/src/DataCollector.h
+++ b/AMDiS/src/DataCollector.h
@@ -92,10 +92,15 @@ namespace AMDiS {
        */          
       DOFVector<int>* getInterpPointInd();
 
+      /** \brief
+       *
+       */
+      DOFVector< ::std::list<WorldVector<double> > >* getInterpPointCoords();
+
       /** \brief
        * Returns list of interpolation point information.
        */
-      ::std::list< ::std::list<int> >* getInterpPoints();
+      ::std::vector< ::std::vector<int> >* getInterpPoints();
 
       /** \brief
        * Returns list of information about periodics.
@@ -237,7 +242,12 @@ namespace AMDiS {
       /** \brief
        * Stores for each simplex the interpolation points.
        */
-      ::std::list< ::std::list<int> > interpPoints_;
+      ::std::vector< ::std::vector<int> > interpPoints_;
+
+      /** \brief
+       *
+       */
+      DOFVector< ::std::list<WorldVector<double> > > *interpPointCoords_;
 
       /** \brief
        * list of coords for each dof
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index bf92f4e0..48bae890 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -71,7 +71,7 @@ namespace AMDiS {
     double c = l[0];
 
     for (int j = 0; j < dimOfWorld; j++)
-      (*ret)[j] = c*v[j];
+      (*ret)[j] = c * v[j];
   
     int vertices =  Global::getGeo(VERTEX, dim);
 
@@ -79,7 +79,7 @@ namespace AMDiS {
       v = (*coords)[i];
       c = l[i];
       for (int j = 0; j < dimOfWorld; j++)
-	(*ret)[j] += c*v[j];
+	(*ret)[j] += c * v[j];
     }
     return ret;
   }
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 3e02ca56..94c42fe9 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -368,15 +368,15 @@ namespace AMDiS {
     void testFlag(const Flag& flag) const;
 
     /** \brief
-     * Returns a pointer to a vector, wich contains the world coordinates
+     * Returns a pointer to a vector, which contains the world coordinates
      * of a point in barycentric coordinates lambda with respect to \ref element.
      * If world is not NULL the world coordinates are stored in this vector.
      * Otherwise the function itself provides memory for this vector. In this
      * case the vector is overwritten during the next call of coordToWorld.
      */
     virtual const WorldVector<double>
-    *coordToWorld(const DimVec<double>& lambda,
-		  WorldVector<double>* world) const;
+      *coordToWorld(const DimVec<double>& lambda,
+		    WorldVector<double>* world) const;
   
 
     /** \brief
diff --git a/AMDiS/src/ElementInfo.h b/AMDiS/src/ElementInfo.h
index d1a15bf8..c4b21de6 100644
--- a/AMDiS/src/ElementInfo.h
+++ b/AMDiS/src/ElementInfo.h
@@ -77,9 +77,15 @@ namespace AMDiS {
        * Element type. Used in 3d.
        */
       unsigned char type;
-      
+
+      /** \brief
+       *
+       */
       int elementRegion;
-      
+
+      /** \brief
+       *
+       */      
       DimVec<int> surfaceRegions;
     };
   
diff --git a/AMDiS/src/FileWriter.cc b/AMDiS/src/FileWriter.cc
index 6b08fb18..a9d30f3a 100644
--- a/AMDiS/src/FileWriter.cc
+++ b/AMDiS/src/FileWriter.cc
@@ -209,13 +209,17 @@ namespace AMDiS {
     }
     
     if (writeParaViewFormat) {
-      VtkWriter::writeFile(&dc, const_cast<char*>( (fn + paraViewFileExt).c_str()));
+      VtkWriter vtkWriter(&dc);
+
+      vtkWriter.writeFile(const_cast<char*>( (fn + paraViewFileExt).c_str()));
       MSG("ParaView file written to %s\n", (fn + paraViewFileExt).c_str());
     }
 
     if (writeParaViewAnimation) {
-      VtkWriter::updateAnimationFile(fn + paraViewFileExt, 
-				     &paraViewAnimationFrames_, 
+      VtkWriter vtkWriter(&dc);
+
+      vtkWriter.updateAnimationFile(fn + paraViewFileExt, 
+				    &paraViewAnimationFrames_, 
 				     const_cast<char*>( (filename + ".pvd").c_str()));
     }
 
diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h
index d2d9c9ea..86bd386e 100644
--- a/AMDiS/src/FiniteElemSpace.h
+++ b/AMDiS/src/FiniteElemSpace.h
@@ -61,10 +61,10 @@ namespace AMDiS {
     /** \brief
      * 
      */
-    static FiniteElemSpace *provideFESpace(DOFAdmin            *admin,
+    static FiniteElemSpace *provideFESpace(DOFAdmin *admin,
 					   const BasisFunction *basFcts,
-					   Mesh                *mesh,
-					   const ::std::string&   name_="");
+					   Mesh *mesh,
+					   const ::std::string& name_ = "");
 
     /** \brief
      * destructor
@@ -74,22 +74,30 @@ namespace AMDiS {
     /** \brief
      * Returns \ref name
      */  
-    inline ::std::string getName() const { return name;};
+    inline ::std::string getName() const { 
+      return name;
+    };
 
     /** \brief
      * Returns \ref admin
      */
-    inline DOFAdmin* getAdmin() const { return admin;};
+    inline DOFAdmin* getAdmin() const { 
+      return admin;
+    };
 
     /** \brief
      * Returns \ref basFcts
      */
-    inline const BasisFunction* getBasisFcts() const { return basFcts;};
+    inline const BasisFunction* getBasisFcts() const { 
+      return basFcts;
+    };
 
     /** \brief
      * Returns \ref mesh
      */
-    inline Mesh* getMesh() const { return mesh; };
+    inline Mesh* getMesh() const { 
+      return mesh; 
+    };
   
   protected:
     /** \brief
@@ -120,7 +128,7 @@ namespace AMDiS {
     /** \brief
      * The Mesh this FiniteElemSpace belongs to
      */
-    Mesh*                 mesh;
+    Mesh* mesh;
 
     /** \brief
      * 
diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc
index be40548b..ec0f6da8 100644
--- a/AMDiS/src/Global.cc
+++ b/AMDiS/src/Global.cc
@@ -297,15 +297,18 @@ namespace AMDiS {
 
   int Global::getGeo(GeoIndex p, int dim) 
   {
-    FUNCNAME("Global::getGeo");
+    FUNCNAME("Global::getGeo()");
+
     initTest();
-    TEST_EXIT((p>=MINPART)&&(p<=MAXPART))
+    TEST_EXIT((p >= MINPART) && (p <= MAXPART))
       ("Calling for invalid geometry value %d\n",p);
-    TEST_EXIT((dim>=0)&&(dim<4))("invalid dim: %d\n", dim);
+    TEST_EXIT((dim >= 0) && (dim < 4))
+      ("invalid dim: %d\n", dim);
 
-    if(p == WORLD) return dimOfWorld;
+    if (p == WORLD) 
+      return dimOfWorld;
 
-    if(dim == 0) {
+    if (dim == 0) {
       switch(p) {
       case PARTS:
 	return 1;
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index 77d8cd81..97cf6571 100644
--- a/AMDiS/src/Global.h
+++ b/AMDiS/src/Global.h
@@ -428,7 +428,7 @@ namespace AMDiS {
      * can get information about the element via Element's getGeo method.
      */
     static const Element *getReferenceElement(int dim) {
-      FUNCNAME("Global::getReferenceElement");
+      FUNCNAME("Global::getReferenceElement()");
       initTest();
       TEST_EXIT((dim > 0) && (dim < 4))("invalid dim\n");
       return referenceElement[dim];
@@ -438,10 +438,13 @@ namespace AMDiS {
      * returns geometrical information. Currently this is only dimOfWorld.
      */
     static inline int getGeo(GeoIndex p) {
-      FUNCNAME("Global::getGeo");
+      FUNCNAME("Global::getGeo()");
+
       initTest();
-      if (WORLD==p) return dimOfWorld; 
-      ERROR_EXIT("Illegal request for geometry data: part=%d!\n",p);
+      if (WORLD == p) 
+	return dimOfWorld; 
+
+      ERROR_EXIT("Illegal request for geometry data: part=%d!\n", p);
       return 0;
     };
 
@@ -467,7 +470,8 @@ namespace AMDiS {
      * calls init if Global is not yet initialized.
      */
     static void initTest() {
-      if(dimOfWorld == 0) init();
+      if (dimOfWorld == 0) 
+	init();
     };
 
   private:
diff --git a/AMDiS/src/MacroWriter.cc b/AMDiS/src/MacroWriter.cc
index 78369f63..6dcdad11 100644
--- a/AMDiS/src/MacroWriter.cc
+++ b/AMDiS/src/MacroWriter.cc
@@ -40,7 +40,8 @@ namespace AMDiS {
 			      Flag traverseFlag,
 			      bool (*writeElem)(ElInfo*))
   {
-    FUNCNAME("MacroWroter::writeFile");
+    FUNCNAME("MacroWroter::writeFile()");
+
     TEST_EXIT(dc)("no data collector\n");
 
     writeElement = writeElem;
@@ -220,453 +221,4 @@ namespace AMDiS {
     file.close();
   }
 
-
-  int MacroWriter::writeMacroOld(const FiniteElemSpace *fe_space, 
-				 const char* name, double time,
-				 int level,
-				 Flag traverseFlag)
-  {
-    FUNCNAME("MacroWriter::writeMacro");
-    TEST_EXIT(fe_space)("no feSpace\n");
-
-    mesh = fe_space->getMesh();
-    dim = mesh->getDim();
-    int dow = Global::getGeo(WORLD);
-
-    bool writeElType = (dim == 3);
-    int vertices = mesh->getGeo(VERTEX);
-
-    // === create filename ===
-    ::std::string filename;
-    if(name) {
-      filename = ::std::string(name);
-    } else {
-      filename = mesh->getName();
-    }
-  
-    // === open file ===
-    if (!(macroFile = fopen(filename.c_str(), "w"))) {
-      ERROR_EXIT("could not open file %s for writing\n", filename.c_str());
-      return(0);
-    }
-
-    // === get admin and fe_space ===
-    DOFAdmin* localAdmin = const_cast<DOFAdmin*>(fe_space->getAdmin());
-
-    //   DOFAdmin* localAdmin = const_cast<DOFAdmin*>(mesh->getVertexAdmin());
-    //   TEST_EXIT(localAdmin)("no dof admin for vertices\n");
-    //   FiniteElemSpace *fe_space = 
-    //     FiniteElemSpace::provideFESpace(localAdmin, NULL, NULL);
-
-    // === create vertex info vector ===
-    vertexInfos = NEW DOFVector< ::std::list<VertexInfo> >(fe_space, "vertex infos");
-
-    // === delete old element list ===
-    elements.resize(0, dim);
-
-    n0 = localAdmin->getNumberOfPreDOFs(VERTEX);
-    nv = ne = 0;
-
-    // create output indices
-    //mesh->traverse(level, traverseFlag, indexFct);
-
-    TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
-    while(elInfo) {
-      if(!writeElement || writeElement(elInfo)) indexFct(elInfo);
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    traverseFlag |= 
-      Mesh::FILL_NEIGH      | 
-      Mesh::FILL_COORDS     |
-      Mesh::FILL_OPP_COORDS |
-      Mesh::FILL_BOUND;
-
-    // === traverse mesh and fill vertex infos and element list ===
-//     mesh->traverse(level, 
-// 		   traverseFlag,
-// 		   traverseFct);
-
-    elInfo = stack.traverseFirst(mesh, level, traverseFlag);
-    while(elInfo) {
-      if(!writeElement || writeElement(elInfo)) traverseFct(elInfo);
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    // === print file header ===
-    fprintf(macroFile, "mesh name: %s\n\n", mesh->getName().c_str());
-    fprintf(macroFile, "time: %e\n\n", time);
-
-    fprintf(macroFile, "DIM: %d\n", dim);
-    fprintf(macroFile, "DIM_OF_WORLD: %d\n\n", dow);
-
-    fprintf(macroFile, "number of vertices: %d\n", nv);
-    fprintf(macroFile, "number of elements: %d\n\n", ne);
-
-    // === print vertex coords and remember global output indices ===
-    fprintf(macroFile, "vertex coordinates:\n");
-
-    DOFVector< ::std::list<VertexInfo> >::Iterator it(vertexInfos, USED_DOFS);
-    int i, counter = 0;
-
-    // for all DOFs
-    for(it.reset(); !it.end(); ++it) {
-      // for all vertex infos of this DOF
-      ::std::list<VertexInfo>::iterator it2;
-      for(it2 = it->begin(); it2 != it->end(); ++it2) {
-	it2->outputIndex = counter++;
-	for(i=0; i < dow; i++) {
-	  fprintf(macroFile, "%e ", it2->coords[i]);
-	}
-	fprintf(macroFile, "\n");
-      }
-    }
-
-    // === print element vertices ===
-    fprintf(macroFile, "\nelement vertices:\n");
-
-    // iterate the element list
-    ::std::list<ElementInfo>::iterator elementIt;
-
-    for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-      // for all vertices
-      for(i=0; i < vertices; i++) {
-	fprintf(macroFile, "%d ", elementIt->vertexInfo[i]->outputIndex);
-      }
-      fprintf(macroFile, "\n");
-    }
-
-    // === print boundaries ===
-    fprintf(macroFile, "\nelement boundaries:\n");
-
-    for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-      // for all vertices
-      for(i=0; i < vertices; i++) {
-	fprintf(macroFile, "%d ", elementIt->boundary[i]);
-      }
-      fprintf(macroFile, "\n");
-    }
-
-    // === print neighbours ===
-    fprintf(macroFile, "\nelement neighbours:\n");
-
-    for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-      // for all vertices
-      for(i=0; i < vertices; i++) {
-	fprintf(macroFile, "%d ", elementIt->neighbour[i]);
-      }
-      fprintf(macroFile, "\n");
-    }
-
-    // === print boundary projections ===
-    fprintf(macroFile, "\nprojections:\n");
-
-    for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-      // for all vertices
-      for(i=0; i < vertices; i++) {
-	if(elementIt->projection[i])
-	  fprintf(macroFile, "%d ", elementIt->projection[i]->getID());
-	else
-	  fprintf(macroFile, "0 ");
-      }
-      fprintf(macroFile, "\n");
-    }
-
-    // === print element regions ===
-    fprintf(macroFile, "\nelement region:\n");
-
-    for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-      fprintf(macroFile, "%d\n", elementIt->elementRegion);
-    }
-
-    // === print surface regions ===
-    fprintf(macroFile, "\nsurface region:\n");
-
-    for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-      for(i=0; i < vertices; i++) {
-	fprintf(macroFile, "%d ", elementIt->surfaceRegions[i]);
-      }
-      fprintf(macroFile, "\n");
-    }
-
-    // === print element types if necessary ===
-    if (dim == 3) {
-      if (writeElType) {
-	fprintf(macroFile, "\nelement type:\n");
-
-	for(elementIt = elements.begin(); elementIt != elements.end(); ++elementIt) {
-	  fprintf(macroFile, "%d ", elementIt->type);
-	}
-      }
-    }
-
-    // === close file ===
-    fclose(macroFile);
-    macroFile = NULL;
-
-    DELETE vertexInfos;
-
-    return(1);
-  }
-
-  int MacroWriter::traverseFct(ElInfo* elInfo)
-  {
-    int i;
-    const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
-    DegreeOfFreedom vertexDOF;
-    WorldVector<double> vertexCoords;
-
-    // increment number of elements
-    //ne++;
-
-    // create ElementInfo
-    ElementInfo elementInfo(dim);
-
-    // read element region
-    ElementData *ed = 
-      elInfo->getElement()->getElementData(ELEMENT_REGION);
-  
-    if(ed) {
-      elementInfo.elementRegion = 
-	dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
-    }
-    else {
-      elementInfo.elementRegion = -1;
-    }
-
-    // read surface regions to element info
-    ed = elInfo->getElement()->getElementData(SURFACE_REGION);
-
-    elementInfo.surfaceRegions.set(-1);
-
-    while(ed) {
-      SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
-      elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
-      ed = ed->getDecorated(SURFACE_REGION);
-    }
-
-    // for all vertices
-    for(i=0; i < dim+1; i++) {
-      // get coords of this vertex
-      vertexCoords = elInfo->getCoord(i);
-
-      // get dof index of this vertex
-      vertexDOF = dof[i][n0];
-
-      // search for coords at this dof
-      ::std::list<VertexInfo>::iterator it =
-	  find((*vertexInfos)[vertexDOF].begin(),
-	       (*vertexInfos)[vertexDOF].end(),
-	       vertexCoords);
-
-      // coords not yet in list?
-      if(it == (*vertexInfos)[vertexDOF].end()) {
-	// create new vertex info
-	VertexInfo newVertexInfo = {vertexCoords, -1};
-      
-	// add new vertex info to list
-	(*vertexInfos)[vertexDOF].push_front(newVertexInfo);
-      
-	// set iterator to new vertex info
-	it = (*vertexInfos)[vertexDOF].begin();
-
-	// increment number of vertices
-	nv++;
-      }
-
-      // fill element info
-      elementInfo.vertexInfo[i] = it;
-      elementInfo.boundary[i] = elInfo->getBoundary(i);
-      elementInfo.projection[i] = elInfo->getProjection(i);
-      elementInfo.neighbour[i] = 
-	elInfo->getNeighbour(i) ?
-	outputIndices[elInfo->getNeighbour(i)->getIndex()] :
-	-1;
-    }
-
-    if(dim == 3) {
-      elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());
-    }
-
-    // remember element info
-    elements.push_back(elementInfo);
-
-    return 0;
-  }
-
-  void MacroWriter::writePeriodicFileOld(Mesh *aMesh, 
-					 const ::std::string filename,
-					 int level,
-					 Flag traverseFlag)
-  {
-    // === open file ===
-    if (!(periodicFile = fopen(filename.c_str(), "w"))) {
-      ERROR_EXIT("could not open file %s for writing\n", filename.c_str());
-    }
-
-    mesh = aMesh;
-    dim = mesh->getDim();
-
-    // === create output indices and count connections ===
-
-    n0 = mesh->getVertexAdmin()->getNumberOfPreDOFs(VERTEX);
-    ne = 0;
-    nc = 0;
-
-    //mesh->traverse(level, traverseFlag, periodicFct1);
-
-    TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
-    while(elInfo) {
-      if(!writeElement || writeElement(elInfo)) periodicFct1(elInfo);
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    nc /= 2;
-
-    fprintf(periodicFile, "associations: %d\n", nc);
-    fprintf(periodicFile, "\nmode  bc  el1 - local vertices <->  el2 - local vertices\n");
-  
-    traverseFlag |=
-      Mesh::FILL_COORDS    |
-      Mesh::FILL_OPP_COORDS|
-      Mesh::FILL_NEIGH     | 
-      Mesh::FILL_BOUND;
-    
-
-    //mesh->traverse(level, 
-    //	   traverseFlag,
-    //	   periodicFct2);
-
-    elInfo = stack.traverseFirst(mesh, level, traverseFlag);
-    while(elInfo) {
-      if(!writeElement || writeElement(elInfo)) periodicFct2(elInfo);
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    outputIndices.clear();
-    periodicConnections.clear();
-
-    // === close periodic file ===
-
-    fclose(periodicFile);
-    periodicFile = NULL;
-  }
-
-  int MacroWriter::periodicFct1(ElInfo *elInfo)
-  {
-    const Element *element = elInfo->getElement();
-
-    outputIndices[element->getIndex()] = ne++;
-
-    LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
-      (element->
-       getElementData()->
-       getElementData(PERIODIC));
-
-    if(ldp) {
-      nc += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
-    }
-
-    periodicConnections.push_back(DimVec<bool>(dim, DEFAULT_VALUE, false));
-
-    return 0;
-  }
-
-  int MacroWriter::periodicFct2(ElInfo *elInfo)
-  {
-    FUNCNAME("MacroWriter::periodicFct2");
-
-    int i, j;
-
-    const Element *element = elInfo->getElement(); 
-
-    LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
-      (element->
-       getElementData()->
-       getElementData(PERIODIC));
-
-    if(ldp) {
-      ::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
-      for(it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
-	  it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
-	  ++it) 
-	{
-	  int outputIndex = outputIndices[elInfo->getElement()->getIndex()];
-	  int neighIndex  = outputIndices[elInfo->
-					  getNeighbour(it->elementSide)->
-					  getIndex()];
-
-	  if(!periodicConnections[outputIndex][it->elementSide]) {
-
-	    fprintf(periodicFile, "%d ", it->periodicMode);
-	    fprintf(periodicFile, "%d ", it->type);
-
-
-	    periodicConnections[outputIndex][it->elementSide] = true;
-	    periodicConnections
-	      [neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
-
-	    int index1, index2, dof1, dof2;
-	    ::std::map<int, int> vertexMap;
-
-	    for(i=0; i < dim; i++) {
-	      index1 = 
-		elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim-1, dim),
-							  it->elementSide,
-							  i);
-	      dof1 = elInfo->getElement()->getDOF(index1, n0);
-
-	      for(j=0; j < dim; j++) {
-		index2 = 
-		  elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim-1, dim),
-							    elInfo->getOppVertex(it->elementSide),
-							    j);
-		dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, n0);
-
-		if((dof1 == dof2) || (mesh->associated(dof1, dof2))) {
-		  vertexMap[index1] = index2;
-		  //MSG("%d ", outputIndices[elInfo->getElement()->getIndex()]); 
-
-		  //	      ::std::cout << index1 << " " << index2 << " ";
-		  break;
-		}
-	      }
-	    }
-
-	    fprintf(periodicFile, "%d ", outputIndex);
-	    for(i=0; i < dim; i++) {
-	      fprintf(periodicFile, "%d ", 
-		      elInfo->
-		      getElement()->
-		      getVertexOfPosition(INDEX_OF_DIM(dim-1, dim),
-					  it->elementSide,
-					  i)
-		      );
-	    }
-	    fprintf(periodicFile, "%d ", neighIndex);
-	    for(i=0; i < dim; i++) {
-	      fprintf(periodicFile, "%d", 
-		      vertexMap[elInfo->
-				getElement()->
-				getVertexOfPosition(INDEX_OF_DIM(dim-1, dim),
-						    it->elementSide,
-						    i)]
-		      );
-	      fprintf(periodicFile, i < dim-1 ? " " : "\n");
-	    }
-	  }
-	}
-    }
-
-    return 0;
-  }
-
-  int MacroWriter::indexFct(ElInfo *elInfo)
-  {
-    outputIndices[elInfo->getElement()->getIndex()] = ne++;  
-    return 0;
-  }
-
 }
diff --git a/AMDiS/src/MacroWriter.h b/AMDiS/src/MacroWriter.h
index 8ea7b4b2..76847333 100644
--- a/AMDiS/src/MacroWriter.h
+++ b/AMDiS/src/MacroWriter.h
@@ -40,141 +40,109 @@ namespace AMDiS {
   class VertexInfo;
   class ElementInfo;
 
-class Mesh;
-class ElInfo;
-class FiniteElemSpace;
-template<typename T> class DOFVector;
-template<typename T> class DimVec;
-
-/** 
- * \ingroup Output
- *
- * \brief
- * Writes the current leaf elements of a mesh as macro triangulation to
- * a text file. Pure static class.
- */
-class MacroWriter
-{
-public:
-  /** \brief
-   * Stores a list of vertex infos for each dof.
-   */
-  static DOFVector< ::std::list<VertexInfo> > *vertexInfos;
-
-  /** \brief
-   * List that stores an ElementInfo for each element.
-   */
-  static ::std::list<ElementInfo> elements;
-
-public:
-  MEMORY_MANAGED(MacroWriter);
-
-  /** \brief
-   * Writes the leaf elements of a Mesh as a macro triangulation to a file.
-   */
-  static int writeMacro(DataCollector *dc,
-			const char *name, 		       
-			double time = 0.0,
-			int level = -1,
-			Flag traverseFlag = Mesh::CALL_LEAF_EL,
-			bool (*writeElem)(ElInfo*) = NULL);
-
-  /** \brief
-   * Init \ref periodicFile for the next macro to be written.
-   */
-  static void writePeriodicFile(DataCollector *dc,
-				const ::std::string filename);
-
-
-  // OBSOLETE
-  // TODO: Remove if MacroWriter::writeMacro is stable.
-  static int writeMacroOld(const FiniteElemSpace *aFESpace, 
-			   const char *name = NULL, 
-			   double time = 0.0,
-			   int level = -1,
-			   Flag traverseFlag = Mesh::CALL_LEAF_EL);
-  
-  // OBSOLETE
-  // TODO: Remove if MacroWriter::writePeriodic is stable.  
-  static void writePeriodicFileOld(Mesh *aMesh, const ::std::string filename,
-				   int level = -1,
-				   Flag traverseFlag = Mesh::CALL_LEAF_EL);
-  
-protected:
-  /** \brief
-   * creates output indices for each element.
-   */
-  static int indexFct(ElInfo *elInfo);
-  
-  /** \brief
-   * Creates vertex infos for 
-   */
-  static int traverseFct(ElInfo* elInfo);
-
-  /** \brief
-   * Used to generate global output indices and count the number of connections.
-   */
-  static int periodicFct1(ElInfo *elInfo);
-
-  /** \brief
-   * Writes periodic file entry for this element.
-   */
-  static int periodicFct2(ElInfo *elInfo);
-
-protected:
-  /** \brief
-   * Mesh that should be written
-   */
-  static Mesh *mesh;
-
-  /** \brief
-   * File to which the mesh should be written
-   */
-  static FILE *macroFile;
-
-  /** \brief
-   * File in which the periodic infos are stored.
-   */
-  static FILE *periodicFile;
-
-  /** \brief
-   * vertex pre-dofs
-   */
-  static int n0;
-
-  /** \brief
-   * Number of vertices.
-   */
-  static int nv;
-
-  /** \brief
-   * Number of elements.
-   */
-  static int ne;
-
-  /** \brief
-   * Number of connections.
-   */
-  static int nc;
-
-  /** \brief
-   * Dimension of \ref mesh
-   */
-  static int dim;
-
-  /** \brief
-   * Maps internal element indices to global output indices.
-   */
-  static ::std::map<int, int> outputIndices;
-
-  /** \brief
-   * periodicConnections[i][j] stores whether the connection at side j of 
-   * the element with output index i has already been written.
-   */
-  static ::std::vector<DimVec<bool> > periodicConnections;
-
-  static bool (*writeElement)(ElInfo*);
-};
+  class Mesh;
+  class ElInfo;
+  class FiniteElemSpace;
+  template<typename T> class DOFVector;
+  template<typename T> class DimVec;
+
+  /** 
+   * \ingroup Output
+   *
+   * \brief
+   * Writes the current leaf elements of a mesh as macro triangulation to
+   * a text file. Pure static class.
+   */
+  class MacroWriter
+  {
+  public:
+    /** \brief
+     * Stores a list of vertex infos for each dof.
+     */
+    static DOFVector< ::std::list<VertexInfo> > *vertexInfos;
+
+    /** \brief
+     * List that stores an ElementInfo for each element.
+     */
+    static ::std::list<ElementInfo> elements;
+
+  public:
+    MEMORY_MANAGED(MacroWriter);
+
+    /** \brief
+     * Writes the leaf elements of a Mesh as a macro triangulation to a file.
+     */
+    static int writeMacro(DataCollector *dc,
+			  const char *name, 		       
+			  double time = 0.0,
+			  int level = -1,
+			  Flag traverseFlag = Mesh::CALL_LEAF_EL,
+			  bool (*writeElem)(ElInfo*) = NULL);
+
+    /** \brief
+     * Init \ref periodicFile for the next macro to be written.
+     */
+    static void writePeriodicFile(DataCollector *dc,
+				  const ::std::string filename);
+
+ 
+  protected:
+    /** \brief
+     * Mesh that should be written
+     */
+    static Mesh *mesh;
+
+    /** \brief
+     * File to which the mesh should be written
+     */
+    static FILE *macroFile;
+
+    /** \brief
+     * File in which the periodic infos are stored.
+     */
+    static FILE *periodicFile;
+
+    /** \brief
+     * vertex pre-dofs
+     */
+    static int n0;
+
+    /** \brief
+     * Number of vertices.
+     */
+    static int nv;
+
+    /** \brief
+     * Number of elements.
+     */
+    static int ne;
+
+    /** \brief
+     * Number of connections.
+     */
+    static int nc;
+
+    /** \brief
+     * Dimension of \ref mesh
+     */
+    static int dim;
+
+    /** \brief
+     * Maps internal element indices to global output indices.
+     */
+    static ::std::map<int, int> outputIndices;
+
+    /** \brief
+     * periodicConnections[i][j] stores whether the connection at side j of 
+     * the element with output index i has already been written.
+     */
+    static ::std::vector<DimVec<bool> > periodicConnections;
+
+    /** \brief
+     *
+     */
+    static bool (*writeElement)(ElInfo*);
+  };
 
 }
 
diff --git a/AMDiS/src/ValueWriter.cc b/AMDiS/src/ValueWriter.cc
index 5772d628..4d602784 100644
--- a/AMDiS/src/ValueWriter.cc
+++ b/AMDiS/src/ValueWriter.cc
@@ -15,184 +15,6 @@ namespace AMDiS {
   DOFVector< ::std::list<WorldVector<double> > >* ValueWriter::dofCoords = NULL;
   int ValueWriter::ni = 0;
 
-  int ValueWriter::markFct(ElInfo *el_info) {
-    int i, j, k;
-    int node_offset, dof_offset, num_dofs, node, dof_index;
-    const DegreeOfFreedom **dof = el_info->getElement()->getDOF();
-    const Mesh *mesh = admin->getMesh();
-    int dim = mesh->getDim();
-
-    /* vertex dofs */
-    dof_offset  = admin->getNumberOfPreDOFs(VERTEX);
-  
-    for (i = 0; i < mesh->getGeo(VERTEX); i++) {
-      (*interpPointInd)[dof[i][dof_offset]] = -2; // mark as vertex
-
-      // get coords of this vertex
-      WorldVector<double> vertexCoords = el_info->getCoord(i);
-
-      // search for coords at this dof
-      ::std::list<WorldVector<double> >::iterator it =
-	  find((*dofCoords)[dof[i][dof_offset]].begin(),
-	       (*dofCoords)[dof[i][dof_offset]].end(),
-	       vertexCoords);
-
-      // coords not yet in list?
-      if(it == (*dofCoords)[dof[i][dof_offset]].end()) {
-	// add new coords to list
-	(*dofCoords)[dof[i][dof_offset]].push_back(vertexCoords);
-      }
-    }
-  
-    for(i = 1; i <= dim; i++) {
-      num_dofs = admin->getNumberOfDOFs(INDEX_OF_DIM(i, dim));
-      node_offset = mesh->getNode(INDEX_OF_DIM(i, dim));
-      dof_offset  = admin->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim));
-      for (j = 0; j < mesh->getGeo(INDEX_OF_DIM(i, dim)); j++) {
-	node = node_offset + j;
-	for(k = 0; k < num_dofs; k++) {
-	  dof_index = dof_offset + k;
-	  if ((*interpPointInd)[dof[node][dof_index]] == -1) {
-	    (*interpPointInd)[dof[node][dof_index]] = -3; 
-	    // mark as interp. point 
-	    ni++;
-	  }
-	}
-      }
-    }
-
-    return 0;
-  }
-
-  int ValueWriter::writeInterpolPoints(ElInfo *el_info)
-  {
-    int     i, j, k;
-    int     node_offset, dof_offset, num_dofs, node, dof_index;
-    const DegreeOfFreedom **dof = el_info->getElement()->getDOF();
-
-    const Mesh *mesh = admin->getMesh();
-    int dim = mesh->getDim();
-
-    for(i = 1; i <= mesh->getDim(); i++) {
-      num_dofs = admin->getNumberOfDOFs(INDEX_OF_DIM(i, dim));
-      node_offset = mesh->getNode(INDEX_OF_DIM(i, dim));
-      dof_offset  = admin->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim));
-      for (j = 0; j < mesh->getGeo(INDEX_OF_DIM(i, dim)); j++) {
-	node = node_offset + j;
-	for(k = 0; k < num_dofs; k++) {
-	  dof_index = dof_offset + k;
-	  fprintf(valueFile, "%d ", (*interpPointInd)[dof[node][dof_index]]);
-	}
-      }
-    }
-
-    fprintf(valueFile, "\n");
-    return 0;
-  }
-
-  void ValueWriter::writeValuesOld(DOFVector<double> *values,
-				   char*              fn,
-				   double             time,
-				   int                level,
-				   Flag               traverseFlag)
-  {
-    FUNCNAME("ValueWriter::writeValues()");
-    TEST_EXIT(values)("no values specified\n");
-
-    ::std::string filename;
-    if(fn) {
-      filename.assign(fn);
-    } else {
-      filename.assign(values->getName());
-    }
-
-    const FiniteElemSpace *feSpace = values->getFESpace();
-    Mesh *mesh  = feSpace->getMesh();
-
-    interpPointInd = NEW DOFVector<int>(feSpace, "interpolation point indices");
-    dofCoords = NEW DOFVector< ::std::list<WorldVector<double> > >(feSpace, "dof coords");
-
-    admin = feSpace->getAdmin();
-
-    valueVec = values;
-
-    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
-
-    for(intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
-      (*intPointIt) = -1;
-    }
-
-    if (!(valueFile = fopen(filename.c_str(), "w")))
-      {
-	ERROR("could not open file %s for writing\n", filename.c_str());
-      }
-
-    ni = 0;
-    mesh->traverse(level, traverseFlag | Mesh::FILL_COORDS, markFct);
-
-    fprintf(valueFile, "mesh name: %s\n\n", 
-	    values->getFESpace()->getMesh()->getName().c_str());
-    fprintf(valueFile, "time: %e\n\n", time);
-    fprintf(valueFile, "number of values: 1\n\n");
-    fprintf(valueFile, "value description: %s\n", values->getName().c_str());
-    fprintf(valueFile, "number of interpolation points: %d\n", ni);
-    fprintf(valueFile, "type: scalar\n");
-    fprintf(valueFile, "interpolation type: lagrange\n");
-    fprintf(valueFile, "interpolation degree: %d\n",
-	    feSpace->getBasisFcts()->getDegree());
-    fprintf(valueFile, "end of description: %s\n\n", values->getName().c_str());
-
-
-    /* ----- write vertex values -----*/
-    DOFVector<double>::Iterator valueIt(valueVec, USED_DOFS);
-    DOFVector< ::std::list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS);
-
-    int i;
-
-    fprintf(valueFile, "vertex values: %s\n", values->getName().c_str());
-
-    for(intPointIt.reset(), valueIt.reset(), coordIt.reset();
-	!intPointIt.end(); 
-	++intPointIt, ++valueIt, ++coordIt) 
-      {
-	if(*intPointIt == -2) {
-	  for(i=0; i < (int) coordIt->size(); i++) {
-	    fprintf(valueFile, "%e\n", *valueIt);
-	  }
-	}
-      }
-
-    fprintf(valueFile, "\n\n");
-
-    /* ----- write interpolation values ----- */
-    fprintf(valueFile, "interpolation values: %s\n", values->getName().c_str());
-
-
-    int nd = 0;
-    for(intPointIt.reset(), valueIt.reset();
-	!intPointIt.end(); 
-	++intPointIt, ++valueIt) 
-      {
-	if(*intPointIt == -3) {
-	  fprintf(valueFile, "%e\n", *valueIt);
-	  *intPointIt = nd++;
-	}
-      }
-
-    fprintf(valueFile, "\n\n");
-
-    /* ----- write interpolation points for each simplex */
-    fprintf(valueFile, "element interpolation points: %s\n", 
-	    values->getName().c_str());
-    mesh->traverse(level, traverseFlag,
-		   writeInterpolPoints);
-    fprintf(valueFile, "\n");
-  
-    fclose(valueFile);
-    DELETE interpPointInd;
-    DELETE dofCoords;
-  }
-
   void ValueWriter::writeValues(DataCollector *dc,
 				const ::std::string filename,
 				double time,
@@ -200,7 +22,8 @@ namespace AMDiS {
 				Flag traverseFlag,
 				bool (*writeElem)(ElInfo*))
   {
-    FUNCNAME("ValueWriter::writeValues2");
+    FUNCNAME("ValueWriter::writeValues()");
+
     TEST_EXIT(dc)("no data collector\n");
 
     ::std::ofstream file;
@@ -226,16 +49,15 @@ namespace AMDiS {
     DOFVector< ::std::list<WorldVector<double> > >::Iterator 
       coordIt(dc->getDofCoords(), USED_DOFS);
 
-    int i;
 
     file << "vertex values: " << dc->getValues()->getName() << ::std::endl;
 
-    for(intPointIt.reset(), valueIt.reset(), coordIt.reset();
-	!intPointIt.end(); 
-	++intPointIt, ++valueIt, ++coordIt) 
-      {
-	if(*intPointIt == -2) {
-	  for(i=0; i < (int) coordIt->size(); i++) {
+    for (intPointIt.reset(), valueIt.reset(), coordIt.reset();
+	 !intPointIt.end(); 
+	 ++intPointIt, ++valueIt, ++coordIt) {
+
+	if (*intPointIt == -2) {
+	  for (int i = 0; i < (int) coordIt->size(); i++) {
 	    file << ::std::scientific << *valueIt << ::std::endl;
 	  }
 	}
@@ -246,11 +68,11 @@ namespace AMDiS {
     /* ----- write interpolation values ----- */
     file << "interpolation values: " << dc->getValues()->getName() << ::std::endl;
 
-    for(intPointIt.reset(), valueIt.reset();
-	!intPointIt.end(); 
-	++intPointIt, ++valueIt) 
-      {
-	if(*intPointIt >= 0) {
+    for (intPointIt.reset(), valueIt.reset();
+	 !intPointIt.end(); 
+	 ++intPointIt, ++valueIt) {
+
+	if (*intPointIt >= 0) {
 	  file << ::std::scientific << *valueIt << ::std::endl;
 	}
       }
@@ -260,12 +82,12 @@ namespace AMDiS {
     /* ----- write interpolation points for each simplex */
     file << "element interpolation points: " << dc->getValues()->getName() << ::std::endl;
 
-    ::std::list< ::std::list<int> >* interpPoints = dc->getInterpPoints();
-    ::std::list< ::std::list<int> >::iterator it1;
-    ::std::list<int>::iterator it2;
+    ::std::vector< ::std::vector<int> >* interpPoints = dc->getInterpPoints();
+    ::std::vector< ::std::vector<int> >::iterator it1;
+    ::std::vector<int>::iterator it2;
 
-    for(it1 = interpPoints->begin(); it1 != interpPoints->end(); ++it1) {
-      for(it2 = it1->begin(); it2 != it1->end(); ++it2) {
+    for (it1 = interpPoints->begin(); it1 != interpPoints->end(); ++it1) {
+      for (it2 = it1->begin(); it2 != it1->end(); ++it2) {
 	file << (*it2) << " ";
       }
       file << ::std::endl;
diff --git a/AMDiS/src/ValueWriter.h b/AMDiS/src/ValueWriter.h
index eb62fcef..db7c01a1 100644
--- a/AMDiS/src/ValueWriter.h
+++ b/AMDiS/src/ValueWriter.h
@@ -32,88 +32,68 @@
 
 namespace AMDiS {
 
-class DOFAdmin;
-class ElInfo;
-
-template<typename T> class DOFVector;
-
-// ============================================================================
-// ===== class ValueWriter ====================================================
-// ============================================================================
-
-/** \ingroup Output
- * \brief
- * ValueWriter is a static class which writes the values of a DOFVector
- * values to an ascii file named values->name.'dat'. This output is done
- * via two leaf-traversals of values->feSpace->mesh. In the first traversal
- * the values at the vertices are printed, in the second these at the
- * interpolation points of each element. For a closer disription of the
- * output format see (...link fehlt noch)
- */
-class ValueWriter
-{
-public:
-  /** \brief
-   * Writes DOFVector values to values->feSpace->mesh.
-   */
-  static void writeValues(DataCollector *dc,
-			  const ::std::string filename,
-			  double time = 0.0,
-			  int level = -1,
-			  Flag traverseFlag = Mesh::CALL_LEAF_EL,
-			  bool (*writeElem)(ElInfo*) = NULL);
-			  
-
-  // OBSOLETE
-  // TODO: Remove if ValueWriter::writeValues is stable.
-  static void writeValuesOld(DOFVector<double> *values,
-			     char*              fn = NULL,
-			     double             time = 0.0,
-			     int level = -1,
-			     Flag traverseFlag = Mesh::CALL_LEAF_EL);
-
-protected:
-  /** \brief
-   * Marks vertices and interpolation points in \ref interpPointInd
-   */
-  static int markFct(ElInfo *el_info);
-
-  /** \brief
-   * Writes indices of interpolation points of element
-   */
-  static int writeInterpolPoints(ElInfo *el_info);
-  
-protected:
-  /** \brief
-   * File to which the values should be written
-   */
-  static FILE *valueFile;
-
-  /** \brief
-   * Global interpolation point indexing
-   */
-  static DOFVector<int> *interpPointInd;
-
-  /** \brief
-   * list of coords for each dof
-   */
-  static DOFVector< ::std::list<WorldVector<double> > > *dofCoords;
-
-  /** \brief
-   * DOFAdmin of values
-   */
-  static const DOFAdmin *admin;
-
-  /** \brief
-   * Pointer to have a global access to values
-   */
-  static DOFVector<double> *valueVec;
-
-  /** \brief
-   * Total number of interpolation points.
+  class DOFAdmin;
+  class ElInfo;
+
+  template<typename T> class DOFVector;
+
+  // ============================================================================
+  // ===== class ValueWriter ====================================================
+  // ============================================================================
+
+  /** \ingroup Output
+   * \brief
+   * ValueWriter is a static class which writes the values of a DOFVector
+   * values to an ascii file named values->name.'dat'. This output is done
+   * via two leaf-traversals of values->feSpace->mesh. In the first traversal
+   * the values at the vertices are printed, in the second these at the
+   * interpolation points of each element. For a closer disription of the
+   * output format see (...link fehlt noch)
    */
-  static int ni;
-};
+  class ValueWriter
+  {
+  public:
+    /** \brief
+     * Writes DOFVector values to values->feSpace->mesh.
+     */
+    static void writeValues(DataCollector *dc,
+			    const ::std::string filename,
+			    double time = 0.0,
+			    int level = -1,
+			    Flag traverseFlag = Mesh::CALL_LEAF_EL,
+			    bool (*writeElem)(ElInfo*) = NULL);
+			   
+  protected:
+    /** \brief
+     * File to which the values should be written
+     */
+    static FILE *valueFile;
+
+    /** \brief
+     * Global interpolation point indexing
+     */
+    static DOFVector<int> *interpPointInd;
+
+    /** \brief
+     * list of coords for each dof
+     */
+    static DOFVector< ::std::list<WorldVector<double> > > *dofCoords;
+
+    /** \brief
+     * DOFAdmin of values
+     */
+    static const DOFAdmin *admin;
+
+    /** \brief
+     * Pointer to have a global access to values
+     */
+    static DOFVector<double> *valueVec;
+
+    /** \brief
+     * Total number of interpolation points.
+     */
+    static int ni;
+  };
 
 }
 
diff --git a/AMDiS/src/VertexInfo.h b/AMDiS/src/VertexInfo.h
index 45306eaf..75f1ae75 100644
--- a/AMDiS/src/VertexInfo.h
+++ b/AMDiS/src/VertexInfo.h
@@ -25,36 +25,37 @@
 #include "FixVec.h"
 
 namespace AMDiS {
+  
   /** \brief
    * Stores coordinates and output index for one vertex.
    */
   class VertexInfo 
-    {
-    public:
-  /** \brief
-   * Coordinates for this vertex.
-   */
-      WorldVector<double> coords;
-      
-      /** \brief
-       * Index for the output file.
-       */
-      int outputIndex;
-      
-      /** \brief
-       * Used to check, whether coords are already stored for a given dof.
-       */ 
-      bool operator==(const WorldVector<double>& c) {
-	return (c == coords);
-      };
-      
-      /** \brief
-       * Used to check, whether coords are already stored for a given dof.
-       */ 
-      bool operator!=(const WorldVector<double>& c) {
-	return (c != coords);
-      };
+  {
+  public:
+    /** \brief
+     * Coordinates for this vertex.
+     */
+    WorldVector<double> coords;
+    
+    /** \brief
+     * Index for the output file.
+     */
+    int outputIndex;
+    
+    /** \brief
+     * Used to check, whether coords are already stored for a given dof.
+     */ 
+    bool operator==(const WorldVector<double>& c) {
+      return (c == coords);
+    };
+    
+    /** \brief
+     * Used to check, whether coords are already stored for a given dof.
+     */ 
+    bool operator!=(const WorldVector<double>& c) {
+      return (c != coords);
     };
+  };
 }
 
 #endif
diff --git a/AMDiS/src/VtkWriter.cc b/AMDiS/src/VtkWriter.cc
index 01181788..1827ed74 100644
--- a/AMDiS/src/VtkWriter.cc
+++ b/AMDiS/src/VtkWriter.cc
@@ -4,65 +4,57 @@
 #include <cmath>
 
 #include "VtkWriter.h"
+#include "Traverse.h"
 #include "DataCollector.h"
 #include "DOFVector.h"
 #include "SurfaceRegion_ED.h"
 #include "ElementRegion_ED.h"
 
 namespace AMDiS { 
-  int VtkWriter::writeFile(::std::vector<DataCollector*> *dc,
-			   const char *name)
-  {
-    ::std::ofstream file;
-    ::std::list<ElementInfo> *elements = (*dc)[0]->getElementInfos();
 
-    DOFVector< ::std::list<VertexInfo> > *vertexInfos = (*dc)[0]->getVertexInfos();
+  int VtkWriter::writeFile(const char *name)
+  {
+    int nVertices = (*dc_)[0]->getNumberVertices();
+    int nElements = (*dc_)[0]->getNumberElements();
+    int vertices = (*dc_)[0]->getMesh()->getGeo(VERTEX);
+    int degree = (*dc_)[0]->getFeSpace()->getBasisFcts()->getDegree();    
+    int dim = (*dc_)[0]->getMesh()->getDim();
 
-    int nv = (*dc)[0]->getNumberVertices();
-    int ne = (*dc)[0]->getNumberElements();
-    int vertices = (*dc)[0]->getMesh()->getGeo(VERTEX);
+    if ((dim == 2) && (degree == 2)) {
+      nVertices += (*dc_)[0]->getNumberInterpPoints();
+      nElements *= 4;
+    }
 
+    ::std::ofstream file;
     file.open(name);
 
     file << "<?xml version=\"1.0\"?>" << ::std::endl;
     file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << ::std::endl;
     file << "  <UnstructuredGrid>" << ::std::endl;
-    file << "    <Piece NumberOfPoints=\"" << nv << "\" NumberOfCells=\"" << ne << "\">" << ::std::endl;
+    file << "    <Piece NumberOfPoints=\"" << nVertices 
+	 << "\" NumberOfCells=\"" << nElements << "\">" << ::std::endl;
     file << "      <Points>" << ::std::endl;
     file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << ::std::endl;
 
-
-    DOFVector< ::std::list<VertexInfo> >::Iterator it(vertexInfos, USED_DOFS);
-    int counter = 0;
-
-    // for all DOFs
-    for (it.reset(); !it.end(); ++it) {
-      // for all vertex infos of this DOF
-      ::std::list<VertexInfo>::iterator it2;
-      for (it2 = it->begin(); it2 != it->end(); ++it2) {
-	it2->outputIndex = counter++;
-	for (int i = 0; i < Global::getGeo(WORLD); i++) {
-	  file << " " << ::std::scientific << it2->coords[i];
-	}
-	for (int i = Global::getGeo(WORLD); i < 3; i++) {
-	  file << " 0.0";
-	}
-	file << ::std::endl;
-      }
+    if ((dim == 2) && (degree == 2)) {
+      writeVertexCoords_dim2_degree2(file);
+    } else {
+      writeVertexCoords(file);
     }
 
     file << "        </DataArray>" << ::std::endl;
     file << "      </Points>" << ::std::endl;
     file << "      <Cells>" << ::std::endl;
-
     file << "        <DataArray type=\"Int32\" Name=\"offsets\">" << ::std::endl;
-    for (int i = 0; i < ne; i++) {
+
+    for (int i = 0; i < nElements; i++) {
       file << " " << (i + 1) * vertices << ::std::endl;
     }
-    file << "        </DataArray>" << ::std::endl;
 
+    file << "        </DataArray>" << ::std::endl;
     file << "        <DataArray type=\"UInt8\" Name=\"types\">" << ::std::endl;
-    for (int i = 0; i < ne; i++) {
+
+    for (int i = 0; i < nElements; i++) {
       switch (vertices) {
       case 2:
 	file << " 3" << ::std::endl;
@@ -77,47 +69,30 @@ namespace AMDiS {
 	break;
       }	   
     }
-    file << "        </DataArray>" << ::std::endl;
 
+    file << "        </DataArray>" << ::std::endl;
     file << "        <DataArray type=\"Int32\" Name=\"connectivity\">" << ::std::endl;
-    ::std::list<ElementInfo>::iterator elementIt;
-    
-    for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
-      // for all vertices
-      for (int i = 0; i < vertices; i++) {
-	file << " " << elementIt->vertexInfo[i]->outputIndex;
-      }
-      file << ::std::endl;
+
+    if ((dim == 2) && (degree == 2)) {
+      writeConnectivity_dim2_degree2(file);
+    } else {
+      writeConnectivity(file);
     }
-    file << "        </DataArray>" << ::std::endl;
-    
+
+    file << "        </DataArray>" << ::std::endl;    
     file << "      </Cells>" << ::std::endl;
     file << "      <PointData>" << ::std::endl;
-
-    /* ----- write vertex values -----*/
-    for (int i = 0; i < static_cast<int>(dc->size()); i++) {
+    
+    for (int i = 0; i < static_cast<int>(dc_->size()); i++) {
       file << "        <DataArray type=\"Float32\" Name=\"value" << i 
 	   << "\" format=\"ascii\">" << ::std::endl;
-
-      DOFVector<int> *interpPointInd = (*dc)[i]->getInterpPointInd();
-      DOFVector<double> *values = (*dc)[i]->getValues();
-      DOFVector< ::std::list<WorldVector<double> > > *dofCoords = (*dc)[i]->getDofCoords();
-
-      DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
-      DOFVector<double>::Iterator valueIt(values, USED_DOFS);
-      DOFVector< ::std::list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS);
-      
-      for (intPointIt.reset(), valueIt.reset(), coordIt.reset();
-	   !intPointIt.end(); 
-	   ++intPointIt, ++valueIt, ++coordIt) 
-	{
-	  if (*intPointIt == -2) {
-	    for (int j = 0; j < static_cast<int>(coordIt->size()); j++) {
-		file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n";
-	    }
-	  }
-	}
       
+      if ((dim == 2) && (degree == 2)) {
+	writeVertexValues_dim2_degree2(file, i);
+      } else {
+	writeVertexValues(file, i);
+      }
+
       file << "        </DataArray>" << ::std::endl;
     }
 
@@ -131,6 +106,137 @@ namespace AMDiS {
     return 0;
   }
 
+
+  void VtkWriter::writeVertexCoords(::std::ofstream &file)
+  {
+    DOFVector< ::std::list<VertexInfo> > *vertexInfos = (*dc_)[0]->getVertexInfos();
+
+    DOFVector< ::std::list<VertexInfo> >::Iterator it(vertexInfos, USED_DOFS);
+    int counter = 0;
+
+    // for all DOFs
+    for (it.reset(); !it.end(); ++it) {
+      // for all vertex infos of this DOF
+      ::std::list<VertexInfo>::iterator it2;
+      for (it2 = it->begin(); it2 != it->end(); ++it2) {
+	it2->outputIndex = counter++;
+	writeCoord(file, it2->coords);
+      }
+    }
+  }
+
+
+  void VtkWriter::writeVertexCoords_dim2_degree2(::std::ofstream &file)
+  {
+    writeVertexCoords(file);    
+
+    DOFVector< ::std::list< WorldVector<double> > > *interpPointCoords = (*dc_)[0]->getInterpPointCoords();
+    DOFVector< ::std::list< WorldVector<double> > >::Iterator it(interpPointCoords, USED_DOFS);
+
+    for (it.reset(); !it.end(); ++it) {
+      ::std::list< WorldVector<double> >::iterator it2;
+      for (it2 = it->begin(); it2 != it->end(); ++it2) {
+	writeCoord(file, *it2);
+      }
+    }
+  }
+
+
+  void VtkWriter::writeVertexValues(::std::ofstream &file, int componentNo)
+  {    
+    DOFVector<int> *interpPointInd = (*dc_)[componentNo]->getInterpPointInd();
+    DOFVector<double> *values = (*dc_)[componentNo]->getValues();
+    DOFVector< ::std::list<WorldVector<double> > > *dofCoords = (*dc_)[componentNo]->getDofCoords();
+    
+    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
+    DOFVector<double>::Iterator valueIt(values, USED_DOFS);
+    DOFVector< ::std::list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS);
+      
+    for (intPointIt.reset(), valueIt.reset(), coordIt.reset();
+	 !intPointIt.end(); 
+	 ++intPointIt, ++valueIt, ++coordIt) {
+
+      if (*intPointIt == -2) {
+	for (int i = 0; i < static_cast<int>(coordIt->size()); i++) {
+	  file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << ::std::endl;
+	}
+      }
+    }        
+  }
+
+
+  void VtkWriter::writeVertexValues_dim2_degree2(::std::ofstream &file, int componentNo)
+  {
+    writeVertexValues(file, componentNo);
+
+    DOFVector<int>::Iterator intPointIt((*dc_)[componentNo]->getInterpPointInd(), USED_DOFS);
+    DOFVector<double>::Iterator valueIt((*dc_)[componentNo]->getValues(), USED_DOFS);
+    DOFVector< ::std::list<WorldVector<double> > >::Iterator 
+      coordIt((*dc_)[componentNo]->getInterpPointCoords(), USED_DOFS);
+
+    for (intPointIt.reset(), valueIt.reset(), coordIt.reset();
+	 !intPointIt.end(); 
+	 ++intPointIt, ++valueIt, ++coordIt) {      
+      
+      if (*intPointIt >= 0) {
+	for (int i = 0; i < static_cast<int>(coordIt->size()); i++) {
+	  file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << ::std::endl;
+	}
+      }
+    }
+  }
+
+
+  void VtkWriter::writeConnectivity(::std::ofstream &file) 
+  {
+    ::std::list<ElementInfo> *elements = (*dc_)[0]->getElementInfos();
+    ::std::list<ElementInfo>::iterator elementIt;
+    int vertices = (*dc_)[0]->getMesh()->getGeo(VERTEX);
+    
+    for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
+      // for all vertices
+      for (int i = 0; i < vertices; i++) {
+	file << " " << elementIt->vertexInfo[i]->outputIndex;
+      }
+      file << ::std::endl;
+    } 
+  }
+
+
+  void VtkWriter::writeConnectivity_dim2_degree2(::std::ofstream &file)
+  {  
+    ::std::list<ElementInfo> *elements = (*dc_)[0]->getElementInfos();
+    ::std::list<ElementInfo>::iterator elementIt;
+
+    ::std::vector< ::std::vector<int> > *interpPoints = (*dc_)[0]->getInterpPoints();
+    ::std::vector< ::std::vector<int> >::iterator pointIt;
+
+    int nVertices = (*dc_)[0]->getNumberVertices();
+
+    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
+	 pointIt != interpPoints->end(); 
+	 ++pointIt, ++elementIt) {
+    
+      file << " " << elementIt->vertexInfo[0]->outputIndex
+	   << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[2] + nVertices << ::std::endl;
+
+      file << " " << elementIt->vertexInfo[2]->outputIndex 
+	   << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[1] + nVertices << ::std::endl;
+
+      file << " " << elementIt->vertexInfo[1]->outputIndex
+	   << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[2] + nVertices << ::std::endl;
+
+      file << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[1] + nVertices 
+	   << " " << (*pointIt)[2] + nVertices << ::std::endl;
+    }
+
+  }
+
+
   int VtkWriter::updateAnimationFile(::std::string valueFilename,
 				     ::std::vector< ::std::string > *paraViewAnimationFrames,
 				     const char *animationFilename)
@@ -161,4 +267,6 @@ namespace AMDiS {
     
     return 0;
   }
+
+
 }
diff --git a/AMDiS/src/VtkWriter.h b/AMDiS/src/VtkWriter.h
index c219f86e..2ae47482 100644
--- a/AMDiS/src/VtkWriter.h
+++ b/AMDiS/src/VtkWriter.h
@@ -22,6 +22,8 @@
 #ifndef AMDIS_VTKWRITER_H
 #define AMDIS_VTKWRITER_H
 
+#include <fstream>
+
 #include "DataCollector.h"
 
 namespace AMDiS {
@@ -29,18 +31,96 @@ namespace AMDiS {
   class VtkWriter
   {
   public:
+    VtkWriter(::std::vector<DataCollector*> *dc,
+	      int level = -1,
+	      Flag traverseFlag = Mesh::CALL_LEAF_EL)
+      : dc_(dc),
+        level_(level),
+        traverseFlag_(traverseFlag)
+    {
+    };  
+
+
     /** \brief
      * Writes a ParaView-VTK file.
      */
-    static int writeFile(::std::vector<DataCollector*> *dc,
-			 const char *name);
+    int writeFile(const char *name);
+
 
     /** \brief
      * Adds a new entry to a ParaView animation file.
      */
-    static int updateAnimationFile(::std::string valueFilename,
-				   ::std::vector< ::std::string > *paraViewAnimationFrames,
-				   const char *animationFilename);
+    int updateAnimationFile(::std::string valueFilename,
+			    ::std::vector< ::std::string > *paraViewAnimationFrames,
+			    const char *animationFilename);
+  protected:
+    /** \brief
+     * 
+     */
+    void writeVertexCoords(::std::ofstream &file);
+
+
+    /** \brief
+     * 
+     */
+    void writeVertexCoords_dim2_degree2(::std::ofstream &file);
+
+
+    /** \brief
+     * 
+     */
+    void writeVertexValues(::std::ofstream &file, int componentNo);
+
+
+    /** \brief
+     * 
+     */
+    void writeVertexValues_dim2_degree2(::std::ofstream &file, int componentNo);
+
+
+    /** \brief
+     * 
+     */
+    void writeConnectivity(::std::ofstream &file);
+
+    /** \brief
+     *
+     */
+    void writeConnectivity_dim2_degree2(::std::ofstream &file);
+
+
+    /** \brief
+     *
+     */
+    inline void writeCoord(::std::ofstream &file, WorldVector<double> coord) {
+      for (int i = 0; i < Global::getGeo(WORLD); i++) {
+	file << " " << ::std::scientific << coord[i];
+      }
+      for (int i = Global::getGeo(WORLD); i < 3; i++) {
+	file << " 0.0";
+      }
+      file << ::std::endl;
+    }
+
+  private:
+    /**
+     *
+     */
+    ::std::vector<DataCollector*> *dc_;
+
+
+    /** \brief
+     * Level information for traversing the mesh.
+     */
+    int level_;
+
+
+    /** \brief
+     * Flags for traversing the mesh.
+     */
+    Flag traverseFlag_;
+
+   
   };
 }
 
-- 
GitLab