Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 523917322b3850ebb35a13d59abc58a39a18dd09)
+++ src/boundary.hpp	(revision 262baeec28914d9a539fe907ea0b7636dcb96672)
@@ -35,5 +35,5 @@
 
 #define DEBUG 1
-#define DoSingleStepOutput 1
+#define DoSingleStepOutput 0
 #define SingleStepWidth 1
 
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 523917322b3850ebb35a13d59abc58a39a18dd09)
+++ src/tesselation.cpp	(revision 262baeec28914d9a539fe907ea0b7636dcb96672)
@@ -558,4 +558,224 @@
 //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
 //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
+  return ost;
+};
+
+// ======================================== Polygons on Boundary =================================
+
+/** Constructor for BoundaryPolygonSet.
+ */
+BoundaryPolygonSet::BoundaryPolygonSet() :
+  Nr(-1)
+{
+  Info FunctionInfo(__func__);
+};
+
+/** Destructor of BoundaryPolygonSet.
+ * Just clears endpoints.
+ * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
+ */
+BoundaryPolygonSet::~BoundaryPolygonSet()
+{
+  Info FunctionInfo(__func__);
+  endpoints.clear();
+  Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl;
+};
+
+/** Calculates the normal vector for this triangle.
+ * Is made unique by comparison with \a OtherVector to point in the other direction.
+ * \param &OtherVector direction vector to make normal vector unique.
+ * \return allocated vector in normal direction
+ */
+Vector * BoundaryPolygonSet::GetNormalVector(const Vector &OtherVector) const
+{
+  Info FunctionInfo(__func__);
+  // get normal vector
+  Vector TemporaryNormal;
+  Vector *TotalNormal = new Vector;
+  PointSet::const_iterator Runner[3];
+  for (int i=0;i<3; i++) {
+    Runner[i] = endpoints.begin();
+    for (int j = 0; j<i; j++) { // go as much further
+      Runner[i]++;
+      if (Runner[i] == endpoints.end()) {
+        eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl;
+        performCriticalExit();
+      }
+    }
+  }
+  TotalNormal->Zero();
+  int counter=0;
+  for (; Runner[2] != endpoints.end(); ) {
+    TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node);
+    for (int i=0;i<3;i++) // increase each of them
+      Runner[i]++;
+    TotalNormal->AddVector(&TemporaryNormal);
+  }
+  TotalNormal->Scale(1./(double)counter);
+
+  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+  if (TotalNormal->ScalarProduct(&OtherVector) > 0.)
+    TotalNormal->Scale(-1.);
+  Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl;
+
+  return TotalNormal;
+};
+
+/** Calculates the center point of the triangle.
+ * Is third of the sum of all endpoints.
+ * \param *center central point on return.
+ */
+void BoundaryPolygonSet::GetCenter(Vector * const center) const
+{
+  Info FunctionInfo(__func__);
+  center->Zero();
+  int counter = 0;
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+    center->AddVector((*Runner)->node->node);
+    counter++;
+  }
+  center->Scale(1./(double)counter);
+}
+
+/** Checks whether the polygons contains all three endpoints of the triangle.
+ * \param *triangle triangle to test
+ * \return true - triangle is contained polygon, false - is not
+ */
+bool BoundaryPolygonSet::ContainsBoundaryTriangle(const BoundaryTriangleSet * const triangle) const
+{
+  Info FunctionInfo(__func__);
+  return ContainsPresentTupel(triangle->endpoints, 3);
+};
+
+/** Checks whether the polygons contains both endpoints of the line.
+ * \param *line line to test
+ * \return true - line is of the triangle, false - is not
+ */
+bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+{
+  return ContainsPresentTupel(line->endpoints, 2);
+};
+
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point point to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryPolygonSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+{
+  Info FunctionInfo(__func__);
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
+    if (point == (*Runner))
+      return true;
+  return false;
+};
+
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point TesselPoint to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryPolygonSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+{
+  Info FunctionInfo(__func__);
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
+    if (point == (*Runner)->node)
+      return true;
+  return false;
+};
+
+/** Checks whether given array of \a *Points coincide with polygons's endpoints.
+ * \param **Points pointer to an array of BoundaryPointSet
+ * \param dim dimension of array
+ * \return true - set of points is contained in polygon, false - is not
+ */
+bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const
+{
+  int counter = 0;
+  for(int i=0;i<dim;i++)
+    if (ContainsBoundaryPoint(Points[i]))
+      counter++;
+
+  if (counter == dim)
+    return true;
+  else
+    return false;
+};
+
+/** Checks whether given PointList coincide with polygons's endpoints.
+ * \param &endpoints PointList
+ * \return true - set of points is contained in polygon, false - is not
+ */
+bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const
+{
+  size_t counter = 0;
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+    if (ContainsBoundaryPoint(*Runner))
+      counter++;
+  }
+
+  if (counter == endpoints.size())
+    return true;
+  else
+    return false;
+};
+
+/** Checks whether given set of \a *Points coincide with polygons's endpoints.
+ * \param *P pointer to BoundaryPolygonSet
+ * \return true - is the very triangle, false - is not
+ */
+bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
+{
+  Info FunctionInfo(__func__);
+  return ContainsPresentTupel((const PointSet)P->endpoints);
+};
+
+/** Gathers all the endpoints' triangles in a unique set.
+ * \return set of all triangles
+ */
+TriangleSet * BoundaryPolygonSet::GetAllTrianglesFromEndpoints() const
+{
+  Info FunctionInfo(__func__);
+  TriangleSet *triangles = new TriangleSet;
+
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
+    for(LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)
+      for(TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++)
+        triangles->insert(Sprinter->second);
+
+  Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl;
+  return triangles;
+};
+
+/** Fills the endpoints of this polygon from the triangles attached to \a *line.
+ * \param *line lines with triangles attached
+ * \return true - polygon contains four points, false - something went wront, polygon does not contain any points
+ */
+bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)
+{
+  for(TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {
+    for (int i=0;i<3;i++)
+      endpoints.insert((Runner->second)->endpoints[i]);
+  }
+
+  if (endpoints.size() != 4)  {// remove all points on error
+    endpoints.clear();
+    return false;
+  } else
+    return true;
+};
+
+/** output operator for BoundaryPolygonSet.
+ * \param &ost output stream
+ * \param &a boundary polygon
+ */
+ostream &operator <<(ostream &ost, const BoundaryPolygonSet &a)
+{
+  ost << "[" << a.Nr << "|";
+  for(PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {
+   ost << (*Runner)->node->Name;
+   Runner++;
+   if (Runner != a.endpoints.end())
+     ost << ",";
+  }
+  ost<< "]";
   return ost;
 };
@@ -3803,2 +4023,103 @@
     TriangleFilesWritten++;
 };
+
+/** Finds all degenerated polygons and calls ReTesselateDegeneratedPolygon()/
+ * \return number of polygons found
+ */
+int Tesselation::CorrectAllDegeneratedPolygons()
+{
+  Info FunctionInfo(__func__);
+
+  /// 1. Find all simply degenerated triangles and sort into a list with each endpoint as key
+  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+
+  /// 2. Go through all lines not contained in degenerated triangles.
+  PolygonList ListofFours;
+  BoundaryPolygonSet *Four = NULL;
+  TriangleSet *T = NULL;
+  for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++) {
+    Four = new BoundaryPolygonSet;
+
+    /// 2a. Get all four endpoints of the two connected triangles.
+    Four->FillPolygonFromTrianglesOfLine((LineRunner->second));
+
+    /// 2b. Get the triangles of all those endpoints
+    T = Four->GetAllTrianglesFromEndpoints();
+
+    /// 2c. Find all pairs of those triangles that contain the four endpoints
+    if (CountTrianglePairContainingPolygon(Four, T) > 1) {
+      /// 2d. Check the number of pairs, if greater 1, we have a degenerated polygon (mark down for later use)
+      ListofFours.push_back(Four);
+    } else
+      delete(Four);
+
+    delete(T);
+  }
+
+  /// 3. Combine all edge-connected degenerated polygons
+  PolygonList::iterator PolygonWalker;  // is the inner iterator
+  PolygonList::iterator PolygonSprinter; // is the inner advanced iterator
+  for (PolygonList::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++) {
+    PolygonWalker = PolygonRunner;
+    PolygonWalker++;
+    PolygonSprinter = PolygonWalker;
+    while (PolygonWalker != ListofFours.end()) {
+      PolygonSprinter++;
+      if (ArePolygonsEdgeConnected((*PolygonRunner), (*PolygonWalker))) { // if connected
+        CombinePolygons((*PolygonRunner), (*PolygonWalker)); // combined and ...
+        ListofFours.erase(PolygonWalker); // ... remove from list
+      }
+      PolygonWalker = PolygonSprinter;
+    }
+  }
+
+  /// 4. Go through all these degenerated polygons
+  Vector NormalVector;
+  for (PolygonList::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++) {
+    /// 4a. Gather all triangles of this polygon
+    T = Four->GetAllTrianglesFromEndpoints();
+
+    /// 4a. Get NormalVector for one side (this is "front")
+    NormalVector.CopyVector(&(*(T->begin()))->NormalVector);
+
+    /// 4b. Remove all triangles whose NormalVector is in opposite direction (i.e. "back")
+    TriangleSet::const_iterator TriangleWalker;  // is the inner iterator
+    TriangleSet::const_iterator TriangleSprinter; // is the inner advanced iterator
+    PolygonWalker = PolygonRunner;
+    PolygonWalker++;
+    PolygonSprinter = PolygonWalker;
+    while (PolygonWalker != ListofFours.end()) {
+      PolygonSprinter++;
+      if ((*TriangleWalker)->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list
+        delete(*TriangleWalker);  // remove triangle
+        T->erase(TriangleWalker);
+      }
+    }
+    /// 4c. Copy all "front" triangles but with inverse NormalVector
+    BoundaryTriangleSet *triangle = NULL;
+    TriangleWalker = T->begin();
+    TriangleSprinter = T->end();
+    TriangleSprinter--; // set onto last "front" triangle
+    while (TriangleWalker != TriangleSprinter) {  // go through all front triangles
+      triangle = new BoundaryTriangleSet(*(*TriangleWalker)); // copy triangle ...
+      triangle->NormalVector.Scale(-1.);
+      T->insert(triangle); // ... and add
+      /// 4d. Add all degenerated triangles to the list of simply degenerated triangles
+      SimplyDegeneratedTriangles->insert(pair <int, int> ((*TriangleWalker)->Nr, triangle->Nr) );
+      TriangleWalker++;
+    }
+    delete(T);  // remove the triangleset
+  }
+
+  /// 5. exit
+  PolygonList::iterator PolygonRunner;
+  while (ListofFours.empty()) {
+    PolygonRunner = ListofFours.begin();
+    delete(*PolygonRunner);
+    ListofFours.erase(PolygonRunner);
+  }
+
+  int counter = SimplyDegeneratedTriangles->size();
+  delete(SimplyDegeneratedTriangles);
+  return counter;
+};
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision 523917322b3850ebb35a13d59abc58a39a18dd09)
+++ src/tesselation.hpp	(revision 262baeec28914d9a539fe907ea0b7636dcb96672)
@@ -50,21 +50,33 @@
 
 #define PointMap map < int, class BoundaryPointSet * >
