Changeset f07f86d

Timestamp:

Apr 20, 2010, 10:34:43 AM (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:

b60a29

Parents:

061b06

git-author:

Frederik Heber <heber@…> (04/20/10 09:26:30)

git-committer:

Frederik Heber <heber@…> (04/20/10 10:34:43)

Message:

Rewrite of tesselation procedure: treatment of degenerated triangles.

The central idea is that when matching open lines are present, then we have to check whether candidate matches desired ThirdPoint and that the calculated OptCenter matches as well. The latter is essential for degenerated triangles, that can only be decerned by their (Other)OptCenter.

• BoundaryTriangleSetCheckDegeneracy() moved back to Tesselation and gets points from TPS.-
• Tesselation:.AddCandidateTriangle()
  • candidates are now in TPS
  • sole parameter is CandidateLine which contains other data to be stored in the new triangle.
  • which is done here not in AddCandidatePolygon() anymore
• Tesselation::AddTesselationLine()
  • new parameter OptCenter to pick the right OpenLine if multiples present
    • functions adapted: Tesselation::RemovePointFromTesselatedSurface(), Tesselation::AddBoundaryPointByDegeneratedTriangle(), Tesselation::CorrectAllDegeneratedPolygons()
  • when an open line is present, we check whether the candidate matches ThirdPoint and OptCenter's match as well
• FIX: CandidateForTesselation::CheckValidity()
  • endpoint and candidate inside check have both been wrong way round (< instead of >)
  • MYEPSILON replaced by HULLEPSILON (there are case with 1.16e-13)
• Tesselation::FindStartingTriangle()
  • if one axis fail for obtaining a starting triangle, the other two are tried as well
  • For this BaseLine is now a pointer, allocated newly for each axis
• Tesselation::AddCandidatePolygon()
  • everything is done in either AddCandidateTriangle() or AddDegeneratedTriangle()
  • only setting TPS to the current three endpoints and calling CheckDegeneracy() remains
• Tesselation::AddDegeneratedTriangle()
  • parameter triangle replaced by CandidateLine
  • no more LineCase et al, but subsequent creation of two triangles and everything done in AddTesselationLine()
• FIX: Tesselation::FindThirdPointForTesselation() - we check whether we actually have candidates before CheckValidity() is called.

Asparagine is now tesselated correctly from radius 1.5 up to 25.

Location:

src

Files:

• tesselation.cpp (modified) (22 diffs)
• tesselation.hpp (modified) (3 diffs)

src/tesselation.cpp

@@ -654 +654 @@
 };
 
