Changeset 71b20e

Timestamp:

Dec 19, 2009, 7:32:24 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:

c15ca2

Parents:

3930eb

Message:

Attempt to fix the embedding.

Basically it would be working, but there was some failures with the FindClosestTriangleToPoint() routines.
We get triangles wrong if we start looking for the closest point. Actually, we should really look at each
triangle and check the distance. Now, we look at least at each line, but code is unfinished and crashes at the end
unexplainedly.

Location:

src

Files:

• boundary.cpp (modified) (7 diffs)
• boundary.hpp (modified) (2 diffs)
• molecule.hpp (modified) (1 diff)
• molecule_pointcloud.cpp (modified) (1 diff)
• tesselation.cpp (modified) (11 diffs)
• tesselation.hpp (modified) (2 diffs)
• tesselationhelpers.cpp (modified) (5 diffs)
• tesselationhelpers.hpp (modified) (1 diff)

src/boundary.cpp

@@ -654 +654 @@
  * \param *out output stream for debugging
  * \param *mol molecule with atoms and bonds
- * \param *&TesselStruct Tesselation with boundary triangles
+ * \param *TesselStruct Tesselation with boundary triangles
  * \param *filename prefix of filename
  * \param *extraSuffix intermediate suffix
  */
-void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
+void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
 {
         Info FunctionInfo(__func__);
@@ -789 +789 @@
  * \param configuration contains box dimensions
  * \param distance[NDIM] distance between filling molecules in each direction
+ * \param epsilon distance to surface which is not filled
  * \param RandAtomDisplacement maximum distance for random displacement per atom
  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
@@ -794 +795 @@
  * \return *mol pointer to new molecule with filled atoms
  */
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double epsilon, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
 {
         Info FunctionInfo(__func__);
@@ -817 +818 @@
   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
-    LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
-    if (TesselStruct[i] == NULL) {
-      Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
-      FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
-    }
+    LCList[i] = new LinkedCell((*ListRunner), 10.); // get linked cell list
+    Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
+    TesselStruct[i] = NULL;
+    FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
     i++;
   }
@@ -835 +835 @@
   FillerDistance.Init(distance[0], distance[1], distance[2]);
   FillerDistance.InverseMatrixMultiplication(M);
-  Log() << Verbose(1) << "INFO: Grid steps are ";
-  for(int i=0;i<NDIM;i++) {
+  for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance.x[i]);
-    Log() << Verbose(1) << N[i];
-    if (i != NDIM-1)
-      Log() << Verbose(1)<< ", ";
-    else
-      Log() << Verbose(1) << "." << endl;
-  }
+  Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;
 
   // go over [0,1]^3 filler grid
@@ -862 +856 @@
             FillIt = false;
           } else {
-            FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
+            const bool FillResult = (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i], epsilon));
+            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." << endl;
+            }
             i++;
           }
         }
 
-        if (FillIt) {
+        if (!FillIt) {
           // fill in Filler
           Log() << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
@@ -931 +931 @@
       }
   Free(&M);
+
+  // output to file
+  TesselStruct[0]->LastTriangle = NULL;
+  StoreTrianglesinFile(Filling, TesselStruct[0], "Tesselated", ".dat");
+
   for (size_t i=0;i<List->ListOfMolecules.size();i++) {
         delete(LCList[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, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double epsilon, 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);
@@ -58 +58 @@
 void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity);
 bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
-void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix);
+void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix);
 double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration);
 

src/molecule.hpp

@@ -110 +110 @@
   TesselPoint *GetPoint() const ;
   TesselPoint *GetTerminalPoint() const ;
+  int GetMaxId() const;
   void GoToNext() const ;
   void GoToPrevious() const ;

src/molecule_pointcloud.cpp

@@ -50 +50 @@
 };
 
