Changeset c15ca2

Timestamp:

Dec 29, 2009, 2:48:49 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:

97498a

Parents:

71b20e

Message:

Changes to FindClosest...ToVector() and use of ..Map,..Set,..List defines in tesselation.cpp.

• Search and Replace for all occurences of defines ..Map,..Set,..List in tesselation.cpp to tighten code
• New function Tesselation::FindClosestBoundaryPointsToVector() - replace Find(Second)ClosestPoint()
• New functions Tesselation::FindClosestBoundaryLineToVector() and Tesselation::FindClosesttriangleToVector()
  • old functions FindClosestTriangleToPoint() and FindClosestTrianglesToPoint() start looking from BoundaryPoints only and thus might find a triangle being closer which we see from the wrong side
• added unit test TesselationInOutsideUnitTest - checks IsInnerPoint and containments
• adapated unit test TesselationUnitTest - checks tesselation only

Location:

src

Files:

• analysis_correlation.cpp (modified) (2 diffs)
• tesselation.cpp (modified) (44 diffs)
• tesselation.hpp (modified) (3 diffs)
• unittests/Makefile.am (modified) (2 diffs)
• unittests/tesselation_insideoutsideunittest.cpp (added)
• unittests/tesselation_insideoutsideunittest.hpp (added)
• unittests/tesselationunittest.cpp (modified) (4 diffs)
• unittests/tesselationunittest.hpp (modified) (2 diffs)

src/analysis_correlation.cpp

@@ -267 +267 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
+          triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
           if (triangle != NULL) {
             distance = DistanceToTrianglePlane(Walker->node, triangle);
@@ -327 +327 @@
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
-                triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
+                triangle = Surface->FindClosestTriangleToVector(&checkX, LC );
                 if (triangle != NULL) {
                   distance = DistanceToTrianglePlane(&checkX, triangle);

src/tesselation.cpp

@@ -1423 +1423 @@
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1438 +1438 @@
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
@@ -2339 +2339 @@
 
   // fill the set of neighbours
-  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  Center.SubtractVector(TurningPoint->node);
-  set<TesselPoint*> SetOfNeighbours;
+  TesselPointSet SetOfNeighbours;
   SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
-  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
+  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
+  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
+    Log() << Verbose(0) << **TesselRunner << endl;
   TesselPointList::iterator Runner = connectedClosestPoints->begin();
   TesselPointList::iterator Sprinter = Runner;
@@ -2357 +2358 @@
     AddTesselationPoint((*Sprinter), 2);
 
-
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2374 +2374 @@
     Runner = Sprinter;
     Sprinter++;
+    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
+    if (Sprinter != connectedClosestPoints->end())
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
   }
   delete(connectedClosestPoints);
@@ -2958 +2961 @@
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  PointMap OrderMap;
+  PointTestPair OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -2979 +2982 @@
 };
 
-/** Finds the triangle that is closest to a given Vector \a *x.
+/** Finds the boundary points that are closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
- */
-list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
-{
-        Info FunctionInfo(__func__);
+ * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
+ */
+DistanceMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
   PointMap::const_iterator FindPoint;
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
@@ -3001 +3004 @@
   Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
 
-  set<BoundaryPointSet *> points;
+  DistanceMap * points = new DistanceMap;
   LC->GetNeighbourBounds(Nlower, Nupper);
   //Log() << Verbose(1) << endl;
@@ -3013 +3016 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end())
-              points.insert(FindPoint->second);
+              points->insert(DistancePair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
           }
         } else {
@@ -3021 +3024 @@
 
   // check whether we found some points
-  if (points.empty()) {
-    Log() << Verbose(1) << "There is no nearest triangle, too far away from the surface." << endl;
+  if (points->empty()) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    delete(points);
+    return NULL;
+  }
+  return points;
+};
+
+/** Finds the boundary line that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return closest BoundaryLineSet or NULL in degenerate case.
+ */
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+
+  // get closest points
+  DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
     return NULL;
   }
@@ -3033 +3055 @@
   Vector Center;
   Vector BaseLine;
-  for (set<BoundaryPointSet *>::iterator Runner = points.begin(); Runner != points.end(); Runner++) {
-    for (LineMap::iterator LineRunner = (*Runner)->lines.begin(); LineRunner != (*Runner)->lines.end(); LineRunner++) {
+  for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       // calculate closest point on line to desired point
       helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
@@ -3067 +3089 @@
     }
   }
+  delete(points);
   // check whether closest line is "too close" :), then it's inside
   if (ClosestLine == NULL) {
@@ -3072 +3095 @@
     return NULL;
   }
