Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselationhelpers.cpp

    rf66195 r658efb  
    66 */
    77
    8 #include <fstream>
    9 
    10 #include "linkedcell.hpp"
    11 #include "tesselation.hpp"
    128#include "tesselationhelpers.hpp"
    13 #include "vector.hpp"
    14 #include "verbose.hpp"
    15 
    16 double DetGet(gsl_matrix *A, int inPlace) {
     9
     10double det_get(gsl_matrix *A, int inPlace) {
    1711  /*
    1812  inPlace = 1 => A is replaced with the LU decomposed copy.
    … …  
    4236};
    4337
    44 void GetSphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
     38void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
    4539{
    4640  gsl_matrix *A = gsl_matrix_calloc(3,3);
    … …  
    5246    gsl_matrix_set(A, i, 2, c.x[i]);
    5347  }
    54   m11 = DetGet(A, 1);
     48  m11 = det_get(A, 1);
    5549
    5650  for(int i=0;i<3;i++) {
    … …  
    5953    gsl_matrix_set(A, i, 2, c.x[i]);
    6054  }
    61   m12 = DetGet(A, 1);
     55  m12 = det_get(A, 1);
    6256
    6357  for(int i=0;i<3;i++) {
    … …  
    6660    gsl_matrix_set(A, i, 2, c.x[i]);
    6761  }
    68   m13 = DetGet(A, 1);
     62  m13 = det_get(A, 1);
    6963
    7064  for(int i=0;i<3;i++) {
    … …  
    7367    gsl_matrix_set(A, i, 2, b.x[i]);
    7468  }
    75   m14 = DetGet(A, 1);
     69  m14 = det_get(A, 1);
    7670
    7771  if (fabs(m11) < MYEPSILON)
    … …  
    107101 * @param Umkreisradius double radius of circumscribing circle
    108102 */
    109 void GetCenterOfSphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
     103void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
    110104    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
    111105{
    … …  
    113107  double Restradius;
    114108  Vector OtherCenter;
    115   cout << Verbose(3) << "Begin of GetCenterOfSphere.\n";
     109  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
    116110  Center->Zero();
    117111  helper.CopyVector(&a);
    … …  
    159153  Center->AddVector(&TempNormal);
    160154  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
    161   GetSphere(&OtherCenter, a, b, c, RADIUS);
     155  get_sphere(&OtherCenter, a, b, c, RADIUS);
    162156  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
    163   cout << Verbose(3) << "End of GetCenterOfSphere.\n";
     157  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
    164158};
    165159
    … …  
    313307    gsl_multimin_fminimizer *s = NULL;
    314308    gsl_vector *ss, *x;
    315     gsl_multimin_function minexFunction;
     309    gsl_multimin_function minex_func;
    316310
    317311    size_t iter = 0;
    … …  
    330324
    331325    /* Initialize method and iterate */
    332     minexFunction.n = NDIM;
    333     minexFunction.f = &MinIntersectDistance;
    334     minexFunction.params = (void *)&par;
     326    minex_func.n = NDIM;
     327    minex_func.f = &MinIntersectDistance;
     328    minex_func.params = (void *)&par;
    335329
    336330    s = gsl_multimin_fminimizer_alloc(T, NDIM);
    337     gsl_multimin_fminimizer_set(s, &minexFunction, x, ss);
     331    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
    338332
    339333    do {
    … …  
    392386
    393387  return result;
    394 };
    395 
    396 /** Gets the angle between a point and a reference relative to the provided center.
    397  * We have two shanks point and reference between which the angle is calculated
    398  * and by scalar product with OrthogonalVector we decide the interval.
    399  * @param point to calculate the angle for
    400  * @param reference to which to calculate the angle
    401  * @param OrthogonalVector points in direction of [pi,2pi] interval
    402  *
    403  * @return angle between point and reference
    404  */
    405 double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
    406 {
    407   if (reference.IsZero())
    408     return M_PI;
    409 
    410   // calculate both angles and correct with in-plane vector
    411   if (point.IsZero())
    412     return M_PI;
    413   double phi = point.Angle(&reference);
    414   if (OrthogonalVector.ScalarProduct(&point) > 0) {
    415     phi = 2.*M_PI - phi;
    416   }
    417 
    418   cout << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
    419 
    420   return phi;
    421388}
    422389
    423 
    424 /** Calculates the volume of a general tetraeder.
    425  * \param *a first vector
    426  * \param *a first vector
    427  * \param *a first vector
    428  * \param *a first vector
    429  * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
    430  */
    431 double CalculateVolumeofGeneralTetraeder(Vector *a, Vector *b, Vector *c, Vector *d)
    432 {
    433   Vector Point, TetraederVector[3];
    434   double volume;
    435 
    436   TetraederVector[0].CopyVector(a);
    437   TetraederVector[1].CopyVector(b);
    438   TetraederVector[2].CopyVector(c);
    439   for (int j=0;j<3;j++)
    440     TetraederVector[j].SubtractVector(d);
    441   Point.CopyVector(&TetraederVector[0]);
    442   Point.VectorProduct(&TetraederVector[1]);
    443   volume = 1./6. * fabs(Point.ScalarProduct(&TetraederVector[2]));
    444   return volume;
    445 };
    446 
    447 
    448 /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
    449  * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
    450  * make it bigger (i.e. closing one (the baseline) and opening two new ones).
    451  * \param TPS[3] nodes of the triangle
    452  * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
    453  */
    454 bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3])
    455 {
    456   bool result = false;
    457   int counter = 0;
    458 
    459   // check all three points
    460   for (int i=0;i<3;i++)
    461     for (int j=i+1; j<3; j++) {
    462       if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
    463         LineMap::iterator FindLine;
    464         pair<LineMap::iterator,LineMap::iterator> FindPair;
    465         FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
    466         for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
    467           // If there is a line with less than two attached triangles, we don't need a new line.
    468           if (FindLine->second->triangles.size() < 2) {
    469             counter++;
    470             break;  // increase counter only once per edge
    471           }
    472         }
    473       } else { // no line
    474         cout << Verbose(1) << "The line between " << *nodes[i] << " and " << *nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
    475         result = true;
    476       }
    477     }
    478   if ((!result) && (counter > 1)) {
    479     cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    480     result = true;
    481   }
    482   return result;
    483 };
    484 
    485 
    486 /** Sort function for the candidate list.
    487  */
    488 bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
    489 {
    490   Vector BaseLineVector, OrthogonalVector, helper;
    491   if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    492     cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    493     //return false;
    494     exit(1);
    495   }
    496   // create baseline vector
    497   BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
    498   BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    499   BaseLineVector.Normalize();
    500 
    501   // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
    502   helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
    503   helper.SubtractVector(candidate1->point->node);
    504   OrthogonalVector.CopyVector(&helper);
    505   helper.VectorProduct(&BaseLineVector);
    506   OrthogonalVector.SubtractVector(&helper);
    507   OrthogonalVector.Normalize();
    508 
    509   // calculate both angles and correct with in-plane vector
    510   helper.CopyVector(candidate1->point->node);
    511   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    512   double phi = BaseLineVector.Angle(&helper);
    513   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    514     phi = 2.*M_PI - phi;
    515   }
    516   helper.CopyVector(candidate2->point->node);
    517   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    518   double psi = BaseLineVector.Angle(&helper);
    519   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    520     psi = 2.*M_PI - psi;
    521   }
    522 
    523   cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
    524   cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
    525 
    526   // return comparison
    527   return phi < psi;
    528 };
    529 
    530 /**
    531  * Finds the point which is second closest to the provided one.
    532  *
    533  * @param Point to which to find the second closest other point
    534  * @param linked cell structure
    535  *
    536  * @return point which is second closest to the provided one
    537  */
    538 TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
    539 {
    540   LinkedNodes *List = NULL;
    541   TesselPoint* closestPoint = NULL;
    542   TesselPoint* secondClosestPoint = NULL;
    543   double distance = 1e16;
    544   double secondDistance = 1e16;
    545   Vector helper;
    546   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
    547 
    548   LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
    549   for(int i=0;i<NDIM;i++) // store indices of this cell
    550     N[i] = LC->n[i];
    551   cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    552 
    553   LC->GetNeighbourBounds(Nlower, Nupper);
    554   //cout << endl;
    555   for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    556     for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    557       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    558         List = LC->GetCurrentCell();
    559         //cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    560         if (List != NULL) {
    561           for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    562             helper.CopyVector(Point);
    563             helper.SubtractVector((*Runner)->node);
    564             double currentNorm = helper. Norm();
    565             if (currentNorm < distance) {
    566               // remember second point
    567               secondDistance = distance;
    568               secondClosestPoint = closestPoint;
    569               // mark down new closest point
    570               distance = currentNorm;
    571               closestPoint = (*Runner);
    572               //cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
    573             }
    574           }
    575         } else {
    576           cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
    577             << LC->n[2] << " is invalid!" << endl;
    578         }
    579       }
    580 
    581   return secondClosestPoint;
    582 };
    583 
    584 /**
    585  * Finds the point which is closest to the provided one.
    586  *
    587  * @param Point to which to find the closest other point
    588  * @param SecondPoint the second closest other point on return, NULL if none found
    589  * @param linked cell structure
    590  *
    591  * @return point which is closest to the provided one, NULL if none found
    592  */
    593 TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
    594 {
    595   LinkedNodes *List = NULL;
    596   TesselPoint* closestPoint = NULL;
    597   SecondPoint = NULL;
    598   double distance = 1e16;
    599   double secondDistance = 1e16;
    600   Vector helper;
    601   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
    602 
    603   LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
    604   for(int i=0;i<NDIM;i++) // store indices of this cell
    605     N[i] = LC->n[i];
    606   cout << Verbose(3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    607 
    608   LC->GetNeighbourBounds(Nlower, Nupper);
    609   //cout << endl;
    610   for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    611     for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    612       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    613         List = LC->GetCurrentCell();
    614         //cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    615         if (List != NULL) {
    616           for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    617             helper.CopyVector(Point);
    618             helper.SubtractVector((*Runner)->node);
    619             double currentNorm = helper. Norm();
    620             if (currentNorm < distance) {
    621               secondDistance = distance;
    622               SecondPoint = closestPoint;
    623               distance = currentNorm;
    624               closestPoint = (*Runner);
    625               //cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    626             } else if (currentNorm < secondDistance) {
    627               secondDistance = currentNorm;
    628               SecondPoint = (*Runner);
    629               //cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
    630             }
    631           }
    632         } else {
    633           cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
    634             << LC->n[2] << " is invalid!" << endl;
    635         }
    636       }
    637 
    638   return closestPoint;
    639 };
    640 
    641 /** Returns the closest point on \a *Base with respect to \a *OtherBase.
    642  * \param *out output stream for debugging
    643  * \param *Base reference line
    644  * \param *OtherBase other base line
    645  * \return Vector on reference line that has closest distance
    646  */
    647 Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
    648 {
    649   // construct the plane of the two baselines (i.e. take both their directional vectors)
    650   Vector Normal;
    651   Vector Baseline, OtherBaseline;
    652   Baseline.CopyVector(Base->endpoints[1]->node->node);
    653   Baseline.SubtractVector(Base->endpoints[0]->node->node);
    654   OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
    655   OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
    656   Normal.CopyVector(&Baseline);
    657   Normal.VectorProduct(&OtherBaseline);
    658   Normal.Normalize();
    659   *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
    660 
    661   // project one offset point of OtherBase onto this plane (and add plane offset vector)
    662   Vector NewOffset;
    663   NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
    664   NewOffset.SubtractVector(Base->endpoints[0]->node->node);
    665   NewOffset.ProjectOntoPlane(&Normal);
    666   NewOffset.AddVector(Base->endpoints[0]->node->node);
    667   Vector NewDirection;
    668   NewDirection.CopyVector(&NewOffset);
    669   NewDirection.AddVector(&OtherBaseline);
    670 
    671   // calculate the intersection between this projected baseline and Base
    672   Vector *Intersection = new Vector;
    673   Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
    674   Normal.CopyVector(Intersection);
    675   Normal.SubtractVector(Base->endpoints[0]->node->node);
    676   *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
    677 
    678   return Intersection;
    679 };
    680 
    681 
    682 /** Creates the objects in a VRML file.
    683  * \param *out output stream for debugging
    684  * \param *vrmlfile output stream for tecplot data
    685  * \param *Tess Tesselation structure with constructed triangles
    686  * \param *mol molecule structure with atom positions
    687  */
    688 void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud)
    689 {
    690   TesselPoint *Walker = NULL;
    691   int i;
    692   Vector *center = cloud->GetCenter(out);
    693   if (vrmlfile != NULL) {
    694     //cout << Verbose(1) << "Writing Raster3D file ... ";
    695     *vrmlfile << "#VRML V2.0 utf8" << endl;
    696     *vrmlfile << "#Created by molecuilder" << endl;
    697     *vrmlfile << "#All atoms as spheres" << endl;
    698     cloud->GoToFirst();
    699     while (!cloud->IsEnd()) {
    700       Walker = cloud->GetPoint();
    701       *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
    702       for (i=0;i<NDIM;i++)
    703         *vrmlfile << Walker->node->x[i]-center->x[i] << " ";
    704       *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    705       cloud->GoToNext();
    706     }
    707 
    708     *vrmlfile << "# All tesselation triangles" << endl;
    709     for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
    710       *vrmlfile << "1" << endl << "  "; // 1 is triangle type
    711       for (i=0;i<3;i++) { // print each node
    712         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    713           *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
    714         *vrmlfile << "\t";
    715       }
    716       *vrmlfile << "1. 0. 0." << endl;  // red as colour
    717       *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
    718     }
    719   } else {
    720     cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
    721   }
    722   delete(center);
    723 };
    724 
    725 /** Writes additionally the current sphere (i.e. the last triangle to file).
    726  * \param *out output stream for debugging
    727  * \param *rasterfile output stream for tecplot data
    728  * \param *Tess Tesselation structure with constructed triangles
    729  * \param *mol molecule structure with atom positions
    730  */
    731 void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
    732 {
    733   Vector helper;
    734   // include the current position of the virtual sphere in the temporary raster3d file
    735   Vector *center = cloud->GetCenter(out);
    736   // make the circumsphere's center absolute again
    737   helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
    738   helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
    739   helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
    740   helper.Scale(1./3.);
    741   helper.SubtractVector(center);
    742   // and add to file plus translucency object
    743   *rasterfile << "# current virtual sphere\n";
    744   *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
    745   *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
    746   *rasterfile << "9\n  terminating special property\n";
    747   delete(center);
    748 };
    749 
    750 /** Creates the objects in a raster3d file (renderable with a header.r3d).
    751  * \param *out output stream for debugging
    752  * \param *rasterfile output stream for tecplot data
    753  * \param *Tess Tesselation structure with constructed triangles
    754  * \param *mol molecule structure with atom positions
    755  */
    756 void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
    757 {
    758   TesselPoint *Walker = NULL;
    759   int i;
    760   Vector *center = cloud->GetCenter(out);
    761   if (rasterfile != NULL) {
    762     //cout << Verbose(1) << "Writing Raster3D file ... ";
    763     *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
    764     *rasterfile << "@header.r3d" << endl;
    765     *rasterfile << "# All atoms as spheres" << endl;
    766     cloud->GoToFirst();
    767     while (!cloud->IsEnd()) {
    768       Walker = cloud->GetPoint();
    769       *rasterfile << "2" << endl << "  ";  // 2 is sphere type
    770       for (i=0;i<NDIM;i++)
    771         *rasterfile << Walker->node->x[i]-center->x[i] << " ";
    772       *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    773       cloud->GoToNext();
    774     }
    775 
    776     *rasterfile << "# All tesselation triangles" << endl;
    777     *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";
    778     for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
    779       *rasterfile << "1" << endl << "  ";  // 1 is triangle type
    780       for (i=0;i<3;i++) {  // print each node
    781         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    782           *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
    783         *rasterfile << "\t";
    784       }
    785       *rasterfile << "1. 0. 0." << endl;  // red as colour
    786       //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
    787     }
    788     *rasterfile << "9\n#  terminating special property\n";
    789   } else {
    790     cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
    791   }
    792   IncludeSphereinRaster3D(out, rasterfile, Tess, cloud);
    793   delete(center);
    794 };
    795 
    796 /** This function creates the tecplot file, displaying the tesselation of the hull.
    797  * \param *out output stream for debugging
    798  * \param *tecplot output stream for tecplot data
    799  * \param N arbitrary number to differentiate various zones in the tecplot format
    800  */
    801 void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N)
    802 {
    803   if ((tecplot != NULL) && (TesselStruct != NULL)) {
    804     // write header
    805     *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    806     *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
    807     *tecplot << "ZONE T=\"" << N << "-";
    808     for (int i=0;i<3;i++)
    809       *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
    810     *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    811     int i=0;
    812     for (cloud->GoToFirst(); !cloud->IsEnd(); cloud->GoToNext(), i++);
    813     int *LookupList = new int[i];
    814     for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)
    815       LookupList[i] = -1;
    816 
    817     // print atom coordinates
    818     *out << Verbose(2) << "The following triangles were created:";
    819     int Counter = 1;
    820     TesselPoint *Walker = NULL;
    821     for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
    822       Walker = target->second->node;
    823       LookupList[Walker->nr] = Counter++;
    824       *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
    825     }
    826     *tecplot << endl;
    827     // print connectivity
    828     for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
    829       *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
    830       *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    831     }
    832     delete[] (LookupList);
    833     *out << endl;
    834   }
    835 };
    836 
    837 /** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
    838  * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
    839  * \param *out output stream for debugging
    840  * \param *TesselStruct pointer to Tesselation structure
    841  */
    842 void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
    843 {
    844   class BoundaryPointSet *point = NULL;
    845   class BoundaryLineSet *line = NULL;
    846 
    847   //*out << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
    848   // calculate remaining concavity
    849   for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
    850     point = PointRunner->second;
    851     *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
    852     point->value = 0;
    853     for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    854       line = LineRunner->second;
    855       //*out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    856       if (!line->CheckConvexityCriterion(out))
    857         point->value += 1;
    858     }
    859   }
    860   //*out << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
    861 };
    862 
    863 
    864 /** Checks whether each BoundaryLineSet in the Tesselation has two triangles.
    865  * \param *out output stream for debugging
    866  * \param *TesselStruct
    867  * \return true - all have exactly two triangles, false - some not, list is printed to screen
    868  */
    869 bool CheckListOfBaselines(ofstream *out, Tesselation *TesselStruct)
    870 {
    871   LineMap::iterator testline;
    872   bool result = false;
    873   int counter = 0;
    874 
    875   *out << Verbose(1) << "Check: List of Baselines with not two connected triangles:" << endl;
    876   for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
    877     if (testline->second->triangles.size() != 2) {
    878       *out << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
    879       counter++;
    880     }
    881   }
    882   if (counter == 0) {
    883     *out << Verbose(1) << "None." << endl;
    884     result = true;
    885   }
    886   return result;
    887 }
    888 
Note: See TracChangeset for help on using the changeset viewer.