+#define PointSet set < class BoundaryPointSet * >
+#define PointList list < class BoundaryPointSet * >
 #define PointPair pair < int, class BoundaryPointSet * >
 #define PointTestPair pair < PointMap::iterator, bool >
+
 #define CandidateList list <class CandidateForTesselation *>
 #define CandidateMap map <class BoundaryLineSet *, class CandidateForTesselation *>
 
 #define LineMap multimap < int, class BoundaryLineSet * >
+#define LineSet set < class BoundaryLineSet * >
+#define LineList list < class BoundaryLineSet * >
 #define LinePair pair < int, class BoundaryLineSet * >
 #define LineTestPair pair < LineMap::iterator, bool >
 
 #define TriangleMap map < int, class BoundaryTriangleSet * >
+#define TriangleSet set < class BoundaryTriangleSet * >
+#define TriangleList list < class BoundaryTriangleSet * >
 #define TrianglePair pair < int, class BoundaryTriangleSet * >
 #define TriangleTestPair pair < TrianglePair::iterator, bool >
 
+#define PolygonMap map < int, class BoundaryPolygonSet * >
+#define PolygonSet set < class BoundaryPolygonSet * >
+#define PolygonList list < class BoundaryPolygonSet * >
+
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
 #define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
 
 #define TesselPointList list <TesselPoint *>