-  list<BoundaryTriangleSet*> *triangles = new list<BoundaryTriangleSet*>;
+  return ClosestLine;
+};
+
+
+/** Finds the triangle that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return BoundaryTriangleSet of nearest triangle or NULL.
+ */
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
+
+        // get closest points
+        DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    return NULL;
+  }
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
+  BoundaryLineSet *ClosestLine = NULL;
+  double MinDistance = -1.;
+  Vector helper;
+  Vector Center;
+  Vector BaseLine;
+  for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+      for (TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++) {
+        // calculate closest point on line to desired point
+        helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+        helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
+        helper.Scale(0.5);
+        Center.CopyVector(x);
+        Center.SubtractVector(&helper);
+        BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+        BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        Center.ProjectOntoPlane(&BaseLine);
+        const double distance = Center.NormSquared();
+        if ((ClosestLine == NULL) || (distance < MinDistance)) {
+          // additionally calculate intersection on line (whether it's on the line section or not)
+          helper.CopyVector(x);
+          helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+          helper.SubtractVector(&Center);
+          const double lengthA = helper.ScalarProduct(&BaseLine);
+          helper.CopyVector(x);
+          helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+          helper.SubtractVector(&Center);
+          const double lengthB = helper.ScalarProduct(&BaseLine);
+          if (lengthB*lengthA < 0) {  // if have different sign
+            ClosestLine = LineRunner->second;
+            MinDistance = distance;
+            Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
+          } else {
+            Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
+          }
+        } else {
+          Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+        }
+      }
+    }
+  }
+  delete(points);
+
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLine == NULL) {
+    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
+    return NULL;
+  }
+  TriangleList * candidates = new TriangleList;
   for (TriangleMap::iterator Runner = ClosestLine->triangles.begin(); Runner != ClosestLine->triangles.end(); Runner++) {
-    triangles->push_back(Runner->second);
-  }
-  return triangles;
+    candidates->push_back(Runner->second);
+  }
+  return candidates;
 };
 
@@ -3085 +3178 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
 {
         Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
-  list<BoundaryTriangleSet*> candidates;
+  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
+  TriangleList candidates;
   Vector Center;
   Vector helper;
@@ -3099 +3192 @@
   // reject all which are not closest
   double MinDistance = -1.;
-  for (list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
     (*Runner)->GetCenter(&Center);
     helper.CopyVector(x);
@@ -3149 +3242 @@
 {
         Info FunctionInfo(__func__);
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
+  class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
   Vector helper;
@@ -3197 +3290 @@
 }
 
-/**
- * Finds the point which is second closest to the provided one.
- *
- * @param Point to which to find the second closest other point
- * @param linked cell structure
- *
- * @return point which is second closest to the provided one
- * @TODO Optimize me! These finds are unnecessary if we have a LC list of BoundaryPoints only!
- */
-TesselPoint* Tesselation::FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const
-{
-  Info FunctionInfo(__func__);
-  TesselPoint* closestPoint = NULL;
-  TesselPoint* secondClosestPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //Log() << Verbose(1) << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        const LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
-            if (PointsOnBoundary.find((*Runner)->nr) != PointsOnBoundary.end()) {
-              if (currentNorm < distance) {
-                // remember second point
-                secondDistance = distance;
-                secondClosestPoint = closestPoint;
-                // mark down new closest point
-                distance = currentNorm;
-                closestPoint = (*Runner);
-                //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
-              }
-            }
-          }
-        } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  return secondClosestPoint;
-};
-
-/**
- * Finds the point which is closest to the provided one.
- *
- * @param Point to which to find the closest other point
- * @param SecondPoint the second closest other point on return, NULL if none found
- * @param linked cell structure
- *
- * @return point which is closest to the provided one, NULL if none found
- * @TODO Optimize me! These finds are unnecessary if we have a LC list of BoundaryPoints only!
- */
-TesselPoint* Tesselation::FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) const
-{
-  Info FunctionInfo(__func__);
-  TesselPoint* closestPoint = NULL;
-  SecondPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //Log() << Verbose(1) << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        const LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper.NormSquared();
-            if (PointsOnBoundary.find((*Runner)->nr) != PointsOnBoundary.end()) {
-              if (currentNorm < distance) {
-                secondDistance = distance;
-                SecondPoint = closestPoint;
-                distance = currentNorm;
-                closestPoint = (*Runner);
-                //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
-              } else if (currentNorm < secondDistance) {
-                secondDistance = currentNorm;
-                SecondPoint = (*Runner);
-                //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
-              }
-            }
-          }
-        } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-  // output
-  if ((closestPoint != NULL) && (SecondPoint != NULL)) {
-    Log() << Verbose(1) << "Closest point is " << *closestPoint << " and second closest is " << *SecondPoint << "." << endl;
-  } else if (closestPoint != NULL) {
-    Log() << Verbose(1) << "Closest point is " << *closestPoint << " and second closest is " << SecondPoint << "." << endl;
-  } else {
-    Log() << Verbose(1) << "No closest point has been found, only " << closestPoint << " and " << SecondPoint << "." << endl;
-  }
-  return closestPoint;
-};
-
-
 /** Gets all points connected to the provided point by triangulation lines.
  *
@@ -3327 +3296 @@
  * @return set of the all points linked to the provided one
  */
