Changeset 57066a

Timestamp:

Sep 28, 2009, 8:00:33 PM (16 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:

e4a379

Parents:

91e7e4a

git-author:

Frederik Heber <heber@…> (09/28/09 18:47:54)

git-committer:

Frederik Heber <heber@…> (09/28/09 20:00:33)

Message:

Various fixes and attempt to get convex hull working.

• Tesselation::FlipBaseline: Now returns pointer to new baseline or NULL if not flipped
• Tesselation::PickFarthestofTwoBaselines() returns volume change due to flipping
• WriteTecplotFile(): ZONE name now contains nodes of last triangle, incorporates upper two changes
• ConvexizeNonconvexEnvelope(): additional dummy string to add single step output of intermediate envelopes
• TesselStruct::RemovePointFromTesselatedSurface(): Now returns volume lost due to point removal
• ConvexizeNonconvexEnvelope() evaluates change in volume due to added/removed points
• FindNonConvexBorder(): last two parameters have changed position, we check for NULL filename before writing, we add straddling points by InsertStraddlingPoints() (is better than AddBoundaryPointByDegeneratedTriangle()), IsInnerAtom() test removed
• new function AddBoundaryPointByDegeneratedTriangle() to add outside points
• StoreTrianglesinFile() contains the triangle storing
• boundary.cpp: typos, Verbose() changes and Verbosity reduced
• new function BoundaryTriangleSet::IsPresentTupel(), checks whether two triangles are the same, used in Tesselation::FindTriangles()
• new attribute Tesselation::LastTriangle contains pointer to last added triangle
• FIX: Tesselation::FindNextSuitableTriangle() checks for existentTriangleCount between 1 and 3 (as is the number of already present lines) and checks for zero NormalVector of created triangle
• Tesselation::FindSecondPointForTesselation(): removed parameter Candidate
• Tesselation::FindClosestTrianglesToPoint() now has a list of triangles from the nearest point which we look through for nearest triangle
• BUGFIX: Tesselation::IsInnerPoint() now takes center of triangle with respect to point to decide whether in- or outside
• FIX: Tesselation::GetPathsOfConnectedPoints() notes down touched triangles as well in order move along the surrounding triangles sequentially
• FIX: Tesselation::RemovePointFromTesselatedSurface(): Huge changes, from the closed path of connected points we remove one by one by adding the respective triangle and try to do it in such an order, that the volume is maximized

• new function Tesselation::FindAllDegeneratedTriangles(), pendant to Tesselation::FindAllDegeneratedLines()
• FIX: new function Tesselation::RemoveDegeneratedTriangles()
• new function Tesselation::AddBoundaryPointByDegeneratedTriangle()

Moved tesselation writing to tesselation.cpp

• WriteVrmlFile(), WriteRaster3dFile(), WriteTecplotFile() moved to tesselation.cpp
• these also don't use molecule anymore, but the PointCloud structure (hence, bonds are no more visible)
• new function TesselStruct::Output() which performs the writing
• new function InsertSphereInRaster3D(), which places a sphere at the position of the last triangle

Location:

src

Files:

• boundary.cpp (modified) (15 diffs)
• boundary.hpp (modified) (2 diffs)
• tesselation.cpp (modified) (59 diffs)
• tesselation.hpp (modified) (8 diffs)
• tesselationhelpers.cpp (modified) (2 diffs)
• tesselationhelpers.hpp (modified) (2 diffs)

src/boundary.cpp

@@ -112 +112 @@
 ;
 
-/** Creates the objects in a VRML file.
- * \param *out output stream for debugging
- * \param *vrmlfile output stream for tecplot data
- * \param *Tess Tesselation structure with constructed triangles
- * \param *mol molecule structure with atom positions
- */
-void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
-{
-  atom *Walker = mol->start;
-  bond *Binder = mol->first;
-  int i;
-  Vector *center = mol->DetermineCenterOfAll(out);
-  if (vrmlfile != NULL) {
-    //cout << Verbose(1) << "Writing Raster3D file ... ";
-    *vrmlfile << "#VRML V2.0 utf8" << endl;
-    *vrmlfile << "#Created by molecuilder" << endl;
-    *vrmlfile << "#All atoms as spheres" << endl;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
-      for (i=0;i<NDIM;i++)
-        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
-      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
-    }
-
-    *vrmlfile << "# All bonds as vertices" << endl;
-    while (Binder->next != mol->last) {
-      Binder = Binder->next;
-      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
-      for (i=0;i<NDIM;i++)
-        *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
-      *vrmlfile << "\t0.03\t";
-      for (i=0;i<NDIM;i++)
-        *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
-      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
-    }
-
-    *vrmlfile << "# All tesselation triangles" << endl;
-    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-      *vrmlfile << "1" << endl << "  "; // 1 is triangle type
-      for (i=0;i<3;i++) { // print each node
-        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
-        *vrmlfile << "\t";
-      }
-      *vrmlfile << "1. 0. 0." << endl;  // red as colour
-      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
-    }
-  } else {
-    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
-  }
-  delete(center);
-};
-
-/** Creates the objects in a raster3d file (renderable with a header.r3d).
- * \param *out output stream for debugging
- * \param *rasterfile output stream for tecplot data
- * \param *Tess Tesselation structure with constructed triangles
- * \param *mol molecule structure with atom positions
- */
-void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
-{
-  atom *Walker = mol->start;
-  bond *Binder = mol->first;
-  int i;
-  Vector *center = mol->DetermineCenterOfAll(out);
-  if (rasterfile != NULL) {
-    //cout << Verbose(1) << "Writing Raster3D file ... ";
-    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
-    *rasterfile << "@header.r3d" << endl;
-    *rasterfile << "# All atoms as spheres" << endl;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
-      for (i=0;i<NDIM;i++)
-        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
-      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
-    }
-
-    *rasterfile << "# All bonds as vertices" << endl;
-    while (Binder->next != mol->last) {
-      Binder = Binder->next;
-      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
-      for (i=0;i<NDIM;i++)
-        *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
-      *rasterfile << "\t0.03\t";
-      for (i=0;i<NDIM;i++)
-        *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
-      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
-    }
-
-    *rasterfile << "# All tesselation triangles" << endl;
-    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
-    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
-      for (i=0;i<3;i++) {  // print each node
-        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
-        *rasterfile << "\t";
-      }
-      *rasterfile << "1. 0. 0." << endl;  // red as colour
-      //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
-    }
-    *rasterfile << "9\n#  terminating special property\n";
-  } else {
-    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
-  }
-  delete(center);
-};
-
-/** This function creates the tecplot file, displaying the tesselation of the hull.
- * \param *out output stream for debugging
- * \param *tecplot output stream for tecplot data
- * \param N arbitrary number to differentiate various zones in the tecplot format
- */
-void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
-{
-  if ((tecplot != NULL) && (TesselStruct != NULL)) {
-    // write header
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int *LookupList = new int[mol->AtomCount];
-    for (int i = 0; i < mol->AtomCount; i++)
-      LookupList[i] = -1;
-
-    // print atom coordinates
-    *out << Verbose(2) << "The following triangles were created:";
-    int Counter = 1;
-    TesselPoint *Walker = NULL;
-    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
-      Walker = target->second->node;
-      LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
-    }
-    *tecplot << endl;
-    // print connectivity
-    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
-      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
-    }
-    delete[] (LookupList);
-    *out << endl;
-  }
-}
-
 
 /** Determines the boundary points of a cluster.
@@ -536 +390 @@
 
         // flip the line
-        if (!mol->TesselStruct->PickFarthestofTwoBaselines(out, line))
+        if (mol->TesselStruct->PickFarthestofTwoBaselines(out, line) == 0.)
           *out << Verbose(1) << "ERROR: Correction of concave baselines failed!" << endl;
-        else
+        else {
+          mol->TesselStruct->FlipBaseline(out, line);
           *out << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
+        }
       }
     }
@@ -644 +500 @@
   class BoundaryLineSet *line = NULL;
   bool Concavity;
+  char dummy[MAXSTRINGSIZE];
   PointMap::iterator PointRunner, PointAdvance;
   LineMap::iterator LineRunner, LineAdvance;
@@ -656 +513 @@
   }
 
-  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-  StoreTrianglesinFile(out, mol, filename, "-first");
-
   // First step: RemovePointFromTesselatedSurface
+  int run = 0;
+  double tmp;
   do {
     Concavity = false;
+    sprintf(dummy, "-first-%d", run);
+    //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+    StoreTrianglesinFile(out, mol, filename, dummy);
+
     PointRunner = TesselStruct->PointsOnBoundary.begin();
     PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
@@ -671 +531 @@
         line = LineRunner->second;
         *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      }
-      if (!line->CheckConvexityCriterion(out)) {
-        *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
-        // remove the point
-        Concavity = true;
-        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+        if (!line->CheckConvexityCriterion(out)) {
+          // remove the point if needed
+          *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
+          volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
+          sprintf(dummy, "-first-%d", ++run);
+          StoreTrianglesinFile(out, mol, filename, dummy);
+          Concavity = true;
+          break;
+        }
       }
       PointRunner = PointAdvance;
     }
 
+    sprintf(dummy, "-second-%d", run);
     //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-    //StoreTrianglesinFile(out, mol, filename, "-second");
+    StoreTrianglesinFile(out, mol, filename, dummy);
 
     // second step: PickFarthestofTwoBaselines
@@ -693 +557 @@
       // take highest of both lines
       if (TesselStruct->IsConvexRectangle(out, line) == NULL) {
-        TesselStruct->PickFarthestofTwoBaselines(out, line);
-        Concavity = true;
+        tmp = TesselStruct->PickFarthestofTwoBaselines(out, line);
+        volume += tmp;
+        if (tmp != 0) {
+          mol->TesselStruct->FlipBaseline(out, line);
+          Concavity = true;
+        }
       }
       LineRunner = LineAdvance;
     }
+    run++;
   } while (Concavity);
-  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-  StoreTrianglesinFile(out, mol, filename, "-third");
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  //StoreTrianglesinFile(out, mol, filename, "-third");
 
   // third step: IsConvexRectangle
@@ -717 +586 @@
       point = TesselStruct->IsConvexRectangle(out, line);
       if (point != NULL)
-        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+        volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
     }
     LineRunner = LineAdvance;
@@ -723 +592 @@
 
   CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-  StoreTrianglesinFile(out, mol, filename, "-fourth");
+  StoreTrianglesinFile(out, mol, filename, "");
 
   // end
+  *out << Verbose(1) << "Volume is " << volume << "." << endl;
   *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
   return volume;
@@ -993 +863 @@
     if ((*ListRunner)->TesselStruct == NULL) {
       *out << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
-      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], NULL, 5.);
+      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], 5., NULL);
     }
     i++;
@@ -1115 +985 @@
  * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
  * \param *LCList atoms in LinkedCell list
+ * \param RADIUS radius of the virtual sphere
  * \param *filename filename prefix for output of vertex data
- * \para RADIUS radius of the virtual sphere
- */
-void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename, const double RADIUS)
+ */
+void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const double RADIUS, const char *filename = NULL)
 {
   int N = 0;
   bool freeLC = false;
-  ofstream *tempstream = NULL;
-  char NumberName[255];
-  int TriangleFilesWritten = 0;
 
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
@@ -1141 +1008 @@
   bool failflag = false;
 
+  // initialise Linked Cell
   if (LCList == NULL) {
     LCList = new LinkedCell(mol, 2.*RADIUS);
@@ -1146 +1014 @@
   }
 
+  // 1. get starting triangle
   mol->TesselStruct->FindStartingTriangle(out, RADIUS, LCList);
 
+  // 2. expand from there
   baseline = mol->TesselStruct->LinesOnBoundary.begin();
-  // the outward most line is dangerous, as we may end up with wrapping up the starting triangle, hence
-  // terminating the algorithm too early.
-  if (baseline != mol->TesselStruct->LinesOnBoundary.end()) // skip first line as it its the outwardmost!
-        baseline++;
+  baseline++; // skip first line
   while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (flag)) {
     if (baseline->second->triangles.size() == 1) {
-      failflag = mol->TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
+      // 3. find next triangle
+      failflag = mol->TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
         cerr << "WARNING: FindNextSuitableTriangle failed." << endl;
+
       // write temporary envelope
-      if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % SingleStepWidth == 0))) { // if we have a new triangle and want to output each new triangle configuration
-        TriangleMap::iterator runner = mol->TesselStruct->TrianglesOnBoundary.end();
-        runner--;
-        class BoundaryTriangleSet *triangle = runner->second;
-        if (triangle != NULL) {
-          sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
-          if (DoTecplotOutput) {
-            string NameofTempFile(filename);
-            NameofTempFile.append(NumberName);
-            for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-            NameofTempFile.erase(npos, 1);
-            NameofTempFile.append(TecplotSuffix);
-            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-            WriteTecplotFile(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
-            tempstream->close();
-            tempstream->flush();
-            delete(tempstream);
-          }
-
-          if (DoRaster3DOutput) {
-            string NameofTempFile(filename);
-            NameofTempFile.append(NumberName);
-            for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-            NameofTempFile.erase(npos, 1);
-            NameofTempFile.append(Raster3DSuffix);
-            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-            WriteRaster3dFile(out, tempstream, mol->TesselStruct, mol);
-    //        // include the current position of the virtual sphere in the temporary raster3d file
-    //        // make the circumsphere's center absolute again
-    //        helper.CopyVector(BaseRay->endpoints[0]->node->node);
-    //        helper.AddVector(BaseRay->endpoints[1]->node->node);
-    //        helper.Scale(0.5);
-    //        (*it)->OptCenter.AddVector(&helper);
-    //        Vector *center = mol->DetermineCenterOfAll(out);
-    //        (*it)->OptCenter.SubtractVector(center);
-    //        delete(center);
-    //        // and add to file plus translucency object
-    //        *tempstream << "# current virtual sphere\n";
-    //        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-    //        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-    //          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-    //          << "\t" << RADIUS << "\t1 0 0\n";
-    //        *tempstream << "9\n  terminating special property\n";
-            tempstream->close();
-            tempstream->flush();
-            delete(tempstream);
-          }
-        }
-        if (DoTecplotOutput || DoRaster3DOutput)
-          TriangleFilesWritten++;
+      if (filename != NULL) {
+        if ((DoSingleStepOutput && ((mol->TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0) || (mol->TesselStruct->TrianglesOnBoundary.size() > 8520 && mol->TesselStruct->TrianglesOnBoundary.size() <= 8530)))) { // if we have a new triangle and want to output each new triangle configuration
+          mol->TesselStruct->Output(out, filename, mol);
+        }
       }
+      baseline = mol->TesselStruct->LinesOnBoundary.end();
     } else {
       //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
@@ -1220 +1042 @@
 
     N++;
-    baseline++;
     if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (flag)) {
       baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      baseline++;
       flag = false;
     }
   }
