Changeset 44fd95 for src/boundary.cpp


Ignore:
Timestamp:
Dec 23, 2008, 11:22:57 AM (16 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
02bfde
Parents:
7c6712
Message:

Cleaned up all debugging output

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    r7c6712 r44fd95  
    22#include "boundary.hpp"
    33
    4 #define DEBUG 0
     4#define DEBUG 1
    55
    66// ======================================== Points on Boundary =================================
    … …  
    14351435    Vector TempNormal, helper;
    14361436    double Restradius;
    1437 
     1437    cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
    14381438    *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
    14391439    Center->Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
    14401440    // Here we calculated center of circumscribing circle, using barycentric coordinates
     1441    cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
    14411442
    14421443    TempNormal.CopyVector(&a);
    … …  
    14521453    TempNormal.Normalize();
    14531454    Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
     1455    cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
    14541456    TempNormal.Scale(Restradius);
     1457    cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
    14551458
    14561459    Center->AddVector(&TempNormal);
     1460    cout << Verbose(4) << "Center of sphere of circumference is " << *Center << ".\n";
     1461    cout << Verbose(3) << "End of Get_center_of_sphere.\n";
    14571462  }
    14581463  ;
    … …  
    14861491    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
    14871492{
    1488 
    1489   cout << Verbose(1) << "Candidate is "<< Candidate->nr << endl;
    1490   cout << "ReferencePoint is " << ReferencePoint.x[0] << " "<< ReferencePoint.x[1] << " "<< ReferencePoint.x[2] << " "<< endl;
     1493  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
     1494  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
     1495  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
     1496  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
     1497  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
     1498  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
    14911499  /* OldNormal is normal vector on the old triangle
    14921500   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
    … …  
    15051513
    15061514
    1507   dif_a.CopyVector(&(a->x));
    1508   dif_a.SubtractVector(&(Candidate->x));
    1509   dif_b.CopyVector(&(b->x));
    1510   dif_b.SubtractVector(&(Candidate->x));
    1511   AngleCheck.CopyVector(&(Candidate->x));
    1512   AngleCheck.SubtractVector(&(a->x));
    1513   AngleCheck.ProjectOntoPlane(Chord);
    1514 
    1515   SideA = dif_b.Norm();
    1516   SideB = dif_a.Norm();
    1517   SideC = Chord->Norm();
    1518   //Chord->Scale(-1);
    1519 
    1520   alpha = Chord->Angle(&dif_a);
    1521   beta = M_PI - Chord->Angle(&dif_b);
    1522   gamma = dif_a.Angle(&dif_b);
    1523 
    15241515  if (a != Candidate and b != Candidate)
    15251516    {
     1517      cout << Verbose(3) << "We have a unique candidate!" << endl;
     1518      dif_a.CopyVector(&(a->x));
     1519      dif_a.SubtractVector(&(Candidate->x));
     1520      dif_b.CopyVector(&(b->x));
     1521      dif_b.SubtractVector(&(Candidate->x));
     1522      AngleCheck.CopyVector(&(Candidate->x));
     1523      AngleCheck.SubtractVector(&(a->x));
     1524      AngleCheck.ProjectOntoPlane(Chord);
     1525
     1526      SideA = dif_b.Norm();
     1527      SideB = dif_a.Norm();
     1528      SideC = Chord->Norm();
     1529      //Chord->Scale(-1);
     1530
     1531      alpha = Chord->Angle(&dif_a);
     1532      beta = M_PI - Chord->Angle(&dif_b);
     1533      gamma = dif_a.Angle(&dif_b);
     1534
     1535      cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
     1536
     1537      if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
     1538        cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
    15261539
    15271540      Umkreisradius = SideA / 2.0 / sin(alpha);
    … …  
    15321545      if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
    15331546        {
     1547          cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
    15341548          sign = AngleCheck.ScalarProduct(direction1);
    1535           sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
    1536 
    1537           Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, sign, alpha, beta, gamma, RADIUS, Umkreisradius);
    1538 
    1539           AngleCheck.CopyVector(&ReferencePoint);
    1540           AngleCheck.Scale(-1);
    1541           cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
    1542           AngleCheck.AddVector(&Mittelpunkt);
    1543           cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
    1544 
    1545           BallAngle = AngleCheck.Angle(OldNormal);
    1546           cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
    1547           cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
    1548           sign = AngleCheck.ScalarProduct(direction1);
    1549           sign /= fabs(sign);
    1550 
    1551 
    1552           if (sign >0)
     1549          if (fabs(sign) < MYEPSILON)
     1550            sign = 0;
     1551          else
     1552            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
     1553          if (sign >= 0)
    15531554            {
     1555            cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
     1556
     1557            Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, sign, alpha, beta, gamma, RADIUS, Umkreisradius);
     1558
     1559            AngleCheck.CopyVector(&ReferencePoint);
     1560            AngleCheck.Scale(-1);
     1561            //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
     1562            AngleCheck.AddVector(&Mittelpunkt);
     1563            //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
     1564            cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
     1565
     1566            BallAngle = AngleCheck.Angle(OldNormal);
     1567            cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
     1568  //          cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
     1569  //          cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
     1570            sign = AngleCheck.ScalarProduct(direction1);
     1571            if (fabs(sign) < MYEPSILON)
     1572              sign = 0;
     1573            else
     1574              sign /= fabs(sign);
     1575
     1576
    15541577              if (Storage[0]< -1.5) // first Candidate at all
    15551578                {
    15561579
    1557                   cout << "Next better candidate is " << *Candidate << " with ";
     1580                  cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
    15581581                  Opt_Candidate = Candidate;
    15591582                  Storage[0] = sign;
    15601583                  Storage[1] = BallAngle;
    1561                   cout << "Angle is " << Storage[1] << ", Halbraum ist "
     1584                  cout << " angle " << Storage[1] << " and Up/Down "
    15621585                    << Storage[0] << endl;
    1563 
    1564 
    15651586                }
    15661587              else
    … …  
    15681589                  if ( Storage[1] > BallAngle)
    15691590                    {
    1570                       cout << "Next better candidate is " << *Candidate << " with ";
     1591                      cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
    15711592                      Opt_Candidate = Candidate;
    15721593                      Storage[0] = sign;
    15731594                      Storage[1] = BallAngle;
    1574                       cout << "Angle is " << Storage[1] << ", Halbraum ist "
     1595                      cout << " angle " << Storage[1] << " and Up/Down "
    15751596                        << Storage[0] << endl;
    15761597                    }
    … …  
    15791600                      if (DEBUG)
    15801601                        {
    1581                           cout << "Looses to better candidate" << endl;
     1602                          cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
    15821603                        }
    15831604                    }
    … …  
    15861607            else
    15871608              {
    1588                 cout << "Refused due to sign which is " << sign << endl;
     1609              if (DEBUG)
     1610                {
     1611                  cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
     1612                }
    15891613              }
    15901614          }
    … …  
    15931617            if (DEBUG)
    15941618              {
    1595                 cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
     1619                cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
    15961620              }
    15971621          }
    … …  
    16011625        if (DEBUG)
    16021626          {
    1603             cout << "identisch mit Ursprungslinie" << endl;
     1627            cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
    16041628          }
    16051629      }
    … …  
    16251649        }
    16261650    }
     1651  cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
    16271652}
    16281653;
    … …  
    17461771          d1->ProjectOntoPlane(&AngleCheckReference);
    17471772          sign = AngleCheck.ScalarProduct(d1);
    1748           sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
     1773          if (fabs(sign) < MYEPSILON)
     1774            sign = 0;
     1775          else
     1776            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
    17491777
    17501778
    … …  
    19041932    const double& RADIUS, int N)
    19051933{
    1906   cout << Verbose(1) << "Looking for next suitable triangle \n";
     1934  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
    19071935  Vector direction1;
    19081936  Vector helper;
    … …  
    19101938  ofstream *tempstream = NULL;
    19111939  char filename[255];
     1940  double tmp;
    19121941  //atom* Walker;
    19131942
    … …  
    19181947  Vector Opt_Mittelpunkt;
    19191948
    1920   cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
     1949  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
     1950
    19211951  helper.CopyVector(&(Line.endpoints[0]->node->x));
    19221952  for (int i = 0; i < 3; i++)
    … …  
    19381968      direction1.Scale(-1);
    19391969    }
     1970  cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
    19401971
    19411972  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
    19421973  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
    1943 
    1944 
    1945     Vector Umkreismittelpunkt, a, b, c;
    1946     double alpha, beta, gamma;
    1947     a.CopyVector(&(T.endpoints[0]->node->x));
    1948     b.CopyVector(&(T.endpoints[1]->node->x));
    1949     c.CopyVector(&(T.endpoints[2]->node->x));
    1950     a.SubtractVector(&(T.endpoints[1]->node->x));
    1951     b.SubtractVector(&(T.endpoints[2]->node->x));
    1952     c.SubtractVector(&(T.endpoints[0]->node->x));
    1953 
    1954     alpha = M_PI - a.Angle(&c);
    1955     beta = M_PI - b.Angle(&a);
    1956     gamma = M_PI - c.Angle(&b);
    1957 
    1958     Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
    1959     cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
    1960     Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
    1961     cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
    1962 
     1974  cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
     1975
     1976
     1977  Vector Umkreismittelpunkt, a, b, c;
     1978  double alpha, beta, gamma;
     1979  a.CopyVector(&(T.endpoints[0]->node->x));
     1980  b.CopyVector(&(T.endpoints[1]->node->x));
     1981  c.CopyVector(&(T.endpoints[2]->node->x));
     1982  a.SubtractVector(&(T.endpoints[1]->node->x));
     1983  b.SubtractVector(&(T.endpoints[2]->node->x));
     1984  c.SubtractVector(&(T.endpoints[0]->node->x));
     1985
     1986  alpha = a.Angle(&c) - M_PI/2.;
     1987  beta = b.Angle(&a);
     1988  gamma = c.Angle(&b) - M_PI/2.;
     1989  if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
     1990    cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
     1991
     1992  cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
     1993  Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
     1994  Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
     1995  cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
     1996  if (DEBUG)
     1997    cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
     1998  tmp = 0;
     1999  for (int i=0;i<NDIM;i++) {
     2000    helper.CopyVector(&T.endpoints[i]->node->x);
     2001    helper.SubtractVector(&Umkreismittelpunkt);
     2002    if (DEBUG)
     2003      cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
     2004    if (tmp == 0) // set on first time for comparison against next ones
     2005      tmp = helper.Norm();
     2006    if (fabs(helper.Norm() - tmp) > MYEPSILON)
     2007      cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
     2008  }
    19632009
    19642010  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
    … …  
    19742020  if ((TrianglesOnBoundaryCount % 10) == 0)
    19752021    {
     2022    cout << Verbose(1) << "Writing temporary non convex hull to file testEnvelope-" << TriangleFilesWritten << "d.dat.\n";
    19762023    sprintf(filename, "testEnvelope-%d.dat", TriangleFilesWritten);
    19772024    tempstream = new ofstream(filename, ios::trunc);
    … …  
    19842031  if (TrianglesOnBoundaryCount >1000 )
    19852032    {
    1986       cout << Verbose(1)
    1987           << "No new Atom found, triangle construction will crash" << endl;
     2033      cerr << Verbose(0)
     2034          << "WARNING: Already more than " << TrianglesOnBoundaryCount-1 << "triangles, construction probably has crashed." << endl;
    19882035      write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
    1989       cout << "This is currently added candidate" << Opt_Candidate << endl;
     2036      cout << Verbose(2) << "This is currently added candidate" << Opt_Candidate << endl;
    19902037    }
    19912038  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
    19922039
    1993   cout << " Optimal candidate is " << *Opt_Candidate << endl;
     2040  cout << Verbose(2) << " Optimal candidate is " << *Opt_Candidate << endl;
    19942041
    19952042  AddTrianglePoint(Opt_Candidate, 0);
    … …  
    20032050  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    20042051  AddTriangleToLines();
    2005   cout << "New triangle with " << *BTS << endl;
    2006   cout << "We have "<< TrianglesOnBoundaryCount << endl;
    2007   cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
    2008 
    2009   BTS->GetNormalVector(BTS->NormalVector);
    2010 
    2011   if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
    2012       (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0))
    2013     {
    2014       BTS->NormalVector.Scale(-1);
    2015     };
    2016 
     2052
     2053  BTS->GetNormalVector(T.NormalVector);
     2054
     2055//  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
     2056//      (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0))
     2057//    {
     2058//      BTS->NormalVector.Scale(-1);
     2059//    };
     2060  cout << Verbose(2) << "New triangle with " << *BTS << "and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
     2061  cout << Verbose(2) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
     2062  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
    20172063}
    20182064;
    … …  
    20222068    molecule* mol, double RADIUS)
    20232069{
    2024   cout << Verbose(1)
    2025       << "Looking for second point of starting triangle, recursive level "
    2026       << RecursionLevel << endl;;
     2070  cout << Verbose(2)
     2071      << "Begin of Find_second_point_for_Tesselation, recursive level "
     2072      << RecursionLevel << endl;
    20272073  int i;
    20282074  Vector AngleCheck;
    20292075  atom* Walker;
    2030   double norm = -1.;
     2076  double norm = -1., angle;
    20312077
    20322078  // check if we only have one unique point yet ...
    20332079  if (a != Candidate)
    20342080    {
     2081    cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
    20352082      AngleCheck.CopyVector(&(Candidate->x));
    20362083      AngleCheck.SubtractVector(&(a->x));
    … …  
    20392086      if (norm < RADIUS)
    20402087        {
    2041           if (AngleCheck.Angle(&Oben) < Storage[0])
     2088        angle = AngleCheck.Angle(&Oben);
     2089          if (angle < Storage[0])
    20422090            {
    20432091              //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
    2044               cout << "Next better candidate is " << *Candidate
    2045                   << " with distance " << norm << ".\n";
     2092              cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
    20462093              Opt_Candidate = Candidate;
    20472094              Storage[0] = AngleCheck.Angle(&Oben);
    … …  
    20502097          else
    20512098            {
    2052               cout << Verbose(1) << "Supposedly looses to a better candidate "
     2099              cout << "Looses with angle " << angle << " to a better candidate "
    20532100                  << *Opt_Candidate << endl;
    20542101            }
    … …  
    20562103      else
    20572104        {
    2058           cout << Verbose(1) << *Candidate << " refused due to Radius " << norm
     2105          cout << "Refused due to Radius " << norm
    20592106              << endl;
    20602107        }
    … …  
    20752122        };
    20762123    };
     2124  cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level "
     2125      << RecursionLevel << endl;
    20772126}
    20782127;
    … …  
    20802129void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
    20812130{
    2082   cout << Verbose(1) << "Looking for starting triangle \n";
     2131  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
    20832132  int i = 0;
    20842133  atom* Walker;
    20852134  atom* FirstPoint;
    20862135  atom* SecondPoint;
    2087   atom* max_index[3];
    2088   double max_coordinate[3];
     2136  atom* max_index[NDIM];
     2137  double max_coordinate[NDIM];
    20892138  Vector Oben;
    20902139  Vector helper;
    … …  
    20992148      max_coordinate[i] = -1;
    21002149    }
    2101   cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount
     2150  cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount
    21022151      << " Atoms \n";
    21032152
    … …  
    21172166    }
    21182167
    2119   cout << Verbose(1) << "Found maximum coordinates. " << endl;
     2168  cout << Verbose(2) << "Found maximum coordinates: ";
     2169  for (int i=0;i<NDIM;i++)
     2170    cout << i << ": " << *max_index[i] << "\t";
     2171  cout << endl;
    21202172  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
    21212173  const int k = 1;
    … …  
    21232175  FirstPoint = max_index[k];
    21242176
    2125   cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << ": "
    2126       << FirstPoint->x.x[0] << endl;
     2177  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
    21272178  double Storage[2];
    21282179  atom* Opt_Candidate = NULL;
    21292180  Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
    21302181  Storage[1] = 999999.; // This will be an angle looking for the third point.
    2131   cout << Verbose(1) << "Number of Bonds: "
    2132       << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
    21332182
    21342183  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
    21352184      Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
    21362185  SecondPoint = Opt_Candidate;
    2137   cout << Verbose(1) << "Found second point is " << *SecondPoint << ".\n";
     2186  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
    21382187
    21392188  helper.CopyVector(&(FirstPoint->x));
    … …  
    21502199  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
    21512200
    2152   cout << Verbose(1) << "Looking for third point candidates \n";
     2201  cout << Verbose(2) << "Looking for third point candidates \n";
    21532202  // look in one direction of baseline for initial candidate
    21542203  Opt_Candidate = NULL;
    … …  
    21572206  CenterOfFirstLine.AddVector(&(SecondPoint->x));
    21582207
     2208  cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
    21592209  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0,
    21602210      &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
    21612211  // look in other direction of baseline for possible better candidate
     2212  cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
    21622213  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, FirstPoint, SecondPoint, 0,
    21632214      &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
    … …  
    21652216
    21662217  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
    2167 
    2168   cout << Verbose(1) << "The found starting triangle consists of "
    2169       << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
    2170       << "." << endl;
    21712218
    21722219  // Finally, we only have to add the found points
    … …  
    21842231  Oben.Scale(-1.);
    21852232  BTS->GetNormalVector(Oben);
     2233  cout << Verbose(0) << "==> The found starting triangle consists of "
     2234      << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
     2235      << " with normal vector " << BTS->NormalVector << ".\n";
     2236  cout << Verbose(1) << "End of Find_starting_triangle\n";
    21862237}
    21872238;
    … …  
    21952246  const double RADIUS = 6.;
    21962247  LineMap::iterator baseline;
     2248  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
     2249  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
     2250
    21972251  Tess->Find_starting_triangle(mol, RADIUS);
    21982252
    21992253  baseline = Tess->LinesOnBoundary.begin();
    2200   while (baseline != Tess->LinesOnBoundary.end())
     2254  while ((baseline != Tess->LinesOnBoundary.end()) || (!flag))
    22012255    {
    22022256      if (baseline->second->TrianglesCount == 1)
    22032257        {
    2204           cout << Verbose(1) << "Begin of Tesselation ... " << endl;
    22052258          Tess->Find_next_suitable_triangle(out, tecplot, mol,
    22062259              *(baseline->second),
    22072260              *(((baseline->second->triangles.begin()))->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
    2208           cout << Verbose(1) << "End of Tesselation ... " << endl;
     2261          flag = true;
    22092262        }
    22102263      else
    22112264        {
    2212           cout << Verbose(1) << "There is a line with "
     2265          cout << Verbose(1) << "Line " << &baseline << " has "
    22132266              << baseline->second->TrianglesCount << " triangles adjacent"
    22142267              << endl;
    … …  
    22162269      N++;
    22172270      baseline++;
    2218     }
     2271      if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
     2272        baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
     2273        flag = false;
     2274      }
     2275    }
     2276  cout << Verbose(1) << "Writing final tecplot file\n";
    22192277  write_tecplot_file(out, tecplot, Tess, mol, -1);
    22202278
    2221 }
    2222 ;
     2279  cout << Verbose(0) << "End of Find_non_convex_border\n";
     2280}
     2281;
Note: See TracChangeset for help on using the changeset viewer.