Changeset 6d4a76 for src/boundary.cpp


Ignore:
Timestamp:
Dec 29, 2008, 12:29:21 PM (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:
e1bc68
Parents:
12298c (diff), 23e09b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'ConcaveHull' of ​ssh://stud64dc-01/home/neuen/workspace/ESPACK into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp

Fixed the incorporation of write_raster3d_file() ...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    r12298c r6d4a76  
    14141414  LineMap::iterator LineWalker;
    14151415  //cout << "Manually checking endpoints for line." << endl;
    1416   if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr)
     1416  if ((a->lines.find(b->node->nr))->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
    14821481 * @param Direction vector indicates up/down
    14831482 * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
    … …  
    14911490 */
    14921491
    1493   void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
     1492  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, Vector* AlternativeDirection,
    14941493      double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
    14951494  {
    14961495    Vector TempNormal, helper;
    14971496    double Restradius;
    1498     cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
     1497
    14991498    *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
    15001499    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";
    15031500    // Here we calculated center of circumscribing circle, using barycentric coordinates
    1504     cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
    15051501
    15061502    TempNormal.CopyVector(&a);
    … …  
    15261522    TempNormal.Normalize();
    15271523    Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
    1528     cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
    15291524    TempNormal.Scale(Restradius);
    1530     cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
    15311525
    15321526    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";
    15351527  }
    15361528  ;
    … …  
    15651557    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
    15661558{
    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;
     1559  //cout << "ReferencePoint is " << ReferencePoint.x[0] << " "<< ReferencePoint.x[1] << " "<< ReferencePoint.x[2] << " "<< endl;
    15731560  /* OldNormal is normal vector on the old triangle
    15741561   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
    … …  
    15861573  atom *Walker; // variable atom point
    15871574
    1588   Vector NewUmkreismittelpunkt;
     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);
    15891592
    15901593
    15911594  if (a != Candidate and b != Candidate and c != Candidate)
    15921595    {
    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";
    16151596
    16161597      Umkreisradius = SideA / 2.0 / sin(alpha);
    … …  
    16211602      if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
    16221603        {
    1623           cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
    1624           cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
     1604          cout << Verbose(1) << "Candidate is "<< *Candidate << endl;
    16251605          sign = AngleCheck.ScalarProduct(direction1);
    16261606          if (fabs(sign)<MYEPSILON)
    … …  
    16401620            {
    16411621              sign /= fabs(sign);
    1642               cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
    16431622            }
    16441623
    1645           Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
     1624
     1625
     1626          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
    16461627
    16471628          AngleCheck.CopyVector(&ReferencePoint);
    … …  
    16501631          AngleCheck.AddVector(&Mittelpunkt);
    16511632          //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;
    16531633
    16541634          BallAngle = AngleCheck.Angle(OldNormal);
    1655           cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
    16561635
    16571636          //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
    16581637          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
    16591638
    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))
     1639          cout << Verbose(1) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
     1640
     1641          if (AngleCheck.ScalarProduct(direction1) >=0)
    16651642            {
    16661643              if (Storage[0]< -1.5) // first Candidate at all
    16671644                {
    16681645
    1669                   cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
     1646                  cout << "Next better candidate is " << *Candidate << " with ";
    16701647                  Opt_Candidate = Candidate;
    16711648                  Storage[0] = sign;
    16721649                  Storage[1] = AlternativeSign;
    16731650                  Storage[2] = BallAngle;
    1674                   cout << " angle " << Storage[2] << " and Up/Down "
     1651                  cout << "Angle is " << Storage[2] << ", Halbraum ist "
    16751652                    << Storage[0] << endl;
     1653
     1654
    16761655                }
    16771656              else
    … …  
    16791658                  if ( Storage[2] > BallAngle)
    16801659                    {
    1681                       cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
     1660                      cout << "Next better candidate is " << *Candidate << " with ";
    16821661                      Opt_Candidate = Candidate;
    16831662                      Storage[0] = sign;
    16841663                      Storage[1] = AlternativeSign;
    16851664                      Storage[2] = BallAngle;
    1686                       cout << " angle " << Storage[2] << " and Up/Down "
     1665                      cout << "Angle is " << Storage[2] << ", Halbraum ist "
    16871666                        << Storage[0] << endl;
    16881667                    }
    16891668                  else
    16901669                    {
    1691                       if (DEBUG)
    1692                         {
    1693                           cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
    1694                         }
     1670                      //if (DEBUG)
     1671                        cout << "Looses to better candidate" << endl;
     1672
    16951673                    }
    16961674                }
    … …  
    16981676            else
    16991677              {
    1700               if (DEBUG)
    1701                 {
    1702                   cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
    1703                 }
     1678                //if (DEBUG)
     1679                  cout << "Refused due to bad direction of ball centre." << endl;
    17041680              }
    17051681          }
    17061682        else
    17071683          {
    1708             if (DEBUG)
    1709               {
    1710                 cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
    1711               }
     1684            //if (DEBUG)
     1685              cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
    17121686          }
    17131687      }
    17141688    else
    17151689      {
    1716         if (DEBUG)
    1717           {
    1718             cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
    1719           }
     1690        //if (DEBUG)
     1691          cout << "identisch mit Ursprungslinie" << endl;
     1692
    17201693      }
    17211694
    … …  
    17401713        }
    17411714    }
    1742   cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
    17431715}
    17441716;
    … …  
    18621834          d1->ProjectOntoPlane(&AngleCheckReference);
    18631835          sign = AngleCheck.ScalarProduct(d1);
    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...
     1836          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
    18681837
    18691838
    … …  
    20221991    const double& RADIUS, int N, const char *tempbasename)
    20231992{
    2024   cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
     1993  cout << Verbose(1) << "Looking for next suitable triangle \n";
    20251994  Vector direction1;
    20261995  Vector helper;
    … …  
    20392008  Vector Opt_Mittelpunkt;
    20402009
    2041   cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
    2042 
     2010  cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
    20432011  helper.CopyVector(&(Line.endpoints[0]->node->x));
    20442012  for (int i = 0; i < 3; i++)
    … …  
    20612029      direction1.Scale(-1);
    20622030    }
    2063   cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
    20642031
    20652032  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
    20662033  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
    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.;
     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
    20822045    alpha = M_PI - a.Angle(&c);
    20832046    beta = M_PI - b.Angle(&a);
    20842047    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";
    20872048
    20882049    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) ;
    … …  
    20932054    cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
    20942055
    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   }
    21122056
    21132057  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
    … …  
    21532097      // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
    21542098
    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)
     2099  cout << " Optimal candidate is " << *Opt_Candidate << endl;
    21582100
    21592101  AddTrianglePoint(Opt_Candidate, 0);
    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";
     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
    21982125}
    21992126;
    … …  
    22032130    molecule* mol, double RADIUS)
    22042131{
    2205   cout << Verbose(2)
    2206       << "Begin of Find_second_point_for_Tesselation, recursive level "
    2207       << RecursionLevel << endl;
     2132  cout << Verbose(1)
     2133      << "Looking for second point of starting triangle, recursive level "
     2134      << RecursionLevel << endl;;
    22082135  int i;
    22092136  Vector AngleCheck;
    22102137  atom* Walker;
    2211   double norm = -1., angle;
     2138  double norm = -1.;
    22122139
    22132140  // check if we only have one unique point yet ...
    22142141  if (a != Candidate)
    22152142    {
    2216     cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
    22172143      AngleCheck.CopyVector(&(Candidate->x));
    22182144      AngleCheck.SubtractVector(&(a->x));
    … …  
    22212147      if (norm < RADIUS)
    22222148        {
    2223         angle = AngleCheck.Angle(&Oben);
    2224           if (angle < Storage[0])
     2149          if (AngleCheck.Angle(&Oben) < Storage[0])
    22252150            {
    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";
     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";
    22282154              Opt_Candidate = Candidate;
    22292155              Storage[0] = AngleCheck.Angle(&Oben);
    … …  
    22322158          else
    22332159            {
    2234               cout << "Looses with angle " << angle << " to a better candidate "
     2160              cout << Verbose(1) << "Supposedly looses to a better candidate "
    22352161                  << *Opt_Candidate << endl;
    22362162            }
    … …  
    22382164      else
    22392165        {
    2240           cout << "Refused due to Radius " << norm
     2166          cout << Verbose(1) << *Candidate << " refused due to Radius " << norm
    22412167              << endl;
    22422168        }
    … …  
    22572183        };
    22582184    };
    2259   cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level "
    2260       << RecursionLevel << endl;
    22612185}
    22622186;
    … …  
    22642188void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
    22652189{
    2266   cout << Verbose(1) << "Begin of Find_starting_triangle\n";
     2190  cout << Verbose(1) << "Looking for starting triangle \n";
    22672191  int i = 0;
    22682192  atom* Walker;
    22692193  atom* FirstPoint;
    22702194  atom* SecondPoint;
    2271   atom* max_index[NDIM];
    2272   double max_coordinate[NDIM];
     2195  atom* max_index[3];
     2196  double max_coordinate[3];
    22732197  Vector Oben;
    22742198  Vector helper;
    … …  
    22832207      max_coordinate[i] = -1;
    22842208    }
    2285   cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount
     2209  cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount
    22862210      << " Atoms \n";
    22872211
    … …  
    23012225    }
    23022226
    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;
     2227  cout << Verbose(1) << "Found maximum coordinates. " << endl;
    23072228  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
    23082229  const int k = 1;
    … …  
    23102231  FirstPoint = max_index[k];
    23112232
    2312   cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
     2233  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << ": "
     2234      << FirstPoint->x.x[0] << endl;
    23132235  double Storage[3];
    23142236  atom* Opt_Candidate = NULL;
    … …  
    23162238  Storage[1] = 999999.; // This will be an angle looking for the third point.
    23172239  Storage[2] = 999999.;
     2240  cout << Verbose(1) << "Number of Bonds: "
     2241      << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
    23182242
    23192243  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
    23202244      Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
    23212245  SecondPoint = Opt_Candidate;
    2322   cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
     2246  cout << Verbose(1) << "Found second point is " << *SecondPoint << ".\n";
    23232247
    23242248  helper.CopyVector(&(FirstPoint->x));
    … …  
    23362260  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
    23372261
    2338   cout << Verbose(2) << "Looking for third point candidates \n";
     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;
    23392264  // look in one direction of baseline for initial candidate
    23402265  Opt_Candidate = NULL;
    … …  
    23432268  CenterOfFirstLine.AddVector(&(SecondPoint->x));
    23442269
    2345   cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
    23462270  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0,
    23472271      &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
    23482272  // look in other direction of baseline for possible better candidate
    2349   cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
    23502273  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0,
    23512274      &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
    … …  
    23532276
    23542277  // 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;
    23552282
    23562283  // Finally, we only have to add the found points
    … …  
    23682295  Oben.Scale(-1.);
    23692296  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";
    23742297}
    23752298;
    … …  
    23882311
    23892312  baseline = Tess->LinesOnBoundary.begin();
    2390   while ((baseline != Tess->LinesOnBoundary.end()) || (flag))
     2313  while (baseline != Tess->LinesOnBoundary.end())
    23912314    {
    23922315      if (baseline->second->TrianglesCount == 1)
    23932316        {
     2317          cout << Verbose(1) << "Begin of Tesselation ... " << endl;
    23942318          Tess->Find_next_suitable_triangle(out, mol,
    23952319              *(baseline->second),
    23962320              *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
    23972321          flag = true;
     2322          cout << Verbose(1) << "End of Tesselation ... " << endl;
    23982323        }
    23992324      else
    24002325        {
    2401           cout << Verbose(1) << "Line " << baseline->second << " has "
     2326          cout << Verbose(1) << "There is a line with "
    24022327              << baseline->second->TrianglesCount << " triangles adjacent"
    24032328              << endl;
    … …  
    24052330      N++;
    24062331      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       }
    24112332    }
    24122333  cout << Verbose(1) << "Writing final tecplot file\n";
    … …  
    24252346    raster.close();
    24262347  }
    2427 
    2428   cout << Verbose(0) << "End of Find_non_convex_border\n";
    2429 }
    2430 ;
    2431 
     2348}
     2349;
Note: See TracChangeset for help on using the changeset viewer.