source: src/Tesselation/boundary.cpp@ 5f7b95

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/** \file boundary.cpp
 *
 * Implementations and super-function for envelopes
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "boundary.hpp"
#include "BoundaryLineSet.hpp"
#include "BoundaryPointSet.hpp"
#include "BoundaryTriangleSet.hpp"
#include "Box.hpp"
#include "CandidateForTesselation.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinkedCell/linkedcell.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
#include "RandomNumbers/RandomNumberGenerator.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "World.hpp"

#include <iostream>
#include <iomanip>

#include<gsl/gsl_poly.h>
#include<time.h>

// ========================================== F U N C T I O N S =================================


/** Determines greatest diameters of a cluster defined by its convex envelope.
 * Looks at lines parallel to one axis and where they intersect on the projected planes
 * \param *out output stream for debugging
 * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
 * \param *mol molecule structure representing the cluster
 * \param *&TesselStruct Tesselation structure with triangles
 * \param IsAngstroem whether we have angstroem or atomic units
 * \return NDIM array of the diameters
 */
double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
{
        //Info FunctionInfo(__func__);
  // get points on boundary of NULL was given as parameter
  bool BoundaryFreeFlag = false;
  double OldComponent = 0.;
  double tmp = 0.;
  double w1 = 0.;
  double w2 = 0.;
  Vector DistanceVector;
  Vector OtherVector;
  int component = 0;
  int Othercomponent = 0;
  Boundaries::const_iterator Neighbour;
  Boundaries::const_iterator OtherNeighbour;
  double *GreatestDiameter = new double[NDIM];

  const Boundaries *BoundaryPoints;
  if (BoundaryPtr == NULL) {
    BoundaryFreeFlag = true;
    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  } else {
    BoundaryPoints = BoundaryPtr;
    LOG(0, "Using given boundary points set.");
  }
  // determine biggest "diameter" of cluster for each axis
  for (int i = 0; i < NDIM; i++)
    GreatestDiameter[i] = 0.;
  for (int axis = 0; axis < NDIM; axis++)
    { // regard each projected plane
      //LOG(1, "Current axis is " << axis << ".");
      for (int j = 0; j < 2; j++)
        { // and for both axis on the current plane
          component = (axis + j + 1) % NDIM;
          Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
          //LOG(1, "Current component is " << component << ", Othercomponent is " << Othercomponent << ".");
          for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
              //LOG(1, "Current runner is " << *(runner->second.second) << ".");
              // seek for the neighbours pair where the Othercomponent sign flips
              Neighbour = runner;
              Neighbour++;
              if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                Neighbour = BoundaryPoints[axis].begin();
              DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
              do { // seek for neighbour pair where it flips
                  OldComponent = DistanceVector[Othercomponent];
                  Neighbour++;
                  if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                    Neighbour = BoundaryPoints[axis].begin();
                  DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
                  //LOG(2, "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << ".");
                } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                  OldComponent) - DistanceVector[Othercomponent] / fabs(
                  DistanceVector[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
              if (runner != Neighbour) {
                  OtherNeighbour = Neighbour;
                  if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
                    OtherNeighbour = BoundaryPoints[axis].end();
                  OtherNeighbour--;
                  //LOG(1, "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << ".");
                  // now we have found the pair: Neighbour and OtherNeighbour
                  OtherVector = (runner->second.second->getPosition()) - (OtherNeighbour->second.second->getPosition());
                  //LOG(1, "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << ".");
                  //LOG(1, "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << ".");
                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                  w1 = fabs(OtherVector[Othercomponent]);
                  w2 = fabs(DistanceVector[Othercomponent]);
                  tmp = fabs((w1 * DistanceVector[component] + w2
                      * OtherVector[component]) / (w1 + w2));
                  // mark if it has greater diameter
                  //LOG(1, "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << ".");
                  GreatestDiameter[component] = (GreatestDiameter[component]
                      > tmp) ? GreatestDiameter[component] : tmp;
                } //else
              //LOG(1, "Saw no sign flip, probably top or bottom node.");
            }
        }
    }
  LOG(0, "RESULT: The biggest diameters are "
      << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
      << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
      : "atomiclength") << ".");

  // free reference lists
  if (BoundaryFreeFlag)
    delete[] (BoundaryPoints);

  return GreatestDiameter;
}
;