-set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
+TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+        TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -3392 +3361 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector PlaneNormal;
   Vector AngleZero;
@@ -3402 +3371 @@
   Vector helper;
   const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
@@ -3413 +3382 @@
   triangles = FindTriangles(TrianglePoints);
   if ((triangles != NULL) && (!triangles->empty())) {
-    for (list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
       PlaneNormal.AddVector(&(*Runner)->NormalVector);
   } else {
@@ -3447 +3416 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -3469 +3438 @@
  * triangle we are looking for.
  *
- * @param *out output stream for debugging
  * @param *SetOfNeighbours all points for which the angle should be calculated
  * @param *Point of which get all connected points
@@ -3475 +3443 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
   Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector center;
   Vector PlaneNormal;
@@ -3493 +3461 @@
 
   // calculate central point
-  for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+  for (TesselPointSet::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -3530 +3498 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
-    Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
-    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (!InserterTest.second) {
+      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      performCriticalExit();
+    }
   }
 
@@ -3552 +3525 @@
  * @return list of the all points linked to the provided one
  */
-list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+list< TesselPointList *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
+  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -3594 +3567 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new list<TesselPoint*>;
+        connectedPath = new TesselPointList;
         triangle = NULL;
         CurrentLine = runner->second;
@@ -3668 +3641 @@
  * @return list of the closed paths
  */
-list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
-  list<TesselPoint*> *newPath = NULL;
+list<TesselPointList *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+  list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
+  TesselPointList *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++) {
+  TesselPointList::iterator CircleRunner;
+  TesselPointList::iterator CircleStart;
+
+  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -3690 +3663 @@
 
     // go through list, look for reappearance of starting Point and create list
-    list<TesselPoint*>::iterator Marker = CircleStart;
+    TesselPointList::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
         Log() << Verbose(1) << count+1 << ". closed path consists of: ";
-        newPath = new list<TesselPoint*>;
-        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        newPath = new TesselPointList;
+        TesselPointList::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3730 +3703 @@
  * \return pointer to allocated list of triangles
  */
-set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-        Info FunctionInfo(__func__);
-  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+        Info FunctionInfo(__func__);
+        TriangleSet *connectedTriangles = new TriangleSet;
 
   if (Point == NULL) {
@@ -3782 +3755 @@
   }
 
-  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  map<class BoundaryTriangleSet *, int> Candidates;
+  TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
     }
   }
@@ -3800 +3773 @@
   count=0;
   NormalVector.Zero();
-  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
-    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
-    RemoveTesselationTriangle(Runner->first);
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
+    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->second);
     count++;
   }
   Log() << 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();
-  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+  TriangleMap::iterator NumberRunner = Candidates.begin();
+  TesselPointList::iterator StartNode, MiddleNode, EndNode;
   double angle;
   double smallestangle;
@@ -3825 +3798 @@
 
       // re-create all triangles by going through connected points list
