Changeset 065e82 for src

Timestamp:

Aug 21, 2009, 12:19:08 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

d04966

Parents:

d4fa23

Message:

FIX: Tesselation::RemovePointFromTesselatedSurface() is working correctly now, heptan is de-tesselated without any memory leaks.

• changed some verbosity of AddTesselation...() and RemoveTesselation...() to ease bug-hunting
• new function Tesselation::GetPresentTriangle() which is similar to CheckPresenceofTriangle() but returns the found triangle.
• new function Tesselation::GetAllTriangles() returns list of all triangles to a given boundary point
• the following functions and changes are all needed to recreate triangles after removal of a boundary point in the correct manner.
  • removed GetNeighboursOnCircleOfConnectedPoints(), just use GetCircleOfConnectedPoints()
  • new function Tesselation::GetAllConnectedPoints() which returns all uniquely connected points (set not list)
  • changed function GetCircleOfConnectedPoints() which works on GetAllConnectedPoints() to bring them into order of their angles
  • new function Tesselation::GetPathsOfConnectedPoints() and Tesselation::GetClosedPathsOfConnectedPoints() which build upon one another, getting the closed paths the connected points represent by looking at the adjacency of the triangles
  • changed Tesselation::RemovePointFromTesselatedSurface(): Triangles are marked for removal, then removed in one loop. We use GetClosedPathsOfConnectedPoints() to obtain the correct sequence and then go through each of the closed paths and re-create the triangles as before

Location:

src

Files:

• tesselation.cpp (modified) (19 diffs)
• tesselation.hpp (modified) (2 diffs)

src/tesselation.cpp

@@ -1179 +1179 @@
   } else {
     delete TPS[n];
-    cout << Verbose(3) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    cout << Verbose(4) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
     TPS[n] = (InsertUnique.first)->second;
   }
@@ -1196 +1196 @@
 
   if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine;
+    LineMap::iterator FindLine = a->lines.find(b->node->nr);
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-
-    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+    cout << Verbose(5) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+
+    for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
       // If there is a line with less than two attached triangles, we don't need a new line.
       if (FindLine->second->triangles.size() < 2) {
         insertNewLine = false;
-        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
+        cout << Verbose(4) << "Using existing line " << *FindLine->second << endl;
 
         BPS[0] = FindLine->second->endpoints[0];
@@ -1231 +1232 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
-  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
+  cout << Verbose(4) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
@@ -1271 +1272 @@
           cout << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
           RemoveTesselationLine(triangle->lines[i]);
-          triangle->lines[i] = NULL;
-      } else
-        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle." << endl;
+      } else {
+        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle: ";
+        for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
+          cout << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
+        cout << endl;
+//        for (int j=0;j<2;j++) {
+//          cout << Verbose(5) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
+//          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
+//            cout << "[" << *(LineRunner->second) << "] \t";
+//          cout << endl;
+//        }
+      }
+      triangle->lines[i] = NULL;  // free'd or not: disconnect
     } else
       cerr << "ERROR: This line " << i << " has already been free'd." << endl;
@@ -1293 +1304 @@
   if (line == NULL)
     return;
-  // get other endpoint number of finding copies of same line
+  // get other endpoint number for finding copies of same line
   if (line->endpoints[1] != NULL)
     Numbers[0] = line->endpoints[1]->Nr;
@@ -1320 +1331 @@
         cout << Verbose(5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         RemoveTesselationPoint(line->endpoints[i]);
-        line->endpoints[i] = NULL;
-      } else
-        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to." << endl;
+      } else {
+        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to: ";
+        for(LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
+          cout << "[" << *(LineRunner->second) << "] \t";
+        cout << endl;
+      }
+      line->endpoints[i] = NULL;  // free'd or not: disconnect
     } else
       cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
@@ -1389 +1404 @@
           }
           // Only one of the triangle lines must be considered for the triangle count.
