Changeset ab1932


Ignore:
Timestamp:
Aug 3, 2009, 8:21:05 PM (15 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
03e57a
Parents:
0dbddc (diff), edb93c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'TesselationRefactoring' into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/molecules.hpp

All of Saskia Metzler's new function were transfered from boundary.cpp to tesselation.cpp and the changes due to TesselPoint, LinkedCell and so on incorporated.

Location:
src
Files:
9 added
22 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    r0dbddc rab1932  
    1 SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp vector.cpp verbose.cpp
    2 HEADER = boundary.hpp defs.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp vector.hpp
     1SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
     2HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
    33
    44bin_PROGRAMS = molecuilder joiner analyzer

  • src/atom.cpp

    r0dbddc rab1932  
    55 */
    66
    7 #include "molecules.hpp"
     7#include "atom.hpp"
    88
    99/************************************* Functions for class atom *************************************/
    … …  
    1313atom::atom()
    1414{
    15   Name = NULL;
     15  father = this;  // generally, father is itself
    1616  previous = NULL;
    1717  next = NULL;
    18   father = this;  // generally, father is itself
    1918  Ancestor = NULL;
    2019  type = NULL;
    2120  sort = NULL;
    2221  FixedIon = 0;
    23   nr = -1;
    2422  GraphNr = -1;
    2523  ComponentNr = NULL;
    … …  
    2927  AdaptiveOrder = 0;
    3028  MaxOrder = false;
     29  // set LCNode::Vector to our Vector
     30  node = &x;
    3131};
    3232
    3333/** Constructor of class atom.
    34  *
    3534 */
    3635atom::atom(atom *pointer)
    … …  
    4140  father = this;  // generally, father is itself
    4241  Ancestor = NULL;
     42  GraphNr = -1;
     43  ComponentNr = NULL;
     44  IsCyclic = false;
     45  SeparationVertex = false;
     46  LowpointNr = -1;
     47  AdaptiveOrder = 0;
     48  MaxOrder = false;
    4349  type = pointer->type;  // copy element of atom
    4450  x.CopyVector(&pointer->x); // copy coordination
    … …  
    5460atom::~atom()
    5561{
    56   Free((void **)&Name, "atom::~atom: *Name");
    5762  Free((void **)&ComponentNr, "atom::~atom: *ComponentNr");
    5863};

  • src/bond.cpp

    r0dbddc rab1932  
    55 */
    66
    7 #include "molecules.hpp"
    8 
     7#include "bond.hpp"
    98
    109/***************************************** Functions for class bond ********************************/

  • src/boundary.cpp

    r0dbddc rab1932  
     1/** \file boundary.hpp
     2 *
     3 * Implementations and super-function for envelopes
     4 */
     5
     6
    17#include "boundary.hpp"
    2 #include "linkedcell.hpp"
    3 #include "molecules.hpp"
    4 #include <gsl/gsl_matrix.h>
    5 #include <gsl/gsl_linalg.h>
    6 #include <gsl/gsl_multimin.h>
    7 #include <gsl/gsl_permutation.h>
    8 
    9 #define DEBUG 1
    10 #define DoSingleStepOutput 0
    11 #define DoTecplotOutput 1
    12 #define DoRaster3DOutput 1
    13 #define DoVRMLOutput 1
    14 #define TecplotSuffix ".dat"
    15 #define Raster3DSuffix ".r3d"
    16 #define VRMLSUffix ".wrl"
    17 #define HULLEPSILON 1e-7
    18 
    19 // ======================================== Points on Boundary =================================
    20 
    21 BoundaryPointSet::BoundaryPointSet()
    22 {
    23   LinesCount = 0;
    24   Nr = -1;
    25 }
    26 ;
    27 
    28 BoundaryPointSet::BoundaryPointSet(atom *Walker)
    29 {
    30   node = Walker;
    31   LinesCount = 0;
    32   Nr = Walker->nr;
    33 }
    34 ;
    35 
    36 BoundaryPointSet::~BoundaryPointSet()
    37 {
    38   cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
    39   if (!lines.empty())
    40     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
    41   node = NULL;
    42 }
    43 ;
    44 
    45 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
    46 {
    47   cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
    48       << endl;
    49   if (line->endpoints[0] == this)
    50     {
    51       lines.insert(LinePair(line->endpoints[1]->Nr, line));
    52     }
    53   else
    54     {
    55       lines.insert(LinePair(line->endpoints[0]->Nr, line));
    56     }
    57   LinesCount++;
    58 }
    59 ;
    60 
    61 ostream &
    62 operator <<(ostream &ost, BoundaryPointSet &a)
    63 {
    64   ost << "[" << a.Nr << "|" << a.node->Name << "]";
    65   return ost;
    66 }
    67 ;
    68 
    69 // ======================================== Lines on Boundary =================================
    70 
    71 BoundaryLineSet::BoundaryLineSet()
    72 {
    73   for (int i = 0; i < 2; i++)
    74     endpoints[i] = NULL;
    75   TrianglesCount = 0;
    76   Nr = -1;
    77 }
    78 ;
    79 
    80 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
    81 {
    82   // set number
    83   Nr = number;
    84   // set endpoints in ascending order
    85   SetEndpointsOrdered(endpoints, Point[0], Point[1]);
    86   // add this line to the hash maps of both endpoints
    87   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
    88   Point[1]->AddLine(this); //
    89   // clear triangles list
    90   TrianglesCount = 0;
    91   cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
    92 }
    93 ;
    94 
    95 BoundaryLineSet::~BoundaryLineSet()
    96 {
    97   int Numbers[2];
    98   Numbers[0] = endpoints[1]->Nr;
    99   Numbers[1] = endpoints[0]->Nr;
    100   for (int i = 0; i < 2; i++) {
    101     cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
    102     // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
    103     pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
    104     for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
    105       if ((*Runner).second == this) {
    106         endpoints[i]->lines.erase(Runner);
    107         break;
    108       }
    109     if (endpoints[i]->lines.empty()) {
    110       cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
    111       if (endpoints[i] != NULL) {
    112         delete(endpoints[i]);
    113         endpoints[i] = NULL;
    114       } else
    115         cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
    116     } else
    117       cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
    118   }
    119   if (!triangles.empty())
    120     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
    121 }
    122 ;
    123 
    124 void
    125 BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
    126 {
    127   cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
    128       << endl;
    129   triangles.insert(TrianglePair(triangle->Nr, triangle));
    130   TrianglesCount++;
    131 }
    132 ;
    133 
    134 /** Checks whether we have a common endpoint with given \a *line.
    135  * \param *line other line to test
    136  * \return true - common endpoint present, false - not connected
    137  */
    138 bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
    139 {
    140   if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
    141     return true;
    142   else
    143     return false;
    144 };
    145 
    146 /** Checks whether the adjacent triangles of a baseline are convex or not.
    147  * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
    148  * If greater/equal M_PI than we are convex.
    149  * \param *out output stream for debugging
    150  * \return true - triangles are convex, false - concave or less than two triangles connected
    151  */
    152 bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
    153 {
    154   Vector BaseLineNormal;
    155   double angle = 0;
    156   // get the two triangles
    157   if (TrianglesCount != 2) {
    158     *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
    159     return false;
    160   }
    161   // have a normal vector on the base line pointing outwards
    162   BaseLineNormal.Zero();
    163   for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
    164     BaseLineNormal.AddVector(&runner->second->NormalVector);
    165   BaseLineNormal.Normalize();
    166   // and calculate the sum of the angles with this normal vector and each of the triangle ones'
    167   for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
    168     angle += BaseLineNormal.Angle(&runner->second->NormalVector);
    169 
    170   if ((angle - M_PI) > -MYEPSILON)
    171     return true;
    172   else
    173     return false;
    174 }
    175 
    176 /** Checks whether point is any of the two endpoints this line contains.
    177  * \param *point point to test
    178  * \return true - point is of the line, false - is not
    179  */
    180 bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
    181 {
    182   for(int i=0;i<2;i++)
    183     if (point == endpoints[i])
    184       return true;
    185   return false;
    186 };
    187 
    188 ostream &
    189 operator <<(ostream &ost, BoundaryLineSet &a)
    190 {
    191   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
    192       << a.endpoints[1]->node->Name << "]";
    193   return ost;
    194 }
    195 ;
    196 
    197 // ======================================== Triangles on Boundary =================================
    198 
    199 
    200 BoundaryTriangleSet::BoundaryTriangleSet()
    201 {
    202   for (int i = 0; i < 3; i++)
    203     {
    204       endpoints[i] = NULL;
    205       lines[i] = NULL;
    206     }
    207   Nr = -1;
    208 }
    209 ;
    210 
    211 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
    212 {
    213   // set number
    214   Nr = number;
    215   // set lines
    216   cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
    217   for (int i = 0; i < 3; i++)
    218     {
    219       lines[i] = line[i];
    220       lines[i]->AddTriangle(this);
    221     }
    222   // get ascending order of endpoints
    223   map<int, class BoundaryPointSet *> OrderMap;
    224   for (int i = 0; i < 3; i++)
    225     // for all three lines
    226     for (int j = 0; j < 2; j++)
    227       { // for both endpoints
    228         OrderMap.insert(pair<int, class BoundaryPointSet *> (
    229             line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
    230         // and we don't care whether insertion fails
    231       }
    232   // set endpoints
    233   int Counter = 0;
    234   cout << Verbose(6) << " with end points ";
    235   for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
    236       != OrderMap.end(); runner++)
    237     {
    238       endpoints[Counter] = runner->second;
    239       cout << " " << *endpoints[Counter];
    240       Counter++;
    241     }
    242   if (Counter < 3)
    243     {
    244       cerr << "ERROR! We have a triangle with only two distinct endpoints!"
    245           << endl;
    246       //exit(1);
    247     }
    248   cout << "." << endl;
    249 }
    250 ;
    251 
    252 BoundaryTriangleSet::~BoundaryTriangleSet()
    253 {
    254   for (int i = 0; i < 3; i++) {
    255     cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
    256     lines[i]->triangles.erase(Nr);
    257     if (lines[i]->triangles.empty()) {
    258       if (lines[i] != NULL) {
    259         cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
    260         delete (lines[i]);
    261         lines[i] = NULL;
    262       } else
    263         cerr << "ERROR: This line " << i << " has already been free'd." << endl;
    264     } else
    265       cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
    266   }
    267 }
    268 ;
    269 
    270 /** Calculates the normal vector for this triangle.
    271  * Is made unique by comparison with \a OtherVector to point in the other direction.
    272  * \param &OtherVector direction vector to make normal vector unique.
    273  */
    274 void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
    275 {
    276   // get normal vector
    277   NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
    278       &endpoints[2]->node->x);
    279 
    280   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
    281   if (NormalVector.Projection(&OtherVector) > 0)
    282     NormalVector.Scale(-1.);
    283 };
    284 
    285 /** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
    286  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
    287  * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
    288  * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
    289  * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
    290  * smaller than the first line.
    291  * \param *out output stream for debugging
    292  * \param *MolCenter offset vector of line
    293  * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
    294  * \param *Intersection intersection on plane on return
    295  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
    296  */
    297 bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
    298 {
    299   Vector CrossPoint;
    300   Vector helper;
    301   int i=0;
    302 
    303   if (Intersection->GetIntersectionWithPlane(out, &NormalVector, &endpoints[0]->node->x, MolCenter, x)) {
    304     *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
    305     return false;
    306   }
    307 
    308   // Calculate cross point between one baseline and the line from the third endpoint to intersection
    309   do {
    310     CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, &endpoints[i%3]->node->x, &endpoints[(i+1)%3]->node->x, &endpoints[(i+2)%3]->node->x, Intersection);
    311     helper.CopyVector(&endpoints[(i+1)%3]->node->x);
    312     helper.SubtractVector(&endpoints[i%3]->node->x);
    313     i++;
    314     if (i>3)
    315       break;
    316   } while (CrossPoint.NormSquared() < MYEPSILON);
    317   if (i>3) {
    318     *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
    319     exit(255);
    320   }
    321   CrossPoint.SubtractVector(&endpoints[i%3]->node->x);
    322 
    323   // check whether intersection is inside or not by comparing length of intersection and length of cross point
    324   if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
    325     return true;
    326   } else { // outside!
    327     Intersection->Zero();
    328     return false;
    329   }
    330 };
    331 
    332 /** Checks whether lines is any of the three boundary lines this triangle contains.
    333  * \param *line line to test
    334  * \return true - line is of the triangle, false - is not
    335  */
    336 bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
    337 {
    338   for(int i=0;i<3;i++)
    339     if (line == lines[i])
    340       return true;
    341   return false;
    342 };
    343 
    344 /** Checks whether point is any of the three endpoints this triangle contains.
    345  * \param *point point to test
    346  * \return true - point is of the triangle, false - is not
    347  */
    348 bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
    349 {
    350   for(int i=0;i<3;i++)
    351     if (point == endpoints[i])
    352       return true;
    353   return false;
    354 };
    355 
    356 ostream &
    357 operator <<(ostream &ost, BoundaryTriangleSet &a)
    358 {
    359   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
    360       << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
    361   return ost;
    362 }
    363 ;
    364 
    365 
    366 // ============================ CandidateForTesselation =============================
    367 
    368 CandidateForTesselation::CandidateForTesselation(
    369   atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
    370 ) {
    371   point = candidate;
    372   BaseLine = line;
    373   OptCenter.CopyVector(&OptCandidateCenter);
    374   OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
    375 }
    376 
    377 CandidateForTesselation::~CandidateForTesselation() {
    378   point = NULL;
    379   BaseLine = NULL;
    380 }
     8
     9#include<gsl/gsl_poly.h>
    38110
    38211// ========================================== F U N C T I O N S =================================
    38312
    384 /** Finds the endpoint two lines are sharing.
    385  * \param *line1 first line
    386  * \param *line2 second line
    387  * \return point which is shared or NULL if none
    388  */
    389 class BoundaryPointSet *
    390 GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
    391 {
    392   class BoundaryLineSet * lines[2] =
    393     { line1, line2 };
    394   class BoundaryPointSet *node = NULL;
    395   map<int, class BoundaryPointSet *> OrderMap;
    396   pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
    397   for (int i = 0; i < 2; i++)
    398     // for both lines
    399     for (int j = 0; j < 2; j++)
    400       { // for both endpoints
    401         OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
    402             lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
    403         if (!OrderTest.second)
    404           { // if insertion fails, we have common endpoint
    405             node = OrderTest.first->second;
    406             cout << Verbose(5) << "Common endpoint of lines " << *line1
    407                 << " and " << *line2 << " is: " << *node << "." << endl;
    408             j = 2;
    409             i = 2;
    410             break;
    411           }
    412       }
    413   return node;
    414 }
    415 ;
    416 
    417 /** Determines the boundary points of a cluster.
    418  * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
    419  * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
    420  * center and first and last point in the triple, it is thrown out.
    421  * \param *out output stream for debugging
    422  * \param *mol molecule structure representing the cluster
    423  */
    424 Boundaries *
    425 GetBoundaryPoints(ofstream *out, molecule *mol)
    426 {
    427   atom *Walker = NULL;
    428   PointMap PointsOnBoundary;
    429   LineMap LinesOnBoundary;
    430   TriangleMap TrianglesOnBoundary;
    431   Vector *MolCenter = mol->DetermineCenterOfAll(out);
    432   Vector helper;
    433 
    434   *out << Verbose(1) << "Finding all boundary points." << endl;
    435   Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
    436   BoundariesTestPair BoundaryTestPair;
    437   Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
    438   double radius, angle;
    439   // 3a. Go through every axis
    440   for (int axis = 0; axis < NDIM; axis++) {
    441     AxisVector.Zero();
    442     AngleReferenceVector.Zero();
    443     AngleReferenceNormalVector.Zero();
    444     AxisVector.x[axis] = 1.;
    445     AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
    446     AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
    447 
    448     *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
    449 
    450     // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    451     Walker = mol->start;
    452     while (Walker->next != mol->end) {
    453       Walker = Walker->next;
    454       Vector ProjectedVector;
    455       ProjectedVector.CopyVector(&Walker->x);
    456       ProjectedVector.SubtractVector(MolCenter);
    457       ProjectedVector.ProjectOntoPlane(&AxisVector);
    458 
    459       // correct for negative side
    460       radius = ProjectedVector.NormSquared();
    461       if (fabs(radius) > MYEPSILON)
    462         angle = ProjectedVector.Angle(&AngleReferenceVector);
    463       else
    464         angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
    465 
    466       //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
    467       if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
    468         angle = 2. * M_PI - angle;
    469       }
    470       *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
    471       BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
    472       if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
    473         *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
    474         *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
    475         *out << Verbose(2) << "New vector: " << *Walker << endl;
    476         double tmp = ProjectedVector.NormSquared();
    477         if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
    478           BoundaryTestPair.first->second.first = tmp;
    479           BoundaryTestPair.first->second.second = Walker;
    480           *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
    481         } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
    482           helper.CopyVector(&Walker->x);
    483           helper.SubtractVector(MolCenter);
    484           tmp = helper.NormSquared();
    485           helper.CopyVector(&BoundaryTestPair.first->second.second->x);
    486           helper.SubtractVector(MolCenter);
    487           if (helper.NormSquared() < tmp) {
    488             BoundaryTestPair.first->second.second = Walker;
    489             *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
    490           } else {
    491             *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
    492           }
    493         } else {
    494           *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
    495         }
    496       }
    497     }
    498     // printing all inserted for debugging
    499     //    {
    500     //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
    501     //      int i=0;
    502     //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    503     //        if (runner != BoundaryPoints[axis].begin())
    504     //          *out << ", " << i << ": " << *runner->second.second;
    505     //        else
    506     //          *out << i << ": " << *runner->second.second;
    507     //        i++;
    508     //      }
    509     //      *out << endl;
    510     //    }
    511     // 3c. throw out points whose distance is less than the mean of left and right neighbours
    512     bool flag = false;
    513     *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
    514     do { // do as long as we still throw one out per round
    515       flag = false;
    516       Boundaries::iterator left = BoundaryPoints[axis].end();
    517       Boundaries::iterator right = BoundaryPoints[axis].end();
    518       for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    519         // set neighbours correctly
    520         if (runner == BoundaryPoints[axis].begin()) {
    521           left = BoundaryPoints[axis].end();
    522         } else {
    523           left = runner;
    524         }
    525         left--;
    526         right = runner;
    527         right++;
    528         if (right == BoundaryPoints[axis].end()) {
    529           right = BoundaryPoints[axis].begin();
    530         }
    531         // check distance
    532 
    533         // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
    534         {
    535           Vector SideA, SideB, SideC, SideH;
    536           SideA.CopyVector(&left->second.second->x);
    537           SideA.SubtractVector(MolCenter);
    538           SideA.ProjectOntoPlane(&AxisVector);
    539           //          *out << "SideA: ";
    540           //          SideA.Output(out);
    541           //          *out << endl;
    542 
    543           SideB.CopyVector(&right->second.second->x);
    544           SideB.SubtractVector(MolCenter);
    545           SideB.ProjectOntoPlane(&AxisVector);
    546           //          *out << "SideB: ";
    547           //          SideB.Output(out);
    548           //          *out << endl;
    549 
    550           SideC.CopyVector(&left->second.second->x);
    551           SideC.SubtractVector(&right->second.second->x);
    552           SideC.ProjectOntoPlane(&AxisVector);
    553           //          *out << "SideC: ";
    554           //          SideC.Output(out);
    555           //          *out << endl;
    556 
    557           SideH.CopyVector(&runner->second.second->x);
    558           SideH.SubtractVector(MolCenter);
    559           SideH.ProjectOntoPlane(&AxisVector);
    560           //          *out << "SideH: ";
    561           //          SideH.Output(out);
    562           //          *out << endl;
    563 
    564           // calculate each length
    565           double a = SideA.Norm();
    566           //double b = SideB.Norm();
    567           //double c = SideC.Norm();
    568           double h = SideH.Norm();
    569           // calculate the angles
    570           double alpha = SideA.Angle(&SideH);
    571           double beta = SideA.Angle(&SideC);
    572           double gamma = SideB.Angle(&SideH);
    573           double delta = SideC.Angle(&SideH);
    574           double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
    575           //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
    576           *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
    577           if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
    578             // throw out point
    579             *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
    580             BoundaryPoints[axis].erase(runner);
    581             flag = true;
    582           }
    583         }
    584       }
    585     } while (flag);
    586   }
    587   delete(MolCenter);
    588   return BoundaryPoints;
    589 }
    590 ;
    59113
    59214/** Determines greatest diameters of a cluster defined by its convex envelope.
    … …  
    60527  bool BoundaryFreeFlag = false;
    60628  Boundaries *BoundaryPoints = BoundaryPtr;
    607   if (BoundaryPoints == NULL)
    608     {
    609       BoundaryFreeFlag = true;
    610       BoundaryPoints = GetBoundaryPoints(out, mol);
    611     }
    612   else
    613     {
    614       *out << Verbose(1) << "Using given boundary points set." << endl;
    615     }
     29  if (BoundaryPoints == NULL) {
     30    BoundaryFreeFlag = true;
     31    BoundaryPoints = GetBoundaryPoints(out, mol);
     32  } else {
     33    *out << Verbose(1) << "Using given boundary points set." << endl;
     34  }
    61635  // determine biggest "diameter" of cluster for each axis
    61736  Boundaries::iterator Neighbour, OtherNeighbour;
    … …  
    735154      for (i=0;i<3;i++) { // print each node
    736155        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    737           *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
     156          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
    738157        *vrmlfile << "\t";
    739158      }
    … …  
    790209      for (i=0;i<3;i++) {  // print each node
    791210        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    792           *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
     211          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
    793212        *rasterfile << "\t";
    794213      }
    … …  
    821240    *out << Verbose(2) << "The following triangles were created:";
    822241    int Counter = 1;
    823     atom *Walker = NULL;
     242    TesselPoint *Walker = NULL;
    824243    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
    825244      Walker = target->second->node;
    826245      LookupList[Walker->nr] = Counter++;
    827       *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
     246      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
    828247    }
    829248    *tecplot << endl;
    … …  
    837256  }
    838257}
     258
     259
     260/** Determines the boundary points of a cluster.
     261 * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
     262 * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
     263 * center and first and last point in the triple, it is thrown out.
     264 * \param *out output stream for debugging
     265 * \param *mol molecule structure representing the cluster
     266 */
     267Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
     268{
     269  atom *Walker = NULL;
     270  PointMap PointsOnBoundary;
     271  LineMap LinesOnBoundary;
     272  TriangleMap TrianglesOnBoundary;
     273  Vector *MolCenter = mol->DetermineCenterOfAll(out);
     274  Vector helper;
     275
     276  *out << Verbose(1) << "Finding all boundary points." << endl;
     277  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
     278  BoundariesTestPair BoundaryTestPair;
     279  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
     280  double radius, angle;
     281  // 3a. Go through every axis
     282  for (int axis = 0; axis < NDIM; axis++) {
     283    AxisVector.Zero();
     284    AngleReferenceVector.Zero();
     285    AngleReferenceNormalVector.Zero();
     286    AxisVector.x[axis] = 1.;
     287    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
     288    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
     289
     290    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
     291
     292    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
     293    Walker = mol->start;
     294    while (Walker->next != mol->end) {
     295      Walker = Walker->next;
     296      Vector ProjectedVector;
     297      ProjectedVector.CopyVector(&Walker->x);
     298      ProjectedVector.SubtractVector(MolCenter);
     299      ProjectedVector.ProjectOntoPlane(&AxisVector);
     300
     301      // correct for negative side
     302      radius = ProjectedVector.NormSquared();
     303      if (fabs(radius) > MYEPSILON)
     304        angle = ProjectedVector.Angle(&AngleReferenceVector);
     305      else
     306        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
     307
     308      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
     309      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
     310        angle = 2. * M_PI - angle;
     311      }
     312      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
     313      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
     314      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
     315        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
     316        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
     317        *out << Verbose(2) << "New vector: " << *Walker << endl;
     318        double tmp = ProjectedVector.NormSquared();
     319        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
     320          BoundaryTestPair.first->second.first = tmp;
     321          BoundaryTestPair.first->second.second = Walker;
     322          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
     323        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
     324          helper.CopyVector(&Walker->x);
     325          helper.SubtractVector(MolCenter);
     326          tmp = helper.NormSquared();
     327          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
     328          helper.SubtractVector(MolCenter);
     329          if (helper.NormSquared() < tmp) {
     330            BoundaryTestPair.first->second.second = Walker;
     331            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
     332          } else {
     333            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
     334          }
     335        } else {
     336          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
     337        }
     338      }
     339    }
     340    // printing all inserted for debugging
     341    //    {
     342    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     343    //      int i=0;
     344    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
     345    //        if (runner != BoundaryPoints[axis].begin())
     346    //          *out << ", " << i << ": " << *runner->second.second;
     347    //        else
     348    //          *out << i << ": " << *runner->second.second;
     349    //        i++;
     350    //      }
     351    //      *out << endl;
     352    //    }
     353    // 3c. throw out points whose distance is less than the mean of left and right neighbours
     354    bool flag = false;
     355    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
     356    do { // do as long as we still throw one out per round
     357      flag = false;
     358      Boundaries::iterator left = BoundaryPoints[axis].end();
     359      Boundaries::iterator right = BoundaryPoints[axis].end();
     360      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
     361        // set neighbours correctly
     362        if (runner == BoundaryPoints[axis].begin()) {
     363          left = BoundaryPoints[axis].end();
     364        } else {
     365          left = runner;
     366        }
     367        left--;
     368        right = runner;
     369        right++;
     370        if (right == BoundaryPoints[axis].end()) {
     371          right = BoundaryPoints[axis].begin();
     372        }
     373        // check distance
     374
     375        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
     376        {
     377          Vector SideA, SideB, SideC, SideH;
     378          SideA.CopyVector(&left->second.second->x);
     379          SideA.SubtractVector(MolCenter);
     380          SideA.ProjectOntoPlane(&AxisVector);
     381          //          *out << "SideA: ";
     382          //          SideA.Output(out);
     383          //          *out << endl;
     384
     385          SideB.CopyVector(&right->second.second->x);
     386          SideB.SubtractVector(MolCenter);
     387          SideB.ProjectOntoPlane(&AxisVector);
     388          //          *out << "SideB: ";
     389          //          SideB.Output(out);
     390          //          *out << endl;
     391
     392          SideC.CopyVector(&left->second.second->x);
     393          SideC.SubtractVector(&right->second.second->x);
     394          SideC.ProjectOntoPlane(&AxisVector);
     395          //          *out << "SideC: ";
     396          //          SideC.Output(out);
     397          //          *out << endl;
     398
     399          SideH.CopyVector(&runner->second.second->x);
     400          SideH.SubtractVector(MolCenter);
     401          SideH.ProjectOntoPlane(&AxisVector);
     402          //          *out << "SideH: ";
     403          //          SideH.Output(out);
     404          //          *out << endl;
     405
     406          // calculate each length
     407          double a = SideA.Norm();
     408          //double b = SideB.Norm();
     409          //double c = SideC.Norm();
     410          double h = SideH.Norm();
     411          // calculate the angles
     412          double alpha = SideA.Angle(&SideH);
     413          double beta = SideA.Angle(&SideC);
     414          double gamma = SideB.Angle(&SideH);
     415          double delta = SideC.Angle(&SideH);
     416          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
     417          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
     418          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
     419          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
     420            // throw out point
     421            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
     422            BoundaryPoints[axis].erase(runner);
     423            flag = true;
     424          }
     425        }
     426      }
     427    } while (flag);
     428  }
     429  delete(MolCenter);
     430  return BoundaryPoints;
     431};
    839432
    840433/** Tesselates the convex boundary by finding all boundary points.
    … …  
    959552      << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
    960553      << endl;
    961   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
    962       != TesselStruct->TrianglesOnBoundary.end(); runner++)
     554  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
    963555    { // go through every triangle, calculate volume of its pyramid with CoG as peak
    964       x.CopyVector(&runner->second->endpoints[0]->node->x);
    965       x.SubtractVector(&runner->second->endpoints[1]->node->x);
    966       y.CopyVector(&runner->second->endpoints[0]->node->x);
    967       y.SubtractVector(&runner->second->endpoints[2]->node->x);
    968       a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
    969           &runner->second->endpoints[1]->node->x));
    970       b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
    971           &runner->second->endpoints[2]->node->x));
    972       c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
    973           &runner->second->endpoints[1]->node->x));
     556      x.CopyVector(runner->second->endpoints[0]->node->node);
     557      x.SubtractVector(runner->second->endpoints[1]->node->node);
     558      y.CopyVector(runner->second->endpoints[0]->node->node);
     559      y.SubtractVector(runner->second->endpoints[2]->node->node);
     560      a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
     561      b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
     562      c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
    974563      G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
    975       x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
    976           &runner->second->endpoints[1]->node->x,
    977           &runner->second->endpoints[2]->node->x);
    978       x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
     564      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
     565      x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
    979566      h = x.Norm(); // distance of CoG to triangle
    980567      PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
    … …  
    1164751        FillIt = true;
    1165752        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    1166           //FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
     753          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
    1167754        }
    1168755
    … …  
    1228815
    1229816
    1230 // =========================================================== class TESSELATION ===========================================
    1231 
    1232 /** Constructor of class Tesselation.
    1233  */
    1234 Tesselation::Tesselation()
    1235 {
    1236   PointsOnBoundaryCount = 0;
    1237   LinesOnBoundaryCount = 0;
    1238   TrianglesOnBoundaryCount = 0;
    1239   TriangleFilesWritten = 0;
    1240 }
    1241 ;
    1242 
    1243 /** Constructor of class Tesselation.
    1244  * We have to free all points, lines and triangles.
    1245  */
    1246 Tesselation::~Tesselation()
    1247 {
    1248   cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
    1249   for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
    1250     if (runner->second != NULL) {
    1251       delete (runner->second);
    1252       runner->second = NULL;
    1253     } else
    1254       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
    1255   }
    1256 }
    1257 ;
    1258 
    1259 /** Gueses first starting triangle of the convex envelope.
    1260  * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
    1261  * \param *out output stream for debugging
    1262  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
    1263  */
    1264 void
    1265 Tesselation::GuessStartingTriangle(ofstream *out)
    1266 {
    1267   // 4b. create a starting triangle
    1268   // 4b1. create all distances
    1269   DistanceMultiMap DistanceMMap;
    1270   double distance, tmp;
    1271   Vector PlaneVector, TrialVector;
    1272   PointMap::iterator A, B, C; // three nodes of the first triangle
    1273   A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
    1274 
    1275   // with A chosen, take each pair B,C and sort
    1276   if (A != PointsOnBoundary.end())
    1277     {
    1278       B = A;
    1279       B++;
    1280       for (; B != PointsOnBoundary.end(); B++)
    1281         {
    1282           C = B;
    1283           C++;
    1284           for (; C != PointsOnBoundary.end(); C++)
    1285             {
    1286               tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
    1287               distance = tmp * tmp;
    1288               tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
    1289               distance += tmp * tmp;
    1290               tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
    1291               distance += tmp * tmp;
    1292               DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
    1293                   PointMap::iterator, PointMap::iterator> (B, C)));
    1294             }
    1295         }
    1296     }
    1297   //    // listing distances
    1298   //    *out << Verbose(1) << "Listing DistanceMMap:";
    1299   //    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
    1300   //      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
    1301   //    }
    1302   //    *out << endl;
    1303   // 4b2. pick three baselines forming a triangle
    1304   // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
    1305   DistanceMultiMap::iterator baseline = DistanceMMap.begin();
    1306   for (; baseline != DistanceMMap.end(); baseline++)
    1307     {
    1308       // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
    1309       // 2. next, we have to check whether all points reside on only one side of the triangle
    1310       // 3. construct plane vector
    1311       PlaneVector.MakeNormalVector(&A->second->node->x,
    1312           &baseline->second.first->second->node->x,
    1313           &baseline->second.second->second->node->x);
    1314       *out << Verbose(2) << "Plane vector of candidate triangle is ";
    1315       PlaneVector.Output(out);
    1316       *out << endl;
    1317       // 4. loop over all points
    1318       double sign = 0.;
    1319       PointMap::iterator checker = PointsOnBoundary.begin();
    1320       for (; checker != PointsOnBoundary.end(); checker++)
    1321         {
    1322           // (neglecting A,B,C)
    1323           if ((checker == A) || (checker == baseline->second.first) || (checker
    1324               == baseline->second.second))
    1325             continue;
    1326           // 4a. project onto plane vector
    1327           TrialVector.CopyVector(&checker->second->node->x);
    1328           TrialVector.SubtractVector(&A->second->node->x);
    1329           distance = TrialVector.Projection(&PlaneVector);
    1330           if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
    1331             continue;
    1332           *out << Verbose(3) << "Projection of " << checker->second->node->Name
    1333               << " yields distance of " << distance << "." << endl;
    1334           tmp = distance / fabs(distance);
    1335           // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
    1336           if ((sign != 0) && (tmp != sign))
    1337             {
    1338               // 4c. If so, break 4. loop and continue with next candidate in 1. loop
    1339               *out << Verbose(2) << "Current candidates: "
    1340                   << A->second->node->Name << ","
    1341                   << baseline->second.first->second->node->Name << ","
    1342                   << baseline->second.second->second->node->Name << " leaves "
    1343                   << checker->second->node->Name << " outside the convex hull."
    1344                   << endl;
    1345               break;
    1346             }
    1347           else
    1348             { // note the sign for later
    1349               *out << Verbose(2) << "Current candidates: "
    1350                   << A->second->node->Name << ","
    1351                   << baseline->second.first->second->node->Name << ","
    1352                   << baseline->second.second->second->node->Name << " leave "
    1353                   << checker->second->node->Name << " inside the convex hull."
    1354                   << endl;
    1355               sign = tmp;
    1356             }
    1357           // 4d. Check whether the point is inside the triangle (check distance to each node
    1358           tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
    1359           int innerpoint = 0;
    1360           if ((tmp < A->second->node->x.DistanceSquared(
    1361               &baseline->second.first->second->node->x)) && (tmp
    1362               < A->second->node->x.DistanceSquared(
    1363                   &baseline->second.second->second->node->x)))
    1364             innerpoint++;
    1365           tmp = checker->second->node->x.DistanceSquared(
    1366               &baseline->second.first->second->node->x);
    1367           if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
    1368               &A->second->node->x)) && (tmp
    1369               < baseline->second.first->second->node->x.DistanceSquared(
    1370                   &baseline->second.second->second->node->x)))
    1371             innerpoint++;
    1372           tmp = checker->second->node->x.DistanceSquared(
    1373               &baseline->second.second->second->node->x);
    1374           if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
    1375               &baseline->second.first->second->node->x)) && (tmp
    1376               < baseline->second.second->second->node->x.DistanceSquared(
    1377                   &A->second->node->x)))
    1378             innerpoint++;
    1379           // 4e. If so, break 4. loop and continue with next candidate in 1. loop
    1380           if (innerpoint == 3)
    1381             break;
    1382         }
    1383       // 5. come this far, all on same side? Then break 1. loop and construct triangle
    1384       if (checker == PointsOnBoundary.end())
    1385         {
    1386           *out << "Looks like we have a candidate!" << endl;
    1387           break;
    1388         }
    1389     }
    1390   if (baseline != DistanceMMap.end())
    1391     {
    1392       BPS[0] = baseline->second.first->second;
    1393       BPS[1] = baseline->second.second->second;
    1394       BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    1395       BPS[0] = A->second;
    1396       BPS[1] = baseline->second.second->second;
    1397       BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    1398       BPS[0] = baseline->second.first->second;
    1399       BPS[1] = A->second;
    1400       BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    1401 
    1402       // 4b3. insert created triangle
    1403       BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    1404       TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    1405       TrianglesOnBoundaryCount++;
    1406       for (int i = 0; i < NDIM; i++)
    1407         {
    1408           LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
    1409           LinesOnBoundaryCount++;
    1410         }
    1411 
    1412       *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
    1413     }
    1414   else
    1415     {
    1416       *out << Verbose(1) << "No starting triangle found." << endl;
    1417       exit(255);
    1418     }
    1419 }
    1420 ;
    1421 
    1422 /** Tesselates the convex envelope of a cluster from a single starting triangle.
    1423  * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
    1424  * 2 triangles. Hence, we go through all current lines:
    1425  * -# if the lines contains to only one triangle
    1426  * -# We search all points in the boundary
    1427  *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
    1428  *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
    1429  *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
    1430  *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
    1431  * \param *out output stream for debugging
    1432  * \param *configuration for IsAngstroem
    1433  * \param *mol the cluster as a molecule structure
    1434  */
    1435 void Tesselation::TesselateOnBoundary(ofstream *out, molecule *mol)
    1436 {
    1437   bool flag;
    1438   PointMap::iterator winner;
    1439   class BoundaryPointSet *peak = NULL;
    1440   double SmallestAngle, TempAngle;
    1441   Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *MolCenter = NULL;
    1442   LineMap::iterator LineChecker[2];
    1443 
    1444   MolCenter = mol->DetermineCenterOfAll(out);
    1445   // create a first tesselation with the given BoundaryPoints
    1446   do {
    1447     flag = false;
    1448     for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
    1449       if (baseline->second->TrianglesCount == 1) {
    1450         // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
    1451         SmallestAngle = M_PI;
    1452 
    1453         // get peak point with respect to this base line's only triangle
    1454         BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
    1455         *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
    1456         for (int i = 0; i < 3; i++)
    1457           if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
    1458             peak = BTS->endpoints[i];
    1459         *out << Verbose(3) << " and has peak " << *peak << "." << endl;
    1460 
    1461         // prepare some auxiliary vectors
    1462         Vector BaseLineCenter, BaseLine;
    1463         BaseLineCenter.CopyVector(&baseline->second->endpoints[0]->node->x);
    1464         BaseLineCenter.AddVector(&baseline->second->endpoints[1]->node->x);
    1465         BaseLineCenter.Scale(1. / 2.); // points now to center of base line
    1466         BaseLine.CopyVector(&baseline->second->endpoints[0]->node->x);
    1467         BaseLine.SubtractVector(&baseline->second->endpoints[1]->node->x);
    1468 
    1469         // offset to center of triangle
    1470         CenterVector.Zero();
    1471         for (int i = 0; i < 3; i++)
    1472           CenterVector.AddVector(&BTS->endpoints[i]->node->x);
    1473         CenterVector.Scale(1. / 3.);
    1474         *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;
    1475 
    1476         // normal vector of triangle
    1477         NormalVector.CopyVector(MolCenter);
    1478         NormalVector.SubtractVector(&CenterVector);
    1479         BTS->GetNormalVector(NormalVector);
    1480         NormalVector.CopyVector(&BTS->NormalVector);
    1481         *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
    1482 
    1483         // vector in propagation direction (out of triangle)
    1484         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
    1485         PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
    1486         TempVector.CopyVector(&CenterVector);
    1487         TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
    1488         //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
    1489         if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
    1490           PropagationVector.Scale(-1.);
    1491         *out << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl;
    1492         winner = PointsOnBoundary.end();
    1493 
    1494         // loop over all points and calculate angle between normal vector of new and present triangle
    1495         for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
    1496           if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
    1497             *out << Verbose(3) << "Target point is " << *(target->second) << ":" << endl;
    1498 
    1499             // first check direction, so that triangles don't intersect
    1500             VirtualNormalVector.CopyVector(&target->second->node->x);
    1501             VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
    1502             VirtualNormalVector.ProjectOntoPlane(&NormalVector);
    1503             TempAngle = VirtualNormalVector.Angle(&PropagationVector);
    1504             *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
    1505             if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
    1506               *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
    1507               continue;
    1508             } else
    1509               *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
    1510 
    1511             // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
    1512             LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
    1513             LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
    1514             if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->TrianglesCount == 2))) {
    1515               *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
    1516               continue;
    1517             }
    1518             if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
    1519               *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
    1520               continue;
    1521             }
    1522 
    1523             // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
    1524             if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
    1525               *out << Verbose(4) << "Current target is peak!" << endl;
    1526               continue;
    1527             }
    1528 
    1529             // check for linear dependence
    1530             TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
    1531             TempVector.SubtractVector(&target->second->node->x);
    1532             helper.CopyVector(&baseline->second->endpoints[1]->node->x);
    1533             helper.SubtractVector(&target->second->node->x);
    1534             helper.ProjectOntoPlane(&TempVector);
    1535             if (fabs(helper.NormSquared()) < MYEPSILON) {
    1536               *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
    1537               continue;
    1538             }
    1539 
    1540             // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
    1541             flag = true;
    1542             VirtualNormalVector.MakeNormalVector(&baseline->second->endpoints[0]->node->x, &baseline->second->endpoints[1]->node->x, &target->second->node->x);
    1543             TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
    1544             TempVector.AddVector(&baseline->second->endpoints[1]->node->x);
    1545             TempVector.AddVector(&target->second->node->x);
    1546             TempVector.Scale(1./3.);
    1547             TempVector.SubtractVector(MolCenter);
    1548             // make it always point outward
    1549             if (VirtualNormalVector.Projection(&TempVector) < 0)
    1550               VirtualNormalVector.Scale(-1.);
    1551             // calculate angle
    1552             TempAngle = NormalVector.Angle(&VirtualNormalVector);
    1553             *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
    1554             if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
    1555               SmallestAngle = TempAngle;
    1556               winner = target;
    1557               *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
    1558             } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
    1559               // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
    1560               helper.CopyVector(&target->second->node->x);
    1561               helper.SubtractVector(&BaseLineCenter);
    1562               helper.ProjectOntoPlane(&BaseLine);
    1563               // ...the one with the smaller angle is the better candidate
    1564               TempVector.CopyVector(&target->second->node->x);
    1565               TempVector.SubtractVector(&BaseLineCenter);
    1566               TempVector.ProjectOntoPlane(&VirtualNormalVector);
    1567               TempAngle = TempVector.Angle(&helper);
    1568               TempVector.CopyVector(&winner->second->node->x);
    1569               TempVector.SubtractVector(&BaseLineCenter);
    1570               TempVector.ProjectOntoPlane(&VirtualNormalVector);
    1571               if (TempAngle < TempVector.Angle(&helper)) {
    1572                 TempAngle = NormalVector.Angle(&VirtualNormalVector);
    1573                 SmallestAngle = TempAngle;
    1574                 winner = target;
    1575                 *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
    1576               } else
    1577                 *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
    1578             } else
    1579               *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
    1580           }
    1581         } // end of loop over all boundary points
    1582 
    1583         // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
    1584         if (winner != PointsOnBoundary.end()) {
    1585           *out << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;
    1586           // create the lins of not yet present
    1587           BLS[0] = baseline->second;
    1588           // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
    1589           LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
    1590           LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
    1591           if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
    1592             BPS[0] = baseline->second->endpoints[0];
    1593             BPS[1] = winner->second;
    1594             BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    1595             LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1]));
    1596             LinesOnBoundaryCount++;
    1597           } else
    1598             BLS[1] = LineChecker[0]->second;
    1599           if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
    1600             BPS[0] = baseline->second->endpoints[1];
    1601             BPS[1] = winner->second;
    1602             BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    1603             LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2]));
    1604             LinesOnBoundaryCount++;
    1605           } else
    1606             BLS[2] = LineChecker[1]->second;
    1607           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    1608           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    1609           TrianglesOnBoundaryCount++;
    1610         } else {
    1611           *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
    1612         }
    1613 
    1614         // 5d. If the set of lines is not yet empty, go to 5. and continue
    1615       } else
    1616         *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
    1617   } while (flag);
    1618 
    1619   // exit
    1620   delete(MolCenter);
    1621 };
    1622 
    1623 /** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
    1624  * \param *out output stream for debugging
    1625  * \param *mol molecule structure with atoms
    1626  * \return true - all straddling points insert, false - something went wrong
    1627  */
    1628 bool Tesselation::InsertStraddlingPoints(ofstream *out, molecule *mol)
    1629 {
    1630   Vector Intersection;
    1631   atom *Walker = mol->start;
    1632   Vector *MolCenter = mol->DetermineCenterOfAll(out);
    1633   while (Walker->next != mol->end) {  // we only have to go once through all points, as boundary can become only bigger
    1634     // get the next triangle
    1635     BTS = FindClosestTriangleToPoint(out, &Walker->x);
    1636     if (BTS == NULL) {
    1637       *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
    1638       return false;
    1639     }
    1640     // get the intersection point
    1641     if (BTS->GetIntersectionInsideTriangle(out, MolCenter, &Walker->x, &Intersection)) {
    1642       // we have the intersection, check whether in- or outside of boundary
    1643       if ((MolCenter->DistanceSquared(&Walker->x) - MolCenter->DistanceSquared(&Intersection)) < -MYEPSILON) {
    1644         // inside, next!
    1645         *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
    1646       } else {
    1647         // outside!
    1648         *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
    1649         class BoundaryLineSet *OldLines[3], *NewLines[3];
    1650         class BoundaryPointSet *OldPoints[3], *NewPoint;
    1651         // store the three old lines and old points
    1652         for (int i=0;i<3;i++) {
    1653           OldLines[i] = BTS->lines[i];
    1654           OldPoints[i] = BTS->endpoints[i];
    1655         }
    1656         // add Walker to boundary points
    1657         AddPoint(Walker);
    1658         if (BPS[0] == NULL)
    1659           NewPoint = BPS[0];
    1660         else
    1661           continue;
    1662         // remove triangle
    1663         TrianglesOnBoundary.erase(BTS->Nr);
    1664         // create three new boundary lines
    1665         for (int i=0;i<3;i++) {
    1666           BPS[0] = NewPoint;
    1667           BPS[1] = OldPoints[i];
    1668           NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    1669           LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
    1670           LinesOnBoundaryCount++;
    1671         }
    1672         // create three new triangle with new point
    1673         for (int i=0;i<3;i++) { // find all baselines
    1674           BLS[0] = OldLines[i];
    1675           int n = 1;
    1676           for (int j=0;j<3;j++) {
    1677             if (NewLines[j]->IsConnectedTo(BLS[0])) {
    1678               if (n>2) {
    1679                 *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
    1680                 return false;
    1681               } else
    1682                 BLS[n++] = NewLines[j];
    1683             }
    1684           }
    1685           // create the triangle
    1686           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    1687           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    1688           TrianglesOnBoundaryCount++;
    1689         }
    1690       }
    1691     } else { // something is wrong with FindClosestTriangleToPoint!
    1692       *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
    1693       return false;
    1694     }
    1695   }
    1696 
    1697   // exit
    1698   delete(MolCenter);
    1699   return true;
    1700 };
    1701 
    1702 /** Goes over all baselines and checks whether adjacent triangles and convex to each other.
    1703  * \param *out output stream for debugging
    1704  * \return true - all baselines were corrected, false - there are still concave pieces
    1705  */
    1706 bool Tesselation::CorrectConcaveBaselines(ofstream *out)
    1707 {
    1708   class BoundaryTriangleSet *triangle[2];
    1709   class BoundaryLineSet *OldLines[4], *NewLine;
    1710   class BoundaryPointSet *OldPoints[2];
    1711   Vector BaseLineNormal;
    1712   class BoundaryLineSet *Base = NULL;
    1713   int OldTriangles[2], OldBaseLine;
    1714   int i;
    1715   for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
    1716     Base = baseline->second;
    1717 
    1718     // check convexity
    1719     if (Base->CheckConvexityCriterion(out)) { // triangles are convex
    1720       *out << Verbose(3) << Base << " has two convex triangles." << endl;
    1721     } else { // not convex!
    1722       // get the two triangles
    1723       i=0;
    1724       for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
    1725         triangle[i++] = runner->second;
    1726       // gather four endpoints and four lines
    1727       for (int j=0;j<4;j++)
    1728         OldLines[j] = NULL;
    1729       for (int j=0;j<2;j++)
    1730         OldPoints[j] = NULL;
    1731       i=0;
    1732       for (int m=0;m<2;m++) { // go over both triangles
    1733         for (int j=0;j<3;j++) { // all of their endpoints and baselines
    1734           if (triangle[m]->lines[j] != Base) // pick not the central baseline
    1735             OldLines[i++] = triangle[m]->lines[j];
    1736           if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
    1737             OldPoints[m] = triangle[m]->endpoints[j];
    1738         }
    1739       }
    1740       if (i<4) {
    1741         *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
    1742         return false;
    1743       }
    1744       for (int j=0;j<4;j++)
    1745         if (OldLines[j] == NULL) {
    1746           *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
    1747           return false;
    1748         }
    1749       for (int j=0;j<2;j++)
    1750         if (OldPoints[j] == NULL) {
    1751           *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
    1752           return false;
    1753         }
    1754 
    1755       // remove triangles
    1756       for (int j=0;j<2;j++) {
    1757         OldTriangles[j] = triangle[j]->Nr;
    1758         TrianglesOnBoundary.erase(OldTriangles[j]);
    1759         triangle[j] = NULL;
    1760       }
    1761 
    1762       // remove baseline
    1763       OldBaseLine = Base->Nr;
    1764       LinesOnBoundary.erase(OldBaseLine);
    1765       Base = NULL;
    1766 
    1767       // construct new baseline (with same number as old one)
    1768       BPS[0] = OldPoints[0];
    1769       BPS[1] = OldPoints[1];
    1770       NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
    1771       LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
    1772 
    1773       // construct new triangles with flipped baseline
    1774       i=-1;
    1775       if (BLS[0]->IsConnectedTo(OldLines[2]))
    1776         i=2;
    1777       if (BLS[0]->IsConnectedTo(OldLines[2]))
    1778         i=3;
    1779       if (i!=-1) {
    1780         BLS[0] = OldLines[0];
    1781         BLS[1] = OldLines[i];
    1782         BLS[2] = NewLine;
    1783         BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
    1784         TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
    1785 
    1786         BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
    1787         BLS[1] = OldLines[1];
    1788         BLS[2] = NewLine;
    1789         BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
    1790         TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
    1791       } else {
    1792         *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
    1793         return false;
    1794       }
    1795     }
    1796   }
    1797   return true;
    1798 };
    1799 
    1800 
    1801 /** States whether point is in- or outside of a tesselated surface.
    1802  * \param *pointer point to be checked
    1803  * \return true - is inside, false - is outside
    1804  */
    1805 bool Tesselation::IsInside(Vector *pointer)
    1806 {
    1807 
    1808   // hier kommt dann Saskias Routine hin...
    1809 
    1810   return true;
    1811 };
    1812 
    1813 /** Finds the closest triangle to a given point.
    1814  * \param *out output stream for debugging
    1815  * \param *x second endpoint of line
    1816  * \return pointer triangle that is closest, NULL if none was found
    1817  */
    1818 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
    1819 {
    1820   class BoundaryTriangleSet *triangle = NULL;
    1821 
    1822   // hier kommt dann Saskias Routine hin...
    1823 
    1824   return triangle;
    1825 };
    1826 
    1827 /** Adds an atom to the tesselation::PointsOnBoundary list.
    1828  * \param *Walker atom to add
    1829  */
    1830 void
    1831 Tesselation::AddPoint(atom *Walker)
    1832 {
    1833   PointTestPair InsertUnique;
    1834   BPS[0] = new class BoundaryPointSet(Walker);
    1835   InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
    1836   if (InsertUnique.second) // if new point was not present before, increase counter
    1837     PointsOnBoundaryCount++;
    1838   else {
    1839     delete(BPS[0]);
    1840     BPS[0] = NULL;
    1841   }
    1842 }
    1843 ;
    1844 
    1845 /** Adds point to Tesselation::PointsOnBoundary if not yet present.
    1846  * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
    1847  * @param Candidate point to add
    1848  * @param n index for this point in Tesselation::TPS array
    1849  */
    1850 void
    1851 Tesselation::AddTrianglePoint(atom* Candidate, int n)
    1852 {
    1853   PointTestPair InsertUnique;
    1854   TPS[n] = new class BoundaryPointSet(Candidate);
    1855   InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
    1856   if (InsertUnique.second) { // if new point was not present before, increase counter
    1857     PointsOnBoundaryCount++;
    1858   } else {
    1859     delete TPS[n];
    1860     cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
    1861     TPS[n] = (InsertUnique.first)->second;
    1862   }
    1863 }
    1864 ;
    1865 
    1866 /** Function tries to add line from current Points in BPS to BoundaryLineSet.
    1867  * If successful it raises the line count and inserts the new line into the BLS,
    1868  * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
    1869  * @param *a first endpoint
    1870  * @param *b second endpoint
    1871  * @param n index of Tesselation::BLS giving the line with both endpoints
    1872  */
    1873 void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
    1874   bool insertNewLine = true;
    1875 
    1876   if (a->lines.find(b->node->nr) != a->lines.end()) {
    1877     LineMap::iterator FindLine;
    1878     pair<LineMap::iterator,LineMap::iterator> FindPair;
    1879     FindPair = a->lines.equal_range(b->node->nr);
    1880 
    1881     for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
    1882       // If there is a line with less than two attached triangles, we don't need a new line.
    1883       if (FindLine->second->TrianglesCount < 2) {
    1884         insertNewLine = false;
    1885         cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
    1886 
    1887         BPS[0] = FindLine->second->endpoints[0];
    1888         BPS[1] = FindLine->second->endpoints[1];
    1889         BLS[n] = FindLine->second;
    1890 
    1891         break;
    1892       }
    1893     }
    1894   }
    1895 
    1896   if (insertNewLine) {
    1897     AlwaysAddTriangleLine(a, b, n);
    1898   }
    1899 }
    1900 ;
    1901 
    1902 /**
    1903  * Adds lines from each of the current points in the BPS to BoundaryLineSet.
    1904  * Raises the line count and inserts the new line into the BLS.
    1905  *
    1906  * @param *a first endpoint
    1907  * @param *b second endpoint
    1908  * @param n index of Tesselation::BLS giving the line with both endpoints
    1909  */
    1910 void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
    1911 {
    1912   cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
    1913   BPS[0] = a;
    1914   BPS[1] = b;
    1915   BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
    1916   // add line to global map
    1917   LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
    1918   // increase counter
    1919   LinesOnBoundaryCount++;
    1920 };
    1921 
    1922 /** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
    1923  * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
    1924  */
    1925 void
    1926 Tesselation::AddTriangle()
    1927 {
    1928   cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
    1929 
    1930   // add triangle to global map
    1931   TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    1932   TrianglesOnBoundaryCount++;
    1933 
    1934   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
    1935 }
    1936 ;
    1937 
    1938 
    1939 double det_get(gsl_matrix *A, int inPlace) {
    1940   /*
    1941   inPlace = 1 => A is replaced with the LU decomposed copy.
    1942   inPlace = 0 => A is retained, and a copy is used for LU.
    1943   */
    1944 
    1945   double det;
    1946   int signum;
    1947   gsl_permutation *p = gsl_permutation_alloc(A->size1);
    1948   gsl_matrix *tmpA;
    1949 
    1950   if (inPlace)
    1951   tmpA = A;
    1952   else {
    1953   gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
    1954   gsl_matrix_memcpy(tmpA , A);
    1955   }
    1956 
    1957 
    1958   gsl_linalg_LU_decomp(tmpA , p , &signum);
    1959   det = gsl_linalg_LU_det(tmpA , signum);
    1960   gsl_permutation_free(p);
    1961   if (! inPlace)
    1962   gsl_matrix_free(tmpA);
    1963 
    1964   return det;
    1965 };
    1966 
    1967 void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
    1968 {
    1969   gsl_matrix *A = gsl_matrix_calloc(3,3);
    1970   double m11, m12, m13, m14;
    1971 
    1972   for(int i=0;i<3;i++) {
    1973     gsl_matrix_set(A, i, 0, a.x[i]);
    1974     gsl_matrix_set(A, i, 1, b.x[i]);
    1975     gsl_matrix_set(A, i, 2, c.x[i]);
    1976   }
    1977   m11 = det_get(A, 1);
    1978 
    1979   for(int i=0;i<3;i++) {
    1980     gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    1981     gsl_matrix_set(A, i, 1, b.x[i]);
    1982     gsl_matrix_set(A, i, 2, c.x[i]);
    1983   }
    1984   m12 = det_get(A, 1);
    1985 
    1986   for(int i=0;i<3;i++) {
    1987     gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    1988     gsl_matrix_set(A, i, 1, a.x[i]);
    1989     gsl_matrix_set(A, i, 2, c.x[i]);
    1990   }
    1991   m13 = det_get(A, 1);
    1992 
    1993   for(int i=0;i<3;i++) {
    1994     gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    1995     gsl_matrix_set(A, i, 1, a.x[i]);
    1996     gsl_matrix_set(A, i, 2, b.x[i]);
    1997   }
    1998   m14 = det_get(A, 1);
    1999 
    2000   if (fabs(m11) < MYEPSILON)
    2001     cerr << "ERROR: three points are colinear." << endl;
    2002 
    2003   center->x[0] =  0.5 * m12/ m11;
    2004   center->x[1] = -0.5 * m13/ m11;
    2005   center->x[2] =  0.5 * m14/ m11;
    2006 
    2007   if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
    2008     cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
    2009 
    2010   gsl_matrix_free(A);
    2011 };
    2012 
    2013 
    2014 
    2015 /**
    2016  * Function returns center of sphere with RADIUS, which rests on points a, b, c
    2017  * @param Center this vector will be used for return
    2018  * @param a vector first point of triangle
    2019  * @param b vector second point of triangle
    2020  * @param c vector third point of triangle
    2021  * @param *Umkreismittelpunkt new cneter point of circumference
    2022  * @param Direction vector indicates up/down
    2023  * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
    2024  * @param Halfplaneindicator double indicates whether Direction is up or down
    2025  * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
    2026  * @param alpha double angle at a
    2027  * @param beta double, angle at b
    2028  * @param gamma, double, angle at c
    2029  * @param Radius, double
    2030  * @param Umkreisradius double radius of circumscribing circle
    2031  */
    2032 void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
    2033     double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
    2034 {
    2035   Vector TempNormal, helper;
    2036   double Restradius;
    2037   Vector OtherCenter;
    2038   cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
    2039   Center->Zero();
    2040   helper.CopyVector(&a);
    2041   helper.Scale(sin(2.*alpha));
    2042   Center->AddVector(&helper);
    2043   helper.CopyVector(&b);
    2044   helper.Scale(sin(2.*beta));
    2045   Center->AddVector(&helper);
    2046   helper.CopyVector(&c);
    2047   helper.Scale(sin(2.*gamma));
    2048   Center->AddVector(&helper);
    2049   //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
    2050   Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    2051   NewUmkreismittelpunkt->CopyVector(Center);
    2052   cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
    2053   // Here we calculated center of circumscribing circle, using barycentric coordinates
    2054   cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
    2055 
    2056   TempNormal.CopyVector(&a);
    2057   TempNormal.SubtractVector(&b);
    2058   helper.CopyVector(&a);
    2059   helper.SubtractVector(&c);
    2060   TempNormal.VectorProduct(&helper);
    2061   if (fabs(HalfplaneIndicator) < MYEPSILON)
    2062     {
    2063       if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
    2064         {
    2065           TempNormal.Scale(-1);
    2066         }
    2067     }
    2068   else
    2069     {
    2070       if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
    2071         {
    2072           TempNormal.Scale(-1);
    2073         }
    2074     }
    2075 
    2076   TempNormal.Normalize();
    2077   Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
    2078   cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
    2079   TempNormal.Scale(Restradius);
    2080   cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
    2081 
    2082   Center->AddVector(&TempNormal);
    2083   cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
    2084   get_sphere(&OtherCenter, a, b, c, RADIUS);
    2085   cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
    2086   cout << Verbose(3) << "End of Get_center_of_sphere.\n";
    2087 };
    2088 
    2089 
    2090 /** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
    2091  * \param *Center new center on return
    2092  * \param *a first point
    2093  * \param *b second point
    2094  * \param *c third point
    2095  */
    2096 void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
    2097 {
    2098   Vector helper;
    2099   double alpha, beta, gamma;
    2100   Vector SideA, SideB, SideC;
    2101   SideA.CopyVector(b);
    2102   SideA.SubtractVector(c);
    2103   SideB.CopyVector(c);
    2104   SideB.SubtractVector(a);
    2105   SideC.CopyVector(a);
    2106   SideC.SubtractVector(b);
    2107   alpha = M_PI - SideB.Angle(&SideC);
    2108   beta = M_PI - SideC.Angle(&SideA);
    2109   gamma = M_PI - SideA.Angle(&SideB);
    2110   //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
    2111   if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
    2112     cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
    2113 
    2114   Center->Zero();
    2115   helper.CopyVector(a);
    2116   helper.Scale(sin(2.*alpha));
    2117   Center->AddVector(&helper);
    2118   helper.CopyVector(b);
    2119   helper.Scale(sin(2.*beta));
    2120   Center->AddVector(&helper);
    2121   helper.CopyVector(c);
    2122   helper.Scale(sin(2.*gamma));
    2123   Center->AddVector(&helper);
    2124   Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    2125 };
    2126 
    2127 /** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
    2128  * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
    2129  * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
    2130  * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
    2131  * \param CircleCenter Center of the parameter circle
    2132  * \param CirclePlaneNormal normal vector to plane of the parameter circle
    2133  * \param CircleRadius radius of the parameter circle
    2134  * \param NewSphereCenter new center of a circumcircle
    2135  * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
    2136  * \param NormalVector normal vector
    2137  * \param SearchDirection search direction to make angle unique on return.
    2138  * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
    2139  */
    2140 double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
    2141 {
    2142   Vector helper;
    2143   double radius, alpha;
    2144 
    2145   helper.CopyVector(&NewSphereCenter);
    2146   // test whether new center is on the parameter circle's plane
    2147   if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    2148     cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
    2149     helper.ProjectOntoPlane(&CirclePlaneNormal);
    2150   }
    2151   radius = helper.ScalarProduct(&helper);
    2152   // test whether the new center vector has length of CircleRadius
    2153   if (fabs(radius - CircleRadius) > HULLEPSILON)
    2154     cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    2155   alpha = helper.Angle(&OldSphereCenter);
    2156   // make the angle unique by checking the halfplanes/search direction
    2157   if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    2158     alpha = 2.*M_PI - alpha;
    2159   //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
    2160   radius = helper.Distance(&OldSphereCenter);
    2161   helper.ProjectOntoPlane(&NormalVector);
    2162   // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    2163   if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    2164     //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
    2165     return alpha;
    2166   } else {
    2167     //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
    2168     return 2.*M_PI;
    2169   }
    2170 };
    2171 
    2172 
    2173 /** Checks whether the triangle consisting of the three atoms is already present.
    2174  * Searches for the points in Tesselation::PointsOnBoundary and checks their
    2175  * lines. If any of the three edges already has two triangles attached, false is
    2176  * returned.
    2177  * \param *out output stream for debugging
    2178  * \param *Candidates endpoints of the triangle candidate
    2179  * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
    2180  *                 triangles exist which is the maximum for three points
    2181  */
    2182 int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
    2183 // TODO: use findTriangles here and return result.size();
    2184   LineMap::iterator FindLine;
    2185   PointMap::iterator FindPoint;
    2186   TriangleMap::iterator FindTriangle;
    2187   int adjacentTriangleCount = 0;
    2188   class BoundaryPointSet *Points[3];
    2189 
    2190   //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
    2191   // builds a triangle point set (Points) of the end points
    2192   for (int i = 0; i < 3; i++) {
    2193     FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
    2194     if (FindPoint != PointsOnBoundary.end()) {
    2195       Points[i] = FindPoint->second;
    2196     } else {
    2197       Points[i] = NULL;
    2198     }
    2199   }
    2200 
    2201   // checks lines between the points in the Points for their adjacent triangles
    2202   for (int i = 0; i < 3; i++) {
    2203     if (Points[i] != NULL) {
    2204       for (int j = i; j < 3; j++) {
    2205         if (Points[j] != NULL) {
    2206           FindLine = Points[i]->lines.find(Points[j]->node->nr);
    2207           if (FindLine != Points[i]->lines.end()) {
    2208             for (; FindLine->first == Points[j]->node->nr; FindLine++) {
    2209               FindTriangle = FindLine->second->triangles.begin();
    2210               for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
    2211                 if ((
    2212                   (FindTriangle->second->endpoints[0] == Points[0])
    2213                     || (FindTriangle->second->endpoints[0] == Points[1])
    2214                     || (FindTriangle->second->endpoints[0] == Points[2])
    2215                   ) && (
    2216                     (FindTriangle->second->endpoints[1] == Points[0])
    2217                     || (FindTriangle->second->endpoints[1] == Points[1])
    2218                     || (FindTriangle->second->endpoints[1] == Points[2])
    2219                   ) && (
    2220                     (FindTriangle->second->endpoints[2] == Points[0])
    2221                     || (FindTriangle->second->endpoints[2] == Points[1])
    2222                     || (FindTriangle->second->endpoints[2] == Points[2])
    2223                   )
    2224                 ) {
    2225                   adjacentTriangleCount++;
    2226                 }
    2227               }
    2228             }
    2229             // Only one of the triangle lines must be considered for the triangle count.
    2230             *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
    2231             return adjacentTriangleCount;
    2232 
    2233           }
    2234         }
    2235       }
    2236     }
    2237   }
    2238 
    2239   *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
    2240   return adjacentTriangleCount;
    2241 };
    2242 
    2243 /** This recursive function finds a third point, to form a triangle with two given ones.
    2244  * Note that this function is for the starting triangle.
    2245  * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
    2246  * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
    2247  * the center of the sphere is still fixed up to a single parameter. The band of possible values
    2248  * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
    2249  * us the "null" on this circle, the new center of the candidate point will be some way along this
    2250  * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
    2251  * by the normal vector of the base triangle that always points outwards by construction.
    2252  * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
    2253  * We construct the normal vector that defines the plane this circle lies in, it is just in the
    2254  * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
    2255  * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
    2256  * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
    2257  * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
    2258  * both.
    2259  * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
    2260  * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
    2261  * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
    2262  * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
    2263  * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
    2264  * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
    2265  * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
    2266  * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
    2267  * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
    2268  * @param BaseLine BoundaryLineSet with the current base line
    2269  * @param ThirdNode third atom to avoid in search
    2270  * @param candidates list of equally good candidates to return
    2271  * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
    2272  * @param RADIUS radius of sphere
    2273  * @param *LC LinkedCell structure with neighbouring atoms
    2274  */
    2275 void Find_third_point_for_Tesselation(
    2276   Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter,
    2277   class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates,
    2278   double *ShortestAngle, const double RADIUS, LinkedCell *LC
    2279 ) {
    2280   Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
    2281   Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
    2282   Vector SphereCenter;
    2283   Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
    2284   Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
    2285   Vector NewNormalVector;   // normal vector of the Candidate's triangle
    2286   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
    2287   LinkedAtoms *List = NULL;
    2288   double CircleRadius; // radius of this circle
    2289   double radius;
    2290   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
    2291   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
    2292   atom *Candidate = NULL;
    2293   CandidateForTesselation *optCandidate = NULL;
    2294 
    2295   cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
    2296 
    2297   //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
    2298 
    2299   // construct center of circle
    2300   CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
    2301   CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
    2302   CircleCenter.Scale(0.5);
    2303 
    2304   // construct normal vector of circle
    2305   CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
    2306   CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
    2307 
    2308   // calculate squared radius atom *ThirdNode,f circle
    2309   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
    2310   if (radius/4. < RADIUS*RADIUS) {
    2311     CircleRadius = RADIUS*RADIUS - radius/4.;
    2312     CirclePlaneNormal.Normalize();
    2313     //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
    2314 
    2315     // test whether old center is on the band's plane
    2316     if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    2317       cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
    2318       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
    2319     }
    2320     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
    2321     if (fabs(radius - CircleRadius) < HULLEPSILON) {
    2322 
    2323       // check SearchDirection
    2324       //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
    2325       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
    2326         cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
    2327       }
    2328 
    2329       // get cell for the starting atom
    2330       if (LC->SetIndexToVector(&CircleCenter)) {
    2331         for(int i=0;i<NDIM;i++) // store indices of this cell
    2332         N[i] = LC->n[i];
    2333         //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    2334       } else {
    2335         cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
    2336         return;
    2337       }
    2338       // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
    2339       //cout << Verbose(2) << "LC Intervals:";
    2340       for (int i=0;i<NDIM;i++) {
    2341         Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
    2342         Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
    2343         //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
    2344       }
    2345       //cout << endl;
    2346       for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    2347         for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    2348           for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    2349             List = LC->GetCurrentCell();
    2350             //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
    2351             if (List != NULL) {
    2352               for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    2353                 Candidate = (*Runner);
    2354 
    2355                 // check for three unique points
    2356                 //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
    2357                 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
    2358 
    2359                   // construct both new centers
    2360                   GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
    2361                   OtherNewSphereCenter.CopyVector(&NewSphereCenter);
    2362 
    2363                   if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
    2364                   && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
    2365                   ) {
    2366                     helper.CopyVector(&NewNormalVector);
    2367                     //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
    2368                     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
    2369                     if (radius < RADIUS*RADIUS) {
    2370                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
    2371                       //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
    2372                       NewSphereCenter.AddVector(&helper);
    2373                       NewSphereCenter.SubtractVector(&CircleCenter);
    2374                       //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
    2375 
    2376                       // OtherNewSphereCenter is created by the same vector just in the other direction
    2377                       helper.Scale(-1.);
    2378                       OtherNewSphereCenter.AddVector(&helper);
    2379                       OtherNewSphereCenter.SubtractVector(&CircleCenter);
    2380                       //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
    2381 
    2382                       alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
    2383                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
    2384                       alpha = min(alpha, Otheralpha);
    2385                       // if there is a better candidate, drop the current list and add the new candidate
    2386                       // otherwise ignore the new candidate and keep the list
    2387                       if (*ShortestAngle > (alpha - HULLEPSILON)) {
    2388                         optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
    2389                         if (fabs(alpha - Otheralpha) > MYEPSILON) {
    2390                           optCandidate->OptCenter.CopyVector(&NewSphereCenter);
    2391                           optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
    2392                         } else {
    2393                           optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
    2394                           optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
    2395                         }
    2396                         // if there is an equal candidate, add it to the list without clearing the list
    2397                         if ((*ShortestAngle - HULLEPSILON) < alpha) {
    2398                           candidates->push_back(optCandidate);
    2399                           cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
    2400                             << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
    2401                         } else {
    2402                           // remove all candidates from the list and then the list itself
    2403                           class CandidateForTesselation *remover = NULL;
    2404                           for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
    2405                             remover = *it;
    2406                             delete(remover);
    2407                           }
    2408                           candidates->clear();
    2409                           candidates->push_back(optCandidate);
    2410                           cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
    2411                             << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
    2412                         }
    2413                         *ShortestAngle = alpha;
    2414                         //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
    2415                       } else {
    2416                         if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
    2417                           //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
    2418                         } else {
    2419                           //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
    2420                         }
    2421                       }
    2422 
    2423                     } else {
    2424                       //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
    2425                     }
    2426                   } else {
    2427                     //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
    2428                   }
    2429                 } else {
    2430                   if (ThirdNode != NULL) {
    2431                     //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
    2432                   } else {
    2433                     //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
    2434                   }
    2435                 }
    2436               }
    2437             }
    2438           }
    2439     } else {
    2440       cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    2441     }
    2442   } else {
    2443     if (ThirdNode != NULL)
    2444       cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
    2445     else
    2446       cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
    2447   }
    2448 
    2449   //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
    2450   if (candidates->size() > 1) {
    2451     candidates->unique();
    2452     candidates->sort(sortCandidates);
    2453   }
    2454 
    2455   cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
    2456 };
    2457 
    2458 struct Intersection {
    2459   Vector x1;
    2460   Vector x2;
    2461   Vector x3;
    2462   Vector x4;
    2463 };
    2464 
    2465 /**
    2466  * Intersection calculation function.
    2467  *
    2468  * @param x to find the result for
    2469  * @param function parameter
    2470  */
    2471 double MinIntersectDistance(const gsl_vector * x, void *params) {
    2472   double retval = 0;
    2473   struct Intersection *I = (struct Intersection *)params;
    2474   Vector intersection;
    2475   Vector SideA,SideB,HeightA, HeightB;
    2476   for (int i=0;i<NDIM;i++)
    2477     intersection.x[i] = gsl_vector_get(x, i);
    2478 
    2479   SideA.CopyVector(&(I->x1));
    2480   SideA.SubtractVector(&I->x2);
    2481   HeightA.CopyVector(&intersection);
    2482   HeightA.SubtractVector(&I->x1);
    2483   HeightA.ProjectOntoPlane(&SideA);
    2484 
    2485   SideB.CopyVector(&I->x3);
    2486   SideB.SubtractVector(&I->x4);
    2487   HeightB.CopyVector(&intersection);
    2488   HeightB.SubtractVector(&I->x3);
    2489   HeightB.ProjectOntoPlane(&SideB);
    2490 
    2491   retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
    2492   //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
    2493 
    2494   return retval;
    2495 };
    2496 
    2497 
    2498 /**
    2499  * Calculates whether there is an intersection between two lines. The first line
    2500  * always goes through point 1 and point 2 and the second line is given by the
    2501  * connection between point 4 and point 5.
    2502  *
    2503  * @param point 1 of line 1
    2504  * @param point 2 of line 1
    2505  * @param point 1 of line 2
    2506  * @param point 2 of line 2
    2507  *
    2508  * @return true if there is an intersection between the given lines, false otherwise
    2509  */
    2510 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
    2511   bool result;
    2512 
    2513   struct Intersection par;
    2514     par.x1.CopyVector(&point1);
    2515     par.x2.CopyVector(&point2);
    2516     par.x3.CopyVector(&point3);
    2517     par.x4.CopyVector(&point4);
    2518 
    2519     const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
    2520     gsl_multimin_fminimizer *s = NULL;
    2521     gsl_vector *ss, *x;
    2522     gsl_multimin_function minex_func;
    2523 
    2524     size_t iter = 0;
    2525     int status;
    2526     double size;
    2527 
    2528     /* Starting point */
    2529     x = gsl_vector_alloc(NDIM);
    2530     gsl_vector_set(x, 0, point1.x[0]);
    2531     gsl_vector_set(x, 1, point1.x[1]);
    2532     gsl_vector_set(x, 2, point1.x[2]);
    2533 
    2534     /* Set initial step sizes to 1 */
    2535     ss = gsl_vector_alloc(NDIM);
    2536     gsl_vector_set_all(ss, 1.0);
    2537 
    2538     /* Initialize method and iterate */
    2539     minex_func.n = NDIM;
    2540     minex_func.f = &MinIntersectDistance;
    2541     minex_func.params = (void *)&par;
    2542 
    2543     s = gsl_multimin_fminimizer_alloc(T, NDIM);
    2544     gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
    2545 
    2546     do {
    2547         iter++;
    2548         status = gsl_multimin_fminimizer_iterate(s);
    2549 
    2550         if (status) {
    2551           break;
    2552         }
    2553 
    2554         size = gsl_multimin_fminimizer_size(s);
    2555         status = gsl_multimin_test_size(size, 1e-2);
    2556 
    2557         if (status == GSL_SUCCESS) {
    2558           cout << Verbose(2) << "converged to minimum" <<  endl;
    2559         }
    2560     } while (status == GSL_CONTINUE && iter < 100);
    2561 
    2562     // check whether intersection is in between or not
    2563   Vector intersection, SideA, SideB, HeightA, HeightB;
    2564   double t1, t2;
    2565   for (int i = 0; i < NDIM; i++) {
    2566     intersection.x[i] = gsl_vector_get(s->x, i);
    2567   }
    2568 
    2569   SideA.CopyVector(&par.x2);
    2570   SideA.SubtractVector(&par.x1);
    2571   HeightA.CopyVector(&intersection);
    2572   HeightA.SubtractVector(&par.x1);
    2573 
    2574   t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
    2575 
    2576   SideB.CopyVector(&par.x4);
    2577   SideB.SubtractVector(&par.x3);
    2578   HeightB.CopyVector(&intersection);
    2579   HeightB.SubtractVector(&par.x3);
    2580 
    2581   t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
    2582 
    2583   cout << Verbose(2) << "Intersection " << intersection << " is at "
    2584     << t1 << " for (" << point1 << "," << point2 << ") and at "
    2585     << t2 << " for (" << point3 << "," << point4 << "): ";
    2586 
    2587   if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    2588     cout << "true intersection." << endl;
    2589     result = true;
    2590   } else {
    2591     cout << "intersection out of region of interest." << endl;
    2592     result = false;
    2593   }
    2594 
    2595   // free minimizer stuff
    2596     gsl_vector_free(x);
    2597     gsl_vector_free(ss);
    2598     gsl_multimin_fminimizer_free(s);
    2599 
    2600   return result;
    2601 }
    2602 
    2603 /** Finds the second point of starting triangle.
    2604  * \param *a first atom
    2605  * \param *Candidate pointer to candidate atom on return
    2606  * \param Oben vector indicating the outside
    2607  * \param Opt_Candidate reference to recommended candidate on return
    2608  * \param Storage[3] array storing angles and other candidate information
    2609  * \param RADIUS radius of virtual sphere
    2610  * \param *LC LinkedCell structure with neighbouring atoms
    2611  */
    2612 void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
    2613 {
    2614   cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
    2615   Vector AngleCheck;
    2616   double norm = -1., angle;
    2617   LinkedAtoms *List = NULL;
    2618   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
    2619 
    2620   if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
    2621     for(int i=0;i<NDIM;i++) // store indices of this cell
    2622       N[i] = LC->n[i];
    2623   } else {
    2624     cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
    2625     return;
    2626   }
    2627   // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
    2628   cout << Verbose(3) << "LC Intervals from [";
    2629   for (int i=0;i<NDIM;i++) {
    2630   cout << " " << N[i] << "<->" << LC->N[i];
    2631   }
    2632   cout << "] :";
    2633   for (int i=0;i<NDIM;i++) {
    2634     Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
    2635     Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
    2636     cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
    2637   }
    2638   cout << endl;
    2639 
    2640 
    2641   for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    2642     for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    2643       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    2644         List = LC->GetCurrentCell();
    2645         //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
    2646         if (List != NULL) {
    2647           for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    2648             Candidate = (*Runner);
    2649             // check if we only have one unique point yet ...
    2650             if (a != Candidate) {
    2651               // Calculate center of the circle with radius RADIUS through points a and Candidate
    2652               Vector OrthogonalizedOben, a_Candidate, Center;
    2653               double distance, scaleFactor;
    2654 
    2655               OrthogonalizedOben.CopyVector(&Oben);
    2656               a_Candidate.CopyVector(&(a->x));
    2657               a_Candidate.SubtractVector(&(Candidate->x));
    2658               OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
    2659               OrthogonalizedOben.Normalize();
    2660               distance = 0.5 * a_Candidate.Norm();
    2661               scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
    2662               OrthogonalizedOben.Scale(scaleFactor);
    2663 
    2664               Center.CopyVector(&(Candidate->x));
    2665               Center.AddVector(&(a->x));
    2666               Center.Scale(0.5);
    2667               Center.AddVector(&OrthogonalizedOben);
    2668 
    2669               AngleCheck.CopyVector(&Center);
    2670               AngleCheck.SubtractVector(&(a->x));
    2671               norm = a_Candidate.Norm();
    2672               // second point shall have smallest angle with respect to Oben vector
    2673               if (norm < RADIUS*2.) {
    2674                 angle = AngleCheck.Angle(&Oben);
    2675                 if (angle < Storage[0]) {
    2676                   //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
    2677                   cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
    2678                   Opt_Candidate = Candidate;
    2679                   Storage[0] = angle;
    2680                   //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
    2681                 } else {
    2682                   //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
    2683                 }
    2684               } else {
    2685                 //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
    2686               }
    2687             } else {
    2688               //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
    2689             }
    2690           }
    2691         } else {
    2692           cout << Verbose(3) << "Linked cell list is empty." << endl;
    2693         }
    2694       }
    2695   cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
    2696 };
    2697 
    2698 /** Finds the starting triangle for find_non_convex_border().
    2699  * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
    2700  * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
    2701  * point are called.
    2702  * \param RADIUS radius of virtual rolling sphere
    2703  * \param *LC LinkedCell structure with neighbouring atoms
    2704  */
    2705 void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
    2706 {
    2707   cout << Verbose(1) << "Begin of Find_starting_triangle\n";
    2708   int i = 0;
    2709   LinkedAtoms *List = NULL;
    2710   atom* FirstPoint = NULL;
    2711   atom* SecondPoint = NULL;
    2712   atom* MaxAtom[NDIM];
    2713   double max_coordinate[NDIM];
    2714   Vector Oben;
    2715   Vector helper;
    2716   Vector Chord;
    2717   Vector SearchDirection;
    2718 
    2719   Oben.Zero();
    2720 
    2721   for (i = 0; i < 3; i++) {
    2722     MaxAtom[i] = NULL;
    2723     max_coordinate[i] = -1;
    2724   }
    2725 
    2726   // 1. searching topmost atom with respect to each axis
    2727   for (int i=0;i<NDIM;i++) { // each axis
    2728     LC->n[i] = LC->N[i]-1; // current axis is topmost cell
    2729     for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
    2730       for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
    2731         List = LC->GetCurrentCell();
    2732         //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
    2733         if (List != NULL) {
    2734           for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
    2735             if ((*Runner)->x.x[i] > max_coordinate[i]) {
    2736               cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl;
    2737               max_coordinate[i] = (*Runner)->x.x[i];
    2738               MaxAtom[i] = (*Runner);
    2739             }
    2740           }
    2741         } else {
    2742           cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
    2743         }
    2744       }
    2745   }
    2746 
    2747   cout << Verbose(2) << "Found maximum coordinates: ";
    2748   for (int i=0;i<NDIM;i++)
    2749     cout << i << ": " << *MaxAtom[i] << "\t";
    2750   cout << endl;
    2751 
    2752   BTS = NULL;
    2753   CandidateList *Opt_Candidates = new CandidateList();
    2754   for (int k=0;k<NDIM;k++) {
    2755     Oben.x[k] = 1.;
    2756     FirstPoint = MaxAtom[k];
    2757     cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
    2758 
    2759     double ShortestAngle;
    2760     atom* Opt_Candidate = NULL;
    2761     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.
    2762 
    2763     Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
    2764     SecondPoint = Opt_Candidate;
    2765     if (SecondPoint == NULL)  // have we found a second point?
    2766       continue;
    2767     else
    2768       cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
    2769 
    2770     helper.CopyVector(&(FirstPoint->x));
    2771     helper.SubtractVector(&(SecondPoint->x));
    2772     helper.Normalize();
    2773     Oben.ProjectOntoPlane(&helper);
    2774     Oben.Normalize();
    2775     helper.VectorProduct(&Oben);
    2776     ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
    2777 
    2778     Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
    2779     Chord.SubtractVector(&(SecondPoint->x));
    2780     double radius = Chord.ScalarProduct(&Chord);
    2781     double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
    2782     helper.CopyVector(&Oben);
    2783     helper.Scale(CircleRadius);
    2784     // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
    2785 
    2786     // look in one direction of baseline for initial candidate
    2787     SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
    2788 
    2789     // adding point 1 and point 2 and the line between them
    2790     AddTrianglePoint(FirstPoint, 0);
    2791     AddTrianglePoint(SecondPoint, 1);
    2792     AddTriangleLine(TPS[0], TPS[1], 0);
    2793 
    2794     //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
    2795     Find_third_point_for_Tesselation(
    2796       Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
    2797     );
    2798     cout << Verbose(1) << "List of third Points is ";
    2799     for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2800         cout << " " << *(*it)->point;
    2801     }
    2802     cout << endl;
    2803 
    2804     for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2805       // add third triangle point
    2806       AddTrianglePoint((*it)->point, 2);
    2807       // add the second and third line
    2808       AddTriangleLine(TPS[1], TPS[2], 1);
    2809       AddTriangleLine(TPS[0], TPS[2], 2);
    2810       // ... and triangles to the Maps of the Tesselation class
    2811       BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    2812       AddTriangle();
    2813       // ... and calculate its normal vector (with correct orientation)
    2814       (*it)->OptCenter.Scale(-1.);
    2815       cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
    2816       BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
    2817       cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
    2818       << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
    2819 
    2820       // if we do not reach the end with the next step of iteration, we need to setup a new first line
    2821       if (it != Opt_Candidates->end()--) {
    2822         FirstPoint = (*it)->BaseLine->endpoints[0]->node;
    2823         SecondPoint = (*it)->point;
    2824         // adding point 1 and point 2 and the line between them
    2825         AddTrianglePoint(FirstPoint, 0);
    2826         AddTrianglePoint(SecondPoint, 1);
    2827         AddTriangleLine(TPS[0], TPS[1], 0);
    2828       }
    2829       cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
    2830     }
    2831     if (BTS != NULL) // we have created one starting triangle
    2832       break;
    2833     else {
    2834       // remove all candidates from the list and then the list itself
    2835       class CandidateForTesselation *remover = NULL;
    2836       for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2837         remover = *it;
    2838         delete(remover);
    2839       }
    2840       Opt_Candidates->clear();
    2841     }
    2842   }
    2843 
    2844   // remove all candidates from the list and then the list itself
    2845   class CandidateForTesselation *remover = NULL;
    2846   for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2847     remover = *it;
    2848     delete(remover);
    2849   }
    2850   delete(Opt_Candidates);
    2851   cout << Verbose(1) << "End of Find_starting_triangle\n";
    2852 };
    2853 
    2854 /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
    2855  * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
    2856  * make it bigger (i.e. closing one (the baseline) and opening two new ones).
    2857  * \param TPS[3] nodes of the triangle
    2858  * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
    2859  */
    2860 bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
    2861 {
    2862   bool result = false;
    2863   int counter = 0;
    2864 
    2865   // check all three points
    2866   for (int i=0;i<3;i++)
    2867     for (int j=i+1; j<3; j++) {
    2868       if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
    2869         LineMap::iterator FindLine;
    2870         pair<LineMap::iterator,LineMap::iterator> FindPair;
    2871         FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
    2872         for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
    2873           // If there is a line with less than two attached triangles, we don't need a new line.
    2874           if (FindLine->second->TrianglesCount < 2) {
    2875             counter++;
    2876             break;  // increase counter only once per edge
    2877           }
    2878         }
    2879       } else { // no line
    2880         cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
    2881         result = true;
    2882       }
    2883     }
    2884   if (counter > 1) {
    2885     cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    2886     result = true;
    2887   }
    2888   return result;
    2889 };
    2890 
    2891 
    2892 /** This function finds a triangle to a line, adjacent to an existing one.
    2893  * @param out output stream for debugging
    2894  * @param *mol molecule with Atom's and Bond's
    2895  * @param Line current baseline to search from
    2896  * @param T current triangle which \a Line is edge of
    2897  * @param RADIUS radius of the rolling ball
    2898  * @param N number of found triangles
    2899  * @param *filename filename base for intermediate envelopes
    2900  * @param *LC LinkedCell structure with neighbouring atoms
    2901  */
    2902 bool Tesselation::Find_next_suitable_triangle(ofstream *out,
    2903     molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
    2904     const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
    2905 {
    2906   cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
    2907   ofstream *tempstream = NULL;
    2908   char NumberName[255];
    2909   bool result = true;
    2910   CandidateList *Opt_Candidates = new CandidateList();
    2911 
    2912   Vector CircleCenter;
    2913   Vector CirclePlaneNormal;
    2914   Vector OldSphereCenter;
    2915   Vector SearchDirection;
    2916   Vector helper;
    2917   atom *ThirdNode = NULL;
    2918   LineMap::iterator testline;
    2919   double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
    2920   double radius, CircleRadius;
    2921 
    2922   cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
    2923   for (int i=0;i<3;i++)
    2924     if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
    2925       ThirdNode = T.endpoints[i]->node;
    2926 
    2927   // construct center of circle
    2928   CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
    2929   CircleCenter.AddVector(&Line.endpoints[1]->node->x);
    2930   CircleCenter.Scale(0.5);
    2931 
    2932   // construct normal vector of circle
    2933   CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
    2934   CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
    2935 
    2936   // calculate squared radius of circle
    2937   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
    2938   if (radius/4. < RADIUS*RADIUS) {
    2939     CircleRadius = RADIUS*RADIUS - radius/4.;
    2940     CirclePlaneNormal.Normalize();
    2941     cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
    2942 
    2943     // construct old center
    2944     GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
    2945     helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
    2946     radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
    2947     helper.Scale(sqrt(RADIUS*RADIUS - radius));
    2948     OldSphereCenter.AddVector(&helper);
    2949     OldSphereCenter.SubtractVector(&CircleCenter);
    2950     //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
    2951 
    2952     // construct SearchDirection
    2953     SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
    2954     helper.CopyVector(&Line.endpoints[0]->node->x);
    2955     helper.SubtractVector(&ThirdNode->x);
    2956     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
    2957       SearchDirection.Scale(-1.);
    2958     SearchDirection.ProjectOntoPlane(&OldSphereCenter);
    2959     SearchDirection.Normalize();
    2960     cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
    2961     if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
    2962       // rotated the wrong way!
    2963       cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
    2964     }
    2965 
    2966     // add third point
    2967     Find_third_point_for_Tesselation(
    2968       T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
    2969       &ShortestAngle, RADIUS, LC
    2970     );
    2971 
    2972   } else {
    2973     cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
    2974   }
    2975 
    2976   if (Opt_Candidates->begin() == Opt_Candidates->end()) {
    2977     cerr << "WARNING: Could not find a suitable candidate." << endl;
    2978     return false;
    2979   }
    2980   cout << Verbose(1) << "Third Points are ";
    2981   for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2982     cout << " " << *(*it)->point;
    2983   }
    2984   cout << endl;
    2985 
    2986   BoundaryLineSet *BaseRay = &Line;
    2987   for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2988     cout << Verbose(1) << " Third point candidate is " << *(*it)->point
    2989     << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
    2990     cout << Verbose(1) << " Baseline is " << *BaseRay << endl;
    2991 
    2992     // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
    2993     atom *AtomCandidates[3];
    2994     AtomCandidates[0] = (*it)->point;
    2995     AtomCandidates[1] = BaseRay->endpoints[0]->node;
    2996     AtomCandidates[2] = BaseRay->endpoints[1]->node;
    2997     int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
    2998 
    2999     BTS = NULL;
    3000     // If there is no triangle, add it regularly.
    3001     if (existentTrianglesCount == 0) {
    3002       AddTrianglePoint((*it)->point, 0);
    3003       AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
    3004       AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
    3005 
    3006       AddTriangleLine(TPS[0], TPS[1], 0);
    3007       AddTriangleLine(TPS[0], TPS[2], 1);
    3008       AddTriangleLine(TPS[1], TPS[2], 2);
    3009 
    3010       BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    3011       AddTriangle();
    3012       (*it)->OptCenter.Scale(-1.);
    3013       BTS->GetNormalVector((*it)->OptCenter);
    3014       (*it)->OptCenter.Scale(-1.);
    3015 
    3016       cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
    3017         << " for this triangle ... " << endl;
    3018       //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
    3019     } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
    3020         AddTrianglePoint((*it)->point, 0);
    3021         AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
    3022         AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
    3023 
    3024         // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
    3025         // i.e. at least one of the three lines must be present with TriangleCount <= 1
    3026         if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
    3027           AddTriangleLine(TPS[0], TPS[1], 0);
    3028           AddTriangleLine(TPS[0], TPS[2], 1);
    3029           AddTriangleLine(TPS[1], TPS[2], 2);
    3030 
    3031           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    3032           AddTriangle();
    3033 
    3034           (*it)->OtherOptCenter.Scale(-1.);
    3035           BTS->GetNormalVector((*it)->OtherOptCenter);
    3036           (*it)->OtherOptCenter.Scale(-1.);
    3037 
    3038           cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
    3039           << " for this triangle ... " << endl;
    3040           cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
    3041         } else {
    3042           cout << Verbose(1) << "WARNING: This triangle consisting of ";
    3043           cout << *(*it)->point << ", ";
    3044           cout << *BaseRay->endpoints[0]->node << " and ";
    3045           cout << *BaseRay->endpoints[1]->node << " ";
    3046           cout << "exists and is not added, as it does not seem helpful!" << endl;
    3047           result = false;
    3048         }
    3049     } else {
    3050       cout << Verbose(1) << "This triangle consisting of ";
    3051       cout << *(*it)->point << ", ";
    3052       cout << *BaseRay->endpoints[0]->node << " and ";
    3053       cout << *BaseRay->endpoints[1]->node << " ";
    3054       cout << "is invalid!" << endl;
    3055       result = false;
    3056     }
    3057 
    3058     if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
    3059       sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
    3060       if (DoTecplotOutput) {
    3061         string NameofTempFile(tempbasename);
    3062         NameofTempFile.append(NumberName);
    3063         for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
    3064         NameofTempFile.erase(npos, 1);
    3065         NameofTempFile.append(TecplotSuffix);
    3066         cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
    3067         tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
    3068         write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
    3069         tempstream->close();
    3070         tempstream->flush();
    3071         delete(tempstream);
    3072       }
    3073 
    3074       if (DoRaster3DOutput) {
    3075         string NameofTempFile(tempbasename);
    3076         NameofTempFile.append(NumberName);
    3077         for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
    3078         NameofTempFile.erase(npos, 1);
    3079         NameofTempFile.append(Raster3DSuffix);
    3080         cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
    3081         tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
    3082         write_raster3d_file(out, tempstream, this, mol);
    3083         // include the current position of the virtual sphere in the temporary raster3d file
    3084         // make the circumsphere's center absolute again
    3085         helper.CopyVector(&BaseRay->endpoints[0]->node->x);
    3086         helper.AddVector(&BaseRay->endpoints[1]->node->x);
    3087         helper.Scale(0.5);
    3088         (*it)->OptCenter.AddVector(&helper);
    3089         Vector *center = mol->DetermineCenterOfAll(out);
    3090         (*it)->OptCenter.SubtractVector(center);
    3091         delete(center);
    3092         // and add to file plus translucency object
    3093         *tempstream << "# current virtual sphere\n";
    3094         *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
    3095         *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
    3096           << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
    3097           << "\t" << RADIUS << "\t1 0 0\n";
    3098         *tempstream << "9\n  terminating special property\n";
    3099         tempstream->close();
    3100         tempstream->flush();
    3101         delete(tempstream);
    3102       }
    3103       if (DoTecplotOutput || DoRaster3DOutput)
    3104         TriangleFilesWritten++;
    3105     }
    3106 
    3107     // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
    3108     BaseRay = BLS[0];
    3109   }
    3110 
    3111   // remove all candidates from the list and then the list itself
    3112   class CandidateForTesselation *remover = NULL;
    3113   for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    3114     remover = *it;
    3115     delete(remover);
    3116   }
    3117   delete(Opt_Candidates);
    3118   cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
    3119   return result;
    3120 };
    3121 
    3122 /**
    3123  * Sort function for the candidate list.
    3124  */
    3125 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
    3126   // TODO: use get angle and remove duplicate code
    3127   Vector BaseLineVector, OrthogonalVector, helper;
    3128         if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    3129           cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    3130           //return false;
    3131           exit(1);
    3132         }
    3133         // create baseline vector
    3134         BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
    3135         BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
    3136         BaseLineVector.Normalize();
    3137 
    3138   // 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!)
    3139   helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x));
    3140   helper.SubtractVector(&(candidate1->point->x));
    3141   OrthogonalVector.CopyVector(&helper);
    3142   helper.VectorProduct(&BaseLineVector);
    3143   OrthogonalVector.SubtractVector(&helper);
    3144   OrthogonalVector.Normalize();
    3145 
    3146   // calculate both angles and correct with in-plane vector
    3147   helper.CopyVector(&(candidate1->point->x));
    3148   helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
    3149   double phi = BaseLineVector.Angle(&helper);
    3150   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    3151     phi = 2.*M_PI - phi;
    3152   }
    3153   helper.CopyVector(&(candidate2->point->x));
    3154   helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
    3155   double psi = BaseLineVector.Angle(&helper);
    3156   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    3157     psi = 2.*M_PI - psi;
    3158   }
    3159 
    3160   cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
    3161   cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
    3162 
    3163   // return comparison
    3164   return phi < psi;
    3165 }
    3166 
    3167 /**
    3168  * Checks whether the provided atom is within the tesselation structure.
    3169  *
    3170  * @param atom of which to check the position
    3171  * @param tesselation structure
    3172  *
    3173  * @return true if the atom is inside the tesselation structure, false otherwise
    3174  */
    3175 bool IsInnerAtom(atom *Atom, class Tesselation *Tess, LinkedCell* LC) {
    3176   if (Tess->LinesOnBoundary.begin() == Tess->LinesOnBoundary.end()) {
    3177     cout << Verbose(0) << "Error: There is no tesselation structure to compare the atom with, "
    3178         << "please create one first.";
    3179     exit(1);
    3180   }
    3181 
    3182   class atom *trianglePoints[3];
    3183   trianglePoints[0] = findClosestAtom(Atom, LC);
    3184   // check whether closest atom is "too close" :), then it's inside
    3185   if (trianglePoints[0]->x.DistanceSquared(&Atom->x) < MYEPSILON)
    3186     return true;
    3187   list<atom*> *connectedClosestAtoms = Tess->getClosestConnectedAtoms(trianglePoints[0], Atom);
    3188   trianglePoints[1] = connectedClosestAtoms->front();
    3189   trianglePoints[2] = connectedClosestAtoms->back();
    3190   for (int i=0;i<3;i++) {
    3191     if (trianglePoints[i] == NULL) {
    3192       cout << Verbose(1) << "IsInnerAtom encounters serious error, point " << i << " not found." << endl;
    3193     }
    3194 
    3195     cout << Verbose(1) << "List of possible atoms:" << endl;
    3196     cout << *trianglePoints[i] << endl;
    3197 
    3198 //    for(list<atom*>::iterator runner = connectedClosestAtoms->begin(); runner != connectedClosestAtoms->end(); runner++)
    3199 //      cout << Verbose(2) << *(*runner) << endl;
    3200   }
    3201   delete(connectedClosestAtoms);
    3202 
    3203   list<BoundaryTriangleSet*> *triangles = Tess->FindTriangles(trianglePoints);
    3204 
    3205   if (triangles->empty()) {
    3206     cout << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
    3207     exit(1);
    3208   }
    3209 
    3210   Vector helper;
    3211   helper.CopyVector(&Atom->x);
    3212 
    3213   // Only in case of degeneration, there will be two different scalar products.
    3214   // If at least one scalar product is positive, the point is considered to be outside.
    3215   bool status = (helper.ScalarProduct(&triangles->front()->NormalVector) < 0)
    3216       && (helper.ScalarProduct(&triangles->back()->NormalVector) < 0);
    3217   delete(triangles);
    3218   return status;
    3219 }
    3220 
    3221 /**
    3222  * Finds the atom which is closest to the provided one.
    3223  *
    3224  * @param atom to which to find the closest other atom
    3225  * @param linked cell structure
    3226  *
    3227  * @return atom which is closest to the provided one
    3228  */
    3229 atom* findClosestAtom(const atom* Atom, LinkedCell* LC) {
    3230   LinkedAtoms *List = NULL;
    3231   atom* closestAtom = NULL;
    3232   double distance = 1e16;
    3233   Vector helper;
    3234   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
    3235 
    3236   LC->SetIndexToVector(&Atom->x); // ignore status as we calculate bounds below sensibly
    3237   for(int i=0;i<NDIM;i++) // store indices of this cell
    3238     N[i] = LC->n[i];
    3239   //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    3240 
    3241   LC->GetNeighbourBounds(Nlower, Nupper);
    3242   //cout << endl;
    3243   for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    3244     for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    3245       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    3246         List = LC->GetCurrentCell();
    3247         //cout << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    3248         if (List != NULL) {
    3249           for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    3250             helper.CopyVector(&Atom->x);
    3251             helper.SubtractVector(&(*Runner)->x);
    3252             double currentNorm = helper. Norm();
    3253             if (currentNorm < distance) {
    3254               distance = currentNorm;
    3255               closestAtom = (*Runner);
    3256             }
    3257           }
    3258         } else {
    3259           cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
    3260             << LC->n[2] << " is invalid!" << endl;
    3261         }
    3262       }
    3263 
    3264   return closestAtom;
    3265 }
    3266 
    3267 /**
    3268  * Gets all atoms connected to the provided atom by triangulation lines.
    3269  *
    3270  * @param atom of which get all connected atoms
    3271  * @param atom to be checked whether it is an inner atom
    3272  *
    3273  * @return list of the two atoms linked to the provided one and closest to the atom to be checked,
    3274  */
    3275 list<atom*> * Tesselation::getClosestConnectedAtoms(atom* Atom, atom* AtomToCheck) {
    3276   list<atom*> connectedAtoms;
    3277   map<double, atom*> anglesOfAtoms;
    3278   map<double, atom*>::iterator runner;
    3279   LineMap::iterator findLines = LinesOnBoundary.begin();
    3280   list<atom*>::iterator listRunner;
    3281   Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
    3282   atom* current;
    3283   bool takeAtom = false;
    3284 
    3285   planeNorm.CopyVector(&Atom->x);
    3286   planeNorm.SubtractVector(&AtomToCheck->x);
    3287   planeNorm.Normalize();
    3288 
    3289   while (findLines != LinesOnBoundary.end()) {
    3290     takeAtom = false;
    3291 
    3292     if (findLines->second->endpoints[0]->Nr == Atom->nr) {
    3293       takeAtom = true;
    3294       current = findLines->second->endpoints[1]->node;
    3295     } else if (findLines->second->endpoints[1]->Nr == Atom->nr) {
    3296       takeAtom = true;
    3297       current = findLines->second->endpoints[0]->node;
    3298     }
    3299 
    3300     if (takeAtom) {
    3301       connectedAtoms.push_back(current);
    3302       currentPoint.CopyVector(&current->x);
    3303       currentPoint.ProjectOntoPlane(&planeNorm);
    3304       center.AddVector(&currentPoint);
    3305     }
    3306 
    3307     findLines++;
    3308   }
    3309 
    3310   cout << "Summed vectors " << center << "; number of atoms " << connectedAtoms.size()
    3311     << "; scale factor " << 1.0/connectedAtoms.size();
    3312 
    3313   center.Scale(1.0/connectedAtoms.size());
    3314   listRunner = connectedAtoms.begin();
    3315 
    3316   cout << " calculated center " << center <<  endl;
    3317 
    3318   // construct one orthogonal vector
    3319   helper.CopyVector(&AtomToCheck->x);
    3320   helper.ProjectOntoPlane(&planeNorm);
    3321   OrthogonalVector.MakeNormalVector(&center, &helper, &(*listRunner)->x);
    3322   while (listRunner != connectedAtoms.end()) {
    3323     double angle = getAngle((*listRunner)->x, AtomToCheck->x, center, OrthogonalVector);
    3324     cout << "Calculated angle " << angle << " for atom " << **listRunner << endl;
    3325     anglesOfAtoms.insert(pair<double, atom*>(angle, (*listRunner)));
    3326     listRunner++;
    3327   }
    3328 
    3329   list<atom*> *result = new list<atom*>;
    3330   runner = anglesOfAtoms.begin();
    3331   cout << "First value is " << *runner->second << endl;
    3332   result->push_back(runner->second);
    3333   runner = anglesOfAtoms.end();
    3334   runner--;
    3335   cout << "Second value is " << *runner->second << endl;
    3336   result->push_back(runner->second);
    3337 
    3338   cout << "List of closest atoms has " << result->size() << " elements, which are "
    3339     << *(result->front()) << " and " << *(result->back()) << endl;
    3340 
    3341   return result;
    3342 }
    3343 
    3344 /**
    3345  * Finds triangles belonging to the three provided atoms.
    3346  *
    3347  * @param atom list, is expected to contain three atoms
    3348  *
    3349  * @return triangles which belong to the provided atoms, will be empty if there are none,
    3350  *         will usually be one, in case of degeneration, there will be two
    3351  */
    3352 list<BoundaryTriangleSet*> *Tesselation::FindTriangles(atom* Points[3]) {
    3353   list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
    3354   LineMap::iterator FindLine;
    3355   PointMap::iterator FindPoint;
    3356   TriangleMap::iterator FindTriangle;
    3357   class BoundaryPointSet *TrianglePoints[3];
    3358 
    3359   for (int i = 0; i < 3; i++) {
    3360     FindPoint = PointsOnBoundary.find(Points[i]->nr);
    3361     if (FindPoint != PointsOnBoundary.end()) {
    3362       TrianglePoints[i] = FindPoint->second;
    3363     } else {
    3364       TrianglePoints[i] = NULL;
    3365     }
    3366   }
    3367 
    3368   // checks lines between the points in the Points for their adjacent triangles
    3369   for (int i = 0; i < 3; i++) {
    3370     if (TrianglePoints[i] != NULL) {
    3371       for (int j = i; j < 3; j++) {
    3372         if (TrianglePoints[j] != NULL) {
    3373           FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr);
    3374           if (FindLine != TrianglePoints[i]->lines.end()) {
    3375             for (; FindLine->first == TrianglePoints[j]->node->nr; FindLine++) {
    3376               FindTriangle = FindLine->second->triangles.begin();
    3377               for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
    3378                 if ((
    3379                   (FindTriangle->second->endpoints[0] == TrianglePoints[0])
    3380                     || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
    3381                     || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
    3382                   ) && (
    3383                     (FindTriangle->second->endpoints[1] == TrianglePoints[0])
    3384                     || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
    3385                     || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
    3386                   ) && (
    3387                     (FindTriangle->second->endpoints[2] == TrianglePoints[0])
    3388                     || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
    3389                     || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
    3390                   )
    3391                 ) {
    3392                   result->push_back(FindTriangle->second);
    3393                 }
    3394               }
    3395             }
    3396             // Is it sufficient to consider one of the triangle lines for this.
    3397             return result;
    3398 
    3399           }
    3400         }
    3401       }
    3402     }
    3403   }
    3404 
    3405   return result;
    3406 }
    3407 
    3408 /**
    3409  * Gets the angle between a point and a reference relative to the provided center.
    3410  *
    3411  * @param point to calculate the angle for
    3412  * @param reference to which to calculate the angle
    3413  * @param center for which to calculate the angle between the vectors
    3414  * @param OrthogonalVector helps discern between [0,pi] and [pi,2pi] in acos()
    3415  *
    3416  * @return angle between point and reference
    3417  */
    3418 double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector) {
    3419   Vector ReferenceVector, helper;
    3420   cout << Verbose(2) << reference << " is our reference " << endl;
    3421   cout << Verbose(2) << center << " is our center " << endl;
    3422   // create baseline vector
    3423   ReferenceVector.CopyVector(&reference);
    3424   ReferenceVector.SubtractVector(&center);
    3425   if (ReferenceVector.IsNull())
    3426     return M_PI;
    3427 
    3428   // calculate both angles and correct with in-plane vector
    3429   helper.CopyVector(&point);
    3430   helper.SubtractVector(&center);
    3431   if (helper.IsNull())
    3432     return M_PI;
    3433   double phi = ReferenceVector.Angle(&helper);
    3434   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    3435     phi = 2.*M_PI - phi;
    3436   }
    3437 
    3438   cout << Verbose(2) << point << " has angle " << phi << endl;
    3439 
    3440   return phi;
    3441 }
    3442 
    3443 /**
    3444  * Checks whether the provided point is within the tesselation structure.
    3445  *
    3446  * This is a wrapper function for IsInnerAtom, so it can be used with a simple
    3447  * vector, too.
    3448  *
    3449  * @param point of which to check the position
    3450  * @param tesselation structure
    3451  *
    3452  * @return true if the point is inside the tesselation structure, false otherwise
    3453  */
    3454 bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC) {
    3455   atom *temp_atom = new atom;
    3456   bool status = false;
    3457   temp_atom->x.CopyVector(&Point);
    3458 
    3459   status = IsInnerAtom(temp_atom, Tess, LC);
    3460   delete(temp_atom);
    3461 
    3462   return status;
    3463 }
    3464 
    3465817/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
    3466818 * \param *out output stream for debugging
    … …  
    3476828  bool freeTess = false;
    3477829  bool freeLC = false;
     830  ofstream *tempstream = NULL;
     831  char NumberName[255];
     832  int TriangleFilesWritten = 0;
     833
    3478834  *out << Verbose(1) << "Entering search for non convex hull. " << endl;
    3479835  if (Tess == NULL) {
    … …  
    3493849  }
    3494850
    3495   Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
     851  Tess->Find_starting_triangle(out, RADIUS, LCList);
    3496852
    3497853  baseline = Tess->LinesOnBoundary.begin();
    3498854  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
    3499855    if (baseline->second->TrianglesCount == 1) {
    3500       failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
     856      failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
    3501857      flag = flag || failflag;
    3502858      if (!failflag)
    … …  
    3514870      flag = false;
    3515871    }
    3516   }
     872
     873    // write temporary envelope
     874    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
     875      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
     876      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
     877      if (DoTecplotOutput) {
     878        string NameofTempFile(filename);
     879        NameofTempFile.append(NumberName);
     880        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
     881        NameofTempFile.erase(npos, 1);
     882        NameofTempFile.append(TecplotSuffix);
     883        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
     884        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
     885        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
     886        tempstream->close();
     887        tempstream->flush();
     888        delete(tempstream);
     889      }
     890
     891      if (DoRaster3DOutput) {
     892        string NameofTempFile(filename);
     893        NameofTempFile.append(NumberName);
     894        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
     895        NameofTempFile.erase(npos, 1);
     896        NameofTempFile.append(Raster3DSuffix);
     897        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
     898        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
     899        write_raster3d_file(out, tempstream, Tess, mol);
     900//        // include the current position of the virtual sphere in the temporary raster3d file
     901//        // make the circumsphere's center absolute again
     902//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
     903//        helper.AddVector(BaseRay->endpoints[1]->node->node);
     904//        helper.Scale(0.5);
     905//        (*it)->OptCenter.AddVector(&helper);
     906//        Vector *center = mol->DetermineCenterOfAll(out);
     907//        (*it)->OptCenter.SubtractVector(center);
     908//        delete(center);
     909//        // and add to file plus translucency object
     910//        *tempstream << "# current virtual sphere\n";
     911//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
     912//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
     913//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
     914//          << "\t" << RADIUS << "\t1 0 0\n";
     915//        *tempstream << "9\n  terminating special property\n";
     916        tempstream->close();
     917        tempstream->flush();
     918        delete(tempstream);
     919      }
     920      if (DoTecplotOutput || DoRaster3DOutput)
     921        TriangleFilesWritten++;
     922    }
     923  }
     924
     925  // write final envelope
    3517926  if (filename != 0) {
    3518927    *out << Verbose(1) << "Writing final tecplot file\n";
    … …  
    3553962  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
    3554963    << "for vector " << x << "." << endl;
    3555   atom* a = Tess->PointsOnBoundary.begin()->second->node;
     964  TesselPoint* a = Tess->PointsOnBoundary.begin()->second->node;
    3556965  cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
    3557   cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
     966  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerPoint(a, Tess, LCList)
    3558967    << "for atom " << a << " (on boundary)." << endl;
    3559   LinkedAtoms *List = NULL;
     968  LinkedNodes *List = NULL;
    3560969  for (int i=0;i<NDIM;i++) { // each axis
    3561970    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
    … …  
    3565974        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
    3566975        if (List != NULL) {
    3567           for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
     976          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
    3568977            if (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
    3569978              a = *Runner;
    … …  
    3577986      }
    3578987  }
    3579   cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
     988  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(a, Tess, LCList)
    3580989    << "for atom " << a << " (inside)." << endl;
    3581990

  • src/boundary.hpp

    r0dbddc rab1932  
    11#ifndef BOUNDARY_HPP_
    22#define BOUNDARY_HPP_
    3 
    4 class BoundaryPointSet;
    5 class BoundaryLineSet;
    6 class BoundaryTriangleSet;
    7 class CandidateForTesselation;
    83
    94// include config.h
    … …  
    149// STL headers
    1510#include <map>
    16 #include <set>
    17 #include <deque>
    1811
    19 #include <gsl/gsl_poly.h>
    20 
     12#include "config.hpp"
    2113#include "linkedcell.hpp"
    2214#include "molecules.hpp"
     15#include "tesselation.hpp"
    2316
    24 template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
    25 {
    26   if (endpoint1->Nr < endpoint2->Nr) {
    27     endpoints[0] = endpoint1;
    28     endpoints[1] = endpoint2;
    29   } else {
    30     endpoints[0] = endpoint2;
    31     endpoints[1] = endpoint1;
    32   }
    33 };
     17#define DEBUG 1
     18#define DoSingleStepOutput 0
     19#define DoTecplotOutput 1
     20#define DoRaster3DOutput 1
     21#define DoVRMLOutput 1
     22#define TecplotSuffix ".dat"
     23#define Raster3DSuffix ".r3d"
     24#define VRMLSUffix ".wrl"
    3425
    35 class BoundaryPointSet {
    36   public:
    37     BoundaryPointSet();
    38     BoundaryPointSet(atom *Walker);
    39     ~BoundaryPointSet();
     26#define DistancePair pair < double, atom* >
     27#define DistanceMap multimap < double, atom* >
     28#define DistanceTestPair pair < DistanceMap::iterator, bool>
    4029
    41     void AddLine(class BoundaryLineSet *line);
    42 
    43     LineMap lines;
    44     int LinesCount;
    45     atom *node;
    46     int Nr;
    47 };
    48 
    49 class BoundaryLineSet {
    50   public:
    51     BoundaryLineSet();
    52     BoundaryLineSet(class BoundaryPointSet *Point[2], int number);
    53     ~BoundaryLineSet();
    54 
    55     void AddTriangle(class BoundaryTriangleSet *triangle);
    56     bool IsConnectedTo(class BoundaryLineSet *line);
    57     bool ContainsBoundaryPoint(class BoundaryPointSet *point);
    58     bool CheckConvexityCriterion(ofstream *out);
    59 
    60     class BoundaryPointSet *endpoints[2];
    61     TriangleMap triangles;
    62     int TrianglesCount;
    63     int Nr;
    64 };
    65 
    66 class BoundaryTriangleSet {
    67   public:
    68     BoundaryTriangleSet();
    69     BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
    70     ~BoundaryTriangleSet();
    71 
    72     void GetNormalVector(Vector &NormalVector);
    73     bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
    74     bool ContainsBoundaryLine(class BoundaryLineSet *line);
    75     bool ContainsBoundaryPoint(class BoundaryPointSet *point);
    76 
    77     class BoundaryPointSet *endpoints[3];
    78     class BoundaryLineSet *lines[3];
    79     Vector NormalVector;
    80     int Nr;
    81 };
    82 
    83 
    84 class CandidateForTesselation {
    85   public :
    86   CandidateForTesselation(atom* candidate, BoundaryLineSet* currentBaseLine, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
    87   ~CandidateForTesselation();
    88   atom *point;
    89   BoundaryLineSet *BaseLine;
    90   Vector OptCenter;
    91   Vector OtherOptCenter;
    92 };
    93 
    94 
    95 class Tesselation {
    96   public:
    97 
    98     Tesselation();
    99     ~Tesselation();
    100 
    101     void TesselateOnBoundary(ofstream *out, molecule *mol);
    102     void GuessStartingTriangle(ofstream *out);
    103     void AddPoint(atom * Walker);
    104     void AddTrianglePoint(atom* Candidate, int n);
    105     void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
    106     void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
    107     void AddTriangle();
    108     void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC);
    109     bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC);
    110     int CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]);
    111     void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol);
    112     class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, Vector *x);
    113     bool IsInside(Vector *pointer);
    114     bool InsertStraddlingPoints(ofstream *out, molecule *mol);
    115     bool CorrectConcaveBaselines(ofstream *out);
    116     list<atom*> *getClosestConnectedAtoms(atom* Atom, atom* AtomToCheck);
    117     list<BoundaryTriangleSet*> *FindTriangles(atom* TrianglePoints[3]);
    118 
    119     PointMap PointsOnBoundary;
    120     LineMap LinesOnBoundary;
    121     TriangleMap TrianglesOnBoundary;
    122     class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
    123     class BoundaryPointSet *BPS[2];
    124     class BoundaryLineSet *BLS[3];
    125     class BoundaryTriangleSet *BTS;
    126     int PointsOnBoundaryCount;
    127     int LinesOnBoundaryCount;
    128     int TrianglesOnBoundaryCount;
    129     int TriangleFilesWritten;
    130 };
    131 
    132 
    133 ostream & operator << (ostream &ost, BoundaryPointSet &a);
    134 ostream & operator << (ostream &ost, BoundaryLineSet &a);
    135 ostream & operator << (ostream &ost, BoundaryTriangleSet &a);
    136 
     30#define Boundaries map <double, DistancePair >
     31#define BoundariesPair pair<double, DistancePair >
     32#define BoundariesTestPair pair< Boundaries::iterator, bool>
    13733
    13834double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration);
    … …  
    14339void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
    14440void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
    145 bool Choose_preferable_third_point(atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, atom *ThirdNode, Tesselation *Tess);
    146 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4);
    147 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2);
    148 bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC);
    149 bool IsInnerAtom(atom *Atom, class Tesselation *Tess, LinkedCell* LC);
    150 atom* findClosestAtom(const atom* Atom, LinkedCell* LC);
    151 double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector);
     41Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
    15242
    15343#endif /*BOUNDARY_HPP_*/

  • src/builder.cpp

    r0dbddc rab1932  
    12291229  }
    12301230  molecules->SimpleMultiAdd(mol, src, N);
    1231   delete[](src);
     1231  delete(src);
     1232
    12321233  // ... and translate back
    12331234  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {

  • src/config.cpp

    r0dbddc rab1932  
    55 */
    66
    7 #include "molecules.hpp"
     7#include "config.hpp"
    88
    99/******************************** Functions for class ConfigFileBuffer **********************/
    … …  
    276276
    277277/** Readin of Thermostat related values from parameter file.
    278  * \param *source parameter file
    279  */
    280 void config::InitThermostats(ifstream *source)
     278 * \param *fb file buffer containing the config file
     279 */
     280void config::InitThermostats(class ConfigFileBuffer *fb)
    281281{
    282282  char *thermo = MallocString(12, "IonsInitRead: thermo");
    … …  
    284284
    285285  // read desired Thermostat from file along with needed additional parameters
    286   if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
     286  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    287287    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
    288288      if (ThermostatImplemented[0] == 1) {
    … …  
    295295      if (ThermostatImplemented[1] == 1) {
    296296        Thermostat = Woodcock;
    297         ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
     297        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
    298298      } else {
    299299        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
    … …  
    303303      if (ThermostatImplemented[2] == 1) {
    304304        Thermostat = Gaussian;
    305         ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
     305        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
    306306      } else {
    307307        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
    … …  
    311311      if (ThermostatImplemented[3] == 1) {
    312312        Thermostat = Langevin;
    313         ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
    314         if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
     313        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
     314        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
    315315          cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
    316316        } else {
    … …  
    324324      if (ThermostatImplemented[4] == 1) {
    325325        Thermostat = Berendsen;
    326         ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
     326        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
    327327      } else {
    328328        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
    … …  
    332332      if (ThermostatImplemented[5] == 1) {
    333333        Thermostat = NoseHoover;
    334         ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
     334        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
    335335        alpha = 0.;
    336336      } else {
    … …  
    704704    return;
    705705  }
     706  file->close();
     707  delete(file);
    706708  RetrieveConfigPathAndName(filename);
    707709
    … …  
    722724  double value[3];
    723725 
    724   InitThermostats(file);
     726  InitThermostats(FileBuffer);
    725727 
    726728  /* Namen einlesen */

  • src/element.cpp

    r0dbddc rab1932  
    55 */
    66
    7 #include "periodentafel.hpp"
     7#include <iomanip>
     8#include <fstream>
     9
     10#include "element.hpp"
    811
    912/************************************* Functions for class element **********************************/

  • src/ellipsoid.cpp

    r0dbddc rab1932  
    66 */
    77
     8#include <gsl/gsl_multimin.h>
     9#include <gsl/gsl_vector.h>
     10
     11#include "boundary.hpp"
    812#include "ellipsoid.hpp"
    913
    … …  
    221225  set<int>::iterator current;
    222226  int index;
    223   atom *Candidate = NULL;
    224   LinkedAtoms *List = NULL;
     227  TesselPoint *Candidate = NULL;
     228  LinkedNodes *List = NULL;
    225229  *out << Verbose(2) << "Begin of PickRandomPointSet" << endl;
    226230
    … …  
    288292//            else
    289293//              *out << Verbose(2) << "Cell is empty ... " << endl;
    290             for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     294            for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    291295              if ((current != PickedAtomNrs.end()) && (*current == index)) {
    292296                Candidate = (*Runner);
    293297                *out << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl;
    294                 x[PointsPicked++].CopyVector(&(Candidate->x));    // we have one more atom picked
     298                x[PointsPicked++].CopyVector(Candidate->node);    // we have one more atom picked
    295299                current++;    // next pre-picked atom
    296300              }
    … …  
    338342      //*out << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
    339343      if (value > threshold) {
    340         x[PointsPicked].CopyVector(&(Runner->second->node->x));
     344        x[PointsPicked].CopyVector(Runner->second->node->node);
    341345        PointsPicked++;
    342346        //*out << "IN." << endl;
    … …  
    375379  Center.Zero();
    376380  for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
    377     Center.AddVector(&Runner->second->node->x);
     381    Center.AddVector(Runner->second->node->node);
    378382  Center.Scale(1./T->PointsOnBoundaryCount);
    379383  *out << Verbose(1) << "Center is at " << Center << "." << endl;

  • src/ellipsoid.hpp

    r0dbddc rab1932  
    1414#endif
    1515
    16 #include <cstdlib>
    17 
    18 #include "boundary.hpp"
     16#include "linkedcell.hpp"
    1917#include "vector.hpp"
    2018

  • src/helpers.hpp

    r0dbddc rab1932  
    88
    99using namespace std;
     10
     11// include config.h
     12#ifdef HAVE_CONFIG_H
     13#include <config.h>
     14#endif
    1015
    1116#include <iostream>
    … …  
    1520#include <math.h>
    1621#include <string>
    17 #include <string.h>
    18 #include <stdio.h>
    19 #include <stdlib.h>
    20 #include <time.h>
    2122
    2223#include "defs.hpp"
    23 
    24 // include config.h
    25 #ifdef HAVE_CONFIG_H
    26 #include <config.h>
    27 #endif
     24#include "verbose.hpp"
    2825
    2926/********************************************** helpful functions *********************************/
    … …  
    275272};
    276273
    277 /************************************* Class Verbose & Binary *******************************/
    278 
    279 /** Verbose is an IO manipulator, that writes tabs according to \a Verbosity level.
    280  */
    281 class Verbose
    282 {
    283   public:
    284     Verbose(int value) : Verbosity(value) { }
    285 
    286     ostream& print (ostream &ost) const;
    287   private:
    288     int Verbosity;
    289 };
    290 
    291 ostream& operator<<(ostream& ost,const Verbose& m);
    292 
    293 /** Binary is an IO manipulator, that writes 0 and 1 according to number \a Binary.
    294  */
    295 class Binary
    296 {
    297   public:
    298     Binary(int value) : BinaryNumber(value) { }
    299 
    300     ostream& print (ostream &ost) const;
    301   private:
    302     int BinaryNumber;
    303 };
    304 
    305 ostream& operator<<(ostream& ost,const Binary& m);
    306 
    307 
    308274
    309275#endif /*HELPERS_HPP_*/

  • src/linkedcell.cpp

    r0dbddc rab1932  
     1/** \file linkedcell.cpp
     2 *
     3 * Function implementations for the class LinkedCell.
     4 *
     5 */
     6
     7
    18#include "linkedcell.hpp"
    29#include "molecules.hpp"
     10#include "tesselation.hpp"
     11
     12// ========================================================= class LinkedCell ===========================================
     13
    314
    415/** Constructor for class LinkedCell.
    … …  
    1627
    1728/** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
    18  * \param *mol molecule structure with all Atom's
     29 * \param *set LCNodeSet class with all LCNode's
    1930 * \param RADIUS edge length of cells
    2031 */
    21 LinkedCell::LinkedCell(molecule *mol, double radius)
    22 {
    23   atom *Walker = NULL;
     32LinkedCell::LinkedCell(PointCloud *set, double radius)
     33{
     34  TesselPoint *Walker = NULL;
    2435
    2536  RADIUS = radius;
    … …  
    3142  min.Zero();
    3243  cout << Verbose(1) << "Begin of LinkedCell" << endl;
    33   if (mol->start->next == mol->end) {
    34     cerr << "ERROR: molecule contains no atoms!" << endl;
     44  if (set->IsEmpty()) {
     45    cerr << "ERROR: set contains no linked cell nodes!" << endl;
    3546    return;
    3647  }
    3748  // 1. find max and min per axis of atoms
    38   Walker = mol->start->next;
    39   for (int i=0;i<NDIM;i++) {
    40     max.x[i] = Walker->x.x[i];
    41     min.x[i] = Walker->x.x[i];
    42   }
    43   while (Walker != mol->end) {
     49  set->GoToFirst();
     50  Walker = set->GetPoint();
     51  for (int i=0;i<NDIM;i++) {
     52    max.x[i] = Walker->node->x[i];
     53    min.x[i] = Walker->node->x[i];
     54  }
     55  set->GoToFirst();
     56  while (!set->IsLast()) {
     57    Walker = set->GetPoint();
    4458    for (int i=0;i<NDIM;i++) {
    45       if (max.x[i] < Walker->x.x[i])
    46         max.x[i] = Walker->x.x[i];
    47       if (min.x[i] > Walker->x.x[i])
    48         min.x[i] = Walker->x.x[i];
     59      if (max.x[i] < Walker->node->x[i])
     60        max.x[i] = Walker->node->x[i];
     61      if (min.x[i] > Walker->node->x[i])
     62        min.x[i] = Walker->node->x[i];
    4963    }
    50     Walker = Walker->next;
     64    set->GoToNext();
    5165  }
    5266  cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
    5367
    54   // 2. find then umber of cells per axis
     68  // 2. find then number of cells per axis
    5569  for (int i=0;i<NDIM;i++) {
    5670    N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
    … …  
    6478    return;
    6579  }
    66   LC = new LinkedAtoms[N[0]*N[1]*N[2]];
     80  LC = new LinkedNodes[N[0]*N[1]*N[2]];
    6781  for (index=0;index<N[0]*N[1]*N[2];index++) {
    6882    LC [index].clear();
    … …  
    7286  // 4. put each atom into its respective cell
    7387  cout << Verbose(2) << "Filling cells ... ";
    74   Walker = mol->start;
    75   while (Walker->next != mol->end) {
    76     Walker = Walker->next;
     88  set->GoToFirst();
     89  while (!set->IsLast()) {
     90    Walker = set->GetPoint();
    7791    for (int i=0;i<NDIM;i++) {
    78       n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS);
     92      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
    7993    }
    8094    index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
    8195    LC[index].push_back(Walker);
    8296    //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
     97    set->GoToNext();
    8398  }
    8499  cout << "done."  << endl;
    … …  
    118133 * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
    119134 */
    120 LinkedAtoms* LinkedCell::GetCurrentCell()
     135LinkedNodes* LinkedCell::GetCurrentCell()
    121136{
    122137  if (CheckBounds()) {
    … …  
    128143};
    129144
    130 /** Calculates the index for a given atom *Walker.
    131  * \param Walker atom to set index to
     145/** Calculates the index for a given LCNode *Walker.
     146 * \param *Walker LCNode to set index tos
    132147 * \return if the atom is also found in this cell - true, else - false
    133148 */
    134 bool LinkedCell::SetIndexToAtom(const atom &Walker)
     149bool LinkedCell::SetIndexToNode(const TesselPoint *Walker)
    135150{
    136151  bool status = false;
    137152  for (int i=0;i<NDIM;i++) {
    138     n[i] = (int)floor((Walker.x.x[i] - min.x[i])/RADIUS);
     153    n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
    139154  }
    140155  index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
    141156  if (CheckBounds()) {
    142     for (LinkedAtoms::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
    143       status = status || ((*Runner) == &Walker);
     157    for (LinkedNodes::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
     158      status = status || ((*Runner) == Walker);
    144159    return status;
    145160  } else {
    146     cerr << Verbose(1) << "ERROR: Atom "<< Walker << " at " << Walker.x << " is out of bounds." << endl;
     161    cerr << Verbose(1) << "ERROR: Node at " << *Walker << " is out of bounds." << endl;
    147162    return false;
    148163  }

  • src/linkedcell.hpp

    r0dbddc rab1932  
     1/*
     2 * linkedcell.hpp
     3 *
     4 *  If the linked cell should be usable, the class has to inherit LCNodeSet and the nodes (containing the Vectors) have to inherit LCNode. This works well
     5 *  for molecule and atom classes.
     6 *
     7 *  Created on: Aug 3, 2009
     8 *      Author: heber
     9 */
     10
    111#ifndef LINKEDCELL_HPP_
    212#define LINKEDCELL_HPP_
     13
     14using namespace std;
    315
    416// include config.h
    … …  
    719#endif
    820
    9 #include "molecules.hpp"
     21#include <list>
    1022
    11 #define LinkedAtoms list <atom *>
     23#include "defs.hpp"
     24#include "helpers.hpp"
     25#include "vector.hpp"
    1226
    13 class LinkedCell{
     27class TesselPoint;
     28class PointCloud;
     29
     30#define LinkedNodes list<TesselPoint *>
     31
     32/** Linked Cell class for containing Vectors in real space efficiently.
     33 */
     34class LinkedCell {
    1435  public:
    1536    Vector max;       // upper boundary
    1637    Vector min;       // lower boundary
    17     LinkedAtoms *LC;  // linked cell list
     38    LinkedNodes *LC;  // linked cell list
    1839    double RADIUS;    // cell edge length
    1940    int N[NDIM];      // number of cells per axis
    … …  
    2243
    2344    LinkedCell();
    24     LinkedCell(molecule *mol, double RADIUS);
     45    LinkedCell(PointCloud *set, double RADIUS);
    2546    ~LinkedCell();
    26     LinkedAtoms* GetCurrentCell();
    27     bool SetIndexToAtom(const atom &Walker);
     47    LinkedNodes* GetCurrentCell();
     48    bool SetIndexToNode(const TesselPoint *Walker);
    2849    bool SetIndexToVector(const Vector *x);
    2950    bool CheckBounds();
    … …  
    3152
    3253    // not implemented yet
    33     bool AddAtom(atom *Walker);
    34     bool DeleteAtom(atom *Walker);
    35     bool MoveAtom(atom *Walker);
     54    bool AddNode(Vector *Walker);
     55    bool DeleteNode(Vector *Walker);
     56    bool MoveNode(Vector *Walker);
    3657};
    3758

  • src/moleculelist.cpp

    r0dbddc rab1932  
    55 */
    66
     7#include "config.hpp"
    78#include "molecules.hpp"
    89

  • src/molecules.cpp

    r0dbddc rab1932  
    55 */
    66
     7#include "config.hpp"
    78#include "molecules.hpp"
    89
    … …  
    4243  end->father = NULL;
    4344  link(start,end);
     45  InternalPointer = start;
    4446  // init bond chain list
    4547  first = new bond(start, end, 1, -1);
    … …  
    8284  delete(end);
    8385  delete(start);
     86};
     87
     88
     89/** Determine center of all atoms.
     90 * \param *out output stream for debugging
     91 * \return pointer to allocated with central coordinates
     92 */
     93Vector *molecule::GetCenter(ofstream *out)
     94{
     95  Vector *center = DetermineCenterOfAll(out);
     96  return center;
     97};
     98
     99/** Return current atom in the list.
     100 * \return pointer to atom or NULL if none present
     101 */
     102TesselPoint *molecule::GetPoint()
     103{
     104  if ((InternalPointer != start) && (InternalPointer != end))
     105    return InternalPointer;
     106  else
     107    return NULL;
     108};
     109
     110/** Return pointer to one after last atom in the list.
     111 * \return pointer to end marker
     112 */
     113TesselPoint *molecule::GetTerminalPoint()
     114{
     115  return end;
     116};
     117
     118/** Go to next atom.
     119 * Stops at last one.
     120 */
     121void molecule::GoToNext()
     122{
     123  if (InternalPointer->next != end)
     124    InternalPointer = InternalPointer->next;
     125};
     126
     127/** Go to previous atom.
     128 * Stops at first one.
     129 */
     130void molecule::GoToPrevious()
     131{
     132  if (InternalPointer->previous != start)
     133    InternalPointer = InternalPointer->previous;
     134};
     135
     136/** Goes to first atom.
     137 */
     138void molecule::GoToFirst()
     139{
     140  InternalPointer = start->next;
     141};
     142
     143/** Goes to last atom.
     144 */
     145void molecule::GoToLast()
     146{
     147  InternalPointer = end->previous;
     148};
     149
     150/** Checks whether we have any atoms in molecule.
     151 * \return true - no atoms, false - not empty
     152 */
     153bool molecule::IsEmpty()
     154{
     155  return (start->next == end);
     156};
     157
     158/** Checks whether we are at the last atom
     159 * \return true - current atom is last one, false - is not last one
     160 */
     161bool molecule::IsLast()
     162{
     163  return (InternalPointer->next == end);
    84164};
    85165

  • src/molecules.hpp

    r0dbddc rab1932  
    1818#include <gsl/gsl_randist.h>
    1919
    20 // STL headers
     20//// STL headers
    2121#include <map>
    2222#include <set>
    … …  
    2525#include <vector>
    2626
    27 #include "helpers.hpp"
     27#include "atom.hpp"
     28#include "bond.hpp"
     29#include "element.hpp"
     30#include "linkedcell.hpp"
    2831#include "parser.hpp"
    2932#include "periodentafel.hpp"
    3033#include "stackclass.hpp"
     34#include "tesselation.hpp"
    3135#include "vector.hpp"
    3236
    33 class atom;
    34 class bond;
    35 class config;
    3637class molecule;
    3738class MoleculeLeafClass;
    3839class MoleculeListClass;
    39 class Verbose;
    40 class Tesselation;
    4140
    4241/******************************** Some definitions for easier reading **********************************/
    … …  
    5049#define GraphTestPair pair<Graph::iterator, bool>
    5150
     51#define MoleculeList list <molecule *>
     52#define MoleculeListTest pair <MoleculeList::iterator, bool>
     53
    5254#define DistancePair pair < double, atom* >
    5355#define DistanceMap multimap < double, atom* >
    5456#define DistanceTestPair pair < DistanceMap::iterator, bool>
    5557
    56 #define Boundaries map <double, DistancePair >
    57 #define BoundariesPair pair<double, DistancePair >
    58 #define BoundariesTestPair pair< Boundaries::iterator, bool>
    59 
    60 #define PointMap map < int, class BoundaryPointSet * >
    61 #define PointPair pair < int, class BoundaryPointSet * >
    62 #define PointTestPair pair < PointMap::iterator, bool >
    63 #define CandidateList list <class CandidateForTesselation *>
    64 
    65 #define LineMap multimap < int, class BoundaryLineSet * >
    66 #define LinePair pair < int, class BoundaryLineSet * >
    67 #define LineTestPair pair < LineMap::iterator, bool >
    68 
    69 #define TriangleMap map < int, class BoundaryTriangleSet * >
    70 #define TrianglePair pair < int, class BoundaryTriangleSet * >
    71 #define TriangleTestPair pair < TrianglePair::iterator, bool >
    72 
    73 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
    74 #define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
    75 
    76 #define MoleculeList list <molecule *>
    77 #define MoleculeListTest pair <MoleculeList::iterator, bool>
    78 
    79 #define LinkedAtoms list <atom *>
     58
     59//#define LinkedAtoms list <atom *>
    8060
    8161/******************************** Some small functions and/or structures **********************************/
    … …  
    125105};
    126106
    127 /** Single atom.
    128  * Class incoporates position, type
    129  */
    130 class atom {
    131   public:
    132     Vector x;       //!< coordinate array of atom, giving position within cell
    133     Vector v;       //!< velocity array of atom
    134     element *type;  //!< pointing to element
    135     atom *previous; //!< previous atom in molecule list
    136     atom *next;     //!< next atom in molecule list
    137     atom *father;   //!< In many-body bond order fragmentations points to originating atom
    138     atom *Ancestor; //!< "Father" in Depth-First-Search
    139     char *Name;      //!< unique name used during many-body bond-order fragmentation
    140     int FixedIon;   //!< config variable that states whether forces act on the ion or not
    141     int *sort;      //!< sort criteria
    142     int nr;         //!< continuous, unique number
    143     int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
    144     int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
    145     int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
    146     bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
    147     bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
    148     unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
    149     bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
    150 
    151   atom();
    152   atom(class atom *pointer);
    153   ~atom();
    154 
    155   bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
    156   bool OutputXYZLine(ofstream *out) const;
    157   atom *GetTrueFather();
    158   bool Compare(const atom &ptr);
    159 
    160   private:
    161 };
    162 
    163 ostream & operator << (ostream &ost, const atom &a);
    164 
    165 /** Bonds between atoms.
    166  * Class incorporates bonds between atoms in a molecule,
    167  * used to derive tge fragments in many-body bond order
    168  * calculations.
    169  */
    170 class bond {
    171   public:
    172     atom *leftatom;    //!< first bond partner
    173     atom *rightatom;  //!< second bond partner
    174     bond *previous; //!< previous atom in molecule list
    175     bond *next;     //!< next atom in molecule list
    176     int HydrogenBond;  //!< Number of hydrogen atoms in the bond
    177     int BondDegree;    //!< single, double, triple, ... bond
    178     int nr;           //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
    179     bool Cyclic;      //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
    180     enum EdgeType Type;//!< whether this is a tree or back edge
    181 
    182   atom * GetOtherAtom(atom *Atom) const;
    183   bond * GetFirstBond();
    184   bond * GetLastBond();
    185 
    186   bool MarkUsed(enum Shading color);
    187   enum Shading IsUsed();
    188   void ResetUsed();
    189   bool Contains(const atom *ptr);
    190   bool Contains(const int nr);
    191 
    192   bond();
    193   bond(atom *left, atom *right);
    194   bond(atom *left, atom *right, int degree);
    195   bond(atom *left, atom *right, int degree, int number);
    196   ~bond();
    197 
    198   private:
    199     enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
    200 };
    201 
    202 
    203 ostream & operator << (ostream &ost, const bond &b);
    204 
    205107#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
    206108enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
    … …  
    210112 * Class incorporates number of types
    211113 */
    212 class molecule {
     114class molecule : public PointCloud {
    213115  public:
    214116    double cell_size[6];//!< cell size
    … …  
    234136    char name[MAXSTRINGSIZE];         //!< arbitrary name
    235137    int IndexNr;        //!< index of molecule in a MoleculeListClass
     138    class Tesselation *TesselStruct;
    236139
    237140  molecule(periodentafel *teil);
    238   ~molecule();
     141  virtual ~molecule();
     142
     143  // re-definition of virtual functions from PointCloud
     144  Vector *GetCenter(ofstream *out);
     145  TesselPoint *GetPoint();
     146  TesselPoint *GetTerminalPoint();
     147  void GoToNext();
     148  void GoToPrevious();
     149  void GoToFirst();
     150  void GoToLast();
     151  bool IsEmpty();
     152  bool IsLast();
    239153
    240154  /// remove atoms from molecule.
    … …  
    350264  private:
    351265  int last_atom;      //!< number given to last atom
     266  atom *InternalPointer;  //!< internal pointer for PointCloud
    352267};
    353268
    … …  
    423338
    424339
    425 /** The config file.
    426  * The class contains all parameters that control a dft run also functions to load and save.
    427  */
    428 class config {
    429   public:
    430     int PsiType;
    431     int MaxPsiDouble;
    432     int PsiMaxNoUp;
    433     int PsiMaxNoDown;
    434     int MaxMinStopStep;
    435     int InitMaxMinStopStep;
    436     int ProcPEGamma;
    437     int ProcPEPsi;
    438     char *configpath;
    439     char *configname;
    440     bool FastParsing;
    441     double Deltat;
    442     string basis;
    443 
    444     char *databasepath;
    445 
    446     int DoConstrainedMD;
    447     int MaxOuterStep;
    448     int Thermostat;
    449     int *ThermostatImplemented;
    450     char **ThermostatNames;
    451     double TempFrequency;
    452     double alpha;
    453     double HooverMass;
    454     double TargetTemp;
    455     int ScaleTempStep;
    456    
    457   private:
    458     char *mainname;
    459     char *defaultpath;
    460     char *pseudopotpath;
    461 
    462     int DoOutVis;
    463     int DoOutMes;
    464     int DoOutNICS;
    465     int DoOutOrbitals;
    466     int DoOutCurrent;
    467     int DoFullCurrent;
    468     int DoPerturbation;
    469     int DoWannier;
    470     int CommonWannier;
    471     double SawtoothStart;
    472     int VectorPlane;
    473     double VectorCut;
    474     int UseAddGramSch;
    475     int Seed;
    476 
    477     int OutVisStep;
    478     int OutSrcStep;
    479     int MaxPsiStep;
    480     double EpsWannier;
    481 
    482     int MaxMinStep;
    483     double RelEpsTotalEnergy;
    484     double RelEpsKineticEnergy;
    485     int MaxMinGapStopStep;
    486     int MaxInitMinStep;
    487     double InitRelEpsTotalEnergy;
    488     double InitRelEpsKineticEnergy;
    489     int InitMaxMinGapStopStep;
    490 
    491     //double BoxLength[NDIM*NDIM];
    492 
    493     double ECut;
    494     int MaxLevel;
    495     int RiemannTensor;
    496     int LevRFactor;
    497     int RiemannLevel;
    498     int Lev0Factor;
    499     int RTActualUse;
    500     int AddPsis;
    501 
    502     double RCut;
    503     int StructOpt;
    504     int IsAngstroem;
    505     int RelativeCoord;
    506     int MaxTypes;
    507 
    508 
    509   int ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical);
    510   int ParseForParameter(int verbose, struct ConfigFileBuffer *FileBuffer, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical);
    511 
    512   public:
    513   config();
    514   ~config();
    515 
    516   int TestSyntax(char *filename, periodentafel *periode, molecule *mol);
    517   void Load(char *filename, periodentafel *periode, molecule *mol);
    518   void LoadOld(char *filename, periodentafel *periode, molecule *mol);
    519   void RetrieveConfigPathAndName(string filename);
    520   bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
    521   bool SaveMPQC(const char *filename, molecule *mol) const;
    522   void Edit();
    523   bool GetIsAngstroem() const;
    524   char *GetDefaultPath() const;
    525   void SetDefaultPath(const char *path);
    526   void InitThermostats(ifstream *source);
    527 };
    528 
    529340#endif /*MOLECULES_HPP_*/
    530341

  • src/parser.hpp

    r0dbddc rab1932  
    5555 
    5656  MatrixContainer();
    57   ~MatrixContainer();
     57  virtual ~MatrixContainer();
    5858 
    5959  bool InitialiseIndices(class MatrixContainer *Matrix = NULL);

  • src/periodentafel.cpp

    r0dbddc rab1932  
    77using namespace std;
    88
     9#include <iomanip>
     10#include <fstream>
     11
     12#include "helpers.hpp"
    913#include "periodentafel.hpp"
     14#include "verbose.hpp"
    1015
    1116/************************************* Functions for class periodentafel ***************************/

  • src/periodentafel.hpp

    r0dbddc rab1932  
    33
    44using namespace std;
    5 
    6 #include "defs.hpp"
    7 #include "helpers.hpp"
    85
    96// include config.h
    … …  
    129#endif
    1310
     11#include <iostream>
     12
     13#include "defs.hpp"
     14#include "element.hpp"
     15
    1416// ====================================== class definitions =========================
    1517
    16 class element;
    17 class periodentafel;
    18 
    19 /** Chemical element.
    20  * Class incorporates data for a certain chemical element to be referenced from atom class.
    21  */
    22 class element {
    23   public:
    24     double mass;    //!< mass in g/mol
    25     double CovalentRadius;  //!< covalent radius
    26     double VanDerWaalsRadius;  //!< can-der-Waals radius
    27     int Z;          //!< atomic number
    28     char name[64];  //!< atom name, i.e. "Hydrogren"
    29     char symbol[3]; //!< short form of the atom, i.e. "H"
    30     char period[8];    //!< period: n quantum number
    31     char group[8];    //!< group: l quantum number
    32     char block[8];    //!< block: l quantum number
    33     element *previous;  //!< previous item in list
    34     element *next;  //!< next element in list
    35     int *sort;      //!< sorc criteria
    36     int No;         //!< number of element set on periodentafel::Output()
    37     double Valence;   //!< number of valence electrons for this element
    38     int NoValenceOrbitals;  //!< number of valence orbitals, used for determining bond degree in molecule::CreateConnectmatrix()
    39     double HBondDistance[NDIM]; //!< distance in Angstrom of this element to hydrogen  (for single, double and triple bonds)
    40     double HBondAngle[NDIM];     //!< typical angle for one, two, three bonded hydrogen (in degrees)
    41 
    42   element();
    43   ~element();
    44 
    45   //> print element entries to screen
    46   bool Output(ofstream *out) const;
    47   bool Checkout(ofstream *out, const int No, const int NoOfAtoms) const;
    48 
    49   private:
    50 };
    5118
    5219/** Periodentafel is a list of all elements sorted by their atomic number.

  • src/stackclass.hpp

    r0dbddc rab1932  
    11#ifndef STACKCLASS_HPP_
    22#define STACKCLASS_HPP_
     3
     4#include "verbose.hpp"
    35
    46template <typename T> class StackClass;

  • src/vector.cpp

    r0dbddc rab1932  
    44 *
    55 */
     6
    67
    78#include "molecules.hpp"

  • src/vector.hpp

    r0dbddc rab1932  
    11#ifndef VECTOR_HPP_
    22#define VECTOR_HPP_
     3
     4using namespace std;
     5
     6#include "helpers.hpp"
    37
    48class Vector;
Note: See TracChangeset for help on using the changeset viewer.