-/** Checks validity of a given sphere of a candidate line.
- * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
- * \param *OtherOptCenter center of the other triangle
- * \param RADIUS radius of sphere
- * \param *LC LinkedCell structure with other atoms
- * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
- */
-bool BoundaryTriangleSet::CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const
-{
-  Info FunctionInfo(__func__);
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
-  bool flag = true;
-
-  // get all points inside the sphere
-  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, OtherOptCenter);
-
-  // remove triangles's endpoints
-  for (int i=0;i<3;i++)
-    ListofPoints->remove(endpoints[i]->node);
-
-  // check for other points
-  if (!ListofPoints->empty()) {
-    Log() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl;
-    flag = false;
-    Log() << Verbose(1) << "External atoms inside of sphere at " << *OtherOptCenter << ":" << endl;
-    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
-  }
-  delete(ListofPoints);
-
-  return flag;
-};
-
-
-
 /** Returns the endpoint which is not contained in the given \a *line.
  * \param *line baseline defining two endpoints
@@ -1078 +1042 @@
   Info FunctionInfo(__func__);
 
-  const double radiusSquared = RADIUS;
+  const double radiusSquared = RADIUS*RADIUS;
   list<const Vector *> VectorList;
   VectorList.push_back(&OptCenter);
@@ -1089 +1053 @@
   // check baseline for OptCenter and OtherOptCenter being on sphere's surface
   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-    for (int i=0;i<2;i++)
-      if (fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared) < MYEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << "." << endl);
+    for (int i=0;i<2;i++) {
+      const double distance = fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared);
+      if (distance > HULLEPSILON) {
+        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << setprecision(13) << distance << "." << endl);
         return false;
       }
+    }
   }
 
@@ -1100 +1066 @@
     const TesselPoint *Walker = *Runner;
     for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-      if (fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared) < MYEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << Walker << " is out of sphere at " << *(*VRunner) << "." << endl);
+      const double distance = fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared);
+      if (distance > HULLEPSILON) {
+        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << Walker << " is out of sphere at " << *(*VRunner) << " by " << setprecision(13) << distance << "." << endl);
         return false;
       }
@@ -1827 +1794 @@
  * If successful it raises the line count and inserts the new line into the BLS,
  * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
+ * @param *OptCenter desired OptCenter if there are more than one candidate line
  * @param *candidate third point of the triangle to be, for checking between multiple open line candidates
  * @param *a first endpoint
@@ -1832 +1800 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddTesselationLine(const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
+void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
   bool insertNewLine = true;
 
@@ -1847 +1815 @@
       // If there is a line with less than two attached triangles, we don't need a new line.
       if (FindLine->second->triangles.size() == 1) {
+        CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
+        if (!Finder->second->pointlist.empty())
+          Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl;
+        else
+          Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl;
         // get open line
-        CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
-        insertNewLine = false;
-        WinningLine = FindLine->second;
-        if ((!Finder->second->pointlist.empty()) && (*(Finder->second->pointlist.begin()) == candidate->node)) { // stop searching if candidate matches
+        if ((!Finder->second->pointlist.empty()) && (*(Finder->second->pointlist.begin()) == candidate->node)
+            && (OptCenter == NULL || *OptCenter == Finder->second->OptCenter)) { // stop searching if candidate matches
+          insertNewLine = false;
+          WinningLine = FindLine->second;
           break;
         }
@@ -2063 +2036 @@
 };
 
+
+/** Checks validity of a given sphere of a candidate line.
+ * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
+ * Note that endpoints are stored in Tesselation::TPS.
+ * \param *OtherOptCenter center of the other triangle
+ * \param RADIUS radius of sphere
+ * \param *LC LinkedCell structure with other atoms
+ * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
+ */
+bool Tesselation::CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const
+{
+  Info FunctionInfo(__func__);
+
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
+  bool flag = true;
+
+  cout << Verbose(1) << "Check by: draw sphere {" << OtherOptCenter->x[0] << " " << OtherOptCenter->x[1] << " " << OtherOptCenter->x[2] << "} radius " << RADIUS << " resolution 30" << endl;
+
+  // get all points inside the sphere
+  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, OtherOptCenter);
+
+  Log() << Verbose(1) << "The following atoms are inside sphere at " << *OtherOptCenter << ":" << endl;
+  for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
+    Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
+
+  // remove triangles's endpoints
+  for (int i=0;i<3;i++)
+    ListofPoints->remove(TPS[i]->node);
+
+  // check for other points
+  if (!ListofPoints->empty()) {
+    Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl;
+    flag = false;
+    Log() << Verbose(1) << "External atoms inside of sphere at " << *OtherOptCenter << ":" << endl;
+    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
+      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
+  }
+  delete(ListofPoints);
+
+  return flag;
+};
+
+
 /** Checks whether the triangle consisting of the three points is already present.
  * Searches for the points in Tesselation::PointsOnBoundary and checks their
@@ -2182 +2198 @@
   TesselPoint* Temporary;
   double maxCoordinate[NDIM];
-  BoundaryLineSet BaseLine;
+  BoundaryLineSet *BaseLine = NULL;
   Vector helper;
   Vector Chord;
@@ -2228 +2244 @@
     NormalVector.Zero();
     NormalVector.x[k] = 1.;
-    BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
-    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
+    BaseLine = new BoundaryLineSet();
+    BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
+    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl;
 
     double ShortestAngle;
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(BaseLine->endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     if (Temporary == NULL)  // have we found a second point?
       continue;
-    BaseLine.endpoints[1] = new BoundaryPointSet(Temporary);
+    BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);
 
     // construct center of circle
-    CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);
-    CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);
+    CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node);
+    CircleCenter.AddVector(BaseLine->endpoints[1]->node->node);
     CircleCenter.Scale(0.5);
 
     // construct normal vector of circle
-    CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node);
-    CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);
+    CirclePlaneNormal.CopyVector(BaseLine->endpoints[0]->node->node);
+    CirclePlaneNormal.SubtractVector(BaseLine->endpoints[1]->node->node);
 
     double radius = CirclePlaneNormal.NormSquared();
@@ -2264 +2281 @@
 
     // adding point 1 and point 2 and add the line between them
-    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
-    Log() << Verbose(0) << "Found second point is at " << *BaseLine.endpoints[1]->node << ".\n";
+    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl;
+    Log() << Verbose(0) << "Found second point is at " << *BaseLine->endpoints[1]->node << ".\n";
 
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    CandidateForTesselation OptCandidates(&BaseLine);
+    CandidateForTesselation OptCandidates(BaseLine);
     FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
     Log() << Verbose(0) << "List of third Points is:" << endl;
@@ -2274 +2291 @@
         Log() << Verbose(0) << " " << *(*it) << endl;
     }
-
-    BTS = NULL;
-    AddCandidatePolygon(OptCandidates, RADIUS, LC);
-//    delete(BaseLine.endpoints[0]);
-//    delete(BaseLine.endpoints[1]);
+    if (!OptCandidates.pointlist.empty()) {
+      BTS = NULL;
+      AddCandidatePolygon(OptCandidates, RADIUS, LC);
+    } else {
+      delete BaseLine;
+      continue;
+    }
 
     if (BTS != NULL) // we have created one starting triangle
@@ -2286 +2305 @@
       OptCandidates.pointlist.clear();
     }
+    delete BaseLine;
   }
 };
@@ -2436 +2456 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
 {
         Info FunctionInfo(__func__);
@@ -2627 +2647 @@
   Sprinter++;
   while(Sprinter != connectedClosestPoints->end()) {
-    // add this triangle
-    AddCandidateTriangle(TurningPoint, *Runner, *Sprinter);
-
-    // store values from tesselation procedure therein
-    BTS->SetTopNode(CandidateLine.T);
-    if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
-      BTS->AngleFromTop = CandidateLine.ShortestAngle;
-    else
-      BTS->AngleFromTop = -1.;
-
-    // create normal vector
-    BTS->GetCenter(&Center);
-    Center.SubtractVector(&CandidateLine.OptCenter);
-    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
-    BTS->GetNormalVector(Center);
-
-    // give some verbose output about the whole procedure
-    if (CandidateLine.T != NULL)
-      Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
-    else
-      Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
-
-    // check whether its degenerate and if so, add the other side triangle as well
-    if (BTS->CheckDegeneracy(&CandidateLine.OtherOptCenter, RADIUS, LC)) {
-      Log() << Verbose(1) << "Triangle " << *BTS << " is degenerated, adding other side." << endl;
-      BoundaryTriangleSet * const triangle = BTS;
-      AddDegeneratedTriangle(BTS, &CandidateLine.OtherOptCenter, RADIUS, LC);
-      // check whether other side triangle has been created
-      if (BTS != NULL) {
-        BTS->SetTopNode(CandidateLine.T);
-        if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
-          BTS->AngleFromTop = CandidateLine.ShortestAngle;
-        else
-          BTS->AngleFromTop = -1.;
-        BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter);
-        // normal vector points other way round
-        BTS->NormalVector.CopyVector(&triangle->NormalVector);
-        BTS->NormalVector.Scale(-1.);
-
-        // give some verbose output about the whole procedure
-        if (CandidateLine.T != NULL)
-          Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
-        else
-          Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
-      }
+    AddTesselationPoint(TurningPoint, 0);
+    AddTesselationPoint(*Runner, 1);
+    AddTesselationPoint(*Sprinter, 2);
+
+    // check whether we add a degenerate or a normal triangle
+    if (CheckDegeneracy(&CandidateLine.OtherOptCenter, RADIUS, LC)) {
+      // add normal and degenerate triangles
+      Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl;
+      BoundaryTriangleSet * const triangle = AddDegeneratedTriangle(CandidateLine, RADIUS, LC);
+    } else {
+      // add this triangle
+      AddCandidateTriangle(CandidateLine);
     }
 
@@ -2683 +2670 @@
 };
 
-/** If a given \a *triangle is degenerated, this adds the "other side".
+/** If a given \a *triangle is degenerated, this adds both sides.
  * i.e. the triangle with same BoundaryPointSet's but NormalVector in opposite direction.
- * \param *triangle degenerated triangle
- * \param *OtherOptCenter center of the sphere for the other triangle
+ * Note that endpoints are stored in Tesselation::TPS
+ * \param CandidateLine CanddiateForTesselation structure for the desired BoundaryLine
  * \param RADIUS radius of sphere
  * \param *LC pointer to LinkedCell structure
  * \return pointer to created other triangle
  */
