Changeset a8c8ba for molecuilder/src

Timestamp:

Jul 9, 2009, 1:33:39 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

6fb785, b65771

Parents:

8bbe43

Message:

FIX: Now tries each of the axis direction to find a starting triangle

• if BTS was not set, as no third point was found, output of NormalVector caused error. This is fixed.
• also, we move in a loop over all three axis direction and try to create a starting triangle.
• Therefore, the candidates in the list have to be free'd. This is done.

File:

• molecuilder/src/boundary.cpp (modified) (4 diffs)

molecuilder/src/boundary.cpp

@@ -2248 +2248 @@
   }
   // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-  cout << Verbose(2) << "LC Intervals from [";
+  cout << Verbose(3) << "LC Intervals from [";
   for (int i=0;i<NDIM;i++) {
   cout << " " << N[i] << "<->" << LC->N[i];
@@ -2296 +2296 @@
                 angle = AngleCheck.Angle(&Oben);
                 if (angle < Storage[0]) {
-                  //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                  cout << Verbose(2) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
+                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
                   Opt_Candidate = Candidate;
                   Storage[0] = angle;
-                  //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
                 } else {
-                  //cout << Verbose(2) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
                 }
               } else {
-                //cout << Verbose(2) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
+                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
               }
             } else {
-              //cout << Verbose(2) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
+              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
             }
           }
         } else {
-          cout << "Linked cell list is empty." << endl;
+          cout << Verbose(3) << "Linked cell list is empty." << endl;
         }
       }
@@ -2371 +2371 @@
     cout << i << ": " << *MaxAtom[i] << "\t";
   cout << endl;
-  const int k = 1;    // arbitrary choice
-  Oben.x[k] = 1.;
-  FirstPoint = MaxAtom[k];
-  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
-
-  double ShortestAngle;
-  atom* Opt_Candidate = NULL;
-  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.
-
-  Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
-  SecondPoint = Opt_Candidate;
-  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-
-  helper.CopyVector(&(FirstPoint->x));
-  helper.SubtractVector(&(SecondPoint->x));
-  helper.Normalize();
-  Oben.ProjectOntoPlane(&helper);
-  Oben.Normalize();
-  helper.VectorProduct(&Oben);
-  ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-
-  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-  Chord.SubtractVector(&(SecondPoint->x));
-  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)
-
-  // look in one direction of baseline for initial candidate
+
+  BTS = NULL;
   CandidateList *Opt_Candidates = new CandidateList();
-  SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
-
-  // adding point 1 and point 2 and the line between them
-  AddTrianglePoint(FirstPoint, 0);
-  AddTrianglePoint(SecondPoint, 1);
-  AddTriangleLine(TPS[0], TPS[1], 0);
-
-  //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
-  Find_third_point_for_Tesselation(
-    Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
-  );
-  cout << Verbose(1) << "List of third Points is ";
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-      cout << " " << *(*it)->point;
-  }
-  cout << endl;
-
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-          // add third triangle point
-          AddTrianglePoint((*it)->point, 2);
-          // add the second and third line
-          AddTriangleLine(TPS[1], TPS[2], 1);
-          AddTriangleLine(TPS[0], TPS[2], 2);
-          // ... and triangles to the Maps of the Tesselation class
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTriangle();
-          // ... and calculate its normal vector (with correct orientation)
-          (*it)->OptCenter.Scale(-1.);
-          cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
-          BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
-          cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
-                << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
-
-          // if we do not reach the end with the next step of iteration, we need to setup a new first line
-          if (it != Opt_Candidates->end()--) {
-                  FirstPoint = (*it)->BaseLine->endpoints[0]->node;
-                  SecondPoint = (*it)->point;
-                  // adding point 1 and point 2 and the line between them
-                  AddTrianglePoint(FirstPoint, 0);
-                  AddTrianglePoint(SecondPoint, 1);
-                  AddTriangleLine(TPS[0], TPS[1], 0);
-          }
-  }
+  for (int k=0;k<NDIM;k++) {
+    Oben.x[k] = 1.;
+    FirstPoint = MaxAtom[k];
+    cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
+
+    double ShortestAngle;
+    atom* Opt_Candidate = NULL;
+    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.
+
+    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    SecondPoint = Opt_Candidate;
+    if (SecondPoint == NULL)  // have we found a second point?
+      continue;
+    else
+      cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+
+    helper.CopyVector(&(FirstPoint->x));
+    helper.SubtractVector(&(SecondPoint->x));
+    helper.Normalize();
+    Oben.ProjectOntoPlane(&helper);
+    Oben.Normalize();
+    helper.VectorProduct(&Oben);
+    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+
+    Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+    Chord.SubtractVector(&(SecondPoint->x));
+    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)
+
+    // look in one direction of baseline for initial candidate
+    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+
+    // adding point 1 and point 2 and the line between them
+    AddTrianglePoint(FirstPoint, 0);
+    AddTrianglePoint(SecondPoint, 1);
+    AddTriangleLine(TPS[0], TPS[1], 0);
+
+    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
+    Find_third_point_for_Tesselation(
+      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
+    );
+    cout << Verbose(1) << "List of third Points is ";
+    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+        cout << " " << *(*it)->point;
+    }
+    cout << endl;
+
+    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+      // add third triangle point
+      AddTrianglePoint((*it)->point, 2);
+      // add the second and third line
+      AddTriangleLine(TPS[1], TPS[2], 1);
+      AddTriangleLine(TPS[0], TPS[2], 2);
+      // ... and triangles to the Maps of the Tesselation class
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      AddTriangle();
+      // ... and calculate its normal vector (with correct orientation)
+      (*it)->OptCenter.Scale(-1.);
+      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
+      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
+      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
+      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
+
+      // if we do not reach the end with the next step of iteration, we need to setup a new first line
+      if (it != Opt_Candidates->end()--) {
+        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
+        SecondPoint = (*it)->point;
+        // adding point 1 and point 2 and the line between them
+        AddTrianglePoint(FirstPoint, 0);
+        AddTrianglePoint(SecondPoint, 1);
+        AddTriangleLine(TPS[0], TPS[1], 0);
+      }
+      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+    }
+    if (BTS != NULL) // we have created one starting triangle
+      break;
+    else {
+      // remove all candidates from the list and then the list itself
+      class CandidateForTesselation *remover = NULL;
+      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+        remover = *it;
+        delete(remover);
+      }
+      Opt_Candidates->clear();
+    }
+  }
+
   // remove all candidates from the list and then the list itself
   class CandidateForTesselation *remover = NULL;
@@ -2452 +2471 @@
   }
   delete(Opt_Candidates);
-  cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
   cout << Verbose(1) << "End of Find_starting_triangle\n";
 };