Diff [bd61b4123d6339083d04e167580d407cbf99c263:271e17edd356e78553feaa53e435fda094d20a46] for / – MoleCuilder

src/Makefile.am

-              rbd61b41
+              r271e17
 ATOMSOURCE = atom.cpp atom_atominfo.cpp atom_bondedparticle.cpp atom_bondedparticleinfo.cpp atom_graphnode.cpp atom_graphnodeinfo.cpp atom_particleinfo.cpp atom_trajectoryparticle.cpp atom_trajectoryparticleinfo.cpp
 ATOMHEADER = atom.hpp atom_atominfo.hpp atom_bondedparticle.hpp atom_bondedparticleinfo.hpp atom_graphnode.hpp atom_graphnodeinfo.hpp atom_particleinfo.hpp atom_trajectoryparticle.hpp atom_trajectoryparticleinfo.hpp
+LINALGSOURCE = gslmatrix.cpp gslvector.cpp linearsystemofequations.cpp
+LINALGHEADER = gslmatrix.hpp gslvector.hpp linearsystemofequations.hpp
 ANALYSISSOURCE = analysis_bonds.cpp analysis_correlation.cpp
 …
 INCLUDES = -I$(top_srcdir)/src/unittests
 noinst_LIBRARIES = libmolecuilder.a
+noinst_LIBRARIES = libmolecuilder.a libgslwrapper.a
 bin_PROGRAMS = molecuilder joiner analyzer
 molecuilderdir = ${bindir}
 libmolecuilder_a_SOURCES = ${SOURCE} ${HEADER}
+libgslwrapper_a_SOURCES = ${LINALGSOURCE} ${LINALGHEADER}
 molecuilder_DATA = elements.db valence.db orbitals.db Hbonddistance.db Hbondangle.db
 molecuilder_LDFLAGS = $(BOOST_LIB)
 molecuilder_SOURCES = builder.cpp
 molecuilder_LDADD = libmolecuilder.a
+molecuilder_LDADD = libmolecuilder.a libgslwrapper.a
 joiner_SOURCES = joiner.cpp datacreator.cpp parser.cpp datacreator.hpp helpers.hpp parser.hpp periodentafel.hpp
 joiner_LDADD = libmolecuilder.a

src/analysis_correlation.cpp

-              rbd61b41
+              r271e17
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
           triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
+          triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
           if (triangle != NULL) {
             distance = DistanceToTrianglePlane(Walker->node, triangle);
 …
+  }
   outmap = new CorrelationToSurfaceMap;
+  double ShortestDistance = 0.;
+  BoundaryTriangleSet *ShortestTriangle = NULL;
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
 …
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
+          ShortestDistance = -1.;
           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
 …
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
+                triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
+                if (triangle != NULL) {
+                  distance = DistanceToTrianglePlane(&checkX, triangle);
+                  outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
+                triangle = Surface->FindClosestTriangleToVector(&checkX, LC);
+                distance = Surface->GetDistanceSquaredToTriangle(checkX, triangle);
+                if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
+                  ShortestDistance = distance;
+                  ShortestTriangle = triangle;
+                }
+          }
+              }
+          // insert
+          ShortestDistance = sqrt(ShortestDistance);
+          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
+          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
+        }
+      }
 …
+{
   Info FunctionInfo(__func__);
   *file << "# BinStart\tCount" << endl;
+  *file << "BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << runner->second << endl;
 …
+{
   Info FunctionInfo(__func__);
   *file << "# BinStart\tAtom1\tAtom2" << endl;
+  *file << "BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
 …
+{
   Info FunctionInfo(__func__);
   *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
+  *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
   for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first;
 …
+{
   Info FunctionInfo(__func__);
   *file << "# BinStart\tTriangle" << endl;
+  *file << "BinStart\tTriangle" << endl;
   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << *(runner->second.second) << endl;
+  }
 };
+    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+  }
+};

src/analysis_correlation.hpp

-              rbd61b41
+              r271e17
 int GetBin ( const double value, const double BinWidth, const double BinStart );
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map );
 void OutputPairCorrelation( ofstream * const file, const BinPairMap * const map );
 void OutputCorrelationToPoint( ofstream * const file, const BinPairMap * const map );
 void OutputCorrelationToSurface( ofstream * const file, const BinPairMap * const map );