-void Tesselation::AddDegeneratedTriangle(BoundaryTriangleSet *&triangle, Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC)
+BoundaryTriangleSet * const  Tesselation::AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC)
 {
   Info FunctionInfo(__func__);
-  TesselPoint *ThirdNode = NULL;
-  BoundaryTriangleSet *otherTriangle = NULL;
-  /// 1. For each of the open sides of \a *triangle find best candidate
+  BoundaryTriangleSet *triangle = NULL;
+  pair<LineMap::iterator,LineMap::iterator> FindPair[3];
+  int LineCase[3];
+  int count[3];
+  Vector Center;
+  CandidateMap::const_iterator CandidateCheck = OpenLines.end();
+
+  /// 1. For each new pair of endpoints, find and count the number of present lines
+  Log() << Verbose(0) << "INFO: Counting present open lines ..." << endl;
   for (int i=0;i<3;i++) {
     BLS[i] = NULL;
-    // go through all possible lines (also degenerate ones) and find ones to be used
-    pair<LineMap::iterator,LineMap::iterator> FindPair = triangle->lines[i]->endpoints[0]->lines.equal_range(triangle->lines[i]->endpoints[1]->Nr);
-    for (LineMap::iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
-      if ((*FindLine).second->triangles.size() == 1) {
-        Log() << Verbose(1) << "Found line " << *(*FindLine).second << " from triangle " << *((*FindLine).second->triangles.begin()->second) << "." << endl;
-        // check whether best Candidate for this line from the degenerate triangle is the third endpoint
-        CandidateForTesselation CandidateLine((*FindLine).second);
-        otherTriangle = ((*FindLine).second->triangles.begin())->second;  // this is unique and present, above we check that there is exactly one
-        ThirdNode = otherTriangle->GetThirdEndpoint((*FindLine).second)->node;
-        FindNextSuitableTriangle(CandidateLine, *triangle, RADIUS, LC);
-
-        /// 2. If best candidate is third endpoint,
-        for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner) {
-          Log() << Verbose(1) << "Found candidate " << *(*Runner) << endl;
-          if ((*Runner) == ThirdNode) { // seems like we found a line to connect to
-            /// sets Tesselation::BLS[i] accordingly and removes the line from Tesselation::OpenLines
-            AddExistingTesselationTriangleLine((*FindLine).second,i);
-            break;
-          }
-        }
+    LineCase[i] = -1;
+    // find already present lines
+    FindPair[i] = TPS[(i+0)%3]->lines.equal_range(TPS[(i+1)%3]->Nr);
+    // count their number
+    count[i] = 0;
+    Log() << Verbose(1) << "For line between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    for (LineMap::iterator FindLine = FindPair[i].first; FindLine != FindPair[i].second; FindLine++)
+      if (FindLine->second->triangles.size() == 1) {
+        Log() << Verbose(1) << "  open line " << *FindLine->second << "." << endl;
+        count[i]++;
       }
-    }
-    /// otherwise we create a new BoundaryLineSet for this Tesselation::BLS[i] and add it to OpenLines
-  }
-  if ((BLS[0] == NULL) && (BLS[1] == NULL) && (BLS[2] == NULL)) {
-    Log() << Verbose(0) << "This degenerate triangle would not share any line with its other side, rejecting." << endl;
-    BTS = NULL;
-  } else {
-    // created new lines for the other egdes of the other side triangle
-    for (int i=0;i<3;i++) {
-      if (BLS[i] == NULL) {
-        AddNewTesselationTriangleLine(triangle->lines[i]->endpoints[0],triangle->lines[i]->endpoints[1],i);
-      }
-    }
-    /// 3. Create the triangle
-    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-    AddTesselationTriangle();
-
-    // add triangle to the CandidateForTesselation of its open lines
-    CandidateMap::iterator Finder = OpenLines.end();
-    for (int i=0;i<3;i++) {
-      Finder = OpenLines.find(BTS->lines[i]);
-      if (Finder != OpenLines.end()) {
-        Log() << Verbose(1) << "INFO: Found open line of this triangle " << *(Finder->first) << endl;
-        if (Finder->second->T == NULL) {
-          Finder->second->T = BTS;
-          Log() << Verbose(1) << "ACCEPT: Setting source triangle of open line " << *(Finder->first) << " to the newly created one." << endl;
-        }
-      } else
-        Log() << Verbose(1) << "REJECT: line " << *BTS->lines[i] << " of this triangle is not found to be open." << endl;
-    }
-  }
+    Log() << Verbose(1) << "  SUM: there are " << count[i] << " line(s) present."  << endl;
+  }
+
+  /// 2. Create or pick the lines for the first triangle
+  // for each count we have a different case:
+  // case 0: no triangles at this line and not closed
+  // case 1: no triangles at new line is closed
+  // case 2: one line with one triangle attached, one will used this line
+  // case 3: one line with two triangles attached, new line will be used by both degenerate triangles
+  // case 4: two lines with one triangle each, each new triangle will pick one
+  Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl;
+  for (int i=0;i<3;i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    switch (count[i]) {
+      case 0:
+        AddNewTesselationTriangleLine(TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      case 1:
+//        // check whether present line has right candidate
+//        CandidateCheck = OpenLines.find(FindPair[i].first->second);
+//        if (CandidateCheck != OpenLines.end()) {
+//          if (!CandidateCheck->second->pointlist.empty()) {
+//            for (TesselPointList::const_iterator Runner = CandidateCheck->second->pointlist.begin(); Runner != CandidateCheck->second->pointlist.end(); ++Runner) {
+//              Log() << Verbose(1) << "Found candidate " << *(*Runner) << " with Candidate center is " << CandidateCheck->second->OptCenter << ", desired opt center is " << CandidateLine.OptCenter << "." << endl;
+//              if (((*Runner) == TPS[(i+2)%3]->node) && (CandidateLine.OptCenter == CandidateCheck->second->OptCenter)) { // seems like we found a line to connect to
+//                AddExistingTesselationTriangleLine(CandidateCheck->first, i);
+//                break;  // have to break as a AddExistingTesselationTriangleLine changes OpenLines
+//              }
+//            }
+//          } else {
+//            Log() << Verbose(1) << "Open line between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << " has no candidate." << endl;
+//          }
+//        } else {
+//          Log() << Verbose(1) << "No open lines present for line between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << "." << endl;
+//        }
+//        if (BLS[i] == NULL) {
+//          AddNewTesselationTriangleLine(TPS[(i+0)%3], TPS[(i+1)%3], i);
+//        }
+        AddTesselationLine(&CandidateLine.OptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      default:
+      case 2:
+        AddTesselationLine(&CandidateLine.OptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+//        DoeLog(0) && (eLog() << Verbose(0) << "There are more than two BoundaryPointSets between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << "." << endl);
+//        performCriticalExit();
+//        break;
+    }
+  }
+
+  /// 3. create the first triangle and NormalVector and so on
+  Log() << Verbose(0) << "INFO: Adding first triangle with center at " << CandidateLine.OptCenter << " ..." << endl;
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTesselationTriangle();
+  // check whether triangle has been created
+  BTS->SetTopNode(CandidateLine.T);
+  if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
+    BTS->AngleFromTop = CandidateLine.ShortestAngle;
+  else
+    BTS->AngleFromTop = -1.;
+  // create normal vector
+  BTS->GetCenter(&Center);
+  Center.SubtractVector(&CandidateLine.OptCenter);
+  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+  BTS->GetNormalVector(Center);
+  // give some verbose output about the whole procedure
+  if (CandidateLine.T != NULL)
+    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+  else
+    Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
+  triangle = BTS;
+
+  /// 4. Gather candidates for each new line
+  Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl;
+  for (int i=0;i<3;i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    CandidateCheck = OpenLines.find(BLS[i]);
+    if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
+      if (CandidateCheck->second->T == NULL)
+        CandidateCheck->second->T = triangle;
+      FindNextSuitableTriangle(*(CandidateCheck->second), *CandidateCheck->second->T, RADIUS, LC);
+    }
+  }
+
+  /// 5. Create or pick the lines for the second triangle
+  Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl;
+  for (int i=0;i<3;i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    switch (count[i]) {
+      case 0:
+        AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      case 1:
+//        // check whether present line has right candidate
+//        CandidateCheck = OpenLines.find(BLS[i]);
+//        BLS[i] = NULL;
+//        if (CandidateCheck == OpenLines.end()) {  // first triangle has closed present line
+//          AddNewTesselationTriangleLine(TPS[(i+0)%3], TPS[(i+1)%3], i);
+//        } else {
+//          for (TesselPointList::const_iterator Runner = CandidateCheck->second->pointlist.begin(); Runner != CandidateCheck->second->pointlist.end(); ++Runner) {
+//            Log() << Verbose(1) << "Found candidate " << *(*Runner) << " with Candidate center is " << CandidateCheck->second->OptCenter << ", desired other opt center is " << CandidateLine.OtherOptCenter << "." << endl;
+//            if (((*Runner) == TPS[(i+2)%3]->node) && (CandidateLine.OtherOptCenter == CandidateCheck->second->OptCenter)) { // seems like we found a line to connect to
+//              AddExistingTesselationTriangleLine(CandidateCheck->first, i);
+//              break;  // have to break as a AddExistingTesselationTriangleLine changes OpenLines
+//            }
+//          }
+//        }
+//        if (BLS[i] == NULL) {
+//          DoeLog(0) && (eLog() << Verbose(0) << "BoundaryLine between " << TPS[(i+0)%3] << " and " << TPS[(i+1)%3] << " has not the correct candidate, i.e. the third endpoint: " << *TPS[(i+2)%3] << endl);
+//          performCriticalExit();
+//        }
+        AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      case 2:
+      default:
+        AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+//        DoeLog(0) && (eLog() << Verbose(0) << "There are more than two BoundaryPointSets between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << "." << endl);
+//        performCriticalExit();
+//        break;
+    }
+  }
+
+  /// 6. create the second triangle and NormalVector and so on
+  Log() << Verbose(0) << "INFO: Adding second triangle with center at " << CandidateLine.OtherOptCenter << " ..." << endl;
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTesselationTriangle();
+  // check whether triangle has been created
+  BTS->SetTopNode(CandidateLine.T);
+  if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
+    BTS->AngleFromTop = CandidateLine.ShortestAngle;
+  else
+    BTS->AngleFromTop = -1.;
+  BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter);
+  // create normal vector in other direction
+  BTS->GetNormalVector(&triangle->NormalVector);
+  BTS->NormalVector.Scale(-1.);
+  // give some verbose output about the whole procedure
+  if (CandidateLine.T != NULL)
+    Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+  else
+    Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
+
+  return triangle;
 };
 
 /** Adds a triangle to the Tesselation structure from three given TesselPoint's.
- * \param *first first TesselPoint
- * \param *second second TesselPoint
- * \param *third third TesselPoint
- */
-void Tesselation::AddCandidateTriangle(TesselPoint *first, TesselPoint *second, TesselPoint *third)
+ * Note that endpoints are in Tesselation::TPS.
+ * \param CandidateLine CandidateForTesselation structure contains other information
+ */
+void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine)
 {
   Info FunctionInfo(__func__);
-  // add the points
-  AddTesselationPoint(first, 0);
-  AddTesselationPoint(second, 1);
-  AddTesselationPoint(third, 2);
+  Vector Center;
 
   // add the lines
-  AddTesselationLine(TPS[2], TPS[0], TPS[1], 0);
-  AddTesselationLine(TPS[1], TPS[0], TPS[2], 1);
-  AddTesselationLine(TPS[0], TPS[1], TPS[2], 2);
+  AddTesselationLine(&CandidateLine.OptCenter, TPS[2], TPS[0], TPS[1], 0);
+  AddTesselationLine(&CandidateLine.OptCenter, TPS[1], TPS[0], TPS[2], 1);
+  AddTesselationLine(&CandidateLine.OptCenter, TPS[0], TPS[1], TPS[2], 2);
 
   // add the triangles
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   AddTesselationTriangle();
+  // store values from tesselation procedure therein
+  BTS->SetTopNode(CandidateLine.T);
+  if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
+    BTS->AngleFromTop = CandidateLine.ShortestAngle;
+  else
+    BTS->AngleFromTop = -1.;
+
+  // create normal vector
+  BTS->GetCenter(&Center);
+  Center.SubtractVector(&CandidateLine.OptCenter);
+  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+  BTS->GetNormalVector(Center);
+
+  // give some verbose output about the whole procedure
+  if (CandidateLine.T != NULL)
+    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+  else
+    Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
 };
 
