source: src/boundary.cpp@ aba92d

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


#include "boundary.hpp"

#include<gsl/gsl_poly.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 IsAngstroem whether we have angstroem or atomic units
 * \return NDIM array of the diameters
 */
double *
GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
    bool IsAngstroem)
{
  // get points on boundary of NULL was given as parameter
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = BoundaryPtr;
  if (BoundaryPoints == NULL) {
    BoundaryFreeFlag = true;
    BoundaryPoints = GetBoundaryPoints(out, mol);
  } else {
    *out << Verbose(1) << "Using given boundary points set." << endl;
  }
  // determine biggest "diameter" of cluster for each axis
  Boundaries::iterator Neighbour, OtherNeighbour;
  double *GreatestDiameter = new double[NDIM];
  for (int i = 0; i < NDIM; i++)
    GreatestDiameter[i] = 0.;
  double OldComponent, tmp, w1, w2;
  Vector DistanceVector, OtherVector;
  int component, Othercomponent;
  for (int axis = 0; axis < NDIM; axis++)
    { // regard each projected plane
      //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
      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;
          //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
              != BoundaryPoints[axis].end(); runner++)
            {
              //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
              // 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.CopyVector(&runner->second.second->x);
              DistanceVector.SubtractVector(&Neighbour->second.second->x);
              do
                { // seek for neighbour pair where it flips
                  OldComponent = DistanceVector.x[Othercomponent];
                  Neighbour++;
                  if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                    Neighbour = BoundaryPoints[axis].begin();
                  DistanceVector.CopyVector(&runner->second.second->x);
                  DistanceVector.SubtractVector(&Neighbour->second.second->x);
                  //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                }
              while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
                  DistanceVector.x[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--;
                  //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
                  // now we have found the pair: Neighbour and OtherNeighbour
                  OtherVector.CopyVector(&runner->second.second->x);
                  OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
                  //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
                  //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                  w1 = fabs(OtherVector.x[Othercomponent]);
                  w2 = fabs(DistanceVector.x[Othercomponent]);
                  tmp = fabs((w1 * DistanceVector.x[component] + w2
                      * OtherVector.x[component]) / (w1 + w2));
                  // mark if it has greater diameter
                  //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
                  GreatestDiameter[component] = (GreatestDiameter[component]
                      > tmp) ? GreatestDiameter[component] : tmp;
                } //else
              //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
            }
        }
    }
  *out << Verbose(0) << "RESULT: The biggest diameters are "
      << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
      << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
      : "atomiclength") << "." << endl;

  // 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
 */
Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
{
  atom *Walker = NULL;
  PointMap PointsOnBoundary;
  LineMap LinesOnBoundary;
  TriangleMap TrianglesOnBoundary;
  Vector *MolCenter = mol->DetermineCenterOfAll(out);
  Vector helper;

  *out << Verbose(1) << "Finding all boundary points." << endl;
  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
  BoundariesTestPair BoundaryTestPair;
  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
  double radius, angle;
  // 3a. Go through every axis
  for (int axis = 0; axis < NDIM; axis++) {
    AxisVector.Zero();
    AngleReferenceVector.Zero();
    AngleReferenceNormalVector.Zero();
    AxisVector.x[axis] = 1.;
    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;

    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;

    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    Walker = mol->start;
    while (Walker->next != mol->end) {
      Walker = Walker->next;
      Vector ProjectedVector;
      ProjectedVector.CopyVector(&Walker->x);
      ProjectedVector.SubtractVector(MolCenter);
      ProjectedVector.ProjectOntoPlane(&AxisVector);

      // correct for negative side
      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

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

          SideB.CopyVector(&right->second.second->x);
          SideB.SubtractVector(MolCenter);
          SideB.ProjectOntoPlane(&AxisVector);
          //          *out << "SideB: ";
          //          SideB.Output(out);
          //          *out << endl;

          SideC.CopyVector(&left->second.second->x);
          SideC.SubtractVector(&right->second.second->x);
          SideC.ProjectOntoPlane(&AxisVector);
          //          *out << "SideC: ";
          //          SideC.Output(out);
          //          *out << endl;

          SideH.CopyVector(&runner->second.second->x);
          SideH.SubtractVector(MolCenter);
          SideH.ProjectOntoPlane(&AxisVector);
          //          *out << "SideH: ";
          //          SideH.Output(out);
          //          *out << endl;

          // calculate each length
          double a = SideA.Norm();
          //double b = SideB.Norm();
          //double c = SideC.Norm();
          double h = SideH.Norm();
          // calculate the angles
          double alpha = SideA.Angle(&SideH);
          double beta = SideA.Angle(&SideC);
          double gamma = SideB.Angle(&SideH);
          double delta = SideC.Angle(&SideH);
          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
          //*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;
          *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;
          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
            // throw out point
            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
            BoundaryPoints[axis].erase(runner);
            flag = true;
          }
        }
      }
    } while (flag);
  }
  delete(MolCenter);
  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 *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList atoms in LinkedCell 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(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename)
{
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

  cout << Verbose(1) << "Begin of FindConvexBorder" << endl;

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

  // 1. Find all points on the boundary
  if (BoundaryPoints == NULL) {
      BoundaryFreeFlag = true;
      BoundaryPoints = GetBoundaryPoints(out, mol);
  } else {
      *out << Verbose(1) << "Using given boundary points set." << endl;
  }

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

  // 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 (!mol->TesselStruct->AddBoundaryPoint(runner->second.second, 0))
          *out << Verbose(3) << "WARNING: Point " << *(runner->second.second) << " is already present!" << endl;

  *out << Verbose(2) << "I found " << mol->TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
  // now we have the whole set of edge points in the BoundaryList

  // listing for debugging
  //  *out << Verbose(1) << "Listing PointsOnBoundary:";
  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
  //    *out << " " << *runner->second;
  //  }
  //  *out << endl;

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

  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
  mol->TesselStruct->TesselateOnBoundary(out, mol);

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

  *out << Verbose(2) << "I created " << mol->TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << mol->TesselStruct->LinesOnBoundary.size() << " lines and " << mol->TesselStruct->PointsOnBoundary.size() << " points." << endl;

  // 4. Store triangles in tecplot file
  if (filename != NULL) {
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append("_intermed");
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append("_intermed");
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
      rasterplot->close();
      delete(rasterplot);
    }
  }

  // 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;
  class BoundaryLineSet *line = NULL;
  do {
    AllConvex = true;
    for (LineMap::iterator LineRunner = mol->TesselStruct->LinesOnBoundary.begin(); LineRunner != mol->TesselStruct->LinesOnBoundary.end(); LineRunner++) {
      line = LineRunner->second;
      *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
      if (!line->CheckConvexityCriterion(out)) {
        *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;

        // flip the line
        if (mol->TesselStruct->PickFarthestofTwoBaselines(out, line) == 0.)
          *out << Verbose(1) << "ERROR: Correction of concave baselines failed!" << endl;
        else {
          mol->TesselStruct->FlipBaseline(out, line);
          *out << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
        }
      }
    }
  } 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 (!mol->TesselStruct->CorrectConcaveTesselPoints(out))
