Changeset af374d

Timestamp:

Nov 28, 2009, 3:07:42 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:

2130dd

Parents:

262bae

git-author:

Frederik Heber <heber@…> (11/28/09 14:58:40)

git-committer:

Frederik Heber <heber@…> (11/28/09 15:07:42)

Message:

Correction of degenerated polygons is working.

• Correction of degenerated polygons:
  • incorporated into FindNonConvexBorder()
  • verbosity added: BoundaryPolygonSet::GetCenter(), BoundaryPolygonSet::ContainsBoundaryLine(), BoundaryPolygonSet::ContainsBoundaryPoint(), BoundaryPolygonSet::ContainsPresentTupel(), Tesselation::AddTesselationLine(), Tesselation::CorrectAllDegeneratedPolygons(), ArePolygonsEdgeConnected()
  • renamed function BoundaryPolygonSet::GetAllTrianglesFromEndpoints() to GetAllContainedTrianglesFromEndpoints() (as it only returns triangles whose endpoints are _all_ in the BoundaryPolygonset)
  • BoundaryPolygonSet::FillPolygonFromTrianglesOfLine() does not strictly return four or NULL endpoints, as it makes subsequent success checks (and verbose output of failure) harder.
  • FIX: Tesselation::RemoveTesselationTriangle() adds line to Tesselation::OpenLines of still attached to other triangle.
  • new struct BoundaryLineSetCompare(), that compares BoundaryLineSet's
  • new struct BoundaryPolygonSetCompare(), that compares BoundaryPolygonSet's
  • FIX: lots of small fixes to Tesselation::CorrectAllDegeneratedPolygons(): new backside triangles were not added to LinesOnBoundary, when joining adjacent polygons triangles were gathered from Four instead of current PolygonRunner, other combination was not added to SimplyDegeneratedTriangles
  • FIX: CountTrianglePairContainingPolygon() - PairTrianglesnodes instead of TriangleNodes was clear()'ed.

• BUGFIX: FindAllDegeneratedLines() did not find the correct lines due to wrong comparison. Now we use a set with the comparator BoundaryLineSetCompare().
• valgrind runs through fine.

Signed-off-by: Frederik Heber <heber@…>

Location:

src

Files:

• boundary.cpp (modified) (1 diff)
• tesselation.cpp (modified) (20 diffs)
• tesselation.hpp (modified) (2 diffs)
• tesselationhelpers.cpp (modified) (3 diffs)

src/boundary.cpp

@@ -1041 +1041 @@
 //  TesselStruct->RemoveDegeneratedTriangles();
 
+  // correct degenerated polygons
+  TesselStruct->CorrectAllDegeneratedPolygons();
+
   // check envelope for consistency
   status = CheckListOfBaselines(TesselStruct);

src/tesselation.cpp

@@ -636 +636 @@
   }
   center->Scale(1./(double)counter);
+  Log() << Verbose(1) << "Center is at " << *center << "." << endl;
 }
 
@@ -654 +655 @@
 bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
 {
+  Info FunctionInfo(__func__);
   return ContainsPresentTupel(line->endpoints, 2);
 };
@@ -664 +666 @@
 {
   Info FunctionInfo(__func__);
-  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
-    if (point == (*Runner))
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+    Log() << Verbose(0) << "Checking against " << **Runner << endl;
+    if (point == (*Runner)) {
+      Log() << Verbose(0) << " Contained." << endl;
       return true;
+    }
+  }
+  Log() << Verbose(0) << " Not contained." << endl;
   return false;
 };
@@ -678 +685 @@
   Info FunctionInfo(__func__);
   for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
-    if (point == (*Runner)->node)
+    if (point == (*Runner)->node) {
+      Log() << Verbose(0) << " Contained." << endl;
       return true;
+    }
+  Log() << Verbose(0) << " Not contained." << endl;
   return false;
 };