+  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
+  mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
+//  mol->GoToFirst();
+//  class TesselPoint *Runner = NULL;
+//  while (!mol->IsEnd()) {
+//    Runner = mol->GetPoint();
+//    *out << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
+//    if (!mol->TesselStruct->IsInnerPoint(out, Runner, LCList)) {
+//      *out << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
+//      mol->TesselStruct->AddBoundaryPointByDegeneratedTriangle(out, Runner, LCList);
+//    } else {
+//      *out << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
+//    }
+//    mol->GoToNext();
+//  }
 
   // Purges surplus triangles.
   mol->TesselStruct->RemoveDegeneratedTriangles();
-
-  // write final envelope
-  if (filename != 0) {
-    *out << Verbose(1) << "Writing final tecplot file\n";
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, -1);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append(Raster3DSuffix);
-      ofstream *raster = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, raster, mol->TesselStruct, mol);
-      raster->close();
-      delete(raster);
-    }
-  } else {
-    cerr << "ERROR: Could definitively not find all necessary triangles!" << endl;
-  }
 
   cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl;
@@ -1264 +1078 @@
     *out << Verbose(2) << "None." << endl;
 
-//  // Tests the IsInnerAtom() function.
-//  Vector x (0, 0, 0);
-//  *out << Verbose(0) << "Point to check: " << x << endl;
-//  *out << Verbose(0) << "Check: IsInnerPoint() returns " << mol->TesselStruct->IsInnerPoint(out, x, LCList)
-//    << "for vector " << x << "." << endl;
-//  TesselPoint* a = mol->TesselStruct->PointsOnBoundary.begin()->second->node;
-//  *out << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
-//  *out << Verbose(0) << "Check: IsInnerAtom() returns " << mol->TesselStruct->IsInnerPoint(out, a, LCList)
-//    << "for atom " << a << " (on boundary)." << endl;
-//  LinkedNodes *List = NULL;
-//  for (int i=0;i<NDIM;i++) { // each axis
-//    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
-//    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
-//      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
-//        List = LCList->GetCurrentCell();
-//        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-//        if (List != NULL) {
-//          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
-//            if (mol->TesselStruct->PointsOnBoundary.find((*Runner)->nr) == mol->TesselStruct->PointsOnBoundary.end()) {
-//              a = *Runner;
-//              i=3;
-//              for (int j=0;j<NDIM;j++)
-//                LCList->n[j] = LCList->N[j];
-//              break;
-//            }
-//          }
-//        }
-//      }
-//  }
-//  *out << Verbose(0) << "Check: IsInnerPoint() returns " << mol->TesselStruct->IsInnerPoint(out, a, LCList)
-//    << "for atom " << a << " (inside)." << endl;
+  // write final envelope
+  CalculateConcavityPerBoundaryPoint(out, mol->TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "");
 
   if (freeLC)
@@ -1300 +1086 @@
   *out << Verbose(0) << "End of FindNonConvexBorder\n";
 };
+
 
 /** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.

src/boundary.hpp

@@ -16 +16 @@
 
 #define DEBUG 1
-#define DoSingleStepOutput 0
+#define DoSingleStepOutput 1
 #define SingleStepWidth 1
-#define DoTecplotOutput 1
-#define DoRaster3DOutput 1
-#define DoVRMLOutput 1
-#define TecplotSuffix ".dat"
-#define Raster3DSuffix ".r3d"
-#define VRMLSUffix ".wrl"
 
 #define DistancePair pair < double, atom* >
@@ -38 +32 @@
 molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
 void FindConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename);
-void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
+void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LC, const double RADIUS, const char *tempbasename);
 double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename);
 void FindNextSuitablePoint(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);

src/tesselation.cpp

@@ -6 +6 @@
  */
 
