Changeset 57066a for src/boundary.cpp

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

File:

• src/boundary.cpp (modified) (15 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.