Changeset 62bb91 for src/tesselation.cpp

Timestamp:

Aug 4, 2009, 1:48:57 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:

ef0e6d

Parents:

03e57a

Message:

Fixes for concave hull creation

• fixed BoundaryLineSet::CheckConvexityCriterion(), we now compare the two heights of the triangles, not the NormalVector which was nonsense
• new functions BoundaryLineSet::GetOtherEndpoint(), BoundaryTriangleSet::GetThirdEndpoint()
• Tesselation::InsertStraddlingPoints() receives LinkedCell pointer and uses FindClosestTrianglesToPoint(), and some verbosity added
• generally replaced [aA]tom by [pPoint]
• Tesselation::CorrectConcaveBaselines() (almost) corrected
• IsInnerPoint is overloaded function which calls FindClosestTriangleToPoint() to decide
• getAngle() receives references instead of copies with exception of OrthogonalVector which is made orthonormal to reference - center (is a FIX)

File:

• src/tesselation.cpp (modified) (33 diffs)

src/tesselation.cpp

@@ -143 +143 @@
 bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
 {
-  Vector BaseLineNormal;
-  double angle = 0;
+  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2];
   // get the two triangles
   if (TrianglesCount != 2) {
@@ -151 +150 @@
   }
   // have a normal vector on the base line pointing outwards
+  *out << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
+  BaseLineCenter.CopyVector(endpoints[0]->node->node);
+  BaseLineCenter.AddVector(endpoints[1]->node->node);
+  BaseLineCenter.Scale(1./2.);
+  BaseLine.CopyVector(endpoints[0]->node->node);
+  BaseLine.SubtractVector(endpoints[1]->node->node);
+  *out << Verbose(3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
+
   BaseLineNormal.Zero();
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
-    BaseLineNormal.AddVector(&runner->second->NormalVector);
-  BaseLineNormal.Normalize();
-  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
-    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
-
+  int i=0;
+  class BoundaryPointSet *node = NULL;
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+    *out << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
+    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
+    node = runner->second->GetThirdEndpoint(this);
+    if (node != NULL) {
+      *out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
+      helper[i].CopyVector(node->node->node);
+      helper[i].SubtractVector(&BaseLineCenter);
+      helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
+      *out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
+      i++;
+    } else {
+      *out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
+      return true;
+    }
+  }
+  *out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
+  double angle = helper[0].Angle(&helper[1]);
+  if (BaseLineNormal.ScalarProduct(&helper[1]) > 0) {
+    angle = 2.*M_PI - angle;
+  }
+  *out << Verbose(3) << "The angle is " << angle << "." << endl;
   if ((angle - M_PI) > -MYEPSILON)
     return true;
@@ -176 +200 @@
   return false;
 };
+
+/** Returns other endpoint of the line.
+ * \param *point other endpoint
+ * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
+ */
+inline class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
+{
+  if (endpoints[0] == point)
+    return endpoints[1];
+  else if (endpoints[1] == point)
+    return endpoints[0];
+  else
+    return NULL;
+};
+
 
 ostream &
@@ -360 +399 @@
             || (endpoints[2] == Points[1])
             || (endpoints[2] == Points[2])
+
           ));
 };
+
+/** Returns the endpoint which is not contained in the given \a *line.
+ * \param *line baseline defining two endpoints
+ * \return pointer third endpoint or NULL if line does not belong to triangle.
+ */
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
+{
+  // sanity check
+  if (!ContainsBoundaryLine(line))
+    return NULL;
+  for(int i=0;i<3;i++)
+    if (!line->ContainsBoundaryPoint(endpoints[i]))
+      return endpoints[i];
+  // actually, that' impossible :)
+  return NULL;
+};
+
+/** Calculates the center point of the triangle.
+ * Is third of the sum of all endpoints.
+ * \param *center central point on return.
+ */
+void BoundaryTriangleSet::GetCenter(Vector *center)
+{
+  center->Zero();
+  for(int i=0;i<3;i++)
+    center->AddVector(endpoints[i]->node->node);
+  center->Scale(1./3.);
+}
 
 ostream &
@@ -809 +877 @@
             BLS[2] = LineChecker[1]->second;
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          BTS->GetCenter(&helper);
+          helper.SubtractVector(Center);
+          helper.Scale(-1);
+          BTS->GetNormalVector(helper);
           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
           TrianglesOnBoundaryCount++;