/** Determines the boundary points of a cluster.
 * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
 * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
 * center and first and last point in the triple, it is thrown out.
 * \param *out output stream for debugging
 * \param *mol molecule structure representing the cluster
 * \param *&TesselStruct pointer to Tesselation structure
 */
Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
{
        //Info FunctionInfo(__func__);
  PointMap PointsOnBoundary;
  LineMap LinesOnBoundary;
  TriangleMap TrianglesOnBoundary;
  Vector MolCenter = mol->DetermineCenterOfAll();
  Vector helper;
  BoundariesTestPair BoundaryTestPair;
  Vector AxisVector;
  Vector AngleReferenceVector;
  Vector AngleReferenceNormalVector;
  Vector ProjectedVector;
  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, Nr)
  double angle = 0.;

  // 3a. Go through every axis
  for (int axis = 0; axis < NDIM; axis++) {
    AxisVector.Zero();
    AngleReferenceVector.Zero();
    AngleReferenceNormalVector.Zero();
    AxisVector[axis] = 1.;
    AngleReferenceVector[(axis + 1) % NDIM] = 1.;
    AngleReferenceNormalVector[(axis + 2) % NDIM] = 1.;

    LOG(1, "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << ".");

    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    // Boundaries stores non-const TesselPoint ref, hence we need iterator here
    for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
      ProjectedVector = (*iter)->getPosition() - (MolCenter);
      ProjectedVector.ProjectOntoPlane(AxisVector);

      // correct for negative side
      const double radius = ProjectedVector.NormSquared();
      if (fabs(radius) > MYEPSILON)
        angle = ProjectedVector.Angle(AngleReferenceVector);
      else
        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues

      //LOG(1, "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << ".");
      if (ProjectedVector.ScalarProduct(AngleReferenceNormalVector) > 0) {
        angle = 2. * M_PI - angle;
      }
    LOG(1, "Inserting " << **iter << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector);
      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, TesselPointDistancePair (radius, (*iter))));
      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
        LOG(2, "Encountered two vectors whose projection onto axis " << axis << " is equal: ");
        LOG(2, "Present vector: " << *BoundaryTestPair.first->second.second);
        LOG(2, "New vector: " << **iter);
        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
          BoundaryTestPair.first->second.second = (*iter);
          LOG(2, "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << ".");
        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
          helper = (*iter)->getPosition() - (MolCenter);
          const double oldhelperNorm = helper.NormSquared();
          helper = BoundaryTestPair.first->second.second->getPosition() - (MolCenter);
          if (helper.NormSquared() < oldhelperNorm) {
            BoundaryTestPair.first->second.second = (*iter);
            LOG(2, "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << ".");
          } else {
            LOG(2, "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << ".");
          }
        } else {
          LOG(2, "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << ".");
        }
      }
    }
    // printing all inserted for debugging
    //    {
    //      std::stringstream output;
    //      output << "Printing list of candidates for axis " << axis << " which we have inserted so far: ";
    //      int i=0;
    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    //        if (runner != BoundaryPoints[axis].begin())
    //          output << ", " << i << ": " << *runner->second.second;
    //        else
    //          output << i << ": " << *runner->second.second;
    //        i++;
    //      }
    //      LOG(1, output.str());
    //    }
    // 3c. throw out points whose distance is less than the mean of left and right neighbours
    bool flag = false;
    LOG(1, "Looking for candidates to kick out by convex condition ... ");
    do { // do as long as we still throw one out per round
      flag = false;
      Boundaries::iterator left = BoundaryPoints[axis].begin();
      Boundaries::iterator right = BoundaryPoints[axis].begin();
      Boundaries::iterator runner = BoundaryPoints[axis].begin();
      bool LoopOnceDone = false;
      while (!LoopOnceDone) {
        runner = right;
        right++;
        // set neighbours correctly
        if (runner == BoundaryPoints[axis].begin()) {
          left = BoundaryPoints[axis].end();
        } else {
          left = runner;
        }
        left--;
        if (right == BoundaryPoints[axis].end()) {
          right = BoundaryPoints[axis].begin();
          LoopOnceDone = true;
        }
        // check distance

        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
        {
          Vector SideA, SideB, SideC, SideH;
          SideA = left->second.second->getPosition() - (MolCenter);
          SideA.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideA: " << SideA);

          SideB = right->second.second->getPosition() -(MolCenter);
          SideB.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideB: " << SideB);

          SideC = left->second.second->getPosition() - right->second.second->getPosition();
          SideC.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideC: " << SideC);

          SideH = runner->second.second->getPosition() -(MolCenter);
          SideH.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideH: " << SideH);

          // calculate each length
          const double a = SideA.Norm();
          //const double b = SideB.Norm();
          //const double c = SideC.Norm();
          const double h = SideH.Norm();
          // calculate the angles
          const double alpha = SideA.Angle(SideH);
          const double beta = SideA.Angle(SideC);
          const double gamma = SideB.Angle(SideH);
          const double delta = SideC.Angle(SideH);
          const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
          //LOG(1, " 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 << ".");
          LOG(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 << ".");
          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
            // throw out point
            LOG(1, "Throwing out " << *runner->second.second << ".");
            BoundaryPoints[axis].erase(runner);
            runner = right;
            flag = true;
          }
        }
      }
    } while (flag);
  }
  return BoundaryPoints;
};