//    *out << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;

  *out << Verbose(2) << "I created " << mol->TesselStruct->TrianglesOnBoundary.size() << " triangles with " << mol->TesselStruct->LinesOnBoundary.size() << " lines and " << mol->TesselStruct->PointsOnBoundary.size() << " points." << endl;

  // 4. Store triangles in tecplot file
  if (filename != NULL) {
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
      rasterplot->close();
      delete(rasterplot);
    }
  }


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

  cout << Verbose(1) << "End of FindConvexBorder" << endl;
};

/** 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(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename)
{
  int i=0;
  char number[MAXSTRINGSIZE];

  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    *out << Verbose(2) << "ERROR: TesselStruct is empty." << endl;
    return false;
  }

  PointMap::iterator PointRunner;
  while (!TesselStruct->PointsOnBoundary.empty()) {
    *out << Verbose(2) << "Remaining points are: ";
    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
      *out << *(PointSprinter->second) << "\t";
      *out << endl;

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

    // store envelope
    sprintf(number, "-%04d", i++);
    StoreTrianglesinFile(out, mol, 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 - VolumeOfConvexEnvelope() - 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(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename)
{
  double volume = 0;
  class BoundaryPointSet *point = NULL;
  class BoundaryLineSet *line = NULL;
  bool Concavity;
  char dummy[MAXSTRINGSIZE];
  PointMap::iterator PointRunner, PointAdvance;
  LineMap::iterator LineRunner, LineAdvance;
  TriangleMap::iterator TriangleRunner, TriangleAdvance;

  *out << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;

  // check whether there is something to work on
  if (TesselStruct == NULL) {
    *out << Verbose(1) << "ERROR: TesselStruct is empty!" << endl;
    return volume;
  }

  // First step: RemovePointFromTesselatedSurface
  int run = 0;
  double tmp;
  do {
    Concavity = false;
    sprintf(dummy, "-first-%d", run);
    //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
    StoreTrianglesinFile(out, mol, filename, dummy);

    PointRunner = TesselStruct->PointsOnBoundary.begin();
    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
    while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
      PointAdvance++;
      point = PointRunner->second;
      *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
        line = LineRunner->second;
        *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
        if (!line->CheckConvexityCriterion(out)) {
          // remove the point if needed
          *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
          volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
          sprintf(dummy, "-first-%d", ++run);
          StoreTrianglesinFile(out, mol, filename, dummy);
          Concavity = true;
          break;
        }
      }
      PointRunner = PointAdvance;
    }

    sprintf(dummy, "-second-%d", run);
    //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
    StoreTrianglesinFile(out, mol, filename, dummy);

    // second step: PickFarthestofTwoBaselines
    LineRunner = TesselStruct->LinesOnBoundary.begin();
    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
      LineAdvance++;
      line = LineRunner->second;
      *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
      // take highest of both lines
      if (TesselStruct->IsConvexRectangle(out, line) == NULL) {
        tmp = TesselStruct->PickFarthestofTwoBaselines(out, line);
        volume += tmp;
        if (tmp != 0) {
          mol->TesselStruct->FlipBaseline(out, line);
          Concavity = true;
        }
      }
      LineRunner = LineAdvance;
    }
    run++;
  } while (Concavity);
  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
  //StoreTrianglesinFile(out, mol, filename, "-third");

  // third step: IsConvexRectangle
//  LineRunner = TesselStruct->LinesOnBoundary.begin();
//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
//    LineAdvance++;
//    line = LineRunner->second;
//    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
//    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
//      //*out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
//    if (!line->CheckConvexityCriterion(out)) {
//      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
//
//      // take highest of both lines
//      point = TesselStruct->IsConvexRectangle(out, line);
//      if (point != NULL)
//        volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
//    }
//    LineRunner = LineAdvance;
//  }

  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
  StoreTrianglesinFile(out, mol, filename, "");

  // end
  *out << Verbose(1) << "Volume is " << volume << "." << endl;
  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
  return volume;
};


/** Determines the volume of a cluster.
 * Determines first the convex envelope, then tesselates it and calculates its volume.
 * \param *out output stream for debugging
 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *configuration needed for path to store convex envelope file
 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
 */
