Changes in / [6d4a76:12298c]


Ignore:
Location:
src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    r6d4a76 r12298c  
    14141414  LineMap::iterator LineWalker;
    14151415  //cout << "Manually checking endpoints for line." << endl;
    1416   if ((a->lines.find(b->node->nr))->first == b->node->nr)
     1416  if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr)
    14171417  //If a line is there, how do I recognize that beyond a shadow of a doubt?
    14181418    {
    … …  
    14791479 * @param b vector second point of triangle
    14801480 * @param c vector third point of triangle
     1481 * @param *Umkreismittelpunkt new cneter point of circumference
    14811482 * @param Direction vector indicates up/down
    14821483 * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
    … …  
    14901491 */
    14911492
    1492   void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, Vector* AlternativeDirection,
     1493  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
    14931494      double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
    14941495  {
    14951496    Vector TempNormal, helper;
    14961497    double Restradius;
    1497 
     1498    cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
    14981499    *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
    14991500    Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
     1501    NewUmkreismittelpunkt->CopyVector(Center);
     1502    cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
    15001503    // Here we calculated center of circumscribing circle, using barycentric coordinates
     1504    cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
    15011505
    15021506    TempNormal.CopyVector(&a);
    … …  
    15221526    TempNormal.Normalize();
    15231527    Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
     1528    cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
    15241529    TempNormal.Scale(Restradius);
     1530    cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
    15251531
    15261532    Center->AddVector(&TempNormal);
     1533    cout << Verbose(4) << "Center of sphere of circumference is " << *Center << ".\n";
     1534    cout << Verbose(3) << "End of Get_center_of_sphere.\n";
    15271535  }
    15281536  ;
    … …  
    15571565    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
    15581566{
    1559   //cout << "ReferencePoint is " << ReferencePoint.x[0] << " "<< ReferencePoint.x[1] << " "<< ReferencePoint.x[2] << " "<< endl;
     1567  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
     1568  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
     1569  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
     1570  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
     1571  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
     1572  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
    15601573  /* OldNormal is normal vector on the old triangle
    15611574   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
    … …  
    15731586  atom *Walker; // variable atom point
    15741587
    1575 
    1576   dif_a.CopyVector(&(a->x));
    1577   dif_a.SubtractVector(&(Candidate->x));
    1578   dif_b.CopyVector(&(b->x));
    1579   dif_b.SubtractVector(&(Candidate->x));
    1580   AngleCheck.CopyVector(&(Candidate->x));
    1581   AngleCheck.SubtractVector(&(a->x));
    1582   AngleCheck.ProjectOntoPlane(Chord);
    1583 
    1584   SideA = dif_b.Norm();
    1585   SideB = dif_a.Norm();
    1586   SideC = Chord->Norm();
    1587   //Chord->Scale(-1);
    1588 
    1589   alpha = Chord->Angle(&dif_a);
    1590   beta = M_PI - Chord->Angle(&dif_b);
    1591   gamma = dif_a.Angle(&dif_b);
     1588  Vector NewUmkreismittelpunkt;
    15921589
    15931590
    15941591  if (a != Candidate and b != Candidate and c != Candidate)
    15951592    {
     1593      cout << Verbose(3) << "We have a unique candidate!" << endl;
     1594      dif_a.CopyVector(&(a->x));
     1595      dif_a.SubtractVector(&(Candidate->x));
     1596      dif_b.CopyVector(&(b->x));
     1597      dif_b.SubtractVector(&(Candidate->x));
     1598      AngleCheck.CopyVector(&(Candidate->x));
     1599      AngleCheck.SubtractVector(&(a->x));
     1600      AngleCheck.ProjectOntoPlane(Chord);
     1601
     1602      SideA = dif_b.Norm();
     1603      SideB = dif_a.Norm();
     1604      SideC = Chord->Norm();
     1605      //Chord->Scale(-1);
     1606
     1607      alpha = Chord->Angle(&dif_a);
     1608      beta = M_PI - Chord->Angle(&dif_b);
     1609      gamma = dif_a.Angle(&dif_b);
     1610
     1611      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;
     1612
     1613      if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
     1614        cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
    15961615
    15971616      Umkreisradius = SideA / 2.0 / sin(alpha);
    … …  
    16021621      if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
    16031622        {
    1604           cout << Verbose(1) << "Candidate is "<< *Candidate << endl;
     1623          cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
     1624          cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
    16051625          sign = AngleCheck.ScalarProduct(direction1);
    16061626          if (fabs(sign)<MYEPSILON)
    … …  
    16201640            {
    16211641              sign /= fabs(sign);
     1642              cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
    16221643            }
    16231644
    1624 
    1625 
    1626           Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
     1645          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
    16271646
    16281647          AngleCheck.CopyVector(&ReferencePoint);
    … …  
    16311650          AngleCheck.AddVector(&Mittelpunkt);
    16321651          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
     1652          cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
    16331653
    16341654          BallAngle = AngleCheck.Angle(OldNormal);
     1655          cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
    16351656
    16361657          //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
    16371658          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
    16381659
    1639           cout << Verbose(1) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
    1640 
    1641           if (AngleCheck.ScalarProduct(direction1) >=0)
     1660          cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
     1661
     1662          NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
     1663
     1664          if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON))
    16421665            {
    16431666              if (Storage[0]< -1.5) // first Candidate at all
    16441667                {
    16451668
    1646                   cout << "Next better candidate is " << *Candidate << " with ";
     1669                  cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
    16471670                  Opt_Candidate = Candidate;
    16481671                  Storage[0] = sign;
    16491672                  Storage[1] = AlternativeSign;
    16501673                  Storage[2] = BallAngle;
    1651                   cout << "Angle is " << Storage[2] << ", Halbraum ist "
     1674                  cout << " angle " << Storage[2] << " and Up/Down "
    16521675                    << Storage[0] << endl;
    1653 
    1654 
    16551676                }
    16561677              else
    … …  
    16581679                  if ( Storage[2] > BallAngle)
    16591680                    {
    1660                       cout << "Next better candidate is " << *Candidate << " with ";
     1681                      cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
    16611682                      Opt_Candidate = Candidate;
    16621683                      Storage[0] = sign;
    16631684                      Storage[1] = AlternativeSign;
    16641685                      Storage[2] = BallAngle;
    1665                       cout << "Angle is " << Storage[2] << ", Halbraum ist "
     1686                      cout << " angle " << Storage[2] << " and Up/Down "
    16661687                        << Storage[0] << endl;
    16671688                    }
    16681689                  else
    16691690                    {
    1670                       //if (DEBUG)
    1671                         cout << "Looses to better candidate" << endl;
    1672 
     1691                      if (DEBUG)
     1692                        {
     1693                          cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
     1694                        }
    16731695                    }
    16741696                }
    … …  
    16761698            else
    16771699              {
    1678                 //if (DEBUG)
    1679                   cout << "Refused due to bad direction of ball centre." << endl;
     1700              if (DEBUG)
     1701                {
     1702                  cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
     1703                }
    16801704              }
    16811705          }
    16821706        else
    16831707          {
    1684             //if (DEBUG)
    1685               cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
     1708            if (DEBUG)
     1709              {
     1710                cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
     1711              }
    16861712          }
    16871713      }
    16881714    else
    16891715      {
    1690         //if (DEBUG)
    1691           cout << "identisch mit Ursprungslinie" << endl;
    1692 
     1716        if (DEBUG)
     1717          {
     1718            cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
     1719          }
    16931720      }
    16941721
    … …  
    17131740        }
    17141741    }
     1742  cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
    17151743}
    17161744;
    … …  
    18341862          d1->ProjectOntoPlane(&AngleCheckReference);
    18351863          sign = AngleCheck.ScalarProduct(d1);
    1836           sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
     1864          if (fabs(sign) < MYEPSILON)
     1865            sign = 0;
     1866          else
     1867            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
    18371868
    18381869
    … …  
    19912022    const double& RADIUS, int N, const char *tempbasename)
    19922023{
    1993   cout << Verbose(1) << "Looking for next suitable triangle \n";
     2024  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
    19942025  Vector direction1;
    19952026  Vector helper;
    … …  
    20082039  Vector Opt_Mittelpunkt;
    20092040
    2010   cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
     2041  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
     2042
    20112043  helper.CopyVector(&(Line.endpoints[0]->node->x));
    20122044  for (int i = 0; i < 3; i++)
    … …  
    20292061      direction1.Scale(-1);
    20302062    }
     2063  cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
    20312064
    20322065  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
    20332066  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
    2034 
    2035 
    2036     Vector Umkreismittelpunkt, a, b, c;
    2037     double alpha, beta, gamma;
    2038     a.CopyVector(&(T.endpoints[0]->node->x));
    2039     b.CopyVector(&(T.endpoints[1]->node->x));
    2040     c.CopyVector(&(T.endpoints[2]->node->x));
    2041     a.SubtractVector(&(T.endpoints[1]->node->x));
    2042     b.SubtractVector(&(T.endpoints[2]->node->x));
    2043     c.SubtractVector(&(T.endpoints[0]->node->x));
    2044 
     2067  cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
     2068
     2069
     2070  Vector Umkreismittelpunkt, a, b, c;
     2071  double alpha, beta, gamma;
     2072  a.CopyVector(&(T.endpoints[0]->node->x));
     2073  b.CopyVector(&(T.endpoints[1]->node->x));
     2074  c.CopyVector(&(T.endpoints[2]->node->x));
     2075  a.SubtractVector(&(T.endpoints[1]->node->x));
     2076  b.SubtractVector(&(T.endpoints[2]->node->x));
     2077  c.SubtractVector(&(T.endpoints[0]->node->x));
     2078
     2079//  alpha = a.Angle(&c) - M_PI/2.;
     2080//  beta = b.Angle(&a);
     2081//  gamma = c.Angle(&b) - M_PI/2.;
    20452082    alpha = M_PI - a.Angle(&c);
    20462083    beta = M_PI - b.Angle(&a);
    20472084    gamma = M_PI - c.Angle(&b);
     2085    if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
     2086      cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
    20482087
    20492088    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) ;
    … …  
    20542093    cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
    20552094
     2095  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;
     2096  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) ;
     2097  Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
     2098  cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
     2099  if (DEBUG)
     2100    cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
     2101  tmp = 0;
     2102  for (int i=0;i<NDIM;i++) {
     2103    helper.CopyVector(&T.endpoints[i]->node->x);
     2104    helper.SubtractVector(&Umkreismittelpunkt);
     2105    if (DEBUG)
     2106      cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
     2107    if (tmp == 0) // set on first time for comparison against next ones
     2108      tmp = helper.Norm();
     2109    if (fabs(helper.Norm() - tmp) > MYEPSILON)
     2110      cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
     2111  }
    20562112
    20572113  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
    … …  
    20972153      // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
    20982154
    2099   cout << " Optimal candidate is " << *Opt_Candidate << endl;
     2155  cout << Verbose(2) << " Optimal candidate is " << *Opt_Candidate << endl;
     2156
     2157  // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
    21002158
    21012159  AddTrianglePoint(Opt_Candidate, 0);
    2102   AddTrianglePoint(Line.endpoints[0]->node, 1);
    2103   AddTrianglePoint(Line.endpoints[1]->node, 2);
    2104 
    2105   AddTriangleLine(TPS[0], TPS[1], 0);
    2106   AddTriangleLine(TPS[0], TPS[2], 1);
    2107   AddTriangleLine(TPS[1], TPS[2], 2);
    2108 
    2109   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    2110   AddTriangleToLines();
    2111   cout << "New triangle with " << *BTS << endl;
    2112   cout << "We have "<< TrianglesOnBoundaryCount << endl;
    2113   cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
    2114 
    2115   BTS->GetNormalVector(BTS->NormalVector);
    2116 
    2117   if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
    2118       (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) ||
    2119       (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) )
    2120 
    2121     {
    2122       BTS->NormalVector.Scale(-1);
    2123     };
    2124 
     2160  LineMap::iterator TryAndError;
     2161  TryAndError = TPS[0]->lines.find(Line.endpoints[0]->node->nr);
     2162  bool flag = true;
     2163  if ((TryAndError != TPS[0]->lines.end()) && (TryAndError->second->TrianglesCount > 1))
     2164    flag = false;
     2165  TryAndError = TPS[0]->lines.find(Line.endpoints[1]->node->nr);
     2166  if ((TryAndError !=  TPS[0]->lines.end()) && (TryAndError->second->TrianglesCount > 1))
     2167    flag = false;
     2168
     2169  if (flag) { // if so, add
     2170    AddTrianglePoint(Line.endpoints[0]->node, 1);
     2171    AddTrianglePoint(Line.endpoints[1]->node, 2);
     2172
     2173    AddTriangleLine(TPS[0], TPS[1], 0);
     2174    AddTriangleLine(TPS[0], TPS[2], 1);
     2175    AddTriangleLine(TPS[1], TPS[2], 2);
     2176
     2177    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     2178    AddTriangleToLines();
     2179
     2180    BTS->GetNormalVector(BTS->NormalVector);
     2181
     2182    if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
     2183        (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) ||
     2184        (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) )
     2185
     2186      {
     2187        BTS->NormalVector.Scale(-1);
     2188      };
     2189    cout << Verbose(2) << "New triangle with " << *BTS << "and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
     2190    cout << Verbose(2) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
     2191  } else { // else, yell and do nothing
     2192    cout << Verbose(2) << "This triangle consisting of ";
     2193    for (int i=0;i<NDIM;i++)
     2194      cout << BLS[i] << "(" << BLS[i]->TrianglesCount << "), ";
     2195    cout << " is invalid!" << endl;
     2196  }
     2197  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
    21252198}
    21262199;
    … …  
    21302203    molecule* mol, double RADIUS)
    21312204{
    2132   cout << Verbose(1)
    2133       << "Looking for second point of starting triangle, recursive level "
    2134       << RecursionLevel << endl;;
     2205  cout << Verbose(2)
     2206      << "Begin of Find_second_point_for_Tesselation, recursive level "
     2207      << RecursionLevel << endl;
    21352208  int i;
    21362209  Vector AngleCheck;
    21372210  atom* Walker;
    2138   double norm = -1.;
     2211  double norm = -1., angle;
    21392212
    21402213  // check if we only have one unique point yet ...
    21412214  if (a != Candidate)
    21422215    {
     2216    cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
    21432217      AngleCheck.CopyVector(&(Candidate->x));
    21442218      AngleCheck.SubtractVector(&(a->x));
    … …  
    21472221      if (norm < RADIUS)
    21482222        {
    2149           if (AngleCheck.Angle(&Oben) < Storage[0])
     2223        angle = AngleCheck.Angle(&Oben);
     2224          if (angle < Storage[0])
    21502225            {
    2151               //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[2]);
    2152               cout << "Next better candidate is " << *Candidate
    2153                   << " with distance " << norm << ".\n";
     2226              //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
     2227              cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
    21542228              Opt_Candidate = Candidate;
    21552229              Storage[0] = AngleCheck.Angle(&Oben);
    … …  
    21582232          else
    21592233            {
    2160               cout << Verbose(1) << "Supposedly looses to a better candidate "
     2234              cout << "Looses with angle " << angle << " to a better candidate "
    21612235                  << *Opt_Candidate << endl;
    21622236            }
    … …  
    21642238      else
    21652239        {
    2166           cout << Verbose(1) << *Candidate << " refused due to Radius " << norm
     2240          cout << "Refused due to Radius " << norm
    21672241              << endl;
    21682242        }
    … …  
    21832257        };
    21842258    };
     2259  cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level "
     2260      << RecursionLevel << endl;
    21852261}
    21862262;
    … …  
    21882264void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
    21892265{
    2190   cout << Verbose(1) << "Looking for starting triangle \n";
     2266  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
    21912267  int i = 0;
    21922268  atom* Walker;
    21932269  atom* FirstPoint;
    21942270  atom* SecondPoint;
    2195   atom* max_index[3];
    2196   double max_coordinate[3];
     2271  atom* max_index[NDIM];
     2272  double max_coordinate[NDIM];
    21972273  Vector Oben;
    21982274  Vector helper;
    … …  
    22072283      max_coordinate[i] = -1;
    22082284    }
    2209   cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount
     2285  cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount
    22102286      << " Atoms \n";
    22112287
    … …  
    22252301    }
    22262302
    2227   cout << Verbose(1) << "Found maximum coordinates. " << endl;
     2303  cout << Verbose(2) << "Found maximum coordinates: ";
     2304  for (int i=0;i<NDIM;i++)
     2305    cout << i << ": " << *max_index[i] << "\t";
     2306  cout << endl;
    22282307  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
    22292308  const int k = 1;
    … …  
    22312310  FirstPoint = max_index[k];
    22322311
    2233   cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << ": "
    2234       << FirstPoint->x.x[0] << endl;
     2312  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
    22352313  double Storage[3];
    22362314  atom* Opt_Candidate = NULL;
    … …  
    22382316  Storage[1] = 999999.; // This will be an angle looking for the third point.
    22392317  Storage[2] = 999999.;
    2240   cout << Verbose(1) << "Number of Bonds: "
    2241       << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
    22422318
    22432319  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
    22442320      Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
    22452321  SecondPoint = Opt_Candidate;
    2246   cout << Verbose(1) << "Found second point is " << *SecondPoint << ".\n";
     2322  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
    22472323
    22482324  helper.CopyVector(&(FirstPoint->x));
    … …  
    22602336  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
    22612337
    2262   cout << Verbose(1) << "Looking for third point candidates \n";
    2263   cout << Verbose(1) << " In direction " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << " " << endl;
     2338  cout << Verbose(2) << "Looking for third point candidates \n";
    22642339  // look in one direction of baseline for initial candidate
    22652340  Opt_Candidate = NULL;
    … …  
    22682343  CenterOfFirstLine.AddVector(&(SecondPoint->x));
    22692344
     2345  cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
    22702346  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0,
    22712347      &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
    22722348  // look in other direction of baseline for possible better candidate
     2349  cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
    22732350  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0,
    22742351      &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
    … …  
    22762353
    22772354  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
    2278 
    2279   cout << Verbose(1) << "The found starting triangle consists of "
    2280       << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
    2281       << "." << endl;
    22822355
    22832356  // Finally, we only have to add the found points
    … …  
    22952368  Oben.Scale(-1.);
    22962369  BTS->GetNormalVector(Oben);
     2370  cout << Verbose(0) << "==> The found starting triangle consists of "
     2371      << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
     2372      << " with normal vector " << BTS->NormalVector << ".\n";
     2373  cout << Verbose(1) << "End of Find_starting_triangle\n";
    22972374}
    22982375;
    … …  
    23112388
    23122389  baseline = Tess->LinesOnBoundary.begin();
    2313   while (baseline != Tess->LinesOnBoundary.end())
     2390  while ((baseline != Tess->LinesOnBoundary.end()) || (flag))
    23142391    {
    23152392      if (baseline->second->TrianglesCount == 1)
    23162393        {
    2317           cout << Verbose(1) << "Begin of Tesselation ... " << endl;
    23182394          Tess->Find_next_suitable_triangle(out, mol,
    23192395              *(baseline->second),
    23202396              *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
    23212397          flag = true;
    2322           cout << Verbose(1) << "End of Tesselation ... " << endl;
    23232398        }
    23242399      else
    23252400        {
    2326           cout << Verbose(1) << "There is a line with "
     2401          cout << Verbose(1) << "Line " << baseline->second << " has "
    23272402              << baseline->second->TrianglesCount << " triangles adjacent"
    23282403              << endl;
    … …  
    23302405      N++;
    23312406      baseline++;
     2407      if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
     2408        baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
     2409        flag = false;
     2410      }
    23322411    }
    23332412  cout << Verbose(1) << "Writing final tecplot file\n";
    … …  
    23462425    raster.close();
    23472426  }
    2348 }
    2349 ;
     2427
     2428  cout << Verbose(0) << "End of Find_non_convex_border\n";
     2429}
     2430;
     2431

  • src/vector.cpp

    r6d4a76 r12298c  
    219219 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
    220220 */
    221 double Vector::Angle(Vector *y) const
     221double Vector::Angle(const Vector *y) const
    222222{
    223223  return acos(this->ScalarProduct(y)/Norm()/y->Norm());
    … …  
    313313};
    314314
    315 ofstream& operator<<(ofstream& ost,Vector& m)
    316 {
    317         m.Output(&ost);
     315/** Prints a 3dim vector to a stream.
     316 * \param ost output stream
     317 * \param v Vector to be printed
     318 * \return output stream
     319 */
     320ostream& operator<<(ostream& ost,Vector& m)
     321{
     322  ost << "(";
     323  for (int i=0;i<NDIM;i++) {
     324    ost << m.x[i];
     325    if (i != 2)
     326      ost << ",";
     327  }
     328  ost << ")";
    318329        return ost;
    319330};

  • src/vector.hpp

    r6d4a76 r12298c  
    2121  double Projection(const Vector *y) const;
    2222  double Norm() const ;
    23   double Angle(Vector *y) const;
     23  double Angle(const Vector *y) const;
    2424
    2525  void AddVector(const Vector *y);
    … …  
    5555};
    5656
    57 ofstream& operator<<(ofstream& ost, Vector& m);
     57ostream & operator << (ostream& ost, Vector &m);
    5858Vector& operator+=(Vector& a, const Vector& b);
    5959Vector& operator*=(Vector& a, const double m);
Note: See TracChangeset for help on using the changeset viewer.