/** Tesselates the convex boundary by finding all boundary points.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's.
 * \param *BoundaryPts set of boundary points to use or NULL
 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList atoms in LinkedCell_deprecated list
 * \param *filename filename prefix for output of vertex data
 * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
 */
void FindConvexBorder(const molecule* mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell_deprecated *LCList, const char *filename)
{
        //Info FunctionInfo(__func__);
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

  if (TesselStruct != NULL) // free if allocated
    delete(TesselStruct);
  TesselStruct = new class Tesselation;

  // 1. Find all points on the boundary
  if (BoundaryPts == NULL) {
    BoundaryFreeFlag = true;
    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  } else {
    BoundaryPoints = BoundaryPts;
    LOG(0, "Using given boundary points set.");
  }

// printing all inserted for debugging
  if (DoLog(1)) {
    for (int axis=0; axis < NDIM; axis++) {
      std::stringstream output;
      output << "Printing list of candidates for axis " << axis << " which we have inserted so far: ";
      int i=0;
      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
        if (runner != BoundaryPoints[axis].begin())
          output << ", " << i << ": " << *runner->second.second;
        else
          output << i << ": " << *runner->second.second;
        i++;
      }
      LOG(1, output.str());
    }
  }

  // 2. fill the boundary point list
  for (int axis = 0; axis < NDIM; axis++)
    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
          LOG(2, "Point " << *(runner->second.second) << " is already present.");

  LOG(0, "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary.");
  // now we have the whole set of edge points in the BoundaryList

  // listing for debugging
  //if (DoLog(1)) {
  //  std::stringstream output;
  //  output << "Listing PointsOnBoundary:";
  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
  //    output << " " << *runner->second;
  //  }
  //  LOG(1, output.str());
  //}

  // 3a. guess starting triangle
  TesselStruct->GuessStartingTriangle();

  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
  PointCloudAdaptor< molecule > cloud(const_cast<molecule *>(mol), mol->name);
  TesselStruct->TesselateOnBoundary(cloud);

  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
  if (!TesselStruct->InsertStraddlingPoints(cloud, LCList))
    ELOG(1, "Insertion of straddling points failed!");

  LOG(0, "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points.");

  // 4. Store triangles in tecplot file
  StoreTrianglesinFile(mol, TesselStruct, filename, "_intermed");

  // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
  bool AllConvex = true;
  class BoundaryLineSet *line = NULL;
  do {
    AllConvex = true;
    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
      line = LineRunner->second;
      LOG(1, "INFO: Current line is " << *line << ".");
      if (!line->CheckConvexityCriterion()) {
        LOG(1, "... line " << *line << " is concave, flipping it.");

        // flip the line
        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
          ELOG(1, "Correction of concave baselines failed!");
        else {
          TesselStruct->FlipBaseline(line);
          LOG(1, "INFO: Correction of concave baselines worked.");
          LineRunner = TesselStruct->LinesOnBoundary.begin(); // LineRunner may have been erase if line was deleted from LinesOnBoundary
        }
      }
    }
  } while (!AllConvex);

  // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
