Changeset 5c7bf8 for src/tesselation.cpp


Ignore:
Timestamp:
Aug 8, 2009, 7:25:28 PM (15 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:
08ef35
Parents:
8bb475f
Message:

BUGFIXes to Find_convex_... and Find_non_convex_... algorithms

  • Find_non_convex_border() - very first baseline of starting triangle is skipped as it is the most outward.
  • BoundaryLineSet::TrianglesCount is removed and replaced by triangles.size(), as it is redundant and we made mistakes in not reducing TrianglesCount on triangle erase.
  • BoundaryLineSet::CheckConvexityCriterion() - extra check when NormalVector point in the very same directions,
  • BoundaryTriangleSet::GetIntersectionInsideTriangle() - CrossPoint return value is checked, i===3 checked not >3!
  • new stream operator for TesselPoint with single argument.
  • class Tesselation inherits PointCloud and implements all functions (including virtual destructor) in order to have PointsOnBoundary parsable by LinkedCell.
  • Tesselation::InsertStraddlingPoints() - LinkedCell is created for every atom just from the BoundaryPoints and not for all (so that only BoundaryPoints can be found), old triangle is now really deleted, lots of verbosity added and fixed.
  • Tesselation::GetCommonEndpoint() - connectedClosestPoint is deleted lateron, SecondPoint argument but not used.
  • Tesselation::FindClosestPoint() - returns second second closest point in new second argument.
  • Tesselation::getCircleOfConnectedPoints() - we find the BoundaryPoint if possible and just look through its lines if found, plane is centered at closest point and normal is between it and the center of all on circle, getAngle() corrections incorporated and vector is always projected onto circle plane.
  • new function Tesselation findSecondClosestPoint().
  • getAngle() - just receives third vector with which it dedices between [0,pi) and [pi,2pi), the vectors are assumed to be relative to one common center already.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselation.cpp

    r8bb475f r5c7bf8  
    6464  for (int i = 0; i < 2; i++)
    6565    endpoints[i] = NULL;
    66   TrianglesCount = 0;
    6766  Nr = -1;
    6867}
    … …  
    7978  Point[1]->AddLine(this); //
    8079  // clear triangles list
    81   TrianglesCount = 0;
    8280  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
    8381}
    … …  
    119117      << endl;
    120118  triangles.insert(TrianglePair(triangle->Nr, triangle));
    121   TrianglesCount++;
    122119}
    123120;
    … …  
    143140bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
    144141{
    145   Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2];
     142  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
    146143  // get the two triangles
    147   if (TrianglesCount != 2) {
    148     *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
     144  if (triangles.size() != 2) {
     145    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
    149146    return false;
    150147  }
     148  // check normal vectors
    151149  // have a normal vector on the base line pointing outwards
    152150  *out << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
    … …  
    159157
    160158  BaseLineNormal.Zero();
     159  NormalCheck.Zero();
     160  double sign = -1.;
    161161  int i=0;
    162162  class BoundaryPointSet *node = NULL;
    163163  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
    164164    *out << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     165    NormalCheck.AddVector(&runner->second->NormalVector);
     166    NormalCheck.Scale(sign);
     167    sign = -sign;
    165168    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
    166169    node = runner->second->GetThirdEndpoint(this);
    … …  
    178181  }
    179182  *out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
    180   double angle = helper[0].Angle(&helper[1]);
    181   if (BaseLineNormal.ScalarProduct(&helper[1]) > 0) {
    182     angle = 2.*M_PI - angle;
    183   }
    184   *out << Verbose(3) << "The angle is " << angle << "." << endl;
     183  if (NormalCheck.NormSquared() < MYEPSILON) {
     184    *out << Verbose(3) << "Normalvectors of both triangles are the same: convex." << endl;
     185    return true;
     186  }
     187  double angle = getAngle(helper[0], helper[1], BaseLineNormal);
    185188  if ((angle - M_PI) > -MYEPSILON)
    186189    return true;
    … …  
    327330  Vector CrossPoint;
    328331  Vector helper;
     332
     333  if (!Intersection->GetIntersectionWithPlane(out, &NormalVector, endpoints[0]->node->node, MolCenter, x)) {
     334    *out << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;
     335    return false;
     336  }
     337
     338  // Calculate cross point between one baseline and the line from the third endpoint to intersection
    329339  int i=0;
    330 
    331   if (Intersection->GetIntersectionWithPlane(out, &NormalVector, endpoints[0]->node->node, MolCenter, x)) {
    332     *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
    333     return false;
    334   }
    335 
    336   // Calculate cross point between one baseline and the line from the third endpoint to intersection
    337340  do {
    338     CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection);
    339     helper.CopyVector(endpoints[(i+1)%3]->node->node);
    340     helper.SubtractVector(endpoints[i%3]->node->node);
    341     i++;
    342     if (i>3)
     341    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
     342      helper.CopyVector(endpoints[(i+1)%3]->node->node);
     343      helper.SubtractVector(endpoints[i%3]->node->node);
     344    } else
     345      i++;
     346    if (i>2)
    343347      break;
    344348  } while (CrossPoint.NormSquared() < MYEPSILON);
    345   if (i>3) {
     349  if (i==3) {
    346350    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
    347351    exit(255);
    … …  
    466470};
    467471
     472/** Prints LCNode to screen.
     473 */
     474ostream & TesselPoint::operator << (ostream &ost)
     475{
     476  ost << "[" << (Name) << "|" << this << "]";
     477  return ost;
     478};
     479
    468480
    469481// =========================================================== class POINTCLOUD ============================================
    … …  
    510522  LinesOnBoundaryCount = 0;
    511523  TrianglesOnBoundaryCount = 0;
     524  InternalPointer = PointsOnBoundary.begin();
    512525}
    513526;
    … …  
    528541}
    529542;
     543
     544/** PointCloud implementation of GetCenter
     545 * Uses PointsOnBoundary and STL stuff.
     546 */   
     547Vector * Tesselation::GetCenter(ofstream *out)
     548{
     549  Vector *Center = new Vector(0.,0.,0.);
     550  int num=0;
     551  for (GoToFirst(); (!IsEnd()); GoToNext()) {
     552    Center->AddVector(GetPoint()->node);
     553    num++;
     554  }
     555  Center->Scale(1./num);
     556  return Center;
     557};
     558
     559/** PointCloud implementation of GoPoint
     560 * Uses PointsOnBoundary and STL stuff.
     561 */   
     562TesselPoint * Tesselation::GetPoint()
     563{
     564  return (InternalPointer->second->node);
     565};
     566
     567/** PointCloud implementation of GetTerminalPoint.
     568 * Uses PointsOnBoundary and STL stuff.
     569 */   
     570TesselPoint * Tesselation::GetTerminalPoint()
     571{
     572  PointMap::iterator Runner = PointsOnBoundary.end();
     573  Runner--;
     574  return (Runner->second->node);
     575};
     576
     577/** PointCloud implementation of GoToNext.
     578 * Uses PointsOnBoundary and STL stuff.
     579 */   
     580void Tesselation::GoToNext()
     581{
     582  if (InternalPointer != PointsOnBoundary.end())
     583    InternalPointer++;
     584};
     585
     586/** PointCloud implementation of GoToPrevious.
     587 * Uses PointsOnBoundary and STL stuff.
     588 */   
     589void Tesselation::GoToPrevious()
     590{
     591  if (InternalPointer != PointsOnBoundary.begin())
     592    InternalPointer--;
     593};
     594
     595/** PointCloud implementation of GoToFirst.
     596 * Uses PointsOnBoundary and STL stuff.
     597 */   
     598void Tesselation::GoToFirst()
     599{
     600  InternalPointer = PointsOnBoundary.begin();
     601};
     602
     603/** PointCloud implementation of GoToLast.
     604 * Uses PointsOnBoundary and STL stuff.
     605 */   
     606void Tesselation::GoToLast()
     607{
     608  InternalPointer = PointsOnBoundary.end();
     609  InternalPointer--;
     610};
     611
     612/** PointCloud implementation of IsEmpty.
     613 * Uses PointsOnBoundary and STL stuff.
     614 */   
     615bool Tesselation::IsEmpty()
     616{
     617  return (PointsOnBoundary.empty());
     618};
     619
     620/** PointCloud implementation of IsLast.
     621 * Uses PointsOnBoundary and STL stuff.
     622 */   
     623bool Tesselation::IsEnd()
     624{
     625  return (InternalPointer == PointsOnBoundary.end());
     626};
     627
    530628
    531629/** Gueses first starting triangle of the convex envelope.
    … …  
    718816    flag = false;
    719817    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
    720       if (baseline->second->TrianglesCount == 1) {
     818      if (baseline->second->triangles.size() == 1) {
    721819        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
    722820        SmallestAngle = M_PI;
    … …  
    783881            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
    784882            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
    785             if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->TrianglesCount == 2))) {
    786               *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
     883            if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) {
     884              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl;
    787885              continue;
    788886            }
    789             if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
    790               *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
     887            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) {
     888              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl;
    791889              continue;
    792890            }
    … …  
    889987        // 5d. If the set of lines is not yet empty, go to 5. and continue
    890988      } else
    891         *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
     989        *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl;
    892990  } while (flag);
    893991
    … …  
    9041002bool Tesselation::InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC)
    9051003{
    906   Vector Intersection;
     1004  Vector Intersection, Normal;
    9071005  TesselPoint *Walker = NULL;
    9081006  Vector *Center = cloud->GetCenter(out);
    … …  
    9131011  cloud->GoToFirst();
    9141012  while (!cloud->IsEnd()) {  // we only have to go once through all points, as boundary can become only bigger
     1013    LinkedCell BoundaryPoints(this, 5.);
    9151014    Walker = cloud->GetPoint();
    9161015    *out << Verbose(2) << "Current point is " << *Walker << "." << endl;
    9171016    // get the next triangle
    918     triangles = FindClosestTrianglesToPoint(out, Walker->node, LC);
     1017    triangles = FindClosestTrianglesToPoint(out, Walker->node, &BoundaryPoints);
    9191018    if (triangles == NULL) {
    9201019      *out << Verbose(1) << "No triangles found, probably a tesselation point itself." << endl;
    … …  
    9241023      BTS = triangles->front();
    9251024    }
    926     *out << Verbose(2) << "Closest triangle is " << BTS << "." << endl;
     1025    *out << Verbose(2) << "Closest triangle is " << *BTS << "." << endl;
    9271026    // get the intersection point
    9281027    if (BTS->GetIntersectionInsideTriangle(out, Center, Walker->node, &Intersection)) {
    … …  
    9311030      if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) {
    9321031        // inside, next!
    933         *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
     1032        *out << Verbose(2) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;
    9341033      } else {
    9351034        // outside!
    936         *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
     1035        *out << Verbose(2) << *Walker << " is outside wrt triangle " << *BTS << "." << endl;
    9371036        class BoundaryLineSet *OldLines[3], *NewLines[3];
    9381037        class BoundaryPointSet *OldPoints[3], *NewPoint;
    … …  
    9421041          OldPoints[i] = BTS->endpoints[i];
    9431042        }
     1043        Normal.CopyVector(&BTS->NormalVector);
    9441044        // add Walker to boundary points
     1045        *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
    9451046        AddPoint(Walker);
    946         if (BPS[0] == NULL)
     1047        if (BPS[0] != NULL)
    9471048          NewPoint = BPS[0];
    9481049        else
    9491050          continue;
    9501051        // remove triangle
     1052        *out << Verbose(2) << "Erasing triangle " << *BTS << "." << endl;
    9511053        TrianglesOnBoundary.erase(BTS->Nr);
     1054        delete(BTS);
    9521055        // create three new boundary lines
    9531056        for (int i=0;i<3;i++) {
    … …  
    9551058          BPS[1] = OldPoints[i];
    9561059          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
     1060          *out << Verbose(3) << "Creating new line " << *NewLines[i] << "." << endl;
    9571061          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
    9581062          LinesOnBoundaryCount++;
    … …  
    9731077          // create the triangle
    9741078          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     1079          Normal.Scale(-1.);
     1080          BTS->GetNormalVector(Normal);
     1081          Normal.Scale(-1.);
     1082          *out << Verbose(2) << "Created new triangle " << *BTS << "." << endl;
    9751083          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    9761084          TrianglesOnBoundaryCount++;
    … …  
    10461154    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
    10471155      // If there is a line with less than two attached triangles, we don't need a new line.
    1048       if (FindLine->second->TrianglesCount < 2) {
     1156      if (FindLine->second->triangles.size() < 2) {
    10491157        insertNewLine = false;
    10501158        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
    … …  
    12471355    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
    12481356
    1249     // adding point 1 and point 2 and the line between them
     1357    // adding point 1 and point 2 and add the line between them
    12501358    AddTrianglePoint(FirstPoint, 0);
    12511359    AddTrianglePoint(SecondPoint, 1);
    … …  
    14641572          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
    14651573          << " for this triangle ... " << endl;
    1466           cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
     1574          cout << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl;
    14671575        } else {
    14681576          cout << Verbose(1) << "WARNING: This triangle consisting of ";
    … …  
    15131621  *out << Verbose(1) << "Begin of CorrectConcaveBaselines" << endl;
    15141622
    1515   for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
     1623  // go through every baseline
     1624  LineMap::iterator Advancer = LinesOnBoundary.begin();
     1625  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline = Advancer) {
     1626    Advancer++;  // as the base line might be removed, we rather have the next one already present
    15161627    Base = baseline->second;
    15171628    *out << Verbose(2) << "Current baseline is " << *Base << " ... " << endl;
     1629    // calculate NormalVector for later use
     1630    BaseLineNormal.Zero();
     1631    if (Base->triangles.size() < 2) {
     1632      *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
     1633      return false;
     1634    }
     1635    for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
     1636      *out << Verbose(4) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;
     1637      BaseLineNormal.AddVector(&(runner->second->NormalVector));
     1638    }
     1639    BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
    15181640    // check convexity
    15191641    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
    … …  
    15611683
    15621684      // remove triangles and baseline removes itself
    1563       m=0;
     1685      *out << Verbose(3) << "INFO: Deleting baseline " << *Base << "." << endl;
    15641686      OldBaseLine = Base->Nr;
    15651687      LinesOnBoundary.erase(OldBaseLine);
     1688      m=0;
    15661689      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
     1690        *out << Verbose(3) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
    15671691        TrianglesOnBoundary.erase(OldTriangles[m++] = runner->second->Nr);
    15681692        delete(runner->second);
    … …  
    15741698      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
    15751699      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
     1700      *out << Verbose(3) << "INFO: Created new baseline " << *NewLine << "." << endl;
    15761701
    15771702      // construct new triangles with flipped baseline
    … …  
    15861711        BLS[2] = NewLine;
    15871712        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
     1713        BTS->GetNormalVector(BaseLineNormal);
    15881714        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
     1715        *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
    15891716
    15901717        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
    … …  
    15921719        BLS[2] = NewLine;
    15931720        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
     1721        BTS->GetNormalVector(BaseLineNormal);
    15941722        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
     1723        *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
    15951724      } else {
    15961725        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
    … …  
    19492078list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC)
    19502079{
    1951   class TesselPoint *trianglePoints[3];
     2080  TesselPoint *trianglePoints[3];
     2081  TesselPoint *SecondPoint = NULL;
    19522082
    19532083  if (LinesOnBoundary.empty()) {
    1954     *out << Verbose(0) << "Error: There is no tesselation structure to compare the point with, " << "please create one first.";
     2084    *out << Verbose(0) << "Error: There is no tesselation structure to compare the point with, please create one first.";
    19552085    return NULL;
    19562086  }
    19572087
    1958   trianglePoints[0] = findClosestPoint(x, LC);
     2088  trianglePoints[0] = findClosestPoint(x, SecondPoint, LC);
     2089 
    19592090  // check whether closest point is "too close" :), then it's inside
     2091  if (trianglePoints[0] == NULL) {
     2092    *out << Verbose(1) << "Is the only point, no one else is closeby." << endl;
     2093    return NULL;
     2094  }
    19602095  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
    19612096    *out << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
    … …  
    19692104      *out << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
    19702105    }
    1971     *out << Verbose(1) << "List of possible points:" << endl;
    1972     *out << *trianglePoints[i] << endl;
    1973   }
     2106    //*out << Verbose(1) << "List of possible points:" << endl;
     2107    //*out << Verbose(2) << *trianglePoints[i] << endl;
     2108  }
     2109
     2110  list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);
     2111
    19742112  delete(connectedClosestPoints);
    1975 
    1976   list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);
    1977 
     2113 
    19782114  if (triangles->empty()) {
    19792115    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
    … …  
    20592195  map<double, TesselPoint*>::iterator runner;
    20602196  list<TesselPoint*>::iterator listRunner;
    2061   Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
     2197  class BoundaryPointSet *ReferencePoint = NULL;
     2198  Vector center, PlaneNormal, currentPoint, OrthogonalVector, helper, AngleZero;
    20622199  TesselPoint* current;
    20632200  bool takePoint = false;
    20642201
    2065   planeNorm.CopyVector(Point->node);
    2066   planeNorm.SubtractVector(Reference);
    2067   planeNorm.Normalize();
    2068 
    2069   LineMap::iterator findLines = LinesOnBoundary.begin();
    2070   while (findLines != LinesOnBoundary.end()) {
     2202  // find the respective boundary point
     2203  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
     2204  if (PointRunner != PointsOnBoundary.end()) {
     2205    ReferencePoint = PointRunner->second;
     2206  } else {
     2207    *out << Verbose(2) << "getCircleOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
     2208    ReferencePoint = NULL;
     2209  }
     2210
     2211  // little trick so that we look just through lines connect to the BoundaryPoint
     2212  // OR fall-back to look through all lines if there is no such BoundaryPoint
     2213  LineMap *Lines = &LinesOnBoundary;
     2214  if (ReferencePoint != NULL)
     2215    Lines = &(ReferencePoint->lines);
     2216  LineMap::iterator findLines = Lines->begin();
     2217  while (findLines != Lines->end()) {
    20712218    takePoint = false;
    20722219
    … …  
    20782225      current = findLines->second->endpoints[0]->node;
    20792226    }
    2080 
     2227   
    20812228    if (takePoint) {
     2229      *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
    20822230      connectedPoints.push_back(current);
    20832231      currentPoint.CopyVector(current->node);
    2084       currentPoint.ProjectOntoPlane(&planeNorm);
     2232      currentPoint.ProjectOntoPlane(&PlaneNormal);
    20852233      center.AddVector(&currentPoint);
    20862234    }
    … …  
    20892237  }
    20902238
    2091   *out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
    2092     << "; scale factor " << 1.0/connectedPoints.size();
    2093 
     2239  //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
     2240  //  << "; scale factor " << 1.0/connectedPoints.size();
     2241
     2242  if (connectedPoints.size() == 0) { // if have not found any points
     2243    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
     2244    return NULL;
     2245  }
    20942246  center.Scale(1.0/connectedPoints.size());
     2247
     2248  *out << Verbose(4) << "INFO: Calculated center of all circle points is " << center << "." << endl;
     2249
     2250  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
     2251  PlaneNormal.CopyVector(Point->node);
     2252  PlaneNormal.SubtractVector(&center);
     2253  PlaneNormal.Normalize();
     2254  *out << Verbose(4) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
     2255
     2256  // construct one orthogonal vector
     2257  AngleZero.CopyVector(Reference);
     2258  AngleZero.SubtractVector(Point->node);
     2259  AngleZero.ProjectOntoPlane(&PlaneNormal);
     2260  *out << Verbose(4) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
     2261  OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
     2262  *out << Verbose(4) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
     2263 
     2264  // go through all connected points and calculate angle
    20952265  listRunner = connectedPoints.begin();
    2096 
    2097   *out << " calculated center " << center <<  endl;
    2098 
    2099   // construct one orthogonal vector
    2100   helper.CopyVector(Reference);
    2101   helper.ProjectOntoPlane(&planeNorm);
    2102   OrthogonalVector.MakeNormalVector(&center, &helper, (*listRunner)->node);
    21032266  while (listRunner != connectedPoints.end()) {
    2104     double angle = getAngle(*((*listRunner)->node), *(Reference), center, OrthogonalVector);
    2105     *out << "Calculated angle " << angle << " for point " << **listRunner << endl;
     2267    helper.CopyVector((*listRunner)->node);
     2268    helper.SubtractVector(Point->node);
     2269    helper.ProjectOntoPlane(&PlaneNormal);
     2270    double angle = getAngle(helper, AngleZero, OrthogonalVector);
     2271    *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
    21062272    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
    21072273    listRunner++;
    … …  
    21102276  list<TesselPoint*> *result = new list<TesselPoint*>;
    21112277  runner = anglesOfPoints.begin();
    2112   *out << "First value is " << *runner->second << endl;
    21132278  result->push_back(runner->second);
    21142279  runner = anglesOfPoints.end();
    21152280  runner--;
    2116   *out << "Second value is " << *runner->second << endl;
    21172281  result->push_back(runner->second);
    21182282
    2119   *out << "List of closest points has " << result->size() << " elements, which are "
     2283  *out << Verbose(2) << "List of closest points has " << result->size() << " elements, which are "
    21202284    << *(result->front()) << " and " << *(result->back()) << endl;
    21212285
    … …  
    21432307        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
    21442308          // If there is a line with less than two attached triangles, we don't need a new line.
    2145           if (FindLine->second->TrianglesCount < 2) {
     2309          if (FindLine->second->triangles.size() < 2) {
    21462310            counter++;
    21472311            break;  // increase counter only once per edge
    … …  
    21492313        }
    21502314      } else { // no line
    2151         cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
     2315        cout << Verbose(1) << "The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
    21522316        result = true;
    21532317      }
    … …  
    22052369};
    22062370
    2207 
    2208 
    22092371/**
    2210  * Finds the point which is closest to the provided one.
     2372 * Finds the point which is second closest to the provided one.
    22112373 *
    2212  * @param Point to which to find the closest other point
     2374 * @param Point to which to find the second closest other point
    22132375 * @param linked cell structure
    22142376 *
    2215  * @return point which is closest to the provided one
    2216  */
    2217 TesselPoint* findClosestPoint(const Vector* Point, LinkedCell* LC)
     2377 * @return point which is second closest to the provided one
     2378 */
     2379TesselPoint* findSecondClosestPoint(const Vector* Point, LinkedCell* LC)
    22182380{
    22192381  LinkedNodes *List = NULL;
    22202382  TesselPoint* closestPoint = NULL;
     2383  TesselPoint* secondClosestPoint = NULL;
    22212384  double distance = 1e16;
     2385  double secondDistance = 1e16;
    22222386  Vector helper;
    22232387  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
    … …  
    22262390  for(int i=0;i<NDIM;i++) // store indices of this cell
    22272391    N[i] = LC->n[i];
    2228   //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     2392  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    22292393
    22302394  LC->GetNeighbourBounds(Nlower, Nupper);
    … …  
    22342398      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    22352399        List = LC->GetCurrentCell();
    2236         //cout << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     2400        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    22372401        if (List != NULL) {
    22382402          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    … …  
    22412405            double currentNorm = helper. Norm();
    22422406            if (currentNorm < distance) {
     2407              // remember second point
     2408              secondDistance = distance;
     2409              secondClosestPoint = closestPoint;
     2410              // mark down new closest point
    22432411              distance = currentNorm;
    22442412              closestPoint = (*Runner);
     2413              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    22452414            }
    22462415          }
    … …  
    22512420      }
    22522421
     2422  return secondClosestPoint;
     2423};
     2424
     2425
     2426
     2427/**
     2428 * Finds the point which is closest to the provided one.
     2429 *
     2430 * @param Point to which to find the closest other point
     2431 * @param SecondPoint the second closest other point on return, NULL if none found
     2432 * @param linked cell structure
     2433 *
     2434 * @return point which is closest to the provided one, NULL if none found
     2435 */
     2436TesselPoint* findClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
     2437{
     2438  LinkedNodes *List = NULL;
     2439  TesselPoint* closestPoint = NULL;
     2440  SecondPoint = NULL;
     2441  double distance = 1e16;
     2442  double secondDistance = 1e16;
     2443  Vector helper;
     2444  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
     2445
     2446  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
     2447  for(int i=0;i<NDIM;i++) // store indices of this cell
     2448    N[i] = LC->n[i];
     2449  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     2450
     2451  LC->GetNeighbourBounds(Nlower, Nupper);
     2452  //cout << endl;
     2453  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
     2454    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
     2455      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
     2456        List = LC->GetCurrentCell();
     2457        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     2458        if (List != NULL) {
     2459          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
     2460            helper.CopyVector(Point);
     2461            helper.SubtractVector((*Runner)->node);
     2462            double currentNorm = helper. Norm();
     2463            if (currentNorm < distance) {
     2464              secondDistance = distance;
     2465              SecondPoint = closestPoint;
     2466              distance = currentNorm;
     2467              closestPoint = (*Runner);
     2468              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
     2469            } else if (currentNorm < secondDistance) {
     2470              secondDistance = currentNorm;
     2471              SecondPoint = (*Runner);
     2472              cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
     2473            }
     2474          }
     2475        } else {
     2476          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
     2477            << LC->n[2] << " is invalid!" << endl;
     2478        }
     2479      }
     2480
    22532481  return closestPoint;
    2254 }
     2482};
    22552483
    22562484
    … …  
    23212549
    23222550/** Gets the angle between a point and a reference relative to the provided center.
    2323  * We have two shanks (point, center) and (reference, center) between which the angle is calculated
     2551 * We have two shanks point and reference between which the angle is calculated
    23242552 * and by scalar product with OrthogonalVector we decide the interval.
    23252553 * @param point to calculate the angle for
    23262554 * @param reference to which to calculate the angle
    2327  * @param center for which to calculate the angle between the vectors
    23282555 * @param OrthogonalVector points in direction of [pi,2pi] interval
    23292556 *
    23302557 * @return angle between point and reference
    23312558 */
    2332 double getAngle(const Vector &point, const Vector &reference, const Vector &center, Vector OrthogonalVector)
    2333 {
    2334   Vector ReferenceVector, helper;
    2335   cout << Verbose(4) << center << " is our center " << endl;
    2336   // create baseline vector
    2337   ReferenceVector.CopyVector(&reference);
    2338   ReferenceVector.SubtractVector(&center);
    2339   OrthogonalVector.MakeNormalVector(&ReferenceVector);
    2340   cout << Verbose(4) << ReferenceVector << " is our reference " << endl;
    2341   if (ReferenceVector.IsNull())
     2559double getAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
     2560{
     2561  if (reference.IsNull())
    23422562    return M_PI;
    23432563
    23442564  // calculate both angles and correct with in-plane vector
    2345   helper.CopyVector(&point);
    2346   helper.SubtractVector(&center);
    2347   if (helper.IsNull())
     2565  if (point.IsNull())
    23482566    return M_PI;
    2349   double phi = ReferenceVector.Angle(&helper);
    2350   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     2567  double phi = point.Angle(&reference);
     2568  if (OrthogonalVector.ScalarProduct(&point) > 0) {
    23512569    phi = 2.*M_PI - phi;
    23522570  }
    23532571
    2354   cout << Verbose(3) << helper << " has angle " << phi << " with respect to reference." << endl;
     2572  cout << Verbose(3) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
    23552573
    23562574  return phi;
Note: See TracChangeset for help on using the changeset viewer.