+#include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "tesselation.hpp"
+#include "tesselationhelpers.hpp"
+#include "vector.hpp"
+
+class molecule;
 
 // ======================================== Points on Boundary =================================
@@ -36 +42 @@
 BoundaryPointSet::~BoundaryPointSet()
 {
-  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  //cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
@@ -66 +72 @@
 ostream & operator <<(ostream &ost, BoundaryPointSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
   return ost;
 }
@@ -124 +130 @@
         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
           if ((*Runner).second == this) {
-            cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+            //cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
             endpoints[i]->lines.erase(Runner);
             break;
           }
       } else { // there's just a single line left
-        if (endpoints[i]->lines.erase(Nr))
-          cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        if (endpoints[i]->lines.erase(Nr)) {
+          //cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        }
       }
       if (endpoints[i]->lines.empty()) {
-        cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+        //cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         if (endpoints[i] != NULL) {
           delete(endpoints[i]);
@@ -167 +174 @@
 
 /** Checks whether the adjacent triangles of a baseline are convex or not.
- * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
+ * We sum the two angles of each height vector with respect to the center of the baseline.
  * If greater/equal M_PI than we are convex.
  * \param *out output stream for debugging
@@ -177 +184 @@
   // get the two triangles
   if (triangles.size() != 2) {
-    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;
     return true;
   }
@@ -200 +207 @@
     NormalCheck.Scale(sign);
     sign = -sign;
-    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
+    if (runner->second->NormalVector.NormSquared() > MYEPSILON)
+      BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
+    else {
+      *out << Verbose(1) << "CRITICAL: Triangle " << *runner->second << " has zero normal vector!" << endl;
+      exit(255);
+    }
     node = runner->second->GetThirdEndpoint(this);
     if (node != NULL) {
-      //*out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
+      *out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
       helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
-      //*out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
+      *out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
     } else {
-      //*out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
+      //*out << Verbose(2) << "ERROR: I cannot find third node in triangle, something's wrong." << endl;
       return true;
     }
   }
-  //*out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
+  *out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
-    *out << Verbose(2) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
+    *out << Verbose(3) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
     return true;
   }
+  BaseLineNormal.Scale(-1.);
   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
   if ((angle - M_PI) > -MYEPSILON) {
-    *out << Verbose(2) << "ACCEPT: Angle is greater than pi: convex." << endl;
+    *out << Verbose(3) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
   } else {
-    *out << Verbose(2) << "REJECT: Angle is less than pi: concave." << endl;
+    *out << Verbose(3) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
   }
@@ -261 +274 @@
 ostream & operator <<(ostream &ost, BoundaryLineSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
   return ost;
 };
@@ -331 +344 @@
   for (int i = 0; i < 3; i++) {
     if (lines[i] != NULL) {
-      if (lines[i]->triangles.erase(Nr))
-        cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      if (lines[i]->triangles.erase(Nr)) {
+        //cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      }
       if (lines[i]->triangles.empty()) {
-          cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+          //cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
           delete (lines[i]);
           lines[i] = NULL;
@@ -340 +354 @@
     }
   }
-  cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  //cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 
@@ -450 +464 @@
 };
 
+/** Checks whether three given \a *Points coincide with triangle's endpoints.
+ * \param *Points[3] pointer to BoundaryPointSet
+ * \return true - is the very triangle, false - is not
+ */
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+{
+  return (((endpoints[0] == T->endpoints[0])
+            || (endpoints[0] == T->endpoints[1])
+            || (endpoints[0] == T->endpoints[2])
+          ) && (
+            (endpoints[1] == T->endpoints[0])
+            || (endpoints[1] == T->endpoints[1])
+            || (endpoints[1] == T->endpoints[2])
+          ) && (
+            (endpoints[2] == T->endpoints[0])
+            || (endpoints[2] == T->endpoints[1])
+            || (endpoints[2] == T->endpoints[2])
+
+          ));
+};
+
 /** Returns the endpoint which is not contained in the given \a *line.
  * \param *line baseline defining two endpoints
@@ -484 +519 @@
 ostream &operator <<(ostream &ost, BoundaryTriangleSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
-      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
+      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
   return ost;
 };
@@ -504 +539 @@
 TesselPoint::~TesselPoint()
 {
-  Free((void **)&Name, "TesselPoint::~TesselPoint: *Name");
 };
 
@@ -511 +545 @@
 ostream & operator << (ostream &ost, const TesselPoint &a)
 {
-  ost << "[" << (a.Name) << "|" << &a << "]";
+  ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
 };
@@ -568 +602 @@
   TrianglesOnBoundaryCount = 0;
   InternalPointer = PointsOnBoundary.begin();
+  LastTriangle = NULL;
+  TriangleFilesWritten = 0;
 }
 ;
@@ -584 +620 @@
       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
   }
+  cout << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -1232 +1269 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
-  cout << Verbose(4) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
+  cout << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
@@ -1252 +1289 @@
   TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
   TrianglesOnBoundaryCount++;
+
+  // set as last new triangle
+  LastTriangle = BTS;
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
@@ -1524 +1564 @@
   CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
+    Oben.Zero();
     Oben.x[k] = 1.;
     FirstPoint = MaxPoint[k];
@@ -1532 +1573 @@
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(FirstPoint, NULL, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(FirstPoint, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     SecondPoint = OptCandidate;
     if (SecondPoint == NULL)  // have we found a second point?
       continue;
-    else
-      cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
 
     helper.CopyVector(FirstPoint->node);
@@ -1559 +1598 @@
 
     // adding point 1 and point 2 and add the line between them
+    cout << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
     AddTesselationPoint(FirstPoint, 0);
+    cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
     AddTesselationPoint(SecondPoint, 1);
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -1632 +1673 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC)
 {
   cout << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
@@ -1667 +1708 @@
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
-    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
     // construct old center
@@ -1754 +1795 @@
         result = false;
       }
-    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
+    } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
         AddTesselationPoint((*it)->point, 0);
         AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
@@ -1795 +1836 @@
     // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
     BaseRay = BLS[0];
+    if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) {
+      *out << Verbose(1) << "CRITICAL: Triangle " << *BTS << " has zero normal vector!" << endl;
+      exit(255);
+    }
+
   }
 
@@ -1813 +1859 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return NULL - concave, otherwise endpoint that makes it concave
+ * \return NULL - convex, otherwise endpoint that makes it concave
  */
 class BoundaryPointSet *Tesselation::IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base)
@@ -1890 +1936 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return true - line was changed, false - same line as before
- */
-bool Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
+ * \return volume change due to flipping (0 - then no flipped occured)
+ */
+double Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
 {
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
+  double volume;
 
   int m=0;
@@ -1915 +1962 @@
   Distance.CopyVector(ClosestPoint[1]);
   Distance.SubtractVector(ClosestPoint[0]);
+
+  // calculate volume
+  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->node, OtherBase->endpoints[0]->node->node, OtherBase->endpoints[1]->node->node, Base->endpoints[0]->node->node);
 
   // delete the temporary other base line and the closest points
@@ -1929 +1979 @@
     if (Base->triangles.size() < 2) {
       *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
-      return false;
+      return 0.;
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
@@ -1939 +1989 @@
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
       *out << Verbose(2) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;
-      FlipBaseline(out, Base);
-      return true;
+      // calculate volume summand as a general tetraeder
+      return volume;
     } else {  // Base higher than OtherBase -> do nothing
       *out << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl;
-      return false;
-    }
-  }
-};
-
-/** Returns the closest point on \a *Base with respect to \a *OtherBase.
- * \param *out output stream for debugging
- * \param *Base reference line
- * \param *OtherBase other base line
- * \return Vector on reference line that has closest distance
- */
-Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
-{
-  // construct the plane of the two baselines (i.e. take both their directional vectors)
-  Vector Normal;
-  Vector Baseline, OtherBaseline;
-  Baseline.CopyVector(Base->endpoints[1]->node->node);
-  Baseline.SubtractVector(Base->endpoints[0]->node->node);
-  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
-  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
-  Normal.CopyVector(&Baseline);
-  Normal.VectorProduct(&OtherBaseline);
-  Normal.Normalize();
-  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
-
-  // project one offset point of OtherBase onto this plane (and add plane offset vector)
-  Vector NewOffset;
-  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
-  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
-  NewOffset.ProjectOntoPlane(&Normal);
-  NewOffset.AddVector(Base->endpoints[0]->node->node);
-  Vector NewDirection;
-  NewDirection.CopyVector(&NewOffset);
-  NewDirection.AddVector(&OtherBaseline);
-
-  // calculate the intersection between this projected baseline and Base
-  Vector *Intersection = new Vector;
-  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
-  Normal.CopyVector(Intersection);
-  Normal.SubtractVector(Base->endpoints[0]->node->node);
-  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
-
-  return Intersection;
+      return 0.;
+    }
+  }
 };
 