+void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map );
+void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map );
+void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map );
 …
 /** Puts given correlation data into bins of given size (histogramming).
  * Note that BinStart and BinEnd are desired to fill the complete range, even where Bins are zero. If this is
+ * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater.
+ * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater. If a
+ * certain start value is desired, give BinStart and a BinEnd that is smaller than BinStart, then only BinEnd will be
+ * calculated automatically, i.e. BinStart = 0. and BinEnd = -1. .
  * Also note that the range is given inclusive, i.e. [ BinStart, BinEnd ].
  * \param *map map of doubles to count
 …
   if (BinStart == BinEnd) { // if same, find range ourselves
     GetMinMax( map, start, end);
+  } else { // if not, initialise range to zero
+  } else if (BinEnd < BinStart) { // if BinEnd smaller, just look for new max
+    GetMinMax( map, start, end);
+    start = BinStart;
+  } else { // else take both values
     start = BinStart;
     end = BinEnd;

src/analyzer.cpp

rbd61b41	r271e17
7	7
8	8	//============================ INCLUDES ===========================
	9
	10	#include <cstring>
9	11
10	12	#include "datacreator.hpp"

src/bondgraph.cpp

-              rbd61b41
+              r271e17
 #include "bondgraph.hpp"
 #include "element.hpp"
+#include "info.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
 …
 bool BondGraph::LoadBondLengthTable(const string &filename)
+{
+  Info FunctionInfo(__func__);
   bool status = true;
   MatrixContainer *TempContainer = NULL;
 …
   // parse in matrix
+  status = TempContainer->ParseMatrix(filename.c_str(), 0, 1, 0);
+  if (status = TempContainer->ParseMatrix(filename.c_str(), 0, 1, 0)) {
+    Log() << Verbose(1) << "Parsing bond length matrix successful." << endl;
+  } else {
+    eLog() << Verbose(1) << "Parsing bond length matrix failed." << endl;
+  }
   // find greatest distance

src/boundary.cpp

-              rbd61b41
+              r271e17
 #include<gsl/gsl_poly.h>
+#include<time.h>
 // ========================================== F U N C T I O N S =================================
 …
  * \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__);
 …
  * \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
  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
 …
  * \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 boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+{
         Info FunctionInfo(__func__);
 …
   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++;
+  }
 …
   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;
+  // initialize seed of random number generator to current time
+  srand ( time(NULL) );
   // go over [0,1]^3 filler grid
 …
           // get linked cell list
           if (TesselStruct[i] == NULL) {
             eLog() << Verbose(1) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
+            eLog() << Verbose(0) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
             FillIt = false;
           } else {
+            FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
+            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 or within boundary." << endl;
+              break;
+            }
             i++;
+          }
 …
+      }
   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

-              rbd61b41
+              r271e17
 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 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);
 …
 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/builder.cpp

-              rbd61b41
+              r271e17
 using namespace std;
+#include <cstring>
 #include "analysis_correlation.hpp"
 …
             Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
             Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
             Log() << Verbose(0) << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
             Log() << Verbose(0) << "\t-F\tFilling Box with water molecules." << endl;
+            Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
+            Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
             Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
             Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
 …
      if (configuration.BG == NULL) {
        configuration.BG = new BondGraph(configuration.GetIsAngstroem());
        if ((BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) {
+       if ((!BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) {
          Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
        } else {
 …
               Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;
               // @TODO rather do the dissection afterwards
               molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration);
+              molecules->DissectMoleculeIntoConnectedSubgraphs(periode, &configuration);
               mol = NULL;
               if (molecules->ListOfMolecules.size() != 0) {
 …
                 int ranges[NDIM] = {1,1,1};
                 CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
+                OutputCorrelationToSurface(&output, surfacemap);
                 BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
                 OutputCorrelation ( &binoutput, binmap );
 …
               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> <epsilon> <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 {
 …
                 Log() << Verbose(1) << "Filling Box with water molecules." << endl;
                 // construct water molecule
                 molecule *filler = new molecule(periode);;
+                molecule *filler = new molecule(periode);
                 molecule *Filling = NULL;
                 atom *second = NULL, *third = NULL;
+//                first = new atom();
+//                first->type = periode->FindElement(5);
+//                first->x.Zero();
+//                filler->AddAtom(first);
                 first = new atom();
                 first->type = periode->FindElement(5);
                 first->x.Zero();
+                first->type = periode->FindElement(1);
+                first->x.Init(0.441, -0.143, 0.);
                 filler->AddAtom(first);
+//                first = new atom();
+//                first->type = periode->FindElement(1);
+//                first->x.Init(0.441, -0.143, 0.);
+//                filler->AddAtom(first);
+//                second = new atom();
+//                second->type = periode->FindElement(1);
+//                second->x.Init(-0.464, 1.137, 0.0);
+//                filler->AddAtom(second);
+//                third = new atom();
+//                third->type = periode->FindElement(8);
+//                third->x.Init(-0.464, 0.177, 0.);
+//                filler->AddAtom(third);
+//                filler->AddBond(first, third, 1);
+//                filler->AddBond(second, third, 1);
+                second = new atom();
+                second->type = periode->FindElement(1);
+                second->x.Init(-0.464, 1.137, 0.0);
+                filler->AddAtom(second);
+                third = new atom();
+                third->type = periode->FindElement(8);
+                third->x.Init(-0.464, 0.177, 0.);
+                filler->AddAtom(third);
+                filler->AddBond(first, third, 1);
+                filler->AddBond(second, third, 1);
                 // call routine
                 double distance[NDIM];

src/config.cpp

-              rbd61b41
+              r271e17
 #include <stdio.h>
+#include <cstring>
 #include "atom.hpp"
 …
     char number1[8];
     char number2[8];
     char *dummy1, *dummy2;
+    const char *dummy1, *dummy2;
     //Log() << Verbose(0) << s1 << "  " << s2 << endl;
     dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
 …
+              }
               line++;
             } while (dummy1 != NULL && (dummy1[0] == '#') || (dummy1[0] == '\n'));
+            } while ((dummy1 != NULL) && ((dummy1[0] == '#') || (dummy1[0] == '\n')));
             dummy = dummy1;
           } else { // simple int, strings or doubles start in the same line

src/helpers.hpp

-              rbd61b41
+              r271e17
   x = y;
   y = tmp;
+};
+/** returns greater of the two values.
+ * \param x first value
+ * \param y second value
+ * \return greater of the two (by operator>())
+ */
+template <typename T> T Max(T x, T y)
+{
+  if (x > y)
+    return x;
+  else return y;
+};
+/** returns smaller of the two values.
+ * \param x first value
+ * \param y second value
+ * \return smaller of the two (by operator<())
+ */
+template <typename T> T Min(T x, T y)
+{
+  if (x < y)
+    return x;
+  else return y;
 };

src/joiner.cpp

rbd61b41	r271e17
7	7
8	8	//============================ INCLUDES ===========================
	9
	10	#include <cstring>
9	11
10	12	#include "datacreator.hpp"

src/memoryallocator.hpp

rbd61b41	r271e17
16	16	#endif
17	17
	18	#include <cstdlib>
18	19	#include <iostream>
19	20	#include <iomanip>

src/memoryusageobserver.cpp

rbd61b41	r271e17
4	4	* This class represents a Singleton for observing memory usage.
5	5	*/
	6
	7	#include <cstdlib>
6	8
7	9	#include "log.hpp"

src/molecule.cpp

-              rbd61b41
+              r271e17
+ *
  */
+#include <cstring>
 #include "atom.hpp"
 …
   else
     molname = filename; // contains no slashes
   char *endname = strchr(molname, '.');
+  const char *endname = strchr(molname, '.');
   if ((endname == NULL) || (endname < molname))
     length = strlen(molname);

src/molecule.hpp

-              rbd61b41
+              r271e17
   TesselPoint *GetPoint() const ;
   TesselPoint *GetTerminalPoint() const ;
+  int GetMaxId() const;
   void GoToNext() const ;
   void GoToPrevious() const ;
 …
   void Enumerate(ofstream *out);
   void Output(ofstream *out);
   void DissectMoleculeIntoConnectedSubgraphs(molecule * const mol, config * const configuration);
+  void DissectMoleculeIntoConnectedSubgraphs(const periodentafel * const periode, config * const configuration);
   int CountAllAtoms() const;

src/molecule_dynamics.cpp

-              rbd61b41
+              r271e17
   // argument minimise the constrained potential in this injective PermutationMap
   Log() << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
+  Log() << Verbose(1) << "Argument minimising the PermutationMap." << endl;
   OldPotential = 1e+10;
   round = 0;
   do {
     Log() << Verbose(2) << "Starting round " << ++round << " ... " << endl;
+    Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl;
     OlderPotential = OldPotential;
     do {

src/molecule_fragmentation.cpp

rbd61b41	r271e17
5	5	* Author: heber
6	6	*/
	7
	8	#include <cstring>
7	9
8	10	#include "atom.hpp"

src/molecule_pointcloud.cpp

-              rbd61b41
+              r271e17
 };
+/** 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/moleculelist.cpp

-              rbd61b41
+              r271e17
+ *
  */
+#include <cstring>
 #include "atom.hpp"
 …
 /** Dissects given \a *mol into connected subgraphs and inserts them as new molecules but with old atoms into \a this.
  * \param *out output stream for debugging
  * \param *mol molecule with atoms to dissect
+ * \param *periode periodentafel
  * \param *configuration config with BondGraph
  */
+void MoleculeListClass::DissectMoleculeIntoConnectedSubgraphs(molecule * const mol, config * const configuration)
+{
+void MoleculeListClass::DissectMoleculeIntoConnectedSubgraphs(const periodentafel * const periode, config * const configuration)
+{
+  molecule *mol = new molecule(periode);
+  atom *Walker = NULL;
+  atom *Advancer = NULL;
+  bond *Binder = NULL;
+  bond *Stepper = NULL;
+  // 0. gather all atoms into single molecule
+  for (MoleculeList::iterator MolRunner = ListOfMolecules.begin(); !ListOfMolecules.empty(); MolRunner = ListOfMolecules.begin()) {
+    // shift all atoms to new molecule
+    Advancer = (*MolRunner)->start->next;
+    while (Advancer != (*MolRunner)->end) {
+      Walker = Advancer;
+      Advancer = Advancer->next;
+      Log() << Verbose(3) << "Re-linking " << *Walker << "..." << endl;
+      unlink(Walker);
+      Walker->father = Walker;
+      mol->AddAtom(Walker);    // counting starts at 1
+    }
+    // remove all bonds
+    Stepper = (*MolRunner)->first->next;
+    while (Stepper != (*MolRunner)->last) {
+      Binder = Stepper;
+      Stepper = Stepper->next;
+      delete(Binder);
+    }
+    // remove the molecule
+    delete(*MolRunner);
+    ListOfMolecules.erase(MolRunner);
+  }
   // 1. dissect the molecule into connected subgraphs
   configuration->BG->ConstructBondGraph(mol);
 …
   int *MolMap = Calloc<int>(mol->AtomCount, "config::Load() - *MolMap");
   MoleculeLeafClass *MolecularWalker = Subgraphs;
   atom *Walker = NULL;
+  Walker = NULL;
   while (MolecularWalker->next != NULL) {
     MolecularWalker = MolecularWalker->next;
 …
+  }
   // 4d. we don't need to redo bonds, as they are connected subgraphs and still maintain their ListOfBonds, but we have to remove them from first..last list
   bond *Binder = mol->first;
+  Binder = mol->first;
   while (mol->first->next != mol->last) {
     Binder = mol->first->next;

src/parser.cpp

-              rbd61b41
+              r271e17
 // ======================================= INCLUDES ==========================================
+#include <cstring>
 #include "helpers.hpp"
 …
   input.open(name, ios::in);
   //Log() << Verbose(0) << "Opening " << name << " ... "  << input << endl;
+  //Log() << Verbose(1) << "Opening " << name << " ... "  << input << endl;
   if (input == NULL) {
     eLog() << Verbose(1) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
 …
+  }
   //Log() << Verbose(0) << line.str() << endl;
   //Log() << Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
+  //Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
   if (ColumnCounter[MatrixNr] == 0) {
     eLog() << Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
 …
+    }
+  }
   //Log() << Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
+  //Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
   if (RowCounter[MatrixNr] == 0) {
     eLog() << Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
 …
       input.getline(filename, 1023);
       stringstream lines(filename);
       //Log() << Verbose(0) << "Matrix at level " << j << ":";// << filename << endl;
+      //Log() << Verbose(2) << "Matrix at level " << j << ":";// << filename << endl;
       for(int k=skipcolumns;k--;)
         lines >> filename;
       for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
         lines >> Matrix[MatrixNr][j][k];
         //Log() << Verbose(0) << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
+        //Log() << Verbose(1) << " " << setprecision(2) << Matrix[MatrixNr][j][k] << endl;
+      }
-      //Log() << Verbose(0) << endl;
       Matrix[MatrixNr][ RowCounter[MatrixNr] ] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
       for(int j=ColumnCounter[MatrixNr];j--;)

src/periodentafel.cpp

rbd61b41	r271e17
9	9	#include <iomanip>
10	10	#include <fstream>
	11	#include <cstring>
11	12
12	13	#include "element.hpp"

src/tesselation.cpp

-              rbd61b41
+              r271e17
  * \param *Walker TesselPoint this boundary point represents
  */
 BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
   LinesCount(0),
   node(Walker),
 …
  * \param *line line to add
  */
 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+{
         Info FunctionInfo(__func__);
 …
  * \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__);
 …
   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;
 …
  * \param *triangle to add
  */
 void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+{
         Info FunctionInfo(__func__);
 …
  * \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__);
 …
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
 bool BoundaryLineSet::CheckConvexityCriterion()
+bool BoundaryLineSet::CheckConvexityCriterion() const
+{
         Info FunctionInfo(__func__);
 …
   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);
 …
  * \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__);
 …
  * \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__);
 …
  * \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)
+{
 …
  * \param &OtherVector direction vector to make normal vector unique.
  */
 void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+{
         Info FunctionInfo(__func__);
 …
 };
 /** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+/** 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
 …
  * \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)
+{
         Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+{
+  Info FunctionInfo(__func__);
   Vector CrossPoint;
   Vector helper;
 …
+  }
+  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
+  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
+  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl;
+    return true;
+  }
   // Calculate cross point between one baseline and the line from the third endpoint to intersection
   int i=0;
 …
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
+      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+      const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared();
+      Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl;
+      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
+        Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl;
+        i=4;
+        break;
+      }
+      i++;
     } else
-      i++;
-    if (i>2)
       break;
   } while (CrossPoint.NormSquared() < MYEPSILON);
+  } while (i<3);
   if (i==3) {
+    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
+  }
+  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+  // check whether intersection is inside or not by comparing length of intersection and length of cross point
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
   } else { // outside!
     Intersection->Zero();
+  } else {
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
     return false;
+  }
+};
+/** 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;
 };
 …
  * \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__);
 …
  * \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__);
 …
  * \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__);
 …
  * \return true - is the very triangle, false - is not
  */
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
+{
+        Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+{
+        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])
 …
  * \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__);
 …
  * \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__);
 …
  * \param *center central point on return.
  */
 void BoundaryTriangleSet::GetCenter(Vector *center)
+void BoundaryTriangleSet::GetCenter(Vector * const center) const
+{
         Info FunctionInfo(__func__);
 …
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
+  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
+}
 …
 TesselPoint::TesselPoint()