@@ -690 +700 @@
 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const
 {
+  Info FunctionInfo(__func__);
   int counter = 0;
-  for(int i=0;i<dim;i++)
-    if (ContainsBoundaryPoint(Points[i]))
+  Log() << Verbose(1) << "Polygon is " << *this << endl;
+  for(int i=0;i<dim;i++) {
+    Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl;
+    if (ContainsBoundaryPoint(Points[i])) {
       counter++;
+    }
+  }
 
   if (counter == dim)
@@ -707 +722 @@
 bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const
 {
+  Info FunctionInfo(__func__);
   size_t counter = 0;
+  Log() << Verbose(1) << "Polygon is " << *this << endl;
   for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+    Log() << Verbose(1) << " Testing endpoint " << **Runner << endl;
     if (ContainsBoundaryPoint(*Runner))
       counter++;
@@ -725 +743 @@
 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
 {
-  Info FunctionInfo(__func__);
   return ContainsPresentTupel((const PointSet)P->endpoints);
 };
@@ -732 +749 @@
  * \return set of all triangles
  */
-TriangleSet * BoundaryPolygonSet::GetAllTrianglesFromEndpoints() const
+TriangleSet * BoundaryPolygonSet::GetAllContainedTrianglesFromEndpoints() const
 {
   Info FunctionInfo(__func__);
+  pair <TriangleSet::iterator, bool> Tester;
   TriangleSet *triangles = new TriangleSet;
 
   for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
     for(LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)
-      for(TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++)
-        triangles->insert(Sprinter->second);
+      for(TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {
+        //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl;
+        if (ContainsBoundaryTriangle(Sprinter->second)) {
+          Tester = triangles->insert(Sprinter->second);
+          if (Tester.second)
+            Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl;
+        }
+      }
 
   Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl;
@@ -748 +772 @@
 /** Fills the endpoints of this polygon from the triangles attached to \a *line.
  * \param *line lines with triangles attached
- * \return true - polygon contains four points, false - something went wront, polygon does not contain any points
+ * \return true - polygon contains endpoints, false - line was NULL
  */
 bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)
 {
+  Info FunctionInfo(__func__);
+  pair <PointSet::iterator, bool> Tester;
+  if (line == NULL)
+    return false;
+  Log() << Verbose(1) << "Filling polygon from line " << *line << endl;
   for(TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {
-    for (int i=0;i<3;i++)
-      endpoints.insert((Runner->second)->endpoints[i]);
-  }
-
-  if (endpoints.size() != 4)  {// remove all points on error
-    endpoints.clear();
-    return false;
-  } else
-    return true;
+    for (int i=0;i<3;i++) {
+      Tester = endpoints.insert((Runner->second)->endpoints[i]);
+      if (Tester.second)
+        Log() << Verbose(1) << "  Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl;
+    }
+  }
+
+  return true;
 };
 
@@ -1577 +1605 @@
         // remove existing line from OpenLines
         CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
-        delete(CandidateLine->second);
-        OpenLines.erase(CandidateLine);
+        if (CandidateLine != OpenLines.end()) {
+          Log() << Verbose(1) << " Removing line from OpenLines." << endl;
+          delete(CandidateLine->second);
+          OpenLines.erase(CandidateLine);
+        } else {
+          eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl;
+        }
 
         break;
@@ -1670 +1703 @@
       } else {
         Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: ";
+        OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL));
         for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
           Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
@@ -2332 +2366 @@
     Sprinter++;
   }
+  delete(connectedClosestPoints);
 };
 
@@ -3707 +3742 @@
 }
 
+struct BoundaryLineSetCompare {
+  bool operator() (const BoundaryLineSet * const a, const BoundaryLineSet * const b) {
+    int lowerNra = -1;
+    int lowerNrb = -1;
+
+    if (a->endpoints[0] < a->endpoints[1])
+      lowerNra = 0;
+    else
+      lowerNra = 1;
+
+    if (b->endpoints[0] < b->endpoints[1])
+      lowerNrb = 0;
+    else
+      lowerNrb = 1;
+
+    if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
+      return true;
+    else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
+      return false;
+    else {  // both lower-numbered endpoints are the same ...
+     if (a->endpoints[(lowerNra+1)%2] < b->endpoints[(lowerNrb+1)%2])
+       return true;
+     else if (a->endpoints[(lowerNra+1)%2] > b->endpoints[(lowerNrb+1)%2])
+       return false;
+    }
+    return false;
+  };
+};
+
+#define UniqueLines set < class BoundaryLineSet *, BoundaryLineSetCompare>
+
 /**
  * Finds all degenerated lines within the tesselation structure.
@@ -3716 +3782 @@
 {
         Info FunctionInfo(__func__);
-  map<int, class BoundaryLineSet *> AllLines;
+        UniqueLines AllLines;
   map<int, int> * DegeneratedLines = new map<int, int>;
 
@@ -3726 +3792 @@
 
   LineMap::iterator LineRunner1;
-  pair<LineMap::iterator, bool> tester;
+  pair< UniqueLines::iterator, bool> tester;
   for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
-    tester = AllLines.insert( pair<int,BoundaryLineSet *> (LineRunner1->second->endpoints[0]->Nr, LineRunner1->second) );
-    if ((!tester.second) && (tester.first->second->endpoints[1]->Nr == LineRunner1->second->endpoints[1]->Nr)) { // found degenerated line
-      DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, tester.first->second->Nr) );
-      DegeneratedLines->insert ( pair<int, int> (tester.first->second->Nr, LineRunner1->second->Nr) );
+    tester = AllLines.insert( LineRunner1->second );
+    if (!tester.second) { // found degenerated line
+      DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr) );
+      DegeneratedLines->insert ( pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr) );
     }
   }
@@ -3739 +3805 @@
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
   map<int,int>::iterator it;
-  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++)
-      Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
+  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
+    const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
+    const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second);
+    if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end())
+      Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl;
+    else
+      eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl;
+  }
 
   return DegeneratedLines;
@@ -4024 +4096 @@
 };
 
+struct BoundaryPolygonSetCompare {
+  bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const {
+    if (s1->endpoints.size() < s2->endpoints.size())
+      return true;
+    else if (s1->endpoints.size() > s2->endpoints.size())
+      return false;
+    else { // equality of number of endpoints
+      PointSet::const_iterator Walker1 = s1->endpoints.begin();
+      PointSet::const_iterator Walker2 = s2->endpoints.begin();
+      while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
+        if ((*Walker1)->Nr < (*Walker2)->Nr)
+          return true;
+        else if ((*Walker1)->Nr > (*Walker2)->Nr)
+          return false;
+        Walker1++;
+        Walker2++;
+      }
+      return false;
+    }
+  }
+};
+
+#define UniquePolygonSet set < BoundaryPolygonSet *, BoundaryPolygonSetCompare>
+
 /** Finds all degenerated polygons and calls ReTesselateDegeneratedPolygon()/
  * \return number of polygons found
@@ -4035 +4131 @@
 
   /// 2. Go through all lines not contained in degenerated triangles.
-  PolygonList ListofFours;
-  BoundaryPolygonSet *Four = NULL;
-  TriangleSet *T = NULL;
+  UniquePolygonSet ListofFours;
   for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++) {
-    Four = new BoundaryPolygonSet;
+    Log() << Verbose(1) << "Current Line is: " << *(LineRunner->second) << endl;
+
+    BoundaryPolygonSet *Four = new BoundaryPolygonSet;
 
     /// 2a. Get all four endpoints of the two connected triangles.
     Four->FillPolygonFromTrianglesOfLine((LineRunner->second));
+    if (Four->endpoints.size() != 4) {
+      Log() << Verbose(0) << "REJECT: There were no four endpoints:" << *Four << endl;
+      delete(Four);
+      continue;
+    }
 
     /// 2b. Get the triangles of all those endpoints
-    T = Four->GetAllTrianglesFromEndpoints();
-
-    /// 2c. Find all pairs of those triangles that contain the four endpoints
-    if (CountTrianglePairContainingPolygon(Four, T) > 1) {
-      /// 2d. Check the number of pairs, if greater 1, we have a degenerated polygon (mark down for later use)
-      ListofFours.push_back(Four);
-    } else
+    TriangleSet *T = Four->GetAllContainedTrianglesFromEndpoints();
+
+    /// 2c. check whether triangles are not already simply degenerated.
+    size_t DegeneratedCount = 0;
+    for (TriangleSet::const_iterator Runner = T->begin(); Runner != T->end(); Runner++) {
+      if (SimplyDegeneratedTriangles->find((*Runner)->Nr) != SimplyDegeneratedTriangles->end())
+        DegeneratedCount++;
+    }
+    Log() << Verbose(1) << "Found " << DegeneratedCount << " degenerated triangles." << endl;
+    if (DegeneratedCount == 0) {
+      /// 2c. Find all pairs of those triangles that contain the four endpoints
+      const int counter = CountTrianglePairContainingPolygon(Four, T);
+      Log() << Verbose(1) << "There are " << counter << " pairs for this polygon of four." << endl;
+      if (counter %3 == 0) {
+        /// 2d. Check the number of pairs, if greater 1, we have a degenerated polygon (mark down for later use)
+        pair < UniquePolygonSet::iterator, bool > Tester = ListofFours.insert(Four);
+        if (Tester.second) {
+          Log() << Verbose(0) << "ACCEPT: We have a degenerated polygon." << *Four << endl;
+        } else {
+          Log() << Verbose(0) << "REJECT: Degenerated polygon already present." << *Four << endl;
+          delete(Four);
+        }
+      } else {
+        Log() << Verbose(0) << "REJECT: Seems to be the edge line of a degenerated polygon only." << *Four << endl;
+        delete(Four);
+      }
+    } else {
+      Log() << Verbose(0) << "REJECT: This polygon contains some simply degenerated triangles." << endl;
       delete(Four);
+    }
 
     delete(T);
@@ -4058 +4181 @@
 
   /// 3. Combine all edge-connected degenerated polygons
-  PolygonList::iterator PolygonWalker;  // is the inner iterator
-  PolygonList::iterator PolygonSprinter; // is the inner advanced iterator
-  for (PolygonList::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++) {
-    PolygonWalker = PolygonRunner;
-    PolygonWalker++;
-    PolygonSprinter = PolygonWalker;
-    while (PolygonWalker != ListofFours.end()) {
-      PolygonSprinter++;
-      if (ArePolygonsEdgeConnected((*PolygonRunner), (*PolygonWalker))) { // if connected
-        CombinePolygons((*PolygonRunner), (*PolygonWalker)); // combined and ...
-        ListofFours.erase(PolygonWalker); // ... remove from list
+  {
+    UniquePolygonSet::iterator PolygonWalker;  // is the inner iterator
+    UniquePolygonSet::iterator PolygonSprinter; // is the inner advanced iterator
+    for (UniquePolygonSet::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++) {
+      const BoundaryPolygonSet *P1 = *PolygonRunner;
+      Log() << Verbose(1) << "Current Polygon is : " << *P1 << endl;
+      PolygonWalker = PolygonRunner;
+      PolygonWalker++;
+      PolygonSprinter = PolygonWalker;
+      while (PolygonWalker != ListofFours.end()) {
+        BoundaryPolygonSet *P2 = *PolygonWalker;
+        Log() << Verbose(1) << "Current Other Polygon is : " << *P2 << endl;
+        PolygonSprinter++;
+        Log() << Verbose(1) << "Are " << *P1 << " and " << *P2 << " edge-connected?" << endl;
+        if (ArePolygonsEdgeConnected((*PolygonRunner), (P2))) { // if connected
+          Log() << Verbose(0) << "Yes!" << endl;
+          CombinePolygons((*PolygonRunner), (P2)); // combined and ...
+          ListofFours.erase(PolygonWalker); // ... remove from list
+          Log() << Verbose(1) << "  New current Polygon is : " << *P1 << endl;
+        } else {
+          Log() << Verbose(0) << "No." << endl;
+        }
+        PolygonWalker = PolygonSprinter;
       }
-      PolygonWalker = PolygonSprinter;
-    }
-  }
+    }
+  }
+
+  Log() << Verbose(0) << "The following degenerated polygons have been found: " << endl;
+  for (UniquePolygonSet::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++)
+    Log() << Verbose(0) << " " << **PolygonRunner << endl;
 
   /// 4. Go through all these degenerated polygons
-  Vector NormalVector;
-  for (PolygonList::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++) {
+  for (UniquePolygonSet::iterator PolygonRunner = ListofFours.begin(); PolygonRunner != ListofFours.end(); PolygonRunner++) {
+    stack <int> TriangleNrs;
+    Vector NormalVector;
     /// 4a. Gather all triangles of this polygon
-    T = Four->GetAllTrianglesFromEndpoints();
-
+    TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();
+
+    TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")
-    NormalVector.CopyVector(&(*(T->begin()))->NormalVector);
-
+    NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
+    Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl;
+    TriangleWalker++;
+    TriangleSet::iterator TriangleSprinter = TriangleWalker; // is the inner advanced iterator
     /// 4b. Remove all triangles whose NormalVector is in opposite direction (i.e. "back")
-    TriangleSet::const_iterator TriangleWalker;  // is the inner iterator
-    TriangleSet::const_iterator TriangleSprinter; // is the inner advanced iterator
-    PolygonWalker = PolygonRunner;
-    PolygonWalker++;
-    PolygonSprinter = PolygonWalker;
-    while (PolygonWalker != ListofFours.end()) {
-      PolygonSprinter++;
-      if ((*TriangleWalker)->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list
-        delete(*TriangleWalker);  // remove triangle
+    BoundaryTriangleSet *triangle = NULL;
+    while (TriangleSprinter != T->end()) {
+      TriangleWalker = TriangleSprinter;
+      triangle = *TriangleWalker;
+      TriangleSprinter++;
+      Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl;
+      if (triangle->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list
+        Log() << Verbose(1) << " Removing ... " << endl;
+        TriangleNrs.push(triangle->Nr);
         T->erase(TriangleWalker);
-      }
+        RemoveTesselationTriangle(triangle);
+      } else
+        Log() << Verbose(1) << " Keeping ... " << endl;
     }
     /// 4c. Copy all "front" triangles but with inverse NormalVector
-    BoundaryTriangleSet *triangle = NULL;
     TriangleWalker = T->begin();
-    TriangleSprinter = T->end();
-    TriangleSprinter--; // set onto last "front" triangle
-    while (TriangleWalker != TriangleSprinter) {  // go through all front triangles
-      triangle = new BoundaryTriangleSet(*(*TriangleWalker)); // copy triangle ...
-      triangle->NormalVector.Scale(-1.);
-      T->insert(triangle); // ... and add
+    while (TriangleWalker != T->end()) {  // go through all front triangles
+      Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << BTS->NormalVector << endl;
+      for (int i = 0; i < 3; i++)
+        AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);
+      AddTesselationLine(TPS[0], TPS[1], 0);
+      AddTesselationLine(TPS[0], TPS[2], 1);
+      AddTesselationLine(TPS[1], TPS[2], 2);
+      BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...
+      AddTesselationTriangle(); // ... and add
+      TriangleNrs.pop();
+      BTS->NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
+      BTS->NormalVector.Scale(-1.);
       /// 4d. Add all degenerated triangles to the list of simply degenerated triangles
-      SimplyDegeneratedTriangles->insert(pair <int, int> ((*TriangleWalker)->Nr, triangle->Nr) );
+      SimplyDegeneratedTriangles->insert(pair <int, int> ((*TriangleWalker)->Nr, BTS->Nr) );
+      SimplyDegeneratedTriangles->insert(pair <int, int> (BTS->Nr, (*TriangleWalker)->Nr) );
       TriangleWalker++;
     }
+    if (!TriangleNrs.empty()) {
+      eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl;
+    }
     delete(T);  // remove the triangleset
   }
 
+  Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
+  map<int,int>::iterator it;
+  for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
+      Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
+
   /// 5. exit
-  PolygonList::iterator PolygonRunner;
-  while (ListofFours.empty()) {
+  UniquePolygonSet::iterator PolygonRunner;
+  while (!ListofFours.empty()) {
     PolygonRunner = ListofFours.begin();
     delete(*PolygonRunner);
@@ -4120 +4278 @@
   }
 
-  int counter = SimplyDegeneratedTriangles->size();
+  const int counter = SimplyDegeneratedTriangles->size();
   delete(SimplyDegeneratedTriangles);
   return counter;

src/tesselation.hpp

@@ -23 +23 @@
 #include <list>
 #include <set>
+#include <stack>
 
 #include "atom_particleinfo.hpp"
@@ -183 +184 @@
     bool ContainsPresentTupel(const BoundaryPolygonSet * const P) const;
     bool ContainsPresentTupel(const PointSet &endpoints) const;
-    TriangleSet * GetAllTrianglesFromEndpoints() const;
+    TriangleSet * GetAllContainedTrianglesFromEndpoints() const;
     bool FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line);
 

src/tesselationhelpers.cpp

@@ -956 +956 @@
   }
 
+  Log() << Verbose(0) << "Polygon is " << *P << endl;
   // create each pair, get the endpoints and check whether *P is contained.
   int counter = 0;
@@ -969 +970 @@
         for (int i=0;i<3;i++)
           PairTrianglenodes.endpoints.insert((*PairWalker)->endpoints[i]);
-        Log() << Verbose(1) << "Current pair of triangles: " << PairTrianglenodes << "." << endl;
-        // now check
-        if (PairTrianglenodes.ContainsPresentTupel(P)) {
-          counter++;
-          Log() << Verbose(1) << "Matches with " << *P << endl;
+        const int size = PairTrianglenodes.endpoints.size();
+        if (size == 4) {
+          Log() << Verbose(0) << " Current pair of triangles: " << **Walker << "," << **PairWalker << " with " << size << " distinct endpoints:" << PairTrianglenodes << endl;
+          // now check
+          if (PairTrianglenodes.ContainsPresentTupel(P)) {
+            counter++;
+            Log() << Verbose(0) << "  ACCEPT: Matches with " << *P << endl;
+          } else {
+            Log() << Verbose(0) << "  REJECT: No match with " << *P << endl;
+          }
         } else {
-          Log() << Verbose(1) << "No match with " << *P << endl;
+          Log() << Verbose(0) << "  REJECT: Less than four endpoints." << endl;
         }
       }
     }
-    PairTrianglenodes.endpoints.clear();
+    Trianglenodes.clear();
   }
   return counter;
@@ -1010 +1016 @@
 {
   Info FunctionInfo(__func__);
-  for(PointSet::iterator Runner = P1->endpoints.begin(); Runner != P1->endpoints.end(); Runner++) {
-    P1->endpoints.insert((*Runner));
+  pair <PointSet::iterator, bool> Tester;
+  for(PointSet::iterator Runner = P2->endpoints.begin(); Runner != P2->endpoints.end(); Runner++) {
+    Tester = P1->endpoints.insert((*Runner));
+    if (Tester.second)
+      Log() << Verbose(0) << "Inserting endpoint " << *(*Runner) << " into first polygon." << endl;
   }
   P2->endpoints.clear();