+#define TesselPointSet set <TesselPoint *>
 
 /********************************************** declarations *******************************/
@@ -147,4 +159,36 @@
 
 ostream & operator << (ostream &ost, const BoundaryTriangleSet &a);
+
+
+// ======================================================== class BoundaryTriangleSet =======================================
+
+/** Set of BoundaryPointSet.
+ * This is just meant as a container for a group of endpoints, extending the node, line, triangle concept. However, this has
+ * only marginally something to do with the tesselation. Hence, there is no incorporation into the bookkeeping of the Tesselation
+ * class (i.e. no allocation, no deletion).
+ * \note we assume that the set of endpoints reside (more or less) on a plane.
+ */
+class BoundaryPolygonSet {
+  public:
+    BoundaryPolygonSet();
+    ~BoundaryPolygonSet();
+
+    Vector * GetNormalVector(const Vector &NormalVector) const;
+    void GetCenter(Vector *center) const;
+    bool ContainsBoundaryLine(const BoundaryLineSet * const line) const;
+    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
+    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
+    bool ContainsBoundaryTriangle(const BoundaryTriangleSet * const point) const;
+    bool ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const;
+    bool ContainsPresentTupel(const BoundaryPolygonSet * const P) const;
+    bool ContainsPresentTupel(const PointSet &endpoints) const;
+    TriangleSet * GetAllTrianglesFromEndpoints() const;
+    bool FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line);
+
+    PointSet endpoints;
+    int Nr;
+};
+
+ostream & operator << (ostream &ost, const BoundaryPolygonSet &a);
 
 // =========================================================== class TESSELPOINT ===========================================