+{
   Info FunctionInfo(__func__);
+  //Info FunctionInfo(__func__);
   node = NULL;
   nr = -1;
 …
 TesselPoint::~TesselPoint()
+{
   Info FunctionInfo(__func__);
+  //Info FunctionInfo(__func__);
 };
 …
 PointCloud::PointCloud()
+{
         Info FunctionInfo(__func__);
+        //Info FunctionInfo(__func__);
 };
 …
 PointCloud::~PointCloud()
+{
         Info FunctionInfo(__func__);
+        //Info FunctionInfo(__func__);
 };
 …
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
  */
+void
+Tesselation::GuessStartingTriangle()
+void Tesselation::GuessStartingTriangle()
+{
         Info FunctionInfo(__func__);
 …
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
   list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
 …
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
     triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
 …
   // fill the set of neighbours
+  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+  Center.SubtractVector(TurningPoint->node);
+  set<TesselPoint*> SetOfNeighbours;
+  TesselPointSet SetOfNeighbours;
   SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
   TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
+  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
+  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
+    Log() << Verbose(0) << **TesselRunner << endl;
   TesselPointList::iterator Runner = connectedClosestPoints->begin();
   TesselPointList::iterator Sprinter = Runner;
 …
     AddTesselationPoint((*Sprinter), 2);
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
 …
     Runner = Sprinter;
     Sprinter++;
+    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
+    if (Sprinter != connectedClosestPoints->end())
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
+  }
   delete(connectedClosestPoints);
 …
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
   map<int, class BoundaryPointSet *> OrderMap;
   pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  PointMap OrderMap;
+  PointTestPair OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
 …
 };
+/** Finds the boundary points that are closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
+ */
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+  PointMap::const_iterator FindPoint;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  if (LinesOnBoundary.empty()) {
+    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
+    return 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;
+  DistanceToPointMap * points = new DistanceToPointMap;
+  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(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
+            }
+          }
+        } 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()) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    delete(points);
+    return NULL;
+  }
+  return points;
+};
+/** Finds the boundary line that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return closest BoundaryLineSet or NULL in degenerate case.
+ */
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+  // get closest points
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: 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 (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->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;
+      }
+    }
+  }
+  delete(points);
+  // 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;
+  }
+  return ClosestLine;
+};
 /** Finds the triangle that is closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
+ * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
+ */
+list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
+  TesselPoint *trianglePoints[3];
+  TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
+  if (LinesOnBoundary.empty()) {
+    eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first.";
+ * \return BoundaryTriangleSet of nearest triangle or NULL.
+ */
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
+        // get closest points
+        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << 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) {
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
+  LineSet ClosestLines;
+  double MinDistance = 1e+16;
+  Vector BaseLineIntersection;
+  Vector Center;
+  Vector BaseLine;
+  Vector BaseLineCenter;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthBase = BaseLine.NormSquared();
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      const double lengthEndA = BaseLineIntersection.NormSquared();
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthEndB = BaseLineIntersection.NormSquared();
+      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
+        const double lengthEnd = Min(lengthEndA, lengthEndB);
+        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
+          ClosestLines.clear();
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = lengthEnd;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
+        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+          ClosestLines.insert(LineRunner->second);
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
+        } else { // line is worse
+          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;
+        }
+      } else { // intersection is closer, calculate
+        // calculate closest point on line to desired point
+        BaseLineIntersection.CopyVector(x);
+        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        Center.CopyVector(&BaseLineIntersection);
+        Center.ProjectOntoPlane(&BaseLine);
+        BaseLineIntersection.SubtractVector(&Center);
+        const double distance = BaseLineIntersection.NormSquared();
+        if (Center.NormSquared() > BaseLine.NormSquared()) {
+          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
+        }
+        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 " << *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;
+        }
+      }
+    }
+  }
+  delete(points);
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLines.empty()) {
     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;
+  TriangleList * candidates = new TriangleList;
+  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
+    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
+    candidates->push_back(Runner->second);
+  }
+  return candidates;
 };
 …
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
+{
         Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
+  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
+  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
+  TriangleList 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;
+  } 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;
+      }
+  // go through all and pick the one with the best alignment to x
+  double MinAlignment = 2.*M_PI;
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+    (*Runner)->GetCenter(&Center);
+    helper.CopyVector(x);
+    helper.SubtractVector(&Center);
+    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
+    if (Alignment < MinAlignment) {
+      result = *Runner;
+      MinAlignment = Alignment;
+      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
+    } else {
+      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
+    }
+  }
   delete(triangles);
   return result;
 };
+/** 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. 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
+{
+        Info FunctionInfo(__func__);
+  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
+ * @return >0 if outside, ==0 if on surface, <0 if inside
+ */
+double Tesselation::GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const
+{
+  Info FunctionInfo(__func__);
   Vector Center;
+  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;
+  }
+  result->GetCenter(&Center);
+  Vector helper;
+  Vector DistanceToCenter;
+  Vector Intersection;
+  double distance = 0.;
+  if (triangle == NULL) {// is boundary point or only point in point cloud?
+    Log() << Verbose(1) << "No triangle given!" << endl;
+    return -1.;
+  } else {
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl;
+  }
+  triangle->GetCenter(&Center);
   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;
+  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
+    Log() << Verbose(1) << Point << " is an inner point." << endl;
+    return true;
+  DistanceToCenter.CopyVector(&Center);
+  DistanceToCenter.SubtractVector(&Point);
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
+  // check whether we are on boundary
+  if (fabs(DistanceToCenter.ScalarProduct(&triangle->NormalVector)) < MYEPSILON) {
+    // calculate whether inside of triangle
+    DistanceToCenter.CopyVector(&Point);
+    Center.CopyVector(&Point);
+    Center.SubtractVector(&triangle->NormalVector); // points towards MolCenter
+    DistanceToCenter.AddVector(&triangle->NormalVector); // points outside
+    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
+    if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
+      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
+      return 0.;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
+      return false;
+    }
   } else {
+    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
+    return false;
+  }
+}
+/** Checks whether the provided TesselPoint is within the tesselation structure.
+ *
+ * @param *Point of which to check the position
+ * @param *LC Linked Cell structure
+ *
+ * @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);
+}
+    // calculate smallest distance
+    distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
+    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
+    // then check direction to boundary
+    if (DistanceToCenter.ScalarProduct(&triangle->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;
+    }
+  }
+};
+/** Calculates distance to a tesselated surface.
+ * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
+ * \param &Point point to calculate distance from
+ * \param *LC needed for finding closest points fast
+ * \return distance squared to closest point on surface
+ */
+double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  BoundaryTriangleSet *triangle = FindClosestTriangleToVector(&Point, LC);
+  const double distance = GetDistanceSquaredToTriangle(Point, triangle);
+  return Min(distance, LC->RADIUS);
+};
 /** Gets all points connected to the provided point by triangulation lines.
 …
  * @return set of the all points linked to the provided one
  */
 set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
         Info FunctionInfo(__func__);
   set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
+TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+        TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
 …
+  }
   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;
 …
  * @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
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+{
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector center;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector PlaneNormal;
   Vector AngleZero;
   Vector OrthogonalVector;
   Vector helper;
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  TriangleList *triangles = NULL;
   if (SetOfNeighbours == NULL) {
 …
   // calculate central point
   for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
   //  << "; scale factor " << 1.0/connectedPoints.size();
   center.Scale(1.0/SetOfNeighbours->size());
   Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
   // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
   PlaneNormal.CopyVector(Point->node);
   PlaneNormal.SubtractVector(&center);
+  triangles = FindTriangles(TrianglePoints);
+  if ((triangles != NULL) && (!triangles->empty())) {
+    for (TriangleList::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();
-  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
   // construct one orthogonal vector
 …
   // go through all connected points and calculate angle
   for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
 …
+}
+/** 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 *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
+ */
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+{
+  Info FunctionInfo(__func__);
+  map<double, TesselPoint*> anglesOfPoints;
+  TesselPointList *connectedCircle = new TesselPointList;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+  if (SetOfNeighbours == NULL) {
+    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    delete(connectedCircle);
+    return NULL;
+  }
+  // check whether there's something to do
+  if (SetOfNeighbours->size() < 3) {
+    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+      connectedCircle->push_back(*TesselRunner);
+    return connectedCircle;
+  }
+  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
+  // calculate central point
+  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
+  TesselB++;
+  TesselC++;
+  TesselC++;
+  int counter = 0;
+  while (TesselC != SetOfNeighbours->end()) {
+    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
+    counter++;
+    TesselA++;
+    TesselB++;
+    TesselC++;
+    PlaneNormal.AddVector(&helper);
+  }
+  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
+  //  << "; scale factor " << counter;
+  PlaneNormal.Scale(1.0/(double)counter);
+//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
+//
+//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
+//  PlaneNormal.CopyVector(Point->node);
+//  PlaneNormal.SubtractVector(&center);
+//  PlaneNormal.Normalize();
+  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
+  // 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
+  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
+  for (TesselPointSet::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);
+    if (angle > M_PI) // the correction is of no use here (and not desired)
+      angle = 2.*M_PI - angle;
+    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (!InserterTest.second) {
+      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      performCriticalExit();
+    }
+  }
+  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 walk along a closed path.
+ *
 …
  * @return list of the all points linked to the provided one
  */
 list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
   list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
   list<TesselPoint*> *connectedPath = NULL;
+  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
 …
       } else if (!LineRunner->second) {
         LineRunner->second = true;
         connectedPath = new list<TesselPoint*>;
+        connectedPath = new TesselPointList;
         triangle = NULL;
         CurrentLine = runner->second;
 …
  * @return list of the closed paths
  */
 list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
         Info FunctionInfo(__func__);
   list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
   list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
   list<TesselPoint*> *connectedPath = NULL;
   list<TesselPoint*> *newPath = NULL;
+ListOfTesselPointList * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+  list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
+  TesselPointList *newPath = NULL;
   int count = 0;
   list<TesselPoint*>::iterator CircleRunner;
   list<TesselPoint*>::iterator CircleStart;
   for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  TesselPointList::iterator CircleRunner;
+  TesselPointList::iterator CircleStart;
+  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 …
     // go through list, look for reappearance of starting Point and create list
     list<TesselPoint*>::iterator Marker = CircleStart;
+    TesselPointList::iterator Marker = CircleStart;
     for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
         Log() << Verbose(1) << count+1 << ". closed path consists of: ";
         newPath = new list<TesselPoint*>;
         list<TesselPoint*>::iterator CircleSprinter = Marker;
+        newPath = new TesselPointList;
+        TesselPointList::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
 …
  * \return pointer to allocated list of triangles
  */
 set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
         Info FunctionInfo(__func__);
   set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+        Info FunctionInfo(__func__);
+        TriangleSet *connectedTriangles = new TriangleSet;
   if (Point == NULL) {
 …
+  }
   list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
   list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
   map<class BoundaryTriangleSet *, int> Candidates;
+  TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
       Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
+    }
+  }
 …
   count=0;
   NormalVector.Zero();
   for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
     Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
     NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
     RemoveTesselationTriangle(Runner->first);
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
+    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->second);
     count++;
+  }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
   list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
   list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
   map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
   list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+  TriangleMap::iterator NumberRunner = Candidates.begin();
+  TesselPointList::iterator StartNode, MiddleNode, EndNode;
   double angle;
   double smallestangle;
 …
       // re-create all triangles by going through connected points list
       list<class BoundaryLineSet *> NewLines;
+      LineList NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
 …
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
         list<class BoundaryLineSet *>::iterator Candidate;
+        LineList::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
           for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
 …
  * 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,
  *         will usually be one, in case of degeneration, there will be two
  */
 list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
         Info FunctionInfo(__func__);
   list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+        TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   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;
+  }
 …
  *         in the list, once as key and once as value
  */
 map<int, int> * Tesselation::FindAllDegeneratedLines()
+IndexToIndex * Tesselation::FindAllDegeneratedLines()
+{
         Info FunctionInfo(__func__);
         UniqueLines AllLines;
   map<int, int> * DegeneratedLines = new map<int, int>;
+  IndexToIndex * DegeneratedLines = new IndexToIndex;
   // sanity check
 …
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
   map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
 …
  *         in the list, once as key and once as value
  */
 map<int, int> * Tesselation::FindAllDegeneratedTriangles()
+{
         Info FunctionInfo(__func__);
   map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
   map<int, int> * DegeneratedTriangles = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
+{
+        Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
+  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
 …
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
   for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
 …
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
   map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
 …
+{
         Info FunctionInfo(__func__);
   map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
   for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
 …
   // 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;
 …
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
   map <int, int> *DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
 …
+  }
   map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
   map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;

src/tesselation.hpp

-              rbd61b41
+              r271e17
 // ======================================================= some template functions =========================================
+#define IndexToIndex map <int, int>
 #define PointMap map < int, class BoundaryPointSet * >
 #define PointSet set < class BoundaryPointSet * >
 …
 #define PolygonList list < class BoundaryPolygonSet * >
+#define DistanceToPointMap multimap <double, class BoundaryPointSet * >
+#define DistanceToPointPair pair <double, class BoundaryPointSet * >
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
 #define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
 …
 #define TesselPointList list <TesselPoint *>
 #define TesselPointSet set <TesselPoint *>
+#define ListOfTesselPointList list<list <TesselPoint *> *>
 /********************************************** declarations *******************************/
 …
   public:
     BoundaryPointSet();
     BoundaryPointSet(TesselPoint * Walker);
+    BoundaryPointSet(TesselPoint * const Walker);
     ~BoundaryPointSet();
     void AddLine(class BoundaryLineSet *line);
+    void AddLine(BoundaryLineSet * const line);
     LineMap lines;
 …
   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];
 …
   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];
 …
   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 {};
 …
     double PickFarthestofTwoBaselines(class BoundaryLineSet *Base);
     class BoundaryPointSet *IsConvexRectangle(class BoundaryLineSet *Base);
     map<int, int> * FindAllDegeneratedTriangles();
     map<int, int> * FindAllDegeneratedLines();
+    IndexToIndex * FindAllDegeneratedTriangles();
+    IndexToIndex * FindAllDegeneratedLines();
     void RemoveDegeneratedTriangles();
     void AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC);
     int CorrectAllDegeneratedPolygons();