-          *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-          return adjacentTriangleCount;
+          //*out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //return adjacentTriangleCount;
         }
       }
@@ -1399 +1414 @@
   *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
   return adjacentTriangleCount;
+};
+
+/** Checks whether the triangle consisting of the three points is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *Candidates endpoints of the triangle candidate
+ * \return NULL - none found or pointer to triangle
+ */
+class BoundaryTriangleSet * Tesselation::GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3])
+{
+  class BoundaryTriangleSet *triangle = NULL;
+  class BoundaryPointSet *Points[3];
+
+  // builds a triangle point set (Points) of the end points
+  for (int i = 0; i < 3; i++) {
+    PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      Points[i] = FindPoint->second;
+    } else {
+      Points[i] = NULL;
+    }
+  }
+
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (Points[i] != NULL) {
+      for (int j = i; j < 3; j++) {
+        if (Points[j] != NULL) {
+          LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
+          for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
+            TriangleMap *triangles = &FindLine->second->triangles;
+            for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
+              if (FindTriangle->second->IsPresentTupel(Points)) {
+                if ((triangle == NULL) || (triangle->Nr > FindTriangle->second->Nr))
+                  triangle = FindTriangle->second;
+              }
+            }
+          }
+          // Only one of the triangle lines must be considered for the triangle count.
+          //*out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //return adjacentTriangleCount;
+        }
+      }
+    }
+  }
+
+  return triangle;
 };
 
@@ -2408 +2472 @@
     return NULL;
   }
-  list<TesselPoint*> *connectedPoints = GetCircleOfConnectedPoints(out, trianglePoints[0]);
-  list<TesselPoint*> *connectedClosestPoints = GetNeighboursOnCircleOfConnectedPoints(out, connectedPoints, trianglePoints[0], x);
-  delete(connectedPoints);
+  list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(out, trianglePoints[0], x);
   trianglePoints[1] = connectedClosestPoints->front();
   trianglePoints[2] = connectedClosestPoints->back();
@@ -2495 +2557 @@
  * @param *Point of which get all connected points
  *
- * @return list of the all points linked to the provided one
- */
-list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point)
-{
-  list<TesselPoint*> *connectedPoints = new list<TesselPoint*>;
+ * @return set of the all points linked to the provided one
+ */
+set<TesselPoint*> * Tesselation::GetAllConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -2509 +2571 @@
     ReferencePoint = PointRunner->second;
   } else {
-    *out << Verbose(2) << "getCircleOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    *out << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
     ReferencePoint = NULL;
   }
 
-  // little trick so that we look just through lines connect to the BoundaryPoint 
+  // little trick so that we look just through lines connect to the BoundaryPoint
   // OR fall-back to look through all lines if there is no such BoundaryPoint
   LineMap *Lines = &LinesOnBoundary;
@@ -2520 +2582 @@
   LineMap::iterator findLines = Lines->begin();
   while (findLines != Lines->end()) {
-    takePoint = false;
-
-    if (findLines->second->endpoints[0]->Nr == Point->nr) {
-      takePoint = true;
-      current = findLines->second->endpoints[1]->node;
-    } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
-      takePoint = true;
-      current = findLines->second->endpoints[0]->node;
-    }
-    
-    if (takePoint) {
-      *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
-      connectedPoints->push_back(current);
-    }
-
-    findLines++;
+   takePoint = false;
+
+   if (findLines->second->endpoints[0]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[1]->node;
+   } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[0]->node;
+   }
+
+   if (takePoint) {
+     *out << Verbose(5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
+     connectedPoints->insert(current);
+   }
+
+   findLines++;
   }
 
@@ -2542 +2604 @@
     return NULL;
   }
+
   return connectedPoints;