@@ -246,4 +290,5 @@
     void RemoveDegeneratedTriangles();
     void AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC);
+    int CorrectAllDegeneratedPolygons();
 
     set<TesselPoint*> * GetAllConnectedPoints(const TesselPoint* const Point) const;
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision 523917322b3850ebb35a13d59abc58a39a18dd09)
+++ src/tesselationhelpers.cpp	(revision 262baeec28914d9a539fe907ea0b7636dcb96672)
@@ -937,2 +937,82 @@
 }
 
+/** Counts the number of triangle pairs that contain the given polygon.
+ * \param *P polygon with endpoints to look for
+ * \param *T set of triangles to create pairs from containing \a *P
+ */
+int CountTrianglePairContainingPolygon(const BoundaryPolygonSet * const P, const TriangleSet * const T)
+{
+  Info FunctionInfo(__func__);
+  // check number of endpoints in *P
+  if (P->endpoints.size() != 4) {
+    eLog() << Verbose(1) << "CountTrianglePairContainingPolygon works only on polygons with 4 nodes!" << endl;
+    return 0;
+  }
+
+  // check number of triangles in *T
+  if (T->size() < 2) {
+    eLog() << Verbose(1) << "Not enough triangles to have pairs!" << endl;
+    return 0;
+  }
+
+  // create each pair, get the endpoints and check whether *P is contained.
+  int counter = 0;
+  PointSet Trianglenodes;
+  class BoundaryPolygonSet PairTrianglenodes;
+  for(TriangleSet::iterator Walker = T->begin(); Walker != T->end(); Walker++) {
+    for (int i=0;i<3;i++)
+      Trianglenodes.insert((*Walker)->endpoints[i]);
+
+    for(TriangleSet::iterator PairWalker = Walker; PairWalker != T->end(); PairWalker++) {
+      if (Walker != PairWalker) { // skip first
+        PairTrianglenodes.endpoints = Trianglenodes;
+        for (int i=0;i<3;i++)
+          PairTrianglenodes.endpoints.insert((*PairWalker)->endpoints[i]);
+        Log() << Verbose(1) << "Current pair of triangles: " << PairTrianglenodes << "." << endl;
+        // now check
+        if (PairTrianglenodes.ContainsPresentTupel(P)) {
+          counter++;
+          Log() << Verbose(1) << "Matches with " << *P << endl;
+        } else {
+          Log() << Verbose(1) << "No match with " << *P << endl;
+        }
+      }
+    }
+    PairTrianglenodes.endpoints.clear();
+  }
+  return counter;
+};
+
+/** Checks whether two give polygons have two or more points in common.
+ * \param *P1 first polygon
+ * \param *P2 second polygon
+ * \return true - are connected, false = are note
+ */
+bool ArePolygonsEdgeConnected(const BoundaryPolygonSet * const P1, const BoundaryPolygonSet * const P2)
+{
+  Info FunctionInfo(__func__);
+  int counter = 0;
+  for(PointSet::const_iterator Runner = P1->endpoints.begin(); Runner != P1->endpoints.end(); Runner++) {
+    if (P2->ContainsBoundaryPoint((*Runner))) {
+      counter++;
+      Log() << Verbose(1) << *(*Runner) << " of second polygon is found in the first one." << endl;
+      return true;
+    }
+  }
+  return false;
+};
+
+/** Combines second into the first and deletes the second.
+ * \param *P1 first polygon, contains all nodes on return
+ * \param *&P2 second polygon, is deleted.
+ */
+void CombinePolygons(BoundaryPolygonSet * const P1, BoundaryPolygonSet * &P2)
+{
+  Info FunctionInfo(__func__);
+  for(PointSet::iterator Runner = P1->endpoints.begin(); Runner != P1->endpoints.end(); Runner++) {
+    P1->endpoints.insert((*Runner));
+  }
+  P2->endpoints.clear();
+  delete(P2);
+};
+
Index: src/tesselationhelpers.hpp
===================================================================
--- src/tesselationhelpers.hpp	(revision 523917322b3850ebb35a13d59abc58a39a18dd09)
+++ src/tesselationhelpers.hpp	(revision 262baeec28914d9a539fe907ea0b7636dcb96672)
@@ -72,4 +72,8 @@
 bool CheckListOfBaselines(const Tesselation * const TesselStruct);
 
+int CountTrianglePairContainingPolygon(const BoundaryPolygonSet * const P, const TriangleSet * const T);
+bool ArePolygonsEdgeConnected(const BoundaryPolygonSet * const P1, const BoundaryPolygonSet * const P2);
+void CombinePolygons(BoundaryPolygonSet * const P1, BoundaryPolygonSet * &P2);
+
 
 #endif /* TESSELATIONHELPERS_HPP_ */