double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
{
  bool IsAngstroem = configuration->GetIsAngstroem();
  double volume = 0.;
  double PyramidVolume = 0.;
  double G, h;
  Vector x, y;
  double a, b, c;

  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
  *out << Verbose(1)
      << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
      << endl;
  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
    { // go through every triangle, calculate volume of its pyramid with CoG as peak
      x.CopyVector(runner->second->endpoints[0]->node->node);
      x.SubtractVector(runner->second->endpoints[1]->node->node);
      y.CopyVector(runner->second->endpoints[0]->node->node);
      y.SubtractVector(runner->second->endpoints[2]->node->node);
      a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
      b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
      c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
      G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
      h = x.Norm(); // distance of CoG to triangle
      PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
      *out << Verbose(2) << "Area of triangle is " << G << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
          << h << " and the volume is " << PyramidVolume << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
      volume += PyramidVolume;
    }
  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
      << endl;

  return volume;
};

/** Stores triangles to file.
 * \param *out output stream for debugging
 * \param *mol molecule with atoms and bonds
 * \param *filename prefix of filename
 * \param *extraSuffix intermediate suffix
 */
void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
{
  // 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(out, tecplot, mol->TesselStruct, mol, 0);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
      rasterplot->close();
      delete(rasterplot);
    }
  }
};

