Changeset f1ef60a


Ignore:
Timestamp:
Nov 10, 2009, 8:59:00 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:
717e0c
Parents:
9a0dc8
Message:

Attempt to fix tesselation for case 1_2-dimethylbenzene.

  • The issue is still the strongly degenerated end of the molecule, that is very flat.
  • new function: Tesselation::HasOtherBaselineBetterCandidate() - checks whether for the two other baselines of a proposed triangle there isn't an even better candidate
  • Tesselation::FindNextSuitableTriangle() uses the above function.
  • FindNonConvexBorder() - instead of jumping to end() if a new line was attempted, we jump to the beginning but only if new line was created succesfully.
  • DoSingleStepOutput reactived for this debugging.
  • Where the newly installed tesselations did only fail for this molecule, now they fail for all and some segfault. Only this runs through cleanly.

Signed-off-by: Frederik Heber <heber@…>

Location:
src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    r9a0dc8 rf1ef60a  
    994994  while ((baseline != TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
    995995    if (baseline->second->triangles.size() == 1) {
     996      CheckListOfBaselines(TesselStruct);
    996997      // 3. find next triangle
    997998      TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
    998999      OneLoopWithoutSuccessFlag = OneLoopWithoutSuccessFlag || TesselationFailFlag;
    9991000      if (!TesselationFailFlag)
    1000         eLog() << Verbose(0) << "WARNING: FindNextSuitableTriangle failed." << endl;
     1001        eLog() << Verbose(2) << "FindNextSuitableTriangle failed." << endl;
    10011002
    10021003      // write temporary envelope
    … …  
    10061007        }
    10071008      }
    1008       baseline = TesselStruct->LinesOnBoundary.end();
    1009       Log() << Verbose(2) << "Baseline set to end." << endl;
     1009      if (TesselationFailFlag) {
     1010        baseline = TesselStruct->LinesOnBoundary.begin();
     1011        OneLoopWithoutSuccessFlag = false;
     1012        Log() << Verbose(2) << "Baseline set to begin." << endl;
     1013      }
    10101014    } else {
    10111015      //Log() << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
    10121016      if (baseline->second->triangles.size() != 2)
    1013         Log() << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
     1017        eLog() << Verbose(1) << "TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
    10141018    }
    10151019

  • src/boundary.hpp

    r9a0dc8 rf1ef60a  
    3535
    3636#define DEBUG 1
    37 #define DoSingleStepOutput 0
     37#define DoSingleStepOutput 1
    3838#define SingleStepWidth 1
    3939

  • src/tesselation.cpp

    r9a0dc8 rf1ef60a  
    12441244;
    12451245
     1246/** Sets point to a present Tesselation::PointsOnBoundary.
     1247 * Tesselation::TPS is set to the existing one or NULL if not found.
     1248 * @param Candidate point to set to
     1249 * @param n index for this point in Tesselation::TPS array
     1250 */
     1251void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const
     1252{
     1253  PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr);
     1254  if (FindPoint != PointsOnBoundary.end())
     1255    TPS[n] = FindPoint->second;
     1256  else
     1257    TPS[n] = NULL;
     1258};
     1259
    12461260/** Function tries to add line from current Points in BPS to BoundaryLineSet.
    12471261 * If successful it raises the line count and inserts the new line into the BLS,
    … …  
    14511465 *                 triangles exist which is the maximum for three points
    14521466 */
    1453 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) {
     1467int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const
     1468{
    14541469  int adjacentTriangleCount = 0;
    14551470  class BoundaryPointSet *Points[3];
    … …  
    14581473  // builds a triangle point set (Points) of the end points
    14591474  for (int i = 0; i < 3; i++) {
    1460     PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
     1475    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
    14611476    if (FindPoint != PointsOnBoundary.end()) {
    14621477      Points[i] = FindPoint->second;
    … …  
    14711486      for (int j = i; j < 3; j++) {
    14721487        if (Points[j] != NULL) {
    1473           LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
     1488          LineMap::const_iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
    14741489          for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
    14751490            TriangleMap *triangles = &FindLine->second->triangles;
    14761491            Log() << Verbose(3) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl;
    1477             for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
     1492            for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
    14781493              if (FindTriangle->second->IsPresentTupel(Points)) {
    14791494                adjacentTriangleCount++;
    … …  
    17001715};
    17011716
     1717/** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates.
     1718 * This is supposed to prevent early closing of the tesselation.
     1719 * \param *BaseRay baseline, i.e. not \a *OptCandidate
     1720 * \param *ThirdNode third point in triangle, not in BoundaryLineSet::endpoints
     1721 * \param ShortestAngle path length on this circle band for the current \a *ThirdNode
     1722 * \param RADIUS radius of sphere
     1723 * \param *LC LinkedCell structure
     1724 * \return true - there is a better candidate (smaller angle than \a ShortestAngle), false - no better TesselPoint candidate found
     1725 */
     1726bool Tesselation::HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const ThirdNode, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const
     1727{
     1728  bool result = false;
     1729  Vector CircleCenter;
     1730  Vector CirclePlaneNormal;
     1731  Vector OldSphereCenter;
     1732  Vector SearchDirection;
     1733  Vector helper;
     1734  TesselPoint *OtherOptCandidate = NULL;
     1735  double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant.
     1736  double radius, CircleRadius;
     1737  BoundaryLineSet *Line = NULL;
     1738  BoundaryTriangleSet *T = NULL;
     1739
     1740  Log() << Verbose(1) << "Begin of HasOtherBaselineBetterCandidate" << endl;
     1741
     1742  // check both other lines
     1743  PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr);
     1744  if (FindPoint != PointsOnBoundary.end()) {
     1745    for (int i=0;i<2;i++) {
     1746      LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr);
     1747      if (FindLine != (FindPoint->second)->lines.end()) {
     1748        Line = FindLine->second;
     1749        Log() << Verbose(1) << "Found line " << *Line << "." << endl;
     1750        if (Line->triangles.size() == 1) {
     1751          T = Line->triangles.begin()->second;
     1752          // construct center of circle
     1753          CircleCenter.CopyVector(Line->endpoints[0]->node->node);
     1754          CircleCenter.AddVector(Line->endpoints[1]->node->node);
     1755          CircleCenter.Scale(0.5);
     1756
     1757          // construct normal vector of circle
     1758          CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node);
     1759          CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node);
     1760
     1761          // calculate squared radius of circle
     1762          radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
     1763          if (radius/4. < RADIUS*RADIUS) {
     1764            CircleRadius = RADIUS*RADIUS - radius/4.;
     1765            CirclePlaneNormal.Normalize();
     1766            //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
     1767
     1768            // construct old center
     1769            GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node);
     1770            helper.CopyVector(&T->NormalVector);  // normal vector ensures that this is correct center of the two possible ones
     1771            radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
     1772            helper.Scale(sqrt(RADIUS*RADIUS - radius));
     1773            OldSphereCenter.AddVector(&helper);
     1774            OldSphereCenter.SubtractVector(&CircleCenter);
     1775            //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
     1776
     1777            // construct SearchDirection
     1778            SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal);
     1779            helper.CopyVector(Line->endpoints[0]->node->node);
     1780            helper.SubtractVector(ThirdNode->node);
     1781            if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
     1782              SearchDirection.Scale(-1.);
     1783            SearchDirection.ProjectOntoPlane(&OldSphereCenter);
     1784            SearchDirection.Normalize();
     1785            Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
     1786            if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
     1787              // rotated the wrong way!
     1788              eLog() << Verbose(0) << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
     1789            }
     1790
     1791            // add third point
     1792            CandidateList *OptCandidates = new CandidateList();
     1793            FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, Line, ThirdNode, OptCandidates, &OtherShortestAngle, RADIUS, LC);
     1794            for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
     1795              if (((*it)->point == BaseRay->endpoints[0]->node) || ((*it)->point == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested
     1796                continue;
     1797              Log() << Verbose(1) << " Third point candidate is " << *(*it)->point
     1798              << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
     1799              Log() << Verbose(1) << " Baseline is " << *BaseRay << endl;
     1800
     1801              // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
     1802              TesselPoint *PointCandidates[3];
     1803              PointCandidates[0] = (*it)->point;
     1804              PointCandidates[1] = BaseRay->endpoints[0]->node;
     1805              PointCandidates[2] = BaseRay->endpoints[1]->node;
     1806              bool check=false;
     1807              int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates);
     1808              // If there is no triangle, add it regularly.
     1809              if (existentTrianglesCount == 0) {
     1810                SetTesselationPoint((*it)->point, 0);
     1811                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
     1812                SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
     1813
     1814                if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
     1815                  OtherOptCandidate = (*it)->point;
     1816                  check = true;
     1817                }
     1818              } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
     1819                SetTesselationPoint((*it)->point, 0);
     1820                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
     1821                SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
     1822
     1823                // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
     1824                // i.e. at least one of the three lines must be present with TriangleCount <= 1
     1825                if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) {
     1826                  OtherOptCandidate = (*it)->point;
     1827                  check = true;
     1828                }
     1829              }
     1830
     1831              if (check) {
     1832                if (ShortestAngle > OtherShortestAngle) {
     1833                  Log() << Verbose(1) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl;
     1834                  result = true;
     1835                  break;
     1836                }
     1837              }
     1838            }
     1839            delete(OptCandidates);
     1840            if (result)
     1841              break;
     1842          } else {
     1843            Log() << Verbose(1) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl;
     1844          }
     1845        } else {
     1846          eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl;
     1847        }
     1848      } else {
     1849        Log() << Verbose(2) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl;
     1850      }
     1851    }
     1852  } else {
     1853    eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl;
     1854  }
     1855
     1856  Log() << Verbose(1) << "End of HasOtherBaselineBetterCandidate" << endl;
     1857
     1858  return result;
     1859};
    17021860
    17031861/** This function finds a triangle to a line, adjacent to an existing one.
    … …  
    17821940  Log() << Verbose(1) << "Third Points are ";
    17831941  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
    1784     Log() << Verbose(0) << " " << *(*it)->point;
    1785   }
    1786   Log() << Verbose(0) << endl;
     1942    Log() << Verbose(1) << " " << *(*it)->point << endl;
     1943  }
    17871944
    17881945  BoundaryLineSet *BaseRay = &Line;
    … …  
    18001957
    18011958    BTS = NULL;
    1802     // If there is no triangle, add it regularly.
    1803     if (existentTrianglesCount == 0) {
    1804       AddTesselationPoint((*it)->point, 0);
    1805       AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
    1806       AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
    1807 
    1808       if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
    1809         AddTesselationLine(TPS[0], TPS[1], 0);
    1810         AddTesselationLine(TPS[0], TPS[2], 1);
    1811         AddTesselationLine(TPS[1], TPS[2], 2);
    1812 
    1813         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    1814         AddTesselationTriangle();
    1815         (*it)->OptCenter.Scale(-1.);
    1816         BTS->GetNormalVector((*it)->OptCenter);
    1817         (*it)->OptCenter.Scale(-1.);
    1818 
    1819         Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
    1820           << " for this triangle ... " << endl;
    1821       //Log() << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
    1822       } else {
    1823         Log() << Verbose(1) << "WARNING: This triangle consisting of ";
    1824         Log() << Verbose(0) << *(*it)->point << ", ";
    1825         Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
    1826         Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
    1827         Log() << Verbose(0) << "exists and is not added, as it does not seem helpful!" << endl;
    1828         result = false;
    1829       }
    1830     } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
     1959    // check for present edges and whether we reach better candidates from them
     1960    if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) {
     1961      result = false;
     1962      break;
     1963    } else {
     1964      // If there is no triangle, add it regularly.
     1965      if (existentTrianglesCount == 0) {
    18311966        AddTesselationPoint((*it)->point, 0);
    18321967        AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
    18331968        AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
    18341969
    1835         // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
    1836         // i.e. at least one of the three lines must be present with TriangleCount <= 1
    1837         if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) {
     1970        if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
    18381971          AddTesselationLine(TPS[0], TPS[1], 0);
    18391972          AddTesselationLine(TPS[0], TPS[2], 1);
    … …  
    18411974
    18421975          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    1843           AddTesselationTriangle();  // add to global map
    1844 
    1845           (*it)->OtherOptCenter.Scale(-1.);
    1846           BTS->GetNormalVector((*it)->OtherOptCenter);
    1847           (*it)->OtherOptCenter.Scale(-1.);
    1848 
    1849           Log() << Verbose(0) << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
    1850           << " for this triangle ... " << endl;
    1851           Log() << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl;
     1976          AddTesselationTriangle();
     1977          (*it)->OptCenter.Scale(-1.);
     1978          BTS->GetNormalVector((*it)->OptCenter);
     1979          (*it)->OptCenter.Scale(-1.);
     1980
     1981          Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
     1982            << " for this triangle ... " << endl;
     1983        //Log() << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
    18521984        } else {
    18531985          Log() << Verbose(1) << "WARNING: This triangle consisting of ";
    … …  
    18581990          result = false;
    18591991        }
    1860     } else {
    1861       Log() << Verbose(1) << "This triangle consisting of ";
    1862       Log() << Verbose(0) << *(*it)->point << ", ";
    1863       Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
    1864       Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
    1865       Log() << Verbose(0) << "is invalid!" << endl;
    1866       result = false;
     1992      } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
     1993          AddTesselationPoint((*it)->point, 0);
     1994          AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
     1995          AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
     1996
     1997          // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
     1998          // i.e. at least one of the three lines must be present with TriangleCount <= 1
     1999          if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) {
     2000            AddTesselationLine(TPS[0], TPS[1], 0);
     2001            AddTesselationLine(TPS[0], TPS[2], 1);
     2002            AddTesselationLine(TPS[1], TPS[2], 2);
     2003
     2004            BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     2005            AddTesselationTriangle();  // add to global map
     2006
     2007            (*it)->OtherOptCenter.Scale(-1.);
     2008            BTS->GetNormalVector((*it)->OtherOptCenter);
     2009            (*it)->OtherOptCenter.Scale(-1.);
     2010
     2011            Log() << Verbose(0) << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
     2012            << " for this triangle ... " << endl;
     2013            Log() << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl;
     2014          } else {
     2015            Log() << Verbose(1) << "WARNING: This triangle consisting of ";
     2016            Log() << Verbose(0) << *(*it)->point << ", ";
     2017            Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
     2018            Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
     2019            Log() << Verbose(0) << "exists and is not added, as it does not seem helpful!" << endl;
     2020            result = false;
     2021          }
     2022      } else {
     2023        Log() << Verbose(1) << "This triangle consisting of ";
     2024        Log() << Verbose(0) << *(*it)->point << ", ";
     2025        Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
     2026        Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
     2027        Log() << Verbose(0) << "is invalid!" << endl;
     2028        result = false;
     2029      }
    18672030    }
    18682031
    … …  
    21772340  }
    21782341  // then go through the current and all neighbouring cells and check the contained points for possible candidates
    2179   Log() << Verbose(3) << "LC Intervals from [";
    2180   for (int i=0;i<NDIM;i++) {
    2181   Log() << Verbose(0) << " " << N[i] << "<->" << LC->N[i];
    2182   }
    2183   Log() << Verbose(0) << "] :";
    21842342  for (int i=0;i<NDIM;i++) {
    21852343    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
    21862344    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
    2187     Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
    2188   }
    2189   Log() << Verbose(0) << endl;
    2190 
     2345  }
     2346  Log() << Verbose(3) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :"
     2347    << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl;
    21912348
    21922349  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    … …  
    22802437 * @param *LC LinkedCell structure with neighbouring points
    22812438 */
    2282 void Tesselation::FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint  *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC)
     2439void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, class BoundaryLineSet *BaseLine, const class TesselPoint  * const ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC) const
    22832440{
    22842441  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers

  • src/tesselation.hpp

    r9a0dc8 rf1ef60a  
    205205
    206206    void AddTesselationPoint(TesselPoint* Candidate, const int n);
     207    void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
    207208    void AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
    208209    void AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
    … …  
    217218    void FindStartingTriangle(const double RADIUS, const LinkedCell *LC);
    218219    void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC);
    219     void FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC);
     220    void FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, class BoundaryLineSet *BaseLine, const class TesselPoint * const ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC) const;
    220221    bool FindNextSuitableTriangle(BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
    221     int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]);
     222    int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]) const;
    222223    class BoundaryTriangleSet * GetPresentTriangle(TesselPoint *Candidates[3]);
    223224
    … …  
    281282
    282283  private:
    283     class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
     284    mutable class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
    284285
    285286    mutable PointMap::const_iterator InternalPointer;
     287
     288    bool HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const OptCandidate, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const;
    286289};
    287290

  • src/tesselationhelpers.cpp

    r9a0dc8 rf1ef60a  
    461461  for (int i=0;i<3;i++)
    462462    for (int j=i+1; j<3; j++) {
    463       if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
     463      if (nodes[i] == NULL) {
     464        Log() << Verbose(1) << "Node nr. " << i << " is not yet present." << endl;
     465        result = true;
     466      } else if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
    464467        LineMap::const_iterator FindLine;
    465468        pair<LineMap::const_iterator,LineMap::const_iterator> FindPair;
Note: See TracChangeset for help on using the changeset viewer.