@@ -1993 +2003 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return true - flipping successful, false - something went awry
- */
-bool Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
+ * \return pointer to allocated new baseline - flipping successful, NULL - something went awry
+ */
+class BoundaryLineSet * Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
 {
   class BoundaryLineSet *OldLines[4], *NewLine;
@@ -2009 +2019 @@
   if (Base->triangles.size() < 2) {
     *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
-    return false;
+    return NULL;
   }
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
@@ -2045 +2055 @@
   if (i<4) {
     *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-    return false;
+    return NULL;
   }
   for (int j=0;j<4;j++)
     if (OldLines[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-      return false;
+      return NULL;
     }
   for (int j=0;j<2;j++)
     if (OldPoints[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-      return false;
+      return NULL;
     }
 
@@ -2099 +2109 @@
   } else {
     *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
-    return false;
+    return NULL;
   }
 
   *out << Verbose(1) << "End of FlipBaseline" << endl;
-  return true;
+  return NewLine;
 };
 
@@ -2109 +2119 @@
 /** Finds the second point of starting triangle.
  * \param *a first node
- * \param *Candidate pointer to candidate node on return
  * \param Oben vector indicating the outside
  * \param OptCandidate reference to recommended candidate on return
@@ -2116 +2125 @@
  * \param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindSecondPointForTesselation(TesselPoint* a, TesselPoint* Candidate, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
+void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
 {
   cout << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;
   Vector AngleCheck;
+  class TesselPoint* Candidate = NULL;
   double norm = -1., angle;
   LinkedNodes *List = NULL;
@@ -2254 +2264 @@
   cout << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl;
 
-  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+  cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
 
   // construct center of circle
@@ -2279 +2289 @@
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
+      //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
@@ -2455 +2466 @@
   TesselPoint *trianglePoints[3];
   TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
 
   if (LinesOnBoundary.empty()) {
@@ -2465 +2477 @@
   // check whether closest point is "too close" :), then it's inside
   if (trianglePoints[0] == NULL) {
-    *out << Verbose(1) << "Is the only point, no one else is closeby." << endl;
+    *out << Verbose(2) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
   if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    *out << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
-    return NULL;
-  }
-  list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(out, trianglePoints[0], x);
-  trianglePoints[1] = connectedClosestPoints->front();
-  trianglePoints[2] = connectedClosestPoints->back();
-  for (int i=0;i<3;i++) {
-    if (trianglePoints[i] == NULL) {
-      *out << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-    }
-    //*out << Verbose(1) << "List of possible points:" << endl;
-    //*out << Verbose(2) << *trianglePoints[i] << endl;
-  }
-
-  list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);
-
-  delete(connectedClosestPoints);
+    *out << Verbose(2) << "Point is right on a tesselation point, no nearest triangle." << endl;
+    PointMap::iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
+    triangles = new list<BoundaryTriangleSet*>;
+    if (PointRunner != PointsOnBoundary.end()) {
+      for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
+        for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
+          triangles->push_back(TriangleRunner->second);
+      triangles->sort();
+      triangles->unique();
+    } else {
+      PointRunner = PointsOnBoundary.find(SecondPoint->nr);
+      trianglePoints[0] = SecondPoint;
+      if (PointRunner != PointsOnBoundary.end()) {
+        for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
+          for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
+            triangles->push_back(TriangleRunner->second);
+        triangles->sort();
+        triangles->unique();
+      } else {
+        *out << Verbose(1) << "ERROR: I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
+        return NULL;
+      }
+    }
+  } else {
+    list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(out, trianglePoints[0], x);
+    trianglePoints[1] = connectedClosestPoints->front();
+    trianglePoints[2] = connectedClosestPoints->back();
+    for (int i=0;i<3;i++) {
+      if (trianglePoints[i] == NULL) {
+        *out << Verbose(1) << "ERROR: IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+      }
+      *out << Verbose(1) << "List of triangle points:" << endl;
+      *out << Verbose(2) << *trianglePoints[i] << endl;
+    }
+
+    triangles = FindTriangles(trianglePoints);
+    *out << Verbose(1) << "List of possible triangles:" << endl;
+    for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      *out << Verbose(2) << **Runner << endl;
+
+    delete(connectedClosestPoints);
+  }
   
   if (triangles->empty()) {
-    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
+    *out << Verbose(0) << "ERROR: There is no nearest triangle. Please check the tesselation structure.";
+    delete(triangles);
     return NULL;
   } else
@@ -2504 +2542 @@
   class BoundaryTriangleSet *result = NULL;
   list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(out, x, LC);
+  Vector Center;
 
   if (triangles == NULL)
     return NULL;
 
-  if (x->ScalarProduct(&triangles->front()->NormalVector) < 0)
-    result = triangles->back();
-  else
+  if (triangles->size() == 1) { // there is no degenerate case
     result = triangles->front();
-
+    *out << Verbose(2) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
+  } else {
+    result = triangles->front();
+    result->GetCenter(&Center);
+    Center.SubtractVector(x);
+    *out << Verbose(2) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
+    if (Center.ScalarProduct(&result->NormalVector) < 0) {
+      result = triangles->back();
+      *out << Verbose(2) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
+      if (Center.ScalarProduct(&result->NormalVector) < 0) {
+        *out << Verbose(1) << "ERROR: Front and back side yield NormalVector in wrong direction!" << endl;
+      }
+    }
+  }
   delete(triangles);
   return result;
@@ -2527 +2577 @@
 {
   class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
-  if (result == NULL)
+  Vector Center;
+
+  if (result == NULL) {// is boundary point or only point in point cloud?
+    *out << Verbose(1) << Point << " is the only point in vicinity." << endl;
+    return false;
+  }
+
+  result->GetCenter(&Center);
+  *out << Verbose(3) << "INFO: Central point of the triangle is " << Center << "." << endl;
+  Center.SubtractVector(&Point);
+  *out << Verbose(3) << "INFO: Vector from center to point to test is " << Center << "." << endl;
+  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
+    *out << Verbose(1) << Point << " is an inner point." << endl;
     return true;
-  if (Point.ScalarProduct(&result->NormalVector) < 0)
-    return true;
-  else
+  } else {
+    *out << Verbose(1) << Point << " is NOT an inner point." << endl;
     return false;
+  }
 }
 
@@ -2544 +2606 @@
 bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
 {
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, Point->node, LC);
-  if (result == NULL)
-    return true;
-  if (Point->node->ScalarProduct(&result->NormalVector) < 0)
-    return true;
-  else
-    return false;
+  return IsInnerPoint(out, *(Point->node), LC);
 }
 
@@ -2705 +2761 @@
   }
 
-  map <class BoundaryLineSet *, bool> Touched;
-  map <class BoundaryLineSet *, bool>::iterator line;
-  for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++)
-    Touched.insert( pair <class BoundaryLineSet *, bool>(runner->second, false) );
+  map <class BoundaryLineSet *, bool> TouchedLine;
+  map <class BoundaryTriangleSet *, bool> TouchedTriangle;
+  map <class BoundaryLineSet *, bool>::iterator LineRunner;
+  map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
+  for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
+    TouchedLine.insert( pair <class BoundaryLineSet *, bool>(Runner->second, false) );
+    for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
+      TouchedTriangle.insert( pair <class BoundaryTriangleSet *, bool>(Sprinter->second, false) );
+  }
   if (!ReferencePoint->lines.empty()) {
     for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++) {
-      line = Touched.find(runner->second);
-      if (line == Touched.end()) {
+      LineRunner = TouchedLine.find(runner->second);
+      if (LineRunner == TouchedLine.end()) {
         *out << Verbose(2) << "ERROR: I could not find " << *runner->second << " in the touched list." << endl;
-      } else if (!line->second) {
-        line->second = true;
+      } else if (!LineRunner->second) {
+        LineRunner->second = true;
         connectedPath = new list<TesselPoint*>;
         triangle = NULL;
@@ -2728 +2789 @@
 
           // find next triangle
-          for (TriangleMap::iterator TriangleRunner = CurrentLine->triangles.begin(); TriangleRunner != CurrentLine->triangles.end(); TriangleRunner++) {
-            if (TriangleRunner->second != triangle) { // look for first triangle not equal to old one
-              triangle = TriangleRunner->second;
-              break;
+          for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
+            *out << Verbose(3) << "INFO: Inspecting triangle " << *Runner->second << "." << endl;
+            if ((Runner->second != triangle)) { // look for first triangle not equal to old one
+              triangle = Runner->second;
+              TriangleRunner = TouchedTriangle.find(triangle);
+              if (TriangleRunner != TouchedTriangle.end()) {
+                if (!TriangleRunner->second) {
+                  TriangleRunner->second = true;
+                  *out << Verbose(3) << "INFO: Connecting triangle is " << *triangle << "." << endl;
+                  break;
+                } else {
+                  *out << Verbose(3) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl;
+                  triangle = NULL;
+                }
+              } else {
+                *out << Verbose(2) << "ERROR: I could not find " << *triangle << " in the touched list." << endl;
+                triangle = NULL;
+              }
             }
           }
+          if (triangle == NULL)
+            break;
           // find next line
           for (int i=0;i<3;i++) {
             if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
               CurrentLine = triangle->lines[i];
-
+              *out << Verbose(3) << "INFO: Connecting line is " << *CurrentLine << "." << endl;
               break;
             }
           }
-          line = Touched.find(CurrentLine);
-          if (line == Touched.end())
+          LineRunner = TouchedLine.find(CurrentLine);
+          if (LineRunner == TouchedLine.end())
             *out << Verbose(2) << "ERROR: I could not find " << *CurrentLine << " in the touched list." << endl;
           else
-            line->second = true;
+            LineRunner->second = true;
           // find next point
           CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
@@ -2753 +2830 @@
         // last point is missing, as it's on start line
         *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
-        connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+        if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
+          connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
 
         ListOfPaths->push_back(connectedPath);
@@ -2853 +2931 @@
 
 /** Removes a boundary point from the envelope while keeping it closed.
- * We create new triangles and remove the old ones connected to the point.
+ * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
+ *  -# a closed path(s) of boundary points surrounding the point to be removed is constructed
+ *  -# on each closed path, we pick three adjacent points, create a triangle with them and subtract the middle point from the path
+ *  -# we advance two points (i.e. the next triangle will start at the ending point of the last triangle) and continue as before
+ *  -# the surface is closed, when the path is empty
+ * Thereby, we (hopefully) make sure that the removed points remains beneath the surface (this is checked via IsInnerPoint eventually).
  * \param *out output stream for debugging
  * \param *point point to be removed
@@ -2861 +2944 @@
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint, TetraederVector[3];
+  Vector OldPoint, NormalVector;
   double volume = 0;
   int count = 0;
@@ -2887 +2970 @@
     count+=LineRunner->second->triangles.size();
   map<class BoundaryTriangleSet *, int> Candidates;
-  for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
@@ -2897 +2980 @@
   // remove all triangles
   count=0;
+  NormalVector.Zero();
   for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
     *out << Verbose(3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
+    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
     RemoveTesselationTriangle(Runner->first);
     count++;
@@ -2907 +2992 @@
   list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
   map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
+  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  double angle;
+  double smallestangle;
+  Vector Point, Reference, OrthogonalVector;
   if (count > 2) {  // less than three triangles, then nothing will be created
     class TesselPoint *TriangleCandidates[3];
@@ -2915 +3004 @@
 
       connectedPath = *ListRunner;
-      // initialize the path to start
-      list<TesselPoint*>::iterator CircleRunner = connectedPath->begin();
-      list<TesselPoint*>::iterator OtherCircleRunner = connectedPath->begin();
-      class TesselPoint *CentralNode = *CircleRunner;
 
       // re-create all triangles by going through connected points list
-      // advance two with CircleRunner and one with OtherCircleRunner
-      CircleRunner++;
-      CircleRunner++;
-      OtherCircleRunner++;
-      cout << Verbose(3) << "INFO: CentralNode is " << *CentralNode << "." << endl;
-      for (; (OtherCircleRunner != connectedPath->end()) && (CircleRunner != connectedPath->end()); (CircleRunner++), (OtherCircleRunner++)) {
-        cout << Verbose(4) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
-        cout << Verbose(4) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
-        *out << Verbose(3) << "INFO: Creating triangle " << CentralNode->Name << ", " << (*OtherCircleRunner)->Name << " and " << (*CircleRunner)->Name << "." << endl;
+      list<class BoundaryLineSet *> NewLines;
+      for (;!connectedPath->empty();) {
+        // search middle node with widest angle to next neighbours
+        EndNode = connectedPath->end();
+        smallestangle = 0.;
+        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+          cout << Verbose(3) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+          // construct vectors to next and previous neighbour
+          StartNode = MiddleNode;
+          if (StartNode == connectedPath->begin())
+            StartNode = connectedPath->end();
+          StartNode--;
+          //cout << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
+          Point.CopyVector((*StartNode)->node);
+          Point.SubtractVector((*MiddleNode)->node);
+          StartNode = MiddleNode;
+          StartNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          //cout << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
+          Reference.CopyVector((*StartNode)->node);
+          Reference.SubtractVector((*MiddleNode)->node);
+          OrthogonalVector.CopyVector((*MiddleNode)->node);
+          OrthogonalVector.SubtractVector(&OldPoint);
+          OrthogonalVector.MakeNormalVector(&Reference);
+          angle = GetAngle(Point, Reference, OrthogonalVector);
+          //if (angle < M_PI)  // no wrong-sided triangles, please?
+            if(fabs(angle - M_PI) < fabs(smallestangle - M_PI)) {  // get straightest angle (i.e. construct those triangles with smallest area first)
+              smallestangle = angle;
+              EndNode = MiddleNode;
+            }
+        }
+        MiddleNode = EndNode;
+        if (MiddleNode == connectedPath->end()) {
+          cout << Verbose(1) << "CRITICAL: Could not find a smallest angle!" << endl;
+          exit(255);
+        }
+        StartNode = MiddleNode;
+        if (StartNode == connectedPath->begin())
+          StartNode = connectedPath->end();
+        StartNode--;
+        EndNode++;
+        if (EndNode == connectedPath->end())
+          EndNode = connectedPath->begin();
+        cout << Verbose(4) << "INFO: StartNode is " << **StartNode << "." << endl;
+        cout << Verbose(4) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+        cout << Verbose(4) << "INFO: EndNode is " << **EndNode << "." << endl;
+        *out << Verbose(3) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl;
+        TriangleCandidates[0] = *StartNode;
+        TriangleCandidates[1] = *MiddleNode;
+        TriangleCandidates[2] = *EndNode;
+        triangle = GetPresentTriangle(out, TriangleCandidates);
+        if (triangle != NULL) {
+          cout << Verbose(1) << "WARNING: New triangle already present, skipping!" << endl;
+          StartNode++;
+          MiddleNode++;
+          EndNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+          continue;
+        }
         *out << Verbose(5) << "Adding new triangle points."<< endl;
-        TriangleCandidates[0] = CentralNode;
-        TriangleCandidates[1] = *OtherCircleRunner;
-        TriangleCandidates[2] = *CircleRunner;
-        triangle = GetPresentTriangle(out, TriangleCandidates);
-        AddTesselationPoint(CentralNode, 0);
-        AddTesselationPoint(*OtherCircleRunner, 1);
-        AddTesselationPoint(*CircleRunner, 2);
+        AddTesselationPoint(*StartNode, 0);
+        AddTesselationPoint(*MiddleNode, 1);
+        AddTesselationPoint(*EndNode, 2);
         *out << Verbose(5) << "Adding new triangle lines."<< endl;
         AddTesselationLine(TPS[0], TPS[1], 0);
         AddTesselationLine(TPS[0], TPS[2], 1);
+        NewLines.push_back(BLS[1]);
         AddTesselationLine(TPS[1], TPS[2], 2);
         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+        BTS->GetNormalVector(NormalVector);
         AddTesselationTriangle();
         // calculate volume summand as a general tetraeder
-        for (int j=0;j<3;j++) {
-          TetraederVector[j].CopyVector(TPS[j]->node->node);
-          TetraederVector[j].SubtractVector(&OldPoint);
+        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->node, TPS[1]->node->node, TPS[2]->node->node, &OldPoint);
+        // advance number
+        count++;
+
+        // prepare nodes for next triangle
+        StartNode = EndNode;
+        cout << Verbose(4) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl;
+        connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+        if (connectedPath->size() == 2) { // we are done
+          connectedPath->remove(*StartNode); // remove the start node
+          connectedPath->remove(*EndNode); // remove the end node
+          break;
+        } else if (connectedPath->size() < 2) { // something's gone wrong!
+          cout << Verbose(1) << "CRITICAL: There are only two endpoints left!" << endl;
+          exit(255);
+        } else {
+          MiddleNode = StartNode;
+          MiddleNode++;
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          EndNode = MiddleNode;
+          EndNode++;
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
         }
-        OldPoint.CopyVector(&TetraederVector[0]);
-        OldPoint.VectorProduct(&TetraederVector[1]);
-        volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
-        // advance number
-        NumberRunner++;
-        count++;
-        if (NumberRunner == Candidates.end())
-          *out << Verbose(2) << "WARN: Maximum of numbers reached!" << endl;
       }
+      // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
+      if (NewLines.size() > 1) {
+        list<class BoundaryLineSet *>::iterator Candidate;
+        class BoundaryLineSet *OtherBase = NULL;
+        double tmp, maxgain;
+        do {
+          maxgain = 0;
+          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+            tmp = PickFarthestofTwoBaselines(out, *Runner);
+            if (maxgain < tmp) {
+              maxgain = tmp;
+              Candidate = Runner;
+            }
+          }
+          if (maxgain != 0) {
+            volume += maxgain;
+            cout << Verbose(3) << "Flipping baseline with highest volume" << **Candidate << "." << endl;
+            OtherBase = FlipBaseline(out, *Candidate);
+            NewLines.erase(Candidate);
+            NewLines.push_back(OtherBase);
+          }
+        } while (maxgain != 0.);
+      }
+
       ListOfClosedPaths->remove(connectedPath);
       delete(connectedPath);
@@ -2971 +3147 @@
     *out << Verbose(1) << "No need to create any triangles." << endl;
   }
-
   delete(ListOfClosedPaths);
 
+  *out << Verbose(1) << "Removed volume is " << volume << "." << endl;
+
   return volume;
 };
 
 
-
-/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
- * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
- * make it bigger (i.e. closing one (the baseline) and opening two new ones).
- * \param TPS[3] nodes of the triangle
- * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
- */
-bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3])
-{
-  bool result = false;
-  int counter = 0;
-
-  // check all three points
-  for (int i=0;i<3;i++)
-    for (int j=i+1; j<3; j++) {
-      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
-        LineMap::iterator FindLine;
-        pair<LineMap::iterator,LineMap::iterator> FindPair;
-        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
-        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-          // If there is a line with less than two attached triangles, we don't need a new line.
-          if (FindLine->second->triangles.size() < 2) {
-            counter++;
-            break;  // increase counter only once per edge
-          }
-        }
-      } else { // no line
-        cout << Verbose(1) << "The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
-        result = true;
-      }
-    }
-  if (counter > 1) {
-    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
-    result = true;
-  }
-  return result;
-};
-
-
-/** Sort function for the candidate list.
- */
-bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
-{
-  Vector BaseLineVector, OrthogonalVector, helper;
-  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
-    //return false;
-    exit(1);
-  }
-  // create baseline vector
-  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
-  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
-  BaseLineVector.Normalize();
-
-  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
-  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
-  helper.SubtractVector(candidate1->point->node);
-  OrthogonalVector.CopyVector(&helper);
-  helper.VectorProduct(&BaseLineVector);
-  OrthogonalVector.SubtractVector(&helper);
-  OrthogonalVector.Normalize();
-
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(candidate1->point->node);
-  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
-  double phi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-  helper.CopyVector(candidate2->point->node);
-  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
-  double psi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    psi = 2.*M_PI - psi;
-  }
-
-  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
-  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
-
-  // return comparison
-  return phi < psi;
-};
-
-/**
- * Finds the point which is second closest to the provided one.
- *
- * @param Point to which to find the second closest other point
- * @param linked cell structure
- *
- * @return point which is second closest to the provided one
- */
-TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
-{
-  LinkedNodes *List = NULL;
-  TesselPoint* closestPoint = NULL;
-  TesselPoint* secondClosestPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
-            if (currentNorm < distance) {
-              // remember second point
-              secondDistance = distance;
-              secondClosestPoint = closestPoint;
-              // mark down new closest point
-              distance = currentNorm;
-              closestPoint = (*Runner);
-              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
-            }
-          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  return secondClosestPoint;
-};
-
-/**
- * Finds the point which is closest to the provided one.
- *
- * @param Point to which to find the closest other point
- * @param SecondPoint the second closest other point on return, NULL if none found
- * @param linked cell structure
- *
- * @return point which is closest to the provided one, NULL if none found
- */
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
-{
-  LinkedNodes *List = NULL;
-  TesselPoint* closestPoint = NULL;
-  SecondPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
-            if (currentNorm < distance) {
-              secondDistance = distance;
-              SecondPoint = closestPoint;
-              distance = currentNorm;
-              closestPoint = (*Runner);
-              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
-            } else if (currentNorm < secondDistance) {
-              secondDistance = currentNorm;
-              SecondPoint = (*Runner);
-              cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
-            }
-          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  return closestPoint;
-};
 
 /**
@@ -3207 +3191 @@
               FindTriangle = FindLine->second->triangles.begin();
               for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
-                  )
-                ) {
+                if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
                   result->push_back(FindTriangle->second);
                 }
@@ -3238 +3209 @@
 
 /**
+ * Finds all degenerated lines within the tesselation structure.
+ *
+ * @return map of keys of degenerated line pairs, each line occurs twice
+ *         in the list, once as key and once as value
+ */
+map<int, int> * Tesselation::FindAllDegeneratedLines()
+{
+  map<int, class BoundaryLineSet *> AllLines;
+  map<int, int> * DegeneratedLines = new map<int, int>;
+
+  // sanity check
+  if (LinesOnBoundary.empty()) {
+    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
+    return DegeneratedLines;
+  }
+
+  LineMap::iterator LineRunner1;
+  pair<LineMap::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) );
+    }
+  }
+
+  AllLines.clear();
+
+  cout << Verbose(1) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
+  map<int,int>::iterator it;
+  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++)
+      cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;
+
+  return DegeneratedLines;
+}
+
+/**
  * Finds all degenerated triangles within the tesselation structure.
  *
@@ -3243 +3251 @@
  *         in the list, once as key and once as value
  */
