Changeset 1e168b for src/tesselation.cpp

Timestamp:

Nov 25, 2009, 12:18:18 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

8725ed

Parents:

791138

Message:

Now the best (in terms of ShortestAngle) baseline is picked among all open ones instead of the next.

• trying to fix the tesselation, we have the following changes
  • in Tesselation::FindNextSuitableTriangle(), the adding of the triangle was split off and shifted to new function Tesselation::AddCandidateTriangle()
  • Tesselation::FindNextSuitableTriangle() now takes CandidateForTesselation instead of a baseline and puts found candidate therein
  • instead of looping through LinesOnBoundary in FindNonConvexBorder() we have a new Map called OpenLines containing all Lines with TriangleCount < 2
  • the loop itself was changed as follows and continues until OpenLines is empty:
    1. Find Candidates for all new OpenLines
    2. Pick smallest ShortestAngle among these
    3. Add this candidate
  • OpenLines is added to in AddAlwaysTesselationLine() and removed from in AllTesselationLine() (existing line case)
  • CandidateForTesselation class has new attributes (Other)ShortestAngle

• Element initialization was done for all constructors in tesselation.cpp

File:

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

src/tesselation.cpp

@@ -22 +22 @@
 /** Constructor of BoundaryPointSet.
  */
-BoundaryPointSet::BoundaryPointSet()
-{
-  LinesCount = 0;
-  Nr = -1;
-  value = 0.;
+BoundaryPointSet::BoundaryPointSet() :
+    LinesCount(0),
+    value(0.),
+    Nr(-1)
+{
 };
 
@@ -85 +85 @@
 /** Constructor of BoundaryLineSet.
  */
-BoundaryLineSet::BoundaryLineSet()
+BoundaryLineSet::BoundaryLineSet() :
+    Nr(-1)
 {
   for (int i = 0; i < 2; i++)
     endpoints[i] = NULL;
-  Nr = -1;
 };
 
@@ -106 +106 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
+  // set skipped to false
+  skipped = false;
   // clear triangles list
   Log() << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
@@ -286 +288 @@
 /** Constructor for BoundaryTriangleSet.
  */
-BoundaryTriangleSet::BoundaryTriangleSet()
+BoundaryTriangleSet::BoundaryTriangleSet() :
+  Nr(-1)
 {
   for (int i = 0; i < 3; i++)
@@ -293 +296 @@
       lines[i] = NULL;
     }
-  Nr = -1;
 };
 