+/** Return the greatest index of all atoms in the list.
+ * \return greatest index
+ */
+int molecule::GetMaxId() const
+{
+  return last_atom;
+};
+
 /** Go to next atom.
  * Stops at last one.

src/tesselation.cpp

@@ -482 +482 @@
 {
         Info FunctionInfo(__func__);
+        Log() << Verbose(1) << "INFO: Checking " << Points[0] << ","  << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl;
   return (((endpoints[0] == Points[0])
             || (endpoints[0] == Points[1])
@@ -2986 +2987 @@
 {
         Info FunctionInfo(__func__);
-  TesselPoint *trianglePoints[3];
-  TesselPoint *SecondPoint = NULL;
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  PointMap::const_iterator FindPoint;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
 
   if (LinesOnBoundary.empty()) {
-    eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first.";
+    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
     return NULL;
   }
-  Log() << Verbose(1) << "Finding closest Tesselpoint to " << *x << " ... " << endl;
-  trianglePoints[0] = FindClosestPoint(x, SecondPoint, LC);
-  
-  // check whether closest point is "too close" :), then it's inside
-  if (trianglePoints[0] == NULL) {
+
+  // gather all points close to the desired one
+  LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  set<BoundaryPointSet *> points;
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //Log() << Verbose(1) << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
+        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            FindPoint = PointsOnBoundary.find((*Runner)->nr);
+            if (FindPoint != PointsOnBoundary.end())
+              points.insert(FindPoint->second);
+          }
+        } else {
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  // check whether we found some points
+  if (points.empty()) {
+    Log() << Verbose(1) << "There is no nearest triangle, too far away from the surface." << endl;
+    return NULL;
+  }
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl;
+  BoundaryLineSet *ClosestLine = NULL;
+  double MinDistance = -1.;
+  Vector helper;
+  Vector Center;
+  Vector BaseLine;
+  for (set<BoundaryPointSet *>::iterator Runner = points.begin(); Runner != points.end(); Runner++) {
+    for (LineMap::iterator LineRunner = (*Runner)->lines.begin(); LineRunner != (*Runner)->lines.end(); LineRunner++) {
+      // calculate closest point on line to desired point
+      helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
+      helper.Scale(0.5);
+      Center.CopyVector(x);
+      Center.SubtractVector(&helper);
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      Center.ProjectOntoPlane(&BaseLine);
+      const double distance = Center.NormSquared();
+      if ((ClosestLine == NULL) || (distance < MinDistance)) {
+        // additionally calculate intersection on line (whether it's on the line section or not)
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthA = helper.ScalarProduct(&BaseLine);
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthB = helper.ScalarProduct(&BaseLine);
+        if (lengthB*lengthA < 0) {  // if have different sign
+          ClosestLine = LineRunner->second;
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
+        }
+      } else {
+        Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+      }
+    }
+  }
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLine == NULL) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
-  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    Log() << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
-    PointMap::const_iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
-    triangles = new list<BoundaryTriangleSet*>;
-    if (PointRunner != PointsOnBoundary.end()) {
-      for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
-        for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
-          triangles->push_back(TriangleRunner->second);
-      triangles->sort();
-      triangles->unique();
-    } else {
-      PointRunner = PointsOnBoundary.find(SecondPoint->nr);
-      trianglePoints[0] = SecondPoint;
-      if (PointRunner != PointsOnBoundary.end()) {
-        for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
-          for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
-            triangles->push_back(TriangleRunner->second);
-        triangles->sort();
-        triangles->unique();
-      } else {
-        eLog() << Verbose(1) << "I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
-        return NULL;
-      }
-    }
-  } else {
-    set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(trianglePoints[0]);
-    TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(connectedPoints, trianglePoints[0], x);
-    delete(connectedPoints);
-    if (connectedClosestPoints != NULL) {
-      trianglePoints[1] = connectedClosestPoints->front();
-      trianglePoints[2] = connectedClosestPoints->back();
-      for (int i=0;i<3;i++) {
-        if (trianglePoints[i] == NULL) {
-          eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-        }
-        //Log() << Verbose(1) << "List of triangle points:" << endl;
-        //Log() << Verbose(2) << *trianglePoints[i] << endl;
-      }
-
-      triangles = FindTriangles(trianglePoints);
-      Log() << Verbose(1) << "List of possible triangles:" << endl;
-      for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-        Log() << Verbose(2) << **Runner << endl;
-
-      delete(connectedClosestPoints);
-    } else {
-      triangles = NULL;
-      eLog() << Verbose(2) << "There is no circle of connected points!" << endl;
-    }
-  }
-  
-  if ((triangles == NULL) || (triangles->empty())) {
-    eLog() << Verbose(1) << "There is no nearest triangle. Please check the tesselation structure.";
-    delete(triangles);
-    return NULL;
-  } else
-    return triangles;
+  list<BoundaryTriangleSet*> *triangles = new list<BoundaryTriangleSet*>;
+  for (TriangleMap::iterator Runner = ClosestLine->triangles.begin(); Runner != ClosestLine->triangles.end(); Runner++) {
+    triangles->push_back(Runner->second);
+  }
+  return triangles;
 };
 
@@ -3072 +3090 @@
   class BoundaryTriangleSet *result = NULL;
   list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
+  list<BoundaryTriangleSet*> candidates;
   Vector Center;
-
-  if (triangles == NULL)
+  Vector helper;
+
+  if ((triangles == NULL) || (triangles->empty()))
     return NULL;
 
-  if (triangles->size() == 1) { // there is no degenerate case
-    result = triangles->front();
-    Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
+  // reject all which are not closest
+  double MinDistance = -1.;
+  for (list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+    (*Runner)->GetCenter(&Center);
+    helper.CopyVector(x);
+    helper.SubtractVector(&Center);
+    helper.ProjectOntoPlane(&(*Runner)->NormalVector);
+    const double distance = helper.NormSquared();
+    if (fabs(distance - MinDistance) < MYEPSILON) {  // degenerate case, keep both
+      candidates.push_back(*Runner);
+    } else if (distance < MinDistance) {
+      candidates.clear();
+      candidates.push_back(*Runner);
+      MinDistance = distance;
+    }
+  }
+  delete(triangles);
+  if (!candidates.empty()) {
+    if (candidates.size() == 1) { // there is no degenerate case
+      result = candidates.front();
+      Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
+    } else {
+      result = candidates.front();
+      result->GetCenter(&Center);
+      Center.SubtractVector(x);
+      Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
+      if (Center.ScalarProduct(&result->NormalVector) < 0) {
+        result = candidates.back();
+        Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
+        if (Center.ScalarProduct(&result->NormalVector) < 0) {
+          eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
+        }
+      }
+    }
   } else {
-    result = triangles->front();
-    result->GetCenter(&Center);
-    Center.SubtractVector(x);
-    Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
-    if (Center.ScalarProduct(&result->NormalVector) < 0) {
-      result = triangles->back();
-      Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
-      if (Center.ScalarProduct(&result->NormalVector) < 0) {
-        eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
-      }
-    }
-  }
-  delete(triangles);
+    result = NULL;
+    Log() << Verbose(0) << "No closest triangle found." << endl;
+  }
   return result;
 };
@@ -3101 +3142 @@
  * @param point of which to check the position
  * @param *LC LinkedCell structure
+ * @param epsilon Distance of \a &Point to Center of triangle (of point outwards) is tested less against this value (standard: -MYEPSILON)
  *
  * @return true if the point is inside the tesselation structure, false otherwise
  */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon) const
 {
         Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
   Vector Center;
+  Vector helper;
 
   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;
+  } else {
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << "." << endl;
   }
 