-}
-
-/** Gets the two neighbouring points with respect to a reference line to the provided point.
+};
+
+
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
  * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
  * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
@@ -2552 +2616 @@
  *
  * @param *out output stream for debugging
- * @param *connectedPoints list of connected points to the central \a *Point
  * @param *Point of which get all connected points
- * @param *Reference Vector to be checked whether it is an inner point
- *
- * @return list of the two points linked to the provided one and closest to the point to be checked,
- */
-list<TesselPoint*> * Tesselation::GetNeighboursOnCircleOfConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference)
+ * @param *Reference Reference vector for zero angle or NULL for no preference
+ * @return list of the all points linked to the provided one
+ */
+list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference)
 {
   map<double, TesselPoint*> anglesOfPoints;
-  map<double, TesselPoint*>::iterator runner;
-  ;
-  Vector center, PlaneNormal, OrthogonalVector, helper, AngleZero;
-
-  if (connectedPoints->size() == 0) { // if have not found any points
-    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
-    return NULL;
-  }
+  set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(out, Point);
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
 
   // calculate central point
-  for (list<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
+  for (set<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -2585 +2646 @@
 
   // construct one orthogonal vector
-  AngleZero.CopyVector(Reference);
+  if (Reference != NULL)
+    AngleZero.CopyVector(Reference);
+  else
+    AngleZero.CopyVector((*connectedPoints->begin())->node);
   AngleZero.SubtractVector(Point->node);
   AngleZero.ProjectOntoPlane(&PlaneNormal);
@@ -2593 +2657 @@
 
   // go through all connected points and calculate angle
-  for (list<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
+  for (set<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
-    *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
+    *out << Verbose(3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
     anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
   }
 
-  list<TesselPoint*> *result = new list<TesselPoint*>;
-  runner = anglesOfPoints.begin();
-  result->push_back(runner->second);
-  runner = anglesOfPoints.end();
-  runner--;
-  result->push_back(runner->second);
-
-  *out << Verbose(2) << "List of closest points has " << result->size() << " elements, which are "
-    << *(result->front()) << " and " << *(result->back()) << endl;
-
-  return result;
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+    connectedCircle->push_back(AngleRunner->second);
+  }
+
+  delete(connectedPoints);
+  return connectedCircle;
 }
+
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
+ *
+ * @param *out output stream for debugging
+ * @param *Point of which get all connected points
+ * @return list of the all points linked to the provided one
+ */
+list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  map<double, TesselPoint*> anglesOfPoints;
+  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+  class BoundaryPointSet *ReferencePoint = NULL;
+  class BoundaryPointSet *CurrentPoint = NULL;
+  class BoundaryTriangleSet *triangle = NULL;
+  class BoundaryLineSet *CurrentLine = NULL;
+  class BoundaryLineSet *StartLine = NULL;
+
+  // find the respective boundary point
+  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
+  if (PointRunner != PointsOnBoundary.end()) {
+    ReferencePoint = PointRunner->second;
+  } else {
+    *out << Verbose(2) << "ERROR: GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    return NULL;
+  }
+
+  map <class BoundaryLineSet *, bool> Touched;
+  map <class BoundaryLineSet *, bool>::iterator line;
+  for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++)
+    Touched.insert( pair <class BoundaryLineSet *, bool>(runner->second, false) );
+  if (!ReferencePoint->lines.empty()) {
+    for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++) {
+      line = Touched.find(runner->second);
+      if (line == Touched.end()) {
+        *out << Verbose(2) << "ERROR: I could not find " << *runner->second << " in the touched list." << endl;
+      } else if (!line->second) {
+        line->second = true;
+        connectedPath = new list<TesselPoint*>;
+        triangle = NULL;
+        CurrentLine = runner->second;
+        StartLine = CurrentLine;
+        CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        *out << Verbose(3)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;
+        do {
+          // push current one
+          *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
+          connectedPath->push_back(CurrentPoint->node);
+
+          // find next triangle
+          for (TriangleMap::iterator TriangleRunner = CurrentLine->triangles.begin(); TriangleRunner != CurrentLine->triangles.end(); TriangleRunner++) {
+            if (TriangleRunner->second != triangle) { // look for first triangle not equal to old one
+              triangle = TriangleRunner->second;
+              break;
+            }
+          }
+          // find next line
+          for (int i=0;i<3;i++) {
+            if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
+              CurrentLine = triangle->lines[i];
+
+              break;
+            }
+          }
+          line = Touched.find(CurrentLine);
+          if (line == Touched.end())
+            *out << Verbose(2) << "ERROR: I could not find " << *CurrentLine << " in the touched list." << endl;
+          else
+            line->second = true;
+          // find next point
+          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+
+        } while (CurrentLine != StartLine);
+        // last point is missing, as it's on start line
+        *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
+        connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+
+        ListOfPaths->push_back(connectedPath);
+      } else {
+        *out << Verbose(3) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl;
+      }
+    }
+  } else {
+    *out << Verbose(1) << "ERROR: There are no lines attached to " << *ReferencePoint << "." << endl;
+  }
+
+  return ListOfPaths;
+}
+
+/** Gets all closed paths on the circle of points connected to the provided point by triangulation lines, if this very point is removed.
+ * From GetPathsOfConnectedPoints() extracts all single loops of intracrossing paths in the list of closed paths.
+ * @param *out output stream for debugging
+ * @param *Point of which get all connected points
+ * @return list of the closed paths
+ */
+list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(out, Point);
+  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPoint*> *newPath = NULL;
+  int count = 0;
+
+
+  list<TesselPoint*>::iterator CircleRunner;
+  list<TesselPoint*>::iterator CircleStart;
+
+  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+    connectedPath = *ListRunner;
+
+    *out << Verbose(2) << "INFO: Current path is " << connectedPath << "." << endl;
+
+    // go through list, look for reappearance of starting Point and count
+    CircleStart = connectedPath->begin();
+
+    // go through list, look for reappearance of starting Point and create list
+    list<TesselPoint*>::iterator Marker = CircleStart;
+    for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
+      if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
+        // we have a closed circle from Marker to new Marker
+        *out << Verbose(3) << count+1 << ". closed path consists of: ";
+        newPath = new list<TesselPoint*>;
+        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        for (; CircleSprinter != CircleRunner; CircleSprinter++) {
+          newPath->push_back(*CircleSprinter);
+          *out << (**CircleSprinter) << " <-> ";
+        }
+        *out << ".." << endl;
+        count++;
+        Marker = CircleRunner;
+
+        // add to list
+        ListofClosedPaths->push_back(newPath);
+      }
+    }
+  }
+  *out << Verbose(3) << "INFO: " << count << " closed additional path(s) have been created." << endl;
+
+  // delete list of paths
+  while (!ListofPaths->empty()) {
+    connectedPath = *(ListofPaths->begin());
+    ListofPaths->remove(connectedPath);
+    delete(connectedPath);
+  }
+  delete(ListofPaths);
+
+  // exit
+  return ListofClosedPaths;
+};
+
+
+/** Gets all belonging triangles for a given BoundaryPointSet.
+ * \param *out output stream for debugging
+ * \param *Point BoundaryPoint
+ * \return pointer to allocated list of triangles
+ */
+set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(ofstream *out, class BoundaryPointSet *Point)
+{
+  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+
+  if (Point == NULL) {
+    *out << Verbose(1) << "ERROR: Point given is NULL." << endl;
+  } else {
+    // go through its lines and insert all triangles
+    for (LineMap::iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
+      for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
+      connectedTriangles->insert(TriangleRunner->second);
+    }
+  }
+
+  return connectedTriangles;
+};
+
 
 /** Removes a boundary point from the envelope while keeping it closed.
@@ -2626 +2863 @@
   Vector OldPoint, TetraederVector[3];
   double volume = 0;
-  int *numbers = NULL;
   int count = 0;
-  int i;
 
   if (point == NULL) {
@@ -2644 +2879 @@
     return 0.;
   }
-  list<TesselPoint*> *CircleofPoints = GetCircleOfConnectedPoints(out, point->node);
-
-  // remove all triangles
+
+  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(out, point->node);
+  list<TesselPoint*> *connectedPath = NULL;
+
+  // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  numbers = new int[count];
-  class BoundaryTriangleSet **Candidates = new BoundaryTriangleSet*[count];
-  i=0;
+  map<class BoundaryTriangleSet *, int> Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates[i] = triangle;
-      numbers[i++] = triangle->Nr;
-    }
-  }
-  for (int j=0;j<i;j++) {
-    RemoveTesselationTriangle(Candidates[j]);
-  }
-  delete[](Candidates);
-  *out << Verbose(1) << i << " triangles were removed." << endl;
-
-  // re-create all triangles by going through connected points list
-  list<TesselPoint*>::iterator CircleRunner = CircleofPoints->begin();
-  list<TesselPoint*>::iterator OtherCircleRunner = CircleofPoints->begin();
-  class TesselPoint *CentralNode = *CircleRunner;
-  // advance two with CircleRunner and one with OtherCircleRunner
-  CircleRunner++;
-  CircleRunner++;
-  OtherCircleRunner++;
-  i=0;
-  cout << Verbose(2) << "INFO: CentralNode is " << *CentralNode << "." << endl;
-  for (; (OtherCircleRunner != CircleofPoints->end()) && (CircleRunner != CircleofPoints->end()); (CircleRunner++), (OtherCircleRunner++)) {
-    cout << Verbose(3) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
-    cout << Verbose(3) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
-    *out << Verbose(4) << "Adding new triangle points."<< endl;
-    AddTesselationPoint(CentralNode, 0);
-    AddTesselationPoint(*OtherCircleRunner, 1);
-    AddTesselationPoint(*CircleRunner, 2);
-    *out << Verbose(4) << "Adding new triangle lines."<< endl;
-    AddTesselationLine(TPS[0], TPS[1], 0);
-    AddTesselationLine(TPS[0], TPS[2], 1);
-    AddTesselationLine(TPS[1], TPS[2], 2);
-    BTS = new class BoundaryTriangleSet(BLS, numbers[i]);
-    TrianglesOnBoundary.insert(TrianglePair(numbers[i], BTS));
-    *out << Verbose(4) << "Created triangle " << *BTS << "." << endl;
-    // calculate volume summand as a general tetraeder
-    for (int j=0;j<3;j++) {
-      TetraederVector[j].CopyVector(TPS[j]->node->node);
-      TetraederVector[j].SubtractVector(&OldPoint);
-    }
-    OldPoint.CopyVector(&TetraederVector[0]);
-    OldPoint.VectorProduct(&TetraederVector[1]);
-    volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
-    // advance number
-    i++;
-    if (i >= count)
-      *out << Verbose(2) << "WARNING: Maximum of numbers reached!" << endl;
-  }
-  *out << Verbose(1) << i << " triangles were created." << endl;
-
-  delete[](numbers);
+      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+    }
+  }
+
+  // remove all triangles
+  count=0;
+  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    *out << Verbose(3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
+    RemoveTesselationTriangle(Runner->first);
+    count++;
+  }
+  *out << Verbose(1) << count << " triangles were removed." << endl;
+
+  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
+  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
+  if (count > 2) {  // less than three triangles, then nothing will be created
+    class TesselPoint *TriangleCandidates[3];
+    count = 0;
+    for ( ; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) {  // go through all closed paths
+      if (ListAdvance != ListOfClosedPaths->end())
+        ListAdvance++;
+
+      connectedPath = *ListRunner;
+      // initialize the path to start
+      list<TesselPoint*>::iterator CircleRunner = connectedPath->begin();
+      list<TesselPoint*>::iterator OtherCircleRunner = connectedPath->begin();
+      class TesselPoint *CentralNode = *CircleRunner;
+
+      // re-create all triangles by going through connected points list
+      // advance two with CircleRunner and one with OtherCircleRunner
+      CircleRunner++;
+      CircleRunner++;
+      OtherCircleRunner++;
+      cout << Verbose(3) << "INFO: CentralNode is " << *CentralNode << "." << endl;
+      for (; (OtherCircleRunner != connectedPath->end()) && (CircleRunner != connectedPath->end()); (CircleRunner++), (OtherCircleRunner++)) {
+        cout << Verbose(4) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
+        cout << Verbose(4) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
+        *out << Verbose(3) << "INFO: Creating triangle " << CentralNode->Name << ", " << (*OtherCircleRunner)->Name << " and " << (*CircleRunner)->Name << "." << endl;
+        *out << Verbose(5) << "Adding new triangle points."<< endl;
+        TriangleCandidates[0] = CentralNode;
+        TriangleCandidates[1] = *OtherCircleRunner;
+        TriangleCandidates[2] = *CircleRunner;
+        triangle = GetPresentTriangle(out, TriangleCandidates);
+        AddTesselationPoint(CentralNode, 0);
+        AddTesselationPoint(*OtherCircleRunner, 1);
+        AddTesselationPoint(*CircleRunner, 2);
+        *out << Verbose(5) << "Adding new triangle lines."<< endl;
+        AddTesselationLine(TPS[0], TPS[1], 0);
+        AddTesselationLine(TPS[0], TPS[2], 1);
+        AddTesselationLine(TPS[1], TPS[2], 2);
+        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+        AddTesselationTriangle();
+        // calculate volume summand as a general tetraeder
+        for (int j=0;j<3;j++) {
+          TetraederVector[j].CopyVector(TPS[j]->node->node);
+          TetraederVector[j].SubtractVector(&OldPoint);
+        }
+        OldPoint.CopyVector(&TetraederVector[0]);
+        OldPoint.VectorProduct(&TetraederVector[1]);
+        volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
+        // advance number
+        NumberRunner++;
+        count++;
+        if (NumberRunner == Candidates.end())
+          *out << Verbose(2) << "WARN: Maximum of numbers reached!" << endl;
+      }
+      ListOfClosedPaths->remove(connectedPath);
+      delete(connectedPath);
+    }
+    *out << Verbose(1) << count << " triangles were created." << endl;
+  } else {
+    while (!ListOfClosedPaths->empty()) {
+      ListRunner = ListOfClosedPaths->begin();
+      connectedPath = *ListRunner;
+      ListOfClosedPaths->remove(connectedPath);
+      delete(connectedPath);
+    }
+    *out << Verbose(1) << "No need to create any triangles." << endl;
+  }
+
+  delete(ListOfClosedPaths);
 
   return volume;
 };
+
+
 
 /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
@@ -3078 +3347 @@
   }
 
-  cout << Verbose(3) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
+  cout << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
 
   return phi;

src/tesselation.hpp

@@ -209 +209 @@
     bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC);
     int CheckPresenceOfTriangle(ofstream *out, class TesselPoint *Candidates[3]);
+    class BoundaryTriangleSet * GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3]);
 
     // convex envelope
@@ -221 +222 @@
     void RemoveDegeneratedTriangles();
 
-    list<TesselPoint*> * GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point);
-    list<TesselPoint*> * GetNeighboursOnCircleOfConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference);
+    set<TesselPoint*> * GetAllConnectedPoints(ofstream *out, TesselPoint* Point);
+    set<BoundaryTriangleSet*> *GetAllTriangles(ofstream *out, class BoundaryPointSet *Point);
+    list<list<TesselPoint*> *> * GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point);
+    list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point);
+    list<TesselPoint*> * GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference = NULL);
     list<BoundaryTriangleSet*> *FindTriangles(TesselPoint* Points[3]);
     list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC);