Changeset 9473f6

Timestamp:

Jan 14, 2010, 5:36:09 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

e9c677

Parents:

7b36fe

Message:

FillBoxWithMolecule now allows boundary between filler molecules and the surface.

• boundary.cpp:
  • FillBoxWithMolecule() has additional parameter boundary and calls Tesselation::GetDistanceSquaredToSurface() to obtain distance to surface.
  • FillBoxWithMolecule(*) verbosity states now inner point or within boundary.

• builder.cpp:
  • ParseCommandLineOptions() in case 'F' takes seven arguments instead of six and calls FillBoxWithMolecule() correctly.

• tesselation.cpp:
  • Make const what needs to be const: BoundaryPointSet, BoundaryLineSet, BoundaryTriangleSet
  • added constructor BoundaryLineSet::BoundaryLineSet() that accept two endpoints instead of an array of two
  • new function BoundaryTriangleSet::GetClosestPointInsideTriangle() calculates closest point in a given triangle to a reference point and returns squared distance.
  • new function Tesselation::GetDistanceSquaredToSurface() that calculates the squared distance whose sign indicates whether reference is inside (-) or outside (+) and uses BoundaryTriangleSet::GetClosestPointInsideTriangle()
  • function Tesselation::IsInnerPoint() renamed to Tesselation::GetDistanceSquaredToSurface().
  • new function Tesselation::IsInnerPoint() just uses Tesselation::GetDistanceSquaredToSurface()

Location:

src

Files:

• boundary.cpp (modified) (3 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (2 diffs)
• tesselation.cpp (modified) (25 diffs)
• tesselation.hpp (modified) (4 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;
 }
 
@@ -3243 +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__);
@@ -3260 +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;
@@ -3284 +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;