//  if (!TesselStruct->CorrectConcaveTesselPoints(out))
//    ELOG(1, "Correction of concave tesselpoints failed!");

  LOG(0, "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points.");

  // 4. Store triangles in tecplot file
  StoreTrianglesinFile(mol, TesselStruct, filename, "");

  // free reference lists
  if (BoundaryFreeFlag)
    delete[] (BoundaryPoints);
};

/** For testing removes one boundary point after another to check for leaks.
 * \param *out output stream for debugging
 * \param *TesselStruct Tesselation containing envelope with boundary points
 * \param *mol molecule
 * \param *filename name of file
 * \return true - all removed, false - something went wrong
 */
bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
{
        //Info FunctionInfo(__func__);
  int i=0;
  char number[MAXSTRINGSIZE];

  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    ELOG(1, "TesselStruct is empty.");
    return false;
  }

  PointMap::iterator PointRunner;
  while (!TesselStruct->PointsOnBoundary.empty()) {
    if (DoLog(1)) {
      std::stringstream output;
      output << "Remaining points are: ";
      for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
        output << *(PointSprinter->second) << "\t";
      LOG(1, output.str());
    }

    PointRunner = TesselStruct->PointsOnBoundary.begin();
    // remove point
    TesselStruct->RemovePointFromTesselatedSurface(PointRunner->second);

    // store envelope
    sprintf(number, "-%04d", i++);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, number);
  }

  return true;
};

/** Creates a convex envelope from a given non-convex one.
 * -# First step, remove concave spots, i.e. singular "dents"
 *   -# We go through all PointsOnBoundary.
 *   -# We CheckConvexityCriterion() for all its lines.
 *   -# If all its lines are concave, it cannot be on the convex envelope.
 *   -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
 *   -# We calculate the additional volume.
 *   -# We go over all lines until none yields a concavity anymore.
 * -# Second step, remove concave lines, i.e. line-shape "dents"
 *   -# We go through all LinesOnBoundary
 *   -# We CheckConvexityCriterion()
 *   -# If it returns concave, we flip the line in this quadrupel of points (abusing the degeneracy of the tesselation)
 *   -# We CheckConvexityCriterion(),
 *   -# if it's concave, we continue
 *   -# if not, we mark an error and stop
 * Note: This routine - for free - calculates the difference in volume between convex and
 * non-convex envelope, as the former is easy to calculate - Tesselation::getVolumeOfConvexEnvelope() - it
 * can be used to compute volumes of arbitrary shapes.
 * \param *out output stream for debugging
 * \param *TesselStruct non-convex envelope, is changed in return!
 * \param *mol molecule
 * \param *filename name of file
 * \return volume difference between the non- and the created convex envelope
 */
