Changeset b998c3 for src

Timestamp:

Dec 16, 2009, 11:53:06 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:

009e37b

Parents:

5e8f82

Message:

Fix attempts of the Tesselation.

One problem still remains:

• degenerated triangles are created more than two times.

The following has been changed:

• BondGraph::~BondGraph() - docu changes
• FindNonConvexBorder() - added CheckListOfBaseLines() and StoreTrianglesInFile() before CorrectDegeneratePolygons() (allows checking of closed tesselation even when correction fails)
• boundary.hpp: SingleStepWidth set to 10
• Tesselation::AddTesselationLine() - FindLine set before if and hence we may use it twice and find() only once
• Tesselation::FindStartingTriangle() - renaming of variables such that they are equal to FindNextSuitableTriangle()
• Tesselation::FindNextSuitableTriangle() - using SphereCenter (is the result of GetCenterofCircumcenter())
• Tesselation::FindThirdPointForTesselation() - some rewrites to make the code clearer
• CandidateForTesselation::IsDegenerated - designates candidate TesselPoint of a degenerated triangle
• BoundaryTriangleSet::SphereCenter - contains the SphereCenter created during Tesselation of the triangle
• FIX: Tesselation::CorrectAllDegeneratedPolygons() - only adds triangles to candidate list that are _not_ simply degenerated (this prevents two(four) singly degenerated but separate triangles to be put into the same polygon, resulting in an uneven number of triangles).
• new function Vector::IsEqualTo() - checks whether vector is equal to another (within MYEPSILON bounds)

Location:

src

Files:

• bondgraph.cpp (modified) (1 diff)
• boundary.cpp (modified) (1 diff)
• boundary.hpp (modified) (1 diff)
• tesselation.cpp (modified) (19 diffs)
• tesselation.hpp (modified) (2 diffs)
• tesselationhelpers.cpp (modified) (2 diffs)
• vector.cpp (modified) (1 diff)
• vector.hpp (modified) (1 diff)

src/bondgraph.cpp