-map<int, int> Tesselation::FindAllDegeneratedTriangles()
-{
-  map<int, int> DegeneratedTriangles;
-
-  // sanity check
-  if (LinesOnBoundary.empty()) {
-    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
-    return DegeneratedTriangles;
-  }
-
-  LineMap::iterator LineRunner1, LineRunner2;
-
-  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
-    for (LineRunner2 = LinesOnBoundary.begin(); LineRunner2 != LinesOnBoundary.end(); ++LineRunner2) {
-      if ((LineRunner1->second != LineRunner2->second)
-        && (LineRunner1->second->endpoints[0] == LineRunner2->second->endpoints[0])
-        && (LineRunner1->second->endpoints[1] == LineRunner2->second->endpoints[1])
-      ) {
-        TriangleMap::iterator TriangleRunner1 = LineRunner1->second->triangles.begin(),
-          TriangleRunner2 = LineRunner2->second->triangles.begin();
-
-        for (; TriangleRunner1 != LineRunner1->second->triangles.end(); ++TriangleRunner1) {
-          for (; TriangleRunner2 != LineRunner2->second->triangles.end(); ++TriangleRunner2) {
-            if ((TriangleRunner1->second != TriangleRunner2->second)
-              && (TriangleRunner1->second->endpoints[0] == TriangleRunner2->second->endpoints[0])
-              && (TriangleRunner1->second->endpoints[1] == TriangleRunner2->second->endpoints[1])
-              && (TriangleRunner1->second->endpoints[2] == TriangleRunner2->second->endpoints[2])
-            ) {
-              DegeneratedTriangles[TriangleRunner1->second->Nr] = TriangleRunner2->second->Nr;
-              DegeneratedTriangles[TriangleRunner2->second->Nr] = TriangleRunner1->second->Nr;
-            }
-          }
+map<int, int> * Tesselation::FindAllDegeneratedTriangles()
+{
+  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
+  map<int, int> * DegeneratedTriangles = new map<int, int>;
+
+  TriangleMap::iterator TriangleRunner1, TriangleRunner2;
+  LineMap::iterator Liner;
+  class BoundaryLineSet *line1 = NULL, *line2 = NULL;
+
+  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+    // run over both lines' triangles
+    Liner = LinesOnBoundary.find(LineRunner->first);
+    if (Liner != LinesOnBoundary.end())
+      line1 = Liner->second;
+    Liner = LinesOnBoundary.find(LineRunner->second);
+    if (Liner != LinesOnBoundary.end())
+      line2 = Liner->second;
+    for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) {
+      for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
+        if ((TriangleRunner1->second != TriangleRunner2->second)
+          && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr) );
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr) );
         }
       }
     }
   }
