Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • molecuilder/src/tesselationhelpers.cpp

    rbe2997 r543ce4  
    88#include <fstream>
    99
    10 #include "info.hpp"
    1110#include "linkedcell.hpp"
    1211#include "log.hpp"
     
    1615#include "verbose.hpp"
    1716
    18 double DetGet(gsl_matrix * const A, const int inPlace)
    19 {
    20         Info FunctionInfo(__func__);
     17double DetGet(gsl_matrix * const A, const int inPlace) {
    2118  /*
    2219  inPlace = 1 => A is replaced with the LU decomposed copy.
     
    4845void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
    4946{
    50         Info FunctionInfo(__func__);
    5147  gsl_matrix *A = gsl_matrix_calloc(3,3);
    5248  double m11, m12, m13, m14;
     
    8177
    8278  if (fabs(m11) < MYEPSILON)
    83     eLog() << Verbose(1) << "three points are colinear." << endl;
     79    eLog() << Verbose(0) << "ERROR: three points are colinear." << endl;
    8480
    8581  center->x[0] =  0.5 * m12/ m11;
     
    8884
    8985  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
    90     eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
     86    eLog() << Verbose(0) << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
    9187
    9288  gsl_matrix_free(A);
     
    115111    const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius)
    116112{
    117         Info FunctionInfo(__func__);
    118113  Vector TempNormal, helper;
    119114  double Restradius;
    120115  Vector OtherCenter;
     116  Log() << Verbose(3) << "Begin of GetCenterOfSphere.\n";
    121117  Center->Zero();
    122118  helper.CopyVector(&a);
     
    132128  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    133129  NewUmkreismittelpunkt->CopyVector(Center);
    134   Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
     130  Log() << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
    135131  // Here we calculated center of circumscribing circle, using barycentric coordinates
    136   Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
     132  Log() << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
    137133
    138134  TempNormal.CopyVector(&a);
     
    158154  TempNormal.Normalize();
    159155  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
    160   Log() << Verbose(1) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
     156  Log() << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
    161157  TempNormal.Scale(Restradius);
    162   Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
     158  Log() << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
    163159
    164160  Center->AddVector(&TempNormal);
    165   Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n";
     161  Log() << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
    166162  GetSphere(&OtherCenter, a, b, c, RADIUS);
    167   Log() << Verbose(1) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
     163  Log() << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
     164  Log() << Verbose(3) << "End of GetCenterOfSphere.\n";
    168165};
    169166
     
    177174void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c)
    178175{
    179         Info FunctionInfo(__func__);
    180176  Vector helper;
    181177  double alpha, beta, gamma;
     
    190186  beta = M_PI - SideC.Angle(&SideA);
    191187  gamma = M_PI - SideA.Angle(&SideB);
    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   }
     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;
    196191
    197192  Center->Zero();
     
    223218double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection)
    224219{
    225         Info FunctionInfo(__func__);
    226220  Vector helper;
    227221  double radius, alpha;
     
    230224  // test whether new center is on the parameter circle's plane
    231225  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    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;
     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;
    233227    helper.ProjectOntoPlane(&CirclePlaneNormal);
    234228  }
     
    236230  // test whether the new center vector has length of CircleRadius
    237231  if (fabs(radius - CircleRadius) > HULLEPSILON)
    238     eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
     232    eLog() << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    239233  alpha = helper.Angle(&OldSphereCenter);
    240234  // make the angle unique by checking the halfplanes/search direction
    241235  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    242236    alpha = 2.*M_PI - alpha;
    243   //Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
     237  //Log() << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
    244238  radius = helper.Distance(&OldSphereCenter);
    245239  helper.ProjectOntoPlane(&NormalVector);
    246240  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    247241  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    248     //Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
     242    //Log() << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
    249243    return alpha;
    250244  } else {
     
    269263double MinIntersectDistance(const gsl_vector * x, void *params)
    270264{
    271         Info FunctionInfo(__func__);
    272265  double retval = 0;
    273266  struct Intersection *I = (struct Intersection *)params;
     
    290283
    291284  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
    292   //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
     285  //Log() << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
    293286
    294287  return retval;
     
    310303bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
    311304{
    312         Info FunctionInfo(__func__);
    313305  bool result;
    314306
     
    358350
    359351        if (status == GSL_SUCCESS) {
    360           Log() << Verbose(1) << "converged to minimum" <<  endl;
     352          Log() << Verbose(2) << "converged to minimum" <<  endl;
    361353        }
    362354    } while (status == GSL_CONTINUE && iter < 100);
     
    383375  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
    384376
    385   Log() << Verbose(1) << "Intersection " << intersection << " is at "
     377  Log() << Verbose(2) << "Intersection " << intersection << " is at "
    386378    << t1 << " for (" << point1 << "," << point2 << ") and at "
    387379    << t2 << " for (" << point3 << "," << point4 << "): ";
    388380
    389381  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    390     Log() << Verbose(1) << "true intersection." << endl;
     382    Log() << Verbose(0) << "true intersection." << endl;
    391383    result = true;
    392384  } else {
    393     Log() << Verbose(1) << "intersection out of region of interest." << endl;
     385    Log() << Verbose(0) << "intersection out of region of interest." << endl;
    394386    result = false;
    395387  }
     
    414406double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector)
    415407{
    416         Info FunctionInfo(__func__);
    417408  if (reference.IsZero())
    418409    return M_PI;
     
    426417  }
    427418
    428   Log() << Verbose(1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
     419  Log() << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
    429420
    430421  return phi;
     
    441432double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
    442433{
    443         Info FunctionInfo(__func__);
    444434  Vector Point, TetraederVector[3];
    445435  double volume;
     
    465455bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3])
    466456{
    467         Info FunctionInfo(__func__);
    468457  bool result = false;
    469458  int counter = 0;
     
    472461  for (int i=0;i<3;i++)
    473462    for (int j=i+1; j<3; j++) {
    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
     463      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
    478464        LineMap::const_iterator FindLine;
    479465        pair<LineMap::const_iterator,LineMap::const_iterator> FindPair;
     
    492478    }
    493479  if ((!result) && (counter > 1)) {
    494     Log() << Verbose(1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
     480    Log() << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    495481    result = true;
    496482  }
     
    499485
    500486
    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 //};
     487/** Sort function for the candidate list.
     488 */
     489bool 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};
    545530
    546531/**
     
    554539TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
    555540{
    556         Info FunctionInfo(__func__);
    557541  TesselPoint* closestPoint = NULL;
    558542  TesselPoint* secondClosestPoint = NULL;
     
    565549  for(int i=0;i<NDIM;i++) // store indices of this cell
    566550    N[i] = LC->n[i];
    567   Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     551  Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    568552
    569553  LC->GetNeighbourBounds(Nlower, Nupper);
    570   //Log() << Verbose(1) << endl;
     554  //Log() << Verbose(0) << endl;
    571555  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    572556    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    573557      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    574558        const LinkedNodes *List = LC->GetCurrentCell();
    575         //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     559        //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    576560        if (List != NULL) {
    577561          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     
    590574          }
    591575        } else {
    592           eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
     576          eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
    593577            << LC->n[2] << " is invalid!" << endl;
    594578        }
     
    609593TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
    610594{
    611         Info FunctionInfo(__func__);
    612595  TesselPoint* closestPoint = NULL;
    613596  SecondPoint = NULL;
     
    620603  for(int i=0;i<NDIM;i++) // store indices of this cell
    621604    N[i] = LC->n[i];
    622   Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     605  Log() << Verbose(3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    623606
    624607  LC->GetNeighbourBounds(Nlower, Nupper);
    625   //Log() << Verbose(1) << endl;
     608  //Log() << Verbose(0) << endl;
    626609  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    627610    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    628611      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    629612        const LinkedNodes *List = LC->GetCurrentCell();
    630         //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     613        //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    631614        if (List != NULL) {
    632615          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     
    639622              distance = currentNorm;
    640623              closestPoint = (*Runner);
    641               //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
     624              //Log() << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    642625            } else if (currentNorm < secondDistance) {
    643626              secondDistance = currentNorm;
    644627              SecondPoint = (*Runner);
    645               //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
     628              //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
    646629            }
    647630          }
    648631        } else {
    649           eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
     632          eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
    650633            << LC->n[2] << " is invalid!" << endl;
    651634        }
     
    653636  // output
    654637  if (closestPoint != NULL) {
    655     Log() << Verbose(1) << "Closest point is " << *closestPoint;
     638    Log() << Verbose(2) << "Closest point is " << *closestPoint;
    656639    if (SecondPoint != NULL)
    657640      Log() << Verbose(0) << " and second closest is " << *SecondPoint;
     
    669652Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase)
    670653{
    671         Info FunctionInfo(__func__);
    672654  // construct the plane of the two baselines (i.e. take both their directional vectors)
    673655  Vector Normal;
     
    680662  Normal.VectorProduct(&OtherBaseline);
    681663  Normal.Normalize();
    682   Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
     664  Log() << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
    683665
    684666  // project one offset point of OtherBase onto this plane (and add plane offset vector)
     
    697679  Normal.CopyVector(Intersection);
    698680  Normal.SubtractVector(Base->endpoints[0]->node->node);
    699   Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
     681  Log() << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
    700682
    701683  return Intersection;
     
    710692double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle)
    711693{
    712         Info FunctionInfo(__func__);
    713694  double distance = 0.;
    714695  if (x == NULL) {
     
    727708void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud)
    728709{
    729         Info FunctionInfo(__func__);
    730710  TesselPoint *Walker = NULL;
    731711  int i;
     
    758738    }
    759739  } else {
    760     eLog() << Verbose(1) << "Given vrmlfile is " << vrmlfile << "." << endl;
     740    eLog() << Verbose(0) << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
    761741  }
    762742  delete(center);
     
    771751void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    772752{
    773         Info FunctionInfo(__func__);
    774753  Vector helper;
    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   }
     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);
    792768};
    793769
     
    800776void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    801777{
    802         Info FunctionInfo(__func__);
    803778  TesselPoint *Walker = NULL;
    804779  int i;
     
    833808    *rasterfile << "9\n#  terminating special property\n";
    834809  } else {
    835     eLog() << Verbose(1) << "Given rasterfile is " << rasterfile << "." << endl;
     810    eLog() << Verbose(0) << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
    836811  }
    837812  IncludeSphereinRaster3D(rasterfile, Tess, cloud);
     
    846821void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N)
    847822{
    848         Info FunctionInfo(__func__);
    849823  if ((tecplot != NULL) && (TesselStruct != NULL)) {
    850824    // write header
    851825    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    852826    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
    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     }
     827    *tecplot << "ZONE T=\"" << N << "-";
     828    for (int i=0;i<3;i++)
     829      *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
    861830    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    862831    int i=0;
     
    867836
    868837    // print atom coordinates
     838    Log() << Verbose(2) << "The following triangles were created:";
    869839    int Counter = 1;
    870840    TesselPoint *Walker = NULL;
     
    876846    *tecplot << endl;
    877847    // print connectivity
    878     Log() << Verbose(1) << "The following triangles were created:" << endl;
    879848    for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
    880       Log() << Verbose(1) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name << endl;
     849      Log() << Verbose(0) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
    881850      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    882851    }
    883852    delete[] (LookupList);
     853    Log() << Verbose(0) << endl;
    884854  }
    885855};
     
    892862void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct)
    893863{
    894         Info FunctionInfo(__func__);
    895864  class BoundaryPointSet *point = NULL;
    896865  class BoundaryLineSet *line = NULL;
    897866
     867  //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
    898868  // calculate remaining concavity
    899869  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     
    903873    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    904874      line = LineRunner->second;
    905       //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     875      //Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    906876      if (!line->CheckConvexityCriterion())
    907877        point->value += 1;
    908878    }
    909879  }
     880  //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
    910881};
    911882
     
    918889bool CheckListOfBaselines(const Tesselation * const TesselStruct)
    919890{
    920         Info FunctionInfo(__func__);
    921891  LineMap::const_iterator testline;
    922892  bool result = false;
     
    926896  for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
    927897    if (testline->second->triangles.size() != 2) {
    928       Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl;
     898      Log() << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
    929899      counter++;
    930900    }
Note: See TracChangeset for help on using the changeset viewer.