Changeset 02bfde

Timestamp:

Dec 23, 2008, 1:06:36 PM (16 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:

12298c, 23e09b, cc2ee5

Parents:

44fd95 (diff), 196a5a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'ConcaveHull' of ​ssh://stud64dc-01/home/neuen/workspace/ESPACK into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp

+ Merge is solved (was Christian's fault :)
+ find_non_convex_hull is checked and working on Heptan boundary

File:

• src/boundary.cpp (modified) (26 diffs)

src/boundary.cpp

@@ -1356 +1356 @@
   LineMap::iterator LineWalker;
   //cout << "Manually checking endpoints for line." << endl;
-  if ((a->lines.find(b->node->nr))->first == b->node->nr)
+  if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr)
   //If a line is there, how do I recognize that beyond a shadow of a doubt?
     {
@@ -1421 +1421 @@
  * @param b vector second point of triangle
  * @param c vector third point of triangle
+ * @param *Umkreismittelpunkt new cneter point of circumference
  * @param Direction vector indicates up/down
+ * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
  * @param Halfplaneindicator double indicates whether Direction is up or down
+ * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
  * @param alpha double angle at a
  * @param beta double, angle at b
@@ -1430 +1433 @@
  */
 
-  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, double HalfplaneIndicator,
-      double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
+      double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
   {
     Vector TempNormal, helper;
@@ -1437 +1440 @@
     cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
     *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
-    Center->Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+    NewUmkreismittelpunkt->CopyVector(Center);
+    cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
     // Here we calculated center of circumscribing circle, using barycentric coordinates
     cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
@@ -1446 +1451 @@
     helper.SubtractVector(&c);
     TempNormal.VectorProduct(&helper);
-    if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+    if (fabs(HalfplaneIndicator) < MYEPSILON)
       {
-        TempNormal.Scale(-1);
+        if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+          {
+            TempNormal.Scale(-1);
+          }
+      }
+    else
+      {
+        if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+          {
+            TempNormal.Scale(-1);
+          }
       }
 
@@ -1474 +1489 @@
  * @param a first point
  * @param b second point
+ * *param c atom old third point of old triangle
  * @param Candidate base point along whose bonds to start looking from
  * @param Parent point to avoid during search as its wrong direction
@@ -1487 +1503 @@
  */
 
-void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* Candidate, atom* Parent,
+void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
     atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
@@ -1509 +1525 @@
   double CurrentEpsilon = 0.1;
   double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-  double BallAngle;
+  double BallAngle, AlternativeSign;
   atom *Walker; // variable atom point
 
-
-  if (a != Candidate and b != Candidate)
+  Vector NewUmkreismittelpunkt;
+
+
+  if (a != Candidate and b != Candidate and c != Candidate)
     {
       cout << Verbose(3) << "We have a unique candidate!" << endl;
@@ -1546 +1564 @@
         {
           cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
+          cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
           sign = AngleCheck.ScalarProduct(direction1);
-          if (fabs(sign) < MYEPSILON)
-            sign = 0;
+          if (fabs(sign)<MYEPSILON)
+            {
+              if (AngleCheck.ScalarProduct(OldNormal)<0)
+                {
+                  sign =0;
+                  AlternativeSign=1;
+                }
+              else
+                {
+                  sign =0;
+                  AlternativeSign=-1;
+                }
+            }
           else
-            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-          if (sign >= 0)
             {
-            cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
-
-            Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, sign, alpha, beta, gamma, RADIUS, Umkreisradius);
-
-            AngleCheck.CopyVector(&ReferencePoint);
-            AngleCheck.Scale(-1);
-            //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-            AngleCheck.AddVector(&Mittelpunkt);
-            //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-            cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
-
-            BallAngle = AngleCheck.Angle(OldNormal);
-            cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-  //          cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
-  //          cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-            sign = AngleCheck.ScalarProduct(direction1);
-            if (fabs(sign) < MYEPSILON)
-              sign = 0;
-            else
               sign /= fabs(sign);
-
-
+              cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
+            }
+
+          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+
+          AngleCheck.CopyVector(&ReferencePoint);
+          AngleCheck.Scale(-1);
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          AngleCheck.AddVector(&Mittelpunkt);
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
+
+          BallAngle = AngleCheck.Angle(OldNormal);
+          cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+
+          //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+          cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+
+          NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
+
+          if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON))
+            {
               if (Storage[0]< -1.5) // first Candidate at all
                 {
@@ -1581 +1612 @@
                   Opt_Candidate = Candidate;
                   Storage[0] = sign;
-                  Storage[1] = BallAngle;
-                  cout << " angle " << Storage[1] << " and Up/Down "
+                  Storage[1] = AlternativeSign;
+                  Storage[2] = BallAngle;
+                  cout << " angle " << Storage[2] << " and Up/Down "
                     << Storage[0] << endl;
                 }
               else
                 {
-                  if ( Storage[1] > BallAngle)
+                  if ( Storage[2] > BallAngle)
                     {
                       cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
                       Opt_Candidate = Candidate;
                       Storage[0] = sign;
-                      Storage[1] = BallAngle;
-                      cout << " angle " << Storage[1] << " and Up/Down "
+                      Storage[1] = AlternativeSign;
+                      Storage[2] = BallAngle;
+                      cout << " angle " << Storage[2] << " and Up/Down "
                         << Storage[0] << endl;
                     }
@@ -1631 +1664 @@
 
 
-  if (RecursionLevel < 7) // Five is the recursion level threshold.
+  if (RecursionLevel < 9) // Seven is the recursion level threshold.
     {
       for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
@@ -1643 +1676 @@
           else
             {
-              Find_next_suitable_point_via_Angle_of_Sphere(a, b, Walker, Candidate, RecursionLevel
+              Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel
                   + 1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS,
                   mol); //call function again
@@ -1940 +1973 @@
   double tmp;
   //atom* Walker;
-
-  double Storage[2];
+  atom* OldThirdPoint;
+
+  double Storage[3];
   Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[1] = 9999999.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[2] = 9999999.;
   atom* Opt_Candidate = NULL;
   Vector Opt_Mittelpunkt;
@@ -1955 +1990 @@
           && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr)
         {
+          OldThirdPoint = T.endpoints[i]->node;
           helper.SubtractVector(&T.endpoints[i]->node->x);
           break;
@@ -1984 +2020 @@
   c.SubtractVector(&(T.endpoints[0]->node->x));
 
-  alpha = a.Angle(&c) - M_PI/2.;
-  beta = b.Angle(&a);
-  gamma = c.Angle(&b) - M_PI/2.;
-  if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-    cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+//  alpha = a.Angle(&c) - M_PI/2.;
+//  beta = b.Angle(&a);
+//  gamma = c.Angle(&b) - M_PI/2.;
+    alpha = M_PI - a.Angle(&c);
+    beta = M_PI - b.Angle(&a);
+    gamma = M_PI - c.Angle(&b);
+    if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
+      cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+
+    Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+    //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
+    cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
 
   cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
@@ -2010 +2056 @@
   cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
 
-  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node,
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
       Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
       &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node,
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
       Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
       &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
 
 
-  if ((TrianglesOnBoundaryCount % 10) == 0)
+      cout << "Letzter Winkel bei " << TrianglesOnBoundaryCount << " Winkel ist " << Storage[2] << endl;
+
+
+  if ((TrianglesOnBoundaryCount % 1) == 0)
     {
     cout << Verbose(1) << "Writing temporary non convex hull to file testEnvelope-" << TriangleFilesWritten << "d.dat.\n";
@@ -2029 +2078 @@
   }
 
-  if (TrianglesOnBoundaryCount >1000 )
+  if (TrianglesOnBoundaryCount >189 and TrianglesOnBoundaryCount < 200 )
     {
       cerr << Verbose(0)
@@ -2040 +2089 @@
   cout << Verbose(2) << " Optimal candidate is " << *Opt_Candidate << endl;
 
+  // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+
   AddTrianglePoint(Opt_Candidate, 0);
-  AddTrianglePoint(Line.endpoints[0]->node, 1);
-  AddTrianglePoint(Line.endpoints[1]->node, 2);
-
-  AddTriangleLine(TPS[0], TPS[1], 0);
-  AddTriangleLine(TPS[0], TPS[2], 1);
-  AddTriangleLine(TPS[1], TPS[2], 2);
-
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  AddTriangleToLines();
-
-  BTS->GetNormalVector(T.NormalVector);
-
-//  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
-//      (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0))
-//    {
-//      BTS->NormalVector.Scale(-1);
-//    };
-  cout << Verbose(2) << "New triangle with " << *BTS << "and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-  cout << Verbose(2) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
+  LineMap::iterator TryAndError;
+  TryAndError = TPS[0]->lines.find(Line.endpoints[0]->node->nr);
+  bool flag = true;
+  if ((TryAndError != TPS[0]->lines.end()) && (TryAndError->second->TrianglesCount > 1))
+    flag = false;
+  TryAndError = TPS[0]->lines.find(Line.endpoints[1]->node->nr);
+  if ((TryAndError !=  TPS[0]->lines.end()) && (TryAndError->second->TrianglesCount > 1))
+    flag = false;
+
+  if (flag) { // if so, add
+    AddTrianglePoint(Line.endpoints[0]->node, 1);
+    AddTrianglePoint(Line.endpoints[1]->node, 2);
+
+    AddTriangleLine(TPS[0], TPS[1], 0);
+    AddTriangleLine(TPS[0], TPS[2], 1);
+    AddTriangleLine(TPS[1], TPS[2], 2);
+
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    AddTriangleToLines();
+
+    BTS->GetNormalVector(BTS->NormalVector);
+
+    if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
+        (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) ||
+        (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) )
+
+      {
+        BTS->NormalVector.Scale(-1);
+      };
+    cout << Verbose(2) << "New triangle with " << *BTS << "and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
+    cout << Verbose(2) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
+  } else { // else, yell and do nothing
+    cout << Verbose(2) << "This triangle consisting of ";
+    for (int i=0;i<NDIM;i++)
+      cout << BLS[i] << "(" << BLS[i]->TrianglesCount << "), ";
+    cout << " is invalid!" << endl;
+  }
   cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
 }
@@ -2065 +2134 @@
 
 void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[2],
+    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[3],
     molecule* mol, double RADIUS)
 {
@@ -2093 +2162 @@
               Opt_Candidate = Candidate;
               Storage[0] = AngleCheck.Angle(&Oben);
-              //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+              //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
             }
           else
@@ -2176 +2245 @@
 
   cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
-  double Storage[2];
+  double Storage[3];
   atom* Opt_Candidate = NULL;
   Storage[0] = 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.
   Storage[1] = 999999.; // This will be an angle looking for the third point.
+  Storage[2] = 999999.;
 
   Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
@@ -2194 +2264 @@
   Storage[0] = -2.; // This will indicate the quadrant.
   Storage[1] = 9999999.; // This will be an angle looking for the third point.
+  Storage[2] = 9999999.;
 
   Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
@@ -2203 +2274 @@
   Opt_Candidate = NULL;
   CenterOfFirstLine.CopyVector(&Chord);
-  CenterOfFirstLine.Scale(-0.5);
+  CenterOfFirstLine.Scale(0.5);
   CenterOfFirstLine.AddVector(&(SecondPoint->x));
 
   cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
-  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0,
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0,
       &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
   // look in other direction of baseline for possible better candidate
   cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
-  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, FirstPoint, SecondPoint, 0,
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0,
       &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
   cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
@@ -2252 +2323 @@
 
   baseline = Tess->LinesOnBoundary.begin();
-  while ((baseline != Tess->LinesOnBoundary.end()) || (!flag))
+  while ((baseline != Tess->LinesOnBoundary.end()) || (flag))
     {
       if (baseline->second->TrianglesCount == 1)
@@ -2263 +2334 @@
       else
         {
-          cout << Verbose(1) << "Line " << &baseline << " has "
+          cout << Verbose(1) << "Line " << baseline->second << " has "
               << baseline->second->TrianglesCount << " triangles adjacent"
               << endl;
@@ -2280 +2351 @@
 }
 ;
+