Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • molecuilder/src/tesselationhelpers.cpp

    r543ce4 rbe2997  
    88#include <fstream>
    99
     10#include "info.hpp"
    1011#include "linkedcell.hpp"
    1112#include "log.hpp"
     
    1516#include "verbose.hpp"
    1617
    17 double DetGet(gsl_matrix * const A, const int inPlace) {
     18double DetGet(gsl_matrix * const A, const int inPlace)
     19{
     20        Info FunctionInfo(__func__);
    1821  /*
    1922  inPlace = 1 => A is replaced with the LU decomposed copy.
     
    4548void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
    4649{
     50        Info FunctionInfo(__func__);
    4751  gsl_matrix *A = gsl_matrix_calloc(3,3);
    4852  double m11, m12, m13, m14;
     
    7781
    7882  if (fabs(m11) < MYEPSILON)
    79     eLog() << Verbose(0) << "ERROR: three points are colinear." << endl;
     83    eLog() << Verbose(1) << "three points are colinear." << endl;
    8084
    8185  center->x[0] =  0.5 * m12/ m11;
     
    8488
    8589  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
    86     eLog() << Verbose(0) << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
     90    eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
    8791
    8892  gsl_matrix_free(A);
     
    111115    const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius)
    112116{
     117        Info FunctionInfo(__func__);
    113118  Vector TempNormal, helper;
    114119  double Restradius;
    115120  Vector OtherCenter;
    116   Log() << Verbose(3) << "Begin of GetCenterOfSphere.\n";
    117121  Center->Zero();
    118122  helper.CopyVector(&a);
     
    128132  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    129133  NewUmkreismittelpunkt->CopyVector(Center);
    130   Log() << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
     134  Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
    131135  // Here we calculated center of circumscribing circle, using barycentric coordinates
    132   Log() << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
     136  Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
    133137
    134138  TempNormal.CopyVector(&a);
     
    154158  TempNormal.Normalize();
    155159  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
    156   Log() << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
     160  Log() << Verbose(1) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
    157161  TempNormal.Scale(Restradius);
    158   Log() << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
     162  Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
    159163
    160164  Center->AddVector(&TempNormal);
    161   Log() << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
     165  Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n";
    162166  GetSphere(&OtherCenter, a, b, c, RADIUS);
    163   Log() << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
    164   Log() << Verbose(3) << "End of GetCenterOfSphere.\n";
     167  Log() << Verbose(1) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
    165168};
    166169
     
    174177void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c)
    175178{
     179        Info FunctionInfo(__func__);
    176180  Vector helper;
    177181  double alpha, beta, gamma;
     
    186190  beta = M_PI - SideC.Angle(&SideA);
    187191  gamma = M_PI - SideA.Angle(&SideB);
    188   //Log() << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
    189   if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
    190     eLog() << Verbose(0) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
     192  //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
     193  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
     194    eLog() << Verbose(1) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
     195  }
    191196
    192197  Center->Zero();
     
    218223double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection)
    219224{
     225        Info FunctionInfo(__func__);
    220226  Vector helper;
    221227  double radius, alpha;
     
    224230  // test whether new center is on the parameter circle's plane
    225231  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    226     eLog() << Verbose(0) << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
     232    eLog() << Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
    227233    helper.ProjectOntoPlane(&CirclePlaneNormal);
    228234  }
     
    230236  // test whether the new center vector has length of CircleRadius
    231237  if (fabs(radius - CircleRadius) > HULLEPSILON)
    232     eLog() << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
     238    eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    233239  alpha = helper.Angle(&OldSphereCenter);
    234240  // make the angle unique by checking the halfplanes/search direction
    235241  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    236242    alpha = 2.*M_PI - alpha;
    237   //Log() << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
     243  //Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
    238244  radius = helper.Distance(&OldSphereCenter);
    239245  helper.ProjectOntoPlane(&NormalVector);
    240246  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    241247  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    242     //Log() << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
     248    //Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
    243249    return alpha;
    244250  } else {
     
    263269double MinIntersectDistance(const gsl_vector * x, void *params)
    264270{
     271        Info FunctionInfo(__func__);
    265272  double retval = 0;
    266273  struct Intersection *I = (struct Intersection *)params;
     
    283290
    284291  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
    285   //Log() << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
     292  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
    286293
    287294  return retval;
     
    303310bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
    304311{
     312        Info FunctionInfo(__func__);
    305313  bool result;
    306314
     
    350358
    351359        if (status == GSL_SUCCESS) {
    352           Log() << Verbose(2) << "converged to minimum" <<  endl;
     360          Log() << Verbose(1) << "converged to minimum" <<  endl;
    353361        }
    354362    } while (status == GSL_CONTINUE && iter < 100);
     
    375383  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
    376384
    377   Log() << Verbose(2) << "Intersection " << intersection << " is at "
     385  Log() << Verbose(1) << "Intersection " << intersection << " is at "
    378386    << t1 << " for (" << point1 << "," << point2 << ") and at "
    379387    << t2 << " for (" << point3 << "," << point4 << "): ";
    380388
    381389  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    382     Log() << Verbose(0) << "true intersection." << endl;
     390    Log() << Verbose(1) << "true intersection." << endl;
    383391    result = true;
    384392  } else {
    385     Log() << Verbose(0) << "intersection out of region of interest." << endl;
     393    Log() << Verbose(1) << "intersection out of region of interest." << endl;
    386394    result = false;
    387395  }
     
    406414double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector)
    407415{
     416        Info FunctionInfo(__func__);
    408417  if (reference.IsZero())
    409418    return M_PI;
     
    417426  }
    418427
    419   Log() << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
     428  Log() << Verbose(1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
    420429
    421430  return phi;
     
    432441double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
    433442{
     443        Info FunctionInfo(__func__);
    434444  Vector Point, TetraederVector[3];
    435445  double volume;
     
    455465bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3])
    456466{
     467        Info FunctionInfo(__func__);
    457468  bool result = false;
    458469  int counter = 0;
     
    461472  for (int i=0;i<3;i++)
    462473    for (int j=i+1; j<3; j++) {
    463       if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
     474      if (nodes[i] == NULL) {
     475        Log() << Verbose(1) << "Node nr. " << i << " is not yet present." << endl;
     476        result = true;
     477      } else if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
    464478        LineMap::const_iterator FindLine;
    465479        pair<LineMap::const_iterator,LineMap::const_iterator> FindPair;
     
    478492    }
    479493  if ((!result) && (counter > 1)) {
    480     Log() << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
     494    Log() << Verbose(1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    481495    result = true;
    482496  }
     
    485499
    486500
    487 /** Sort function for the candidate list.
    488  */
    489 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
    490 {
    491   Vector BaseLineVector, OrthogonalVector, helper;
    492   if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    493     Log() << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    494     //return false;
    495     exit(1);
    496   }
    497   // create baseline vector
    498   BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
    499   BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    500   BaseLineVector.Normalize();
    501 
    502   // 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!)
    503   helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
    504   helper.SubtractVector(candidate1->point->node);
    505   OrthogonalVector.CopyVector(&helper);
    506   helper.VectorProduct(&BaseLineVector);
    507   OrthogonalVector.SubtractVector(&helper);
    508   OrthogonalVector.Normalize();
    509 
    510   // calculate both angles and correct with in-plane vector
    511   helper.CopyVector(candidate1->point->node);
    512   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    513   double phi = BaseLineVector.Angle(&helper);
    514   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    515     phi = 2.*M_PI - phi;
    516   }
    517   helper.CopyVector(candidate2->point->node);
    518   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    519   double psi = BaseLineVector.Angle(&helper);
    520   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    521     psi = 2.*M_PI - psi;
    522   }
    523 
    524   Log() << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
    525   Log() << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
    526 
    527   // return comparison
    528   return phi < psi;
    529 };
     501///** Sort function for the candidate list.
     502// */
     503//bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
     504//{
     505//      Info FunctionInfo(__func__);
     506//  Vector BaseLineVector, OrthogonalVector, helper;
     507//  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
     508//    eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
     509//    //return false;
     510//    exit(1);
     511//  }
     512//  // create baseline vector
     513//  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
     514//  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     515//  BaseLineVector.Normalize();
     516//
     517//  // 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!)
     518//  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
     519//  helper.SubtractVector(candidate1->point->node);
     520//  OrthogonalVector.CopyVector(&helper);
     521//  helper.VectorProduct(&BaseLineVector);
     522//  OrthogonalVector.SubtractVector(&helper);
     523//  OrthogonalVector.Normalize();
     524//
     525//  // calculate both angles and correct with in-plane vector
     526//  helper.CopyVector(candidate1->point->node);
     527//  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     528//  double phi = BaseLineVector.Angle(&helper);
     529//  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     530//    phi = 2.*M_PI - phi;
     531//  }
     532//  helper.CopyVector(candidate2->point->node);
     533//  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     534//  double psi = BaseLineVector.Angle(&helper);
     535//  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     536//    psi = 2.*M_PI - psi;
     537//  }
     538//
     539//  Log() << Verbose(1) << *candidate1->point << " has angle " << phi << endl;
     540//  Log() << Verbose(1) << *candidate2->point << " has angle " << psi << endl;
     541//
     542//  // return comparison
     543//  return phi < psi;
     544//};
    530545
    531546/**
     
    539554TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
    540555{
     556        Info FunctionInfo(__func__);
    541557  TesselPoint* closestPoint = NULL;
    542558  TesselPoint* secondClosestPoint = NULL;
     
    549565  for(int i=0;i<NDIM;i++) // store indices of this cell
    550566    N[i] = LC->n[i];
    551   Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     567  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    552568
    553569  LC->GetNeighbourBounds(Nlower, Nupper);
    554   //Log() << Verbose(0) << endl;
     570  //Log() << Verbose(1) << endl;
    555571  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    556572    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    557573      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    558574        const LinkedNodes *List = LC->GetCurrentCell();
    559         //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     575        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    560576        if (List != NULL) {
    561577          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     
    574590          }
    575591        } else {
    576           eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
     592          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
    577593            << LC->n[2] << " is invalid!" << endl;
    578594        }
     
    593609TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
    594610{
     611        Info FunctionInfo(__func__);
    595612  TesselPoint* closestPoint = NULL;
    596613  SecondPoint = NULL;
     
    603620  for(int i=0;i<NDIM;i++) // store indices of this cell
    604621    N[i] = LC->n[i];
    605   Log() << Verbose(3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     622  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    606623
    607624  LC->GetNeighbourBounds(Nlower, Nupper);
    608   //Log() << Verbose(0) << endl;
     625  //Log() << Verbose(1) << endl;
    609626  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    610627    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    611628      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    612629        const LinkedNodes *List = LC->GetCurrentCell();
    613         //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     630        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    614631        if (List != NULL) {
    615632          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     
    622639              distance = currentNorm;
    623640              closestPoint = (*Runner);
    624               //Log() << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
     641              //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    625642            } else if (currentNorm < secondDistance) {
    626643              secondDistance = currentNorm;
    627644              SecondPoint = (*Runner);
    628               //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
     645              //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
    629646            }
    630647          }
    631648        } else {
    632           eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
     649          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
    633650            << LC->n[2] << " is invalid!" << endl;
    634651        }
     
    636653  // output
    637654  if (closestPoint != NULL) {
    638     Log() << Verbose(2) << "Closest point is " << *closestPoint;
     655    Log() << Verbose(1) << "Closest point is " << *closestPoint;
    639656    if (SecondPoint != NULL)
    640657      Log() << Verbose(0) << " and second closest is " << *SecondPoint;
     
    652669Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase)
    653670{
     671        Info FunctionInfo(__func__);
    654672  // construct the plane of the two baselines (i.e. take both their directional vectors)
    655673  Vector Normal;
     
    662680  Normal.VectorProduct(&OtherBaseline);
    663681  Normal.Normalize();
    664   Log() << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
     682  Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
    665683
    666684  // project one offset point of OtherBase onto this plane (and add plane offset vector)
     
    679697  Normal.CopyVector(Intersection);
    680698  Normal.SubtractVector(Base->endpoints[0]->node->node);
    681   Log() << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
     699  Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
    682700
    683701  return Intersection;
     
    692710double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle)
    693711{
     712        Info FunctionInfo(__func__);
    694713  double distance = 0.;
    695714  if (x == NULL) {
     
    708727void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud)
    709728{
     729        Info FunctionInfo(__func__);
    710730  TesselPoint *Walker = NULL;
    711731  int i;
     
    738758    }
    739759  } else {
    740     eLog() << Verbose(0) << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
     760    eLog() << Verbose(1) << "Given vrmlfile is " << vrmlfile << "." << endl;
    741761  }
    742762  delete(center);
     
    751771void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    752772{
     773        Info FunctionInfo(__func__);
    753774  Vector helper;
    754   // include the current position of the virtual sphere in the temporary raster3d file
    755   Vector *center = cloud->GetCenter();
    756   // make the circumsphere's center absolute again
    757   helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
    758   helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
    759   helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
    760   helper.Scale(1./3.);
    761   helper.SubtractVector(center);
    762   // and add to file plus translucency object
    763   *rasterfile << "# current virtual sphere\n";
    764   *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
    765   *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
    766   *rasterfile << "9\n  terminating special property\n";
    767   delete(center);
     775
     776  if (Tess->LastTriangle != NULL) {
     777    // include the current position of the virtual sphere in the temporary raster3d file
     778    Vector *center = cloud->GetCenter();
     779    // make the circumsphere's center absolute again
     780    helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
     781    helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
     782    helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
     783    helper.Scale(1./3.);
     784    helper.SubtractVector(center);
     785    // and add to file plus translucency object
     786    *rasterfile << "# current virtual sphere\n";
     787    *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
     788    *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
     789    *rasterfile << "9\n  terminating special property\n";
     790    delete(center);
     791  }
    768792};
    769793
     
    776800void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    777801{
     802        Info FunctionInfo(__func__);
    778803  TesselPoint *Walker = NULL;
    779804  int i;
     
    808833    *rasterfile << "9\n#  terminating special property\n";
    809834  } else {
    810     eLog() << Verbose(0) << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
     835    eLog() << Verbose(1) << "Given rasterfile is " << rasterfile << "." << endl;
    811836  }
    812837  IncludeSphereinRaster3D(rasterfile, Tess, cloud);
     
    821846void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N)
    822847{
     848        Info FunctionInfo(__func__);
    823849  if ((tecplot != NULL) && (TesselStruct != NULL)) {
    824850    // write header
    825851    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    826852    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
    827     *tecplot << "ZONE T=\"" << N << "-";
    828     for (int i=0;i<3;i++)
    829       *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
     853    *tecplot << "ZONE T=\"";
     854    if (N < 0) {
     855      *tecplot << cloud->GetName();
     856    } else {
     857      *tecplot << N << "-";
     858      for (int i=0;i<3;i++)
     859        *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
     860    }
    830861    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    831862    int i=0;
     
    836867
    837868    // print atom coordinates
    838     Log() << Verbose(2) << "The following triangles were created:";
    839869    int Counter = 1;
    840870    TesselPoint *Walker = NULL;
     
    846876    *tecplot << endl;
    847877    // print connectivity
     878    Log() << Verbose(1) << "The following triangles were created:" << endl;
    848879    for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
    849       Log() << Verbose(0) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
     880      Log() << Verbose(1) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name << endl;
    850881      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    851882    }
    852883    delete[] (LookupList);
    853     Log() << Verbose(0) << endl;
    854884  }
    855885};
     
    862892void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct)
    863893{
     894        Info FunctionInfo(__func__);
    864895  class BoundaryPointSet *point = NULL;
    865896  class BoundaryLineSet *line = NULL;
    866897
    867   //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
    868898  // calculate remaining concavity
    869899  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     
    873903    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    874904      line = LineRunner->second;
    875       //Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     905      //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    876906      if (!line->CheckConvexityCriterion())
    877907        point->value += 1;
    878908    }
    879909  }
    880   //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
    881910};
    882911
     
    889918bool CheckListOfBaselines(const Tesselation * const TesselStruct)
    890919{
     920        Info FunctionInfo(__func__);
    891921  LineMap::const_iterator testline;
    892922  bool result = false;
     
    896926  for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
    897927    if (testline->second->triangles.size() != 2) {
    898       Log() << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
     928      Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl;
    899929      counter++;
    900930    }
Note: See TracChangeset for help on using the changeset viewer.