@@ -300 +302 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
+  Nr(number)
 {
   // set number
-  Nr = number;
   // set lines
-  Log() << Verbose(5) << "New triangle " << Nr << ":" << endl;
+  Log() << Verbose(5) << "New triangle " << Nr << ":";
   for (int i = 0; i < 3; i++)
     {
@@ -579 +581 @@
 PointCloud::PointCloud()
 {
-
 };
 
@@ -586 +587 @@
 PointCloud::~PointCloud()
 {
-
 };
 
@@ -593 +593 @@
 /** Constructor of class CandidateForTesselation.
  */
-CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) {
-  point = candidate;
-  BaseLine = line;
+CandidateForTesselation::CandidateForTesselation (BoundaryLineSet* line) :
+  point(NULL),
+  BaseLine(line),
+  ShortestAngle(2.*M_PI),
+  OtherShortestAngle(2.*M_PI)
+{
+};
+
+
+/** Constructor of class CandidateForTesselation.
+ */
+CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
+    point(candidate),
+    BaseLine(line),
+    ShortestAngle(2.*M_PI),
+    OtherShortestAngle(2.*M_PI)
+{
   OptCenter.CopyVector(&OptCandidateCenter);
   OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
@@ -607 +621 @@
 };
 
+/** output operator for CandidateForTesselation.
+ * \param &ost output stream
+ * \param &a boundary line
+ */
+ostream & operator <<(ostream &ost, const  CandidateForTesselation &a)
+{
+  ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with ";
+  if (a.point == NULL)
+    ost << "no candidate.";
+  else
+    ost << "candidate " << *(a.point) << " at angle " << (a.ShortestAngle)<< ".";
+
+  return ost;
+};
+
+
 // =========================================================== class TESSELATION ===========================================
 
 /** Constructor of class Tesselation.
  */
-Tesselation::Tesselation()
-{
-  PointsOnBoundaryCount = 0;
-  LinesOnBoundaryCount = 0;
-  TrianglesOnBoundaryCount = 0;
-  InternalPointer = PointsOnBoundary.begin();
-  LastTriangle = NULL;
-  TriangleFilesWritten = 0;
+Tesselation::Tesselation() :
+  PointsOnBoundaryCount(0),
+  LinesOnBoundaryCount(0),
+  TrianglesOnBoundaryCount(0),
+  LastTriangle(NULL),
+  TriangleFilesWritten(0),
+  InternalPointer(PointsOnBoundary.begin())
+{
 }
 ;
@@ -1283 +1313 @@
         BLS[n] = FindLine->second;
 
+        // remove existing line from OpenLines
+        CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
+        delete(CandidateLine->second);
+        OpenLines.erase(CandidateLine);
+
         break;
       }
@@ -1304 +1339 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
 {
-  Log() << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
+  Log() << Verbose(4) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
@@ -1312 +1347 @@
   // increase counter
   LinesOnBoundaryCount++;
+  // also add to open lines
+  CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
+  OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
 };
 
@@ -1680 +1718 @@
       << *(*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 != OptCandidates->end()--) {
-        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
-        SecondPoint = (*it)->point;
-        // adding point 1 and point 2 and the line between them
-        AddTesselationPoint(FirstPoint, 0);
-        AddTesselationPoint(SecondPoint, 1);
-        AddTesselationLine(TPS[0], TPS[1], 0);
-      }
+//      // if we do not reach the end with the next step of iteration, we need to setup a new first line
+//      if (it != OptCandidates->end()--) {
+//        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
+//        SecondPoint = (*it)->point;
+//        // adding point 1 and point 2 and the line between them
+//        AddTesselationPoint(FirstPoint, 0);
+//        AddTesselationPoint(SecondPoint, 1);
+//        AddTesselationLine(TPS[0], TPS[1], 0);
+//      }
       Log() << Verbose(2) << "Projection is " << BTS->NormalVector.ScalarProduct(&Oben) << "." << endl;
     }
@@ -1860 +1898 @@
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out output stream for debugging
- * @param Line current baseline to search from
+ * @param CandidateLine current cadndiate baseline to search from
  * @param T current triangle which \a Line is edge of
  * @param RADIUS radius of the rolling ball
@@ -1866 +1904 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
 {
   Log() << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
@@ -1882 +1920 @@
   double radius, CircleRadius;
 
-  Log() << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+  Log() << Verbose(1) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl;
   for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.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;
 
   // construct center of circle
-  CircleCenter.CopyVector(Line.endpoints[0]->node->node);
-  CircleCenter.AddVector(Line.endpoints[1]->node->node);
+  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+  CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node);
   CircleCenter.Scale(0.5);
 
   // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(Line.endpoints[0]->node->node);
-  CirclePlaneNormal.SubtractVector(Line.endpoints[1]->node->node);
+  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // calculate squared radius of circle
@@ -1906 +1944 @@
     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 = Line.endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
+    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
     helper.Scale(sqrt(RADIUS*RADIUS - radius));
     OldSphereCenter.AddVector(&helper);
@@ -1914 +1952 @@
     // construct SearchDirection
     SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(Line.endpoints[0]->node->node);
+    helper.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
     helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
@@ -1927 +1965 @@
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, OptCandidates, &ShortestAngle, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine.BaseLine, ThirdNode, OptCandidates, &ShortestAngle, RADIUS, LC);
 
   } else {
-    Log() << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+    Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl;
   }
 
@@ -1942 +1980 @@
   }
 
-  BoundaryLineSet *BaseRay = &Line;
+  BoundaryLineSet *BaseRay = CandidateLine.BaseLine;
   for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
     Log() << Verbose(1) << " Third point candidate is " << *(*it)->point