-      list<class BoundaryLineSet *> NewLines;
+      LineList NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -3931 +3904 @@
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
-        list<class BoundaryLineSet *>::iterator Candidate;
+        LineList::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -3983 +3956 @@
  *         will usually be one, in case of degeneration, there will be two
  */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-        Info FunctionInfo(__func__);
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+        TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
@@ -4114 +4087 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedLines()
+IndexToIndex * Tesselation::FindAllDegeneratedLines()
 {
         Info FunctionInfo(__func__);
         UniqueLines AllLines;
-  map<int, int> * DegeneratedLines = new map<int, int>;
+  IndexToIndex * DegeneratedLines = new IndexToIndex;
 
   // sanity check
@@ -4139 +4112 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
@@ -4158 +4131 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedTriangles()
-{
-        Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
-  map<int, int> * DegeneratedTriangles = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
+{
+        Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
+  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -4168 +4141 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -4189 +4162 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
@@ -4203 +4176 @@
 {
         Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -4463 +4436 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
-  map <int, int> *DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4615 +4588 @@
   }
 
-  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;

src/tesselation.hpp

@@ -52 +52 @@
 // ======================================================= some template functions =========================================
 
+#define IndexToIndex map <int, int>
+
 #define PointMap map < int, class BoundaryPointSet * >
 #define PointSet set < class BoundaryPointSet * >
@@ -76 +78 @@
 #define PolygonSet set < class BoundaryPolygonSet * >
 #define PolygonList list < class BoundaryPolygonSet * >
+
+#define DistanceMap multimap <double, class BoundaryPointSet * >
+#define DistancePair pair <double, class BoundaryPointSet * >
 
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
@@ -310 +315 @@
     bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
-    TesselPoint* FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const;
-    TesselPoint* FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) const;
+    DistanceMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
+    BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const;
+    TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
+    BoundaryTriangleSet* FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const;
 
     // print for debugging

src/unittests/Makefile.am

@@ -4 +4 @@
 AM_CXXFLAGS = $(CPPUNIT_CFLAGS)
 
-TESTS = ActOnAllUnitTest AnalysisBondsUnitTests AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest BondGraphUnitTest InfoUnitTest ListOfBondsUnitTest LogUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest StackClassUnitTest VectorUnitTest
+TESTS = ActOnAllUnitTest AnalysisBondsUnitTests AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest BondGraphUnitTest InfoUnitTest ListOfBondsUnitTest LogUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest StackClassUnitTest TesselationUnitTest TesselationInOutsideUnitTest VectorUnitTest
 check_PROGRAMS = $(TESTS)
 noinst_PROGRAMS = $(TESTS)
@@ -44 +44 @@
 StackClassUnitTest_LDADD = ../libmolecuilder.a
 
+TesselationUnitTest_SOURCES = tesselationunittest.cpp tesselationunittest.hpp
+TesselationUnitTest_LDADD = ../libmolecuilder.a
+
+TesselationInOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
+TesselationInOutsideUnitTest_LDADD = ../libmolecuilder.a
+
 VectorUnitTest_SOURCES = vectorunittest.cpp vectorunittest.hpp 
 VectorUnitTest_LDADD = ../libmolecuilder.a

src/unittests/tesselationunittest.cpp

@@ -30 +30 @@
   class TesselPoint *Walker;
   Walker = new TesselPoint;
-  Walker->node = new Vector(1., 0., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(1., 0., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "1");
   Walker->nr = 1;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
-  Walker->node = new Vector(-1., 1., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(-1., 1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "2");
   Walker->nr = 2;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
-  Walker->node = new Vector(-1., -1., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(-1., -1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "3");
   Walker->nr = 3;
@@ -49 +49 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 0., 1.);
-  Walker->Name = new char[3];
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "4");
   Walker->nr = 4;
@@ -59 +59 @@
   // create tesselation
   TesselStruct = new Tesselation;
-  TesselStruct->PointsOnBoundary.clear();
-  TesselStruct->LinesOnBoundary.clear();
-  TesselStruct->TrianglesOnBoundary.clear();
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->PointsOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->LinesOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->TrianglesOnBoundary.empty() );
   TesselStruct->FindStartingTriangle(SPHERERADIUS, LinkedList);
-  bool flag = false;
 
-  LineMap::iterator baseline = TesselStruct->LinesOnBoundary.begin();
-  while (baseline != TesselStruct->LinesOnBoundary.end()) {
-    if (baseline->second->triangles.size() == 1) {
-      flag = TesselStruct->FindNextSuitableTriangle(*(baseline->second), *(((baseline->second->triangles.begin()))->second), SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+  bool OneLoopWithoutSuccessFlag = true;
+  bool TesselationFailFlag = false;
+  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
+    // 2a. fill all new OpenLines
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      baseline = Runner->second;
+      if (baseline->pointlist.empty()) {
+        T = (((baseline->BaseLine->triangles.begin()))->second);
+        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+      }
     }
-    baseline++;
-    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (flag)) {
-      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-      flag = false;
+
+    // 2b. search for smallest ShortestAngle among all candidates
+    double ShortestAngle = 4.*M_PI;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+      }
+    }
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      OneLoopWithoutSuccessFlag = false;
+    else {
+      TesselStruct->AddCandidateTriangle(*baseline);
     }
   }
@@ -84 +100 @@
   delete(TesselStruct);
   for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
-    delete[]((*Runner)->Name);
     delete((*Runner)->node);
     delete(*Runner);
   }
   Corners.clear();
-};
-
-/** UnitTest for Tesselation::IsInnerPoint()
- */
-void TesselationTest::IsInnerPointTest()
-{
-  // true inside points
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.5,0.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.5,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.5), LinkedList) );
-
-  // corners
-  for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
-    CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint((*Runner), LinkedList) );
-
-  // true outside points
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,5.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,0.,5.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(1.,1.,1.), LinkedList) );
-
-  // tricky point, there are three equally close triangles
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(5.,0.,0.), LinkedList) );
-
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
 };
 

src/unittests/tesselationunittest.hpp

@@ -21 +21 @@
 {
     CPPUNIT_TEST_SUITE( TesselationTest) ;
-    CPPUNIT_TEST ( IsInnerPointTest );
     CPPUNIT_TEST ( GetAllTrianglesTest );
     CPPUNIT_TEST ( ContainmentTest );
@@ -29 +28 @@
       void setUp();
       void tearDown();
-      void IsInnerPointTest();
       void GetAllTrianglesTest();
       void ContainmentTest();