double ConvexizeNonconvexEnvelope(
    Tesselation *&TesselStruct,
    const molecule * const mol,
    const char * const filename,
    bool DebugOutputEveryStep)
{
        //Info FunctionInfo(__func__);
  double volume = 0;
  class BoundaryPointSet *point = NULL;
  class BoundaryLineSet *line = NULL;
  bool Concavity = false;
  char dummy[MAXSTRINGSIZE];
  PointMap::iterator PointRunner;
  PointMap::iterator PointAdvance;
  LineMap::iterator LineRunner;
  LineMap::iterator LineAdvance;
  TriangleMap::iterator TriangleRunner;
  TriangleMap::iterator TriangleAdvance;
  int run = 0;

  // check whether there is something to work on
  if (TesselStruct == NULL) {
    ELOG(1, "TesselStruct is empty!");
    return volume;
  }

  LOG(1, "INFO: Making tesselated surface with " << TesselStruct->TrianglesOnBoundaryCount
      << " convex ...");

  // first purge all degenerate triangles
  TesselStruct->RemoveDegeneratedTriangles();

  do {
    Concavity = false;

    if (DebugOutputEveryStep) {
      sprintf(dummy, "-%d", run++);
      //CalculateConcavityPerBoundaryPoint(TesselStruct);
      LOG(1, "INFO: Writing " << run << "th tesselation file.");
      StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
    }

    // first step: remove all full-concave point
    PointRunner = TesselStruct->PointsOnBoundary.begin();
    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
    while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
      PointAdvance++;
      point = PointRunner->second;
      LOG(2, "INFO: Current point is " << *point << ".");
      // check that at least a single line is concave
      LineMap::iterator LineRunner = point->lines.begin();
      for (; LineRunner != point->lines.end(); LineRunner++) {
        const BoundaryLineSet * line = LineRunner->second;
        LOG(3, "INFO: Current line of point " << *point << " is " << *line << ".");
        if (!line->CheckConvexityCriterion())
          break;
      }
      // remove the point if needed
      if (LineRunner != point->lines.end()) {
        const double tmp = TesselStruct->RemoveFullConcavePointFromTesselatedSurface(point);
        if (tmp > 0.) {
          volume += tmp;
          Concavity = true;
        }
      }
      PointRunner = PointAdvance;
    }

    if (DebugOutputEveryStep) {
      sprintf(dummy, "-%d", run++);
      //CalculateConcavityPerBoundaryPoint(TesselStruct);
      LOG(1, "INFO: Writing " << run << "th tesselation file.");
      StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
    }

    // second step: flip baselines, i.e. add general tetraeder at concave lines
    // when the tetraeder does not intersect with other already present triangles
    LineRunner = TesselStruct->LinesOnBoundary.begin();
    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
    std::map<double, std::pair<BoundaryLineSet *, double> > GainMap;
    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
      LineAdvance++;
      line = LineRunner->second;
      if (!line->CheckConvexityCriterion()) {
        LOG(2, "INFO: concave line is " << *line << ".");
        // gather the other points
        BoundaryPointSet *BPS[4];
        int m = 0;
        {
          for (TriangleMap::iterator runner = line->triangles.begin(); runner != line->triangles.end(); runner++)
            for (int j = 0; j < 3; j++) // all of their endpoints and baselines
              if (!line->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
                BPS[m++] = runner->second->endpoints[j];
        }
        BPS[2] = line->endpoints[0];
        BPS[3] = line->endpoints[1];
        LOG(3, "DEBUG: other line would consist of " << *BPS[0] << " and "
            << *BPS[1] << ".");

        // check for already present (third) side of the tetraeder as we then
        // would create a degenerate triangle
        bool TetraederSidePresent = false;
        {
          class TesselPoint *TriangleCandidates[3];
          TriangleCandidates[0] = BPS[0]->node;
          TriangleCandidates[1] = BPS[1]->node;
          TriangleCandidates[2] = BPS[2]->node;
          if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {
            LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","
                << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");
            TetraederSidePresent = true;
          }
          TriangleCandidates[2] = BPS[3]->node;
          if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {
            LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","
                << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");
            TetraederSidePresent = true;
          }
        }

        if ((BPS[0] != BPS[1]) && (m == 2) && (!TetraederSidePresent)) {
          // check whether all adjacent triangles do not intersect with new line
          bool no_line_intersects = true;
          Vector Intersection;
          TriangleSet triangles;
          TriangleSet *firsttriangles = TesselStruct->GetAllTriangles(line->endpoints[0]);
          TriangleSet *secondtriangles = TesselStruct->GetAllTriangles(line->endpoints[1]);
          triangles.insert( firsttriangles->begin(), firsttriangles->end() );
          triangles.insert( secondtriangles->begin(), secondtriangles->end() );
          delete firsttriangles;
          delete secondtriangles;
          for (TriangleSet::const_iterator triangleiter = triangles.begin();
              triangleiter != triangles.end(); ++triangleiter) {
            const BoundaryTriangleSet * triangle = *triangleiter;
            bool line_intersects = triangle->GetIntersectionInsideTriangle(
                BPS[0]->node->getPosition(),
                BPS[1]->node->getPosition(),
                Intersection);
            // switch result when coinciding with endpoint
            bool concave_adjacent_line = false;
            bool intersection_is_endnode = false;
            for (int j=0;j<2;++j) {
              if (Intersection.DistanceSquared(BPS[j]->node->getPosition()) < MYEPSILON) {
                intersection_is_endnode = true;
                // check whether its an adjacent triangle and if it's concavely connected
                // only then are we in danger of cutting through it and need to check
                // sign of normal vector and intersecting line
                for (int i=0;i<2;++i)
                  for (int lineindex=0;lineindex < 3;++lineindex)
                    if ((triangle->lines[lineindex]->ContainsBoundaryPoint(line->endpoints[i]))
                        && (triangle->lines[lineindex]->ContainsBoundaryPoint(BPS[j]))) {
                      concave_adjacent_line = !triangle->lines[lineindex]->CheckConvexityCriterion();
                  }
                if (concave_adjacent_line) {
                  const Vector intersector =
                      BPS[(j+1)%2]->node->getPosition() - Intersection;
                  if (triangle->NormalVector.ScalarProduct(intersector) >= -MYEPSILON) {
                    LOG(4, "ACCEPT: Intersection coincides with first endpoint "
                        << *BPS[j] << ".");
                    line_intersects = false;
                  } else {
                    LOG(4, "REJECT: Intersection ends on wrong side of triangle.");
                  }
                } else {
                  LOG(4, "ACCEPT: Intersection coincides with first endpoint "
                      << *BPS[j] << ".");
                  line_intersects = false;
                }
              }
            }
            // if we have an intersection, check that it is within either
            // endpoint, i.e. check that scalar product between vectors going
            // from intersction to either endpoint has negative sign (both
            // vectors point in opposite directions)
            if (!intersection_is_endnode && line_intersects) {
              const Vector firstvector = BPS[0]->node->getPosition() - Intersection;
              const Vector secondvector = BPS[1]->node->getPosition() - Intersection;
              if (firstvector.ScalarProduct(secondvector) >= 0)
                line_intersects = false;
            }
            no_line_intersects &= !line_intersects;
          }

          if (no_line_intersects) {
            // calculate the volume
            const double tmp = line->CalculateConvexity();
            const double gain =
              CalculateVolumeofGeneralTetraeder(
                  BPS[0]->node->getPosition(),
                  BPS[1]->node->getPosition(),
                  BPS[2]->node->getPosition(),
                  BPS[3]->node->getPosition());

            GainMap.insert(std::make_pair(tmp, std::make_pair(line,gain) ));
            LOG(2, "DEBUG: Adding concave line " << *line << " with gain of "
                << gain << ".");
          } else {
            // if 2 or 3 don't
            LOG(2, "DEBUG: We don't added concave line " << *line
                << " as other line intersects with adjacent triangles.");
          }
        }
      }
      LineRunner = LineAdvance;
    }
    // flip line with most gain
    if (!GainMap.empty()) {
      line = GainMap.begin()->second.first;
      const double tmp = GainMap.begin()->second.second;
      volume += tmp;

//      GainMap.clear();

      // and flip the line
      LOG(1, "INFO: Flipping current most concave line " << *line << " with gain of "
          << tmp << ".");
      TesselStruct->FlipBaseline(line);
      Concavity = true;
    }
  } while ((Concavity)); // && (run < 100)

  CalculateConcavityPerBoundaryPoint(TesselStruct);
  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");

  // end
  LOG(0, "RESULT: Added volume in convexization is " << volume << ".");
  return volume;
};