+    set<TesselPoint*> * GetAllConnectedPoints(const TesselPoint* const Point) const;
+    set<BoundaryTriangleSet*> *GetAllTriangles(const BoundaryPointSet * const Point) const;
+    list<list<TesselPoint*> *> * GetPathsOfConnectedPoints(const TesselPoint* const Point) const;
+    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;
+    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;
+    TesselPointSet * GetAllConnectedPoints(const TesselPoint* const Point) const;
+    TriangleSet * GetAllTriangles(const BoundaryPointSet * const Point) const;
+    ListOfTesselPointList * GetPathsOfConnectedPoints(const TesselPoint* const Point) const;
+    ListOfTesselPointList * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
+    TesselPointList * GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
+    TesselPointList * GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
+    class BoundaryPointSet * GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
+    TriangleList * FindTriangles(const TesselPoint* const Points[3]) const;
+    TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
+    BoundaryTriangleSet * FindClosestTriangleToVector(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;
+    double GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) 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;
+    BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const;
     // print for debugging

src/tesselationhelpers.cpp

-              rbd61b41
+              r271e17
   if (fabs(HalfplaneIndicator) < MYEPSILON)
+    {
       if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(AlternativeDirection) >0 && AlternativeIndicator <0))
+        {
           TempNormal.Scale(-1);
 …
   else
+    {
       if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+      if (((TempNormal.ScalarProduct(Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(Direction)>0) && (HalfplaneIndicator<0)))
+        {
           TempNormal.Scale(-1);
 …
  * @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__);
 …
  * @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__);
 …
             helper.CopyVector(Point);
             helper.SubtractVector((*Runner)->node);
             double currentNorm = helper. Norm();
+            double currentNorm = helper.NormSquared();
             if (currentNorm < distance) {
               secondDistance = distance;
 …
+    }
     *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

-              rbd61b41
+              r271e17
 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);

src/unittests/AnalysisCorrelationToPointUnitTest.cpp

rbd61b41	r271e17
11	11	#include <cppunit/extensions/TestFactoryRegistry.h>
12	12	#include <cppunit/ui/text/TestRunner.h>
	13
	14	#include <cstring>
13	15
14	16	#include "analysis_correlation.hpp"

src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

rbd61b41	r271e17
11	11	#include <cppunit/extensions/TestFactoryRegistry.h>
12	12	#include <cppunit/ui/text/TestRunner.h>
	13
	14	#include <cstring>
13	15
14	16	#include "analysis_correlation.hpp"

src/unittests/AnalysisPairCorrelationUnitTest.cpp

rbd61b41	r271e17
11	11	#include <cppunit/extensions/TestFactoryRegistry.h>
12	12	#include <cppunit/ui/text/TestRunner.h>
	13
	14	#include <cstring>
13	15
14	16	#include "analysis_correlation.hpp"

src/unittests/Makefile.am

-              rbd61b41
+              r271e17
 AM_CXXFLAGS = $(CPPUNIT_CFLAGS)
+TESTS = ActOnAllUnitTest AnalysisBondsUnitTests AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest BondGraphUnitTest InfoUnitTest ListOfBondsUnitTest LogUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest StackClassUnitTest VectorUnitTest
+TESTS = \
+  ActOnAllUnitTest \
+  AnalysisBondsUnitTests \
+  AnalysisCorrelationToPointUnitTest \
+  AnalysisCorrelationToSurfaceUnitTest \
+  AnalysisPairCorrelationUnitTest \
+  BondGraphUnitTest \
+  GSLMatrixSymmetricUnitTest \
+  GSLMatrixUnitTest \
+  GSLVectorUnitTest \
+  InfoUnitTest \
+  LinearSystemOfEquationsUnitTest \
+  ListOfBondsUnitTest \
+  LogUnitTest \
+  MemoryUsageObserverUnitTest \
+  MemoryAllocatorUnitTest \
+  StackClassUnitTest \
+  TesselationUnitTest \
+  Tesselation_BoundaryTriangleUnitTest \
+  Tesselation_InOutsideUnitTest \
+  VectorUnitTest
 check_PROGRAMS = $(TESTS)
 noinst_PROGRAMS = $(TESTS)
 ActOnAllUnitTest_SOURCES = ../test/ActOnAllTest.hpp ActOnAllUnitTest.cpp ActOnAllUnitTest.hpp
 ActOnAllUnitTest_LDADD = ../libmolecuilder.a
+ActOnAllUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 AnalysisBondsUnitTests_SOURCES = analysisbondsunittest.cpp analysisbondsunittest.hpp
 AnalysisBondsUnitTests_LDADD = ../libmolecuilder.a
+AnalysisBondsUnitTests_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 AnalysisCorrelationToPointUnitTest_SOURCES = analysis_correlation.hpp AnalysisCorrelationToPointUnitTest.cpp AnalysisCorrelationToPointUnitTest.hpp
 AnalysisCorrelationToPointUnitTest_LDADD = ../libmolecuilder.a
+AnalysisCorrelationToPointUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 AnalysisCorrelationToSurfaceUnitTest_SOURCES = analysis_correlation.hpp AnalysisCorrelationToSurfaceUnitTest.cpp AnalysisCorrelationToSurfaceUnitTest.hpp
 AnalysisCorrelationToSurfaceUnitTest_LDADD = ../libmolecuilder.a
+AnalysisCorrelationToSurfaceUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 AnalysisPairCorrelationUnitTest_SOURCES = analysis_correlation.hpp AnalysisPairCorrelationUnitTest.cpp AnalysisPairCorrelationUnitTest.hpp
 AnalysisPairCorrelationUnitTest_LDADD = ../libmolecuilder.a
+AnalysisPairCorrelationUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 BondGraphUnitTest_SOURCES = bondgraphunittest.cpp bondgraphunittest.hpp
+BondGraphUnitTest_LDADD = ../libmolecuilder.a
+BondGraphUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+GSLMatrixSymmetricUnitTest_SOURCES = gslmatrixsymmetricunittest.cpp gslmatrixsymmetricunittest.hpp
+GSLMatrixSymmetricUnitTest_LDADD = ../libgslwrapper.a
+GSLMatrixUnitTest_SOURCES = gslmatrixunittest.cpp gslmatrixunittest.hpp
+GSLMatrixUnitTest_LDADD = ../libgslwrapper.a
+GSLVectorUnitTest_SOURCES = gslvectorunittest.cpp gslvectorunittest.hpp
+GSLVectorUnitTest_LDADD = ../libgslwrapper.a
 InfoUnitTest_SOURCES = infounittest.cpp infounittest.hpp
+InfoUnitTest_LDADD = ../libmolecuilder.a
+InfoUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+LinearSystemOfEquationsUnitTest_SOURCES = linearsystemofequationsunittest.cpp linearsystemofequationsunittest.hpp
+LinearSystemOfEquationsUnitTest_LDADD = ../libgslwrapper.a ../libmolecuilder.a
 ListOfBondsUnitTest_SOURCES = listofbondsunittest.cpp listofbondsunittest.hpp
 ListOfBondsUnitTest_LDADD = ../libmolecuilder.a
+ListOfBondsUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 LogUnitTest_SOURCES = logunittest.cpp logunittest.hpp
 LogUnitTest_LDADD = ../libmolecuilder.a
+LogUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 MemoryAllocatorUnitTest_SOURCES = memoryallocatorunittest.cpp memoryallocatorunittest.hpp
 MemoryAllocatorUnitTest_LDADD = ../libmolecuilder.a
+MemoryAllocatorUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 MemoryUsageObserverUnitTest_SOURCES = memoryusageobserverunittest.cpp memoryusageobserverunittest.hpp
 MemoryUsageObserverUnitTest_LDADD = ../libmolecuilder.a
+MemoryUsageObserverUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 StackClassUnitTest_SOURCES = stackclassunittest.cpp stackclassunittest.hpp
+StackClassUnitTest_LDADD = ../libmolecuilder.a
+StackClassUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+TesselationUnitTest_SOURCES = tesselationunittest.cpp tesselationunittest.hpp
+TesselationUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+Tesselation_BoundaryTriangleUnitTest_SOURCES = tesselation_boundarytriangleunittest.cpp tesselation_boundarytriangleunittest.hpp
+Tesselation_BoundaryTriangleUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+Tesselation_InOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
+Tesselation_InOutsideUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 VectorUnitTest_SOURCES = vectorunittest.cpp vectorunittest.hpp
 VectorUnitTest_LDADD = ../libmolecuilder.a
+VectorUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a

src/unittests/analysisbondsunittest.cpp

rbd61b41	r271e17
14	14	#include <iostream>
15	15	#include <stdio.h>
	16	#include <cstring>
16	17
17	18	#include "analysis_bonds.hpp"

src/unittests/bondgraphunittest.cpp

rbd61b41	r271e17
14	14	#include <iostream>
15	15	#include <stdio.h>
	16	#include <cstring>
16	17
17	18	#include "atom.hpp"

src/unittests/listofbondsunittest.cpp

rbd61b41	r271e17
11	11	#include <cppunit/extensions/TestFactoryRegistry.h>
12	12	#include <cppunit/ui/text/TestRunner.h>
	13
	14	#include <cstring>
13	15
14	16	#include "listofbondsunittest.hpp"

src/unittests/tesselationunittest.cpp

-              rbd61b41
+              r271e17
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+#include <cstring>
 #include "defs.hpp"
 …
   class TesselPoint *Walker;
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 0., 0.);
   Walker->Name = new char[3];