@@ -3118 +3163 @@
   Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
   Center.SubtractVector(&Point);
-  Log() << Verbose(2) << "INFO: Vector from center to point to test is " << Center << "." << endl;
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << Center << "." << endl;
   if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
     Log() << Verbose(1) << Point << " is an inner point." << endl;
     return true;
   } else {
+    for (int i=0;i<3;i++) {
+      helper.CopyVector(result->endpoints[i]->node->node);
+      helper.SubtractVector(&Point);
+      if (Center.NormSquared() > helper.NormSquared())
+        Center.CopyVector(&helper);
+    }
+    if (Center.NormSquared() < epsilon*epsilon) {
+      Log() << Verbose(1) << Point << " is inner point, being sufficiently close (wrt " << epsilon << ") to boundary." << endl;
+      return true;
+    }
     Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
     return false;
@@ -3132 +3187 @@
  * @param *Point of which to check the position
  * @param *LC Linked Cell structure
+ * @param epsilon Distance of \a &Point to Center of triangle is tested greater against this value (standard: -MYEPSILON)
  *
  * @return true if the point is inside the tesselation structure, false otherwise
  */
-bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const
-{
-        Info FunctionInfo(__func__);
-  return IsInnerPoint(*(Point->node), LC);
+bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon) const
+{
+        Info FunctionInfo(__func__);
+  return IsInnerPoint(*(Point->node), LC, epsilon);
 }
+
+/**
+ * Finds the point which is second closest to the provided one.
+ *
+ * @param Point to which to find the second closest other point
+ * @param linked cell structure
+ *
+ * @return point which is second closest to the provided one
+ * @TODO Optimize me! These finds are unnecessary if we have a LC list of BoundaryPoints only!
+ */
+TesselPoint* Tesselation::FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  TesselPoint* closestPoint = NULL;
+  TesselPoint* secondClosestPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //Log() << Verbose(1) << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
+        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (PointsOnBoundary.find((*Runner)->nr) != PointsOnBoundary.end()) {
+              if (currentNorm < distance) {
+                // remember second point
+                secondDistance = distance;
+                secondClosestPoint = closestPoint;
+                // mark down new closest point
+                distance = currentNorm;
+                closestPoint = (*Runner);
+                //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
+              }
+            }
+          }
+        } else {
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  return secondClosestPoint;
+};
+
+/**
+ * Finds the point which is closest to the provided one.
+ *
+ * @param Point to which to find the closest other point
+ * @param SecondPoint the second closest other point on return, NULL if none found
+ * @param linked cell structure
+ *
+ * @return point which is closest to the provided one, NULL if none found
+ * @TODO Optimize me! These finds are unnecessary if we have a LC list of BoundaryPoints only!
+ */
+TesselPoint* Tesselation::FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  TesselPoint* closestPoint = NULL;
+  SecondPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //Log() << Verbose(1) << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
+        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper.NormSquared();
+            if (PointsOnBoundary.find((*Runner)->nr) != PointsOnBoundary.end()) {
+              if (currentNorm < distance) {
+                secondDistance = distance;
+                SecondPoint = closestPoint;
+                distance = currentNorm;
+                closestPoint = (*Runner);
+                //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
+              } else if (currentNorm < secondDistance) {
+                secondDistance = currentNorm;
+                SecondPoint = (*Runner);
+                //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
+              }
+            }
+          }
+        } else {
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  // output
+  if ((closestPoint != NULL) && (SecondPoint != NULL)) {
+    Log() << Verbose(1) << "Closest point is " << *closestPoint << " and second closest is " << *SecondPoint << "." << endl;
+  } else if (closestPoint != NULL) {
+    Log() << Verbose(1) << "Closest point is " << *closestPoint << " and second closest is " << SecondPoint << "." << endl;
+  } else {
+    Log() << Verbose(1) << "No closest point has been found, only " << closestPoint << " and " << SecondPoint << "." << endl;
+  }
+  return closestPoint;
+};
+
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3191 +3371 @@
   }
 
-  if (connectedPoints->size() == 0) { // if have not found any points
+  if (connectedPoints->empty()) { // if have not found any points
     eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
     return NULL;
@@ -3212 +3392 @@
  * @return list of the all points linked to the provided one
  */
+list<TesselPoint*> * Tesselation::GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+{
+        Info FunctionInfo(__func__);
+  map<double, TesselPoint*> anglesOfPoints;
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  list<BoundaryTriangleSet*> *triangles = NULL;
+
+  if (SetOfNeighbours == NULL) {
+    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    delete(connectedCircle);
+    return NULL;
+  }
+
+  // calculate central point
+  triangles = FindTriangles(TrianglePoints);
+  if ((triangles != NULL) && (!triangles->empty())) {
+    for (list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      PlaneNormal.AddVector(&(*Runner)->NormalVector);
+  } else {
+    eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
+    performCriticalExit();
+  }
+  PlaneNormal.Scale(1.0/triangles->size());
+  Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl;
+  PlaneNormal.Normalize();
+
+  // construct one orthogonal vector
+  if (Reference != NULL) {
+    AngleZero.CopyVector(Reference);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+  }
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
+    Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
+    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
+      performCriticalExit();
+    }
+  }
+  Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
+  if (AngleZero.NormSquared() > MYEPSILON)
+    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+  else
+    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+  Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
+
+  // go through all connected points and calculate angle
+  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+    helper.CopyVector((*listRunner)->node);
+    helper.SubtractVector(Point->node);
+    helper.ProjectOntoPlane(&PlaneNormal);
+    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
+    Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
+    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+  }
+
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+    connectedCircle->push_back(AngleRunner->second);
+  }
+
+  return connectedCircle;
+}
+
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
+ * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
+ * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
+ * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
+ * triangle we are looking for.
+ *
+ * @param *out output stream for debugging
+ * @param *SetOfNeighbours all points for which the angle should be calculated
+ * @param *Point of which get all connected points
+ * @param *Reference Reference vector for zero angle or NULL for no preference
+ * @return list of the all points linked to the provided one
+ */
 list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
   list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
