Changeset d74077 for src/tesselation.cpp


Ignore:
Timestamp:
Jul 31, 2010, 3:23:10 PM (14 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:
8f4df1
Parents:
5fbaeb
Message:

Member variable Vector and element of class atom are now private.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselation.cpp

    r5fbaeb rd74077  
    1111#include <iomanip>
    1212
     13#include "BoundaryPointSet.hpp"
     14#include "BoundaryLineSet.hpp"
     15#include "BoundaryTriangleSet.hpp"
     16#include "BoundaryPolygonSet.hpp"
     17#include "TesselPoint.hpp"
     18#include "CandidateForTesselation.hpp"
     19#include "PointCloud.hpp"
    1320#include "helpers.hpp"
    1421#include "info.hpp"
    … …  
    2835class molecule;
    2936
    30 // ======================================== Points on Boundary =================================
    31 
    32 /** Constructor of BoundaryPointSet.
    33  */
    34 BoundaryPointSet::BoundaryPointSet() :
    35   LinesCount(0), value(0.), Nr(-1)
    36 {
    37   Info FunctionInfo(__func__);
    38   DoLog(1) && (Log() << Verbose(1) << "Adding noname." << endl);
    39 }
    40 ;
    41 
    42 /** Constructor of BoundaryPointSet with Tesselpoint.
    43  * \param *Walker TesselPoint this boundary point represents
    44  */
    45 BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
    46   LinesCount(0), node(Walker), value(0.), Nr(Walker->nr)
    47 {
    48   Info FunctionInfo(__func__);
    49   DoLog(1) && (Log() << Verbose(1) << "Adding Node " << *Walker << endl);
    50 }
    51 ;
    52 
    53 /** Destructor of BoundaryPointSet.
    54  * Sets node to NULL to avoid removing the original, represented TesselPoint.
    55  * \note When removing point from a class Tesselation, use RemoveTesselationPoint()
    56  */
    57 BoundaryPointSet::~BoundaryPointSet()
    58 {
    59   Info FunctionInfo(__func__);
    60   //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
    61   if (!lines.empty())
    62     DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
    63   node = NULL;
    64 }
    65 ;
    66 
    67 /** Add a line to the LineMap of this point.
    68  * \param *line line to add
    69  */
    70 void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
    71 {
    72   Info FunctionInfo(__func__);
    73   DoLog(1) && (Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." << endl);
    74   if (line->endpoints[0] == this) {
    75     lines.insert(LinePair(line->endpoints[1]->Nr, line));
    76   } else {
    77     lines.insert(LinePair(line->endpoints[0]->Nr, line));
    78   }
    79   LinesCount++;
    80 }
    81 ;
    82 
    83 /** output operator for BoundaryPointSet.
    84  * \param &ost output stream
    85  * \param &a boundary point
    86  */
    87 ostream & operator <<(ostream &ost, const BoundaryPointSet &a)
    88 {
    89   ost << "[" << a.Nr << "|" << a.node->getName() << " at " << *a.node->node << "]";
    90   return ost;
    91 }
    92 ;
    93 
    94 // ======================================== Lines on Boundary =================================
    95 
    96 /** Constructor of BoundaryLineSet.
    97  */
    98 BoundaryLineSet::BoundaryLineSet() :
    99   Nr(-1)
    100 {
    101   Info FunctionInfo(__func__);
    102   for (int i = 0; i < 2; i++)
    103     endpoints[i] = NULL;
    104 }
    105 ;
    106 
    107 /** Constructor of BoundaryLineSet with two endpoints.
    108  * Adds line automatically to each endpoints' LineMap
    109  * \param *Point[2] array of two boundary points
    110  * \param number number of the list
    111  */
    112 BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
    113 {
    114   Info FunctionInfo(__func__);
    115   // set number
    116   Nr = number;
    117   // set endpoints in ascending order
    118   SetEndpointsOrdered(endpoints, Point[0], Point[1]);
    119   // add this line to the hash maps of both endpoints
    120   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
    121   Point[1]->AddLine(this); //
    122   // set skipped to false
    123   skipped = false;
    124   // clear triangles list
    125   DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl);
    126 }
    127 ;
    128 
    129 /** Constructor of BoundaryLineSet with two endpoints.
    130  * Adds line automatically to each endpoints' LineMap
    131  * \param *Point1 first boundary point
    132  * \param *Point2 second boundary point
    133  * \param number number of the list
    134  */
    135 BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number)
    136 {
    137   Info FunctionInfo(__func__);
    138   // set number
    139   Nr = number;
    140   // set endpoints in ascending order
    141   SetEndpointsOrdered(endpoints, Point1, Point2);
    142   // add this line to the hash maps of both endpoints
    143   Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
    144   Point2->AddLine(this); //
    145   // set skipped to false
    146   skipped = false;
    147   // clear triangles list
    148   DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl);
    149 }
    150 ;
    151 
    152 /** Destructor for BoundaryLineSet.
    153  * Removes itself from each endpoints' LineMap, calling RemoveTrianglePoint() when point not connected anymore.
    154  * \note When removing lines from a class Tesselation, use RemoveTesselationLine()
    155  */
    156 BoundaryLineSet::~BoundaryLineSet()
    157 {
    158   Info FunctionInfo(__func__);
    159   int Numbers[2];
    160 
    161   // get other endpoint number of finding copies of same line
    162   if (endpoints[1] != NULL)
    163     Numbers[0] = endpoints[1]->Nr;
    164   else
    165     Numbers[0] = -1;
    166   if (endpoints[0] != NULL)
    167     Numbers[1] = endpoints[0]->Nr;
    168   else
    169     Numbers[1] = -1;
    170 
    171   for (int i = 0; i < 2; i++) {
    172     if (endpoints[i] != NULL) {
    173       if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
    174         pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
    175         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
    176           if ((*Runner).second == this) {
    177             //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
    178             endpoints[i]->lines.erase(Runner);
    179             break;
    180           }
    181       } else { // there's just a single line left
    182         if (endpoints[i]->lines.erase(Nr)) {
    183           //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
    184         }
    185       }
    186       if (endpoints[i]->lines.empty()) {
    187         //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
    188         if (endpoints[i] != NULL) {
    189           delete (endpoints[i]);
    190           endpoints[i] = NULL;
    191         }
    192       }
    193     }
    194   }
    195   if (!triangles.empty())
    196     DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
    197 }
    198 ;
    199 
    200 /** Add triangle to TriangleMap of this boundary line.
    201  * \param *triangle to add
    202  */
    203 void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
    204 {
    205   Info FunctionInfo(__func__);
    206   DoLog(0) && (Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl);
    207   triangles.insert(TrianglePair(triangle->Nr, triangle));
    208 }
    209 ;
    210 
    211 /** Checks whether we have a common endpoint with given \a *line.
    212  * \param *line other line to test
    213  * \return true - common endpoint present, false - not connected
    214  */
    215 bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
    216 {
    217   Info FunctionInfo(__func__);
    218   if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
    219     return true;
    220   else
    221     return false;
    222 }
    223 ;
    224 
    225 /** Checks whether the adjacent triangles of a baseline are convex or not.
    226  * We sum the two angles of each height vector with respect to the center of the baseline.
    227  * If greater/equal M_PI than we are convex.
    228  * \param *out output stream for debugging
    229  * \return true - triangles are convex, false - concave or less than two triangles connected
    230  */
    231 bool BoundaryLineSet::CheckConvexityCriterion() const
    232 {
    233   Info FunctionInfo(__func__);
    234   double angle = CalculateConvexity();
    235   if (angle > -MYEPSILON) {
    236     DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
    237     return true;
    238   } else {
    239     DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
    240     return false;
    241   }
    242 }
    243 
    244 
    245 /** Calculates the angle between two triangles with respect to their normal vector.
    246  * We sum the two angles of each height vector with respect to the center of the baseline.
    247  * \return angle > 0 then convex, if < 0 then concave
    248  */
    249 double BoundaryLineSet::CalculateConvexity() const
    250 {
    251   Info FunctionInfo(__func__);
    252   Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
    253   // get the two triangles
    254   if (triangles.size() != 2) {
    255     DoeLog(0) && (eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
    256     return true;
    257   }
    258   // check normal vectors
    259   // have a normal vector on the base line pointing outwards
    260   //Log() << Verbose(0) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
    261   BaseLineCenter = (1./2.)*((*endpoints[0]->node->node) + (*endpoints[1]->node->node));
    262   BaseLine = (*endpoints[0]->node->node) - (*endpoints[1]->node->node);
    263 
    264   //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
    265 
    266   BaseLineNormal.Zero();
    267   NormalCheck.Zero();
    268   double sign = -1.;
    269   int i = 0;
    270   class BoundaryPointSet *node = NULL;
    271   for (TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
    272     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
    273     NormalCheck += runner->second->NormalVector;
    274     NormalCheck *= sign;
    275     sign = -sign;
    276     if (runner->second->NormalVector.NormSquared() > MYEPSILON)
    277       BaseLineNormal = runner->second->NormalVector;   // yes, copy second on top of first
    278     else {
    279       DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
    280     }
    281     node = runner->second->GetThirdEndpoint(this);
    282     if (node != NULL) {
    283       //Log() << Verbose(0) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
    284       helper[i] = (*node->node->node) - BaseLineCenter;
    285       helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
    286       //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
    287       i++;
    288     } else {
    289       DoeLog(1) && (eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
    290       return true;
    291     }
    292   }
    293   //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
    294   if (NormalCheck.NormSquared() < MYEPSILON) {
    295     DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl);
    296     return true;
    297   }
    298   BaseLineNormal.Scale(-1.);
    299   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
    300   return (angle - M_PI);
    301 }
    302 
    303 /** Checks whether point is any of the two endpoints this line contains.
    304  * \param *point point to test
    305  * \return true - point is of the line, false - is not
    306  */
    307 bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
    308 {
    309   Info FunctionInfo(__func__);
    310   for (int i = 0; i < 2; i++)
    311     if (point == endpoints[i])
    312       return true;
    313   return false;
    314 }
    315 ;
    316 
    317 /** Returns other endpoint of the line.
    318  * \param *point other endpoint
    319  * \return NULL - if endpoint not contained in BoundaryLineSet::lines, or pointer to BoundaryPointSet otherwise
    320  */
    321 class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
    322 {
    323   Info FunctionInfo(__func__);
    324   if (endpoints[0] == point)
    325     return endpoints[1];
    326   else if (endpoints[1] == point)
    327     return endpoints[0];
    328   else
    329     return NULL;
    330 }
    331 ;
    332 
    333 /** Returns other triangle of the line.
    334  * \param *point other endpoint
    335  * \return NULL - if triangle not contained in BoundaryLineSet::triangles, or pointer to BoundaryTriangleSet otherwise
    336  */
    337 class BoundaryTriangleSet *BoundaryLineSet::GetOtherTriangle(const BoundaryTriangleSet * const triangle) const
    338 {
    339   Info FunctionInfo(__func__);
    340   if (triangles.size() == 2) {
    341     for (TriangleMap::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); ++TriangleRunner)
    342       if (TriangleRunner->second != triangle)
    343         return TriangleRunner->second;
    344   }
    345   return NULL;
    346 }
    347 ;
    348 
    349 /** output operator for BoundaryLineSet.
    350  * \param &ost output stream
    351  * \param &a boundary line
    352  */
    353 ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
    354 {
    355   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->getName() << " at " << *a.endpoints[1]->node->node << "]";
    356   return ost;
    357 }
    358 ;
    359 
    360 // ======================================== Triangles on Boundary =================================
    361 
    362 /** Constructor for BoundaryTriangleSet.
    363  */
    364 BoundaryTriangleSet::BoundaryTriangleSet() :
    365   Nr(-1)
    366 {
    367   Info FunctionInfo(__func__);
    368   for (int i = 0; i < 3; i++) {
    369     endpoints[i] = NULL;
    370     lines[i] = NULL;
    371   }
    372 }
    373 ;
    374 
    375 /** Constructor for BoundaryTriangleSet with three lines.
    376  * \param *line[3] lines that make up the triangle
    377  * \param number number of triangle
    378  */
    379 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
    380   Nr(number)
    381 {
    382   Info FunctionInfo(__func__);
    383   // set number
    384   // set lines
    385   for (int i = 0; i < 3; i++) {
    386     lines[i] = line[i];
    387     lines[i]->AddTriangle(this);
    388   }
    389   // get ascending order of endpoints
    390   PointMap OrderMap;
    391   for (int i = 0; i < 3; i++) {
    392     // for all three lines
    393     for (int j = 0; j < 2; j++) { // for both endpoints
    394       OrderMap.insert(pair<int, class BoundaryPointSet *> (line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
    395       // and we don't care whether insertion fails
    396     }
    397   }
    398   // set endpoints
    399   int Counter = 0;
    400   DoLog(0) && (Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl);
    401   for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
    402     endpoints[Counter] = runner->second;
    403     DoLog(0) && (Log() << Verbose(0) << " " << *endpoints[Counter] << endl);
    404     Counter++;
    405   }
    406   ASSERT(Counter >= 3,"We have a triangle with only two distinct endpoints!");
    407 };
    408 
    409 
    410 /** Destructor of BoundaryTriangleSet.
    411  * Removes itself from each of its lines' LineMap and removes them if necessary.
    412  * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
    413  */
    414 BoundaryTriangleSet::~BoundaryTriangleSet()
    415 {
    416   Info FunctionInfo(__func__);
    417   for (int i = 0; i < 3; i++) {
    418     if (lines[i] != NULL) {
    419       if (lines[i]->triangles.erase(Nr)) {
    420         //Log() << Verbose(0) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
    421       }
    422       if (lines[i]->triangles.empty()) {
    423         //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
    424         delete (lines[i]);
    425         lines[i] = NULL;
    426       }
    427     }
    428   }
    429   //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl;
    430 }
    431 ;
    432 
    433 /** Calculates the normal vector for this triangle.
    434  * Is made unique by comparison with \a OtherVector to point in the other direction.
    435  * \param &OtherVector direction vector to make normal vector unique.
    436  */
    437 void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
    438 {
    439   Info FunctionInfo(__func__);
    440   // get normal vector
    441   NormalVector = Plane(*(endpoints[0]->node->node),
    442                        *(endpoints[1]->node->node),
    443                        *(endpoints[2]->node->node)).getNormal();
    444 
    445   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
    446   if (NormalVector.ScalarProduct(OtherVector) > 0.)
    447     NormalVector.Scale(-1.);
    448   DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl);
    449 }
    450 ;
    451 
    452 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
    453  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
    454  * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
    455  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
    456  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
    457  * the first two basepoints) or not.
    458  * \param *out output stream for debugging
    459  * \param *MolCenter offset vector of line
    460  * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
    461  * \param *Intersection intersection on plane on return
    462  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
    463  */
    464 
    465 bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
    466 {
    467   Info FunctionInfo(__func__);
    468   Vector CrossPoint;
    469   Vector helper;
    470 
    471   try {
    472     Line centerLine = makeLineThrough(*MolCenter, *x);
    473     *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);
    474 
    475     DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
    476     DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
    477     DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
    478 
    479     if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
    480       DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
    481       return true;
    482     }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
    483       DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
    484       return true;
    485     }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
    486       DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
    487       return true;
    488     }
    489     // Calculate cross point between one baseline and the line from the third endpoint to intersection
    490     int i = 0;
    491     do {
    492       Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
    493       Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
    494       CrossPoint = line1.getIntersection(line2);
    495       helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
    496       CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
    497       const double s = CrossPoint.ScalarProduct(helper)/helper.NormSquared();
    498       DoLog(1) && (Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl);
    499       if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
    500         DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
    501         return false;
    502       }
    503       i++;
    504     } while (i < 3);
    505     DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
    506     return true;
    507   }
    508   catch (MathException &excp) {
    509     Log() << Verbose(1) << excp;
    510     DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
    511     return false;
    512   }
    513 }
    514 ;
    515 
    516 /** Finds the point on the triangle to the point \a *x.
    517  * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point.
    518  * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the
    519  * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down.
    520  * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
    521  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
    522  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
    523  * the first two basepoints) or not.
    524  * \param *x point
    525  * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
    526  * \return Distance squared between \a *x and closest point inside triangle
    527  */
    528 double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
    529 {
    530   Info FunctionInfo(__func__);
    531   Vector Direction;
    532 
    533   // 1. get intersection with plane
    534   DoLog(1) && (Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl);
    535   GetCenter(&Direction);
    536   try {
    537     Line l = makeLineThrough(*x, Direction);
    538     *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
    539   }
    540   catch (MathException &excp) {
    541     (*ClosestPoint) = (*x);
    542   }
    543 
    544   // 2. Calculate in plane part of line (x, intersection)
    545   Vector InPlane = (*x) - (*ClosestPoint); // points from plane intersection to straight-down point
    546   InPlane.ProjectOntoPlane(NormalVector);
    547   InPlane += *ClosestPoint;
    548 
    549   DoLog(2) && (Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl);
    550   DoLog(2) && (Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl);
    551   DoLog(2) && (Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl);
    552 
    553   // Calculate cross point between one baseline and the desired point such that distance is shortest
    554   double ShortestDistance = -1.;
    555   bool InsideFlag = false;
    556   Vector CrossDirection[3];
    557   Vector CrossPoint[3];
    558   Vector helper;
    559   for (int i = 0; i < 3; i++) {
    560     // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
    561     Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
    562     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
    563     Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
    564     CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
    565     CrossDirection[i] = CrossPoint[i] - InPlane;
    566     CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
    567     const double s = CrossPoint[i].ScalarProduct(Direction)/Direction.NormSquared();
    568     DoLog(2) && (Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl);
    569     if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
    570           CrossPoint[i] += (*endpoints[i%3]->node->node);  // make cross point absolute again
    571       DoLog(2) && (Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i % 3]->node->node << " and " << *endpoints[(i + 1) % 3]->node->node << "." << endl);
    572       const double distance = CrossPoint[i].DistanceSquared(*x);
    573       if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
    574         ShortestDistance = distance;
    575         (*ClosestPoint) = CrossPoint[i];
    576       }
    577     } else
    578       CrossPoint[i].Zero();
    579   }
    580   InsideFlag = true;
    581   for (int i = 0; i < 3; i++) {
    582     const double sign = CrossDirection[i].ScalarProduct(CrossDirection[(i + 1) % 3]);
    583     const double othersign = CrossDirection[i].ScalarProduct(CrossDirection[(i + 2) % 3]);
    584 
    585     if ((sign > -MYEPSILON) && (othersign > -MYEPSILON)) // have different sign
    586       InsideFlag = false;
    587   }
    588   if (InsideFlag) {
    589     (*ClosestPoint) = InPlane;
    590     ShortestDistance = InPlane.DistanceSquared(*x);
    591   } else { // also check endnodes
    592     for (int i = 0; i < 3; i++) {
    593       const double distance = x->DistanceSquared(*endpoints[i]->node->node);
    594       if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
    595         ShortestDistance = distance;
    596         (*ClosestPoint) = (*endpoints[i]->node->node);
    597       }
    598     }
    599   }
    600   DoLog(1) && (Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl);
    601   return ShortestDistance;
    602 }
    603 ;
    604 
    605 /** Checks whether lines is any of the three boundary lines this triangle contains.
    606  * \param *line line to test
    607  * \return true - line is of the triangle, false - is not
    608  */
    609 bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
    610 {
    611   Info FunctionInfo(__func__);
    612   for (int i = 0; i < 3; i++)
    613     if (line == lines[i])
    614       return true;
    615   return false;
    616 }
    617 ;
    618 
    619 /** Checks whether point is any of the three endpoints this triangle contains.
    620  * \param *point point to test
    621  * \return true - point is of the triangle, false - is not
    622  */
    623 bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
    624 {
    625   Info FunctionInfo(__func__);
    626   for (int i = 0; i < 3; i++)
    627     if (point == endpoints[i])
    628       return true;
    629   return false;
    630 }
    631 ;
    632 
    633 /** Checks whether point is any of the three endpoints this triangle contains.
    634  * \param *point TesselPoint to test
    635  * \return true - point is of the triangle, false - is not
    636  */
    637 bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
    638 {
    639   Info FunctionInfo(__func__);
    640   for (int i = 0; i < 3; i++)
    641     if (point == endpoints[i]->node)
    642       return true;
    643   return false;
    644 }
    645 ;
    646 
    647 /** Checks whether three given \a *Points coincide with triangle's endpoints.
    648  * \param *Points[3] pointer to BoundaryPointSet
    649  * \return true - is the very triangle, false - is not
    650  */
    651 bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
    652 {
    653   Info FunctionInfo(__func__);
    654   DoLog(1) && (Log() << Verbose(1) << "INFO: Checking " << Points[0] << "," << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl);
    655   return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])
    656 
    657   ));
    658 }
    659 ;
    660 
    661 /** Checks whether three given \a *Points coincide with triangle's endpoints.
    662  * \param *Points[3] pointer to BoundaryPointSet
    663  * \return true - is the very triangle, false - is not
    664  */
    665 bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
    666 {
    667   Info FunctionInfo(__func__);
    668   return (((endpoints[0] == T->endpoints[0]) || (endpoints[0] == T->endpoints[1]) || (endpoints[0] == T->endpoints[2])) && ((endpoints[1] == T->endpoints[0]) || (endpoints[1] == T->endpoints[1]) || (endpoints[1] == T->endpoints[2])) && ((endpoints[2] == T->endpoints[0]) || (endpoints[2] == T->endpoints[1]) || (endpoints[2] == T->endpoints[2])
    669 
    670   ));
    671 }
    672 ;
    673 
    674 /** Returns the endpoint which is not contained in the given \a *line.
    675  * \param *line baseline defining two endpoints
    676  * \return pointer third endpoint or NULL if line does not belong to triangle.
    677  */
    678 class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
    679 {
    680   Info FunctionInfo(__func__);
    681   // sanity check
    682   if (!ContainsBoundaryLine(line))
    683     return NULL;
    684   for (int i = 0; i < 3; i++)
    685     if (!line->ContainsBoundaryPoint(endpoints[i]))
    686       return endpoints[i];
    687   // actually, that' impossible :)
    688   return NULL;
    689 }
    690 ;
    691 
    692 /** Returns the baseline which does not contain the given boundary point \a *point.
    693  * \param *point endpoint which is neither endpoint of the desired line
    694  * \return pointer to desired third baseline
    695  */
    696 class BoundaryLineSet *BoundaryTriangleSet::GetThirdLine(const BoundaryPointSet * const point) const
    697 {
    698   Info FunctionInfo(__func__);
    699   // sanity check
    700   if (!ContainsBoundaryPoint(point))
    701     return NULL;
    702   for (int i = 0; i < 3; i++)
    703     if (!lines[i]->ContainsBoundaryPoint(point))
    704       return lines[i];
    705   // actually, that' impossible :)
    706   return NULL;
    707 }
    708 ;
    709 
    710 /** Calculates the center point of the triangle.
    711  * Is third of the sum of all endpoints.
    712  * \param *center central point on return.
    713  */
    714 void BoundaryTriangleSet::GetCenter(Vector * const center) const
    715 {
    716   Info FunctionInfo(__func__);
    717   center->Zero();
    718   for (int i = 0; i < 3; i++)
    719     (*center) += (*endpoints[i]->node->node);
    720   center->Scale(1. / 3.);
    721   DoLog(1) && (Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl);
    722 }
    723 
    724 /**
    725  * gets the Plane defined by the three triangle Basepoints
    726  */
    727 Plane BoundaryTriangleSet::getPlane() const{
    728   ASSERT(endpoints[0] && endpoints[1] && endpoints[2], "Triangle not fully defined");
    729 
    730   return Plane(*endpoints[0]->node->node,
    731                *endpoints[1]->node->node,
    732                *endpoints[2]->node->node);
    733 }
    734 
    735 Vector BoundaryTriangleSet::getEndpoint(int i) const{
    736   ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
    737 
    738   return *endpoints[i]->node->node;
    739 }
    740 
    741 string BoundaryTriangleSet::getEndpointName(int i) const{
    742   ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
    743 
    744   return endpoints[i]->node->getName();
    745 }
    746 
    747 /** output operator for BoundaryTriangleSet.
    748  * \param &ost output stream
    749  * \param &a boundary triangle
    750  */
    751 ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a)
    752 {
    753   ost << "[" << a.Nr << "|" << a.getEndpointName(0) << "," << a.getEndpointName(1) << "," << a.getEndpointName(2) << "]";
    754   //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
    755   //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
    756   return ost;
    757 }
    758 ;
    759 
    760 // ======================================== Polygons on Boundary =================================
    761 
    762 /** Constructor for BoundaryPolygonSet.
    763  */
    764 BoundaryPolygonSet::BoundaryPolygonSet() :
    765   Nr(-1)
    766 {
    767   Info FunctionInfo(__func__);
    768 }
    769 ;
    770 
    771 /** Destructor of BoundaryPolygonSet.
    772  * Just clears endpoints.
    773  * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
    774  */
    775 BoundaryPolygonSet::~BoundaryPolygonSet()
    776 {
    777   Info FunctionInfo(__func__);
    778   endpoints.clear();
    779   DoLog(1) && (Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl);
    780 }
    781 ;
    782 
    783 /** Calculates the normal vector for this triangle.
    784  * Is made unique by comparison with \a OtherVector to point in the other direction.
    785  * \param &OtherVector direction vector to make normal vector unique.
    786  * \return allocated vector in normal direction
    787  */
    788 Vector * BoundaryPolygonSet::GetNormalVector(const Vector &OtherVector) const
    789 {
    790   Info FunctionInfo(__func__);
    791   // get normal vector
    792   Vector TemporaryNormal;
    793   Vector *TotalNormal = new Vector;
    794   PointSet::const_iterator Runner[3];
    795   for (int i = 0; i < 3; i++) {
    796     Runner[i] = endpoints.begin();
    797     for (int j = 0; j < i; j++) { // go as much further
    798       Runner[i]++;
    799       if (Runner[i] == endpoints.end()) {
    800         DoeLog(0) && (eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
    801         performCriticalExit();
    802       }
    803     }
    804   }
    805   TotalNormal->Zero();
    806   int counter = 0;
    807   for (; Runner[2] != endpoints.end();) {
    808     TemporaryNormal = Plane(*((*Runner[0])->node->node),
    809                             *((*Runner[1])->node->node),
    810                             *((*Runner[2])->node->node)).getNormal();
    811     for (int i = 0; i < 3; i++) // increase each of them
    812       Runner[i]++;
    813     (*TotalNormal) += TemporaryNormal;
    814   }
    815   TotalNormal->Scale(1. / (double) counter);
    816 
    817   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
    818   if (TotalNormal->ScalarProduct(OtherVector) > 0.)
    819     TotalNormal->Scale(-1.);
    820   DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl);
    821 
    822   return TotalNormal;
    823 }
    824 ;
    825 
    826 /** Calculates the center point of the triangle.
    827  * Is third of the sum of all endpoints.
    828  * \param *center central point on return.
    829  */
    830 void BoundaryPolygonSet::GetCenter(Vector * const center) const
    831 {
    832   Info FunctionInfo(__func__);
    833   center->Zero();
    834   int counter = 0;
    835   for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
    836     (*center) += (*(*Runner)->node->node);
    837     counter++;
    838   }
    839   center->Scale(1. / (double) counter);
    840   DoLog(1) && (Log() << Verbose(1) << "Center is at " << *center << "." << endl);
    841 }
    842 
    843 /** Checks whether the polygons contains all three endpoints of the triangle.
    844  * \param *triangle triangle to test
    845  * \return true - triangle is contained polygon, false - is not
    846  */
    847 bool BoundaryPolygonSet::ContainsBoundaryTriangle(const BoundaryTriangleSet * const triangle) const
    848 {
    849   Info FunctionInfo(__func__);
    850   return ContainsPresentTupel(triangle->endpoints, 3);
    851 }
    852 ;
    853 
    854 /** Checks whether the polygons contains both endpoints of the line.
    855  * \param *line line to test
    856  * \return true - line is of the triangle, false - is not
    857  */
    858 bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
    859 {
    860   Info FunctionInfo(__func__);
    861   return ContainsPresentTupel(line->endpoints, 2);
    862 }
    863 ;
    864 
    865 /** Checks whether point is any of the three endpoints this triangle contains.
    866  * \param *point point to test
    867  * \return true - point is of the triangle, false - is not
    868  */
    869 bool BoundaryPolygonSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
    870 {
    871   Info FunctionInfo(__func__);
    872   for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
    873     DoLog(0) && (Log() << Verbose(0) << "Checking against " << **Runner << endl);
    874     if (point == (*Runner)) {
    875       DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
    876       return true;
    877     }
    878   }
    879   DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
    880   return false;
    881 }
    882 ;
    883 
    884 /** Checks whether point is any of the three endpoints this triangle contains.
    885  * \param *point TesselPoint to test
    886  * \return true - point is of the triangle, false - is not
    887  */
    888 bool BoundaryPolygonSet::ContainsBoundaryPoint(const TesselPoint * const point) const
    889 {
    890   Info FunctionInfo(__func__);
    891   for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
    892     if (point == (*Runner)->node) {
    893       DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
    894       return true;
    895     }
    896   DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
    897   return false;
    898 }
    899 ;
    900 
    901 /** Checks whether given array of \a *Points coincide with polygons's endpoints.
    902  * \param **Points pointer to an array of BoundaryPointSet
    903  * \param dim dimension of array
    904  * \return true - set of points is contained in polygon, false - is not
    905  */
    906 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const
    907 {
    908   Info FunctionInfo(__func__);
    909   int counter = 0;
    910   DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
    911   for (int i = 0; i < dim; i++) {
    912     DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl);
    913     if (ContainsBoundaryPoint(Points[i])) {
    914       counter++;
    915     }
    916   }
    917 
    918   if (counter == dim)
    919     return true;
    920   else
    921     return false;
    922 }
    923 ;
    924 
    925 /** Checks whether given PointList coincide with polygons's endpoints.
    926  * \param &endpoints PointList
    927  * \return true - set of points is contained in polygon, false - is not
    928  */
    929 bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const
    930 {
    931   Info FunctionInfo(__func__);
    932   size_t counter = 0;
    933   DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
    934   for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
    935     DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << **Runner << endl);
    936     if (ContainsBoundaryPoint(*Runner))
    937       counter++;
    938   }
    939 
    940   if (counter == endpoints.size())
    941     return true;
    942   else
    943     return false;
    944 }
    945 ;
    946 
    947 /** Checks whether given set of \a *Points coincide with polygons's endpoints.
    948  * \param *P pointer to BoundaryPolygonSet
    949  * \return true - is the very triangle, false - is not
    950  */
    951 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
    952 {
    953   return ContainsPresentTupel((const PointSet) P->endpoints);
    954 }
    955 ;
    956 
    957 /** Gathers all the endpoints' triangles in a unique set.
    958  * \return set of all triangles
    959  */
    960 TriangleSet * BoundaryPolygonSet::GetAllContainedTrianglesFromEndpoints() const
    961 {
    962   Info FunctionInfo(__func__);
    963   pair<TriangleSet::iterator, bool> Tester;
    964   TriangleSet *triangles = new TriangleSet;
    965 
    966   for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
    967     for (LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)
    968       for (TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {
    969         //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl;
    970         if (ContainsBoundaryTriangle(Sprinter->second)) {
    971           Tester = triangles->insert(Sprinter->second);
    972           if (Tester.second)
    973             DoLog(0) && (Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl);
    974         }
    975       }
    976 
    977   DoLog(1) && (Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl);
    978   return triangles;
    979 }
    980 ;
    981 
    982 /** Fills the endpoints of this polygon from the triangles attached to \a *line.
    983  * \param *line lines with triangles attached
    984  * \return true - polygon contains endpoints, false - line was NULL
    985  */
    986 bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)
    987 {
    988   Info FunctionInfo(__func__);
    989   pair<PointSet::iterator, bool> Tester;
    990   if (line == NULL)
    991     return false;
    992   DoLog(1) && (Log() << Verbose(1) << "Filling polygon from line " << *line << endl);
    993   for (TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {
    994     for (int i = 0; i < 3; i++) {
    995       Tester = endpoints.insert((Runner->second)->endpoints[i]);
    996       if (Tester.second)
    997         DoLog(1) && (Log() << Verbose(1) << "  Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl);
    998     }
    999   }
    1000 
    1001   return true;
    1002 }
    1003 ;
    1004 
    1005 /** output operator for BoundaryPolygonSet.
    1006  * \param &ost output stream
    1007  * \param &a boundary polygon
    1008  */
    1009 ostream &operator <<(ostream &ost, const BoundaryPolygonSet &a)
    1010 {
    1011   ost << "[" << a.Nr << "|";
    1012   for (PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {
    1013     ost << (*Runner)->node->getName();
    1014     Runner++;
    1015     if (Runner != a.endpoints.end())
    1016       ost << ",";
    1017   }
    1018   ost << "]";
    1019   return ost;
    1020 }
    1021 ;
    1022 
    1023 // =========================================================== class TESSELPOINT ===========================================
    1024 
    1025 /** Constructor of class TesselPoint.
    1026  */
    1027 TesselPoint::TesselPoint()
    1028 {
    1029   //Info FunctionInfo(__func__);
    1030   node = NULL;
    1031   nr = -1;
    1032 }
    1033 ;
    1034 
    1035 /** Destructor for class TesselPoint.
    1036  */
    1037 TesselPoint::~TesselPoint()
    1038 {
    1039   //Info FunctionInfo(__func__);
    1040 }
    1041 ;
    1042 
    1043 /** Prints LCNode to screen.
    1044  */
    1045 ostream & operator <<(ostream &ost, const TesselPoint &a)
    1046 {
    1047   ost << "[" << a.getName() << "|" << *a.node << "]";
    1048   return ost;
    1049 }
    1050 ;
    1051 
    1052 /** Prints LCNode to screen.
    1053  */
    1054 ostream & TesselPoint::operator <<(ostream &ost)
    1055 {
    1056   Info FunctionInfo(__func__);
    1057   ost << "[" << (nr) << "|" << this << "]";
    1058   return ost;
    1059 }
    1060 ;
    1061 
    1062 // =========================================================== class POINTCLOUD ============================================
    1063 
    1064 /** Constructor of class PointCloud.
    1065  */
    1066 PointCloud::PointCloud()
    1067 {
    1068   //Info FunctionInfo(__func__);
    1069 }
    1070 ;
    1071 
    1072 /** Destructor for class PointCloud.
    1073  */
    1074 PointCloud::~PointCloud()
    1075 {
    1076   //Info FunctionInfo(__func__);
    1077 }
    1078 ;
    1079 
    1080 // ============================ CandidateForTesselation =============================
    1081 
    1082 /** Constructor of class CandidateForTesselation.
    1083  */
    1084 CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) :
    1085   BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
    1086 {
    1087   Info FunctionInfo(__func__);
    1088 }
    1089 ;
    1090 
    1091 /** Constructor of class CandidateForTesselation.
    1092  */
    1093 CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
    1094   BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
    1095 {
    1096         Info FunctionInfo(__func__);
    1097   OptCenter = OptCandidateCenter;
    1098   OtherOptCenter = OtherOptCandidateCenter;
    1099 };
    1100 
    1101 
    1102 /** Destructor for class CandidateForTesselation.
    1103  */
    1104 CandidateForTesselation::~CandidateForTesselation()
    1105 {
    1106 }
    1107 ;
    1108 
    1109 /** Checks validity of a given sphere of a candidate line.
    1110  * Sphere must touch all candidates and the baseline endpoints and there must be no other atoms inside.
    1111  * \param RADIUS radius of sphere
    1112  * \param *LC LinkedCell structure with other atoms
    1113  * \return true - sphere is valid, false - sphere contains other points
    1114  */
    1115 bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell *LC) const
    1116 {
    1117   Info FunctionInfo(__func__);
    1118 
    1119   const double radiusSquared = RADIUS * RADIUS;
    1120   list<const Vector *> VectorList;
    1121   VectorList.push_back(&OptCenter);
    1122   //VectorList.push_back(&OtherOptCenter);  // don't check the other (wrong) center
    1123 
    1124   if (!pointlist.empty())
    1125     DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
    1126   else
    1127     DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
    1128   // check baseline for OptCenter and OtherOptCenter being on sphere's surface
    1129   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
    1130     for (int i = 0; i < 2; i++) {
    1131       const double distance = fabs((*VRunner)->DistanceSquared(*BaseLine->endpoints[i]->node->node) - radiusSquared);
    1132       if (distance > HULLEPSILON) {
    1133         DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
    1134         return false;
    1135       }
    1136     }
    1137   }
    1138 
    1139   // check Candidates for OptCenter and OtherOptCenter being on sphere's surface
    1140   for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
    1141     const TesselPoint *Walker = *Runner;
    1142     for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
    1143       const double distance = fabs((*VRunner)->DistanceSquared(*Walker->node) - radiusSquared);
    1144       if (distance > HULLEPSILON) {
    1145         DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
    1146         return false;
    1147       } else {
    1148         DoLog(1) && (Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl);
    1149       }
    1150     }
    1151   }
    1152 
    1153   DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
    1154   bool flag = true;
    1155   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
    1156     // get all points inside the sphere
    1157     TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));
    1158 
    1159     DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << (*VRunner) << ":" << endl);
    1160     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    1161       DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(*(*VRunner)) << "." << endl);
    1162 
    1163     // remove baseline's endpoints and candidates
    1164     for (int i = 0; i < 2; i++) {
    1165       DoLog(1) && (Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl);
    1166       ListofPoints->remove(BaseLine->endpoints[i]->node);
    1167     }
    1168     for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
    1169       DoLog(1) && (Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl);
    1170       ListofPoints->remove(*Runner);
    1171     }
    1172     if (!ListofPoints->empty()) {
    1173       DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
    1174       flag = false;
    1175       DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
    1176       for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    1177         DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << " at distance " << setprecision(13) << (*Runner)->node->distance(*(*VRunner)) << setprecision(6) << "." << endl);
    1178 
    1179       // check with animate_sphere.tcl VMD script
    1180       if (ThirdPoint != NULL) {
    1181         DoeLog(1) && (eLog() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
    1182       } else {
    1183         DoeLog(1) && (eLog() << Verbose(1) << "Check by: ... missing third point ..." << endl);
    1184         DoeLog(1) && (eLog() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
    1185       }
    1186     }
    1187     delete (ListofPoints);
    1188 
    1189   }
    1190   return flag;
    1191 }
    1192 ;
    1193 
    1194 /** output operator for CandidateForTesselation.
    1195  * \param &ost output stream
    1196  * \param &a boundary line
    1197  */
    1198 ostream & operator <<(ostream &ost, const CandidateForTesselation &a)
    1199 {
    1200   ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->getName() << "," << a.BaseLine->endpoints[1]->node->getName() << "] with ";
    1201   if (a.pointlist.empty())
    1202     ost << "no candidate.";
    1203   else {
    1204     ost << "candidate";
    1205     if (a.pointlist.size() != 1)
    1206       ost << "s ";
    1207     else
    1208       ost << " ";
    1209     for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)
    1210       ost << *(*Runner) << " ";
    1211     ost << " at angle " << (a.ShortestAngle) << ".";
    1212   }
    1213 
    1214   return ost;
    1215 }
    1216 ;
    1217 
    1218 // =========================================================== class TESSELATION ===========================================
    1219 
    122037/** Constructor of class Tesselation.
    122138 */
    … …  
    125471  int num = 0;
    125572  for (GoToFirst(); (!IsEnd()); GoToNext()) {
    1256     (*Center) += (*GetPoint()->node);
     73    (*Center) += (GetPoint()->getPosition());
    125774    num++;
    125875  }
    … …  
    1337154      C++;
    1338155      for (; C != PointsOnBoundary.end(); C++) {
    1339         tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
     156        tmp = A->second->node->DistanceSquared(B->second->node->getPosition());
    1340157        distance = tmp * tmp;
    1341         tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
     158        tmp = A->second->node->DistanceSquared(C->second->node->getPosition());
    1342159        distance += tmp * tmp;
    1343         tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
     160        tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
    1344161        distance += tmp * tmp;
    1345162        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
    … …  
    1360177    // 2. next, we have to check whether all points reside on only one side of the triangle
    1361178    // 3. construct plane vector
    1362     PlaneVector = Plane(*A->second->node->node,
    1363                         *baseline->second.first->second->node->node,
    1364                         *baseline->second.second->second->node->node).getNormal();
     179    PlaneVector = Plane(A->second->node->getPosition(),
     180                        baseline->second.first->second->node->getPosition(),
     181                        baseline->second.second->second->node->getPosition()).getNormal();
    1365182    DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
    1366183    // 4. loop over all points
    … …  
    1372189        continue;
    1373190      // 4a. project onto plane vector
    1374       TrialVector = (*checker->second->node->node);
    1375       TrialVector.SubtractVector(*A->second->node->node);
     191      TrialVector = (checker->second->node->getPosition() - A->second->node->getPosition());
    1376192      distance = TrialVector.ScalarProduct(PlaneVector);
    1377193      if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
    … …  
    1389205      }
    1390206      // 4d. Check whether the point is inside the triangle (check distance to each node
    1391       tmp = checker->second->node->node->DistanceSquared(*A->second->node->node);
     207      tmp = checker->second->node->DistanceSquared(A->second->node->getPosition());
    1392208      int innerpoint = 0;
    1393       if ((tmp < A->second->node->node->DistanceSquared(*baseline->second.first->second->node->node)) && (tmp < A->second->node->node->DistanceSquared(*baseline->second.second->second->node->node)))
     209      if ((tmp < A->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < A->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
    1394210        innerpoint++;
    1395       tmp = checker->second->node->node->DistanceSquared(*baseline->second.first->second->node->node);
    1396       if ((tmp < baseline->second.first->second->node->node->DistanceSquared(*A->second->node->node)) && (tmp < baseline->second.first->second->node->node->DistanceSquared(*baseline->second.second->second->node->node)))
     211      tmp = checker->second->node->DistanceSquared(baseline->second.first->second->node->getPosition());
     212      if ((tmp < baseline->second.first->second->node->DistanceSquared(A->second->node->getPosition())) && (tmp < baseline->second.first->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
    1397213        innerpoint++;
    1398       tmp = checker->second->node->node->DistanceSquared(*baseline->second.second->second->node->node);
    1399       if ((tmp < baseline->second.second->second->node->node->DistanceSquared(*baseline->second.first->second->node->node)) && (tmp < baseline->second.second->second->node->node->DistanceSquared(*A->second->node->node)))
     214      tmp = checker->second->node->DistanceSquared(baseline->second.second->second->node->getPosition());
     215      if ((tmp < baseline->second.second->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < baseline->second.second->second->node->DistanceSquared(A->second->node->getPosition())))
    1400216        innerpoint++;
    1401217      // 4e. If so, break 4. loop and continue with next candidate in 1. loop
    … …  
    1478294        // prepare some auxiliary vectors
    1479295        Vector BaseLineCenter, BaseLine;
    1480         BaseLineCenter = 0.5 * ((*baseline->second->endpoints[0]->node->node) +
    1481                                 (*baseline->second->endpoints[1]->node->node));
    1482         BaseLine = (*baseline->second->endpoints[0]->node->node) - (*baseline->second->endpoints[1]->node->node);
     296        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) +
     297                                (baseline->second->endpoints[1]->node->getPosition()));
     298        BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition());
    1483299
    1484300        // offset to center of triangle
    … …  
    1498314        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
    1499315        PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
    1500         TempVector = CenterVector - (*baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
     316        TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
    1501317        //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
    1502318        if (PropagationVector.ScalarProduct(TempVector) > 0) // make sure normal propagation vector points outward from baseline
    … …  
    1511327
    1512328            // first check direction, so that triangles don't intersect
    1513             VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
     329            VirtualNormalVector = (target->second->node->getPosition()) - BaseLineCenter;
    1514330            VirtualNormalVector.ProjectOntoPlane(NormalVector);
    1515331            TempAngle = VirtualNormalVector.Angle(PropagationVector);
    … …  
    1540356
    1541357            // check for linear dependence
    1542             TempVector = (*baseline->second->endpoints[0]->node->node) - (*target->second->node->node);
    1543             helper = (*baseline->second->endpoints[1]->node->node) - (*target->second->node->node);
     358            TempVector = (baseline->second->endpoints[0]->node->getPosition()) - (target->second->node->getPosition());
     359            helper = (baseline->second->endpoints[1]->node->getPosition()) - (target->second->node->getPosition());
    1544360            helper.ProjectOntoPlane(TempVector);
    1545361            if (fabs(helper.NormSquared()) < MYEPSILON) {
    … …  
    1550366            // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
    1551367            flag = true;
    1552             VirtualNormalVector = Plane(*(baseline->second->endpoints[0]->node->node),
    1553                                         *(baseline->second->endpoints[1]->node->node),
    1554                                         *(target->second->node->node)).getNormal();
    1555             TempVector = (1./3.) * ((*baseline->second->endpoints[0]->node->node) +
    1556                                     (*baseline->second->endpoints[1]->node->node) +
    1557                                     (*target->second->node->node));
     368            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()),
     369                                        (baseline->second->endpoints[1]->node->getPosition()),
     370                                        (target->second->node->getPosition())).getNormal();
     371            TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) +
     372                                    (baseline->second->endpoints[1]->node->getPosition()) +
     373                                    (target->second->node->getPosition()));
    1558374            TempVector -= (*Center);
    1559375            // make it always point outward
    … …  
    1569385            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
    1570386              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
    1571               helper = (*target->second->node->node) - BaseLineCenter;
     387              helper = (target->second->node->getPosition()) - BaseLineCenter;
    1572388              helper.ProjectOntoPlane(BaseLine);
    1573389              // ...the one with the smaller angle is the better candidate
    1574               TempVector = (*target->second->node->node) - BaseLineCenter;
     390              TempVector = (target->second->node->getPosition()) - BaseLineCenter;
    1575391              TempVector.ProjectOntoPlane(VirtualNormalVector);
    1576392              TempAngle = TempVector.Angle(helper);
    1577               TempVector = (*winner->second->node->node) - BaseLineCenter;
     393              TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
    1578394              TempVector.ProjectOntoPlane(VirtualNormalVector);
    1579395              if (TempAngle < TempVector.Angle(helper)) {
    … …  
    1614430            BLS[2] = LineChecker[1]->second;
    1615431          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    1616           BTS->GetCenter(&helper);
     432          BTS->GetCenter(helper);
    1617433          helper -= (*Center);
    1618434          helper *= -1;
    … …  
    1662478    DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
    1663479    // get the next triangle
    1664     triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     480    triangles = FindClosestTrianglesToVector(Walker->getPosition(), BoundaryPoints);
    1665481    if (triangles != NULL)
    1666482      BTS = triangles->front();
    … …  
    1676492    DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
    1677493    // get the intersection point
    1678     if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) {
     494    if (BTS->GetIntersectionInsideTriangle(*Center, Walker->getPosition(), Intersection)) {
    1679495      DoLog(0) && (Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl);
    1680496      // we have the intersection, check whether in- or outside of boundary
    1681       if ((Center->DistanceSquared(*Walker->node) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
     497      if ((Center->DistanceSquared(Walker->getPosition()) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
    1682498        // inside, next!
    1683499        DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
    … …  
    2086902  DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
    2087903  for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    2088     DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
     904    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
    2089905
    2090906  // remove triangles's endpoints
    … …  
    2102918    DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
    2103919    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    2104       DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
     920      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
    2105921  }
    2106922  delete (ListofPoints);
    … …  
    22571073        if (List != NULL) {
    22581074          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    2259             if ((*Runner)->node->at(map[0]) > maxCoordinate[map[0]]) {
    2260               DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);
    2261               maxCoordinate[map[0]] = (*Runner)->node->at(map[0]);
     1075            if ((*Runner)->at(map[0]) > maxCoordinate[map[0]]) {
     1076              DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << (*Runner)->getPosition() << "." << endl);
     1077              maxCoordinate[map[0]] = (*Runner)->at(map[0]);
    22621078              MaxPoint[map[0]] = (*Runner);
    22631079            }
    … …  
    22951111
    22961112    // construct center of circle
    2297     CircleCenter = 0.5 * ((*BaseLine->endpoints[0]->node->node) + (*BaseLine->endpoints[1]->node->node));
     1113    CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
    22981114
    22991115    // construct normal vector of circle
    2300     CirclePlaneNormal = (*BaseLine->endpoints[0]->node->node) - (*BaseLine->endpoints[1]->node->node);
     1116    CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
    23011117
    23021118    double radius = CirclePlaneNormal.NormSquared();
    … …  
    25151331
    25161332  // construct center of circle
    2517   CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
    2518                         (*CandidateLine.BaseLine->endpoints[1]->node->node));
     1333  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
     1334                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
    25191335
    25201336  // construct normal vector of circle
    2521   CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
    2522                       (*CandidateLine.BaseLine->endpoints[1]->node->node);
     1337  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
     1338                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
    25231339
    25241340  // calculate squared radius of circle
    … …  
    25361352    // construct SearchDirection and an "outward pointer"
    25371353    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
    2538     helper = CircleCenter - (*ThirdPoint->node->node);
     1354    helper = CircleCenter - (ThirdPoint->node->getPosition());
    25391355    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
    25401356      SearchDirection.Scale(-1.);
    … …  
    26111427  for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
    26121428    SetOfNeighbours.insert(*Runner);
    2613   TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
     1429  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition());
    26141430
    26151431  DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
    … …  
    27121528
    27131529  // create normal vector
    2714   BTS->GetCenter(&Center);
     1530  BTS->GetCenter(Center);
    27151531  Center -= CandidateLine.OptCenter;
    27161532  BTS->SphereCenter = CandidateLine.OptCenter;
    … …  
    27911607
    27921608  // create normal vector
    2793   BTS->GetCenter(&Center);
     1609  BTS->GetCenter(Center);
    27941610  Center.SubtractVector(*OptCenter);
    27951611  BTS->SphereCenter = *OptCenter;
    … …  
    28371653  Vector DistanceToIntersection[2], BaseLine;
    28381654  double distance[2];
    2839   BaseLine = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
     1655  BaseLine = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
    28401656  for (int i = 0; i < 2; i++) {
    2841     DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
     1657    DistanceToIntersection[i] = (*ClosestPoint) - (Base->endpoints[i]->node->getPosition());
    28421658    distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
    28431659  }
    … …  
    29191735
    29201736  // calculate volume
    2921   volume = CalculateVolumeofGeneralTetraeder(*Base->endpoints[1]->node->node, *OtherBase->endpoints[0]->node->node, *OtherBase->endpoints[1]->node->node, *Base->endpoints[0]->node->node);
     1737  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->getPosition(), OtherBase->endpoints[0]->node->getPosition(), OtherBase->endpoints[1]->node->getPosition(), Base->endpoints[0]->node->getPosition());
    29221738
    29231739  // delete the temporary other base line and the closest points
    … …  
    31251941
    31261942              OrthogonalizedOben = Oben;
    3127               aCandidate = (*a->node) - (*Candidate->node);
     1943              aCandidate = (a->getPosition()) - (Candidate->getPosition());
    31281944              OrthogonalizedOben.ProjectOntoPlane(aCandidate);
    31291945              OrthogonalizedOben.Normalize();
    … …  
    31321948              OrthogonalizedOben.Scale(scaleFactor);
    31331949
    3134               Center = 0.5 * ((*Candidate->node) + (*a->node));
     1950              Center = 0.5 * ((Candidate->getPosition()) + (a->getPosition()));
    31351951              Center += OrthogonalizedOben;
    31361952
    3137               AngleCheck = Center - (*a->node);
     1953              AngleCheck = Center - (a->getPosition());
    31381954              norm = aCandidate.Norm();
    31391955              // second point shall have smallest angle with respect to Oben vector
    … …  
    32202036
    32212037  // construct center of circle
    3222   CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
    3223                         (*CandidateLine.BaseLine->endpoints[1]->node->node));
     2038  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
     2039                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
    32242040
    32252041  // construct normal vector of circle
    3226   CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
    3227                       (*CandidateLine.BaseLine->endpoints[1]->node->node);
     2042  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
     2043                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
    32282044
    32292045  RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
    … …  
    32522068
    32532069      // get cell for the starting point
    3254       if (LC->SetIndexToVector(&CircleCenter)) {
     2070      if (LC->SetIndexToVector(CircleCenter)) {
    32552071        for (int i = 0; i < NDIM; i++) // store indices of this cell
    32562072          N[i] = LC->n[i];
    … …  
    32822098
    32832099                  // find center on the plane
    3284                   GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
     2100                  GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
    32852101                  DoLog(1) && (Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl);
    32862102
    32872103                  try {
    3288                     NewNormalVector = Plane(*(CandidateLine.BaseLine->endpoints[0]->node->node),
    3289                                             *(CandidateLine.BaseLine->endpoints[1]->node->node),
    3290                                             *(Candidate->node)).getNormal();
     2104                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
     2105                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
     2106                                            (Candidate->getPosition())).getNormal();
    32912107                    DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
    3292                     radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(NewPlaneCenter);
     2108                    radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter);
    32932109                    DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
    32942110                    DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
    32952111                    DoLog(1) && (Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl);
    32962112                    if (radius < RADIUS * RADIUS) {
    3297                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(NewPlaneCenter);
     2113                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
    32982114                      if (fabs(radius - otherradius) < HULLEPSILON) {
    32992115                        // construct both new centers
    … …  
    34222238 * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
    34232239 */
    3424 DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
     2240DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector &x, const LinkedCell* LC) const
    34252241{
    34262242  Info FunctionInfo(__func__);
    … …  
    34502266            FindPoint = PointsOnBoundary.find((*Runner)->nr);
    34512267            if (FindPoint != PointsOnBoundary.end()) {
    3452               points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(*x), FindPoint->second));
     2268              points->insert(DistanceToPointPair(FindPoint->second->node->DistanceSquared(x), FindPoint->second));
    34532269              DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
    34542270            }
    … …  
    34742290 * \return closest BoundaryLineSet or NULL in degenerate case.
    34752291 */
    3476 BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
     2292BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector &x, const LinkedCell* LC) const
    34772293{
    34782294  Info FunctionInfo(__func__);
    … …  
    34852301
    34862302  // for each point, check its lines, remember closest
    3487   DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl);
     2303  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << x << " ... " << endl);
    34882304  BoundaryLineSet *ClosestLine = NULL;
    34892305  double MinDistance = -1.;
    … …  
    34942310    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
    34952311      // calculate closest point on line to desired point
    3496       helper = 0.5 * ((*(LineRunner->second)->endpoints[0]->node->node) +
    3497                       (*(LineRunner->second)->endpoints[1]->node->node));
    3498       Center = (*x) - helper;
    3499       BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
    3500                  (*(LineRunner->second)->endpoints[1]->node->node);
     2312      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) +
     2313                      ((LineRunner->second)->endpoints[1]->node->getPosition()));
     2314      Center = (x) - helper;
     2315      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
     2316                 ((LineRunner->second)->endpoints[1]->node->getPosition());
    35012317      Center.ProjectOntoPlane(BaseLine);
    35022318      const double distance = Center.NormSquared();
    35032319      if ((ClosestLine == NULL) || (distance < MinDistance)) {
    35042320        // additionally calculate intersection on line (whether it's on the line section or not)
    3505         helper = (*x) - (*(LineRunner->second)->endpoints[0]->node->node) - Center;
     2321        helper = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition()) - Center;
    35062322        const double lengthA = helper.ScalarProduct(BaseLine);
    3507         helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
     2323        helper = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition()) - Center;
    35082324        const double lengthB = helper.ScalarProduct(BaseLine);
    35092325        if (lengthB * lengthA < 0) { // if have different sign
    … …  
    35342350 * \return BoundaryTriangleSet of nearest triangle or NULL.
    35352351 */
    3536 TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
     2352TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector &x, const LinkedCell* LC) const
    35372353{
    35382354  Info FunctionInfo(__func__);
    … …  
    35452361
    35462362  // for each point, check its lines, remember closest
    3547   DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl);
     2363  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << x << " ... " << endl);
    35482364  LineSet ClosestLines;
    35492365  double MinDistance = 1e+16;
    … …  
    35552371    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
    35562372
    3557       BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
    3558                  (*(LineRunner->second)->endpoints[1]->node->node);
     2373      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
     2374                 ((LineRunner->second)->endpoints[1]->node->getPosition());
    35592375      const double lengthBase = BaseLine.NormSquared();
    35602376
    3561       BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
     2377      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition());
    35622378      const double lengthEndA = BaseLineIntersection.NormSquared();
    35632379
    3564       BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
     2380      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
    35652381      const double lengthEndB = BaseLineIntersection.NormSquared();
    35662382
    … …  
    35802396      } else { // intersection is closer, calculate
    35812397        // calculate closest point on line to desired point
    3582         BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
     2398        BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
    35832399        Center = BaseLineIntersection;
    35842400        Center.ProjectOntoPlane(BaseLine);
    … …  
    36212437 * \return list of BoundaryTriangleSet of nearest triangles or NULL.
    36222438 */
    3623 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
     2439class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector &x, const LinkedCell* LC) const
    36242440{
    36252441  Info FunctionInfo(__func__);
    … …  
    36362452  double MinAlignment = 2. * M_PI;
    36372453  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
    3638     (*Runner)->GetCenter(&Center);
    3639     helper = (*x) - Center;
     2454    (*Runner)->GetCenter(Center);
     2455    helper = (x) - Center;
    36402456    const double Alignment = helper.Angle((*Runner)->NormalVector);
    36412457    if (Alignment < MinAlignment) {
    … …  
    36632479{
    36642480  Info FunctionInfo(__func__);
    3665   TriangleIntersectionList Intersections(&Point, this, LC);
     2481  TriangleIntersectionList Intersections(Point, this, LC);
    36662482
    36672483  return Intersections.IsInside();
    … …  
    37032519  }
    37042520
    3705   triangle->GetCenter(&Center);
     2521  triangle->GetCenter(Center);
    37062522  DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
    37072523  DistanceToCenter = Center - Point;
    … …  
    37142530    Center = Point - triangle->NormalVector; // points towards MolCenter
    37152531    DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);
    3716     if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
     2532    if (triangle->GetIntersectionInsideTriangle(Center, DistanceToCenter, Intersection)) {
    37172533      DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);
    37182534      return 0.;
    … …  
    37232539  } else {
    37242540    // calculate smallest distance
    3725     distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
     2541    distance = triangle->GetClosestPointInsideTriangle(Point, Intersection);
    37262542    DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);
    37272543
    … …  
    37472563{
    37482564  Info FunctionInfo(__func__);
    3749   TriangleIntersectionList Intersections(&Point, this, LC);
     2565  TriangleIntersectionList Intersections(Point, this, LC);
    37502566
    37512567  return Intersections.GetSmallestDistance();
    … …  
    37622578{
    37632579  Info FunctionInfo(__func__);
    3764   TriangleIntersectionList Intersections(&Point, this, LC);
     2580  TriangleIntersectionList Intersections(Point, this, LC);
    37652581
    37662582  return Intersections.GetClosestTriangle();
    … …  
    38392655 * @return list of the all points linked to the provided one
    38402656 */
    3841 TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
     2657TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
    38422658{
    38432659  Info FunctionInfo(__func__);
    … …  
    38712687
    38722688  // construct one orthogonal vector
    3873   if (Reference != NULL) {
    3874     AngleZero = (*Reference) - (*Point->node);
    3875     AngleZero.ProjectOntoPlane(PlaneNormal);
    3876   }
    3877   if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
    3878     DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
    3879     AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
     2689  AngleZero = (Reference) - (Point->getPosition());
     2690  AngleZero.ProjectOntoPlane(PlaneNormal);
     2691  if ((AngleZero.NormSquared() < MYEPSILON)) {
     2692    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
     2693    AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
    38802694    AngleZero.ProjectOntoPlane(PlaneNormal);
    38812695    if (AngleZero.NormSquared() < MYEPSILON) {
    … …  
    38932707  // go through all connected points and calculate angle
    38942708  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
    3895     helper = (*(*listRunner)->node) - (*Point->node);
     2709    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
    38962710    helper.ProjectOntoPlane(PlaneNormal);
    38972711    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
    … …  
    39152729 * @param *SetOfNeighbours all points for which the angle should be calculated
    39162730 * @param *Point of which get all connected points
    3917  * @param *Reference Reference vector for zero angle or NULL for no preference
     2731 * @param *Reference Reference vector for zero angle or (0,0,0) for no preference
    39182732 * @return list of the all points linked to the provided one
    39192733 */
    3920 TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
     2734TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
    39212735{
    39222736  Info FunctionInfo(__func__);
    … …  
    39422756  }
    39432757
    3944   DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl);
     2758  DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << Reference << "." << endl);
    39452759  // calculate central point
    39462760  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
    … …  
    39522766  int counter = 0;
    39532767  while (TesselC != SetOfNeighbours->end()) {
    3954     helper = Plane(*((*TesselA)->node),
    3955                    *((*TesselB)->node),
    3956                    *((*TesselC)->node)).getNormal();
     2768    helper = Plane(((*TesselA)->getPosition()),
     2769                   ((*TesselB)->getPosition()),
     2770                   ((*TesselC)->getPosition())).getNormal();
    39572771    DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);
    39582772    counter++;
    … …  
    39742788
    39752789  // construct one orthogonal vector
    3976   if (Reference != NULL) {
    3977     AngleZero = (*Reference) - (*Point->node);
     2790  if (!Reference.IsZero()) {
     2791    AngleZero = (Reference) - (Point->getPosition());
    39782792    AngleZero.ProjectOntoPlane(PlaneNormal);
    39792793  }
    3980   if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
    3981     DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
    3982     AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
     2794  if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON )) {
     2795    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
     2796    AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
    39832797    AngleZero.ProjectOntoPlane(PlaneNormal);
    39842798    if (AngleZero.NormSquared() < MYEPSILON) {
    … …  
    39972811  pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
    39982812  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
    3999     helper = (*(*listRunner)->node) - (*Point->node);
     2813    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
    40002814    helper.ProjectOntoPlane(PlaneNormal);
    40012815    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
    … …  
    42423056
    42433057  // copy old location for the volume
    4244   OldPoint = (*point->node->node);
     3058  OldPoint = (point->node->getPosition());
    42453059
    42463060  // get list of connected points
    … …  
    43053119          StartNode--;
    43063120          //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
    4307           Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
     3121          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
    43083122          StartNode = MiddleNode;
    43093123          StartNode++;
    … …  
    43113125            StartNode = connectedPath->begin();
    43123126          //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
    4313           Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
    4314           OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
     3127          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
     3128          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
    43153129          OrthogonalVector.MakeNormalTo(Reference);
    43163130          angle = GetAngle(Point, Reference, OrthogonalVector);
    … …  
    43673181        AddTesselationTriangle();
    43683182        // calculate volume summand as a general tetraeder
    4369         volume += CalculateVolumeofGeneralTetraeder(*TPS[0]->node->node, *TPS[1]->node->node, *TPS[2]->node->node, OldPoint);
     3183        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
    43703184        // advance number
    43713185        count++;
    … …  
    47523566  // find nearest boundary point
    47533567  class TesselPoint *BackupPoint = NULL;
    4754   class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
     3568  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->getPosition(), BackupPoint, LC);
    47553569  class BoundaryPointSet *NearestBoundaryPoint = NULL;
    47563570  PointMap::iterator PointRunner;
    … …  
    47733587  class BoundaryLineSet *BestLine = NULL;
    47743588  for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
    4775     BaseLine = (*Runner->second->endpoints[0]->node->node) -
    4776                (*Runner->second->endpoints[1]->node->node);
    4777     CenterToPoint = 0.5 * ((*Runner->second->endpoints[0]->node->node) +
    4778                            (*Runner->second->endpoints[1]->node->node));
    4779     CenterToPoint -= (*point->node);
     3589    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) -
     3590               (Runner->second->endpoints[1]->node->getPosition());
     3591    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) +
     3592                           (Runner->second->endpoints[1]->node->getPosition()));
     3593    CenterToPoint -= (point->getPosition());
    47803594    angle = CenterToPoint.Angle(BaseLine);
    47813595    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
Note: See TracChangeset for help on using the changeset viewer.