@@ -3163 +3273 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
+void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
 {
         Info FunctionInfo(__func__);
@@ -3344 +3454 @@
   }
 
-  if (!CandidateLine.CheckValidity(RADIUS, LC)) {
+  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
     DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
     performCriticalExit();
@@ -4349 +4459 @@
         AddTesselationPoint(*EndNode, 2);
         Log() << Verbose(3) << "Adding new triangle lines."<< endl;
-        AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-        AddTesselationLine(NULL, TPS[0], TPS[2], 1);
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
         NewLines.push_back(BLS[1]);
-        AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
         BTS->GetNormalVector(NormalVector);
@@ -4801 +4911 @@
   AddTesselationPoint(point, 2);
   Log() << Verbose(2) << "Adding new triangle lines."<< endl;
-  AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-  AddTesselationLine(NULL, TPS[0], TPS[2], 1);
-  AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
@@ -4816 +4926 @@
   AddTesselationPoint(point, 2);
   Log() << Verbose(2) << "Adding new triangle lines."<< endl;
-  AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-  AddTesselationLine(NULL, TPS[0], TPS[2], 1);
-  AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
@@ -5051 +5161 @@
       for (int i = 0; i < 3; i++)
         AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);
-      AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-      AddTesselationLine(NULL, TPS[0], TPS[2], 1);
-      AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+      AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+      AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+      AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
       if (TriangleNrs.empty())
         DoeLog(0) && (eLog()<< Verbose(0) << "No more free triangle numbers!" << endl);