@@ -1969 +2007 @@
 
         if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
-          AddTesselationLine(TPS[0], TPS[1], 0);
-          AddTesselationLine(TPS[0], TPS[2], 1);
-          AddTesselationLine(TPS[1], TPS[2], 2);
-
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTesselationTriangle();
-          (*it)->OptCenter.Scale(-1.);
-          BTS->GetNormalVector((*it)->OptCenter);
-          (*it)->OptCenter.Scale(-1.);
-
-          Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-            << " for this triangle ... " << endl;
-        //Log() << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
+          CandidateLine.point = (*it)->point;
+          CandidateLine.OptCenter.CopyVector(&((*it)->OptCenter));
+          CandidateLine.OtherOptCenter.CopyVector(&((*it)->OtherOptCenter));
+          CandidateLine.ShortestAngle = ShortestAngle;
         } else {
-          eLog() << Verbose(2) << "This triangle consisting of ";
-          Log() << Verbose(0) << *(*it)->point << ", ";
-          Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
-          Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
-          Log() << Verbose(0) << "exists and is not added, as it does not seem helpful!" << endl;
+//          eLog() << Verbose(2) << "This triangle consisting of ";
+//          Log() << Verbose(0) << *(*it)->point << ", ";
+//          Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
+//          Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
+//          Log() << Verbose(0) << "exists and is not added, as it 0x80000000006fc150(does not seem helpful!" << endl;
           result = false;
         }
@@ -1997 +2026 @@
           // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
           // i.e. at least one of the three lines must be present with TriangleCount <= 1
-          if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) {
-            AddTesselationLine(TPS[0], TPS[1], 0);
-            AddTesselationLine(TPS[0], TPS[2], 1);
-            AddTesselationLine(TPS[1], TPS[2], 2);
-
-            BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-            AddTesselationTriangle();  // add to global map
-
-            (*it)->OtherOptCenter.Scale(-1.);
-            BTS->GetNormalVector((*it)->OtherOptCenter);
-            (*it)->OtherOptCenter.Scale(-1.);
-
-            eLog() << Verbose(2) << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-            Log() << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl;
+          if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS) || CandidateLine.BaseLine->skipped) {
+            CandidateLine.point = (*it)->point;
+            CandidateLine.OptCenter.CopyVector(&(*it)->OptCenter);
+            CandidateLine.OtherOptCenter.CopyVector(&(*it)->OtherOptCenter);
+            CandidateLine.ShortestAngle = ShortestAngle+2.*M_PI;
+
           } else {
-            eLog() << Verbose(2) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl;
+//            eLog() << Verbose(2) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl;
             result = false;
           }
       } else {
-        Log() << Verbose(1) << "This triangle consisting of ";
-        Log() << Verbose(0) << *(*it)->point << ", ";
-        Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
-        Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
-        Log() << Verbose(0) << "is invalid!" << endl;
+//        Log() << Verbose(1) << "This triangle consisting of ";
+//        Log() << Verbose(0) << *(*it)->point << ", ";
+//        Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
+//        Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
+//        Log() << Verbose(0) << "is invalid!" << endl;
         result = false;
       }
@@ -2043 +2064 @@
   Log() << Verbose(0) << "End of FindNextSuitableTriangle\n";
   return result;
+};
+
+/** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.
+ * \param &CandidateLine triangle to add
+ */
+void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine)
+{
+  Vector Center;
+  Log() << Verbose(2) << "BaseLine is :" << *(CandidateLine.BaseLine) << " with candidate " << *(CandidateLine.point) << "." << endl;
+  // add the points
+  AddTesselationPoint(CandidateLine.point, 0);
+  AddTesselationPoint(CandidateLine.BaseLine->endpoints[0]->node, 1);
+  AddTesselationPoint(CandidateLine.BaseLine->endpoints[1]->node, 2);
+
+  Center.CopyVector(&CandidateLine.OptCenter);
+  // add the lines
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+
+  // add the triangles
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTesselationTriangle();
+  Center.Scale(-1.);
+  BTS->GetNormalVector(Center);
+
+  Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
 };