+  Walker->node = new Vector(1., 0., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "1");
   Walker->nr = 1;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 1., 0.);
   Walker->Name = new char[3];
+  Walker->node = new Vector(-1., 1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "2");
   Walker->nr = 2;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., -1., 0.);
   Walker->Name = new char[3];
+  Walker->node = new Vector(-1., -1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "3");
   Walker->nr = 3;
 …
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 0., 1.);
   Walker->Name = new char[3];
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "4");
   Walker->nr = 4;
 …
   // create tesselation
   TesselStruct = new Tesselation;
   TesselStruct->PointsOnBoundary.clear();
   TesselStruct->LinesOnBoundary.clear();
   TesselStruct->TrianglesOnBoundary.clear();
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->PointsOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->LinesOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->TrianglesOnBoundary.empty() );
   TesselStruct->FindStartingTriangle(SPHERERADIUS, LinkedList);
-  bool flag = false;
+  LineMap::iterator baseline = TesselStruct->LinesOnBoundary.begin();
+  while (baseline != TesselStruct->LinesOnBoundary.end()) {
+    if (baseline->second->triangles.size() == 1) {
+      flag = TesselStruct->FindNextSuitableTriangle(*(baseline->second), *(((baseline->second->triangles.begin()))->second), SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+  bool OneLoopWithoutSuccessFlag = true;
+  bool TesselationFailFlag = false;
+  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
+    // 2a. fill all new OpenLines
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      baseline = Runner->second;
+      if (baseline->pointlist.empty()) {
+        T = (((baseline->BaseLine->triangles.begin()))->second);
+        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+      }
+    }
+    baseline++;
+    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (flag)) {
+      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      flag = false;
+    // 2b. search for smallest ShortestAngle among all candidates
+    double ShortestAngle = 4.*M_PI;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+      }
+    }
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      OneLoopWithoutSuccessFlag = false;
+    else {
+      TesselStruct->AddCandidateTriangle(*baseline);
+    }
+  }
 …
   delete(TesselStruct);
   for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
-    delete[]((*Runner)->Name);
     delete((*Runner)->node);
     delete(*Runner);
+  }
   Corners.clear();
+};
+/** UnitTest for Tesselation::IsInnerPoint()
+ */
+void TesselationTest::IsInnerPointTest()
+{
+  // true inside points
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.), LinkedList) );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.5,0.,0.), LinkedList) );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.5,0.), LinkedList) );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.5), LinkedList) );
+  // corners
+  for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
+    CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint((*Runner), LinkedList) );
+  // true outside points
+  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,5.,0.), LinkedList) );
+  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,0.,5.), LinkedList) );
+  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(1.,1.,1.), LinkedList) );
+  // tricky point, there are three equally close triangles
+  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(5.,0.,0.), LinkedList) );
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
 };

src/unittests/tesselationunittest.hpp

-              rbd61b41
+              r271e17
+{
     CPPUNIT_TEST_SUITE( TesselationTest) ;
-    CPPUNIT_TEST ( IsInnerPointTest );
     CPPUNIT_TEST ( GetAllTrianglesTest );
     CPPUNIT_TEST ( ContainmentTest );
 …
       void setUp();
       void tearDown();
-      void IsInnerPointTest();
       void GetAllTrianglesTest();
       void ContainmentTest();

src/vector.cpp

-              rbd61b41
+              r271e17
 #include "defs.hpp"
 #include "helpers.hpp"
+#include "memoryallocator.hpp"
+#include "info.hpp"
+#include "gslmatrix.hpp"
 #include "leastsquaremin.hpp"
 #include "log.hpp"
+#include "memoryallocator.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_vector.h>
 /************************************ Functions for class vector ************************************/
 …
  * \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)