-
-  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles.size() << " triangles." << endl;
+  delete(DegeneratedLines);
+
+  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
   map<int,int>::iterator it;
-  for (it = DegeneratedTriangles.begin(); it != DegeneratedTriangles.end(); it++)
+  for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;
 
@@ -3294 +3294 @@
 void Tesselation::RemoveDegeneratedTriangles()
 {
-  map<int, int> DegeneratedTriangles = FindAllDegeneratedTriangles();
-
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles.begin();
-    TriangleKeyRunner != DegeneratedTriangles.end(); ++TriangleKeyRunner
+  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  TriangleMap::iterator finder;
+  BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
+  int count  = 0;
+
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+    TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
-    BoundaryTriangleSet *triangle = TrianglesOnBoundary.find(TriangleKeyRunner->first)->second,
-      *partnerTriangle = TrianglesOnBoundary.find(TriangleKeyRunner->second)->second;
+    finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
+    if (finder != TrianglesOnBoundary.end())
+      triangle = finder->second;
+    else
+      break;
+    finder = TrianglesOnBoundary.find(TriangleKeyRunner->second);
+    if (finder != TrianglesOnBoundary.end())
+      partnerTriangle = finder->second;
+    else
+      break;
 
     bool trianglesShareLine = false;
@@ -3312 +3323 @@
       && (triangle->endpoints[0]->LinesCount > 2)
     ) {
+      // check whether we have to fix lines
+      BoundaryTriangleSet *Othertriangle = NULL;
+      BoundaryTriangleSet *OtherpartnerTriangle = NULL;
+      TriangleMap::iterator TriangleRunner;
+      for (int i = 0; i < 3; ++i)
+        for (int j = 0; j < 3; ++j)
+          if (triangle->lines[i] != partnerTriangle->lines[j]) {
+            // get the other two triangles
+            for (TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); ++TriangleRunner)
+              if (TriangleRunner->second != triangle) {
+                Othertriangle = TriangleRunner->second;
+              }
+            for (TriangleRunner = partnerTriangle->lines[i]->triangles.begin(); TriangleRunner != partnerTriangle->lines[i]->triangles.end(); ++TriangleRunner)
+              if (TriangleRunner->second != partnerTriangle) {
+                OtherpartnerTriangle = TriangleRunner->second;
+              }
+            /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
+            // the line of triangle receives the degenerated ones
+            triangle->lines[i]->triangles.erase(Othertriangle->Nr);
+            triangle->lines[i]->triangles.insert( TrianglePair( partnerTriangle->Nr, partnerTriangle) );
+            for (int k=0;k<3;k++)
+              if (triangle->lines[i] == Othertriangle->lines[k]) {
+                Othertriangle->lines[k] = partnerTriangle->lines[j];
+                break;
+              }
+            // the line of partnerTriangle receives the non-degenerated ones
+            partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr);
+            partnerTriangle->lines[j]->triangles.insert( TrianglePair( Othertriangle->Nr, Othertriangle) );
+            partnerTriangle->lines[j] = triangle->lines[i];
+          }
+
+      // erase the pair
+      count += (int) DegeneratedTriangles->erase(triangle->Nr);
       cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
       RemoveTesselationTriangle(triangle);