@@ -35 +35 @@
 /** Parses the bond lengths in a given file and puts them int a matrix form.
  * Allocates \a MatrixContainer for BondGraph::BondLengthMatrix, using MatrixContainer::ParseMatrix(),
- * but only if parsing is successfull. Otherwise variable is left as NULL.
+ * but only if parsing is successful. Otherwise variable is left as NULL.
  * \param *out output stream for debugging
  * \param filename file with bond lengths to parse

src/boundary.cpp

@@ -1041 +1041 @@
 //  TesselStruct->RemoveDegeneratedTriangles();
 
+  // check envelope for consistency
+  status = CheckListOfBaselines(TesselStruct);
+
+  // store before correction
+  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
+
   // correct degenerated polygons
   TesselStruct->CorrectAllDegeneratedPolygons();

src/boundary.hpp

@@ -36 +36 @@
 #define DEBUG 1
 #define DoSingleStepOutput 0
-#define SingleStepWidth 1
+#define SingleStepWidth 10
 
 #define DistancePair pair < double, atom* >

src/tesselation.cpp

@@ -879 +879 @@
 CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
     BaseLine(line),
+    IsDegenerated(false),
     ShortestAngle(2.*M_PI),
     OtherShortestAngle(2.*M_PI)
@@ -1587 +1588 @@
   bool insertNewLine = true;
 
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  if (FindLine != a->lines.end()) {
+    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
 
     for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
@@ -1915 +1917 @@
   double maxCoordinate[NDIM];
   BoundaryLineSet BaseLine;
-  Vector Oben;
   Vector helper;
   Vector Chord;
   Vector SearchDirection;
-
-  Oben.Zero();
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector SphereCenter;
+  Vector NormalVector;
+
+  NormalVector.Zero();
 
   for (i = 0; i < 3; i++) {
@@ -1955 +1960 @@
   BTS = NULL;
   for (int k=0;k<NDIM;k++) {
-    Oben.Zero();
-    Oben.x[k] = 1.;
+    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;
@@ -1963 +1968 @@
     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, Oben, 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);
 
-    helper.CopyVector(BaseLine.endpoints[0]->node->node);
-    helper.SubtractVector(BaseLine.endpoints[1]->node->node);
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
+    // construct center of circle
+    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);
+
+    double radius = CirclePlaneNormal.NormSquared();
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+
+    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
+    NormalVector.Normalize();
     ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
 
-    Chord.CopyVector(BaseLine.endpoints[0]->node->node); // bring into calling function
-    Chord.SubtractVector(BaseLine.endpoints[1]->node->node);
-    double radius = Chord.ScalarProduct(&Chord);
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
-    helper.CopyVector(&Oben);
-    helper.Scale(CircleRadius);
-    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+    SphereCenter.CopyVector(&NormalVector);
+    SphereCenter.Scale(CircleRadius);
+    SphereCenter.AddVector(&CircleCenter);
+    // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -1993 +2003 @@
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
     CandidateForTesselation OptCandidates(&BaseLine);
-    FindThirdPointForTesselation(Oben, SearchDirection, helper, OptCandidates, NULL, RADIUS, LC);
+    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
     Log() << Verbose(0) << "List of third Points is:" << endl;
     for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
@@ -2167 +2177 @@
   Vector CircleCenter;
   Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
+  Vector RelativeSphereCenter;
   Vector SearchDirection;
   Vector helper;
@@ -2174 +2184 @@
   double radius, CircleRadius;
 
-  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl;
   for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node))
+    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) {
       ThirdNode = T.endpoints[i]->node;
+      break;
+    }
+  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl;
 
   // construct center of circle
@@ -2191 +2203 @@
   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
   if (radius/4. < RADIUS*RADIUS) {
+    // construct relative sphere center with now known CircleCenter
+    RelativeSphereCenter.CopyVector(&T.SphereCenter);
+    RelativeSphereCenter.SubtractVector(&CircleCenter);
+
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
     Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
-    // construct old center
-    GetCenterofCircumcircle(&OldSphereCenter, *T.endpoints[0]->node->node, *T.endpoints[1]->node->node, *T.endpoints[2]->node->node);
-    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl;
+
+    // construct SearchDirection and an "outward pointer"
+    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
     Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
@@ -2219 +2226 @@
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
 
   } else {
@@ -2350 +2357 @@
     AddTesselationPoint((*Sprinter), 2);
 
-    Center.CopyVector(&CandidateLine.OptCenter);
+
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2359 +2366 @@
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     AddTesselationTriangle();
-    Center.Scale(-1.);
+    BTS->GetCenter(&Center);
+    Center.SubtractVector(&CandidateLine.OptCenter);
+    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
     BTS->GetNormalVector(Center);
 
@@ -2766 +2775 @@
   Vector NewNormalVector;   // normal vector of the Candidate's triangle
   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
+  Vector RelativeOldSphereCenter;
+  Vector NewPlaneCenter;
   double CircleRadius; // radius of this circle
   double radius;
+  double otherradius;
   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+  bool IsDegenerated;
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
   TesselPoint *Candidate = NULL;
+  TesselPoint *CandidateTriangle[3];
 
   Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
@@ -2783 +2797 @@
   CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+
+  CandidateTriangle[0] = CandidateLine.BaseLine->endpoints[0]->node;
+  CandidateTriangle[1] = CandidateLine.BaseLine->endpoints[1]->node;
+
   // calculate squared radius TesselPoint *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
+  radius = CirclePlaneNormal.NormSquared()/4.;
+  if (radius < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius;
     CirclePlaneNormal.Normalize();
-    //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
     // test whether old center is on the band's plane
-    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      eLog() << Verbose(1) << "Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
-    }
-    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+      Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;
 
       // check SearchDirection
-      //Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
@@ -2832 +2852 @@
 
                 // check for three unique points
-                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << *(Candidate->node) << "." << endl;
+                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl;
                 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){
 
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+                  // find center on the plane
+                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
+                  Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
 
                   if ((NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
                   ) {
-                    helper.CopyVector(&NewNormalVector);
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                    Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl;
                     if (radius < RADIUS*RADIUS) {
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
+                      if (fabs(radius - otherradius) > HULLEPSILON) {
+                        eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;
+                      }
+                      // construct both new centers
+                      NewSphereCenter.CopyVector(&NewPlaneCenter);
+                      OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
+                      helper.CopyVector(&NewNormalVector);
                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
                       NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-
                       // OtherNewSphereCenter is created by the same vector just in the other direction
                       helper.Scale(-1.);
                       OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
 
@@ -2861 +2888 @@
                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                       alpha = min(alpha, Otheralpha);
-                      // if there is a better candidate, drop the current list and add the new candidate
-                      // otherwise ignore the new candidate and keep the list
-                      if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
-                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                          CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
-                          CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+
+                      CandidateTriangle[2] = Candidate;
+                      // the idea of the IsDegenerated flag is not to put a penalty on degenerated triangles, but to push them to the
+                      // very end of the Tesselation::OpenLines list.
+                      IsDegenerated = (CheckPresenceOfTriangle(CandidateTriangle) >= 3);
+                      if (!IsDegenerated && CandidateLine.IsDegenerated) {  // if current is not, but old one was, comparison would be unfair
+                        // if there is a better candidate, drop the current list and add the new candidate
+                        // otherwise ignore the new candidate and keep the list
+                        if (CandidateLine.ShortestAngle-2.*M_PI > (alpha - HULLEPSILON)) {
+                          if (fabs(alpha - Otheralpha) > MYEPSILON) {
+                            CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
+                            CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                          } else {
+                            CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
+                            CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
+                          }
+                          // if there is an equal candidate, add it to the list without clearing the list
+                          if ((CandidateLine.ShortestAngle-2.*M_PI - HULLEPSILON) < alpha) {
+                            CandidateLine.pointlist.push_back(Candidate);
+                            Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with "
+                              << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                          } else {
+                            // remove all candidates from the list and then the list itself
+                            CandidateLine.pointlist.clear();
+                            CandidateLine.pointlist.push_back(Candidate);
+                            Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with "
+                              << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                          }
+                          CandidateLine.ShortestAngle = alpha;
+                          CandidateLine.IsDegenerated = IsDegenerated;
+                          if (IsDegenerated)
+                            CandidateLine.ShortestAngle += 2.*M_PI;
+                          Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl;
                         } else {
-                          CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
-                          CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
+                          if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
+                            Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;
+                          } else {
+                            Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                          }
                         }
-                        // if there is an equal candidate, add it to the list without clearing the list
-                        if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
-                          CandidateLine.pointlist.push_back(Candidate);
-                          Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with "
-                            << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                      } else {
+                        // if there is a better candidate, drop the current list and add the new candidate
+                        // otherwise ignore the new candidate and keep the list
+                        if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
+                          if (fabs(alpha - Otheralpha) > MYEPSILON) {
+                            CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
+                            CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                          } else {
+                            CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
+                            CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
+                          }
+                          // if there is an equal candidate, add it to the list without clearing the list
+                          if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
+                            CandidateLine.pointlist.push_back(Candidate);
+                            Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with "
+                              << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                          } else {
+                            // remove all candidates from the list and then the list itself
+                            CandidateLine.pointlist.clear();
+                            CandidateLine.pointlist.push_back(Candidate);
+                            Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with "
+                              << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                          }
+                          CandidateLine.ShortestAngle = alpha;
+                          Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl;
                         } else {
-                          // remove all candidates from the list and then the list itself
-                          CandidateLine.pointlist.clear();
-                          CandidateLine.pointlist.push_back(Candidate);
-                          Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with "
-                            << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
-                        }
-                        CandidateLine.ShortestAngle = alpha;
-                        Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl;
-                      } else {
-                        if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
-                          Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;
-                        } else {
-                          Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                          if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
+                            Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;
+                          } else {
+                            Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                          }
                         }
                       }
@@ -4128 +4197 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
+  map <int, int> *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4138 +4208 @@
     for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)
       for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
-        TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
-        if (TriangleInsertionTester.second)
-          Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
+          TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
+          if (TriangleInsertionTester.second)
+            Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        } else {
+          Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl;
+        }
       }
     // check whether there are two that are parallel
@@ -4213 +4287 @@
     /// 4a. Gather all triangles of this polygon
     TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();
+
+    if (T->size() == 2) {
+      Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;
+      delete(T);
+      continue;
+    }
 
     TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator

src/tesselation.hpp

@@ -158 +158 @@
     class BoundaryLineSet *lines[3];
     Vector NormalVector;
+    Vector SphereCenter;
     int Nr;
 };
@@ -245 +246 @@
   Vector OptCenter;
   Vector OtherOptCenter;
+  bool IsDegenerated;
   double ShortestAngle;
   double OtherShortestAngle;

src/tesselationhelpers.cpp

@@ -226 +226 @@
   Vector helper;
   double radius, alpha;
-
-  helper.CopyVector(&NewSphereCenter);
+  Vector RelativeOldSphereCenter;
+  Vector RelativeNewSphereCenter;
+
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
+  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
+  helper.CopyVector(&RelativeNewSphereCenter);
   // test whether new center is on the parameter circle's plane
   if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
@@ -233 +239 @@
     helper.ProjectOntoPlane(&CirclePlaneNormal);
   }
-  radius = helper.ScalarProduct(&helper);
+  radius = helper.NormSquared();
   // test whether the new center vector has length of CircleRadius
   if (fabs(radius - CircleRadius) > HULLEPSILON)
     eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-  alpha = helper.Angle(&OldSphereCenter);
+  alpha = helper.Angle(&RelativeOldSphereCenter);
   // make the angle unique by checking the halfplanes/search direction
   if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
     alpha = 2.*M_PI - alpha;
-  //Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&OldSphereCenter);
+  Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&RelativeOldSphereCenter);
   helper.ProjectOntoPlane(&NormalVector);
   // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
   if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
-    //Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
     return alpha;
   } else {
-    //Log() << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    Log() << Verbose(1) << "INFO: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "." << endl;
     return 2.*M_PI;
   }

src/vector.cpp

@@ -480 +480 @@
   else
     return false;
+};
+
+/** Checks whether vector is normal to \a *normal.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool Vector::IsEqualTo(const Vector * const a) const
+{
+  bool status = true;
+  for (int i=0;i<NDIM;i++) {
+    if (fabs(x[i] - a->x[i]) > MYEPSILON)
+      status = false;
+  }
+  return status;
 };
 

src/vector.hpp

@@ -42 +42 @@
   bool IsOne() const;
   bool IsNormalTo(const Vector * const normal) const;
+  bool IsEqualTo(const Vector * const a) const;
 
   void AddVector(const Vector * const y);