/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
 * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
 * \param *out output stream for debugging
 * \param *configuration needed for path to store convex envelope file
 * \param *mol molecule structure representing the cluster
 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
 * \param celldensity desired average density in final cell
 */
void
PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
    double ClusterVolume, double celldensity)
{
  // transform to PAS
  mol->PrincipalAxisSystem(out, true);

  // some preparations beforehand
  bool IsAngstroem = configuration->GetIsAngstroem();
  Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
  class Tesselation *TesselStruct = NULL;
  LinkedCell LCList(mol, 10.);
  FindConvexBorder(out, mol, &LCList, NULL);
  double clustervolume;
  if (ClusterVolume == 0)
    clustervolume = VolumeOfConvexEnvelope(out, TesselStruct, configuration);
  else
    clustervolume = ClusterVolume;
  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
  Vector BoxLengths;
  int repetition[NDIM] =
    { 1, 1, 1 };
  int TotalNoClusters = 1;
  for (int i = 0; i < NDIM; i++)
    TotalNoClusters *= repetition[i];

  // sum up the atomic masses
  double totalmass = 0.;
  atom *Walker = mol->start;
  while (Walker->next != mol->end)
    {
      Walker = Walker->next;
      totalmass += Walker->type->mass;
    }
  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
      << totalmass << " atomicmassunit." << endl;

  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
      << totalmass / clustervolume << " atomicmassunit/"
      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;

  // solve cubic polynomial
  *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
      << endl;
  double cellvolume;
  if (IsAngstroem)
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
        / clustervolume)) / (celldensity - 1);
  else
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
        / clustervolume)) / (celldensity - 1);
  *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
      << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
      : "atomiclength") << "^3." << endl;

  double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
      * GreatestDiameter[1] * GreatestDiameter[2]);
  *out << Verbose(1)
      << "Minimum volume of the convex envelope contained in a rectangular box is "
      << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
      : "atomiclength") << "^3." << endl;
  if (minimumvolume > cellvolume)
    {
      cerr << Verbose(0)
          << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
          << endl;
      cout << Verbose(0)
          << "Setting Box dimensions to minimum possible, the greatest diameters."
          << endl;
      for (int i = 0; i < NDIM; i++)
        BoxLengths.x[i] = GreatestDiameter[i];
      mol->CenterEdge(out, &BoxLengths);
    }
  else
    {
      BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
          * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
      BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
          * GreatestDiameter[1] + repetition[0] * repetition[2]
          * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
          * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
      BoxLengths.x[2] = minimumvolume - cellvolume;
      double x0 = 0., x1 = 0., x2 = 0.;
      if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
          BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
        *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
            << " ." << endl;
      else
        {
          *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
              << " and " << x1 << " and " << x2 << " ." << endl;
          x0 = x2; // sorted in ascending order
        }

      cellvolume = 1;
      for (int i = 0; i < NDIM; i++)
        {
          BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
          cellvolume *= BoxLengths.x[i];
        }

      // set new box dimensions
      *out << Verbose(0) << "Translating to box with these boundaries." << endl;
      mol->SetBoxDimension(&BoxLengths);
      mol->CenterInBox((ofstream *) &cout);
    }
  // update Box of atoms by boundary
  mol->SetBoxDimension(&BoxLengths);
  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
      << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
      << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
}
;