+      count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
       cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
       RemoveTesselationTriangle(partnerTriangle);
-      DegeneratedTriangles.erase(DegeneratedTriangles.find(partnerTriangle->Nr));
     } else {
       cout << Verbose(1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
@@ -3323 +3367 @@
     }
   }
+  delete(DegeneratedTriangles);
+
+  cout << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
 }
 
-/** Gets the angle between a point and a reference relative to the provided center.
- * We have two shanks point and reference between which the angle is calculated
- * and by scalar product with OrthogonalVector we decide the interval.
- * @param point to calculate the angle for
- * @param reference to which to calculate the angle
- * @param OrthogonalVector points in direction of [pi,2pi] interval
- *
- * @return angle between point and reference
- */
-double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
-{
-  if (reference.IsZero())
-    return M_PI;
-
-  // calculate both angles and correct with in-plane vector
-  if (point.IsZero())
-    return M_PI;
-  double phi = point.Angle(&reference);
-  if (OrthogonalVector.ScalarProduct(&point) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-
-  cout << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
-
-  return phi;
-}
-
+/** Adds an outside Tesselpoint to the envelope via (two) degenerated triangles.
+ * We look for the closest point on the boundary, we look through its connected boundary lines and
+ * seek the one with the minimum angle between its center point and the new point and this base line.
+ * We open up the line by adding a degenerated triangle, whose other side closes the base line again.
+ * \param *out output stream for debugging
+ * \param *point point to add
+ * \param *LC Linked Cell structure to find nearest point
+ */
+void Tesselation::AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC)
+{
+  *out << Verbose(2) << "Begin of AddBoundaryPointByDegeneratedTriangle" << endl;
+
+  // find nearest boundary point
+  class TesselPoint *BackupPoint = NULL;
+  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class BoundaryPointSet *NearestBoundaryPoint = NULL;
+  PointMap::iterator PointRunner;
+
+  if (NearestPoint == point)
+    NearestPoint = BackupPoint;
+  PointRunner = PointsOnBoundary.find(NearestPoint->nr);
+  if (PointRunner != PointsOnBoundary.end()) {
+    NearestBoundaryPoint = PointRunner->second;
+  } else {
+    *out << Verbose(1) << "ERROR: I cannot find the boundary point." << endl;
+    return;
+  }
+  *out << Verbose(2) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl;
+
+  // go through its lines and find the best one to split
+  Vector CenterToPoint;
+  Vector BaseLine;
+  double angle, BestAngle = 0.;
+  class BoundaryLineSet *BestLine = NULL;
+  for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
+    BaseLine.CopyVector(Runner->second->endpoints[0]->node->node);
+    BaseLine.SubtractVector(Runner->second->endpoints[1]->node->node);
+    CenterToPoint.CopyVector(Runner->second->endpoints[0]->node->node);
+    CenterToPoint.AddVector(Runner->second->endpoints[1]->node->node);
+    CenterToPoint.Scale(0.5);
+    CenterToPoint.SubtractVector(point->node);
+    angle = CenterToPoint.Angle(&BaseLine);
+    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
+      BestAngle = angle;
+      BestLine = Runner->second;
+    }
+  }
+
+  // remove one triangle from the chosen line
+  class BoundaryTriangleSet *TempTriangle = (BestLine->triangles.begin())->second;
+  BestLine->triangles.erase(TempTriangle->Nr);
+  int nr = -1;
+  for (int i=0;i<3; i++) {
+    if (TempTriangle->lines[i] == BestLine) {
+      nr = i;
+      break;
+    }
+  }
+
+  // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
+  *out << Verbose(5) << "Adding new triangle points."<< endl;
+  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
+  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
+  AddTesselationPoint(point, 2);
+  *out << Verbose(5) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  BTS->GetNormalVector(TempTriangle->NormalVector);
+  BTS->NormalVector.Scale(-1.);
+  *out << Verbose(3) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl;
+  AddTesselationTriangle();
+
+  // create other side of this triangle and close both new sides of the first created triangle
+  *out << Verbose(5) << "Adding new triangle points."<< endl;
+  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
+  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
+  AddTesselationPoint(point, 2);
+  *out << Verbose(5) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  BTS->GetNormalVector(TempTriangle->NormalVector);
+  *out << Verbose(3) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;
+  AddTesselationTriangle();
+
+  // add removed triangle to the last open line of the second triangle
+  for (int i=0;i<3;i++) { // look for the same line as BestLine (only it's its degenerated companion)
+    if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
+      if (BestLine == BTS->lines[i]){
+        *out << Verbose(1) << "CRITICAL: BestLine is same as found line, something's wrong here!" << endl;
+        exit(255);
+      }
+      BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) );
+      TempTriangle->lines[nr] = BTS->lines[i];
+      break;
+    }
+  }
+
+  // exit
+  *out << Verbose(2) << "End of AddBoundaryPointByDegeneratedTriangle" << endl;
+};
+
+/** Writes the envelope to file.
+ * \param *out otuput stream for debugging
+ * \param *filename basename of output file
+ * \param *cloud PointCloud structure with all nodes
+ */
+void Tesselation::Output(ofstream *out, const char *filename, PointCloud *cloud)
+{
+  ofstream *tempstream = NULL;
+  string NameofTempFile;
+  char NumberName[255];
+
+  if (LastTriangle != NULL) {
+    sprintf(NumberName, "-%04d-%s_%s_%s", (int)TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);
+    if (DoTecplotOutput) {
+      string NameofTempFile(filename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(TecplotSuffix);
+      *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      WriteTecplotFile(out, tempstream, this, cloud, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+
+    if (DoRaster3DOutput) {
+      string NameofTempFile(filename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(Raster3DSuffix);
+      *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      WriteRaster3dFile(out, tempstream, this, cloud);
+      IncludeSphereinRaster3D(out, tempstream, this, cloud);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+  }
+  if (DoTecplotOutput || DoRaster3DOutput)
+    TriangleFilesWritten++;
+};

src/tesselation.hpp

@@ -34 +34 @@
 class PointCloud;
 class Tesselation;
+
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 1
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
 
 // ======================================================= some template functions =========================================
@@ -115 +122 @@
 
     void GetNormalVector(Vector &NormalVector);
+    void GetCenter(Vector *center);
     bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
     bool ContainsBoundaryLine(class BoundaryLineSet *line);
@@ -120 +128 @@
     class BoundaryPointSet *GetThirdEndpoint(class BoundaryLineSet *line);
     bool IsPresentTupel(class BoundaryPointSet *Points[3]);
-    void GetCenter(Vector *center);
+    bool IsPresentTupel(class BoundaryTriangleSet *T);
 
     class BoundaryPointSet *endpoints[3];
@@ -200 +208 @@
     void RemoveTesselationPoint(class BoundaryPointSet *point);
 
-    bool IsInside(Vector *pointer);
     class BoundaryPointSet *GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2);
 
     // concave envelope
     void FindStartingTriangle(ofstream *out, const double RADIUS, class LinkedCell *LC);
-    void FindSecondPointForTesselation(class TesselPoint* a, class TesselPoint* Candidate, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, class LinkedCell *LC);
+    void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, class LinkedCell *LC);
     void FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, class LinkedCell *LC);
-    bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC);
+    bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC);
     int CheckPresenceOfTriangle(ofstream *out, class TesselPoint *Candidates[3]);
     class BoundaryTriangleSet * GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3]);
@@ -216 +223 @@
     bool InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC);
     double RemovePointFromTesselatedSurface(ofstream *out, class BoundaryPointSet *point);
-    bool FlipBaseline(ofstream *out, class BoundaryLineSet *Base);
-    bool PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base);
+    class BoundaryLineSet * FlipBaseline(ofstream *out, class BoundaryLineSet *Base);
+    double PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base);
     class BoundaryPointSet *IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base);
-    map<int, int> FindAllDegeneratedTriangles();
+    map<int, int> * FindAllDegeneratedTriangles();
+    map<int, int> * FindAllDegeneratedLines();
     void RemoveDegeneratedTriangles();
+    void AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC);
 
     set<TesselPoint*> * GetAllConnectedPoints(ofstream *out, TesselPoint* Point);
@@ -239 +248 @@
     void PrintAllBoundaryTriangles(ofstream *out);
 
+    // store envelope in file
+    void Output(ofstream *out, const char *filename, PointCloud *cloud);
 
     PointMap PointsOnBoundary;
@@ -261 +272 @@
     class BoundaryLineSet *BLS[3];
     class BoundaryTriangleSet *BTS;
+    class BoundaryTriangleSet *LastTriangle;
+    int TriangleFilesWritten;
 
   private:
@@ -268 +281 @@
 };
 
-bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3]);
-bool SortCandidates(class CandidateForTesselation* candidate1, class CandidateForTesselation* candidate2);
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC);
-TesselPoint* FindSecondClosestPoint(const Vector*, LinkedCell*);
-double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector);
-Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase);
 
 #endif /* TESSELATION_HPP_ */

src/tesselationhelpers.cpp

@@ -388 +388 @@
 };
 
+/** Gets the angle between a point and a reference relative to the provided center.
+ * We have two shanks point and reference between which the angle is calculated
+ * and by scalar product with OrthogonalVector we decide the interval.
+ * @param point to calculate the angle for
+ * @param reference to which to calculate the angle
+ * @param OrthogonalVector points in direction of [pi,2pi] interval
+ *
+ * @return angle between point and reference
+ */
+double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
+{
+  if (reference.IsZero())
+    return M_PI;
+
+  // calculate both angles and correct with in-plane vector
+  if (point.IsZero())
+    return M_PI;
+  double phi = point.Angle(&reference);
+  if (OrthogonalVector.ScalarProduct(&point) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+
+  cout << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
+
+  return phi;
+}
+
 
 /** Calculates the volume of a general tetraeder.
@@ -412 +439 @@
 };
 
+
+/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
+ * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
+ * make it bigger (i.e. closing one (the baseline) and opening two new ones).
+ * \param TPS[3] nodes of the triangle
+ * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
+ */
+bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3])
+{
+  bool result = false;
+  int counter = 0;
+
+  // check all three points
+  for (int i=0;i<3;i++)
+    for (int j=i+1; j<3; j++) {
+      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
+        LineMap::iterator FindLine;
+        pair<LineMap::iterator,LineMap::iterator> FindPair;
+        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
+        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+          // If there is a line with less than two attached triangles, we don't need a new line.
+          if (FindLine->second->triangles.size() < 2) {
+            counter++;
+            break;  // increase counter only once per edge
+          }
+        }
+      } else { // no line
+        cout << Verbose(1) << "The line between " << *nodes[i] << " and " << *nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
+        result = true;
+      }
+    }
+  if ((!result) && (counter > 1)) {
+    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
+    result = true;
+  }
+  return result;
+};
+
+
+/** Sort function for the candidate list.
+ */
+bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
+{
+  Vector BaseLineVector, OrthogonalVector, helper;
+  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
+    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
+    //return false;
+    exit(1);
+  }
+  // create baseline vector
+  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
+  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  BaseLineVector.Normalize();
+
+  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
+  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
+  helper.SubtractVector(candidate1->point->node);
+  OrthogonalVector.CopyVector(&helper);
+  helper.VectorProduct(&BaseLineVector);
+  OrthogonalVector.SubtractVector(&helper);
+  OrthogonalVector.Normalize();
+
+  // calculate both angles and correct with in-plane vector
+  helper.CopyVector(candidate1->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  double phi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  helper.CopyVector(candidate2->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  double psi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    psi = 2.*M_PI - psi;
+  }
+
+  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
+  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
+
+  // return comparison
+  return phi < psi;
+};
+
+/**
+ * Finds the point which is second closest to the provided one.
+ *
+ * @param Point to which to find the second closest other point
+ * @param linked cell structure
+ *
+ * @return point which is second closest to the provided one
+ */
+TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
+{
+  LinkedNodes *List = NULL;
+  TesselPoint* closestPoint = NULL;
+  TesselPoint* secondClosestPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //cout << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              // remember second point
+              secondDistance = distance;
+              secondClosestPoint = closestPoint;
+              // mark down new closest point
+              distance = currentNorm;
+              closestPoint = (*Runner);
+              //cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  return secondClosestPoint;
+};
+
+/**
+ * Finds the point which is closest to the provided one.
+ *
+ * @param Point to which to find the closest other point
+ * @param SecondPoint the second closest other point on return, NULL if none found
+ * @param linked cell structure
+ *
+ * @return point which is closest to the provided one, NULL if none found
+ */
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
+{
+  LinkedNodes *List = NULL;
+  TesselPoint* closestPoint = NULL;
+  SecondPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //cout << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              secondDistance = distance;
+              SecondPoint = closestPoint;
+              distance = currentNorm;
+              closestPoint = (*Runner);
+              //cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
+            } else if (currentNorm < secondDistance) {
+              secondDistance = currentNorm;
+              SecondPoint = (*Runner);
+              //cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  return closestPoint;
+};
+
+/** Returns the closest point on \a *Base with respect to \a *OtherBase.
+ * \param *out output stream for debugging
+ * \param *Base reference line
+ * \param *OtherBase other base line
+ * \return Vector on reference line that has closest distance
+ */
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
+{
+  // construct the plane of the two baselines (i.e. take both their directional vectors)
+  Vector Normal;
+  Vector Baseline, OtherBaseline;
+  Baseline.CopyVector(Base->endpoints[1]->node->node);
+  Baseline.SubtractVector(Base->endpoints[0]->node->node);
+  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
+  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
+  Normal.CopyVector(&Baseline);
+  Normal.VectorProduct(&OtherBaseline);
+  Normal.Normalize();
+  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
+
+  // project one offset point of OtherBase onto this plane (and add plane offset vector)
+  Vector NewOffset;
+  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
+  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
+  NewOffset.ProjectOntoPlane(&Normal);
+  NewOffset.AddVector(Base->endpoints[0]->node->node);
+  Vector NewDirection;
+  NewDirection.CopyVector(&NewOffset);
+  NewDirection.AddVector(&OtherBaseline);
+
+  // calculate the intersection between this projected baseline and Base
+  Vector *Intersection = new Vector;
+  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+  Normal.CopyVector(Intersection);
+  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
+
+  return Intersection;
+};
+
+
+/** Creates the objects in a VRML file.
+ * \param *out output stream for debugging
+ * \param *vrmlfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud)
+{
+  TesselPoint *Walker = NULL;
+  int i;
+  Vector *center = cloud->GetCenter(out);
+  if (vrmlfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *vrmlfile << "#VRML V2.0 utf8" << endl;
+    *vrmlfile << "#Created by molecuilder" << endl;
+    *vrmlfile << "#All atoms as spheres" << endl;
+    cloud->GoToFirst();
+    while (!cloud->IsEnd()) {
+      Walker = cloud->GetPoint();
+      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Walker->node->x[i]-center->x[i] << " ";
+      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+      cloud->GoToNext();
+    }
+
+    *vrmlfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *vrmlfile << "1" << endl << "  "; // 1 is triangle type
+      for (i=0;i<3;i++) { // print each node
+        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
+        *vrmlfile << "\t";
+      }
+      *vrmlfile << "1. 0. 0." << endl;  // red as colour
+      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+  }
+  delete(center);
+};
+
+/** Writes additionally the current sphere (i.e. the last triangle to file).
+ * \param *out output stream for debugging
+ * \param *rasterfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
+{
+  Vector helper;
+  // include the current position of the virtual sphere in the temporary raster3d file
+  Vector *center = cloud->GetCenter(out);
+  // make the circumsphere's center absolute again
+  helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
+  helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
+  helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
+  helper.Scale(1./3.);
+  helper.SubtractVector(center);
+  // and add to file plus translucency object
+  *rasterfile << "# current virtual sphere\n";
+  *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+  *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
+  *rasterfile << "9\n  terminating special property\n";
+  delete(center);
+};
+
+/** Creates the objects in a raster3d file (renderable with a header.r3d).
+ * \param *out output stream for debugging
+ * \param *rasterfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
+{
+  TesselPoint *Walker = NULL;
+  int i;
+  Vector *center = cloud->GetCenter(out);
+  if (rasterfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+    *rasterfile << "@header.r3d" << endl;
+    *rasterfile << "# All atoms as spheres" << endl;
+    cloud->GoToFirst();
+    while (!cloud->IsEnd()) {
+      Walker = cloud->GetPoint();
+      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Walker->node->x[i]-center->x[i] << " ";
+      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+      cloud->GoToNext();
+    }
+
+    *rasterfile << "# All tesselation triangles" << endl;
+    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
+      for (i=0;i<3;i++) {  // print each node
+        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
+        *rasterfile << "\t";
+      }
+      *rasterfile << "1. 0. 0." << endl;  // red as colour
+      //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
+    }
+    *rasterfile << "9\n#  terminating special property\n";
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  IncludeSphereinRaster3D(out, rasterfile, Tess, cloud);
+  delete(center);
+};
+
+/** This function creates the tecplot file, displaying the tesselation of the hull.
+ * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ * \param N arbitrary number to differentiate various zones in the tecplot format
+ */
+void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N)
+{
+  if ((tecplot != NULL) && (TesselStruct != NULL)) {
+    // write header
+    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
+    *tecplot << "ZONE T=\"" << N << "-";
+    for (int i=0;i<3;i++)
+      *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
+    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+    int i=0;
+    for (cloud->GoToFirst(); !cloud->IsEnd(); cloud->GoToNext(), i++);
+    int *LookupList = new int[i];
+    for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)
+      LookupList[i] = -1;
+
+    // print atom coordinates
+    *out << Verbose(2) << "The following triangles were created:";
+    int Counter = 1;
+    TesselPoint *Walker = NULL;
+    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
+      Walker = target->second->node;
+      LookupList[Walker->nr] = Counter++;
+      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
+    }
+    *tecplot << endl;
+    // print connectivity
+    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
+      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
+      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+    }
+    delete[] (LookupList);
+    *out << endl;
+  }
+};

src/tesselationhelpers.hpp

@@ -25 +25 @@
 
 #include "defs.hpp"
+#include "tesselation.hpp"
 #include "vector.hpp"
 
@@ -37 +38 @@
 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4);
 double CalculateVolumeofGeneralTetraeder(Vector *a, Vector *b, Vector *c, Vector *d);
+double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector);
+
+bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3]);
+bool SortCandidates(class CandidateForTesselation* candidate1, class CandidateForTesselation* candidate2);
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC);
+TesselPoint* FindSecondClosestPoint(const Vector*, LinkedCell*);
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase);
+
+void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N);
+void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud);
+void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud);
+void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud);