Changeset 474961 for src

Timestamp:

Apr 16, 2010, 12:26:16 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:

d5fea7

Parents:

ebb50e

Message:

New approach to degenerated triangles: Recognize on creation and add both sides at once.

• FindNonConvexBorder(): Approach with constructing a tree of smallest ShortestAngle has bee abandoned.

• new function BoundaryTriangleSet::CheckDegeneracy(): Similar to CandidateForTesselation::CheckValidity()
  • checks by looking whether points are inside the sphere at OtherOptCenter
• new function AddCandidateTriangle(): refactored out of AddCandidatePolygon()
• new function AddDegenerateTriangle(): basically just checks which lines to take, newly created or existing ones and which
• AddCandidatePolygon(): (renamed from AddCandidateTriangle())
  • Now needs radius of sphere and LinkedCell structure (for the degenerated searches), adapted where called
  • calls AddCandidateTriangle() for one triangle in the polygon, subsequently calls CheckDegeneracy()
  • if appropriate, calls AddDegenerateTriangle() for the other side triangle
• rename AlwaysAddTesselationTriangleLine() -> AddNewTesselationTriangleLine()
• new function AddExistingTesselationTriangleLine(), refactored out of AddTriangleLine()

Location:

src

Files:

• boundary.cpp (modified) (1 diff)
• tesselation.cpp (modified) (11 diffs)
• tesselation.hpp (modified) (2 diffs)
• unittests/tesselation_insideoutsideunittest.cpp (modified) (1 diff)
• unittests/tesselationunittest.cpp (modified) (1 diff)

src/boundary.cpp

@@ -1014 +1014 @@
     double ShortestAngle = 4.*M_PI;
     for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
-      if (Runner->second->T->AngleFromTop < Runner->second->ShortestAngle) {
-        Log() << Verbose(1) << "INFO: Checking candidate is " << *Runner->second->BaseLine << " with point " << *(*Runner->second->pointlist.begin()) << " and angle " << Runner->second->ShortestAngle << endl;
-        if (Runner->second->ShortestAngle < ShortestAngle) {
-          baseline = Runner->second;
-          ShortestAngle = baseline->ShortestAngle;
-          Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl;
-        }
-      } else {
-        Log() << Verbose(1) << "Rejecting openline " << *Runner->second->BaseLine << " with point " << *(*Runner->second->pointlist.begin()) << " and angle " << Runner->second->ShortestAngle << " against angle from top " << Runner->second->T->AngleFromTop << "." << endl;
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+        Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl;
       }
     }
-    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
-      for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
-        if (Runner->second->ShortestAngle < ShortestAngle) {
-          baseline = Runner->second;
-          ShortestAngle = baseline->ShortestAngle;
-          Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl;
-        }
-      }
     // 2e. if we found one, add candidate
     if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
       OneLoopWithoutSuccessFlag = false;
     else {
-      TesselStruct->AddCandidatePolygon(*baseline);
+      TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
     }
 

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
@@ -1104 +1140 @@
   return flag;
 };
-
-
 
 /** output operator for CandidateForTesselation.
@@ -1811 +1845 @@
       if (FindLine->second->triangles.size() < 2) {
         insertNewLine = false;
-        Log() << Verbose(0) << "Using existing line " << *FindLine->second << endl;
-
-        BPS[0] = FindLine->second->endpoints[0];
-        BPS[1] = FindLine->second->endpoints[1];
-        BLS[n] = FindLine->second;
-
-        // remove existing line from OpenLines
-        CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
-        if (CandidateLine != OpenLines.end()) {
-          Log() << Verbose(1) << " Removing line from OpenLines." << endl;
-          delete(CandidateLine->second);
-          OpenLines.erase(CandidateLine);
-        } else {
-          DoeLog(1) && (eLog()<< Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
-        }
-
+        AddExistingTesselationTriangleLine(FindLine->second, n);
         break;
       }
@@ -1833 +1852 @@
 
   if (insertNewLine) {
-    AlwaysAddTesselationTriangleLine(a, b, n);
+    AddNewTesselationTriangleLine(a, b, n);
   }
 }
@@ -1846 +1865 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+void Tesselation::AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
 {
         Info FunctionInfo(__func__);
@@ -1860 +1879 @@
   CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
   OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
+};
+
+/** Uses an existing line for a new triangle.
+ * Sets Tesselation::BLS[\a n] and removes the lines from Tesselation::OpenLines.
+ * \param *FindLine the line to add
+ * \param n index of the line to set in Tesselation::BLS
+ */
+void Tesselation::AddExistingTesselationTriangleLine(class BoundaryLineSet *Line, int n)
+{
+  Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Using existing line " << *Line << endl;
+
+  // set endpoints and line
+  BPS[0] = Line->endpoints[0];
+  BPS[1] = Line->endpoints[1];
+  BLS[n] = Line;
+
+  // remove existing line from OpenLines
+  CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
+  if (CandidateLine != OpenLines.end()) {
+    Log() << Verbose(1) << " Removing line from OpenLines." << endl;
+    delete(CandidateLine->second);
+    OpenLines.erase(CandidateLine);
+  } else {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
+  }
 };
 
@@ -2222 +2267 @@
 
     BTS = NULL;
-    AddCandidatePolygon(OptCandidates);
+    AddCandidatePolygon(OptCandidates, RADIUS, LC);
 //    delete(BaseLine.endpoints[0]);
 //    delete(BaseLine.endpoints[1]);
@@ -2544 +2589 @@
 /** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.
  * \param CandidateLine triangle to add
- * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in AddTesselationLine()
- */
-void Tesselation::AddCandidatePolygon(CandidateForTesselation CandidateLine)
+ * \param RADIUS Radius of sphere
+ * \param *LC LinkedCell structure
+ * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in
+ * AddTesselationLine() in AddCandidateTriangle()
+ */
+void Tesselation::AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC)
 {
   Info FunctionInfo(__func__);
@@ -2570 +2618 @@
   Sprinter++;
   while(Sprinter != connectedClosestPoints->end()) {
-    // add the points
-    AddTesselationPoint(TurningPoint, 0);
-    AddTesselationPoint((*Runner), 1);
-    AddTesselationPoint((*Sprinter), 2);
-
-    // 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();
-    BTS->GetCenter(&Center);
+    // 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
@@ -2589 +2627 @@
     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;
+      }
+    }
+
     Runner = Sprinter;
     Sprinter++;
@@ -2604 +2672 @@
   }
   delete(connectedClosestPoints);
+};
+
+/** If a given \a *triangle is degenerated, this adds the "other side".
+ * 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
+ * \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)
+{
+  Info FunctionInfo(__func__);
+  TesselPoint *ThirdNode = NULL;
+  BoundaryTriangleSet *otherTriangle = NULL;
+  /// 1. For each of the open sides of \a *triangle find best candidate
+  for (int i=0;i<3;i++) {
+    BLS[i] = NULL;
+    // go through all possible lines (also degenerate ones)
+    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) {
+        // 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;
+          }
+        }
+      }
+    }
+    /// 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 {
+    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 with opposite NormalVector and other center
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    AddTesselationTriangle();
+  }
+};
+
+/** 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)
+{
+  Info FunctionInfo(__func__);
+  // add the points
+  AddTesselationPoint(first, 0);
+  AddTesselationPoint(second, 1);
+  AddTesselationPoint(third, 2);
+
+  // 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();
 };
 

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 +284 @@
     void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
     void AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
-    void AlwaysAddTesselationTriangleLine(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 AddCandidatePolygon(CandidateForTesselation CandidateLine);
+    void AddCandidateTriangle(TesselPoint *first, TesselPoint *second, TesselPoint *third);
+    void AddDegeneratedTriangle(BoundaryTriangleSet *&triangle, Vector *OtherOptCenter, 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);

src/unittests/tesselation_insideoutsideunittest.cpp

@@ -123 +123 @@
       OneLoopWithoutSuccessFlag = false;
     else {
-      TesselStruct->AddCandidatePolygon(*baseline);
+      TesselStruct->AddCandidatePolygon(*baseline, SPHERERADIUS, LinkedList);
     }
   }

src/unittests/tesselationunittest.cpp

@@ -91 +91 @@
       OneLoopWithoutSuccessFlag = false;
     else {
-      TesselStruct->AddCandidatePolygon(*baseline);
+      TesselStruct->AddCandidatePolygon(*baseline, SPHERERADIUS, LinkedList);
     }
   }