/** Fills the empty space of the simulation box with water/
 * \param *out output stream for debugging
 * \param *List list of molecules already present in box
 * \param *TesselStruct contains tesselated surface
 * \param *filler molecule which the box is to be filled with
 * \param configuration contains box dimensions
 * \param distance[NDIM] distance between filling molecules in each direction
 * \param RandAtomDisplacement maximum distance for random displacement per atom
 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
 * \param DoRandomRotation true - do random rotiations, false - don't
 * \return *mol pointer to new molecule with filled atoms
 */
molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
{
  molecule *Filling = new molecule(filler->elemente);
  Vector CurrentPosition;
  int N[NDIM];
  int n[NDIM];
  double *M =  filler->ReturnFullMatrixforSymmetric(filler->cell_size);
  double Rotations[NDIM*NDIM];
  Vector AtomTranslations;
  Vector FillerTranslations;
  Vector FillerDistance;
  double FillIt = false;
  atom *Walker = NULL;
  bond *Binder = NULL;
  int i;
  LinkedCell *LCList[List->ListOfMolecules.size()];

  *out << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;

  i=0;
  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    *out << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
    LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
    if ((*ListRunner)->TesselStruct == NULL) {
      *out << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], 5., NULL);
    }
    i++;
  }

  // Center filler at origin
  filler->CenterOrigin(out);
  filler->Center.Zero();

  filler->CountAtoms(out);
  atom * CopyAtoms[filler->AtomCount];
  int nr = 0;

  // calculate filler grid in [0,1]^3
  FillerDistance.Init(distance[0], distance[1], distance[2]);
  FillerDistance.InverseMatrixMultiplication(M);
  *out << Verbose(1) << "INFO: Grid steps are ";
  for(int i=0;i<NDIM;i++) {
    N[i] = (int) ceil(1./FillerDistance.x[i]);
    *out << N[i];
    if (i != NDIM-1)
      *out<< ", ";
    else
      *out << "." << endl;
  }

  // go over [0,1]^3 filler grid
  for (n[0] = 0; n[0] < N[0]; n[0]++)
    for (n[1] = 0; n[1] < N[1]; n[1]++)
      for (n[2] = 0; n[2] < N[2]; n[2]++) {
        // calculate position of current grid vector in untransformed box
        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
        CurrentPosition.MatrixMultiplication(M);
        *out << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
        // Check whether point is in- or outside
        FillIt = true;
        i=0;
        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
          // get linked cell list
          if ((*ListRunner)->TesselStruct == NULL) {
            *out << Verbose(1) << "ERROR: TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
            FillIt = false;
          } else
            FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInnerPoint(out, CurrentPosition, LCList[i++]));
        }

        if (FillIt) {
          // fill in Filler
          *out << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;

          // create molecule random translation vector ...
          for (int i=0;i<NDIM;i++)
            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
          *out << Verbose(3) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;

          // go through all atoms
          nr=0;
          Walker = filler->start;
          while (Walker->next != filler->end) {
            Walker = Walker->next;
            // copy atom ...
            CopyAtoms[Walker->nr] = new atom(Walker);

            // create atomic random translation vector ...
            for (int i=0;i<NDIM;i++)
              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);

            // ... and rotation matrix
            if (DoRandomRotation) {
              double phi[NDIM];
              for (int i=0;i<NDIM;i++) {
                phi[i] = rand()/(RAND_MAX/(2.*M_PI));
              }

              Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
              Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
              Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
              Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
              Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
              Rotations[7] =               sin(phi[1])                                                ;
              Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
              Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
              Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
            }

            // ... and put at new position
            if (DoRandomRotation)
              CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
            CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
            CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
            CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);

            // insert into Filling

            // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
            *out << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
            Filling->AddAtom(CopyAtoms[Walker->nr]);
          }

          // go through all bonds and add as well
          Binder = filler->first;
          while(Binder->next != filler->last) {
            Binder = Binder->next;
            *out << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
          }
        } else {
          // leave empty
          *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
        }
      }
  *out << Verbose(0) << "End of FillBoxWithMolecule" << endl;

  return Filling;
};


/** 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 *Tess Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList atoms in LinkedCell list
 * \param RADIUS radius of the virtual sphere
 * \param *filename filename prefix for output of vertex data
 */
void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const double RADIUS, const char *filename = NULL)
{
  bool freeLC = false;

  *out << Verbose(1) << "Entering search for non convex hull. " << endl;
  if (mol->TesselStruct == NULL) {
    *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
    mol->TesselStruct = new Tesselation;
  } else {
    delete(mol->TesselStruct);
    *out << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
    mol->TesselStruct = new Tesselation;
  }
  LineMap::iterator baseline;
  LineMap::iterator testline;
  *out << Verbose(0) << "Begin of FindNonConvexBorder\n";
  bool OneLoopWithoutSuccessFlag = false;  // marks whether we went once through all baselines without finding any without two triangles
  bool TesselationFailFlag = false;

  // initialise Linked Cell
  if (LCList == NULL) {
    LCList = new LinkedCell(mol, 2.*RADIUS);
    freeLC = true;
  }

  // 1. get starting triangle
  mol->TesselStruct->FindStartingTriangle(out, RADIUS, LCList);

  // 2. expand from there
  baseline = mol->TesselStruct->LinesOnBoundary.begin();
  baseline++; // skip first line
  while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
    if (baseline->second->triangles.size() == 1) {
      // 3. find next triangle
      TesselationFailFlag = mol->TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
      OneLoopWithoutSuccessFlag = OneLoopWithoutSuccessFlag || TesselationFailFlag;
      if (!TesselationFailFlag)
        cerr << "WARNING: FindNextSuitableTriangle failed." << endl;

      // write temporary envelope
      if (filename != NULL) {
        if ((DoSingleStepOutput && ((mol->TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
          mol->TesselStruct->Output(out, filename, mol);
        }
      }
      baseline = mol->TesselStruct->LinesOnBoundary.end();
      *out << Verbose(2) << "Baseline set to end." << endl;
    } else {
      //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
      if (baseline->second->triangles.size() != 2)
        *out << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
    }

    if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
      baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
      OneLoopWithoutSuccessFlag = false;
    }
    baseline++;
  }
  // check envelope for consistency
  CheckListOfBaselines(out, mol->TesselStruct);

  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
  mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
//  mol->GoToFirst();
//  class TesselPoint *Runner = NULL;
//  while (!mol->IsEnd()) {
//    Runner = mol->GetPoint();
//    *out << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
//    if (!mol->TesselStruct->IsInnerPoint(out, Runner, LCList)) {
//      *out << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
//      mol->TesselStruct->AddBoundaryPointByDegeneratedTriangle(out, Runner, LCList);
//    } else {
//      *out << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
//    }
//    mol->GoToNext();
//  }

  // Purges surplus triangles.
  mol->TesselStruct->RemoveDegeneratedTriangles();

  // check envelope for consistency
  CheckListOfBaselines(out, mol->TesselStruct);

  // write final envelope
  CalculateConcavityPerBoundaryPoint(out, mol->TesselStruct);
  StoreTrianglesinFile(out, mol, filename, "");

  if (freeLC)
    delete(LCList);
  *out << Verbose(0) << "End of FindNonConvexBorder\n";
};


/** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.
 * \param *out output stream for debugging
 * \param *srcmol molecule to embed into
 * \return *Vector new center of \a *srcmol for embedding relative to \a this
 */
Vector* molecule::FindEmbeddingHole(ofstream *out, molecule *srcmol)
{
  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;
};