/** Stores triangles to file.
 * \param *out output stream for debugging
 * \param *mol molecule with atoms and bonds
 * \param *TesselStruct Tesselation with boundary triangles
 * \param *filename prefix of filename
 * \param *extraSuffix intermediate suffix
 */
void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
{
        //Info FunctionInfo(__func__);
  PointCloudAdaptor< molecule > cloud(const_cast<molecule *>(mol), mol->name);
  // 4. Store triangles in tecplot file
  if (filename != NULL) {
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      WriteTecplotFile(tecplot, TesselStruct, cloud, -1);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      WriteRaster3dFile(rasterplot, TesselStruct, cloud);
      rasterplot->close();
      delete(rasterplot);
    }
  }
};

/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's
 * \param *&TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *&LCList atoms in LinkedCell_deprecated list
 * \param RADIUS radius of the virtual sphere
 * \param *filename filename prefix for output of vertex data
 * \return true - tesselation successful, false - tesselation failed
 */
bool FindNonConvexBorder(molecule* const mol, Tesselation *&TesselStruct, const LinkedCell_deprecated *&LCList, const double RADIUS, const char *filename = NULL)
{
        //Info FunctionInfo(__func__);
  bool freeLC = false;
  bool status = false;
  CandidateForTesselation *baseline = NULL;
  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
//  bool TesselationFailFlag = false;

  mol->getAtomCount();

  if (TesselStruct == NULL) {
    LOG(1, "Allocating Tesselation struct ...");
    TesselStruct= new Tesselation;
  } else {
    delete(TesselStruct);
    LOG(1, "Re-Allocating Tesselation struct ...");
    TesselStruct = new Tesselation;
  }

  // initialise Linked Cell
  PointCloudAdaptor< molecule > cloud(mol, mol->name);
  if (LCList == NULL) {
    LCList = new LinkedCell_deprecated(cloud, 2.*RADIUS);
    freeLC = true;
  }

  // 1. get starting triangle
  if (!TesselStruct->FindStartingTriangle(RADIUS, LCList)) {
    ELOG(0, "No valid starting triangle found.");
    //performCriticalExit();
  }
  if (filename != NULL) {
    if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
      TesselStruct->Output(filename, cloud);
    }
  }

  // 2. expand from there
  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
    (cerr << "There are " <<  TesselStruct->TrianglesOnBoundary.size() << " triangles and " << TesselStruct->OpenLines.size() << " open lines to scan for candidates." << endl);
    // 2a. print OpenLines without candidates
    LOG(1, "There are the following open lines to scan for a candidates:");
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
      if (Runner->second->pointlist.empty())
        LOG(1, " " << *(Runner->second));

    // 2b. find best candidate for each OpenLine
    const bool TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);
    ASSERT( TesselationFailFlag,
        "FindNonConvexBorder() - at least one open line without candidate exists.");

    // 2c. print OpenLines with candidates again
    LOG(1, "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:");
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
      LOG(1, " " << *(Runner->second));

    // 2d. search for smallest ShortestAngle among all candidates
    double ShortestAngle = 4.*M_PI;
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
      if (Runner->second->ShortestAngle < ShortestAngle) {
        baseline = Runner->second;
        ShortestAngle = baseline->ShortestAngle;
        LOG(1, "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle);
      }
    }
    // 2e. if we found one, add candidate
    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
      OneLoopWithoutSuccessFlag = false;
    else {
      TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
    }

    // 2f. write temporary envelope
    if (filename != NULL) {
      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
        TesselStruct->Output(filename, cloud);
      }
    }
  }
//  // check envelope for consistency
//  status = CheckListOfBaselines(TesselStruct);
//
//  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
//  //->InsertStraddlingPoints(mol, LCList);
//  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
//  class TesselPoint *Runner = NULL;
//    Runner = *iter;
//    LOG(1, "Checking on " << Runner->Name << " ... ");
//    if (!->IsInnerPoint(Runner, LCList)) {
//      LOG(2, Runner->Name << " is outside of envelope, adding via degenerated triangles.");
//      ->AddBoundaryPointByDegeneratedTriangle(Runner, LCList);
//    } else {
//      LOG(2, Runner->Name << " is inside of or on envelope.");
//    }
//  }

//  // Purges surplus triangles.
//  TesselStruct->RemoveDegeneratedTriangles();
//
//  // check envelope for consistency
//  status = CheckListOfBaselines(TesselStruct);

//  cout << "before correction" << endl;

  // store before correction
  StoreTrianglesinFile(mol, TesselStruct, filename, "");

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

  // write final envelope
  CalculateConcavityPerBoundaryPoint(TesselStruct);
//  cout << "after correction" << endl;
  StoreTrianglesinFile(mol, TesselStruct, filename, "");

  if (freeLC)
    delete(LCList);

  return status;
};


/** Finds a hole of sufficient size in \a *mols to embed \a *srcmol into it.
 * \param *out output stream for debugging
 * \param *mols molecules in the domain to embed in between
 * \param *srcmol embedding molecule
 * \return *Vector new center of \a *srcmol for embedding relative to \a this
 */
Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
{
        //Info FunctionInfo(__func__);
  Vector *Center = new Vector;
  Center->Zero();
  // calculate volume/shape of \a *srcmol

  // find embedding holes

  // if more than one, let user choose

  // return embedding center
  return Center;
};