@@ -824 +896 @@
 };
 
-/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
+/** Inserts all points outside of the tesselated surface into it by adding new triangles.
  * \param *out output stream for debugging
  * \param *cloud cluster of points
+ * \param *LC LinkedCell structure to find nearest point quickly
  * \return true - all straddling points insert, false - something went wrong
  */
-bool Tesselation::InsertStraddlingPoints(ofstream *out, PointCloud *cloud)
+bool Tesselation::InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC)
 {
   Vector Intersection;
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter(out);
+  list<BoundaryTriangleSet*> *triangles = NULL;
+
+  *out << Verbose(1) << "Begin of InsertStraddlingPoints" << endl;
 
   cloud->GoToFirst();
   while (!cloud->IsLast()) {  // we only have to go once through all points, as boundary can become only bigger
     Walker = cloud->GetPoint();
+    *out << Verbose(2) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    BTS = FindClosestTriangleToPoint(out, Walker->node);
-    if (BTS == NULL) {
-      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
-      return false;
-    }
+    triangles = FindClosestTrianglesToPoint(out, Walker->node, LC);
+    if (triangles == NULL) {
+      *out << Verbose(1) << "No triangles found, probably a tesselation point itself." << endl;
+      cloud->GoToNext();
+      continue;
+    } else {
+      BTS = triangles->front();
+    }
+    *out << Verbose(2) << "Closest triangle is " << BTS << "." << endl;
     // get the intersection point
     if (BTS->GetIntersectionInsideTriangle(out, Center, Walker->node, &Intersection)) {
+      *out << Verbose(2) << "We have an intersection at " << Intersection << "." << endl;
       // we have the intersection, check whether in- or outside of boundary
       if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) {
@@ -904 +986 @@
   // exit
   delete(Center);
+  *out << Verbose(1) << "End of InsertStraddlingPoints" << endl;
   return true;
 };
 
-/** Adds an atom to the tesselation::PointsOnBoundary list.
- * \param *Walker atom to add
+/** Adds an point to the tesselation::PointsOnBoundary list.
+ * \param *Walker point to add
  */
 void
@@ -1020 +1103 @@
 
 
-/** Checks whether the triangle consisting of the three atoms is already present.
+/** 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
@@ -1075 +1158 @@
 
 /** Finds the starting triangle for find_non_convex_border().
- * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
+ * Looks at the outermost point per axis, then Find_second_point_for_Tesselation()
  * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
  * point are called.
@@ -1089 +1172 @@
   TesselPoint* FirstPoint = NULL;
   TesselPoint* SecondPoint = NULL;
-  TesselPoint* MaxAtom[NDIM];
+  TesselPoint* MaxPoint[NDIM];
   double max_coordinate[NDIM];
   Vector Oben;
@@ -1099 +1182 @@
 
   for (i = 0; i < 3; i++) {
-    MaxAtom[i] = NULL;
+    MaxPoint[i] = NULL;
     max_coordinate[i] = -1;
   }
 
-  // 1. searching topmost atom with respect to each axis
+  // 1. searching topmost point with respect to each axis
   for (int i=0;i<NDIM;i++) { // each axis
     LC->n[i] = LC->N[i]-1; // current axis is topmost cell
@@ -1115 +1198 @@
               cout << Verbose(2) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
               max_coordinate[i] = (*Runner)->node->x[i];
-              MaxAtom[i] = (*Runner);
+              MaxPoint[i] = (*Runner);
             }
           }
@@ -1126 +1209 @@
   cout << Verbose(2) << "Found maximum coordinates: ";
   for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *MaxAtom[i] << "\t";
+    cout << i << ": " << *MaxPoint[i] << "\t";
   cout << endl;
 
@@ -1133 +1216 @@
   for (int k=0;k<NDIM;k++) {
     Oben.x[k] = 1.;
-    FirstPoint = MaxAtom[k];
+    FirstPoint = MaxPoint[k];
     cout << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
 
@@ -1238 +1321 @@
  * @param RADIUS radius of the rolling ball
  * @param N number of found triangles
- * @param *LC LinkedCell structure with neighbouring atoms
+ * @param *LC LinkedCell structure with neighbouring points
  */
 bool Tesselation::Find_next_suitable_triangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
@@ -1327 +1410 @@
 
     // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-    TesselPoint *AtomCandidates[3];
-    AtomCandidates[0] = (*it)->point;
-    AtomCandidates[1] = BaseRay->endpoints[0]->node;
-    AtomCandidates[2] = BaseRay->endpoints[1]->node;
-    int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
+    TesselPoint *PointCandidates[3];
+    PointCandidates[0] = (*it)->point;
+    PointCandidates[1] = BaseRay->endpoints[0]->node;
+    PointCandidates[2] = BaseRay->endpoints[1]->node;
+    int existentTrianglesCount = CheckPresenceOfTriangle(out, PointCandidates);
 
     BTS = NULL;
@@ -1414 +1497 @@
 bool Tesselation::CorrectConcaveBaselines(ofstream *out)
 {
-  class BoundaryTriangleSet *triangle[2];
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
@@ -1420 +1502 @@
   class BoundaryLineSet *Base = NULL;
   int OldTriangles[2], OldBaseLine;
-  int i;
+  int i,m;
+
+  *out << Verbose(1) << "Begin of CorrectConcaveBaselines" << endl;
+
   for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
     Base = baseline->second;
-
+    *out << Verbose(2) << "Current baseline is " << *Base << " ... " << endl;
     // check convexity
     if (Base->CheckConvexityCriterion(out)) { // triangles are convex
-      *out << Verbose(3) << Base << " has two convex triangles." << endl;
+      *out << Verbose(2) << "... has two convex triangles." << endl;
     } else { // not convex!
+      *out << Verbose(2) << "... has two concave triangles!" << endl;
       // get the two triangles
-      i=0;
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-        triangle[i++] = runner->second;
       // gather four endpoints and four lines
       for (int j=0;j<4;j++)
@@ -1438 +1521 @@
         OldPoints[j] = NULL;
       i=0;
-      for (int m=0;m<2;m++) { // go over both triangles
-        for (int j=0;j<3;j++) { // all of their endpoints and baselines
-          if (triangle[m]->lines[j] != Base) // pick not the central baseline
-            OldLines[i++] = triangle[m]->lines[j];
-          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
-            OldPoints[m] = triangle[m]->endpoints[j];
-        }
-      }
+      m=0;
+      *out << Verbose(3) << "The four old lines are: ";
+      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+        for (int j=0;j<3;j++) // all of their endpoints and baselines
+          if (runner->second->lines[j] != Base) { // pick not the central baseline
+            OldLines[i++] = runner->second->lines[j];
+            *out << *runner->second->lines[j] << "\t";
+          }
+      *out << endl;
+      *out << Verbose(3) << "The two old points are: ";
+      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+        for (int j=0;j<3;j++) // all of their endpoints and baselines
+          if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
+            OldPoints[m++] = runner->second->endpoints[j];
+            *out << *runner->second->endpoints[j] << "\t";
+          }
+      *out << endl;
       if (i<4) {
         *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
@@ -1461 +1553 @@
         }
 
-      // remove triangles
-      for (int j=0;j<2;j++) {
-        OldTriangles[j] = triangle[j]->Nr;
-        TrianglesOnBoundary.erase(OldTriangles[j]);
-        triangle[j] = NULL;
-      }
-
-      // remove baseline
+      // remove triangles and baseline removes itself
+      m=0;
       OldBaseLine = Base->Nr;
       LinesOnBoundary.erase(OldBaseLine);
-      Base = NULL;
+      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+        TrianglesOnBoundary.erase(OldTriangles[m++] = runner->second->Nr);
+        delete(runner->second);
+      }
 
       // construct new baseline (with same number as old one)
@@ -1481 +1570 @@
       // construct new triangles with flipped baseline
       i=-1;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
+      if (OldLines[0]->IsConnectedTo(OldLines[2]))
         i=2;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
+      if (OldLines[0]->IsConnectedTo(OldLines[3]))
         i=3;
       if (i!=-1) {
@@ -1503 +1592 @@
     }
   }
+  *out << Verbose(1) << "End of CorrectConcaveBaselines" << endl;
   return true;
 };
-
-
-/** States whether point is in- or outside of a tesselated surface.
- * \param *pointer point to be checked
- * \return true - is inside, false - is outside
- */
-bool Tesselation::IsInside(Vector *pointer)
-{
-
-  // hier kommt dann Saskias Routine hin...
-
-  return true;
-};
-
-
-/** Finds the closest triangle to a given point.
- * \param *out output stream for debugging
- * \param *x second endpoint of line
- * \return pointer triangle that is closest, NULL if none was found
- */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
-{
-  class BoundaryTriangleSet *triangle = NULL;
-
-  // hier kommt dann Saskias Routine hin...
-
-  return triangle;
-};
-
 
 /** Finds the second point of starting triangle.
@@ -1542 +1603 @@
  * \param Storage[3] array storing angles and other candidate information
  * \param RADIUS radius of virtual sphere
- * \param *LC LinkedCell structure with neighbouring atoms
+ * \param *LC LinkedCell structure with neighbouring points
  */
 void Tesselation::Find_second_point_for_Tesselation(TesselPoint* a, TesselPoint* Candidate, Vector Oben, TesselPoint*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
@@ -1552 +1613 @@
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
 
-  if (LC->SetIndexToNode(a)) {  // get cell for the starting atom
+  if (LC->SetIndexToNode(a)) {  // get cell for the starting point
     for(int i=0;i<NDIM;i++) // store indices of this cell
       N[i] = LC->n[i];
   } else {
-    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
+    cerr << "ERROR: Point " << *a << " is not found in cell " << LC->index << "." << endl;
     return;
   }
-  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+  // then go through the current and all neighbouring cells and check the contained points for possible candidates
   cout << Verbose(3) << "LC Intervals from [";
   for (int i=0;i<NDIM;i++) {
@@ -1657 +1718 @@
  * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
  * @param BaseLine BoundaryLineSet with the current base line
- * @param ThirdNode third atom to avoid in search
+ * @param ThirdNode third point to avoid in search
  * @param candidates list of equally good candidates to return
  * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
  * @param RADIUS radius of sphere
- * @param *LC LinkedCell structure with neighbouring atoms
+ * @param *LC LinkedCell structure with neighbouring points
  */
 void Tesselation::Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint  *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
@@ -1714 +1775 @@
       }
 
-      // get cell for the starting atom
+      // get cell for the starting point
       if (LC->SetIndexToVector(&CircleCenter)) {
         for(int i=0;i<NDIM;i++) // store indices of this cell
@@ -1723 +1784 @@
         return;
       }
-      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+      // then go through the current and all neighbouring cells and check the contained points for possible candidates
       //cout << Verbose(2) << "LC Intervals:";
       for (int i=0;i<NDIM;i++) {
@@ -1874 +1935 @@
 };
 
+/** 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 list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
+ */
+list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC)
+{
+  class TesselPoint *trianglePoints[3];
+
+  if (LinesOnBoundary.empty()) {
+    *out << Verbose(0) << "Error: There is no tesselation structure to compare the point with, " << "please create one first.";
+    return NULL;
+  }
+
+  trianglePoints[0] = findClosestPoint(x, LC);
+  // check whether closest point is "too close" :), then it's inside
+  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
+    *out << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
+    return NULL;
+  }
+  list<TesselPoint*> *connectedClosestPoints = getCircleOfConnectedPoints(out, trianglePoints[0], x);
+  trianglePoints[1] = connectedClosestPoints->front();
+  trianglePoints[2] = connectedClosestPoints->back();
+  for (int i=0;i<3;i++) {
+    if (trianglePoints[i] == NULL) {
+      *out << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+    }
+    *out << Verbose(1) << "List of possible points:" << endl;
+    *out << *trianglePoints[i] << endl;
+  }
+  delete(connectedClosestPoints);
+
+  list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);
+
+  if (triangles->empty()) {
+    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
+    return NULL;
+  } else
+    return triangles;
+};
+
+/** Finds closest triangle to a point.
+ * This basically just takes care of the degenerate case, which is not handled in FindClosestTrianglesToPoint().
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return list of BoundaryTriangleSet of nearest triangles or NULL.
+ */
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x, LinkedCell* LC)
+{
+  class BoundaryTriangleSet *result = NULL;
+  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(out, x, LC);
+
+  if (triangles == NULL)
+    return NULL;
+
+  if (x->ScalarProduct(&triangles->front()->NormalVector) < 0)
+    result = triangles->back();
+  else
+    result = triangles->front();
+
+  delete(triangles);
+  return result;
+};
+
+/** Checks whether the provided Vector is within the tesselation structure.
+ *
+ * @param point of which to check the position
+ * @param *LC LinkedCell structure
+ *
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(ofstream *out, Vector Point, LinkedCell* LC)
+{
+  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
+  if (result == NULL)
+    return true;
+  if (Point.ScalarProduct(&result->NormalVector) < 0)
+    return true;
+  else
+    return false;
+}
+
+/** Checks whether the provided TesselPoint is within the tesselation structure.
+ *
+ * @param *Point of which to check the position
+ * @param *LC Linked Cell structure
+ *
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
+{
+  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, Point->node, LC);
+  if (result == NULL)
+    return true;
+  if (Point->node->ScalarProduct(&result->NormalVector) < 0)
+    return true;
+  else
+    return false;
+}
+
+/** Gets all points connected to the provided point by triangulation lines.
+ * 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
+ * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
+ * triangle we are looking for.
+ *
+ * @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::getCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector* Reference)
+{
+  list<TesselPoint*> connectedPoints;
+  map<double, TesselPoint*> anglesOfPoints;
+  map<double, TesselPoint*>::iterator runner;
+  list<TesselPoint*>::iterator listRunner;
+  Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
+  TesselPoint* current;
+  bool takePoint = false;
+
+  planeNorm.CopyVector(Point->node);
+  planeNorm.SubtractVector(Reference);
+  planeNorm.Normalize();
+
+  LineMap::iterator findLines = LinesOnBoundary.begin();
+  while (findLines != LinesOnBoundary.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) {
+      connectedPoints.push_back(current);
+      currentPoint.CopyVector(current->node);
+      currentPoint.ProjectOntoPlane(&planeNorm);
+      center.AddVector(&currentPoint);
+    }
+
+    findLines++;
+  }
+
+  *out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
+    << "; scale factor " << 1.0/connectedPoints.size();
+
+  center.Scale(1.0/connectedPoints.size());
+  listRunner = connectedPoints.begin();
+
+  *out << " calculated center " << center <<  endl;
+
+  // construct one orthogonal vector
+  helper.CopyVector(Reference);
+  helper.ProjectOntoPlane(&planeNorm);
+  OrthogonalVector.MakeNormalVector(&center, &helper, (*listRunner)->node);
+  while (listRunner != connectedPoints.end()) {
+    double angle = getAngle(*((*listRunner)->node), *(Reference), center, OrthogonalVector);
+    *out << "Calculated angle " << angle << " for point " << **listRunner << endl;
+    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    listRunner++;
+  }
+
+  list<TesselPoint*> *result = new list<TesselPoint*>;
+  runner = anglesOfPoints.begin();
+  *out << "First value is " << *runner->second << endl;
+  result->push_back(runner->second);
+  runner = anglesOfPoints.end();
+  runner--;
+  *out << "Second value is " << *runner->second << endl;
+  result->push_back(runner->second);
+
+  *out << "List of closest points has " << result->size() << " elements, which are "
+    << *(result->front()) << " and " << *(result->back()) << endl;
+
+  return result;
+}
+
 /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
  * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
@@ -1957 +2199 @@
 
 
+
 /**
- * Checks whether the provided atom is within the tesselation structure.
+ * Finds the point which is closest to the provided one.
  *
- * @param *Atom of which to check the position
- * @param tesselation structure
- *
- * @return true if the atom is inside the tesselation structure, false otherwise
- */
-bool IsInnerPoint(TesselPoint *Atom, class Tesselation *Tess, LinkedCell* LC) 
-{
-  if (Tess->LinesOnBoundary.begin() == Tess->LinesOnBoundary.end()) {
-    cout << Verbose(0) << "Error: There is no tesselation structure to compare the point with, "
-        << "please create one first.";
-    exit(1);
-  }
-
-  class TesselPoint *trianglePoints[3];
-  trianglePoints[0] = findClosestAtom(Atom, LC);
-  // check whether closest atom is "too close" :), then it's inside
-  if (trianglePoints[0]->node->DistanceSquared(Atom->node) < MYEPSILON)
-    return true;
-  list<TesselPoint*> *connectedClosestPoints = Tess->getClosestConnectedAtoms(trianglePoints[0], Atom);
-  trianglePoints[1] = connectedClosestPoints->front();
-  trianglePoints[2] = connectedClosestPoints->back();
-  for (int i=0;i<3;i++) {
-    if (trianglePoints[i] == NULL) {
-      cout << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-    }
-
-    cout << Verbose(1) << "List of possible atoms:" << endl;
-    cout << *trianglePoints[i] << endl;
-
-//    for(list<TesselPoint*>::iterator runner = connectedClosestPoints->begin(); runner != connectedClosestPoints->end(); runner++)
-//      cout << Verbose(2) << *(*runner) << endl;
-  }
-  delete(connectedClosestPoints);
-
-  list<BoundaryTriangleSet*> *triangles = Tess->FindTriangles(trianglePoints);
-
-  if (triangles->empty()) {
-    cout << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
-    exit(1);
-  }
-
-  Vector helper;
-  helper.CopyVector(Atom->node);
-
-  // Only in case of degeneration, there will be two different scalar products.
-  // If at least one scalar product is positive, the point is considered to be outside.
-  bool status = (helper.ScalarProduct(&triangles->front()->NormalVector) < 0)
-      && (helper.ScalarProduct(&triangles->back()->NormalVector) < 0);
-  delete(triangles);
-  return status;
-}
-
-/**
- * Finds the atom which is closest to the provided one.
- *
- * @param atom to which to find the closest other atom
+ * @param Point to which to find the closest other point
  * @param linked cell structure
  *
- * @return atom which is closest to the provided one
- */
-TesselPoint* findClosestAtom(const TesselPoint* Atom, LinkedCell* LC) 
+ * @return point which is closest to the provided one
+ */
+TesselPoint* findClosestPoint(const Vector* Point, LinkedCell* LC)
 {
   LinkedNodes *List = NULL;
-  TesselPoint* closestAtom = NULL;
+  TesselPoint* closestPoint = NULL;
   double distance = 1e16;
   Vector helper;
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
 
-  LC->SetIndexToVector(Atom->node); // ignore status as we calculate bounds below sensibly
+  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];
@@ -2042 +2230 @@
         if (List != NULL) {
           for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Atom->node);
+            helper.CopyVector(Point);
             helper.SubtractVector((*Runner)->node);
             double currentNorm = helper. Norm();
             if (currentNorm < distance) {
               distance = currentNorm;
-              closestAtom = (*Runner);
+              closestPoint = (*Runner);
             }
           }
@@ -2056 +2244 @@
       }
 
-  return closestAtom;
-}
+  return closestPoint;
+}
+
 
 /**
- * Gets all atoms connected to the provided atom by triangulation lines.
+ * Finds triangles belonging to the three provided points.
  *
- * @param atom of which get all connected atoms
- * @param atom to be checked whether it is an inner atom
+ * @param *Points[3] list, is expected to contain three points
  *
- * @return list of the two atoms linked to the provided one and closest to the atom to be checked,
- */
-list<TesselPoint*> * Tesselation::getClosestConnectedAtoms(TesselPoint* Atom, TesselPoint* AtomToCheck) 
-{
-  list<TesselPoint*> connectedAtoms;
-  map<double, TesselPoint*> anglesOfAtoms;
-  map<double, TesselPoint*>::iterator runner;
-  LineMap::iterator findLines = LinesOnBoundary.begin();
-  list<TesselPoint*>::iterator listRunner;
-  Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
-  TesselPoint* current;
-  bool takeAtom = false;
-
-  planeNorm.CopyVector(Atom->node);
-  planeNorm.SubtractVector(AtomToCheck->node);
-  planeNorm.Normalize();
-
-  while (findLines != LinesOnBoundary.end()) {
-    takeAtom = false;
-
-    if (findLines->second->endpoints[0]->Nr == Atom->nr) {
-      takeAtom = true;
-      current = findLines->second->endpoints[1]->node;
-    } else if (findLines->second->endpoints[1]->Nr == Atom->nr) {
-      takeAtom = true;
-      current = findLines->second->endpoints[0]->node;
-    }
-
-    if (takeAtom) {
-      connectedAtoms.push_back(current);
-      currentPoint.CopyVector(current->node);
-      currentPoint.ProjectOntoPlane(&planeNorm);
-      center.AddVector(&currentPoint);
-    }
-
-    findLines++;
-  }
-
-  cout << "Summed vectors " << center << "; number of atoms " << connectedAtoms.size()
-    << "; scale factor " << 1.0/connectedAtoms.size();
-
-  center.Scale(1.0/connectedAtoms.size());
-  listRunner = connectedAtoms.begin();
-
-  cout << " calculated center " << center <<  endl;
-
-  // construct one orthogonal vector
-  helper.CopyVector(AtomToCheck->node);
-  helper.ProjectOntoPlane(&planeNorm);
-  OrthogonalVector.MakeNormalVector(&center, &helper, (*listRunner)->node);
-  while (listRunner != connectedAtoms.end()) {
-    double angle = getAngle(*((*listRunner)->node), *(AtomToCheck->node), center, OrthogonalVector);
-    cout << "Calculated angle " << angle << " for atom " << **listRunner << endl;
-    anglesOfAtoms.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
-    listRunner++;
-  }
-
-  list<TesselPoint*> *result = new list<TesselPoint*>;
-  runner = anglesOfAtoms.begin();
-  cout << "First value is " << *runner->second << endl;
-  result->push_back(runner->second);
-  runner = anglesOfAtoms.end();
-  runner--;
-  cout << "Second value is " << *runner->second << endl;
-  result->push_back(runner->second);
-
-  cout << "List of closest atoms has " << result->size() << " elements, which are "
-    << *(result->front()) << " and " << *(result->back()) << endl;
-
-  return result;
-}
-
-/**
- * Finds triangles belonging to the three provided atoms.
- *
- * @param atom list, is expected to contain three atoms
- *
- * @return triangles which belong to the provided atoms, will be empty if there are none,
+ * @return triangles which belong to the provided points, will be empty if there are none,
  *         will usually be one, in case of degeneration, there will be two
  */
