Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 018741f1886afdc512cdde2e7bb4e7e5d4f1e041)
+++ src/boundary.cpp	(revision 86234be7f2f4f7062c6ddf83c803a72d008d6945)
@@ -8,7 +8,7 @@
 
 #define DEBUG 1
-#define DoSingleStepOutput 1
+#define DoSingleStepOutput 0
 #define DoTecplotOutput 1
-#define DoRaster3DOutput 1
+#define DoRaster3DOutput 0
 #define DoVRMLOutput 1
 #define TecplotSuffix ".dat"
@@ -45,6 +45,5 @@
 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
-  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
-      << endl;
+  cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
   if (line->endpoints[0] == this)
     {
@@ -100,5 +99,11 @@
   for (int i = 0; i < 2; i++) {
     cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-    endpoints[i]->lines.erase(Numbers[i]);
+    // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
+    pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
+    for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
+      if ((*Runner).second == this) {
+        endpoints[i]->lines.erase(Runner);
+        break;
+      }
     if (endpoints[i]->lines.empty()) {
       cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
@@ -196,6 +201,6 @@
     lines[i]->triangles.erase(Nr);
     if (lines[i]->triangles.empty()) {
-      cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
       if (lines[i] != NULL) {
+        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
         delete (lines[i]);
         lines[i] = NULL;
@@ -203,5 +208,5 @@
         cerr << "ERROR: This line " << i << " has already been free'd." << endl;
     } else
-      cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
   }
 }
@@ -1471,15 +1476,11 @@
   TPS[n] = new class BoundaryPointSet(Candidate);
   InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
-  if (InsertUnique.second) // if new point was not present before, increase counter
-    {
-      PointsOnBoundaryCount++;
-    }
-  else
-    {
-      delete TPS[n];
-      cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node)
-          << " gibt's schon in der PointMap." << endl;
-      TPS[n] = (InsertUnique.first)->second;
-    }
+  if (InsertUnique.second) { // if new point was not present before, increase counter
+    PointsOnBoundaryCount++;
+  } else {
+    delete TPS[n];
+    cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    TPS[n] = (InsertUnique.first)->second;
+  }
 }
 ;
@@ -1504,6 +1505,5 @@
       if (FindLine->second->TrianglesCount < 2) {
         insertNewLine = false;
-        cout << Verbose(2)
-        << "Using existing line " << *FindLine->second << endl;
+        cout << Verbose(2) << "Using existing line " << *FindLine->second << endl;
 
         BPS[0] = FindLine->second->endpoints[0];
@@ -1532,12 +1532,11 @@
 void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
-  cout << Verbose(2)
-    << "Adding line which has not been used before between "
-    << *(a->node) << " and " << *(b->node) << "." << endl;
+  cout << Verbose(2) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
-  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-
+  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
+  // add line to global map
   LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
+  // increase counter
   LinesOnBoundaryCount++;
 };
@@ -1547,20 +1546,13 @@
  */
 void
-Tesselation::AddTriangleToLines()
-{
-
-  cout << Verbose(1) << "Adding triangle to its lines" << endl;
-  int i = 0;
+Tesselation::AddTriangle()
+{
+  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
+
+  // add triangle to global map
   TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
   TrianglesOnBoundaryCount++;
 
-  /*
-   * this is apparently done when constructing triangle
-
-   for (i=0; i<3; i++)
-   {
-   BLS[i]->AddTriangle(BTS);
-   }
-   */
+  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
 }
 ;
@@ -2868,5 +2860,5 @@
 	  // ... and triangles to the Maps of the Tesselation class
 	  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-	  AddTriangleToLines();
+	  AddTriangle();
 	  // ... and calculate its normal vector (with correct orientation)
 	  (*it)->OptCenter.Scale(-1.);
@@ -2921,5 +2913,5 @@
     }
   if (counter > 1) {
-    cout << Verbose(2) << "INFO: Degenerate triangle is ok, there is only at most " << counter << " new line needed." << endl;
+    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
     result = true;
   }
@@ -3050,5 +3042,5 @@
 
       BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangleToLines();
+      AddTriangle();
       (*it)->OptCenter.Scale(-1.);
       BTS->GetNormalVector((*it)->OptCenter);
@@ -3071,6 +3063,5 @@
 
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          //TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-          AddTriangleToLines();
+          AddTriangle();
 
           (*it)->OtherOptCenter.Scale(-1.);
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 018741f1886afdc512cdde2e7bb4e7e5d4f1e041)
+++ src/boundary.hpp	(revision 86234be7f2f4f7062c6ddf83c803a72d008d6945)
@@ -99,5 +99,5 @@
     void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
     void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AddTriangleToLines();
+    void AddTriangle();
     void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC);
     bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC);