@@ -3715 +3978 @@
  * Finds triangles belonging to the three provided points.
  *
- * @param *Points[3] list, is expected to contain three points
+ * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
  *
  * @return triangles which belong to the provided points, will be empty if there are none,
@@ -3727 +3990 @@
   TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
+  size_t NoOfWildcards = 0;
 
   for (int i = 0; i < 3; i++) {
-    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      TrianglePoints[i] = FindPoint->second;
+    if (Points[i] == NULL) {
+      NoOfWildcards++;
+      TrianglePoints[i] = NULL;
     } else {
-      TrianglePoints[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (TrianglePoints[i] != NULL) {
-      for (int j = i+1; j < 3; j++) {
-        if (TrianglePoints[j] != NULL) {
-          for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
-              (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
-              FindLine++) {
-            for (FindTriangle = FindLine->second->triangles.begin();
-                FindTriangle != FindLine->second->triangles.end();
-                FindTriangle++) {
-              if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
-                result->push_back(FindTriangle->second);
+      PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
+      if (FindPoint != PointsOnBoundary.end()) {
+        TrianglePoints[i] = FindPoint->second;
+      } else {
+        TrianglePoints[i] = NULL;
+      }
+    }
+  }
+
+  switch (NoOfWildcards) {
+    case 0: // checks lines between the points in the Points for their adjacent triangles
+      for (int i = 0; i < 3; i++) {
+        if (TrianglePoints[i] != NULL) {
+          for (int j = i+1; j < 3; j++) {
+            if (TrianglePoints[j] != NULL) {
+              for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
+                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
+                  FindLine++) {
+                for (FindTriangle = FindLine->second->triangles.begin();
+                    FindTriangle != FindLine->second->triangles.end();
+                    FindTriangle++) {
+                  if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+                    result->push_back(FindTriangle->second);
+                  }
+                }
               }
+              // Is it sufficient to consider one of the triangle lines for this.
+              return result;
             }
           }
-          // Is it sufficient to consider one of the triangle lines for this.
-          return result;
         }
       }
-    }
+      break;
+    case 1: // copy all triangles of the respective line
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] == NULL)
+          break;
+      for (FindLine = TrianglePoints[(i+1)%3]->lines.find(TrianglePoints[(i+2)%3]->node->nr); // is a multimap!
+          (FindLine != TrianglePoints[(i+1)%3]->lines.end()) && (FindLine->first == TrianglePoints[(i+2)%3]->node->nr);
+          FindLine++) {
+        for (FindTriangle = FindLine->second->triangles.begin();
+            FindTriangle != FindLine->second->triangles.end();
+            FindTriangle++) {
+          if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+            result->push_back(FindTriangle->second);
+          }
+        }
+      }
+      break;
+    }
+    case 2: // copy all triangles of the respective point
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] != NULL)
+          break;
+      for (LineMap::const_iterator line = TrianglePoints[i]->lines.begin(); line != TrianglePoints[i]->lines.end(); line++)
+        for (TriangleMap::const_iterator triangle = line->second->triangles.begin(); triangle != line->second->triangles.end(); triangle++)
+          result->push_back(triangle->second);
+      result->sort();
+      result->unique();
+      break;
+    }
+    case 3: // copy all triangles
+    {
+      for (TriangleMap::const_iterator triangle = TrianglesOnBoundary.begin(); triangle != TrianglesOnBoundary.end(); triangle++)
+        result->push_back(triangle->second);
+      break;
+    }
+    default:
+      eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl;
+      performCriticalExit();
+      break;
   }
 
@@ -3982 +4296 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;

src/tesselation.hpp

@@ -226 +226 @@
   virtual TesselPoint *GetPoint() const { return NULL; };
   virtual TesselPoint *GetTerminalPoint() const { return NULL; };
+  virtual int GetMaxId() const { return 0; };
   virtual void GoToNext() const {};
   virtual void GoToPrevious() const {};
@@ -301 +302 @@
     list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
     list<TesselPoint*> * GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
+    list<TesselPoint*> * GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
     class BoundaryPointSet *GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
     list<BoundaryTriangleSet*> *FindTriangles(const TesselPoint* const Points[3]) const;
     list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const;
     class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
-    bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
-    bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const;
+    bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
+    bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
+    TesselPoint* FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const;
+    TesselPoint* FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) const;
 
     // print for debugging

src/tesselationhelpers.cpp

@@ -558 +558 @@
  * @return point which is second closest to the provided one
  */
-TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
+TesselPoint* FindSecondClosestTesselPoint(const Vector* Point, const LinkedCell* const LC)
 {
         Info FunctionInfo(__func__);
@@ -613 +613 @@
  * @return point which is closest to the provided one, NULL if none found
  */
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
+TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
 {
         Info FunctionInfo(__func__);
@@ -639 +639 @@
             helper.CopyVector(Point);
             helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
+            double currentNorm = helper.NormSquared();
             if (currentNorm < distance) {
               secondDistance = distance;
@@ -809 +809 @@
   TesselPoint *Walker = NULL;
   int i;
-  Vector *center = cloud->GetCenter();
+  Vector *center = new Vector(); //cloud->GetCenter();
   if (rasterfile != NULL) {
     //Log() << Verbose(1) << "Writing Raster3D file ... ";
@@ -866 +866 @@
     }
     *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int i=0;
-    for (cloud->GoToFirst(); !cloud->IsEnd(); cloud->GoToNext(), i++);
+    int i=cloud->GetMaxId();
     int *LookupList = new int[i];
     for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)

src/tesselationhelpers.hpp

@@ -59 +59 @@
 bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3]);
 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation *candidate2);
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC);
-TesselPoint* FindSecondClosestPoint(const Vector*, const LinkedCell* const LC);
+TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC);
+TesselPoint* FindSecondClosestTesselPoint(const Vector*, const LinkedCell* const LC);
 Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase);