Changes in / [124e14:76102e]

Location:

src

Files:

• boundary.cpp (modified) (3 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (2 diffs)
• tesselation.cpp (modified) (29 diffs)
• tesselation.hpp (modified) (4 diffs)
• unittests/Makefile.am (modified) (2 diffs)
• unittests/tesselation_boundarytriangleunittest.cpp (deleted)
• unittests/tesselation_boundarytriangleunittest.hpp (deleted)
• vector.cpp (modified) (5 diffs)

src/boundary.cpp

@@ -789 +789 @@
  * \param configuration contains box dimensions
  * \param distance[NDIM] distance between filling molecules in each direction
- * \param boundary length of boundary zone between molecule and filling mollecules
  * \param epsilon distance to surface which is not filled
  * \param RandAtomDisplacement maximum distance for random displacement per atom
@@ -796 +795 @@
  * \return *mol pointer to new molecule with filled atoms
  */
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
 {
         Info FunctionInfo(__func__);
@@ -857 +856 @@
             FillIt = false;
           } else {
-            const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
-            FillIt = FillIt && (distance > boundary*boundary);
-            if (FillIt) {
+            const bool FillResult = (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
+            FillIt = FillIt && FillResult;
+            if (FillResult) {
               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
             } else {
-              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
-              break;
+              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point." << endl;
             }
             i++;

src/boundary.hpp

@@ -49 +49 @@
 
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
 void FindConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);

src/builder.cpp

@@ -1738 +1738 @@
               if (argptr+6 >=argc) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
+                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
                 performCriticalExit();
               } else {
@@ -1769 +1769 @@
                 for (int i=0;i<NDIM;i++)
                   distance[i] = atof(argv[argptr+i]);
-                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
+                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
                 if (Filling != NULL) {
                   Filling->ActiveFlag = false;

src/tesselation.cpp

@@ -35 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
   LinesCount(0),
   node(Walker),
@@ -61 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -105 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
 {
         Info FunctionInfo(__func__);
@@ -115 +115 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
-  // set skipped to false
-  skipped = false;
-  // clear triangles list
-  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
-};
-
-/** Constructor of BoundaryLineSet with two endpoints.
- * Adds line automatically to each endpoints' LineMap
- * \param *Point1 first boundary point
- * \param *Point2 second boundary point
- * \param number number of the list
- */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number)
-{
-  Info FunctionInfo(__func__);
-  // set number
-  Nr = number;
-  // set endpoints in ascending order
-  SetEndpointsOrdered(endpoints, Point1, Point2);
-  // add this line to the hash maps of both endpoints
-  Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
-  Point2->AddLine(this); //
   // set skipped to false
   skipped = false;
@@ -193 +171 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
         Info FunctionInfo(__func__);
@@ -204 +182 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -219 +197 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion() const
+bool BoundaryLineSet::CheckConvexityCriterion()
 {
         Info FunctionInfo(__func__);
@@ -243 +221 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
@@ -286 +264 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -299 +277 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -339 +317 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
   Nr(number)
 {
@@ -398 +376 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
         Info FunctionInfo(__func__);
@@ -412 +390 @@
 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
@@ -422 +400 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
 {
   Info FunctionInfo(__func__);
@@ -474 +452 @@
 };
 
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
- * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
- * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
- * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
- * the first two basepoints) or not.
- * \param *x point
- * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
- * \return Distance squared between \a *x and closest point inside triangle
- */
-double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
-{
-  Info FunctionInfo(__func__);
-  Vector Direction;
-
-  // 1. get intersection with plane
-  Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
-  GetCenter(&Direction);
-  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
-    ClosestPoint->CopyVector(x);
-  }
-
-  // 2. Calculate in plane part of line (x, intersection)
-  Vector InPlane;
-  InPlane.CopyVector(x);
-  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
-  InPlane.ProjectOntoPlane(&NormalVector);
-  InPlane.AddVector(ClosestPoint);
-
-  Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
-  Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl;
-  Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl;
-
-  // Calculate cross point between one baseline and the desired point such that distance is shortest
-  double ShortestDistance = -1.;
-  bool InsideFlag = false;
-  Vector CrossDirection[3];
-  Vector CrossPoint[3];
-  Vector helper;
-  for (int i=0;i<3;i++) {
-    // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
-    Direction.CopyVector(endpoints[(i+1)%3]->node->node);
-    Direction.SubtractVector(endpoints[i%3]->node->node);
-    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
-    CrossDirection[i].CopyVector(&CrossPoint[i]);
-    CrossDirection[i].SubtractVector(&InPlane);
-    CrossPoint[i].SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-    const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared();
-    Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl;
-    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
-      CrossPoint[i].AddVector(endpoints[i%3]->node->node);  // make cross point absolute again
-      Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl;
-      const double distance = CrossPoint[i].DistanceSquared(x);
-      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
-        ShortestDistance = distance;
-        ClosestPoint->CopyVector(&CrossPoint[i]);
-      }
-    } else
-      CrossPoint[i].Zero();
-  }
-  InsideFlag = true;
-  for (int i=0;i<3;i++) {
-    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]);
-    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);;
-    if ((sign > -MYEPSILON) && (othersign > -MYEPSILON))  // have different sign
-      InsideFlag = false;
-  }
-  if (InsideFlag) {
-    ClosestPoint->CopyVector(&InPlane);
-    ShortestDistance = InPlane.DistanceSquared(x);
-  } else {  // also check endnodes
-    for (int i=0;i<3;i++) {
-      const double distance = x->DistanceSquared(endpoints[i]->node->node);
-      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
-        ShortestDistance = distance;
-        ClosestPoint->CopyVector(endpoints[i]->node->node);
-      }
-    }
-  }
-  Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
-  return ShortestDistance;
-};
-
 /** Checks whether lines is any of the three boundary lines this triangle contains.
  * \param *line line to test
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -575 +469 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -588 +482 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
 {
         Info FunctionInfo(__func__);
@@ -601 +495 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
 {
         Info FunctionInfo(__func__);
@@ -624 +518 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
 {
         Info FunctionInfo(__func__);
@@ -646 +540 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -663 +557 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector * const center) const
+void BoundaryTriangleSet::GetCenter(Vector *center)
 {
         Info FunctionInfo(__func__);
@@ -670 +564 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
-  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -3242 +3135 @@
   // for each point, check its lines, remember closest
   Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
-  LineSet ClosestLines;
+  BoundaryLineSet *ClosestLine = NULL;
   double MinDistance = 1e+16;
   Vector BaseLineIntersection;
@@ -3264 +3157 @@
 
       if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
-        const double lengthEnd = Min(lengthEndA, lengthEndB);
-        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
-          ClosestLines.clear();
-          ClosestLines.insert(LineRunner->second);
-          MinDistance = lengthEnd;
-          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
-        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
-          ClosestLines.insert(LineRunner->second);
-          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
-        } else { // line is worse
+        if ((lengthEndA < MinDistance) && (lengthEndA < lengthEndB)) {
+          ClosestLine = LineRunner->second;
+          MinDistance = lengthEndA;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEndA << "." << endl;
+        } else if (lengthEndB < MinDistance) {
+          ClosestLine = LineRunner->second;
+          MinDistance = lengthEndB;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is closer with " << lengthEndB << "." << endl;
+        } else {
           Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
         }
@@ -3287 +3179 @@
           eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
         }
-        if ((ClosestLines.empty()) || (distance < MinDistance)) {
-          ClosestLines.insert(LineRunner->second);
+        if ((ClosestLine == NULL) || (distance < MinDistance)) {
+          ClosestLine = LineRunner->second;
           MinDistance = distance;
-          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
         } else {
           Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
@@ -3300 +3192 @@
 
   // check whether closest line is "too close" :), then it's inside
-  if (ClosestLines.empty()) {
+  if (ClosestLine == NULL) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
   TriangleList * candidates = new TriangleList;
-  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
-    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
+  for (TriangleMap::iterator Runner = ClosestLine->triangles.begin(); Runner != ClosestLine->triangles.end(); Runner++) {
     candidates->push_back(Runner->second);
   }
@@ -3350 +3241 @@
 };
 
-/** Checks whether the provided Vector is within the Tesselation structure.
- * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
- * @param point of which to check the position
- * @param *LC LinkedCell structure
- *
- * @return true if the point is inside the Tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
-{
-  return (GetDistanceSquaredToSurface(Point, LC) > MYEPSILON);
-}
-
-/** Returns the distance to the surface given by the tesselation.
+/** Checks whether the provided Vector is within the tesselation structure.
  * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
- * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
- * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
- * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
- * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
- *  -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
- *  -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
- *  -# If inside, take it to calculate closest distance
- *  -# If not, take intersection with BoundaryLine as distance
- *
- * @note distance is squared despite it still contains a sign to determine in-/outside!
+ * towards or away from the given \a &Point.
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @return >0 if outside, ==0 if on surface, <0 if inside (Note that distance can be at most LinkedCell::RADIUS.)
- */
-double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
 {
   Info FunctionInfo(__func__);
@@ -3388 +3258 @@
   Vector DistanceToCenter;
   Vector Intersection;
-  double distance = 0.;
 
   if (result == NULL) {// is boundary point or only point in point cloud?
     Log() << Verbose(1) << Point << " is the only point in vicinity." << endl;
-    return LC->RADIUS;
+    return false;
   } else {
     Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
@@ -3413 +3282 @@
     if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
       Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
-      return 0.;
+      return true;
     } else {
       Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
       return false;
     }
+  }
+
+  // then check direction to boundary
+  if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
+    Log() << Verbose(1) << Point << " is an inner point." << endl;
+    return true;
   } else {
-    // calculate smallest distance
-    distance = result->GetClosestPointInsideTriangle(&Point, &Intersection);
-    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
-    distance = Min(distance, (LC->RADIUS*LC->RADIUS));
-
-    // then check direction to boundary
-    if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
-      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
-      return -distance;
-    } else {
-      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
-      return +distance;
-    }
-  }
-};
+    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
+    return false;
+  }
+}
 
 /** Gets all points connected to the provided point by triangulation lines.

src/tesselation.hpp

@@ -106 +106 @@
   public:
     BoundaryPointSet();
-    BoundaryPointSet(TesselPoint * const Walker);
+    BoundaryPointSet(TesselPoint * Walker);
     ~BoundaryPointSet();
 
-    void AddLine(BoundaryLineSet * const line);
+    void AddLine(class BoundaryLineSet *line);
 
     LineMap lines;
@@ -125 +125 @@
   public:
     BoundaryLineSet();
-    BoundaryLineSet(BoundaryPointSet * const Point[2], const int number);
-    BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number);
+    BoundaryLineSet(class BoundaryPointSet *Point[2], const int number);
     ~BoundaryLineSet();
 
-    void AddTriangle(BoundaryTriangleSet * const triangle);
-    bool IsConnectedTo(const BoundaryLineSet * const line) const;
-    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
-    bool CheckConvexityCriterion() const;
-    class BoundaryPointSet *GetOtherEndpoint(const BoundaryPointSet * const point) const;
+    void AddTriangle(class BoundaryTriangleSet *triangle);
+    bool IsConnectedTo(class BoundaryLineSet *line);
+    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
+    bool CheckConvexityCriterion();
+    class BoundaryPointSet *GetOtherEndpoint(class BoundaryPointSet *);
 
     class BoundaryPointSet *endpoints[2];
@@ -148 +147 @@
   public:
     BoundaryTriangleSet();
-    BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number);
+    BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
     ~BoundaryTriangleSet();
 
-    void GetNormalVector(const Vector &NormalVector);
-    void GetCenter(Vector * const center) const;
-    bool GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const;
-    double GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const;
-    bool ContainsBoundaryLine(const BoundaryLineSet * const line) const;
-    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
-    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
-    class BoundaryPointSet *GetThirdEndpoint(const BoundaryLineSet * const line) const;
-    bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
-    bool IsPresentTupel(const BoundaryTriangleSet * const T) const;
+    void GetNormalVector(Vector &NormalVector);
+    void GetCenter(Vector *center);
+    bool GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection);
+    bool ContainsBoundaryLine(class BoundaryLineSet *line);
+    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
+    bool ContainsBoundaryPoint(class TesselPoint *point);
+    class BoundaryPointSet *GetThirdEndpoint(class BoundaryLineSet *line);
+    bool IsPresentTupel(class BoundaryPointSet *Points[3]);
+    bool IsPresentTupel(class BoundaryTriangleSet *T);
 
     class BoundaryPointSet *endpoints[3];
@@ -315 +313 @@
     class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
     bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
-    double GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
     DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;

src/unittests/Makefile.am

@@ -22 +22 @@
   StackClassUnitTest \
   TesselationUnitTest \
-  Tesselation_BoundaryTriangleUnitTest \
-  Tesselation_InOutsideUnitTest \
+  TesselationInOutsideUnitTest \
   VectorUnitTest
   
@@ -80 +79 @@
 TesselationUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
-Tesselation_BoundaryTriangleUnitTest_SOURCES = tesselation_boundarytriangleunittest.cpp tesselation_boundarytriangleunittest.hpp
-Tesselation_BoundaryTriangleUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
-
-Tesselation_InOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
-Tesselation_InOutsideUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+TesselationInOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
+TesselationInOutsideUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 VectorUnitTest_SOURCES = vectorunittest.cpp vectorunittest.hpp 

src/vector.cpp

@@ -222 +222 @@
  * \param *Origin first vector of line
  * \param *LineVector second vector of line
- * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
+ * \return true -  \a this contains intersection point on return, false - line is parallel to plane
  */
 bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
@@ -235 +235 @@
   Direction.Normalize();
   Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
-  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
-  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
-    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
+  if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
+    eLog() << Verbose(2) << "Line is parallel to plane, no intersection." << endl;
     return false;
   }
@@ -244 +243 @@
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
-  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
-    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
+  if (factor < MYEPSILON) { // Origin is in-plane
+    Log() << Verbose(1) << "Origin of line is in-plane, simple." << endl;
     CopyVector(Origin);
     return true;
@@ -259 +258 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
+    Log() << Verbose(1) << "INFO: Intersection at " << *this << " is good." << endl;
     return true;
   } else {
@@ -354 +353 @@
   Vector parallel;
   double factor = 0.;
+  double pfactor = 0.;
   if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
     parallel.CopyVector(Line1a);