@@ -2202 +2313 @@
 }
 
-/**
- * Gets the angle between a point and a reference relative to the provided center.
- *
+/** Gets the angle between a point and a reference relative to the provided center.
+ * We have two shanks (point, center) and (reference, center) between which the angle is calculated
+ * and by scalar product with OrthogonalVector we decide the interval.
  * @param point to calculate the angle for
  * @param reference to which to calculate the angle
  * @param center for which to calculate the angle between the vectors
- * @param OrthogonalVector helps discern between [0,pi] and [pi,2pi] in acos()
+ * @param OrthogonalVector points in direction of [pi,2pi] interval
  *
  * @return angle between point and reference
  */
-double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector) 
+double getAngle(const Vector &point, const Vector &reference, const Vector &center, Vector OrthogonalVector)
 {
   Vector ReferenceVector, helper;
-  cout << Verbose(2) << reference << " is our reference " << endl;
-  cout << Verbose(2) << center << " is our center " << endl;
+  cout << Verbose(4) << center << " is our center " << endl;
   // create baseline vector
   ReferenceVector.CopyVector(&reference);
   ReferenceVector.SubtractVector(&center);
+  OrthogonalVector.MakeNormalVector(&ReferenceVector);
+  cout << Verbose(4) << ReferenceVector << " is our reference " << endl;
   if (ReferenceVector.IsNull())
     return M_PI;
@@ -2233 +2345 @@
   }
 
-  cout << Verbose(2) << point << " has angle " << phi << endl;
+  cout << Verbose(3) << helper << " has angle " << phi << " with respect to reference." << endl;
 
   return phi;
 }
 
-/**
- * Checks whether the provided point is within the tesselation structure.
- *
- * This is a wrapper function for IsInnerAtom, so it can be used with a simple
- * vector, too.
- *
- * @param point of which to check the position
- * @param tesselation structure
- *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC) 
-{
-  TesselPoint *temp_atom = new TesselPoint;
-  bool status = false;
-  temp_atom->node->CopyVector(&Point);
-
-  status = IsInnerPoint(temp_atom, Tess, LC);
-  delete(temp_atom);
-
-  return status;
-}
-