src/tesselation.hpp

@@ -163 +163 @@
     bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
     bool IsPresentTupel(const BoundaryTriangleSet * const T) const;
-    bool CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const;
     void SetTopNode(const BoundaryTriangleSet * const topnode);
 
@@ -283 +282 @@
     void AddTesselationPoint(TesselPoint* Candidate, const int n);
     void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
-    void AddTesselationLine(const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
+    void AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddExistingTesselationTriangleLine(class BoundaryLineSet *FindLine, int n);
     void AddTesselationTriangle();
     void AddTesselationTriangle(const int nr);
-    void AddCandidateTriangle(TesselPoint *first, TesselPoint *second, TesselPoint *third);
-    void AddDegeneratedTriangle(BoundaryTriangleSet *&triangle, Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC);
+    void AddCandidateTriangle(CandidateForTesselation &CandidateLine);
+    BoundaryTriangleSet * const AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC);
     void AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC);
     void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
     void RemoveTesselationLine(class BoundaryLineSet *line);
     void RemoveTesselationPoint(class BoundaryPointSet *point);
+    bool CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const;
 
 
@@ -299 +299 @@
     void FindStartingTriangle(const double RADIUS, const LinkedCell *LC);
     void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC);
-    void FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
-    bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
+    void FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
+    bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
     int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]) const;
     class BoundaryTriangleSet * GetPresentTriangle(TesselPoint *Candidates[3]);