Diff [0af58ff89def082f4f9f59de50b950e56eddc91a:551a58e29f3a436c0de09a3afdffe11e4563f894] for / – MoleCuilder

.gitignore

r0af58f	r551a58
25	25
26	26	# ctags files
27		tags
	27	tags
28	28
29	29	# exclude automake generated stuff

molecuilder/src/Makefile.am

-              r0af58f
+              r551a58
 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
 …
 FORCE:
 $(srcdir)/.git-version: FORCE
         @if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe --tags HEAD) > .git-version-t 2>/dev/null \
+        @if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe HEAD) > .git-version-t 2>/dev/null \
           && ! diff .git-version-t $(srcdir)/.git-version >/dev/null 2>&1; then \
           mv -f .git-version-t $(srcdir)/.git-version; \

molecuilder/src/analysis_correlation.cpp

-              r0af58f
+              r551a58
 #include "analysis_correlation.hpp"
 #include "element.hpp"
+#include "info.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
 …
 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
+{
+  Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
 …
 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
+{
+  Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
 …
 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
+{
+  Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
 …
 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
+{
+  Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
 …
 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
+{
+  Info FunctionInfo(__func__);
   CorrelationToSurfaceMap *outmap = NULL;
   double distance = 0;
 …
         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);
 …
 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
+{
+  Info FunctionInfo(__func__);
   CorrelationToSurfaceMap *outmap = NULL;
   double distance = 0;
 …
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
                 triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
+                triangle = Surface->FindClosestTriangleToVector(&checkX, LC );
                 if (triangle != NULL) {
                   distance = DistanceToTrianglePlane(&checkX, triangle);
 …
 double GetBin ( const double value, const double BinWidth, const double BinStart )
+{
+  Info FunctionInfo(__func__);
   double bin =(double) (floor((value - BinStart)/BinWidth));
   return (bin*BinWidth+BinStart);
 …
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
+{
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
 …
 void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
+{
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
 …
 void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
+{
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
   for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
 …
 void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
+{
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tTriangle" << endl;
   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {

molecuilder/src/analyzer.cpp

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

molecuilder/src/bondgraph.cpp

r0af58f	r551a58
35	35	/** Parses the bond lengths in a given file and puts them int a matrix form.
36	36	* Allocates \a MatrixContainer for BondGraph::BondLengthMatrix, using MatrixContainer::ParseMatrix(),
37		* but only if parsing is successfull. Otherwise variable is left as NULL.
	37	* but only if parsing is successful. Otherwise variable is left as NULL.
38	38	* \param *out output stream for debugging
39	39	* \param filename file with bond lengths to parse

molecuilder/src/boundary.cpp

-              r0af58f
+              r551a58
  * \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;
   // 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]);
 …
 //  TesselStruct->RemoveDegeneratedTriangles();
+  // check envelope for consistency
+  status = CheckListOfBaselines(TesselStruct);
+  // store before correction
+  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
   // correct degenerated polygons
   TesselStruct->CorrectAllDegeneratedPolygons();

molecuilder/src/boundary.hpp

-              r0af58f
+              r551a58
 #define DEBUG 1
 #define DoSingleStepOutput 0
 #define SingleStepWidth 1
+#define SingleStepWidth 10
 #define DistancePair pair < double, atom* >
 …
 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);

molecuilder/src/builder.cpp

-              r0af58f
+              r551a58
 using namespace std;
+#include <cstring>
 #include "analysis_correlation.hpp"
 …
             Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
             Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
             Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
+            Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
             Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
             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 <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;
+            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
             Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
             Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
 …
                 else {
                   Log() << Verbose(2) << "File found and parsed." << endl;
-                  // @TODO rather do the dissection afterwards
-//                  mol->SetNameFromFilename(argv[argptr]);
-//                  molecules->ListOfMolecules.remove(mol);
-//                  molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration);
-//                  delete(mol);
-//                  if (molecules->ListOfMolecules.size() != 0) {
-//                    for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
-//                      if ((*ListRunner)->ActiveFlag) {
-//                        mol = *ListRunner;
-//                        break;
-//                      }
-//                  }
                   configPresent = present;
+                }
 …
               //argptr+=1;
               break;
+            case 'I':
+              Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;
+              // @TODO rather do the dissection afterwards
+              molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration);
+              mol = NULL;
+              if (molecules->ListOfMolecules.size() != 0) {
+                for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+                  if ((*ListRunner)->ActiveFlag) {
+                    mol = *ListRunner;
+                    break;
+                  }
+              }
+              if (mol == NULL) {
+                mol = *(molecules->ListOfMolecules.begin());
+                mol->ActiveFlag = true;
+              }
+              break;
             case 'C':
               if (ExitFlag == 0) ExitFlag = 1;
 …
                 performCriticalExit();
               } else {
-                SaveFlag = false;
                 ofstream output(argv[argptr+1]);
                 ofstream binoutput(argv[argptr+2]);
 …
                 counter = 0;
                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                   Actives[counter] = (*BigFinder)->ActiveFlag;
+                  Actives[counter++] = (*BigFinder)->ActiveFlag;
                   (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
+                }
 …
                 binoutput.close();
                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
                   (*BigFinder)->ActiveFlag = Actives[counter];
+                  (*BigFinder)->ActiveFlag = Actives[counter++];
                 Free(&Actives);
                 delete(LCList);
 …
             case 'F':
               if (ExitFlag == 0) ExitFlag = 1;
               if (argptr+5 >=argc) {
+              if (argptr+6 >=argc) {
                 ExitFlag = 255;
                 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
+                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
                 performCriticalExit();
               } else {
 …
                 molecule *filler = new molecule(periode);;
                 molecule *Filling = NULL;
                 atom *second = NULL, *third = NULL;
+//                atom *second = NULL, *third = NULL;
                 first = new atom();
                 first->type = periode->FindElement(1);
                 first->x.Init(0.441, -0.143, 0.);
+                first->type = periode->FindElement(5);
+                first->x.Zero();
                 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);
+//                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);
                 // call routine
                 double distance[NDIM];
                 for (int i=0;i<NDIM;i++)
                   distance[i] = atof(argv[argptr+i]);
                 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
+                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
                 if (Filling != NULL) {
+                  Filling->ActiveFlag = false;
                   molecules->insert(Filling);
+                }
 …
                 if (volume != -1)
                   ExitFlag = 255;
                   eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
+                  eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;
                   performCriticalExit();
               } else {

molecuilder/src/config.cpp

-              r0af58f
+              r551a58
 #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"
 …
   // don't do this here ...
   //MolList->DissectMoleculeIntoConnectedSubgraphs(mol,this);
+  delete(mol);
+  //delete(mol);
   delete(FileBuffer);
 };
 …
+              }
               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

molecuilder/src/helpers.hpp

-              r0af58f
+              r551a58
   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;
 };

molecuilder/src/joiner.cpp

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

molecuilder/src/memoryallocator.hpp

r0af58f	r551a58
16	16	#endif
17	17
	18	#include <cstdlib>
18	19	#include <iostream>
19	20	#include <iomanip>

molecuilder/src/memoryusageobserver.cpp

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

molecuilder/src/molecule.cpp

-              r0af58f
+              r551a58
+ *
  */
+#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);

molecuilder/src/molecule.hpp

r0af58f	r551a58
110	110	TesselPoint *GetPoint() const ;
111	111	TesselPoint *GetTerminalPoint() const ;
	112	int GetMaxId() const;
112	113	void GoToNext() const ;
113	114	void GoToPrevious() const ;

molecuilder/src/molecule_dynamics.cpp

-              r0af58f
+              r551a58
   // 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 {

molecuilder/src/molecule_fragmentation.cpp

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

molecuilder/src/molecule_geometry.cpp

-              r0af58f
+              r551a58
+{
   int Num = 0;
   atom *ptr = start->next;  // start at first in list
+  atom *ptr = start;  // start at first in list
   Center.Zero();
   if (ptr != end) {   //list not empty?
+  if (ptr->next != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
       ptr = ptr->next;

molecuilder/src/molecule_graph.cpp

-              r0af58f
+              r551a58
   bool status = true;
   if (ReferenceStack->IsEmpty()) {
+    eLog() << Verbose(0) << "ReferenceStack is empty!" << endl;
+    performCriticalExit();
+    Log() << Verbose(1) << "ReferenceStack is empty!" << endl;
     return false;
+  }

molecuilder/src/molecule_pointcloud.cpp

-              r0af58f
+              r551a58
 };
+/** 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.

molecuilder/src/moleculelist.cpp

-              r0af58f
+              r551a58
+ *
  */
+#include <cstring>
 #include "atom.hpp"
 …
   input.open(line.c_str());
   if (input == NULL) {
+    eLog() << Verbose(0) << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+    performCriticalExit();
+    Log() << Verbose(1) << endl << "Unable to open " << line << ", is the directory correct?" << endl;
     return false;
+  }

molecuilder/src/parser.cpp

-              r0af58f
+              r551a58
 // ======================================= INCLUDES ==========================================
+#include <cstring>
 #include "helpers.hpp"
 …
   //Log() << Verbose(0) << "Opening " << name << " ... "  << input << endl;
   if (input == NULL) {
     eLog() << Verbose(0) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
     performCriticalExit();
+    eLog() << Verbose(1) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
+    //performCriticalExit();
     return false;
+  }

molecuilder/src/periodentafel.cpp

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

molecuilder/src/tesselation.cpp

-              r0af58f
+              r551a58
  * \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 *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))) {
 …
   bool insertNewLine = true;
+  if (a->lines.find(b->node->nr) != a->lines.end()) {
+    LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  if (FindLine != a->lines.end()) {
+    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
     for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
 …
   double maxCoordinate[NDIM];
   BoundaryLineSet BaseLine;
-  Vector Oben;
   Vector helper;
   Vector Chord;
   Vector SearchDirection;
+  Oben.Zero();
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector SphereCenter;
+  Vector NormalVector;
+  NormalVector.Zero();
   for (i = 0; i < 3; i++) {
 …
   BTS = NULL;
   for (int k=0;k<NDIM;k++) {
     Oben.Zero();
     Oben.x[k] = 1.;
+    NormalVector.Zero();
+    NormalVector.x[k] = 1.;
     BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
 …
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
     FindSecondPointForTesselation(BaseLine.endpoints[0]->node, Oben, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     if (Temporary == NULL)  // have we found a second point?
       continue;
     BaseLine.endpoints[1] = new BoundaryPointSet(Temporary);
+    helper.CopyVector(BaseLine.endpoints[0]->node->node);
+    helper.SubtractVector(BaseLine.endpoints[1]->node->node);
+    helper.Normalize();
+    Oben.ProjectOntoPlane(&helper);
+    Oben.Normalize();
+    helper.VectorProduct(&Oben);
+    // construct center of circle
+    CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);
+    CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);
+    CircleCenter.Scale(0.5);
+    // construct normal vector of circle
+    CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node);
+    CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);
+    double radius = CirclePlaneNormal.NormSquared();
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
+    NormalVector.Normalize();
     ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+    Chord.CopyVector(BaseLine.endpoints[0]->node->node); // bring into calling function
+    Chord.SubtractVector(BaseLine.endpoints[1]->node->node);
+    double radius = Chord.ScalarProduct(&Chord);
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+    helper.CopyVector(&Oben);
+    helper.Scale(CircleRadius);
+    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+    SphereCenter.CopyVector(&NormalVector);
+    SphereCenter.Scale(CircleRadius);
+    SphereCenter.AddVector(&CircleCenter);
+    // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
     // look in one direction of baseline for initial candidate
     SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
     // adding point 1 and point 2 and add the line between them
 …
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
     CandidateForTesselation OptCandidates(&BaseLine);
     FindThirdPointForTesselation(Oben, SearchDirection, helper, OptCandidates, NULL, RADIUS, LC);
+    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
     Log() << Verbose(0) << "List of third Points is:" << endl;
     for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
 …
   Vector CircleCenter;
   Vector CirclePlaneNormal;
   Vector OldSphereCenter;
+  Vector RelativeSphereCenter;
   Vector SearchDirection;
   Vector helper;
 …
   double radius, CircleRadius;
-  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl;
   for (int i=0;i<3;i++)
     if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node))
+    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) {
       ThirdNode = T.endpoints[i]->node;
+      break;
+    }
+  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl;
   // construct center of circle
 …
   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
   if (radius/4. < RADIUS*RADIUS) {
+    // construct relative sphere center with now known CircleCenter
+    RelativeSphereCenter.CopyVector(&T.SphereCenter);
+    RelativeSphereCenter.SubtractVector(&CircleCenter);
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
     Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    // construct old center
+    GetCenterofCircumcircle(&OldSphereCenter, *T.endpoints[0]->node->node, *T.endpoints[1]->node->node, *T.endpoints[2]->node->node);
+    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
+    helper.Scale(sqrt(RADIUS*RADIUS - radius));
+    OldSphereCenter.AddVector(&helper);
+    OldSphereCenter.SubtractVector(&CircleCenter);
+    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+    // construct SearchDirection
+    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
+    helper.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl;
+    // construct SearchDirection and an "outward pointer"
+    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
     Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
     if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
 …
     // add third point
     FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
   } else {
 …
   // 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);
-    Center.CopyVector(&CandidateLine.OptCenter);
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
 …
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     AddTesselationTriangle();
+    Center.Scale(-1.);
+    BTS->GetCenter(&Center);
+    Center.SubtractVector(&CandidateLine.OptCenter);
+    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
     BTS->GetNormalVector(Center);
 …
     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);
 …
   Vector NewNormalVector;   // normal vector of the Candidate's triangle
   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
+  Vector RelativeOldSphereCenter;
+  Vector NewPlaneCenter;
   double CircleRadius; // radius of this circle
   double radius;
+  double otherradius;
   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
 …
   CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
   // calculate squared radius TesselPoint *ThirdNode,f circle
   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
   if (radius/4. < RADIUS*RADIUS) {
     CircleRadius = RADIUS*RADIUS - radius/4.;
+  radius = CirclePlaneNormal.NormSquared()/4.;
+  if (radius < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius;
     CirclePlaneNormal.Normalize();
     //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
     // test whether old center is on the band's plane
     if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
       eLog() << Verbose(1) << "Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
       //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+      Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;
       // check SearchDirection
       //Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+      }
 …
                 // check for three unique points
                 Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << *(Candidate->node) << "." << endl;
+                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl;
                 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){
                   // construct both new centers
                   GetCenterofCircumcircle(&NewSphereCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
                   OtherNewSphereCenter.CopyVector(&NewSphereCenter);
                   if ((NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node))
                   && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  // find center on the plane
+                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
+                  Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
+                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)
+                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
                   ) {
-                    helper.CopyVector(&NewNormalVector);
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                    Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl;
                     if (radius < RADIUS*RADIUS) {
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
+                      if (fabs(radius - otherradius) > HULLEPSILON) {
+                        eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;
+                      }
+                      // construct both new centers
+                      NewSphereCenter.CopyVector(&NewPlaneCenter);
+                      OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
+                      helper.CopyVector(&NewNormalVector);
                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
                       Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
                       NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
                       // OtherNewSphereCenter is created by the same vector just in the other direction
                       helper.Scale(-1.);
                       OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
 …
                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                       alpha = min(alpha, Otheralpha);
                       // if there is a better candidate, drop the current list and add the new candidate
                       // otherwise ignore the new candidate and keep the list
 …
+                        }
+                      }
                     } else {
                       Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
 …
   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 (Note that distance can be at most LinkedCell::RADIUS.)
+ */
+double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
+  Vector helper;
+  Vector DistanceToCenter;
+  Vector Intersection;
+  double distance = 0.;
   if (result == NULL) {// is boundary point or only point in point cloud?
     Log() << Verbose(1) << Point << " is the only point in vicinity." << endl;
+    return false;
+    return LC->RADIUS;
+  } else {
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
+  }
   result->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(&result->NormalVector)) < MYEPSILON) {
+    // calculate whether inside of triangle
+    DistanceToCenter.CopyVector(&Point);
+    Center.CopyVector(&Point);
+    Center.SubtractVector(&result->NormalVector); // points towards MolCenter
+    DistanceToCenter.AddVector(&result->NormalVector); // points outside
+    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
+    if (result->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 = result->GetClosestPointInsideTriangle(&Point, &Intersection);
+    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
+    distance = Min(distance, (LC->RADIUS*LC->RADIUS));
+    // then check direction to boundary
+    if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
+      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
+      return -distance;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
+      return +distance;
+    }
+  }
+};
 /** Gets all points connected to the provided point by triangulation lines.
 …
  * @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
+list< TesselPointList *> * 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;
+list<TesselPointList *> * 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
+  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
 …
     for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)
       for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
+        TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
+        if (TriangleInsertionTester.second)
+          Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
+          TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
+          if (TriangleInsertionTester.second)
+            Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        } else {
+          Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl;
+        }
+      }
     // check whether there are two that are parallel
 …
     /// 4a. Gather all triangles of this polygon
     TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();
+    // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do.
+    if (T->size() == 2) {
+      Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;
+      delete(T);
+      continue;
+    }
+    // check whether number is even
+    // If this case occurs, we have to think about it!
+    // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
+    // connections to either polygon ...
+    if (T->size() % 2 != 0) {
+      eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl;
+      performCriticalExit();
+    }
     TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
 …
+  }
   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;

molecuilder/src/tesselation.hpp

-              r0af58f
+              r551a58
 // ======================================================= some template functions =========================================
+#define IndexToIndex map <int, int>
 #define PointMap map < int, class BoundaryPointSet * >
 #define PointSet set < class BoundaryPointSet * >
 …
 #define PolygonSet set < class BoundaryPolygonSet * >
 #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> >
 …
   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];
     class BoundaryLineSet *lines[3];
     Vector NormalVector;
+    Vector SphereCenter;
     int Nr;
 };
 …
   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 {};
 …
     list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
     list<TesselPoint*> * GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
+    list<TesselPoint*> * GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
     class BoundaryPointSet *GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
     list<BoundaryTriangleSet*> *FindTriangles(const TesselPoint* const Points[3]) const;
 …
     class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
     bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
     bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const;
+    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;
+    TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
+    BoundaryTriangleSet* FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const;
     // print for debugging

molecuilder/src/tesselationhelpers.cpp

-              r0af58f
+              r551a58
   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);
 …
   Vector helper;
   double radius, alpha;
+  helper.CopyVector(&NewSphereCenter);
+  Vector RelativeOldSphereCenter;
+  Vector RelativeNewSphereCenter;
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
+  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
+  helper.CopyVector(&RelativeNewSphereCenter);
   // test whether new center is on the parameter circle's plane
   if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
 …
     helper.ProjectOntoPlane(&CirclePlaneNormal);
+  }
   radius = helper.ScalarProduct(&helper);
+  radius = helper.NormSquared();
   // test whether the new center vector has length of CircleRadius
   if (fabs(radius - CircleRadius) > HULLEPSILON)
     eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
   alpha = helper.Angle(&OldSphereCenter);
+  alpha = helper.Angle(&RelativeOldSphereCenter);
   // make the angle unique by checking the halfplanes/search direction
   if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
     alpha = 2.*M_PI - alpha;
   //Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
   radius = helper.Distance(&OldSphereCenter);
+  Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&RelativeOldSphereCenter);
   helper.ProjectOntoPlane(&NormalVector);
   // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
   if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
     //Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
     return alpha;
   } else {
     //Log() << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    Log() << Verbose(1) << "INFO: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "." << endl;
     return 2.*M_PI;
+  }
 …
  * @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++)

molecuilder/src/tesselationhelpers.hpp

-              r0af58f
+              r551a58
 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);

molecuilder/src/unittests/AnalysisCorrelationToPointUnitTest.cpp

r0af58f	r551a58
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"

molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

r0af58f	r551a58
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"

molecuilder/src/unittests/AnalysisPairCorrelationUnitTest.cpp

r0af58f	r551a58
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"

molecuilder/src/unittests/Makefile.am

-              r0af58f
+              r551a58
 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

molecuilder/src/unittests/analysisbondsunittest.cpp

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

molecuilder/src/unittests/bondgraphunittest.cpp

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

molecuilder/src/unittests/listofbondsunittest.cpp

r0af58f	r551a58
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"

molecuilder/src/unittests/tesselationunittest.cpp

-              r0af58f
+              r551a58
 #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();
 };

molecuilder/src/unittests/tesselationunittest.hpp

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

molecuilder/src/vector.cpp

-              r0af58f
+              r551a58
 #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;
 };
 …
   else
     return false;
+};
+/** Checks whether vector is normal to \a *normal.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool Vector::IsEqualTo(const Vector * const a) const
+{
+  bool status = true;
+  for (int i=0;i<NDIM;i++) {
+    if (fabs(x[i] - a->x[i]) > MYEPSILON)
+      status = false;
+  }
+  return status;
 };

molecuilder/src/vector.hpp

r0af58f	r551a58
42	42	bool IsOne() const;
43	43	bool IsNormalTo(const Vector * const normal) const;
	44	bool IsEqualTo(const Vector * const a) const;
44	45
45	46	void AddVector(const Vector * const y);

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

-              r0af58f
+              r551a58
 TITLE = "3D CONVEX SHELL"
 VARIABLES = "X" "Y" "Z" "U"
 ZONE T="0- H08_ H09_ H11", N=8, E=12, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
+ZONE T="test", N=8, E=12, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
 .78209 2.64589 2.64589 0
 .78209 2.64589 4.42589 0
 …
 3 4
 4 7
+6 7
 4 5
+7 8
+2 8
+6 7
+5 6
 3 5
-5 6
-7 8
-2 3
 5 8
 6 8
+7 8
+7 8
+2 7
+2 3

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

r0af58f	r551a58
34	34	1.37419 -0.89433 -0.89 0.12489 1.24767 -0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
35	35	1
36		~~2.26414 0.364321 -4.44089e-16 0.12489 1.24767 -0.89 0.12489 1.24767 0.89 1. 0. 0.~~
37		1
38		~~1.37419 -0.89433 -0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
39		1
40		~~1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
41		1
42	36	0.12489 1.24767 -0.89 -2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 1. 0. 0.
43	37	1
44		~~1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16~~ 0.12489 1.24767 0.89 1. 0. 0.
	38	2.26414 0.364321 -4.44089e-16 0.12489 1.24767 -0.89 0.12489 1.24767 0.89 1. 0. 0.
45	39	1
46	40	0.12489 1.24767 -0.89 0.12489 1.24767 0.89 -2.01426 0.364321 -4.44089e-16 1. 0. 0.
47	41	1
48		-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
49		1
50		1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.
	42	1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 0.12489 1.24767 0.89 1. 0. 0.
51	43	1
52	44	1.37419 -0.89433 0.89 0.12489 1.24767 0.89 -1.12431 -0.89433 0.89 1. 0. 0.
53	45	1
54	46	0.12489 1.24767 0.89 -2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 0.89 1. 0. 0.
	47	1
	48	-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
	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
	52	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
	53	1
	54	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.
55	55	9
56	56	# terminating special property
…	…
59	59	25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0
60	60	2
61		~~-1.00456 0.23922 0.593333~~ 5 1 0 0
	61	1.67084 -0.47478 -8.88178e-16 5 1 0 0
62	62	9
63	63	terminating special property

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

-              r0af58f
+              r551a58
 TITLE = "3D CONVEX SHELL"
 VARIABLES = "X" "Y" "Z" "U"
 ZONE T="0- H08_ H09_ H11", N=8, E=12, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
+ZONE T="test", N=8, E=12, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
 .78209 2.64589 2.64589 0
 .78209 2.64589 4.42589 0
 …
 3 4
 4 7
+6 7
 4 5
+7 8
+2 8
+6 7
+5 6
 3 5
-5 6
-7 8
-2 3
 5 8
 6 8
+7 8
+7 8
+2 7
+2 3

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

r0af58f	r551a58
34	34	1.37419 -0.89433 -0.89 0.12489 1.24767 -0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
35	35	1
36		~~2.26414 0.364321 -4.44089e-16 0.12489 1.24767 -0.89 0.12489 1.24767 0.89 1. 0. 0.~~
37		1
38		~~1.37419 -0.89433 -0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
39		1
40		~~1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
41		1
42	36	0.12489 1.24767 -0.89 -2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 1. 0. 0.
43	37	1
44		~~1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16~~ 0.12489 1.24767 0.89 1. 0. 0.
	38	2.26414 0.364321 -4.44089e-16 0.12489 1.24767 -0.89 0.12489 1.24767 0.89 1. 0. 0.
45	39	1
46	40	0.12489 1.24767 -0.89 0.12489 1.24767 0.89 -2.01426 0.364321 -4.44089e-16 1. 0. 0.
47	41	1
48		-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
49		1
50		1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.
	42	1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 0.12489 1.24767 0.89 1. 0. 0.
51	43	1
52	44	1.37419 -0.89433 0.89 0.12489 1.24767 0.89 -1.12431 -0.89433 0.89 1. 0. 0.
53	45	1
54	46	0.12489 1.24767 0.89 -2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 0.89 1. 0. 0.
	47	1
	48	-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
	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
	52	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
	53	1
	54	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.
55	55	9
56	56	# terminating special property
…	…
59	59	25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0
60	60	2
61		~~-1.00456 0.23922 0.593333~~ 5 1 0 0
	61	1.67084 -0.47478 -8.88178e-16 5 1 0 0
62	62	9
63	63	terminating special property

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

-              r0af58f
+              r551a58
 TITLE = "3D CONVEX SHELL"
 VARIABLES = "X" "Y" "Z" "U"
 ZONE T="0- H08_ H09_ H11", N=8, E=12, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
+ZONE T="test", N=8, E=12, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
 .78209 2.64589 2.64589 0
 .78209 2.64589 4.42589 0
 …
 3 4
 4 7
+6 7
 4 5
+7 8
+2 8
+6 7
+5 6
 3 5
-5 6
-7 8
-2 3
 5 8
 6 8
+7 8
+7 8
+2 7
+2 3

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

r0af58f	r551a58
34	34	1.37419 -0.89433 -0.89 0.12489 1.24767 -0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
35	35	1
36		~~2.26414 0.364321 -4.44089e-16 0.12489 1.24767 -0.89 0.12489 1.24767 0.89 1. 0. 0.~~
37		1
38		~~1.37419 -0.89433 -0.89 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
39		1
40		~~1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 0.89 1. 0. 0.~~
41		1
42	36	0.12489 1.24767 -0.89 -2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 1. 0. 0.
43	37	1
44		~~1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16~~ 0.12489 1.24767 0.89 1. 0. 0.
	38	2.26414 0.364321 -4.44089e-16 0.12489 1.24767 -0.89 0.12489 1.24767 0.89 1. 0. 0.
45	39	1
46	40	0.12489 1.24767 -0.89 0.12489 1.24767 0.89 -2.01426 0.364321 -4.44089e-16 1. 0. 0.
47	41	1
48		-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
49		1
50		1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.
	42	1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 0.12489 1.24767 0.89 1. 0. 0.
51	43	1
52	44	1.37419 -0.89433 0.89 0.12489 1.24767 0.89 -1.12431 -0.89433 0.89 1. 0. 0.
53	45	1
54	46	0.12489 1.24767 0.89 -2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 0.89 1. 0. 0.
	47	1
	48	-2.01426 0.364321 -4.44089e-16 -1.12431 -0.89433 -0.89 -1.12431 -0.89433 0.89 1. 0. 0.
	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
	52	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 -1.12431 -0.89433 -0.89 1. 0. 0.
	53	1
	54	1.37419 -0.89433 -0.89 1.37419 -0.89433 0.89 2.26414 0.364321 -4.44089e-16 1. 0. 0.
55	55	9
56	56	# terminating special property
…	…
59	59	25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0
60	60	2
61		~~-1.00456 0.23922 0.593333~~ 5 1 0 0
	61	1.67084 -0.47478 -8.88178e-16 5 1 0 0
62	62	9
63	63	terminating special property

molecuilder/tests/regression/Tesselation/3/post/NonConvexEnvelope.dat

-              r0af58f
+              r551a58
 TITLE = "3D CONVEX SHELL"
 VARIABLES = "X" "Y" "Z" "U"
 ZONE T="0- H49_ H50_ H51", N=44, E=86, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
+ZONE T="test", N=44, E=86, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
 .9077 1.1106 0.1214 1
 .3612 -3.628 1.323 1
 …
 .8131 1.4776 2.5103 0
 .9137 2.2936 1.3739 0
 .159 2.5738 1.2698 2
+.159 2.5738 1.2698 5
 .6606 -0.4593 2.1396 2
 .2007 -1.4419 0.7311 4
 -3.3002 2.3589 0.0094 5
+-3.3002 2.3589 0.0094 8
 -4.377 1.6962 -1.2433 3
 .2593 1.4547 -1.7445 0
 …
 .2727 1.6068 1.2828 2
 -6.2394 4.6427 0.0632 0
 -4.4738 4.5591 -0.1458 1
+-4.4738 4.5591 -0.1458 3
 -5.5506 3.8964 -1.3985 0
 -6.7081 0.9923 0.6224 2
 …
 32 44
 32 35
-33 44
 34 35
 32 35
+32 33
+23 35
+17 35
+10 17
+8 10
+8 35
+4 35
+29 35
+34 35
+3 4
+4 29
+15 29
+28 29
+29 34
+7 15
+14 15
+15 28
+25 28
+28 37
+34 37
+34 37
+37 44
+26 37
+26 27
+27 33
 33 44
-34 37
-34 35
-23 35
-23 33
-27 33
 37 44
-37 44
-34 37
-29 34
-29 35
-18 23
-17 35
-27 33
-26 27
-26 37
-28 37
-28 29
-4 29
-4 35
-17 18
-10 17
-8 35
-18 27
 25 27
-25 28
-15 28
-15 29
-3 4
-10 17
-18 31
-8 10
-18 30
-27 30
 27 30
-14 15
-7 15
-3 5
-10 12
-19 31
-30 31
-10 12
 14 30
+24 30
+30 36
+36 39
+30 39
+7 24
 7 14
-5 7
-5 12
-12 13
-13 19
-19 31
-31 39
-24 30
-7 24
-7 11
-12 22
-13 21
-21 31
-30 39
-39 40
-30 36
 11 24
+11 22
+21 22
+31 43
+40 43
+42 43
+42 43
+41 43
+21 41
+20 24
 24 41
 36 41
 41 42
+20 22
+11 22
+7 11
+5 7
+30 39
+27 30
+18 30
+18 27
+27 33
+23 33
+38 39
 38 42
+30 31
+31 39
+39 40
+18 31
+17 18
+18 23
+19 31
+13 19
+19 31
+21 31
+31 43
+40 43
+12 13
+10 12
+10 17
+13 21
+21 22
+12 22
+5 12
+3 5
+10 12
+21 22
+21 41
+41 43
+42 43
+32 33
+33 44
+42 43
 40 42
 39 40
-38 39
-36 39
-30 39
-20 24
-20 22
-21 22

molecuilder/tests/regression/Tesselation/3/post/NonConvexEnvelope.r3d

-              r0af58f
+              r551a58
 .33693 1.04442 0.0886712     5.68853 1.38852 -1.77723        5.55043 -0.952879 -0.490929     1. 0. 0.
+.33693 1.04442 0.0886712     6.17523 -0.969279 1.17127       5.55043 -0.952879 -0.490929     1. 0. 0.
+.62023 1.25552 -2.86813      5.68853 1.38852 -1.77723        5.55043 -0.952879 -0.490929     1. 0. 0.
+.76663 -0.360379 -2.99373    3.62023 1.25552 -2.86813        5.55043 -0.952879 -0.490929     1. 0. 0.
+.19193 -1.51408 -0.867629    1.76663 -0.360379 -2.99373      5.55043 -0.952879 -0.490929     1. 0. 0.
+.19193 -1.51408 -0.867629    0.450934 -2.72908 -2.23353      1.76663 -0.360379 -2.99373      1. 0. 0.
+.917634 -3.66448 -0.484829   2.19193 -1.51408 -0.867629      0.450934 -2.72908 -2.23353      1. 0. 0.
+.917634 -3.66448 -0.484829   2.19193 -1.51408 -0.867629      5.55043 -0.952879 -0.490929     1. 0. 0.
+.917634 -3.66448 -0.484829   1.92773 -2.57738 0.498071       5.55043 -0.952879 -0.490929     1. 0. 0.
+.92773 -2.57738 0.498071     3.62993 -1.50808 0.698371       5.55043 -0.952879 -0.490929     1. 0. 0.
+.62993 -1.50808 0.698371     6.17523 -0.969279 1.17127       5.55043 -0.952879 -0.490929     1. 0. 0.
+.790434 -3.69418 1.29027     0.917634 -3.66448 -0.484829     1.92773 -2.57738 0.498071       1. 0. 0.
+.790434 -3.69418 1.29027     1.92773 -2.57738 0.498071       3.62993 -1.50808 0.698371       1. 0. 0.
+.790434 -3.69418 1.29027     2.06173 -1.39848 1.79787        3.62993 -1.50808 0.698371       1. 0. 0.
+.06173 -1.39848 1.79787      4.08983 -0.525479 2.10687       3.62993 -1.50808 0.698371       1. 0. 0.
+.08983 -0.525479 2.10687     3.62993 -1.50808 0.698371       6.17523 -0.969279 1.17127       1. 0. 0.
+.790434 -3.69418 1.29027     -0.287266 -1.67078 2.48017      2.06173 -1.39848 1.79787        1. 0. 0.
+  -0.287266 -1.67078 2.48017    1.46453 0.112321 2.50927        2.06173 -1.39848 1.79787        1. 0. 0.
+.46453 0.112321 2.50927      2.06173 -1.39848 1.79787        4.08983 -0.525479 2.10687       1. 0. 0.
+.46453 0.112321 2.50927      3.24233 1.41142 2.47757         4.08983 -0.525479 2.10687       1. 0. 0.
+.24233 1.41142 2.47757       4.08983 -0.525479 2.10687       5.70193 1.54062 1.25007         1. 0. 0.
+.08983 -0.525479 2.10687     6.17523 -0.969279 1.17127       5.70193 1.54062 1.25007         1. 0. 0.
+.33693 1.04442 0.0886712     6.17523 -0.969279 1.17127       5.70193 1.54062 1.25007         1. 0. 0.
+.33693 1.04442 0.0886712     5.70193 1.54062 1.25007         7.56833 1.97852 -0.00632877     1. 0. 0.
+.24233 1.41142 2.47757       4.34293 2.22742 1.34117         5.70193 1.54062 1.25007         1. 0. 0.
+.24233 1.41142 2.47757       4.34293 2.22742 1.34117         2.58823 2.50762 1.23707         1. 0. 0.
+.34293 2.22742 1.34117       2.58823 2.50762 1.23707         5.11553 2.70122 -0.710229       1. 0. 0.
+.34293 2.22742 1.34117       5.11553 2.70122 -0.710229       7.56833 1.97852 -0.00632877     1. 0. 0.
+.34293 2.22742 1.34117       5.70193 1.54062 1.25007         7.56833 1.97852 -0.00632877     1. 0. 0.
+.46453 0.112321 2.50927      3.24233 1.41142 2.47757         2.58823 2.50762 1.23707         1. 0. 0.
+.46453 0.112321 2.50927      2.58823 2.50762 1.23707         -2.87097 2.29272 -0.0233288     1. 0. 0.
+  -0.759066 -0.265179 1.48487   1.46453 0.112321 2.50927        -2.87097 2.29272 -0.0233288     1. 0. 0.
+  -0.759066 -0.265179 1.48487   -3.66127 0.565021 1.57007       -2.87097 2.29272 -0.0233288     1. 0. 0.
+  -3.66127 0.565021 1.57007     -2.87097 2.29272 -0.0233288     -4.83487 2.76522 1.41487        1. 0. 0.
+  -2.87097 2.29272 -0.0233288   -4.83487 2.76522 1.41487        -4.04457 4.49292 -0.178529      1. 0. 0.
+.58823 2.50762 1.23707       -2.87097 2.29272 -0.0233288     -4.04457 4.49292 -0.178529      1. 0. 0.
+  -0.759066 -0.265179 1.48487   -0.287266 -1.67078 2.48017      -3.66127 0.565021 1.57007       1. 0. 0.
+  -0.759066 -0.265179 1.48487   -0.287266 -1.67078 2.48017      1.46453 0.112321 2.50927        1. 0. 0.
+  -0.287266 -1.67078 2.48017    -2.45927 -1.63678 1.72157       -3.66127 0.565021 1.57007       1. 0. 0.
+  -2.45927 -1.63678 1.72157     -4.75167 -1.93408 0.935971      -3.66127 0.565021 1.57007       1. 0. 0.
+  -4.75167 -1.93408 0.935971    -3.66127 0.565021 1.57007       -6.27887 0.926121 0.589671      1. 0. 0.
+  -3.66127 0.565021 1.57007     -4.83487 2.76522 1.41487        -6.27887 0.926121 0.589671      1. 0. 0.
+  -4.83487 2.76522 1.41487      -6.27887 0.926121 0.589671      -7.11497 2.49352 0.479071       1. 0. 0.
+  -2.45927 -1.63678 1.72157     -4.75167 -1.93408 0.935971      -3.90927 -3.49908 0.839771      1. 0. 0.
+  -1.52417 -3.64138 0.503471    -2.45927 -1.63678 1.72157       -3.90927 -3.49908 0.839771      1. 0. 0.
+  -1.52417 -3.64138 0.503471    -0.287266 -1.67078 2.48017      -2.45927 -1.63678 1.72157       1. 0. 0.
+.790434 -3.69418 1.29027     -1.52417 -3.64138 0.503471      -0.287266 -1.67078 2.48017      1. 0. 0.
+.58823 2.50762 1.23707       -2.87097 2.29272 -0.0233288     -4.04457 4.49292 -0.178529      1. 0. 0.
+.15633 1.11082 0.326671      2.58823 2.50762 1.23707         -2.87097 2.29272 -0.0233288     1. 0. 0.
+.15633 1.11082 0.326671      1.03283 1.01972 -2.14533        -2.87097 2.29272 -0.0233288     1. 0. 0.
+.15633 1.11082 0.326671      1.03283 1.01972 -2.14533        2.58823 2.50762 1.23707         1. 0. 0.
+.03283 1.01972 -2.14533      2.58823 2.50762 1.23707         5.11553 2.70122 -0.710229       1. 0. 0.
+.03283 1.01972 -2.14533      3.62023 1.25552 -2.86813        5.11553 2.70122 -0.710229       1. 0. 0.
+  -4.83487 2.76522 1.41487      -5.81017 4.57652 0.0304712      -4.04457 4.49292 -0.178529      1. 0. 0.
+  -4.83487 2.76522 1.41487      -5.81017 4.57652 0.0304712      -7.11497 2.49352 0.479071       1. 0. 0.
+.03283 1.01972 -2.14533      -2.87097 2.29272 -0.0233288     -3.94777 1.63002 -1.27603       1. 0. 0.
+  -2.87097 2.29272 -0.0233288   -3.94777 1.63002 -1.27603       -4.04457 4.49292 -0.178529      1. 0. 0.
+  -3.94777 1.63002 -1.27603     -4.04457 4.49292 -0.178529      -5.12137 3.83022 -1.43123       1. 0. 0.
+  -0.573766 -1.42458 -2.91753   1.03283 1.01972 -2.14533        -3.94777 1.63002 -1.27603       1. 0. 0.
+  -0.573766 -1.42458 -2.91753   1.76663 -0.360379 -2.99373      1.03283 1.01972 -2.14533        1. 0. 0.
+.76663 -0.360379 -2.99373    1.03283 1.01972 -2.14533        3.62023 1.25552 -2.86813        1. 0. 0.
+  -0.573766 -1.42458 -2.91753   -2.77417 -0.570279 -1.12083     -3.94777 1.63002 -1.27603       1. 0. 0.
+  -0.573766 -1.42458 -2.91753   -2.49667 -2.18078 -1.79993      -2.77417 -0.570279 -1.12083     1. 0. 0.
+  -2.49667 -2.18078 -1.79993    -2.77417 -0.570279 -1.12083     -3.94777 1.63002 -1.27603       1. 0. 0.
+  -2.49667 -2.18078 -1.79993    -4.17327 -2.53828 -0.635229     -3.94777 1.63002 -1.27603       1. 0. 0.
+  -4.17327 -2.53828 -0.635229   -3.94777 1.63002 -1.27603       -6.42617 1.74722 -0.982629      1. 0. 0.
+  -3.94777 1.63002 -1.27603     -5.12137 3.83022 -1.43123       -6.42617 1.74722 -0.982629      1. 0. 0.
+  -0.573766 -1.42458 -2.91753   -1.62867 -3.74268 -1.79493      -2.49667 -2.18078 -1.79993      1. 0. 0.
+  -0.573766 -1.42458 -2.91753   0.450934 -2.72908 -2.23353      -1.62867 -3.74268 -1.79493      1. 0. 0.
+  -0.573766 -1.42458 -2.91753   0.450934 -2.72908 -2.23353      1.76663 -0.360379 -2.99373      1. 0. 0.
+  -1.62867 -3.74268 -1.79493    -2.49667 -2.18078 -1.79993      -4.17327 -2.53828 -0.635229     1. 0. 0.
+  -1.62867 -3.74268 -1.79493    -4.17327 -2.53828 -0.635229     -3.90927 -3.49908 0.839771      1. 0. 0.
+  -1.52417 -3.64138 0.503471    -1.62867 -3.74268 -1.79493      -3.90927 -3.49908 0.839771      1. 0. 0.
+.917634 -3.66448 -0.484829   -1.52417 -3.64138 0.503471      -1.62867 -3.74268 -1.79493      1. 0. 0.
+.790434 -3.69418 1.29027     0.917634 -3.66448 -0.484829     -1.52417 -3.64138 0.503471      1. 0. 0.
+.917634 -3.66448 -0.484829   0.450934 -2.72908 -2.23353      -1.62867 -3.74268 -1.79493      1. 0. 0.
+  -4.75167 -1.93408 0.935971    -4.17327 -2.53828 -0.635229     -3.90927 -3.49908 0.839771      1. 0. 0.
+  -4.75167 -1.93408 0.935971    -4.17327 -2.53828 -0.635229     -6.27887 0.926121 0.589671      1. 0. 0.
+  -4.17327 -2.53828 -0.635229   -6.27887 0.926121 0.589671      -6.42617 1.74722 -0.982629      1. 0. 0.
+  -6.27887 0.926121 0.589671    -7.11497 2.49352 0.479071       -6.42617 1.74722 -0.982629      1. 0. 0.
+.62023 1.25552 -2.86813      5.68853 1.38852 -1.77723        5.11553 2.70122 -0.710229       1. 0. 0.
 .68853 1.38852 -1.77723      5.11553 2.70122 -0.710229       7.56833 1.97852 -0.00632877     1. 0. 0.
-.33693 1.04442 0.0886712     6.17523 -0.969279 1.17127       5.55043 -0.952879 -0.490929     1. 0. 0.
-.62023 1.25552 -2.86813      5.68853 1.38852 -1.77723        5.55043 -0.952879 -0.490929     1. 0. 0.
-.62023 1.25552 -2.86813      5.68853 1.38852 -1.77723        5.11553 2.70122 -0.710229       1. 0. 0.
-.34293 2.22742 1.34117       5.11553 2.70122 -0.710229       7.56833 1.97852 -0.00632877     1. 0. 0.
-.33693 1.04442 0.0886712     6.17523 -0.969279 1.17127       5.70193 1.54062 1.25007         1. 0. 0.
-.62993 -1.50808 0.698371     6.17523 -0.969279 1.17127       5.55043 -0.952879 -0.490929     1. 0. 0.
-.76663 -0.360379 -2.99373    3.62023 1.25552 -2.86813        5.55043 -0.952879 -0.490929     1. 0. 0.
-.03283 1.01972 -2.14533      3.62023 1.25552 -2.86813        5.11553 2.70122 -0.710229       1. 0. 0.
-.34293 2.22742 1.34117       2.58823 2.50762 1.23707         5.11553 2.70122 -0.710229       1. 0. 0.
-.34293 2.22742 1.34117       5.70193 1.54062 1.25007         7.56833 1.97852 -0.00632877     1. 0. 0.
-.33693 1.04442 0.0886712     5.70193 1.54062 1.25007         7.56833 1.97852 -0.00632877     1. 0. 0.
-.08983 -0.525479 2.10687     6.17523 -0.969279 1.17127       5.70193 1.54062 1.25007         1. 0. 0.
-.08983 -0.525479 2.10687     3.62993 -1.50808 0.698371       6.17523 -0.969279 1.17127       1. 0. 0.
-.92773 -2.57738 0.498071     3.62993 -1.50808 0.698371       5.55043 -0.952879 -0.490929     1. 0. 0.
-.76663 -0.360379 -2.99373    1.03283 1.01972 -2.14533        3.62023 1.25552 -2.86813        1. 0. 0.
-.19193 -1.51408 -0.867629    1.76663 -0.360379 -2.99373      5.55043 -0.952879 -0.490929     1. 0. 0.
-.03283 1.01972 -2.14533      2.58823 2.50762 1.23707         5.11553 2.70122 -0.710229       1. 0. 0.
-.24233 1.41142 2.47757       4.34293 2.22742 1.34117         2.58823 2.50762 1.23707         1. 0. 0.
-.24233 1.41142 2.47757       4.34293 2.22742 1.34117         5.70193 1.54062 1.25007         1. 0. 0.
-.24233 1.41142 2.47757       4.08983 -0.525479 2.10687       5.70193 1.54062 1.25007         1. 0. 0.
-.06173 -1.39848 1.79787      4.08983 -0.525479 2.10687       3.62993 -1.50808 0.698371       1. 0. 0.
-.790434 -3.69418 1.29027     1.92773 -2.57738 0.498071       3.62993 -1.50808 0.698371       1. 0. 0.
-.917634 -3.66448 -0.484829   1.92773 -2.57738 0.498071       5.55043 -0.952879 -0.490929     1. 0. 0.
-  -0.573766 -1.42458 -2.91753   1.76663 -0.360379 -2.99373      1.03283 1.01972 -2.14533        1. 0. 0.
-.19193 -1.51408 -0.867629    0.450934 -2.72908 -2.23353      1.76663 -0.360379 -2.99373      1. 0. 0.
-.917634 -3.66448 -0.484829   2.19193 -1.51408 -0.867629      5.55043 -0.952879 -0.490929     1. 0. 0.
-.15633 1.11082 0.326671      1.03283 1.01972 -2.14533        2.58823 2.50762 1.23707         1. 0. 0.
-.46453 0.112321 2.50927      3.24233 1.41142 2.47757         2.58823 2.50762 1.23707         1. 0. 0.
-.46453 0.112321 2.50927      3.24233 1.41142 2.47757         4.08983 -0.525479 2.10687       1. 0. 0.
-.46453 0.112321 2.50927      2.06173 -1.39848 1.79787        4.08983 -0.525479 2.10687       1. 0. 0.
-.790434 -3.69418 1.29027     2.06173 -1.39848 1.79787        3.62993 -1.50808 0.698371       1. 0. 0.
-.790434 -3.69418 1.29027     0.917634 -3.66448 -0.484829     1.92773 -2.57738 0.498071       1. 0. 0.
-  -0.573766 -1.42458 -2.91753   0.450934 -2.72908 -2.23353      1.76663 -0.360379 -2.99373      1. 0. 0.
-  -0.573766 -1.42458 -2.91753   1.03283 1.01972 -2.14533        -3.94777 1.63002 -1.27603       1. 0. 0.
-.917634 -3.66448 -0.484829   2.19193 -1.51408 -0.867629      0.450934 -2.72908 -2.23353      1. 0. 0.
-.15633 1.11082 0.326671      1.03283 1.01972 -2.14533        -2.87097 2.29272 -0.0233288     1. 0. 0.
-.15633 1.11082 0.326671      2.58823 2.50762 1.23707         -2.87097 2.29272 -0.0233288     1. 0. 0.
-.46453 0.112321 2.50927      2.58823 2.50762 1.23707         -2.87097 2.29272 -0.0233288     1. 0. 0.
-  -0.287266 -1.67078 2.48017    1.46453 0.112321 2.50927        2.06173 -1.39848 1.79787        1. 0. 0.
-.790434 -3.69418 1.29027     -0.287266 -1.67078 2.48017      2.06173 -1.39848 1.79787        1. 0. 0.
-.790434 -3.69418 1.29027     0.917634 -3.66448 -0.484829     -1.52417 -3.64138 0.503471      1. 0. 0.
-  -0.573766 -1.42458 -2.91753   0.450934 -2.72908 -2.23353      -1.62867 -3.74268 -1.79493      1. 0. 0.
-  -0.573766 -1.42458 -2.91753   -2.77417 -0.570279 -1.12083     -3.94777 1.63002 -1.27603       1. 0. 0.
-.03283 1.01972 -2.14533      -2.87097 2.29272 -0.0233288     -3.94777 1.63002 -1.27603       1. 0. 0.
-.917634 -3.66448 -0.484829   0.450934 -2.72908 -2.23353      -1.62867 -3.74268 -1.79493      1. 0. 0.
-  -0.759066 -0.265179 1.48487   1.46453 0.112321 2.50927        -2.87097 2.29272 -0.0233288     1. 0. 0.
-  -0.759066 -0.265179 1.48487   -0.287266 -1.67078 2.48017      1.46453 0.112321 2.50927        1. 0. 0.
-.790434 -3.69418 1.29027     -1.52417 -3.64138 0.503471      -0.287266 -1.67078 2.48017      1. 0. 0.
-.917634 -3.66448 -0.484829   -1.52417 -3.64138 0.503471      -1.62867 -3.74268 -1.79493      1. 0. 0.
-  -0.573766 -1.42458 -2.91753   -1.62867 -3.74268 -1.79493      -2.49667 -2.18078 -1.79993      1. 0. 0.
-  -0.573766 -1.42458 -2.91753   -2.49667 -2.18078 -1.79993      -2.77417 -0.570279 -1.12083     1. 0. 0.
-  -2.49667 -2.18078 -1.79993    -2.77417 -0.570279 -1.12083     -3.94777 1.63002 -1.27603       1. 0. 0.
-  -2.87097 2.29272 -0.0233288   -3.94777 1.63002 -1.27603       -4.04457 4.49292 -0.178529      1. 0. 0.
-  -0.759066 -0.265179 1.48487   -3.66127 0.565021 1.57007       -2.87097 2.29272 -0.0233288     1. 0. 0.
-  -0.759066 -0.265179 1.48487   -0.287266 -1.67078 2.48017      -3.66127 0.565021 1.57007       1. 0. 0.
-  -1.52417 -3.64138 0.503471    -0.287266 -1.67078 2.48017      -2.45927 -1.63678 1.72157       1. 0. 0.
-  -1.52417 -3.64138 0.503471    -1.62867 -3.74268 -1.79493      -3.90927 -3.49908 0.839771      1. 0. 0.
-  -1.62867 -3.74268 -1.79493    -2.49667 -2.18078 -1.79993      -4.17327 -2.53828 -0.635229     1. 0. 0.
-  -2.49667 -2.18078 -1.79993    -4.17327 -2.53828 -0.635229     -3.94777 1.63002 -1.27603       1. 0. 0.
-.58823 2.50762 1.23707       -2.87097 2.29272 -0.0233288     -4.04457 4.49292 -0.178529      1. 0. 0.
-  -3.94777 1.63002 -1.27603     -4.04457 4.49292 -0.178529      -5.12137 3.83022 -1.43123       1. 0. 0.
-  -3.66127 0.565021 1.57007     -2.87097 2.29272 -0.0233288     -4.83487 2.76522 1.41487        1. 0. 0.
-  -0.287266 -1.67078 2.48017    -2.45927 -1.63678 1.72157       -3.66127 0.565021 1.57007       1. 0. 0.
-  -1.52417 -3.64138 0.503471    -2.45927 -1.63678 1.72157       -3.90927 -3.49908 0.839771      1. 0. 0.
-  -1.62867 -3.74268 -1.79493    -4.17327 -2.53828 -0.635229     -3.90927 -3.49908 0.839771      1. 0. 0.
-  -4.17327 -2.53828 -0.635229   -3.94777 1.63002 -1.27603       -6.42617 1.74722 -0.982629      1. 0. 0.
-  -3.94777 1.63002 -1.27603     -5.12137 3.83022 -1.43123       -6.42617 1.74722 -0.982629      1. 0. 0.
   -5.12137 3.83022 -1.43123     -7.11497 2.49352 0.479071       -6.42617 1.74722 -0.982629      1. 0. 0.
-  -6.27887 0.926121 0.589671    -7.11497 2.49352 0.479071       -6.42617 1.74722 -0.982629      1. 0. 0.
-  -4.17327 -2.53828 -0.635229   -6.27887 0.926121 0.589671      -6.42617 1.74722 -0.982629      1. 0. 0.
-  -4.75167 -1.93408 0.935971    -4.17327 -2.53828 -0.635229     -6.27887 0.926121 0.589671      1. 0. 0.
-  -4.75167 -1.93408 0.935971    -3.66127 0.565021 1.57007       -6.27887 0.926121 0.589671      1. 0. 0.
-  -3.66127 0.565021 1.57007     -4.83487 2.76522 1.41487        -6.27887 0.926121 0.589671      1. 0. 0.
-  -4.83487 2.76522 1.41487      -6.27887 0.926121 0.589671      -7.11497 2.49352 0.479071       1. 0. 0.
-  -4.83487 2.76522 1.41487      -5.81017 4.57652 0.0304712      -7.11497 2.49352 0.479071       1. 0. 0.
   -5.81017 4.57652 0.0304712    -5.12137 3.83022 -1.43123       -7.11497 2.49352 0.479071       1. 0. 0.
   -5.81017 4.57652 0.0304712    -4.04457 4.49292 -0.178529      -5.12137 3.83022 -1.43123       1. 0. 0.
-  -4.83487 2.76522 1.41487      -5.81017 4.57652 0.0304712      -4.04457 4.49292 -0.178529      1. 0. 0.
-  -2.87097 2.29272 -0.0233288   -4.83487 2.76522 1.41487        -4.04457 4.49292 -0.178529      1. 0. 0.
-.58823 2.50762 1.23707       -2.87097 2.29272 -0.0233288     -4.04457 4.49292 -0.178529      1. 0. 0.
-  -2.45927 -1.63678 1.72157     -4.75167 -1.93408 0.935971      -3.66127 0.565021 1.57007       1. 0. 0.
-  -2.45927 -1.63678 1.72157     -4.75167 -1.93408 0.935971      -3.90927 -3.49908 0.839771      1. 0. 0.
-  -4.75167 -1.93408 0.935971    -4.17327 -2.53828 -0.635229     -3.90927 -3.49908 0.839771      1. 0. 0.
 #  terminating special property
 …
 .0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0
   -4.27807 -2.65715 0.380171    5       1 0 0
+  -4.99203 4.29989 -0.526429    5       1 0 0
   terminating special property

molecuilder/tests/testsuite.at

-              r0af58f
+              r551a58
 AT_CHECK([pwd],[ignore],[ignore])
 AT_CHECK([../../molecuilder -v], 0, [stdout], [ignore])
 AT_CHECK([fgrep molecuilder stdout], 0, [ignore], [ignore])
+AT_CHECK([fgrep olecuilder stdout], 0, [ignore], [ignore])
 AT_CHECK([../../molecuilder -h], 0, [stdout], [ignore])
 AT_CHECK([fgrep "Give this help screen" stdout], 0, [ignore], [ignore])
 AT_CHECK([../../molecuilder -e], 0, [ignore], [stderr])
+AT_CHECK([../../molecuilder -e], 255, [ignore], [stderr])
 AT_CHECK([fgrep "Not enough or invalid arguments" stderr], 0, [ignore], [ignore])
 AT_CHECK([../../molecuilder test.conf], 0, [stdout], [stderr])
 …
 AT_SETUP([Simple configuration - invalid commands on present configs])
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/7/pre/test.conf .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -t -s -b -E -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "Not enough or invalid" stderr], 0, [9
+], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -t], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -s -b -E -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -b -E -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -E -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "CRITICAL: Not enough or invalid" stderr], 0, [ignore], [ignore])
 AT_CLEANUP
 …
 AT_SETUP([Graph - DFS analysis])
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Graph/1/pre/test.conf .], 0)
 AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -D 2.], 0, [stdout], [stderr])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -vvv -D 2.], 0, [stdout], [stderr])
 AT_CHECK([fgrep -c "No rings were detected in the molecular structure." stdout], 0, [1
 ], [ignore])
 …
 AT_CHECK([file=ConvexEnvelope.dat; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/2/post/$file], 0, [ignore], [ignore])
 AT_CHECK([file=ConvexEnvelope.r3d; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/2/post/$file], 0, [ignore], [ignore])
 AT_CHECK([fgrep "RESULT: The summed volume is 16.401577 angstrom^3" stdout], 0, [ignore], [ignore])
+AT_CHECK([fgrep "tesselated volume area is 16.4016 angstrom^3" stdout], 0, [ignore], [ignore])
 AT_CHECK([diff ConvexEnvelope.dat NonConvexEnvelope.dat], 0, [ignore], [ignore])
 AT_CLEANUP
 …
 #AT_CHECK([file=ConvexEnvelope.dat; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/4/post/$file], 0, [ignore], [ignore])
 #AT_CHECK([file=ConvexEnvelope.r3d; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/4/post/$file], 0, [ignore], [ignore])
 #AT_CHECK([fgrep "RESULT: The summed volume is 16.401577 angstrom^3" stdout], 0, [ignore], [ignore])
+#AT_CHECK([fgrep "tesselated volume area is 16.4016 angstrom^3" stdout], 0, [ignore], [ignore])
 #AT_CLEANUP

pcp/src/Makefile.am

r0af58f	r551a58
7	7	FORCE:
8	8	$(srcdir)/.git-version: FORCE
9		@if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe ~~--tags~~ HEAD) > .git-version-t 2>/dev/null \
	9	@if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe HEAD) > .git-version-t 2>/dev/null \
10	10	&& ! diff .git-version-t $(srcdir)/.git-version >/dev/null 2>&1; then \
11	11	mv -f .git-version-t $(srcdir)/.git-version; \

util/src/Makefile.am

-              r0af58f
+              r551a58
         LinearlyInterpolateBetweenPCPConfigs.py \
         OutputRaster3DHeader.py \
-        LinearlyInterpolateBetweenPCPConfigs.py \
         PairCorrelationTIP3toResiduum.py \
         PairCorrelationToCluster.py \
 …
         ReMapDBondFileFromXYZ.py \
         ReMapDBondFile.py \
-        RemoveConvexPart.py \
-        RemoveRandomAtoms.py \
         RemoveConvexPart.py \
         RemoveRandomAtoms.py \
 …
 FORCE:
 $(srcdir)/.git-version: FORCE
         @if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe --tags HEAD) > .git-version-t 2>/dev/null \
+        @if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe HEAD) > .git-version-t 2>/dev/null \
           && ! diff .git-version-t $(srcdir)/.git-version >/dev/null 2>&1; then \
           mv -f .git-version-t $(srcdir)/.git-version; \

Context Navigation

Changes in / [0af58f:551a58]

Legend:

Download in other formats: