Changeset 856098 for src

Timestamp:

Nov 28, 2009, 3:45:59 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, fad93c

Parents:

262bae

git-author:

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

git-committer:

Frederik Heber <heber@…> (11/28/09 15:45:59)

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, we did not check whether polygon to test was really planar (thx UnitTest)
  • 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;
+  UniquePolygonSet TestedPolygons;
   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();
+
+    /// 2c1. Check whether all triangles have (besides sign) same normal vector)
+    TriangleSet::const_iterator Runner = T->begin();
+    Vector NormalVector;
+    NormalVector.CopyVector(&(*Runner)->NormalVector);
+    Runner++;
+    bool IsPlanar = true;
+    for (; Runner != T->end(); Runner++) {
+      if (fabs(fabs(NormalVector.ScalarProduct(&(*Runner)->NormalVector)) - 1.) > MYEPSILON)  // check whether all NormalVectors are parallel
+        IsPlanar = false;
+    }
+
+    /// 2c2. 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 (IsPlanar) {
+      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);
+      }
+    } else {
+      Log() << Verbose(0) << "REJECT: This polygon does not lay on a plane." << endl;
       delete(Four);
+    }
 
     delete(T);
@@ -4058 +4199 @@
 
   /// 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 +4296 @@
   }
 
-  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();