+{
+  Info FunctionInfo(__func__);
   double factor;
   Vector Direction, helper;
 …
   Direction.SubtractVector(Origin);
   Direction.Normalize();
+  //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
+  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;
+  }
 …
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
   if (factor < MYEPSILON) { // Origin is in-plane
     //Log() << Verbose(2) << "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;
 …
   Direction.Scale(factor);
   CopyVector(Origin);
   //Log() << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
   AddVector(&Direction);
 …
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
     //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
+    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
     return true;
   } else {
 …
 /** Calculates the intersection of the two lines that are both on the same plane.
+ * We construct auxiliary plane with its vector normal to one line direction and the PlaneNormal, then a vector
+ * from the first line's offset onto the plane. Finally, scale by factor is 1/cos(angle(line1,line2..)) = 1/SP(...), and
+ * project onto the first line's direction and add its offset.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
  * \param *out output stream for debugging
  * \param *Line1a first vector of first line
 …
 bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
+{
+  bool result = true;
+  Vector Direction, OtherDirection;
+  Vector AuxiliaryNormal;
+  Vector Distance;
+  const Vector *Normal = NULL;
+  Vector *ConstructedNormal = NULL;
+  bool FreeNormal = false;
+  // construct both direction vectors
+  Zero();
+  Direction.CopyVector(Line1b);
+  Direction.SubtractVector(Line1a);
+  if (Direction.IsZero())
+  Info FunctionInfo(__func__);
+  GSLMatrix *M = new GSLMatrix(4,4);
+  M->SetAll(1.);
+  for (int i=0;i<3;i++) {
+    M->Set(0, i, Line1a->x[i]);
+    M->Set(1, i, Line1b->x[i]);
+    M->Set(2, i, Line2a->x[i]);
+    M->Set(3, i, Line2b->x[i]);
+  }
+  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  //for (int i=0;i<4;i++) {
+  //  for (int j=0;j<4;j++)
+  //    cout << "\t" << M->Get(i,j);
+  //  cout << endl;
+  //}
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
     return false;
+  OtherDirection.CopyVector(Line2b);
+  OtherDirection.SubtractVector(Line2a);
+  if (OtherDirection.IsZero())
+  }
+  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
+  // constuct a,b,c
+  Vector a;
+  Vector b;
+  Vector c;
+  Vector d;
+  a.CopyVector(Line1b);
+  a.SubtractVector(Line1a);
+  b.CopyVector(Line2b);
+  b.SubtractVector(Line2a);
+  c.CopyVector(Line2a);
+  c.SubtractVector(Line1a);
+  d.CopyVector(Line2b);
+  d.SubtractVector(Line1b);
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
+   Zero();
+   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   return false;
+  }
+  // check for parallelity
+  Vector parallel;
+  double factor = 0.;
+  if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+    parallel.CopyVector(Line1a);
+    parallel.SubtractVector(Line2a);
+    factor = parallel.ScalarProduct(&a)/a.Norm();
+    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+      CopyVector(Line2a);
+      Log() << Verbose(1) << "Lines conincide." << endl;
+      return true;
+    } else {
+      parallel.CopyVector(Line1a);
+      parallel.SubtractVector(Line2b);
+      factor = parallel.ScalarProduct(&a)/a.Norm();
+      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+        CopyVector(Line2b);
+        Log() << Verbose(1) << "Lines conincide." << endl;
+        return true;
+      }
+    }
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    Zero();
     return false;
+  Direction.Normalize();
+  OtherDirection.Normalize();
+  //Log() << Verbose(4) << "INFO: Normalized Direction " << Direction << " and OtherDirection " << OtherDirection << "." << endl;
+  if (fabs(OtherDirection.ScalarProduct(&Direction) - 1.) < MYEPSILON) { // lines are parallel
+    if ((Line1a == Line2a) || (Line1a == Line2b))
+      CopyVector(Line1a);
+    else if ((Line1b == Line2b) || (Line1b == Line2b))
+        CopyVector(Line1b);
+    else
+      return false;
+    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
+    return true;
+  } else {
+    // check whether we have a plane normal vector
+    if (PlaneNormal == NULL) {
+      ConstructedNormal = new Vector;
+      ConstructedNormal->MakeNormalVector(&Direction, &OtherDirection);
+      Normal = ConstructedNormal;
+      FreeNormal = true;
+    } else
+      Normal = PlaneNormal;
+    AuxiliaryNormal.MakeNormalVector(&OtherDirection, Normal);
+    //Log() << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
+    Distance.CopyVector(Line2a);
+    Distance.SubtractVector(Line1a);
+    //Log() << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
+    if (Distance.IsZero()) {
+      // offsets are equal, match found
+      CopyVector(Line1a);
+      result = true;
+    } else {
+      CopyVector(Distance.Projection(&AuxiliaryNormal));
+      //Log() << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
+      double factor = Direction.ScalarProduct(&AuxiliaryNormal);
+      //Log() << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
+      Scale(1./(factor*factor));
+      //Log() << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
+      CopyVector(Projection(&Direction));
+      //Log() << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
+      if (this->IsZero())
+        result = false;
+      else
+        result = true;
+      AddVector(Line1a);
+    }
+    if (FreeNormal)
+      delete(ConstructedNormal);
+  }
+  if (result)
+    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
+  return result;
+  }
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1.CopyVector(&c);
+  temp1.VectorProduct(&b);
+  temp2.CopyVector(&a);
+  temp2.VectorProduct(&b);
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+  // construct intersection
+  CopyVector(&a);
+  Scale(s);
+  AddVector(Line1a);
+  Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
+  return true;
 };

tests/regression/Tesselation/1/post/NonConvexEnvelope.dat

rbd61b41	r271e17
20	20	5 6 8
21	21	6 7 8
	22	2 7 8
22	23	1 2 7
23		~~2 7 8~~
24	24	1 2 3

tests/regression/Tesselation/1/post/NonConvexEnvelope.r3d

rbd61b41	r271e17
48	48	-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
49	49	1
	50	1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
	51	1
50	52	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
51		1
52		~~1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
53	53	1
54	54	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.

tests/regression/Tesselation/2/post/ConvexEnvelope.dat

rbd61b41	r271e17
20	20	5 6 8
21	21	6 7 8
	22	2 7 8
22	23	1 2 7
23		~~2 7 8~~
24	24	1 2 3

tests/regression/Tesselation/2/post/ConvexEnvelope.r3d

rbd61b41	r271e17
48	48	-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
49	49	1
	50	1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
	51	1
50	52	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
51		1
52		~~1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
53	53	1
54	54	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.

tests/regression/Tesselation/2/post/NonConvexEnvelope.dat

rbd61b41	r271e17
20	20	5 6 8
21	21	6 7 8
	22	2 7 8
22	23	1 2 7
23		~~2 7 8~~
24	24	1 2 3

tests/regression/Tesselation/2/post/NonConvexEnvelope.r3d

rbd61b41	r271e17
48	48	-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
49	49	1
	50	1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
	51	1
50	52	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
51		1
52		~~1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
53	53	1
54	54	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.

Context Navigation

Changes in / [bd61b41:271e17]

Legend:

Download in other formats: