source: src/boundary.cpp@ edb93c

/** \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;
}
;

/** Creates the objects in a VRML file.
 * \param *out output stream for debugging
 * \param *vrmlfile output stream for tecplot data
 * \param *Tess Tesselation structure with constructed triangles
 * \param *mol molecule structure with atom positions
 */
void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
{
  atom *Walker = mol->start;
  bond *Binder = mol->first;
  int i;
  Vector *center = mol->DetermineCenterOfAll(out);
  if (vrmlfile != NULL) {
    //cout << Verbose(1) << "Writing Raster3D file ... ";
    *vrmlfile << "#VRML V2.0 utf8" << endl;
    *vrmlfile << "#Created by molecuilder" << endl;
    *vrmlfile << "#All atoms as spheres" << endl;
    while (Walker->next != mol->end) {
      Walker = Walker->next;
      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
      for (i=0;i<NDIM;i++)
        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    }

    *vrmlfile << "# All bonds as vertices" << endl;
    while (Binder->next != mol->last) {
      Binder = Binder->next;
      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
      for (i=0;i<NDIM;i++)
        *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
      *vrmlfile << "\t0.03\t";
      for (i=0;i<NDIM;i++)
        *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
    }

    *vrmlfile << "# All tesselation triangles" << endl;
    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
      *vrmlfile << "1" << endl << "  "; // 1 is triangle type
      for (i=0;i<3;i++) { // print each node
        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
        *vrmlfile << "\t";
      }
      *vrmlfile << "1. 0. 0." << endl;  // red as colour
      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
    }
  } else {
    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
  }
  delete(center);
};

/** Creates the objects in a raster3d file (renderable with a header.r3d).
 * \param *out output stream for debugging
 * \param *rasterfile output stream for tecplot data
 * \param *Tess Tesselation structure with constructed triangles
 * \param *mol molecule structure with atom positions
 */
void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
{
  atom *Walker = mol->start;
  bond *Binder = mol->first;
  int i;
  Vector *center = mol->DetermineCenterOfAll(out);
  if (rasterfile != NULL) {
    //cout << Verbose(1) << "Writing Raster3D file ... ";
    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
    *rasterfile << "@header.r3d" << endl;
    *rasterfile << "# All atoms as spheres" << endl;
    while (Walker->next != mol->end) {
      Walker = Walker->next;
      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
      for (i=0;i<NDIM;i++)
        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    }

    *rasterfile << "# All bonds as vertices" << endl;
    while (Binder->next != mol->last) {
      Binder = Binder->next;
      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
      for (i=0;i<NDIM;i++)
        *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
      *rasterfile << "\t0.03\t";
      for (i=0;i<NDIM;i++)
        *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
    }

    *rasterfile << "# All tesselation triangles" << endl;
    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
      for (i=0;i<3;i++) {  // print each node
        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
        *rasterfile << "\t";
      }
      *rasterfile << "1. 0. 0." << endl;  // red as colour
      //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
    }
    *rasterfile << "9\n  terminating special property\n";
  } else {
    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
  }
  delete(center);
};

/** This function creates the tecplot file, displaying the tesselation of the hull.
 * \param *out output stream for debugging
 * \param *tecplot output stream for tecplot data
 * \param N arbitrary number to differentiate various zones in the tecplot format
 */
void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
{
  if (tecplot != NULL) {
    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundaryCount << ", E=" << TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    int *LookupList = new int[mol->AtomCount];
    for (int i = 0; i < mol->AtomCount; i++)
      LookupList[i] = -1;

    // print atom coordinates
    *out << Verbose(2) << "The following triangles were created:";
    int Counter = 1;
    TesselPoint *Walker = NULL;
    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
      Walker = target->second->node;
      LookupList[Walker->nr] = Counter++;
      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
    }
    *tecplot << endl;
    // print connectivity
    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    }
    delete[] (LookupList);
    *out << endl;
  }
}


/** Determines the boundary points of a cluster.
 * 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.Projection(&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 Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
{
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

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

  if (TesselStruct != NULL) // free if allocated
    delete(TesselStruct);
  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++)
        TesselStruct->AddPoint(runner->second.second);

  *out << Verbose(2) << "I found " << 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
  TesselStruct->GuessStartingTriangle(out);

  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
  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 (!TesselStruct->InsertStraddlingPoints(out, mol))
    *out << Verbose(1) << "Insertion of straddling points failed!" << endl;

  // 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
  if (!TesselStruct->CorrectConcaveBaselines(out))
    *out << Verbose(1) << "Correction of concave baselines failed!" << endl;

  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << 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());
      write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      write_raster3d_file(out, rasterplot, TesselStruct, mol);
      rasterplot->close();
      delete(rasterplot);
    }
  }

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

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


/** 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->Projection(&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;
}
;

/** 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.);
  Find_convex_border(out, mol, TesselStruct, &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, *Runner = NULL;
  bond *Binder = NULL;

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

  // calculate filler grid in [0,1]^3
  FillerDistance.Init(distance[0], distance[1], distance[2]);
  FillerDistance.InverseMatrixMultiplication(M);
  for(int i=0;i<NDIM;i++)
    N[i] = ceil(1./FillerDistance.x[i]);

  // 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(n[0], n[1], n[2]);
        CurrentPosition.MatrixMultiplication(M);
        // Check whether point is in- or outside
        FillIt = true;
        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
        }

        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.);

          // go through all atoms
          Walker = filler->start;
          while (Walker != filler->end) {
            Walker = Walker->next;
            // copy atom ...
            *Runner = 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)
              Runner->x.MatrixMultiplication(Rotations);
            Runner->x.AddVector(&AtomTranslations);
            Runner->x.AddVector(&FillerTranslations);
            Runner->x.AddVector(&CurrentPosition);
            // insert into Filling
            Filling->AddAtom(Runner);
          }

          // go through all bonds and add as well
          Binder = filler->first;
          while(Binder != filler->last) {
            Binder = Binder->next;
          }
        } else {
          // leave empty
          *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << 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 *filename filename prefix for output of vertex data
 * \para RADIUS radius of the virtual sphere
 */
void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
{
  int N = 0;
  bool freeTess = false;
  bool freeLC = false;
  ofstream *tempstream = NULL;
  char NumberName[255];
  int TriangleFilesWritten = 0;

  *out << Verbose(1) << "Entering search for non convex hull. " << endl;
  if (Tess == NULL) {
    *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
    Tess = new Tesselation;
    freeTess = true;
  }
  LineMap::iterator baseline;
  LineMap::iterator testline;
  *out << Verbose(0) << "Begin of Find_non_convex_border\n";
  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
  bool failflag = false;

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

  Tess->Find_starting_triangle(out, RADIUS, LCList);

  baseline = Tess->LinesOnBoundary.begin();
  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
    if (baseline->second->TrianglesCount == 1) {
      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.
      flag = flag || failflag;
      if (!failflag)
        cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
    } else {
      //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
      if (baseline->second->TrianglesCount != 2)
        cout << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
    }

    N++;
    baseline++;
    if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
      baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
      flag = false;
    }

    // write temporary envelope
    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
      if (DoTecplotOutput) {
        string NameofTempFile(filename);
        NameofTempFile.append(NumberName);
        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
        NameofTempFile.erase(npos, 1);
        NameofTempFile.append(TecplotSuffix);
        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
        tempstream->close();
        tempstream->flush();
        delete(tempstream);
      }

      if (DoRaster3DOutput) {
        string NameofTempFile(filename);
        NameofTempFile.append(NumberName);
        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
        NameofTempFile.erase(npos, 1);
        NameofTempFile.append(Raster3DSuffix);
        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
        write_raster3d_file(out, tempstream, Tess, mol);
//        // include the current position of the virtual sphere in the temporary raster3d file
//        // make the circumsphere's center absolute again
//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
//        helper.AddVector(BaseRay->endpoints[1]->node->node);
//        helper.Scale(0.5);
//        (*it)->OptCenter.AddVector(&helper);
//        Vector *center = mol->DetermineCenterOfAll(out);
//        (*it)->OptCenter.SubtractVector(center);
//        delete(center);
//        // and add to file plus translucency object
//        *tempstream << "# current virtual sphere\n";
//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
//          << "\t" << RADIUS << "\t1 0 0\n";
//        *tempstream << "9\n  terminating special property\n";
        tempstream->close();
        tempstream->flush();
        delete(tempstream);
      }
      if (DoTecplotOutput || DoRaster3DOutput)
        TriangleFilesWritten++;
    }
  }

  // write final envelope
  if (filename != 0) {
    *out << Verbose(1) << "Writing final tecplot file\n";
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      write_tecplot_file(out, tecplot, Tess, mol, -1);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(Raster3DSuffix);
      ofstream *raster = new ofstream(OutputName.c_str());
      write_raster3d_file(out, raster, Tess, mol);
      raster->close();
      delete(raster);
    }
  } else {
    cerr << "ERROR: Could definitively not find all necessary triangles!" << endl;
  }

  cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl;
  int counter = 0;
  for (testline = Tess->LinesOnBoundary.begin(); testline != Tess->LinesOnBoundary.end(); testline++) {
    if (testline->second->TrianglesCount != 2) {
      cout << Verbose(2) << *testline->second << "\t" << testline->second->TrianglesCount << endl;
      counter++;
    }
  }
  if (counter == 0)
    cout << Verbose(2) << "None." << endl;

  if (freeTess)
    delete(Tess);
  if (freeLC)
    delete(LCList);
  *out << Verbose(0) << "End of Find_non_convex_border\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;
};