Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselationhelpers.cpp

    re138de r299554  
    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    DoeLog(1) && (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    DoeLog(1) && (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  DoLog(1) && (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  DoLog(1) && (Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n");
    133137
    134138  TempNormal.CopyVector(&a);
    … …  
    139143  if (fabs(HalfplaneIndicator) < MYEPSILON)
    140144    {
    141       if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
     145      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(AlternativeDirection) >0 && AlternativeIndicator <0))
    142146        {
    143147          TempNormal.Scale(-1);
    … …  
    146150  else
    147151    {
    148       if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
     152      if (((TempNormal.ScalarProduct(Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(Direction)>0) && (HalfplaneIndicator<0)))
    149153        {
    150154          TempNormal.Scale(-1);
    … …  
    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  DoLog(1) && (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  DoLog(1) && (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  DoLog(1) && (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  DoLog(1) && (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    DoeLog(2) && (eLog()<< Verbose(2) << "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;
    222 
    223   helper.CopyVector(&NewSphereCenter);
     228  Vector RelativeOldSphereCenter;
     229  Vector RelativeNewSphereCenter;
     230
     231  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
     232  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
     233  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
     234  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
     235  helper.CopyVector(&RelativeNewSphereCenter);
    224236  // test whether new center is on the parameter circle's plane
    225237  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;
     238    DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl);
    227239    helper.ProjectOntoPlane(&CirclePlaneNormal);
    228240  }
    229   radius = helper.ScalarProduct(&helper);
     241  radius = helper.NormSquared();
    230242  // test whether the new center vector has length of CircleRadius
    231243  if (fabs(radius - CircleRadius) > HULLEPSILON)
    232     eLog() << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    233   alpha = helper.Angle(&OldSphereCenter);
     244    DoeLog(1) && (eLog()<< Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
     245  alpha = helper.Angle(&RelativeOldSphereCenter);
    234246  // make the angle unique by checking the halfplanes/search direction
    235247  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    236248    alpha = 2.*M_PI - alpha;
    237   //Log() << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
    238   radius = helper.Distance(&OldSphereCenter);
     249  DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl);
     250  radius = helper.Distance(&RelativeOldSphereCenter);
    239251  helper.ProjectOntoPlane(&NormalVector);
    240252  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    241253  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;
     254    DoLog(1) && (Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl);
    243255    return alpha;
    244256  } else {
    245     //Log() << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
     257    DoLog(1) && (Log() << Verbose(1) << "INFO: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "." << endl);
    246258    return 2.*M_PI;
    247259  }
    … …  
    263275double MinIntersectDistance(const gsl_vector * x, void *params)
    264276{
     277        Info FunctionInfo(__func__);
    265278  double retval = 0;
    266279  struct Intersection *I = (struct Intersection *)params;
    … …  
    283296
    284297  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
    285   //Log() << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
     298  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
    286299
    287300  return retval;
    … …  
    303316bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
    304317{
     318        Info FunctionInfo(__func__);
    305319  bool result;
    306320
    … …  
    350364
    351365        if (status == GSL_SUCCESS) {
    352           Log() << Verbose(2) << "converged to minimum" <<  endl;
     366          DoLog(1) && (Log() << Verbose(1) << "converged to minimum" <<  endl);
    353367        }
    354368    } while (status == GSL_CONTINUE && iter < 100);
    … …  
    375389  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
    376390
    377   Log() << Verbose(2) << "Intersection " << intersection << " is at "
     391  Log() << Verbose(1) << "Intersection " << intersection << " is at "
    378392    << t1 << " for (" << point1 << "," << point2 << ") and at "
    379393    << t2 << " for (" << point3 << "," << point4 << "): ";
    380394
    381395  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    382     Log() << Verbose(0) << "true intersection." << endl;
     396    DoLog(1) && (Log() << Verbose(1) << "true intersection." << endl);
    383397    result = true;
    384398  } else {
    385     Log() << Verbose(0) << "intersection out of region of interest." << endl;
     399    DoLog(1) && (Log() << Verbose(1) << "intersection out of region of interest." << endl);
    386400    result = false;
    387401  }
    … …  
    406420double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector)
    407421{
     422        Info FunctionInfo(__func__);
    408423  if (reference.IsZero())
    409424    return M_PI;
    … …  
    417432  }
    418433
    419   Log() << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
     434  DoLog(1) && (Log() << Verbose(1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl);
    420435
    421436  return phi;
    … …  
    432447double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
    433448{
     449        Info FunctionInfo(__func__);
    434450  Vector Point, TetraederVector[3];
    435451  double volume;
    … …  
    455471bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3])
    456472{
     473        Info FunctionInfo(__func__);
    457474  bool result = false;
    458475  int counter = 0;
    … …  
    461478  for (int i=0;i<3;i++)
    462479    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
     480      if (nodes[i] == NULL) {
     481        DoLog(1) && (Log() << Verbose(1) << "Node nr. " << i << " is not yet present." << endl);
     482        result = true;
     483      } else if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
    464484        LineMap::const_iterator FindLine;
    465485        pair<LineMap::const_iterator,LineMap::const_iterator> FindPair;
    … …  
    473493        }
    474494      } else { // no line
    475         Log() << Verbose(1) << "The line between " << *nodes[i] << " and " << *nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
     495        DoLog(1) && (Log() << Verbose(1) << "The line between " << *nodes[i] << " and " << *nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl);
    476496        result = true;
    477497      }
    478498    }
    479499  if ((!result) && (counter > 1)) {
    480     Log() << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
     500    DoLog(1) && (Log() << Verbose(1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl);
    481501    result = true;
    482502  }
    … …  
    485505
    486506
    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 };
     507///** Sort function for the candidate list.
     508// */
     509//bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
     510//{
     511//      Info FunctionInfo(__func__);
     512//  Vector BaseLineVector, OrthogonalVector, helper;
     513//  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
     514//    DoeLog(1) && (eLog()<< Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl);
     515//    //return false;
     516//    exit(1);
     517//  }
     518//  // create baseline vector
     519//  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
     520//  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     521//  BaseLineVector.Normalize();
     522//
     523//  // 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!)
     524//  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
     525//  helper.SubtractVector(candidate1->point->node);
     526//  OrthogonalVector.CopyVector(&helper);
     527//  helper.VectorProduct(&BaseLineVector);
     528//  OrthogonalVector.SubtractVector(&helper);
     529//  OrthogonalVector.Normalize();
     530//
     531//  // calculate both angles and correct with in-plane vector
     532//  helper.CopyVector(candidate1->point->node);
     533//  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     534//  double phi = BaseLineVector.Angle(&helper);
     535//  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     536//    phi = 2.*M_PI - phi;
     537//  }
     538//  helper.CopyVector(candidate2->point->node);
     539//  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     540//  double psi = BaseLineVector.Angle(&helper);
     541//  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     542//    psi = 2.*M_PI - psi;
     543//  }
     544//
     545//  Log() << Verbose(1) << *candidate1->point << " has angle " << phi << endl;
     546//  Log() << Verbose(1) << *candidate2->point << " has angle " << psi << endl;
     547//
     548//  // return comparison
     549//  return phi < psi;
     550//};
    530551
    531552/**
    … …  
    537558 * @return point which is second closest to the provided one
    538559 */
    539 TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
    540 {
     560TesselPoint* FindSecondClosestTesselPoint(const Vector* Point, const LinkedCell* const LC)
     561{
     562        Info FunctionInfo(__func__);
    541563  TesselPoint* closestPoint = NULL;
    542564  TesselPoint* secondClosestPoint = NULL;
    … …  
    549571  for(int i=0;i<NDIM;i++) // store indices of this cell
    550572    N[i] = LC->n[i];
    551   Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     573  DoLog(1) && (Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl);
    552574
    553575  LC->GetNeighbourBounds(Nlower, Nupper);
    554   //Log() << Verbose(0) << endl;
     576  //Log() << Verbose(1) << endl;
    555577  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    556578    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    557579      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    558         const LinkedNodes *List = LC->GetCurrentCell();
    559         //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     580        const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
     581        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    560582        if (List != NULL) {
    561           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     583          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    562584            helper.CopyVector(Point);
    563585            helper.SubtractVector((*Runner)->node);
    … …  
    574596          }
    575597        } else {
    576           eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
     598          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
    577599            << LC->n[2] << " is invalid!" << endl;
    578600        }
    … …  
    591613 * @return point which is closest to the provided one, NULL if none found
    592614 */
    593 TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
    594 {
     615TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
     616{
     617        Info FunctionInfo(__func__);
    595618  TesselPoint* closestPoint = NULL;
    596619  SecondPoint = NULL;
    … …  
    603626  for(int i=0;i<NDIM;i++) // store indices of this cell
    604627    N[i] = LC->n[i];
    605   Log() << Verbose(3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     628  DoLog(1) && (Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl);
    606629
    607630  LC->GetNeighbourBounds(Nlower, Nupper);
    608   //Log() << Verbose(0) << endl;
     631  //Log() << Verbose(1) << endl;
    609632  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    610633    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    611634      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    612         const LinkedNodes *List = LC->GetCurrentCell();
    613         //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     635        const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
     636        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    614637        if (List != NULL) {
    615           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     638          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    616639            helper.CopyVector(Point);
    617640            helper.SubtractVector((*Runner)->node);
    618             double currentNorm = helper. Norm();
     641            double currentNorm = helper.NormSquared();
    619642            if (currentNorm < distance) {
    620643              secondDistance = distance;
    … …  
    622645              distance = currentNorm;
    623646              closestPoint = (*Runner);
    624               //Log() << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
     647              //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    625648            } else if (currentNorm < secondDistance) {
    626649              secondDistance = currentNorm;
    627650              SecondPoint = (*Runner);
    628               //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
     651              //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
    629652            }
    630653          }
    631654        } else {
    632           eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
     655          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
    633656            << LC->n[2] << " is invalid!" << endl;
    634657        }
    … …  
    636659  // output
    637660  if (closestPoint != NULL) {
    638     Log() << Verbose(2) << "Closest point is " << *closestPoint;
     661    DoLog(1) && (Log() << Verbose(1) << "Closest point is " << *closestPoint);
    639662    if (SecondPoint != NULL)
    640       Log() << Verbose(0) << " and second closest is " << *SecondPoint;
    641     Log() << Verbose(0) << "." << endl;
     663      DoLog(0) && (Log() << Verbose(0) << " and second closest is " << *SecondPoint);
     664    DoLog(0) && (Log() << Verbose(0) << "." << endl);
    642665  }
    643666  return closestPoint;
    … …  
    652675Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase)
    653676{
     677        Info FunctionInfo(__func__);
    654678  // construct the plane of the two baselines (i.e. take both their directional vectors)
    655679  Vector Normal;
    … …  
    662686  Normal.VectorProduct(&OtherBaseline);
    663687  Normal.Normalize();
    664   Log() << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
     688  DoLog(1) && (Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl);
    665689
    666690  // project one offset point of OtherBase onto this plane (and add plane offset vector)
    … …  
    679703  Normal.CopyVector(Intersection);
    680704  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;
     705  DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl);
    682706
    683707  return Intersection;
    … …  
    692716double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle)
    693717{
     718        Info FunctionInfo(__func__);
    694719  double distance = 0.;
    695720  if (x == NULL) {
    … …  
    708733void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud)
    709734{
     735        Info FunctionInfo(__func__);
    710736  TesselPoint *Walker = NULL;
    711737  int i;
    … …  
    738764    }
    739765  } else {
    740     eLog() << Verbose(0) << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
     766    DoeLog(1) && (eLog()<< Verbose(1) << "Given vrmlfile is " << vrmlfile << "." << endl);
    741767  }
    742768  delete(center);
    … …  
    751777void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    752778{
     779        Info FunctionInfo(__func__);
    753780  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);
     781
     782  if (Tess->LastTriangle != NULL) {
     783    // include the current position of the virtual sphere in the temporary raster3d file
     784    Vector *center = cloud->GetCenter();
     785    // make the circumsphere's center absolute again
     786    helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
     787    helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
     788    helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
     789    helper.Scale(1./3.);
     790    helper.SubtractVector(center);
     791    // and add to file plus translucency object
     792    *rasterfile << "# current virtual sphere\n";
     793    *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
     794    *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
     795    *rasterfile << "9\n  terminating special property\n";
     796    delete(center);
     797  }
    768798};
    769799
    … …  
    776806void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    777807{
     808        Info FunctionInfo(__func__);
    778809  TesselPoint *Walker = NULL;
    779810  int i;
    … …  
    808839    *rasterfile << "9\n#  terminating special property\n";
    809840  } else {
    810     eLog() << Verbose(0) << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
     841    DoeLog(1) && (eLog()<< Verbose(1) << "Given rasterfile is " << rasterfile << "." << endl);
    811842  }
    812843  IncludeSphereinRaster3D(rasterfile, Tess, cloud);
    … …  
    821852void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N)
    822853{
     854        Info FunctionInfo(__func__);
    823855  if ((tecplot != NULL) && (TesselStruct != NULL)) {
    824856    // write header
    825857    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    826858    *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;
     859    *tecplot << "ZONE T=\"";
     860    if (N < 0) {
     861      *tecplot << cloud->GetName();
     862    } else {
     863      *tecplot << N << "-";
     864      if (TesselStruct->LastTriangle != NULL) {
     865        for (int i=0;i<3;i++)
     866          *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
     867      } else {
     868        *tecplot << "none";
     869      }
     870    }
    830871    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    831     int i=0;
    832     for (cloud->GoToFirst(); !cloud->IsEnd(); cloud->GoToNext(), i++);
     872    int i=cloud->GetMaxId();
    833873    int *LookupList = new int[i];
    834874    for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)
    … …  
    836876
    837877    // print atom coordinates
    838     Log() << Verbose(2) << "The following triangles were created:";
    839878    int Counter = 1;
    840879    TesselPoint *Walker = NULL;
    … …  
    846885    *tecplot << endl;
    847886    // print connectivity
     887    DoLog(1) && (Log() << Verbose(1) << "The following triangles were created:" << endl);
    848888    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;
     889      DoLog(1) && (Log() << Verbose(1) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name << endl);
    850890      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    851891    }
    852892    delete[] (LookupList);
    853     Log() << Verbose(0) << endl;
    854893  }
    855894};
    … …  
    862901void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct)
    863902{
     903        Info FunctionInfo(__func__);
    864904  class BoundaryPointSet *point = NULL;
    865905  class BoundaryLineSet *line = NULL;
    866906
    867   //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
    868907  // calculate remaining concavity
    869908  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
    870909    point = PointRunner->second;
    871     Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
     910    DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
    872911    point->value = 0;
    873912    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    874913      line = LineRunner->second;
    875       //Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     914      //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    876915      if (!line->CheckConvexityCriterion())
    877916        point->value += 1;
    878917    }
    879918  }
    880   //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
    881919};
    882920
    … …  
    889927bool CheckListOfBaselines(const Tesselation * const TesselStruct)
    890928{
     929        Info FunctionInfo(__func__);
    891930  LineMap::const_iterator testline;
    892931  bool result = false;
    893932  int counter = 0;
    894933
    895   Log() << Verbose(1) << "Check: List of Baselines with not two connected triangles:" << endl;
     934  DoLog(1) && (Log() << Verbose(1) << "Check: List of Baselines with not two connected triangles:" << endl);
    896935  for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
    897936    if (testline->second->triangles.size() != 2) {
    898       Log() << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
     937      DoLog(2) && (Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl);
    899938      counter++;
    900939    }
    901940  }
    902941  if (counter == 0) {
    903     Log() << Verbose(1) << "None." << endl;
     942    DoLog(1) && (Log() << Verbose(1) << "None." << endl);
    904943    result = true;
    905944  }
    … …  
    907946}
    908947
     948/** Counts the number of triangle pairs that contain the given polygon.
     949 * \param *P polygon with endpoints to look for
     950 * \param *T set of triangles to create pairs from containing \a *P
     951 */
     952int CountTrianglePairContainingPolygon(const BoundaryPolygonSet * const P, const TriangleSet * const T)
     953{
     954  Info FunctionInfo(__func__);
     955  // check number of endpoints in *P
     956  if (P->endpoints.size() != 4) {
     957    DoeLog(1) && (eLog()<< Verbose(1) << "CountTrianglePairContainingPolygon works only on polygons with 4 nodes!" << endl);
     958    return 0;
     959  }
     960
     961  // check number of triangles in *T
     962  if (T->size() < 2) {
     963    DoeLog(1) && (eLog()<< Verbose(1) << "Not enough triangles to have pairs!" << endl);
     964    return 0;
     965  }
     966
     967  DoLog(0) && (Log() << Verbose(0) << "Polygon is " << *P << endl);
     968  // create each pair, get the endpoints and check whether *P is contained.
     969  int counter = 0;
     970  PointSet Trianglenodes;
     971  class BoundaryPolygonSet PairTrianglenodes;
     972  for(TriangleSet::iterator Walker = T->begin(); Walker != T->end(); Walker++) {
     973    for (int i=0;i<3;i++)
     974      Trianglenodes.insert((*Walker)->endpoints[i]);
     975
     976    for(TriangleSet::iterator PairWalker = Walker; PairWalker != T->end(); PairWalker++) {
     977      if (Walker != PairWalker) { // skip first
     978        PairTrianglenodes.endpoints = Trianglenodes;
     979        for (int i=0;i<3;i++)
     980          PairTrianglenodes.endpoints.insert((*PairWalker)->endpoints[i]);
     981        const int size = PairTrianglenodes.endpoints.size();
     982        if (size == 4) {
     983          DoLog(0) && (Log() << Verbose(0) << " Current pair of triangles: " << **Walker << "," << **PairWalker << " with " << size << " distinct endpoints:" << PairTrianglenodes << endl);
     984          // now check
     985          if (PairTrianglenodes.ContainsPresentTupel(P)) {
     986            counter++;
     987            DoLog(0) && (Log() << Verbose(0) << "  ACCEPT: Matches with " << *P << endl);
     988          } else {
     989            DoLog(0) && (Log() << Verbose(0) << "  REJECT: No match with " << *P << endl);
     990          }
     991        } else {
     992          DoLog(0) && (Log() << Verbose(0) << "  REJECT: Less than four endpoints." << endl);
     993        }
     994      }
     995    }
     996    Trianglenodes.clear();
     997  }
     998  return counter;
     999};
     1000
     1001/** Checks whether two give polygons have two or more points in common.
     1002 * \param *P1 first polygon
     1003 * \param *P2 second polygon
     1004 * \return true - are connected, false = are note
     1005 */
     1006bool ArePolygonsEdgeConnected(const BoundaryPolygonSet * const P1, const BoundaryPolygonSet * const P2)
     1007{
     1008  Info FunctionInfo(__func__);
     1009  int counter = 0;
     1010  for(PointSet::const_iterator Runner = P1->endpoints.begin(); Runner != P1->endpoints.end(); Runner++) {
     1011    if (P2->ContainsBoundaryPoint((*Runner))) {
     1012      counter++;
     1013      DoLog(1) && (Log() << Verbose(1) << *(*Runner) << " of second polygon is found in the first one." << endl);
     1014      return true;
     1015    }
     1016  }
     1017  return false;
     1018};
     1019
     1020/** Combines second into the first and deletes the second.
     1021 * \param *P1 first polygon, contains all nodes on return
     1022 * \param *&P2 second polygon, is deleted.
     1023 */
     1024void CombinePolygons(BoundaryPolygonSet * const P1, BoundaryPolygonSet * &P2)
     1025{
     1026  Info FunctionInfo(__func__);
     1027  pair <PointSet::iterator, bool> Tester;
     1028  for(PointSet::iterator Runner = P2->endpoints.begin(); Runner != P2->endpoints.end(); Runner++) {
     1029    Tester = P1->endpoints.insert((*Runner));
     1030    if (Tester.second)
     1031      DoLog(0) && (Log() << Verbose(0) << "Inserting endpoint " << *(*Runner) << " into first polygon." << endl);
     1032  }
     1033  P2->endpoints.clear();
     1034  delete(P2);
     1035};
     1036
Note: See TracChangeset for help on using the changeset viewer.