Changes in / [76102e:124e14]

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 (added)
• unittests/tesselation_boundarytriangleunittest.hpp (added)
• 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
@@ -795 +796 @@
  * \return *mol pointer to new molecule with filled atoms
  */
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
 {
         Info FunctionInfo(__func__);
@@ -856 +857 @@
             FillIt = false;
           } else {
-            const bool FillResult = (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
-            FillIt = FillIt && FillResult;
-            if (FillResult) {
+            const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
+            FillIt = FillIt && (distance > boundary*boundary);
+            if (FillIt) {
               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
             } else {
-              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point." << endl;
+              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
+              break;
             }
             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 RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, 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> <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> <boundary> <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]), atoi(argv[argptr+5]));
+                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
                 if (Filling != NULL) {
                   Filling->ActiveFlag = false;

src/tesselation.cpp

@@ -35 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
   LinesCount(0),
   node(Walker),
@@ -61 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
 {
         Info FunctionInfo(__func__);
@@ -105 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const 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;
@@ -171 +193 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
 {
         Info FunctionInfo(__func__);
@@ -182 +204 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
 {
         Info FunctionInfo(__func__);
@@ -197 +219 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion()
+bool BoundaryLineSet::CheckConvexityCriterion() const
 {
         Info FunctionInfo(__func__);
@@ -221 +243 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::const_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);
@@ -264 +286 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
         Info FunctionInfo(__func__);
@@ -277 +299 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
 {
         Info FunctionInfo(__func__);
@@ -317 +339 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
   Nr(number)
 {
@@ -376 +398 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
 {
         Info FunctionInfo(__func__);
@@ -390 +412 @@
 /** 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
- * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * 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
@@ -400 +422 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
 {
   Info FunctionInfo(__func__);
@@ -452 +474 @@
 };
 
+/** 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(class BoundaryLineSet *line)
+bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
 {
         Info FunctionInfo(__func__);
@@ -469 +575 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
         Info FunctionInfo(__func__);
@@ -482 +588 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
 {
         Info FunctionInfo(__func__);
@@ -495 +601 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
 {
         Info FunctionInfo(__func__);
@@ -518 +624 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
 {
         Info FunctionInfo(__func__);
@@ -540 +646 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
 {
         Info FunctionInfo(__func__);
@@ -557 +663 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector *center)
+void BoundaryTriangleSet::GetCenter(Vector * const center) const
 {
         Info FunctionInfo(__func__);
@@ -564 +670 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
+  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -3135 +3242 @@
   // for each point, check its lines, remember closest
   Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
-  BoundaryLineSet *ClosestLine = NULL;
+  LineSet ClosestLines;
   double MinDistance = 1e+16;
   Vector BaseLineIntersection;
@@ -3157 +3264 @@
 
       if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
-        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 {
+        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
           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;
         }
@@ -3179 +3287 @@
           eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
         }
-        if ((ClosestLine == NULL) || (distance < MinDistance)) {
-          ClosestLine = LineRunner->second;
+        if ((ClosestLines.empty()) || (distance < MinDistance)) {
+          ClosestLines.insert(LineRunner->second);
           MinDistance = distance;
-          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " 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;
@@ -3192 +3300 @@
 
   // check whether closest line is "too close" :), then it's inside
-  if (ClosestLine == NULL) {
+  if (ClosestLines.empty()) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
   TriangleList * candidates = new TriangleList;
-  for (TriangleMap::iterator Runner = ClosestLine->triangles.begin(); Runner != ClosestLine->triangles.end(); Runner++) {
+  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
+    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
     candidates->push_back(Runner->second);
   }
@@ -3241 +3350 @@
 };
 
-/** Checks whether the provided Vector is within the tesselation structure.
+/** 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.
  * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
- * towards or away from the given \a &Point.
+ * 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!
  *
  * @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 >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
 {
   Info FunctionInfo(__func__);
@@ -3258 +3388 @@
   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 false;
+    return LC->RADIUS;
   } else {
     Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
@@ -3282 +3413 @@
     if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
       Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
-      return true;
+      return 0.;
     } 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 {
-    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
-    return false;
-  }
-}
+    // 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;
+    }
+  }
+};
 
 /** Gets all points connected to the provided point by triangulation lines.

src/tesselation.hpp

@@ -106 +106 @@
   public:
     BoundaryPointSet();
-    BoundaryPointSet(TesselPoint * Walker);
+    BoundaryPointSet(TesselPoint * const Walker);
     ~BoundaryPointSet();
 
-    void AddLine(class BoundaryLineSet *line);
+    void AddLine(BoundaryLineSet * const line);
 
     LineMap lines;
@@ -125 +125 @@
   public:
     BoundaryLineSet();
-    BoundaryLineSet(class BoundaryPointSet *Point[2], const int number);
+    BoundaryLineSet(BoundaryPointSet * const Point[2], const int number);
+    BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number);
     ~BoundaryLineSet();
 
-    void AddTriangle(class BoundaryTriangleSet *triangle);
-    bool IsConnectedTo(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-    bool CheckConvexityCriterion();
-    class BoundaryPointSet *GetOtherEndpoint(class BoundaryPointSet *);
+    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;
 
     class BoundaryPointSet *endpoints[2];
@@ -147 +148 @@
   public:
     BoundaryTriangleSet();
-    BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
+    BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number);
     ~BoundaryTriangleSet();
 
-    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);
+    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;
 
     class BoundaryPointSet *endpoints[3];
@@ -313 +315 @@
     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 \
-  TesselationInOutsideUnitTest \
+  Tesselation_BoundaryTriangleUnitTest \
+  Tesselation_InOutsideUnitTest \
   VectorUnitTest
   
@@ -79 +80 @@
 TesselationUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
-TesselationInOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
-TesselationInOutsideUnitTest_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
 
 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
+ * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-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 (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
-    eLog() << Verbose(2) << "Line is parallel to plane, no intersection." << endl;
+  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
+    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
     return false;
   }
@@ -243 +244 @@
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
-  if (factor < MYEPSILON) { // Origin is in-plane
-    Log() << Verbose(1) << "Origin of line is in-plane, simple." << endl;
+  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
+    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
     CopyVector(Origin);
     return true;
@@ -258 +259 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    Log() << Verbose(1) << "INFO: Intersection at " << *this << " is good." << endl;
+    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
     return true;
   } else {
@@ -353 +354 @@
   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);