Diff [8d6d31e159e6d6b40179835f02717acb3398e53f:8f4df1ef77448922e1b977a3e138bcc2a4a4dad0] for / – MoleCuilder

src/Actions/AnalysisAction/MolecularVolumeAction.cpp

r8d6d31	r8f4df1
60	60	LCList = new LinkedCell(mol, 10.);
61	61	config * const configuration = World::getInstance().getConfig();
62		~~Boundaries *BoundaryPoints = NULL;~~
63	62	//FindConvexBorder(mol, BoundaryPoints, TesselStruct, LCList, argv[argptr]);
64	63	FindNonConvexBorder(mol, TesselStruct, LCList, 5., NULL);

src/Actions/AnalysisAction/PairCorrelationAction.cpp

r8d6d31	r8f4df1
70	70	double BinStart = 0.;
71	71	double BinWidth = 0.;
72		~~molecule *Boundary = NULL;~~
73	72	string outputname;
74	73	string binoutputname;

src/Actions/AnalysisAction/PointCorrelationAction.cpp

r8d6d31	r8f4df1
71	71	double BinStart = 0.;
72	72	double BinWidth = 0.;
73		~~molecule *Boundary = NULL;~~
74	73	string outputname;
75	74	string binoutputname;

src/Actions/AnalysisAction/PrincipalAxisSystemAction.cpp

-              r8d6d31
+              r8f4df1
     // sum up inertia tensor
     for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       Vector x = (*iter)->x;
+      Vector x = (*iter)->getPosition();
       x -= *CenterOfGravity;
+      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
+      const double mass = (*iter)->getType()->mass;
+      InertiaTensor.at(0,0) += mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += mass*(x[0]*x[0] + x[1]*x[1]);
+    }
     // print InertiaTensor for debugging

src/Actions/AnalysisAction/SurfaceCorrelationAction.cpp

r8d6d31	r8f4df1
111	111	// find biggest molecule
112	112	std::vector<molecule*> molecules = World::getInstance().getSelectedMolecules();
113		~~int counter = molecules.size();~~
114	113	LCList = new LinkedCell(Boundary, LCWidth);
115	114	FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);

src/Actions/AtomAction/AddAction.cpp

-              r8d6d31
+              r8f4df1
   // execute action
   atom * first = World::getInstance().createAtom();
   first->type = elemental;
   first->x = position;
   DoLog(1) && (Log() << Verbose(1) << "Adding new atom with element " << first->type->name << " at " << (first->x) << "." << endl);
+  first->setType(elemental);
+  first->setPosition(position);
+  DoLog(1) && (Log() << Verbose(1) << "Adding new atom with element " << first->getType()->name << " at " << (first->getPosition()) << "." << endl);
   // TODO: remove when all of World's atoms are stored.
   std::vector<molecule *> molecules = World::getInstance().getAllMolecules();

src/Actions/AtomAction/ChangeElementAction.cpp

r8d6d31	r8f4df1
61	61	mol = first->getMolecule();
62	62	first->removeFromMolecule(); // remove atom
63		first->~~type = elemental~~;
	63	first->setType(elemental);
64	64	mol->AddAtom(first); // add atom to ensure correctness of formula
65	65	}

src/Actions/FragmentationAction/SubgraphDissectionAction.cpp

-              r8d6d31
+              r8f4df1
 #include "Actions/FragmentationAction/SubgraphDissectionAction.hpp"
 #include "Actions/ActionRegistry.hpp"
+#include "Descriptors/MoleculeDescriptor.hpp"
 #include "atom.hpp"
+#include "bond.hpp"
+#include "bondgraph.hpp"
 #include "config.hpp"
 #include "Helpers/Log.hpp"
 #include "molecule.hpp"
-#include "Descriptors/MoleculeDescriptor.hpp"
 #include "stackclass.hpp"
 #include "World.hpp"
 …
   // @TODO rather do the dissection afterwards
   MoleculeListClass *molecules = World::getInstance().getMolecules();
+  molecules->DissectMoleculeIntoConnectedSubgraphs(World::getInstance().getPeriode(), World::getInstance().getConfig());
+  config * const configuration = World::getInstance().getConfig();
+  // 0a. remove all present molecules
+  vector<molecule *> allmolecules = World::getInstance().getAllMolecules();
+  for (vector<molecule *>::iterator MolRunner = allmolecules.begin(); MolRunner != allmolecules.end(); ++MolRunner) {
+    molecules->erase(*MolRunner);
+    World::getInstance().destroyMolecule(*MolRunner);
+  }
+  // 0b. remove all bonds and construct a molecule with all atoms
+  molecule *mol = World::getInstance().createMolecule();
+  vector <atom *> allatoms = World::getInstance().getAllAtoms();
+  for(vector<atom *>::iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) {
+    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
+      delete(*BondRunner);
+    mol->AddAtom(*AtomRunner);
+  }
+  // 1. create the bond structure of the single molecule
+  if (configuration->BG != NULL) {
+    if (!configuration->BG->ConstructBondGraph(mol)) {
+      World::getInstance().destroyMolecule(mol);
+      DoeLog(1) && (eLog()<< Verbose(1) << "There are no bonds." << endl);
+      return Action::failure;
+    }
+  } else {
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no BondGraph class present to create bonds." << endl);
+    return Action::failure;
+  }
+  // 2. scan for connected subgraphs
+  MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
+  class StackClass<bond *> *BackEdgeStack = NULL;
+  Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
+  delete(BackEdgeStack);
+  if ((Subgraphs == NULL) || (Subgraphs->next == NULL)) {
+    World::getInstance().destroyMolecule(mol);
+    DoeLog(1) && (eLog()<< Verbose(1) << "There are no atoms." << endl);
+    return Action::failure;
+  }
+  // 3. dissect (the following construct is needed to have the atoms not in the order of the DFS, but in
+  // the original one as parsed in)
+  // TODO: Optimize this, when molecules just contain pointer list of global atoms!
+  // 4a. create array of molecules to fill
+  const int MolCount = Subgraphs->next->Count();
+  char number[MAXSTRINGSIZE];
+  molecule **moleculelist = new molecule *[MolCount];
+  MoleculeLeafClass *MolecularWalker = Subgraphs;
+  for (int i=0;i<MolCount;i++) {
+    MolecularWalker = MolecularWalker->next;
+    moleculelist[i] = World::getInstance().createMolecule();
+    moleculelist[i]->ActiveFlag = true;
+    strncpy(moleculelist[i]->name, mol->name, MAXSTRINGSIZE);
+    if (MolCount > 1) {
+      sprintf(number, "-%d", i+1);
+      strncat(moleculelist[i]->name, number, MAXSTRINGSIZE - strlen(mol->name) - 1);
+    }
+    DoLog(1) && (Log() << Verbose(1) << "MolName is " << moleculelist[i]->name << ", id is " << moleculelist[i]->getId() << endl);
+    for (molecule::iterator iter = MolecularWalker->Leaf->begin(); iter != MolecularWalker->Leaf->end(); ++iter) {
+      DoLog(1) && (Log() << Verbose(1) << **iter << endl);
+    }
+    molecules->insert(moleculelist[i]);
+  }
+  // 4b. create and fill map of which atom is associated to which connected molecule (note, counting starts at 1)
+  int FragmentCounter = 0;
+  map<int, atom *> AtomToFragmentMap;
+  MolecularWalker = Subgraphs;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    for (molecule::iterator iter = MolecularWalker->Leaf->begin(); !MolecularWalker->Leaf->empty(); iter = MolecularWalker->Leaf->begin()) {
+      atom * Walker = *iter;
+      DoLog(1) && (Log() << Verbose(1) << "Re-linking " << Walker << "..." << endl);
+      MolecularWalker->Leaf->erase(iter);
+      moleculelist[FragmentCounter]->AddAtom(Walker);    // counting starts at 1
+    }
+    FragmentCounter++;
+  }
+  // TODO: When DepthFirstSearchAnalysis does not use AddCopyAtom() anymore, we don't need to delete all original atoms.
+  // 4d. destroy the original molecule
+  for (molecule::iterator AtomRunner = mol->begin(); !mol->empty(); AtomRunner = mol->begin())
+    World::getInstance().destroyAtom(*AtomRunner);
+  World::getInstance().destroyMolecule(mol);
+  // 4d. we don't need to redo bonds, as they are connected subgraphs and still maintain their ListOfBonds, but we have to remove them from first..last list
+  // TODO: check whether this is really not needed anymore
+  // 4e. free Leafs
+  MolecularWalker = Subgraphs;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    delete(MolecularWalker->previous);
+  }
+  delete(MolecularWalker);
+  delete[](moleculelist);
+  DoLog(1) && (Log() << Verbose(1) << "I scanned " << FragmentCounter << " molecules." << endl);
   return Action::success;
+}

src/Actions/MoleculeAction/RotateAroundOriginByAngleAction.cpp

r8d6d31	r8f4df1
91	91
92	92	for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
93		((iter)->node) = RotationAxis.rotateVector(((iter)->node), alpha);
	93	(iter)->setPosition(RotationAxis.rotateVector((iter)->getPosition(), alpha));
94	94	}
95	95	}
…	…
109	109
110	110	for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
111		((iter)->node) = RotationAxis.rotateVector(((iter)->node), -state->alpha);
	111	(iter)->setPosition(RotationAxis.rotateVector((iter)->getPosition(), -state->alpha));
112	112	}
113	113	}
…	…
127	127
128	128	for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
129		((iter)->node) = RotationAxis.rotateVector(((iter)->node), state->alpha);
	129	(iter)->setPosition(RotationAxis.rotateVector((iter)->getPosition(), state->alpha));
130	130	}
131	131	}

src/Actions/MoleculeAction/RotateAroundSelfByAngleAction.cpp

-              r8d6d31
+              r8f4df1
   for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     *((*iter)->node) = RotationAxis.rotateVector(*((*iter)->node), alpha);
+    (*iter)->setPosition(RotationAxis.rotateVector((*iter)->getPosition(), alpha));
+  }
   DoLog(0) && (Log() << Verbose(0) << "done." << endl);
 …
   for (molecule::iterator iter = state->mol->begin(); iter != state->mol->end(); ++iter) {
     *((*iter)->node) = RotationAxis.rotateVector(*((*iter)->node), -state->alpha);
+    (*iter)->setPosition(RotationAxis.rotateVector((*iter)->getPosition(), -state->alpha));
+  }
 …
   for (molecule::iterator iter = state->mol->begin(); iter != state->mol->end(); ++iter) {
     *((*iter)->node) = RotationAxis.rotateVector(*((*iter)->node), state->alpha);
+    (*iter)->setPosition(RotationAxis.rotateVector((*iter)->getPosition(), state->alpha));
+  }

src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp

-              r8d6d31
+              r8f4df1
     // sum up inertia tensor
     for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       Vector x = (*iter)->x;
+      Vector x = (*iter)->getPosition();
       x -= *CenterOfGravity;
+      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
+      const double mass = (*iter)->getType()->mass;
+      InertiaTensor.at(0,0) += mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += mass*(x[0]*x[0] + x[1]*x[1]);
+    }
     // print InertiaTensor for debugging
 …
     for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       *((*iter)->node) -= *CenterOfGravity;
       *((*iter)->node) = RotationAxis.rotateVector(*((*iter)->node), alpha);
       *((*iter)->node) += *CenterOfGravity;
+      *(*iter) -= *CenterOfGravity;
+      (*iter)->setPosition(RotationAxis.rotateVector((*iter)->getPosition(), alpha));
+      *(*iter) += *CenterOfGravity;
+    }
     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
 …
     // sum up inertia tensor
     for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       Vector x = (*iter)->x;
+      Vector x = (*iter)->getPosition();
       x -= *CenterOfGravity;
+      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
+      const double mass = (*iter)->getType()->mass;
+      InertiaTensor.at(0,0) += mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += mass*(x[0]*x[0] + x[1]*x[1]);
       // print InertiaTensor for debugging
       DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << InertiaTensor << endl);

src/Actions/TesselationAction/ConvexEnvelopeAction.cpp

r8d6d31	r8f4df1
84	84	DoLog(1) && (Log() << Verbose(1) << "Storing tecplot non-convex data in " << filenameNonConvex << "." << endl);
85	85	LCList = new LinkedCell(mol, 100.);
86		~~Boundaries *BoundaryPoints = NULL;~~
87	86	//FindConvexBorder(mol, BoundaryPoints, TesselStruct, LCList, argv[argptr]);
88	87	// TODO: Beide Funktionen sollten streams anstelle des Filenamen benutzen, besser fuer unit tests

src/Actions/WorldAction/AddEmptyBoundaryAction.cpp

-              r8d6d31
+              r8f4df1
   ASSERT(AllAtoms.size() > 0, "There must be atoms present for AddingEmptyBoundary.");
   vector<atom *>::iterator AtomRunner = AllAtoms.begin();
   Min = (*AtomRunner)->x;
   Max = (*AtomRunner)->x;
+  Min = (*AtomRunner)->getPosition();
+  Max = (*AtomRunner)->getPosition();
   for (; AtomRunner != AllAtoms.end(); ++AtomRunner) {
     for (int i=0;i<NDIM;i++) {
       if ((*AtomRunner)->x[i] > Max[i])
         Max[i] = (*AtomRunner)->x[i];
       if ((*AtomRunner)->x[i] < Min[i])
         Min[i] = (*AtomRunner)->x[i];
+      if ((*AtomRunner)->at(i) > Max[i])
+        Max[i] = (*AtomRunner)->at(i);
+      if ((*AtomRunner)->at(i) < Min[i])
+        Min[i] = (*AtomRunner)->at(i);
+    }
+  }
 …
   AtomRunner = AllAtoms.begin();
   for (; AtomRunner != AllAtoms.end(); ++AtomRunner)
     (*AtomRunner)->x -= Min - boundary;
+    *(*AtomRunner) -= Min - boundary;
   return Action::success;
+}

src/Actions/WorldAction/CenterOnEdgeAction.cpp

-              r8d6d31
+              r8f4df1
   ASSERT(AllAtoms.size() > 0, "For CenteronEdge atoms must be present.");
   vector<atom *>::iterator AtomRunner = AllAtoms.begin();
   Min = (*AtomRunner)->x;
   Max = (*AtomRunner)->x;
+  Min = (*AtomRunner)->getPosition();
+  Max = (*AtomRunner)->getPosition();
   for (; AtomRunner != AllAtoms.end(); ++AtomRunner) {
     for (int i=0;i<NDIM;i++) {
       if ((*AtomRunner)->x[i] > Max[i])
         Max[i] = (*AtomRunner)->x[i];
       if ((*AtomRunner)->x[i] < Min[i])
         Min[i] = (*AtomRunner)->x[i];
+      if ((*AtomRunner)->at(i) > Max[i])
+        Max[i] = (*AtomRunner)->at(i);
+      if ((*AtomRunner)->at(i) < Min[i])
+        Min[i] = (*AtomRunner)->at(i);
+    }
+  }
 …
   // translate all atoms, such that Min is aty (0,0,0)
   for (vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner)
     (*AtomRunner)->x -= Min;
+    *(*AtomRunner) -= Min;
   return Action::success;

src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp

-              r8d6d31
+              r8f4df1
   vector<molecule *> molecules = World::getInstance().getAllMolecules();
   for (vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner) {
+    if (point.DistanceSquared((*AtomRunner)->x) > radius*radius) { // distance to first above radius ...
+    if ((*AtomRunner)->DistanceSquared(point) > radius*radius) { // distance to first above radius ...
+//      cout << "Removing " << (*AtomRunner)->getType()->symbol << (*AtomRunner)->getId() << " at " << (*AtomRunner)->getPosition() << " as distance is " << sqrt((*AtomRunner)->DistanceSquared(point)) << endl;
       // TODO: This is not necessary anymore when atoms are completely handled by World (create/destroy and load/save)
       for (vector<molecule *>::iterator iter = molecules.begin();iter != molecules.end();++iter)
 …
       World::getInstance().destroyAtom(*AtomRunner);
+    }
+//    else {
+//      cout << "Keeping" << (*AtomRunner)->getType()->symbol << (*AtomRunner)->getId() << " at " << (*AtomRunner)->getPosition() << " as distance is " << sqrt((*AtomRunner)->DistanceSquared(point)) << endl;
+//    }
+  }
   return Action::success;

src/Actions/WorldAction/RepeatBoxAction.cpp

-              r8d6d31
+              r8f4df1
 #include <iostream>
 #include <string>
+#include <vector>
 using namespace std;
 …
   molecule *newmol = NULL;
   Vector ** vectors = NULL;
+  std::vector<Vector> vectors;
   for (n[0] = 0; n[0] < Repeater[0]; n[0]++) {
     y[0] = n[0];
 …
           if (count != 0) {  // if there is more than none
             Elements = new const element *[count];
             vectors = new Vector *[count];
+            vectors.resize(count);
             j = 0;
             for(molecule::iterator AtomRunner = mol->begin(); AtomRunner != mol->end(); ++AtomRunner) {
               Elements[j] = (*AtomRunner)->type;
               vectors[j] = &(*AtomRunner)->x;
+              Elements[j] = (*AtomRunner)->getType();
+              vectors[j] = (*AtomRunner)->getPosition();
               j++;
+            }
 …
             for (int k=count;k--;) { // go through every atom of the original cell
               Walker = World::getInstance().createAtom(); // create a new body
               Walker->x = (*vectors[k]) + x;
               Walker->type = Elements[k];  // insert original element
+              Walker->setPosition((vectors[k]) + x);
+              Walker->setType(Elements[k]);  // insert original element
               cout << "new atom is " << *Walker << endl;
               newmol->AddAtom(Walker);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
 …
             // free memory
             delete[](Elements);
-            delete[](vectors);
           } else {
             DoLog(1) && (Log() << Verbose(1) << "\t ... is empty." << endl);

src/Actions/WorldAction/ScaleBoxAction.cpp

r8d6d31	r8f4df1
60	60	vector<atom*> AllAtoms = World::getInstance().getAllAtoms();
61	61	for(vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner) {
62		(*AtomRunner)->x.ScaleAll(x);
	62	(*AtomRunner)->ScaleAll(x);
63	63	}
64	64

src/AtomSet.hpp

-              r8d6d31
+              r8f4df1
     workOnNodePointer(Function &_f) : f(_f){}
     void operator()(atom *atom){
+      if(atom->node)
+        *atom->node = f(*atom->node);
+      atom->setPosition(f(atom->getPosition()));
+    }
     Function &f;
 …
 inline void AtomSetMixin<Set>::translate(const Vector &translater){
   BOOST_FOREACH(atom *atom,*this){
+    if(atom->node){
+      *(atom->node) += translater;
+    }
+    *(atom) += translater;
+  }
+}

src/Legacy/oldmenu.cpp

-              r8d6d31
+              r8f4df1
         first = World::getInstance().createAtom();
         std::vector<element *> elements;
         dialog->queryVector("Please enter coordinates: ",&first->x,World::getInstance().getDomain(), false);
+        dialog->queryVector("Please enter coordinates: ",&first->AtomicPosition,World::getInstance().getDomain(), false);
         dialog->queryElement("Please choose element: ",&elements);
         if(dialog->display()){
           if (elements.size() == 1) {
             first->type = elements.at(0);
+            first->setType(elements.at(0));
             mol->AddAtom(first);  // add to molecule
           } else {
 …
           auto_ptr<Dialog> dialog(UIFactory::getInstance().makeDialog());
           dialog->queryVector("Enter reference coordinates.",&x,World::getInstance().getDomain(), true);
           dialog->queryVector("Enter relative coordinates.",&first->x,World::getInstance().getDomain(), false);
+          dialog->queryVector("Enter relative coordinates.",&first->AtomicPosition,World::getInstance().getDomain(), false);
           if((aborted = !dialog->display())){
             continue;
+          }
           first->x += x;
         } while (!aborted && !(valid = mol->CheckBounds((const Vector *)&first->x)));
+          first->AtomicPosition += x;
+        } while (!aborted && !(valid = mol->CheckBounds((const Vector *)&first->AtomicPosition)));
         if(!aborted){
           first->type = periode->AskElement();  // give type
+          first->AtomicElement = periode->AskElement();  // give type
           mol->AddAtom(first);  // add to molecule
+        }
 …
           auto_ptr<Dialog> dialog(UIFactory::getInstance().makeDialog());
           second = mol->AskAtom("Enter atom number: ");
           dialog->queryVector("Enter relative coordinates.",&first->x,World::getInstance().getDomain(), false);
+          dialog->queryVector("Enter relative coordinates.",&first->AtomicPosition,World::getInstance().getDomain(), false);
           dialog->display();
           for (int i=NDIM;i--;) {
             first->x[i] += second->x[i];
+            first->AtomicPosition[i] += second->AtomicPosition[i];
+          }
         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
         first->type = periode->AskElement();  // give type
+        } while (!(valid = mol->CheckBounds((const Vector *)&first->AtomicPosition)));
+        first->AtomicElement = periode->AskElement();  // give type
         mol->AddAtom(first);  // add to molecule
+      }
 …
         do {
           if (!valid) {
             eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl;
+            eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->AtomicPosition << endl;
+          }
           Log() << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
 …
           */
           // calc axis vector
           x= second->x - third->x;
+          x= second->AtomicPosition - third->AtomicPosition;
           x.Normalize();
           Log() << Verbose(0) << "x: " << x << endl;
           z = Plane(second->x,third->x,fourth->x).getNormal();
+          z = Plane(second->AtomicPosition,third->AtomicPosition,fourth->AtomicPosition).getNormal();
           Log() << Verbose(0) << "z: " << z << endl;
           y = Plane(x,z,0).getNormal();
 …
           // rotate vector around first angle
           first->x = x;
           first->x = Line(zeroVec,z).rotateVector(first->x,b - M_PI);
           Log() << Verbose(0) << "Rotated vector: " << first->x << endl,
+          first->AtomicPosition = x;
+          first->AtomicPosition = Line(zeroVec,z).rotateVector(first->AtomicPosition,b - M_PI);
+          Log() << Verbose(0) << "Rotated vector: " << first->AtomicPosition << endl,
           // remove the projection onto the rotation plane of the second angle
           n = y;
           n.Scale(first->x.ScalarProduct(y));
+          n.Scale(first->AtomicPosition.ScalarProduct(y));
           Log() << Verbose(0) << "N1: " << n << endl;
           first->x -= n;
           Log() << Verbose(0) << "Subtracted vector: " << first->x << endl;
+          first->AtomicPosition -= n;
+          Log() << Verbose(0) << "Subtracted vector: " << first->AtomicPosition << endl;
           n = z;
           n.Scale(first->x.ScalarProduct(z));
+          n.Scale(first->AtomicPosition.ScalarProduct(z));
           Log() << Verbose(0) << "N2: " << n << endl;
           first->x -= n;
           Log() << Verbose(0) << "2nd subtracted vector: " << first->x << endl;
+          first->AtomicPosition -= n;
+          Log() << Verbose(0) << "2nd subtracted vector: " << first->AtomicPosition << endl;
           // rotate another vector around second angle
 …
           // add the two linear independent vectors
           first->x += n;
           first->x.Normalize();
           first->x.Scale(a);
           first->x += second->x;
           Log() << Verbose(0) << "resulting coordinates: " << first->x << endl;
         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
         first->type = periode->AskElement();  // give type
+          first->AtomicPosition += n;
+          first->AtomicPosition.Normalize();
+          first->AtomicPosition.Scale(a);
+          first->AtomicPosition += second->AtomicPosition;
+          Log() << Verbose(0) << "resulting coordinates: " << first->AtomicPosition << endl;
+        } while (!(valid = mol->CheckBounds((const Vector *)&first->AtomicPosition)));
+        first->AtomicElement = periode->AskElement();  // give type
         mol->AddAtom(first);  // add to molecule
         break;
 …
           if (j != -1) {
             second = mol->FindAtom(j);
             atoms[i++] = &(second->x);
+            atoms[i++] = &(second->AtomicPosition);
+          }
         } while ((j != -1) && (i<128));
         if (i >= 2) {
           LSQdistance(first->x,(const Vector **)atoms, i);
           Log() << Verbose(0) << first->x;
           first->type = periode->AskElement();  // give type
+          LSQdistance(first->AtomicPosition,(const Vector **)atoms, i);
+          Log() << Verbose(0) << first->AtomicPosition;
+          first->AtomicElement = periode->AskElement();  // give type
           mol->AddAtom(first);  // add to molecule
         } else {
 …
       third = mol->AskAtom("Enter third atom: ");
       n = Plane(first->x,second->x,third->x).getNormal();
+      n = Plane(first->AtomicPosition,second->AtomicPosition,third->AtomicPosition).getNormal();
       break;
     case 'b': // normal vector of mirror plane
 …
       second = mol->AskAtom("Enter second atom: ");
       n = first->x - second->x;
+      n = first->AtomicPosition - second->AtomicPosition;
       n.Normalize();
       break;
 …
       third = mol->AskAtom("Enter third atom: ");
       n = Plane(first->x,second->x,third->x).getNormal();
+      n = Plane(first->AtomicPosition,second->AtomicPosition,third->AtomicPosition).getNormal();
       break;
     case 'b': // normal vector of mirror plane
 …
       second = mol->AskAtom("Enter second atom: ");
       n = first->x - second->x;
+      n = first->AtomicPosition - second->AtomicPosition;
       n.Normalize();
       break;
 …
         for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
           runner = iter++;
           if ((*runner)->x.DistanceSquared((*runner)->x) > tmp1*tmp1) // distance to first above radius ...
+          if ((*runner)->AtomicPosition.DistanceSquared((*runner)->AtomicPosition) > tmp1*tmp1) // distance to first above radius ...
             mol->RemoveAtom((*runner));
+        }
 …
       for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
         runner = iter++;
         if (((*runner)->x[axis] < tmp1) || ((*runner)->x[axis] > tmp2)) {// out of boundary ...
+        if (((*runner)->AtomicPosition[axis] < tmp1) || ((*runner)->AtomicPosition[axis] > tmp2)) {// out of boundary ...
           //Log() << Verbose(0) << "Atom " << *(*runner) << " with " << (*runner)->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
           mol->RemoveAtom((*runner));
 …
       for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
         Z = (*iter)->type->Z;
+        Z = (*iter)->getType()->Z;
         tmp1 = 0.;
         if (first != (*iter)) {
           x = first->x - (*iter)->x;
+          x = first->AtomicPosition - (*iter)->AtomicPosition;
           tmp1 = x.Norm();
+        }
 …
+      }
       for (int i=MAX_ELEMENTS;i--;)
         if (min[i] != 0.) Log() << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
+        if (min[i] != 0.) Log() << Verbose(0) << "Minimum Bond length between " << first->getType()->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
       break;
 …
       for (int i=NDIM;i--;)
         min[i] = 0.;
       x = first->x - second->x;
+      x = first->AtomicPosition - second->AtomicPosition;
       tmp1 = x.Norm();
       Log() << Verbose(1) << "Distance vector is " << x << "." << "/n"
 …
       third  = mol->AskAtom("Enter last atom: ");
       tmp1 = tmp2 = tmp3 = 0.;
       x = first->x - second->x;
       y = third->x - second->x;
+      x = first->AtomicPosition - second->AtomicPosition;
+      y = third->AtomicPosition - second->AtomicPosition;
       Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
       Log() << Verbose(0) << (acos(x.ScalarProduct(y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
 …
         minBond = 0.;
         for (int i=NDIM;i--;)
           minBond += (first->x[i]-second->x[i])*(first->x[i] - second->x[i]);
+          minBond += (first->AtomicPosition[i]-second->AtomicPosition[i])*(first->AtomicPosition[i] - second->AtomicPosition[i]);
         minBond = sqrt(minBond);
         Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;
+        Log() << Verbose(0) << "Current Bond length between " << first->getType()->name << " Atom " << first->nr << " and " << second->getType()->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;
         Log() << Verbose(0) << "Enter new bond length [a.u.]: ";
         cin >> bond;
         for (int i=NDIM;i--;) {
           second->x[i] -= (second->x[i]-first->x[i])/minBond*(minBond-bond);
+          second->AtomicPosition[i] -= (second->AtomicPosition[i]-first->AtomicPosition[i])/minBond*(minBond-bond);
+        }
         //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
 …
         cin >> Z;
         first->setType(Z);
         Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
+        Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->getType()->name << "." << endl;
+      }
       break;
 …
       j = 0;
       for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
         Elements[j] = (*iter)->type;
         vectors[j] = &(*iter)->x;
+        Elements[j] = (*iter)->AtomicElement;
+        vectors[j] = &(*iter)->AtomicPosition;
         j++;
+      }
 …
         for (int k=count;k--;) { // go through every atom of the original cell
           first = World::getInstance().createAtom(); // create a new body
           first->x = (*vectors[k]) + x;  // use coordinate of original atom
           first->type = Elements[k];  // insert original element
+          first->AtomicPosition = (*vectors[k]) + x;  // use coordinate of original atom
+          first->AtomicElement = Elements[k];  // insert original element
           mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
+        }

src/LinearAlgebra/Makefile.am

-              r8d6d31
+              r8f4df1
   Plane.cpp \
   Space.cpp \
+  Vector.cpp
+  vector_ops.cpp \
+  Vector.cpp \
+  VectorInterface.cpp
 LINALGHEADER = \
 …
   Plane.hpp \
   Space.hpp \
+  Vector.hpp
+  vector_ops.hpp \
+  Vector.hpp \
+  VectorInterface.hpp

src/LinearAlgebra/Vector.cpp

-              r8d6d31
+              r8f4df1
   gsl_vector_set(content->content,1,x2);
   gsl_vector_set(content->content,2,x3);
+};
+/** Constructor of class vector.
+ */
+Vector::Vector(const double x[3])
+{
+  content = new VectorContent();
+  gsl_vector_set(content->content,0,x[0]);
+  gsl_vector_set(content->content,1,x[1]);
+  gsl_vector_set(content->content,2,x[2]);
 };

src/LinearAlgebra/Vector.hpp

r8d6d31	r8f4df1
36	36	Vector();
37	37	Vector(const double x1, const double x2, const double x3);
	38	Vector(const double x[3]);
38	39	Vector(const Vector& src);
39	40	virtual ~Vector();

src/Makefile.am

-              r8d6d31
+              r8f4df1
 QTUI_DEFS =
+TESSELATIONSOURCE = \
+  BoundaryLineSet.cpp \
+  BoundaryPointSet.cpp \
+  BoundaryPolygonSet.cpp \
+  BoundaryTriangleSet.cpp \
+  CandidateForTesselation.cpp \
+  PointCloud.cpp \
+  tesselation.cpp \
+  tesselationhelpers.cpp \
+  TesselPoint.cpp
+TESSELATIONHEADER = \
+  BoundaryLineSet.hpp \
+  BoundaryPointSet.hpp \
+  BoundaryPolygonSet.hpp \
+  BoundaryTriangleSet.hpp \
+  CandidateForTesselation.hpp \
+  PointCloud.hpp \
+  tesselation.hpp \
+  tesselationhelpers.hpp \
+  TesselPoint.hpp
 MOLECUILDERSOURCE = \
   ${ANALYSISSOURCE} \
 …
   ${ATOMSOURCE} \
   ${PATTERNSOURCE} \
-  ${PARSERSOURCE} \
   ${SHAPESOURCE} \
   ${DESCRIPTORSOURCE} \
+  ${TESSELATIONSOURCE} \
   bond.cpp \
   bondgraph.cpp \
 …
   parser.cpp \
   periodentafel.cpp \
-  tesselation.cpp \
-  tesselationhelpers.cpp \
   ThermoStatContainer.cpp \
   triangleintersectionlist.cpp \
   UIElements/UIFactory.cpp \
-  vector_ops.cpp \
   World.cpp
 …
   ${ACTIONSHEADER} \
   ${ATOMHEADER} \
-  ${PARSERHEADER} \
   ${PATTERNHEADER} \
   ${SHAPEHEADER} \
   ${DESCRIPTORHEADER} \
+  ${TESSELATIONHEADER} \
   bond.hpp \
   bondgraph.hpp \
 …
   periodentafel.hpp \
   stackclass.hpp \
-  tesselation.hpp \
-  tesselationhelpers.hpp \
   ThermoStatContainer.hpp \
   triangleintersectionlist.hpp \
   UIElements/UIFactory.hpp \
-  vector_ops.hpp \
   World.hpp

src/Parser/MpqcParser.cpp

-              r8d6d31
+              r8f4df1
   // calculate center
   for (vector<atom *>::iterator runner = allatoms.begin();runner != allatoms.end(); ++runner)
     center += (*runner)->x;
+    center += (*runner)->getPosition();
   center.Scale(1./allatoms.size());
 …
   // calculate center
   for (vector<atom *>::iterator runner = allatoms.begin();runner != allatoms.end(); ++runner)
     center += (*runner)->x;
+    center += (*runner)->getPosition();
   center.Scale(1./allatoms.size());

src/Parser/PcpParser.cpp

-              r8d6d31
+              r8f4df1
   PlaneWaveSpecifics.MaxPsiDouble = PlaneWaveSpecifics.PsiMaxNoDown = PlaneWaveSpecifics.PsiMaxNoUp = PlaneWaveSpecifics.PsiType = 0;
   for (vector<atom *>::iterator runner = allatoms.begin(); runner != allatoms.end(); ++runner) {
     PlaneWaveSpecifics.MaxPsiDouble += (*runner)->type->NoValenceOrbitals;
+    PlaneWaveSpecifics.MaxPsiDouble += (*runner)->getType()->NoValenceOrbitals;
+  }
   cout << PlaneWaveSpecifics.MaxPsiDouble << endl;
 …
   // insert all found elements into the map
   for (vector<atom *>::iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
     Inserter = PresentElements.insert(pair<int, int>((*AtomRunner)->type->Z, 1));
+    Inserter = PresentElements.insert(pair<int, int>((*AtomRunner)->getType()->Z, 1));
     if (!Inserter.second) // increase if present
       Inserter.first->second += 1;
 …
   int nr = 0;
   for (vector<atom *>::iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
     Inserter = ZtoCountMap.insert( pair<int, int>((*AtomRunner)->type->Z, 1) );
+    Inserter = ZtoCountMap.insert( pair<int, int>((*AtomRunner)->getType()->Z, 1) );
     if (!Inserter.second)
       Inserter.first->second += 1;
     const int Z = (*AtomRunner)->type->Z;
+    const int Z = (*AtomRunner)->getType()->Z;
     *file << "Ion_Type" << ZtoIndexMap[Z] << "_" << ZtoCountMap[Z] << "\t"  << fixed << setprecision(9) << showpoint;
     *file << (*AtomRunner)->x[0] << "\t" << (*AtomRunner)->x[1] << "\t" << (*AtomRunner)->x[2];
+    *file << (*AtomRunner)->at(0) << "\t" << (*AtomRunner)->at(1) << "\t" << (*AtomRunner)->at(2);
     *file << "\t" << (*AtomRunner)->FixedIon;
     if ((*AtomRunner)->v.Norm() > MYEPSILON)
       *file << "\t" << scientific << setprecision(6) << (*AtomRunner)->v[0] << "\t" << (*AtomRunner)->v[1] << "\t" << (*AtomRunner)->v[2] << "\t";
+    if ((*AtomRunner)->AtomicVelocity.Norm() > MYEPSILON)
+      *file << "\t" << scientific << setprecision(6) << (*AtomRunner)->AtomicVelocity[0] << "\t" << (*AtomRunner)->AtomicVelocity[1] << "\t" << (*AtomRunner)->AtomicVelocity[2] << "\t";
     *file << " # molecule nr " << nr++ << endl;
+  }

src/Parser/TremoloParser.cpp

-              r8d6d31
+              r8f4df1
       case TremoloKey::x :
         // for the moment, assume there are always three dimensions
         *file << currentAtom->x[0] << "\t";
         *file << currentAtom->x[1] << "\t";
         *file << currentAtom->x[2] << "\t";
+        *file << currentAtom->at(0) << "\t";
+        *file << currentAtom->at(1) << "\t";
+        *file << currentAtom->at(2) << "\t";
         break;
       case TremoloKey::u :
         // for the moment, assume there are always three dimensions
         *file << currentAtom->v[0] << "\t";
         *file << currentAtom->v[1] << "\t";
         *file << currentAtom->v[2] << "\t";
+        *file << currentAtom->AtomicVelocity[0] << "\t";
+        *file << currentAtom->AtomicVelocity[1] << "\t";
+        *file << currentAtom->AtomicVelocity[2] << "\t";
         break;
       case TremoloKey::Type :
 …
   string word;
   int oldId;
+  double tmp;
   lineStream << line;
 …
       case TremoloKey::x :
         // for the moment, assume there are always three dimensions
+        lineStream >> newAtom->x[0];
+        lineStream >> newAtom->x[1];
+        lineStream >> newAtom->x[2];
+        for (int i=0;i<NDIM;i++) {
+          lineStream >> tmp;
+          newAtom->set(i, tmp);
+        }
         break;
       case TremoloKey::u :
         // for the moment, assume there are always three dimensions
         lineStream >> newAtom->v[0];
         lineStream >> newAtom->v[1];
         lineStream >> newAtom->v[2];
+        lineStream >> newAtom->AtomicVelocity[0];
+        lineStream >> newAtom->AtomicVelocity[1];
+        lineStream >> newAtom->AtomicVelocity[2];
         break;
       case TremoloKey::Type :

src/Parser/XyzParser.cpp

-              r8d6d31
+              r8f4df1
   int numberOfAtoms;
   char commentBuffer[512], type[3];
+  double tmp;
   // the first line tells number of atoms, the second line is always a comment
 …
   for (int i = 0; i < numberOfAtoms; i++) {
     newAtom = World::getInstance().createAtom();
+    *file >> type >> ws >> newAtom->x[0] >> ws >> newAtom->x[1] >> ws >> newAtom->x[2];
+    *file >> type;
+    for (int j=0;j<NDIM;j++) {
+      *file >> tmp;
+      newAtom->set(j, tmp);
+    }
     newAtom->setType(World::getInstance().getPeriode()->FindElement(type));
     newmol->AddAtom(newAtom);
 …
   vector<atom*> atoms = World::getInstance().getAllAtoms();
   for(vector<atom*>::iterator it = atoms.begin(); it != atoms.end(); it++) {
     *file << noshowpoint << (*it)->getType()->symbol << "\t" << (*it)->x[0] << "\t" << (*it)->x[1] << "\t" << (*it)->x[2] << endl;
+    *file << noshowpoint << (*it)->getType()->symbol << "\t" << (*it)->at(0) << "\t" << (*it)->at(1) << "\t" << (*it)->at(2) << endl;
+  }
+}

src/UIElements/QT4/QTDialog.cpp

-              r8d6d31
+              r8f4df1
     parent(_parent)
+{
-  const Matrix& M = World::getInstance().getDomain().getM();
-  const char *coords[3] = {"x:","y:","z:"};
   mainLayout= new QHBoxLayout();
   titleLabel = new QLabel(QString(getTitle().c_str()));
 …
     parent(_parent)
+{
-  const Matrix& M = World::getInstance().getDomain().getM();
-  const char *coords[3] = {"x:","y:","z:"};
   mainLayout= new QHBoxLayout();
   titleLabel = new QLabel(QString(getTitle().c_str()));

src/UIElements/TextUI/TextDialog.cpp

r8d6d31	r8f4df1
525	525	}
526	526	}
	527
527	528	return true;
528	529	}

src/analysis_bonds.cpp

-              r8d6d31
+              r8f4df1
   int AtomNo = 0;
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     if ((*iter)->type == type1)
+    if ((*iter)->getType() == type1)
       for (BondList::const_iterator BondRunner = (*iter)->ListOfBonds.begin(); BondRunner != (*iter)->ListOfBonds.end(); BondRunner++)
         if ((*BondRunner)->GetOtherAtom((*iter))->type == type2) {
+        if ((*BondRunner)->GetOtherAtom((*iter))->getType() == type2) {
           const double distance = (*BondRunner)->GetDistanceSquared();
           if (Min > distance)
 …
  * \return angle between \a *first and \a *second, both relative to origin at \a *origin.
  */
 double CalculateAngle(Vector *first, Vector *central, Vector *second)
+double CalculateAngle(const Vector &first, const Vector &central, const Vector &second)
+{
   Vector OHBond;
   Vector OOBond;
   OHBond = (*first) - (*central);
   OOBond = (*second) - (*central);
+  OHBond = first - central;
+  OOBond = second - central;
   const double angle = OHBond.Angle(OOBond);
   return angle;
 …
   // check angle
   if (CalculateAngle(&Hydrogen->x, &Oxygen->x, &OtherOxygen->x) < M_PI*(30./180.)) {
+  if (CalculateAngle(Hydrogen->getPosition(), Oxygen->getPosition(), OtherOxygen->getPosition()) < M_PI*(30./180.)) {
     return true;
   } else {
 …
         molecule::iterator Runner = (*MolRunner)->begin();
         for(;Runner!=(*MolRunner)->end();++Runner){
           if (((*Walker)->type->Z  == 8) && ((*Runner)->type->Z  == 8)) {
+          if (((*Walker)->getType()->Z  == 8) && ((*Runner)->getType()->Z  == 8)) {
             // check distance
             const double distance = (*Runner)->x.DistanceSquared((*Walker)->x);
+            const double distance = (*Runner)->DistanceSquared(*(*Walker));
             if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means  different atoms
               // on other atom(Runner) we check for bond to interface element and
 …
                 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
                 // if hydrogen, check angle to be greater(!) than 30 degrees
                 if (OtherAtom->type->Z == 1) {
                   const double angle = CalculateAngle(&OtherAtom->x, &(*Runner)->x, &(*Walker)->x);
+                if (OtherAtom->getType()->Z == 1) {
+                  const double angle = CalculateAngle(OtherAtom->getPosition(), (*Runner)->getPosition(), (*Walker)->getPosition());
                   OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
                   Otherangle += angle;
                   OtherHydrogens++;
+                }
                 InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement);
                 Interface2Flag = Interface2Flag || (OtherAtom->type == Interface2Element);
+                InterfaceFlag = InterfaceFlag || (OtherAtom->getType() == InterfaceElement);
+                Interface2Flag = Interface2Flag || (OtherAtom->getType() == Interface2Element);
+              }
               DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl);
 …
                 for (BondList::const_iterator BondRunner = (*Walker)->ListOfBonds.begin(); BondRunner != (*Walker)->ListOfBonds.end(); BondRunner++) {
                   atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
                   if (OtherAtom->type->Z == 1) {
+                  if (OtherAtom->getType()->Z == 1) {
                     // check angle
                     if (CheckHydrogenBridgeBondAngle(*Walker, OtherAtom, *Runner)) {
                       DoLog(1) && (Log() << Verbose(1) << (*Walker)->getName() << ", " << OtherAtom->getName() << " and " << (*Runner)->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &(*Walker)->x, &(*Runner)->x)*(180./M_PI) << "." << endl);
+                      DoLog(1) && (Log() << Verbose(1) << (*Walker)->getName() << ", " << OtherAtom->getName() << " and " << (*Runner)->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(OtherAtom->getPosition(), (*Walker)->getPosition(), (*Runner)->getPosition())*(180./M_PI) << "." << endl);
                       count++;
                       break;
 …
     for(;Walker!=(*MolWalker)->end();++Walker){
       atom * theAtom = *Walker;
       if ((theAtom->type == first) || (theAtom->type == second)) {  // first element matches
+      if ((theAtom->getType() == first) || (theAtom->getType() == second)) {  // first element matches
         for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
           if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (theAtom->nr < OtherAtom->nr)) {
+          if (((OtherAtom->getType() == first) || (OtherAtom->getType() == second)) && (theAtom->nr < OtherAtom->nr)) {
             count++;
             DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl);
 …
     for(;Walker!=(*MolWalker)->end();++Walker){
       atom *theAtom = *Walker;
       if (theAtom->type == second) {  // first element matches
+      if (theAtom->getType() == second) {  // first element matches
         for (int i=0;i<2;i++)
           MatchFlag[i] = false;
 …
           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
           for (int i=0;i<2;i++)
             if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) {
+            if ((!MatchFlag[i]) && (OtherAtom->getType() == ElementArray[i])) {
               MatchFlag[i] = true;
               break;  // each bonding atom can match at most one element we are looking for

src/analysis_correlation.cpp

-              r8d6d31
+              r8f4df1
 #include <iomanip>
+#include "BoundaryTriangleSet.hpp"
 #include "analysis_correlation.hpp"
 #include "element.hpp"
 …
           if ((*iter)->getId() < (*runner)->getId()){
             for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
               if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
                 distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
+              if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
+                distance = domain.periodicDistance((*iter)->getPosition(),(*runner)->getPosition());
                 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
 …
     for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
       DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
       periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
+      periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
       // go through every range in xyz and get distance
       for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
 …
                   if ((*iter)->getId() < (*runner)->getId()){
                     for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                       if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
                         periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
+                      if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
+                        periodicOtherX = FullInverseMatrix * ((*runner)->getPosition()); // x now in [0,1)^3
                         // go through every range in xyz and get distance
                         for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
 …
       DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
       for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
         if ((*type == NULL) || ((*iter)->type == *type)) {
           distance = domain.periodicDistance(*(*iter)->node,*point);
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          distance = domain.periodicDistance((*iter)->getPosition(),*point);
           DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
 …
       DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
       for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
         if ((*type == NULL) || ((*iter)->type == *type)) {
           periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
           // go through every range in xyz and get distance
           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
 …
       DoLog(3) && (Log() << Verbose(3) << "\tCurrent atom is " << *(*iter) << "." << endl);
       for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
         if ((*type == NULL) || ((*iter)->type == *type)) {
           TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          TriangleIntersectionList Intersections((*iter)->getPosition(),Surface,LC);
           distance = Intersections.GetSmallestDistance();
           triangle = Intersections.GetClosestTriangle();
 …
       DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
       for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
         if ((*type == NULL) || ((*iter)->type == *type)) {
           periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
           // go through every range in xyz and get distance
           ShortestDistance = -1.;
 …
               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
                 checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                 TriangleIntersectionList Intersections(&checkX,Surface,LC);
+                TriangleIntersectionList Intersections(checkX,Surface,LC);
                 distance = Intersections.GetSmallestDistance();
                 triangle = Intersections.GetClosestTriangle();
 …
     *file << runner->first;
     for (int i=0;i<NDIM;i++)
       *file << "\t" << setprecision(8) << (runner->second.first->node->at(i) - runner->second.second->at(i));
+      *file << "\t" << setprecision(8) << (runner->second.first->at(i) - runner->second.second->at(i));
     *file << endl;
+  }
 …
   if (!map->empty())
     for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
+      *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+    }
+};
+      *file << setprecision(8) << runner->first << "\t";
+      *file << *(runner->second.first) << "\t";
+      *file << *(runner->second.second) << endl;
+    }
+};

src/atom.cpp

-              r8d6d31
+              r8f4df1
 atom::atom() :
   father(this), sort(&nr), mol(0)
+{
+  node = &x;  // TesselPoint::x can only be referenced from here
+};
+{};
 /** Constructor of class atom.
 …
     ParticleInfo(pointer),father(pointer), sort(&nr)
+{
   type = pointer->type;  // copy element of atom
   x = pointer->x; // copy coordination
   v = pointer->v; // copy velocity
+  setType(pointer->getType());  // copy element of atom
+  setPosition(pointer->getPosition()); // copy coordination
+  AtomicVelocity = pointer->AtomicVelocity; // copy velocity
   FixedIon = pointer->FixedIon;
-  node = &x;
   mol = 0;
 };
 …
   res->father = this;
   res->sort = &res->nr;
   res->type = type;
   res->x = this->x;
   res->v = this->v;
+  res->setType(getType());
+  res->setPosition(this->getPosition());
+  res->AtomicVelocity = this->AtomicVelocity;
   res->FixedIon = FixedIon;
-  res->node = &x;
   res->mol = 0;
   World::getInstance().registerAtom(res);
 …
 bool atom::IsInShape(const Shape& shape) const
+{
   return shape.isInside(*node);
+  return shape.isInside(getPosition());
 };
 …
   if (out != NULL) {
     *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
     *out << x[0] << "\t" << x[1] << "\t" << x[2];
+    *out << at(0) << "\t" << at(1) << "\t" << at(2);
     *out << "\t" << FixedIon;
     if (v.Norm() > MYEPSILON)
       *out << "\t" << scientific << setprecision(6) << v[0] << "\t" << v[1] << "\t" << v[2] << "\t";
+    if (AtomicVelocity.Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << AtomicVelocity[0] << "\t" << AtomicVelocity[1] << "\t" << AtomicVelocity[2] << "\t";
     if (comment != NULL)
       *out << " # " << comment << endl;
 …
 bool atom::OutputArrayIndexed(ostream * const out,const enumeration<const element*> &elementLookup, int *AtomNo, const char *comment) const
+{
+  AtomNo[type->Z]++;  // increment number
+  if (out != NULL) {
+    cout << "Looking for atom with element " << *type << endl;
+    ASSERT(elementLookup.there.find(type)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration");
+    *out << "Ion_Type" << elementLookup.there.find(type)->second << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
+    *out << x[0] << "\t" << x[1] << "\t" << x[2];
+  AtomNo[getType()->Z]++;  // increment number
+  if (out != NULL) {
+    const element *elemental = getType();
+    cout << "Looking for atom with element " << *elemental << endl;
+    ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration");
+    *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->Z] << "\t"  << fixed << setprecision(9) << showpoint;
+    *out << at(0) << "\t" << at(1) << "\t" << at(2);
     *out << "\t" << FixedIon;
     if (v.Norm() > MYEPSILON)
       *out << "\t" << scientific << setprecision(6) << v[0] << "\t" << v[1] << "\t" << v[2] << "\t";
+    if (AtomicVelocity.Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << AtomicVelocity[0] << "\t" << AtomicVelocity[1] << "\t" << AtomicVelocity[2] << "\t";
     if (comment != NULL)
       *out << " # " << comment << endl;
 …
+{
   if (out != NULL) {
     *out << type->symbol << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << "\t" << endl;
+    *out << getType()->symbol << "\t" << at(0) << "\t" << at(1) << "\t" << at(2) << "\t" << endl;
     return true;
   } else
 …
 bool atom::OutputTrajectory(ofstream * const out, const int *ElementNo, int *AtomNo, const int step) const
+{
   AtomNo[type->Z]++;
   if (out != NULL) {
     *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
+  AtomNo[getType()->Z]++;
+  if (out != NULL) {
+    *out << "Ion_Type" << ElementNo[getType()->Z] << "_" << AtomNo[getType()->Z] << "\t"  << fixed << setprecision(9) << showpoint;
     *out << Trajectory.R.at(step)[0] << "\t" << Trajectory.R.at(step)[1] << "\t" << Trajectory.R.at(step)[2];
     *out << "\t" << FixedIon;
 …
+{
   if (out != NULL) {
     *out << type->symbol << "\t";
+    *out << getType()->symbol << "\t";
     *out << Trajectory.R.at(step)[0] << "\t";
     *out << Trajectory.R.at(step)[1] << "\t";
 …
 void atom::OutputMPQCLine(ostream * const out, const Vector *center, int *AtomNo = NULL) const
+{
+  *out << "\t\t" << type->symbol << " [ " << x[0]-center->at(0) << "\t" << x[1]-center->at(1) << "\t" << x[2]-center->at(2) << " ]" << endl;
+  Vector recentered(getPosition());
+  recentered -= *center;
+  *out << "\t\t" << getType()->symbol << " [ " << recentered[0] << "\t" << recentered[1] << "\t" << recentered[2] << " ]" << endl;
   if (AtomNo != NULL)
     *AtomNo++;
 …
 double atom::DistanceSquaredToVector(const Vector &origin) const
+{
   return origin.DistanceSquared(x);
+  return DistanceSquared(origin);
 };
 …
 double atom::DistanceToVector(const Vector &origin) const
+{
   return origin.distance(x);
+  return distance(origin);
 };
 …
 };
+std::ostream & atom::operator << (std::ostream &ost) const
+{
+  ParticleInfo::operator<<(ost);
+  ost << "," << getPosition();
+  return ost;
+}
+std::ostream & operator << (std::ostream &ost, const atom &a)
+{
+  a.ParticleInfo::operator<<(ost);
+  ost << "," << a.getPosition();
+  return ost;
+}
 bool operator < (atom &a, atom &b)

src/atom.hpp

-              r8d6d31
+              r8f4df1
 #include <vector>
+#include "Helpers/helpers.hpp"
 #include "atom_atominfo.hpp"
 #include "atom_bondedparticle.hpp"
 …
 #include "atom_particleinfo.hpp"
 #include "atom_trajectoryparticle.hpp"
 #include "tesselation.hpp"
+#include "TesselPoint.hpp"
 #include "types.hpp"
 …
  * Class incorporates position, type
  */
 class atom : public TesselPoint, public TrajectoryParticle, public GraphNode, public BondedParticle, public virtual ParticleInfo, public virtual AtomInfo {
+class atom : public TrajectoryParticle, public GraphNode, public BondedParticle, public TesselPoint {
   friend atom* NewAtom(atomId_t);
   friend void  DeleteAtom(atom*);
 …
    void removeFromMolecule();
+   // Output operator
+   std::ostream & operator << (std::ostream &ost) const;
   protected:
 …
 /**
+ * Global output operator for class atom.
+ */
+std::ostream & operator << (std::ostream &ost, const atom &_atom);
+/**
  * internal method used by the world. Do not use if you don't know what you are doing.
  * You might get burned...

src/atom_atominfo.cpp

-              r8d6d31
+              r8f4df1
 #include "periodentafel.hpp"
 #include "World.hpp"
+#include "element.hpp"
 #include "atom_atominfo.hpp"
 /** Constructor of class AtomInfo.
  */
+AtomInfo::AtomInfo() : type(NULL) {};
+AtomInfo::AtomInfo() : AtomicElement(NULL) {};
+/** Copy constructor of class AtomInfo.
+ */
+AtomInfo::AtomInfo(const AtomInfo &_atom) : AtomicPosition(_atom.AtomicPosition), AtomicElement(_atom.AtomicElement) {};
+AtomInfo::AtomInfo(const VectorInterface &_v) : AtomicPosition(_v.getPosition()), AtomicElement(NULL)
+{};
 /** Destructor of class AtomInfo.
 …
 };
+const element *AtomInfo::getType(){
+  return type;
+const element *AtomInfo::getType() const
+{
+  return AtomicElement;
+}
+const double& AtomInfo::operator[](size_t i) const
+{
+  return AtomicPosition[i];
+}
+const double& AtomInfo::at(size_t i) const
+{
+  return AtomicPosition.at(i);
+}
+void AtomInfo::set(size_t i, const double value)
+{
+  AtomicPosition.at(i) = value;
+}
+const Vector& AtomInfo::getPosition() const
+{
+  return AtomicPosition;
+}
 void AtomInfo::setType(const element* _type) {
   type = _type;
+  AtomicElement = _type;
+}
 void AtomInfo::setType(int Z) {
+void AtomInfo::setType(const int Z) {
   const element *elem = World::getInstance().getPeriode()->FindElement(Z);
   setType(elem);
+}
+void AtomInfo::setPosition(const Vector& _vector)
+{
+  AtomicPosition = _vector;
+  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
+}
+const VectorInterface& AtomInfo::operator+=(const Vector& b)
+{
+  AtomicPosition += b;
+  return *this;
+}
+const VectorInterface& AtomInfo::operator-=(const Vector& b)
+{
+  AtomicPosition -= b;
+  return *this;
+}
+Vector const AtomInfo::operator+(const Vector& b) const
+{
+  Vector a(AtomicPosition);
+  a += b;
+  return a;
+}
+Vector const AtomInfo::operator-(const Vector& b) const
+{
+  Vector a(AtomicPosition);
+  a -= b;
+  return a;
+}
+double AtomInfo::distance(const Vector &point) const
+{
+  return AtomicPosition.distance(point);
+}
+double AtomInfo::DistanceSquared(const Vector &y) const
+{
+  return AtomicPosition.DistanceSquared(y);
+}
+double AtomInfo::distance(const VectorInterface &_atom) const
+{
+  return _atom.distance(AtomicPosition);
+}
+double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
+{
+  return _atom.DistanceSquared(AtomicPosition);
+}
+VectorInterface &AtomInfo::operator=(const Vector& _vector)
+{
+  AtomicPosition = _vector;
+  return *this;
+}
+void AtomInfo::ScaleAll(const double *factor)
+{
+  AtomicPosition.ScaleAll(factor);
+}
+void AtomInfo::ScaleAll(const Vector &factor)
+{
+  AtomicPosition.ScaleAll(factor);
+}
+void AtomInfo::Scale(const double factor)
+{
+  AtomicPosition.Scale(factor);
+}
+void AtomInfo::Zero()
+{
+  AtomicPosition.Zero();
+}
+void AtomInfo::One(const double one)
+{
+  AtomicPosition.One(one);
+}
+void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
+{
+  AtomicPosition.LinearCombinationOfVectors(x1,x2,x3,factors);
+}
+const AtomInfo& operator*=(AtomInfo& a, const double m)
+{
+  a.Scale(m);
+  return a;
+}
+AtomInfo const operator*(const AtomInfo& a, const double m)
+{
+  AtomInfo copy(a);
+  copy *= m;
+  return copy;
+}
+AtomInfo const operator*(const double m, const AtomInfo& a)
+{
+  AtomInfo copy(a);
+  copy *= m;
+  return copy;
+}
+std::ostream & AtomInfo::operator << (std::ostream &ost) const
+{
+  return (ost << getPosition());
+}
+std::ostream & operator << (std::ostream &ost, const AtomInfo &a)
+{
+  ost << a;
+  return ost;
+}

src/atom_atominfo.hpp

-              r8d6d31
+              r8f4df1
 #include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/VectorInterface.hpp"
 /****************************************** forward declarations *****************************/
+class AtomInfo;
 class element;
+class Matrix;
 /********************************************** declarations *******************************/
+class AtomInfo {
+class AtomInfo : public VectorInterface {
 public:
+  Vector x;       //!< coordinate vector of atom, giving last position within cell
+  Vector v;       //!< velocity vector of atom, giving last velocity within cell
+  Vector F;       //!< Force vector of atom, giving last force within cell
+  const element *type;  //!< pointing to element
+  Vector AtomicVelocity;       //!< velocity vector of atom, giving last velocity within cell
+  Vector AtomicForce;       //!< Force vector of atom, giving last force within cell
   AtomInfo();
+  ~AtomInfo();
+  AtomInfo(const AtomInfo &_atom);
+  AtomInfo(const VectorInterface &_v);
+  virtual ~AtomInfo();
   const element *getType();
+  const element *getType() const;
   void setType(const element *);
+  void setType(int);
+  void setType(const int);
+  ///// manipulation of the atomic position
+  // Accessors ussually come in pairs... and sometimes even more than that
+  const double& operator[](size_t i) const;
+  const double& at(size_t i) const;
+  void set(size_t i, const double value);
+  const Vector& getPosition() const;
+  // Assignment operator
+  void setPosition(const Vector& _vector);
+  class VectorInterface &operator=(const Vector& _vector);
+  // operators for mathematical operations
+  const VectorInterface& operator+=(const Vector& b);
+  const VectorInterface& operator-=(const Vector& b);
+  Vector const operator+(const Vector& b) const;
+  Vector const operator-(const Vector& b) const;
+  void Zero();
+  void One(const double one);
+  void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
+  double distance(const Vector &point) const;
+  double DistanceSquared(const Vector &y) const;
+  double distance(const VectorInterface &_atom) const;
+  double DistanceSquared(const VectorInterface &_atom) const;
+  void ScaleAll(const double *factor);
+  void ScaleAll(const Vector &factor);
+  void Scale(const double factor);
+  std::ostream & operator << (std::ostream &ost) const;
 private:
+  Vector AtomicPosition;       //!< coordinate vector of atom, giving last position within cell
+  const element *AtomicElement;  //!< pointing to element
 };
+std::ostream & operator << (std::ostream &ost, const AtomInfo &a);
+const AtomInfo& operator*=(AtomInfo& a, const double m);
+AtomInfo const operator*(const AtomInfo& a, const double m);
+AtomInfo const operator*(const double m, const AtomInfo& a);
 #endif /* ATOM_ATOMINFO_HPP_ */

src/atom_bondedparticle.cpp

-              r8d6d31
+              r8f4df1
   NoBonds = CountBonds();
   //Log() << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
   if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
+  if ((int)(getType()->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
     for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
       OtherWalker = (*Runner)->GetOtherAtom(this);
       OtherNoBonds = OtherWalker->CountBonds();
       //Log() << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << OtherNoBonds << "?" << endl;
       if ((int)(OtherWalker->type->NoValenceOrbitals) > OtherNoBonds) { // check if possible candidate
+      if ((int)(OtherWalker->getType()->NoValenceOrbitals) > OtherNoBonds) { // check if possible candidate
         if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
           CandidateBond = (*Runner);
 …
 };
+std::ostream & BondedParticle::operator << (std::ostream &ost) const
+{
+  ParticleInfo::operator<<(ost);
+  ost << "," << getPosition();
+  return ost;
+}
+std::ostream & operator << (std::ostream &ost, const BondedParticle &a)
+{
+  a.ParticleInfo::operator<<(ost);
+  ost << "," << a.getPosition();
+  return ost;
+}

src/atom_bondedparticle.hpp

-              r8d6d31
+              r8f4df1
   void OutputOrder(ofstream *file) const;
+  std::ostream & operator << (std::ostream &ost) const;
 private:
 };
+std::ostream & operator << (std::ostream &ost, const BondedParticle &a);
 #endif /* ATOM_BONDEDPARTICLE_HPP_ */

src/atom_bondedparticleinfo.hpp

r8d6d31	r8f4df1
36	36
37	37	BondedParticleInfo();
38		~BondedParticleInfo();
	38	virtual ~BondedParticleInfo();
39	39
40	40	private:

src/atom_graphnodeinfo.hpp

r8d6d31	r8f4df1
37	37
38	38	GraphNodeInfo();
39		~GraphNodeInfo();
	39	virtual ~GraphNodeInfo();
40	40	private:
41	41

src/atom_particleinfo.cpp

-              r8d6d31
+              r8f4df1
 ostream & operator << (ostream &ost, const ParticleInfo &a)
+{
   ost << "[" << a.getName() << "|" << &a << "]";
+  ost << a.getName();
   return ost;
 };
 …
 ostream & ParticleInfo::operator << (ostream &ost) const
+{
   ost << "[" << name << "|" << this << "]";
+  ost << getName();
   return ost;
 };

src/atom_particleinfo.hpp

r8d6d31	r8f4df1
32	32	ParticleInfo();
33	33	ParticleInfo(ParticleInfo*);
34		~ParticleInfo();
	34	virtual ~ParticleInfo();
35	35
36	36	const std::string& getName() const;

src/atom_trajectoryparticle.cpp

-              r8d6d31
+              r8f4df1
+{
   for (int i=NDIM;i--;)
     *temperature += type->mass * Trajectory.U.at(step)[i]* Trajectory.U.at(step)[i];
+    *temperature += getType()->mass * Trajectory.U.at(step)[i]* Trajectory.U.at(step)[i];
 };
 …
   for(int d=0;d<NDIM;d++) {
     Trajectory.U.at(Step)[d] -= CoGVelocity->at(d);
     *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step)[d] * Trajectory.U.at(Step)[d];
+    *ActualTemp += 0.5 * getType()->mass * Trajectory.U.at(Step)[d] * Trajectory.U.at(Step)[d];
+  }
 };
 …
     Trajectory.R.at(NextStep)[d] = Trajectory.R.at(NextStep-1)[d];
     Trajectory.R.at(NextStep)[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
     Trajectory.R.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep)[d]/type->mass);     // F = m * a and s =
+    Trajectory.R.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep)[d]/getType()->mass);     // F = m * a and s =
+  }
   // Update U
   for (int d=0; d<NDIM; d++) {
     Trajectory.U.at(NextStep)[d] = Trajectory.U.at(NextStep-1)[d];
     Trajectory.U.at(NextStep)[d] += configuration->Deltat * (Trajectory.F.at(NextStep)[d]+Trajectory.F.at(NextStep-1)[d]/type->mass); // v = F/m * t
+    Trajectory.U.at(NextStep)[d] += configuration->Deltat * (Trajectory.F.at(NextStep)[d]+Trajectory.F.at(NextStep-1)[d]/getType()->mass); // v = F/m * t
+  }
   // Update R (and F)
 …
 void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity ) const
+{
   *TotalMass += type->mass;  // sum up total mass
+  *TotalMass += getType()->mass;  // sum up total mass
   for(int d=0;d<NDIM;d++) {
     TotalVelocity->at(d) += Trajectory.U.at(Step)[d]*type->mass;
+    TotalVelocity->at(d) += Trajectory.U.at(Step)[d]*getType()->mass;
+  }
 };
 …
     for (int d=0; d<NDIM; d++) {
       U[d] *= ScaleTempFactor;
       *ekin += 0.5*type->mass * U[d]*U[d];
+      *ekin += 0.5*getType()->mass * U[d]*U[d];
+    }
 };
 …
     for (int d=0; d<NDIM; d++) {
       *G += U[d] * F[d];
       *E += U[d]*U[d]*type->mass;
+      *E += U[d]*U[d]*getType()->mass;
+    }
 };
 …
   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
       U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
       *ekin += type->mass * U[d]*U[d];
+      U[d] += configuration->Deltat/getType()->mass * ( (G_over_E) * (U[d]*getType()->mass) );
+      *ekin += getType()->mass * U[d]*U[d];
+    }
 };
 …
 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
+{
   double sigma  = sqrt(configuration->Thermostats->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+  double sigma  = sqrt(configuration->Thermostats->TargetTemp/getType()->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
   Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
 …
+    }
     for (int d=0; d<NDIM; d++)
       *ekin += 0.5*type->mass * U[d]*U[d];
+      *ekin += 0.5*getType()->mass * U[d]*U[d];
+  }
 };
 …
     for (int d=0; d<NDIM; d++) {
       U[d] *= sqrt(1+(configuration->Deltat/configuration->Thermostats->TempFrequency)*(ScaleTempFactor-1));
       *ekin += 0.5*type->mass * U[d]*U[d];
+      *ekin += 0.5*getType()->mass * U[d]*U[d];
+    }
+  }
 …
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
       *delta_alpha += U[d]*U[d]*type->mass;
+      *delta_alpha += U[d]*U[d]*getType()->mass;
+    }
+  }
 …
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
         U[d] += configuration->Deltat/type->mass * (configuration->Thermostats->alpha * (U[d] * type->mass));
         *ekin += (0.5*type->mass) * U[d]*U[d];
+        U[d] += configuration->Deltat/getType()->mass * (configuration->Thermostats->alpha * (U[d] * getType()->mass));
+        *ekin += (0.5*getType()->mass) * U[d]*U[d];
+      }
+  }
 };
+std::ostream & TrajectoryParticle::operator << (std::ostream &ost) const
+{
+  ParticleInfo::operator<<(ost);
+  ost << "," << getPosition();
+  return ost;
+}
+std::ostream & operator << (std::ostream &ost, const TrajectoryParticle &a)
+{
+  a.ParticleInfo::operator<<(ost);
+  ost << "," << a.getPosition();
+  return ost;
+}

src/atom_trajectoryparticle.hpp

-              r8d6d31
+              r8f4df1
   void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
+  std::ostream & operator << (std::ostream &ost) const;
 private:
 };
+ostream & operator << (ostream &ost, const TrajectoryParticle &a);
 #endif /* ATOM_TRAJECTORYPARTICLE_HPP_ */

src/atom_trajectoryparticleinfo.hpp

r8d6d31	r8f4df1
36	36
37	37	TrajectoryParticleInfo();
38		~TrajectoryParticleInfo();
	38	virtual ~TrajectoryParticleInfo();
39	39
40	40	private:

src/bond.cpp

-              r8d6d31
+              r8f4df1
 #include "Helpers/MemDebug.hpp"
+#include "Helpers/Log.hpp"
 #include "Helpers/Verbose.hpp"
 #include "atom.hpp"
 …
+{
   if ((left != NULL) && (right != NULL)) {
     if ((left->type != NULL) && (left->type->Z == 1))
+    if ((left->getType() != NULL) && (left->getType()->Z == 1))
       HydrogenBond++;
     if ((right->type != NULL) && (right->type->Z == 1))
+    if ((right->getType() != NULL) && (right->getType()->Z == 1))
       HydrogenBond++;
+  }
 …
 double bond::GetDistance() const
+{
   return (leftatom->node->distance(*rightatom->node));
+  return (leftatom->distance(*rightatom));
 };
 …
 double bond::GetDistanceSquared() const
+{
   return (leftatom->node->DistanceSquared(*rightatom->node));
+  return (leftatom->DistanceSquared(*rightatom));
 };

src/bondgraph.cpp

-              r8d6d31
+              r8f4df1
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     if ((*iter)->type->CovalentRadius > max_distance)
       max_distance = (*iter)->type->CovalentRadius;
+    if ((*iter)->getType()->CovalentRadius > max_distance)
+      max_distance = (*iter)->getType()->CovalentRadius;
+  }
   max_distance *= 2.;
 …
 void BondGraph::CovalentMinMaxDistance(BondedParticle * const Walker, BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem)
+{
   MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
+  MinDistance = OtherWalker->getType()->CovalentRadius + Walker->getType()->CovalentRadius;
   MinDistance *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
   MaxDistance = MinDistance + BONDTHRESHOLD;
 …
     CovalentMinMaxDistance(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
   } else {
     MinDistance = GetBondLength(Walker->type->Z-1, OtherWalker->type->Z-1);
+    MinDistance = GetBondLength(Walker->getType()->Z-1, OtherWalker->getType()->Z-1);
     MinDistance *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
     MaxDistance = MinDistance + BONDTHRESHOLD;

src/boundary.cpp

-              r8d6d31
+              r8f4df1
 #include "Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp"
+#include "BoundaryPointSet.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "CandidateForTesselation.hpp"
+//#include "TesselPoint.hpp"
 #include "World.hpp"
 #include "atom.hpp"
 …
               if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                 Neighbour = BoundaryPoints[axis].begin();
               DistanceVector = runner->second.second->x - Neighbour->second.second->x;
+              DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
               do { // seek for neighbour pair where it flips
                   OldComponent = DistanceVector[Othercomponent];
 …
                   if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                     Neighbour = BoundaryPoints[axis].begin();
                   DistanceVector = runner->second.second->x - Neighbour->second.second->x;
+                  DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
                   //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
 …
                   //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
                   // now we have found the pair: Neighbour and OtherNeighbour
                   OtherVector = runner->second.second->x - OtherNeighbour->second.second->x;
+                  OtherVector = (runner->second.second->getPosition()) - (OtherNeighbour->second.second->getPosition());
                   //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
                   //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
 …
     // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
     for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       ProjectedVector = (*iter)->x - (*MolCenter);
+      ProjectedVector = (*iter)->getPosition() - (*MolCenter);
       ProjectedVector.ProjectOntoPlane(AxisVector);
 …
           DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
         } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
+          helper = (*iter)->x;
+          helper -= *MolCenter;
+          helper = (*iter)->getPosition() - (*MolCenter);
           const double oldhelperNorm = helper.NormSquared();
           helper = BoundaryTestPair.first->second.second->x - (*MolCenter);
+          helper = BoundaryTestPair.first->second.second->getPosition() - (*MolCenter);
           if (helper.NormSquared() < oldhelperNorm) {
             BoundaryTestPair.first->second.second = (*iter);
 …
+        {
           Vector SideA, SideB, SideC, SideH;
           SideA = left->second.second->x - (*MolCenter);
+          SideA = left->second.second->getPosition() - (*MolCenter);
           SideA.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideA: " << SideA << endl;
           SideB = right->second.second->x -(*MolCenter);
+          SideB = right->second.second->getPosition() -(*MolCenter);
           SideB.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideB: " << SideB << endl;
           SideC = left->second.second->x - right->second.second->x;
+          SideC = left->second.second->getPosition() - right->second.second->getPosition();
           SideC.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideC: " << SideC << endl;
           SideH = runner->second.second->x -(*MolCenter);
+          SideH = runner->second.second->getPosition() -(*MolCenter);
           SideH.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideH: " << SideH << endl;
 …
   // sum up the atomic masses
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       totalmass += (*iter)->type->mass;
+      totalmass += (*iter)->getType()->mass;
+  }
   DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl);
 …
           // ... and put at new position
           Inserter = (*iter)->x;
+          Inserter = (*iter)->getPosition();
           if (DoRandomRotation)
             Inserter *= Rotations;
 …
             // copy atom ...
             CopyAtoms[(*iter)->nr] = (*iter)->clone();
             CopyAtoms[(*iter)->nr]->x = Inserter;
+            (*CopyAtoms[(*iter)->nr]).setPosition(Inserter);
             Filling->AddAtom(CopyAtoms[(*iter)->nr]);
             DoLog(4) && (Log() << Verbose(4) << "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->nr]->x) << "." << endl);
+            DoLog(1) && (Log() << Verbose(1) << "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->nr]->getPosition()) << "." << endl);
           } else {
             DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl);

src/config.cpp

-              r8d6d31
+              r8f4df1
     map<int, atom *> LinearList;
     atom *neues = NULL;
+    Vector position;
     if (!FastParsing) {
       // parse in trajectories
 …
               AtomList[i][j] = neues;
               LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
               neues->type = elementhash[i]; // find element type
+              neues->setType(elementhash[i]); // find element type
             } else
               neues = AtomList[i][j];
             status = (status &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x[0], 1, (repetition == 0) ? critical : optional) &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x[1], 1, (repetition == 0) ? critical : optional) &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x[2], 1, (repetition == 0) ? critical : optional) &&
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &position[0], 1, (repetition == 0) ? critical : optional) &&
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], 1, (repetition == 0) ? critical : optional) &&
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], 1, (repetition == 0) ? critical : optional) &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
+            if (!status) break;
+            if (!status)
+              break;
+            neues ->setPosition(position);
             // check size of vectors
 …
             // put into trajectories list
             for (int d=0;d<NDIM;d++)
               neues->Trajectory.R.at(repetition)[d] = neues->x[d];
+              neues->Trajectory.R.at(repetition)[d] = neues->at(d);
             // parse velocities if present
             if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v[0], 1,optional))
               neues->v[0] = 0.;
             if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v[1], 1,optional))
               neues->v[1] = 0.;
             if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v[2], 1,optional))
               neues->v[2] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->AtomicVelocity[0], 1,optional))
+              neues->AtomicVelocity[0] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->AtomicVelocity[1], 1,optional))
+              neues->AtomicVelocity[1] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->AtomicVelocity[2], 1,optional))
+              neues->AtomicVelocity[2] = 0.;
             for (int d=0;d<NDIM;d++)
               neues->Trajectory.U.at(repetition)[d] = neues->v[d];
+              neues->Trajectory.U.at(repetition)[d] = neues->AtomicVelocity[d];
             // parse forces if present
 …
             AtomList[i][j] = neues;
             LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
             neues->type = elementhash[i]; // find element type
+            neues->setType(elementhash[i]); // find element type
           } else
             neues = AtomList[i][j];
           // then parse for each atom the coordinates as often as present
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x[0], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x[1], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x[2], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &position[0], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], repetition,critical);
+          neues->setPosition(position);
           ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
           if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v[0], repetition,optional))
             neues->v[0] = 0.;
           if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v[1], repetition,optional))
             neues->v[1] = 0.;
           if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v[2], repetition,optional))
             neues->v[2] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->AtomicVelocity[0], repetition,optional))
+            neues->AtomicVelocity[0] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->AtomicVelocity[1], repetition,optional))
+            neues->AtomicVelocity[1] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->AtomicVelocity[2], repetition,optional))
+            neues->AtomicVelocity[2] = 0.;
           // here we don't care if forces are present (last in trajectories is always equal to current position)
           neues->type = elementhash[i]; // find element type
+          neues->setType(elementhash[i]); // find element type
           mol->AddAtom(neues);
+        }
 …
         istringstream input2(zeile);
         atom *neues = World::getInstance().createAtom();
+        input2 >> neues->x[0]; // x
+        input2 >> neues->x[1]; // y
+        input2 >> neues->x[2]; // z
+        double tmp;
+        for (int j=0;j<NDIM;j++) {
+          input2 >> tmp;
+          neues->set(j,tmp);
+        }
         input2 >> l;
         neues->type = elementhash[No]; // find element type
+        neues->setType(elementhash[No]); // find element type
         mol->AddAtom(neues);
+      }
 …
     AtomNo = 0;
     for (molecule::const_iterator iter = (*MolRunner)->begin(); iter != (*MolRunner)->end(); ++iter) {
       sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
       elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
+      sprintf(name, "%2s%2d",(*iter)->getType()->symbol, elementNo[(*iter)->getType()->Z]);
+      elementNo[(*iter)->getType()->Z] = (elementNo[(*iter)->getType()->Z]+1) % 100;   // confine to two digits
       fprintf(f,
              "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
 …
              'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
              MolNo,         /* residue sequence number */
              (*iter)->node->at(0),                 /* position X in Angstroem */
              (*iter)->node->at(1),                 /* position Y in Angstroem */
              (*iter)->node->at(2),                 /* position Z in Angstroem */
              (double)(*iter)->type->Valence,         /* occupancy */
              (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
+             (*iter)->at(0),                 /* position X in Angstroem */
+             (*iter)->at(1),                 /* position Y in Angstroem */
+             (*iter)->at(2),                 /* position Z in Angstroem */
+             (double)(*iter)->getType()->Valence,         /* occupancy */
+             (double)(*iter)->getType()->NoValenceOrbitals,          /* temperature factor */
              "0",            /* segment identifier */
              (*iter)->type->symbol,    /* element symbol */
+             (*iter)->getType()->symbol,    /* element symbol */
              "0");           /* charge */
       AtomNo++;
 …
   AtomNo = 0;
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
     elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
+    sprintf(name, "%2s%2d",(*iter)->getType()->symbol, elementNo[(*iter)->getType()->Z]);
+    elementNo[(*iter)->getType()->Z] = (elementNo[(*iter)->getType()->Z]+1) % 100;   // confine to two digits
     fprintf(f,
            "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
 …
            'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
 ,         /* residue sequence number */
            (*iter)->node->at(0),                 /* position X in Angstroem */
            (*iter)->node->at(1),                 /* position Y in Angstroem */
            (*iter)->node->at(2),                 /* position Z in Angstroem */
            (double)(*iter)->type->Valence,         /* occupancy */
            (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
+           (*iter)->at(0),                 /* position X in Angstroem */
+           (*iter)->at(1),                 /* position Y in Angstroem */
+           (*iter)->at(2),                 /* position Z in Angstroem */
+           (double)(*iter)->getType()->Valence,         /* occupancy */
+           (double)(*iter)->getType()->NoValenceOrbitals,          /* temperature factor */
            "0",            /* segment identifier */
            (*iter)->type->symbol,    /* element symbol */
+           (*iter)->getType()->symbol,    /* element symbol */
            "0");           /* charge */
     AtomNo++;
 …
     *output << mol->name << "\t";
     *output << 0 << "\t";
     *output << (*iter)->node->at(0) << "\t" << (*iter)->node->at(1) << "\t" << (*iter)->node->at(2) << "\t";
     *output << static_cast<double>((*iter)->type->Valence) << "\t";
     *output << (*iter)->type->symbol << "\t";
+    *output << (*iter)->at(0) << "\t" << (*iter)->at(1) << "\t" << (*iter)->at(2) << "\t";
+    *output << static_cast<double>((*iter)->getType()->Valence) << "\t";
+    *output << (*iter)->getType()->symbol << "\t";
     for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
       *output << (*runner)->GetOtherAtom(*iter)->nr << "\t";
 …
         *output << (*MolWalker)->name << "\t";
         *output << MolCounter+1 << "\t";
         *output << (*iter)->node->at(0) << "\t" << (*iter)->node->at(1) << "\t" << (*iter)->node->at(2) << "\t";
         *output << (double)(*iter)->type->Valence << "\t";
         *output << (*iter)->type->symbol << "\t";
+        *output << (*iter)->at(0) << "\t" << (*iter)->at(1) << "\t" << (*iter)->at(2) << "\t";
+        *output << (double)(*iter)->getType()->Valence << "\t";
+        *output << (*iter)->getType()->symbol << "\t";
         for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
           *output << LocalNotoGlobalNoMap[ (*runner)->GetOtherAtom((*iter))->getId() ] << "\t";

src/ellipsoid.cpp

-              r8d6d31
+              r8f4df1
 #include <set>
+#include "BoundaryPointSet.hpp"
 #include "boundary.hpp"
 #include "ellipsoid.hpp"
 …
               if ((current != PickedAtomNrs.end()) && (*current == index)) {
                 Candidate = (*Runner);
                 DoLog(2) && (Log() << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl);
                 x[PointsPicked++] = *Candidate->node;    // we have one more atom picked
+                DoLog(2) && (Log() << Verbose(2) << "Current picked node is " << (*Runner)->getName() << " with index " << index << "." << endl);
+                x[PointsPicked++] = Candidate->getPosition();    // we have one more atom picked
                 current++;    // next pre-picked atom
+              }
 …
       //Log() << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
       if (value > threshold) {
         x[PointsPicked] = (*Runner->second->node->node);
+        x[PointsPicked] = (Runner->second->node->getPosition());
         PointsPicked++;
         //Log() << Verbose(0) << "IN." << endl;
 …
   Center.Zero();
   for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
     Center += (*Runner->second->node->node);
+    Center += (Runner->second->node->getPosition());
   Center.Scale(1./T->PointsOnBoundaryCount);
   DoLog(1) && (Log() << Verbose(1) << "Center is at " << Center << "." << endl);

src/linkedcell.cpp

-              r8d6d31
+              r8f4df1
   Walker = set->GetPoint();
   for (int i=0;i<NDIM;i++) {
     max[i] = Walker->node->at(i);
     min[i] = Walker->node->at(i);
+    max[i] = Walker->at(i);
+    min[i] = Walker->at(i);
+  }
   set->GoToFirst();
 …
     Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
       if (max[i] < Walker->node->at(i))
         max[i] = Walker->node->at(i);
       if (min[i] > Walker->node->at(i))
         min[i] = Walker->node->at(i);
+      if (max[i] < Walker->at(i))
+        max[i] = Walker->at(i);
+      if (min[i] > Walker->at(i))
+        min[i] = Walker->at(i);
+    }
     set->GoToNext();
 …
     Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
       n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
+      n[i] = static_cast<int>(floor((Walker->at(i) - min[i])/RADIUS));
+    }
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
 …
   LinkedNodes::iterator Runner = set->begin();
   for (int i=0;i<NDIM;i++) {
     max[i] = (*Runner)->node->at(i);
     min[i] = (*Runner)->node->at(i);
+    max[i] = (*Runner)->at(i);
+    min[i] = (*Runner)->at(i);
+  }
   for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
     Walker = *Runner;
     for (int i=0;i<NDIM;i++) {
       if (max[i] < Walker->node->at(i))
         max[i] = Walker->node->at(i);
       if (min[i] > Walker->node->at(i))
         min[i] = Walker->node->at(i);
+      if (max[i] < Walker->at(i))
+        max[i] = Walker->at(i);
+      if (min[i] > Walker->at(i))
+        min[i] = Walker->at(i);
+    }
+  }
 …
     Walker = *Runner;
     for (int i=0;i<NDIM;i++) {
       n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
+      n[i] = static_cast<int>(floor((Walker->at(i) - min[i])/RADIUS));
+    }
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
 …
  * \return Vector is inside bounding box - true, else - false
  */
 bool LinkedCell::SetIndexToVector(const Vector * const x) const
+bool LinkedCell::SetIndexToVector(const Vector & x) const
+{
   for (int i=0;i<NDIM;i++)
     n[i] = (int)floor((x->at(i) - min[i])/RADIUS);
+    n[i] = (int)floor((x.at(i) - min[i])/RADIUS);
   return CheckBounds();
 …
   bool status = false;
   for (int i=0;i<NDIM;i++) {
     n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
+    n[i] = static_cast<int>(floor((Walker->at(i) - min[i])/RADIUS));
+  }
   index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
 …
       Walker = *Runner;
       //Log() << Verbose(1) << "Current neighbour is at " << *Walker->node << "." << endl;
       if ((center->DistanceSquared(*Walker->node) - radiusSquared) < MYEPSILON) {
+      if ((Walker->DistanceSquared(*center) - radiusSquared) < MYEPSILON) {
         TesselList->push_back(Walker);
+      }

src/linkedcell.hpp

r8d6d31	r8f4df1
60	60	const LinkedCell::LinkedNodes* GetRelativeToCurrentCell(const int relative[NDIM])const ;
61	61	bool SetIndexToNode(const TesselPoint * const Walker)const ;
62		bool SetIndexToVector(const Vector * const x)const ;
	62	bool SetIndexToVector(const Vector &x)const ;
63	63	double SetClosestIndexToOutsideVector(const Vector * const x) const;
64	64	bool CheckBounds()const ;

src/molecule.cpp

-              r8d6d31
+              r8f4df1
   atomIds.erase( atom->getId() );
   atoms.remove( atom );
   formula-=atom->type;
+  formula-=atom->getType();
   atom->removeFromMolecule();
   return iter;
 …
     atomIds.erase( key->getId() );
     atoms.remove( key );
     formula-=key->type;
+    formula-=key->getType();
     key->removeFromMolecule();
+  }
 …
   if (res.second) { // push atom if went well
     atoms.push_back(key);
     formula+=key->type;
+    formula+=key->getType();
     return pair<iterator,bool>(molecule::iterator(--end()),res.second);
   } else {
 …
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
     if (pointer->type != NULL) {
       formula += pointer->type;
       if (pointer->type->Z != 1)
+    if (pointer->getType() != NULL) {
+      formula += pointer->getType();
+      if (pointer->getType()->Z != 1)
         NoNonHydrogen++;
       if(pointer->getName() == "Unknown"){
         stringstream sstr;
         sstr << pointer->type->symbol << pointer->nr+1;
+        sstr << pointer->getType()->symbol << pointer->nr+1;
         pointer->setName(sstr.str());
+      }
 …
   if (pointer != NULL) {
     atom *walker = pointer->clone();
     formula += walker->type;
+    formula += walker->getType();
     walker->setName(pointer->getName());
     walker->nr = last_atom++;  // increase number within molecule
     insert(walker);
     if ((pointer->type != NULL) && (pointer->type->Z != 1))
+    if ((pointer->getType() != NULL) && (pointer->getType()->Z != 1))
       NoNonHydrogen++;
     retval=walker;
 …
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
   InBondvector = TopReplacement->x - TopOrigin->x;
+  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
   bondlength = InBondvector.Norm();
 …
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
       l = TopReplacement->x[i] - TopOrigin->x[i];
+      l = TopReplacement->at(i) - TopOrigin->at(i);
       if (fabs(l) > BondDistance) { // is component greater than bond distance
         Orthovector1[i] = (l < 0) ? -1. : +1.;
 …
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
   ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
+  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
+  BondRescale = TopOrigin->getType()->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
     DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
 …
     case 1:
       FirstOtherAtom = World::getInstance().createAtom();    // new atom
       FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
       FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);  // element is Hydrogen
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
       if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
+      if (TopReplacement->getType()->Z == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
         BondRescale = bondlength;
 …
+      }
       InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
+      FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
+      FirstOtherAtom->x += InBondvector;  // ... and add distance vector to replacement atom
+      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
 …
         // determine the plane of these two with the *origin
         try {
           Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal();
+          Orthovector1 =Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
+        }
         catch(LinearDependenceException &excp){
 …
       FirstOtherAtom = World::getInstance().createAtom();
       SecondOtherAtom = World::getInstance().createAtom();
       FirstOtherAtom->type = elemente->FindElement(1);
       SecondOtherAtom->type = elemente->FindElement(1);
       FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);
+      SecondOtherAtom->setType(1);
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
       SecondOtherAtom->v = TopReplacement->v; // copy velocity
+      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       bondangle = TopOrigin->type->HBondAngle[1];
+      bondangle = TopOrigin->getType()->HBondAngle[1];
       if (bondangle == -1) {
         DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
 …
 //      Log() << Verbose(0) << endl;
 //      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
       FirstOtherAtom->x.Zero();
       SecondOtherAtom->x.Zero();
+      FirstOtherAtom->Zero();
+      SecondOtherAtom->Zero();
       for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
         FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle));
         SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
+      }
       FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
       SecondOtherAtom->x *= BondRescale;
+        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
+        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
+      }
+      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
+      SecondOtherAtom->Scale(BondRescale);
       //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
+      for(int i=NDIM;i--;) { // and make relative to origin atom
+        FirstOtherAtom->x[i] += TopOrigin->x[i];
+        SecondOtherAtom->x[i] += TopOrigin->x[i];
+      }
+      *FirstOtherAtom += TopOrigin->getPosition();
+      *SecondOtherAtom += TopOrigin->getPosition();
       // ... and add to molecule
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 …
       SecondOtherAtom = World::getInstance().createAtom();
       ThirdOtherAtom = World::getInstance().createAtom();
       FirstOtherAtom->type = elemente->FindElement(1);
       SecondOtherAtom->type = elemente->FindElement(1);
       ThirdOtherAtom->type = elemente->FindElement(1);
       FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);
+      SecondOtherAtom->setType(1);
+      ThirdOtherAtom->setType(1);
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
       SecondOtherAtom->v = TopReplacement->v; // copy velocity
+      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
       ThirdOtherAtom->v = TopReplacement->v; // copy velocity
+      ThirdOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
 …
       // create correct coordination for the three atoms
       alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
+      alpha = (TopOrigin->getType()->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
       l = BondRescale;        // desired bond length
       b = 2.*l*sin(alpha);    // base length of isosceles triangle
 …
       factors[1] = f;
       factors[2] = 0.;
       FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g;
       SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[2] = -g;
       ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       // rescale each to correct BondDistance
 …
       // and relative to *origin atom
       FirstOtherAtom->x += TopOrigin->x;
       SecondOtherAtom->x += TopOrigin->x;
       ThirdOtherAtom->x += TopOrigin->x;
+      *FirstOtherAtom += TopOrigin->getPosition();
+      *SecondOtherAtom += TopOrigin->getPosition();
+      *ThirdOtherAtom += TopOrigin->getPosition();
       // ... and add to molecule
 …
     *item >> x[1];
     *item >> x[2];
     Walker->type = elemente->FindElement(shorthand);
     if (Walker->type == NULL) {
+    Walker->setType(elemente->FindElement(shorthand));
+    if (Walker->getType() == NULL) {
       DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
       Walker->type = elemente->FindElement(1);
+      Walker->setType(1);
+    }
     if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
 …
       Walker->Trajectory.F.resize(MDSteps+10);
+    }
+    Walker->setPosition(Vector(x));
     for(j=NDIM;j--;) {
-      Walker->x[j] = x[j];
       Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
       Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
 …
   atom1->RegisterBond(Binder);
   atom2->RegisterBond(Binder);
   if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
+  if ((atom1->getType() != NULL) && (atom1->getType()->Z != 1) && (atom2->getType() != NULL) && (atom2->getType()->Z != 1))
     NoNonBonds++;
 …
   ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
   OBSERVE;
   formula-=pointer->type;
+  formula-=pointer->getType();
   RemoveBonds(pointer);
   erase(pointer);
 …
   if (pointer == NULL)
     return false;
   formula-=pointer->type;
+  formula-=pointer->getType();
   erase(pointer);
   return true;
 …
   for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
     if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
+    if ((*iter)->getType()->Z != 1) // count non-hydrogen atoms whilst at it
       NoNonHydrogen++;
     stringstream sstr;
     sstr << (*iter)->type->symbol << (*iter)->nr+1;
+    sstr << (*iter)->getType()->symbol << (*iter)->nr+1;
     (*iter)->setName(sstr.str());
     DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);

src/molecule.hpp

-              r8d6d31
+              r8f4df1
   bool IsEmpty() const ;
   bool IsEnd() const ;
-  // templates for allowing global manipulation of all vectors
-  template <typename res> void ActOnAllVectors( res (Vector::*f)() ) const;
-  template <typename res> void ActOnAllVectors( res (Vector::*f)() const) const;
-  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t ) const;
-  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t ) const;
-  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T&), T &t ) const;
-  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const;
-  template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const;
-  template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const;
-  template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const;
-  template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const;
   // templates for allowing global manipulation of molecule with each atom as single argument

src/molecule_dynamics.cpp

r8d6d31	r8f4df1
513	513	Sprinter = mol->AddCopyAtom((*iter));
514	514	for (int n=NDIM;n--;) {
515		Sprinter->~~x[n] = (iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(iter)->nr]->Trajectory.R.at(endstep)[n] - (iter)->Trajectory.R.at(startstep)[n])((double)step/(double)MaxSteps~~);
	515	Sprinter->set(n, (iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(iter)->nr]->Trajectory.R.at(endstep)[n] - (iter)->Trajectory.R.at(startstep)[n])((double)step/(double)MaxSteps));
516	516	// add to Trajectories
517	517	//Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;

src/molecule_fragmentation.cpp

-              r8d6d31
+              r8f4df1
       for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     #ifdef ADDHYDROGEN
         if ((*iter)->type->Z != 1) // skip hydrogen
+        if ((*iter)->getType()->Z != 1) // skip hydrogen
     #endif
+        {
 …
     for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
   #ifdef ADDHYDROGEN
       if ((*iter)->type->Z != 1) // skip hydrogen
+      if ((*iter)->getType()->Z != 1) // skip hydrogen
   #endif
+      {
 …
   for (KeySet::iterator runner = Leaf->begin(); runner != Leaf->end(); runner++) {
     Runner = FindAtom((*runner));
     if (Runner->type->Z != 1) { // skip all those added hydrogens when re-filling snake stack
+    if (Runner->getType()->Z != 1) { // skip all those added hydrogens when re-filling snake stack
       if (ShortestPathList[(*runner)] > SP) {  // remove the oldest one with longest shortest path
         SP = ShortestPathList[(*runner)];
 …
     ++iter;
 #ifdef ADDHYDROGEN
     while ((iter != Leaf->end()) && ((*iter)->type->Z == 1)){ // skip added hydrogen
+    while ((iter != Leaf->end()) && ((*iter)->getType()->Z == 1)){ // skip added hydrogen
       iter++;
+    }
 …
         if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
   #ifdef ADDHYDROGEN
          && (OtherWalker->type->Z != 1)
+         && (OtherWalker->getType()->Z != 1)
   #endif
                                                               ) {  // skip hydrogens and restrict to fragment
 …
         Binder = (*BondRunner);
         for (int i=NDIM;i--;) {
           tmp = fabs(Binder->leftatom->x[i] - Binder->rightatom->x[i]);
+          tmp = fabs(Binder->leftatom->at(i) - Binder->rightatom->at(i));
           //Log() << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
           if (tmp > BondDistance) {
 …
       // create translation vector from their periodically modified distance
       for (int i=NDIM;i--;) {
         tmp = Binder->leftatom->x[i] - Binder->rightatom->x[i];
+        tmp = Binder->leftatom->at(i) - Binder->rightatom->at(i);
         if (fabs(tmp) > BondDistance)
           Translationvector[i] = (tmp < 0) ? +1. : -1.;
 …
         //Log() << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
         ColorList[Walker->nr] = black;    // mark as explored
         Walker->x += Translationvector; // translate
+        *Walker += Translationvector; // translate
         for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
           if ((*Runner) != Binder) {

src/molecule_geometry.cpp

-              r8d6d31
+              r8f4df1
   // go through all atoms
+  ActOnAllVectors( &Vector::SubtractVector, *Center);
+  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter -= *Center;
+    *iter -= *CenterBox;
+  }
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 …
   Box &domain = World::getInstance().getDomain();
+  // go through all atoms
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 …
   if (iter != end()) { //list not empty?
     for (int i=NDIM;i--;) {
       max->at(i) = (*iter)->x[i];
       min->at(i) = (*iter)->x[i];
+      max->at(i) = (*iter)->at(i);
+      min->at(i) = (*iter)->at(i);
+    }
     for (; iter != end(); ++iter) {// continue with second if present
       //(*iter)->Output(1,1,out);
       for (int i=NDIM;i--;) {
         max->at(i) = (max->at(i) < (*iter)->x[i]) ? (*iter)->x[i] : max->at(i);
         min->at(i) = (min->at(i) > (*iter)->x[i]) ? (*iter)->x[i] : min->at(i);
+        max->at(i) = (max->at(i) < (*iter)->at(i)) ? (*iter)->at(i) : max->at(i);
+        min->at(i) = (min->at(i) > (*iter)->at(i)) ? (*iter)->at(i) : min->at(i);
+      }
+    }
 …
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
       Center += (*iter)->x;
+      Center += (*iter)->getPosition();
+    }
     Center.Scale(-1./(double)Num); // divide through total number (and sign for direction)
 …
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
       (*a) += (*iter)->x;
+      (*a) += (*iter)->getPosition();
+    }
     a->Scale(1./(double)Num); // divide through total mass (and sign for direction)
 …
   if (iter != end()) {   //list not empty?
     for (; iter != end(); ++iter) {  // continue with second if present
       Num += (*iter)->type->mass;
       tmp = (*iter)->type->mass * (*iter)->x;
+      Num += (*iter)->getType()->mass;
+      tmp = (*iter)->getType()->mass * (*iter)->getPosition();
       (*a) += tmp;
+    }
 …
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
     (*iter)->x.ScaleAll(*factor);
+    (*iter)->ScaleAll(*factor);
+  }
 };
 …
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j) += (*trans);
     (*iter)->x += (*trans);
+    *(*iter) += (*trans);
+  }
 };
 …
   // go through all atoms
+  ActOnAllVectors( &Vector::AddVector, *trans);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter += *trans;
+  }
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 …
     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 #ifdef ADDHYDROGEN
       if ((*iter)->type->Z != 1) {
+      if ((*iter)->getType()->Z != 1) {
 #endif
         Testvector = inversematrix * (*iter)->x;
+        Testvector = inversematrix * (*iter)->getPosition();
         Translationvector.Zero();
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
          if ((*iter)->nr < (*Runner)->GetOtherAtom((*iter))->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
               tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
+              tmp = (*iter)->at(j) - (*Runner)->GetOtherAtom(*iter)->at(j);
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
 …
         // now also change all hydrogens
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
           if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
             Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
+          if ((*Runner)->GetOtherAtom((*iter))->getType()->Z == 1) {
+            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->getPosition();
             Testvector += Translationvector;
             Testvector *= matrix;
 …
   DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     tmp = (*iter)->x[0];
     (*iter)->x[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
     (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
+    tmp = (*iter)->at(0);
+    (*iter)->set(0,  cos(alpha) * tmp + sin(alpha) * (*iter)->at(2));
+    (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
     for (int j=0;j<MDSteps;j++) {
       tmp = (*iter)->Trajectory.R.at(j)[0];
 …
   DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     tmp = (*iter)->x[1];
     (*iter)->x[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
     (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
+    tmp = (*iter)->at(1);
+    (*iter)->set(1,  cos(alpha) * tmp + sin(alpha) * (*iter)->at(2));
+    (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
     for (int j=0;j<MDSteps;j++) {
       tmp = (*iter)->Trajectory.R.at(j)[1];
 …
   // go through all atoms
   for (molecule::const_iterator iter = par->mol->begin(); iter != par->mol->end(); ++iter) {
     if ((*iter)->type == ((struct lsq_params *)params)->type) { // for specific type
       c = (*iter)->x - a;
+    if ((*iter)->getType() == ((struct lsq_params *)params)->type) { // for specific type
+      c = (*iter)->getPosition() - a;
       t = c.ScalarProduct(b);           // get direction parameter
       d = t*b;       // and create vector

src/molecule_graph.cpp

-              r8d6d31
+              r8f4df1
                           //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                           (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
                           const double distance = domain.periodicDistanceSquared(OtherWalker->x,Walker->x);
+                          const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
                           const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
 //                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
 …
         OtherAtom = (*Runner)->GetOtherAtom(Walker);
 #ifdef ADDHYDROGEN
         if (OtherAtom->type->Z != 1) {
+        if (OtherAtom->getType()->Z != 1) {
 #endif
         DoLog(2) && (Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "." << endl);

src/molecule_template.hpp

-              r8d6d31
+              r8f4df1
 /********************************************** declarations *******************************/
-// ================== Acting on all Vectors ========================== //
-// zero arguments
-template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() ) const
+    {
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)();
+  }
-};
-template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() const ) const
+    {
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)();
+  }
-};
-// one argument
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t ) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t);
+  }
-};
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t ) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t);
+  }
-};
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&), T &t ) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t);
+  }
-};
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t);
+  }
-};
-// two arguments
-template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t, u);
+  }
-};
-template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t, u);
+  }
-};
-// three arguments
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t, u, v);
+  }
-};
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const
+{
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    (((*iter)->node)->*f)(t, u, v);
+  }
-};
 // ========================= Summing over each Atoms =================================== //
 …
   int inc = 1;
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     (*Setor) (&array[((*iter)->type->*index)], &inc);
+    (*Setor) (&array[((*iter)->getType()->*index)], &inc);
+  }
 };
 …
+{
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     (*Setor) (&array[((*iter)->type->*index)], &value);
+    (*Setor) (&array[((*iter)->getType()->*index)], &value);
+  }
 };
 …
+{
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     (*Setor) (&array[((*iter)->type->*index)], value);
+    (*Setor) (&array[((*iter)->getType()->*index)], value);
+  }
 };

src/moleculelist.cpp

-              r8d6d31
+              r8f4df1
       size=0.;
       for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
         counts[(*iter)->type->getNumber()]++;
         if ((*iter)->x.DistanceSquared(Origin) > size)
           size = (*iter)->x.DistanceSquared(Origin);
+        counts[(*iter)->getType()->getNumber()]++;
+        if ((*iter)->DistanceSquared(Origin) > size)
+          size = (*iter)->DistanceSquared(Origin);
+      }
       // output Index, Name, number of atoms, chemical formula
 …
   for (molecule::const_iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
     DoLog(2) && (Log() << Verbose(2) << "INFO: Current Walker is " << **iter << "." << endl);
     if (!TesselStruct->IsInnerPoint((*iter)->x, LCList)) {
+    if (!TesselStruct->IsInnerPoint((*iter)->getPosition(), LCList)) {
       CopyAtoms[(*iter)->nr] = (*iter)->clone();
       mol->AddAtom(CopyAtoms[(*iter)->nr]);
 …
     for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
       //Log() << Verbose(1) << "(*iter): " << *(*iter) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
       if (((*iter)->type->Z == 1) && (((*iter)->father == NULL)
           || ((*iter)->father->type->Z != 1))) { // if it's a hydrogen
+      if (((*iter)->getType()->Z == 1) && (((*iter)->father == NULL)
+          || ((*iter)->father->getType()->Z != 1))) { // if it's a hydrogen
         for (molecule::const_iterator runner = (*ListRunner)->begin(); runner != (*ListRunner)->end(); ++runner) {
           //Log() << Verbose(2) << "Runner: " << *(*runner) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
           Binder = *((*runner)->ListOfBonds.begin());
           if (((*runner)->type->Z == 1) && ((*runner)->nr > (*iter)->nr) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
+          if (((*runner)->getType()->Z == 1) && ((*runner)->nr > (*iter)->nr) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
             // 4. evaluate the morse potential for each matrix component and add up
             distance = (*runner)->x.distance((*iter)->x);
+            distance = (*runner)->distance(*(*iter));
             //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *(*runner) << "<= " << distance << "=>" << *(*iter) << ":" << endl;
             for (int k = 0; k < a; k++) {
 …
         if ((*ListRunner)->hasElement((*elemIter).first)) { // if this element got atoms
           for(molecule::iterator atomIter = (*ListRunner)->begin(); atomIter !=(*ListRunner)->end();++atomIter){
             if ((*atomIter)->type->getNumber() == (*elemIter).first) {
+            if ((*atomIter)->getType()->getNumber() == (*elemIter).first) {
               if (((*atomIter)->GetTrueFather() != NULL) && ((*atomIter)->GetTrueFather() != (*atomIter))) {// if there is a rea
                 //Log() << Verbose(0) << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
 …
         if (AtomMask[Father->nr]) // apply mask
 #ifdef ADDHYDROGEN
           if ((*iter)->type->Z != 1) // skip hydrogen
+          if ((*iter)->getType()->Z != 1) // skip hydrogen
 #endif
           RootStack[FragmentCounter].push_front((*iter)->nr);

src/tesselation.cpp

-              r8d6d31
+              r8f4df1
 #include "Helpers/helpers.hpp"
 #include "Helpers/Info.hpp"
+#include "BoundaryPointSet.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "BoundaryPolygonSet.hpp"
+#include "TesselPoint.hpp"
+#include "CandidateForTesselation.hpp"
+#include "PointCloud.hpp"
 #include "linkedcell.hpp"
 #include "Helpers/Log.hpp"
 …
 #include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/Line.hpp"
 #include "vector_ops.hpp"
+#include "LinearAlgebra/vector_ops.hpp"
 #include "Helpers/Verbose.hpp"
 #include "LinearAlgebra/Plane.hpp"
 …
 class molecule;
-// ======================================== Points on Boundary =================================
-/** Constructor of BoundaryPointSet.
- */
-BoundaryPointSet::BoundaryPointSet() :
-  LinesCount(0), value(0.), Nr(-1)
+{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding noname." << endl);
+}
+;
-/** Constructor of BoundaryPointSet with Tesselpoint.
- * \param *Walker TesselPoint this boundary point represents
- */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
-  LinesCount(0), node(Walker), value(0.), Nr(Walker->nr)
+{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding Node " << *Walker << endl);
+}
+;
-/** Destructor of BoundaryPointSet.
- * Sets node to NULL to avoid removing the original, represented TesselPoint.
- * \note When removing point from a class Tesselation, use RemoveTesselationPoint()
- */
-BoundaryPointSet::~BoundaryPointSet()
+{
-  Info FunctionInfo(__func__);
-  //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
-  if (!lines.empty())
-    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
-  node = NULL;
+}
+;
-/** Add a line to the LineMap of this point.
- * \param *line line to add
- */
-void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." << endl);
-  if (line->endpoints[0] == this) {
-    lines.insert(LinePair(line->endpoints[1]->Nr, line));
-  } else {
-    lines.insert(LinePair(line->endpoints[0]->Nr, line));
+  }
-  LinesCount++;
+}
+;
-/** output operator for BoundaryPointSet.
- * \param &ost output stream
- * \param &a boundary point
- */
-ostream & operator <<(ostream &ost, const BoundaryPointSet &a)
+{
-  ost << "[" << a.Nr << "|" << a.node->getName() << " at " << *a.node->node << "]";
-  return ost;
+}
+;
-// ======================================== Lines on Boundary =================================
-/** Constructor of BoundaryLineSet.
- */
-BoundaryLineSet::BoundaryLineSet() :
-  Nr(-1)
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 2; i++)
-    endpoints[i] = NULL;
+}
+;
-/** Constructor of BoundaryLineSet with two endpoints.
- * Adds line automatically to each endpoints' LineMap
- * \param *Point[2] array of two boundary points
- * \param number number of the list
- */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
+{
-  Info FunctionInfo(__func__);
-  // set number
-  Nr = number;
-  // set endpoints in ascending order
-  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
-  // add this line to the hash maps of both endpoints
-  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
-  DoLog(0) && (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;
-  // clear triangles list
-  DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl);
+}
+;
-/** Destructor for BoundaryLineSet.
- * Removes itself from each endpoints' LineMap, calling RemoveTrianglePoint() when point not connected anymore.
- * \note When removing lines from a class Tesselation, use RemoveTesselationLine()
- */
-BoundaryLineSet::~BoundaryLineSet()
+{
-  Info FunctionInfo(__func__);
-  int Numbers[2];
-  // get other endpoint number of finding copies of same line
-  if (endpoints[1] != NULL)
-    Numbers[0] = endpoints[1]->Nr;
-  else
-    Numbers[0] = -1;
-  if (endpoints[0] != NULL)
-    Numbers[1] = endpoints[0]->Nr;
-  else
-    Numbers[1] = -1;
-  for (int i = 0; i < 2; i++) {
-    if (endpoints[i] != NULL) {
-      if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
-        pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
-        for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
-          if ((*Runner).second == this) {
-            //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-            endpoints[i]->lines.erase(Runner);
-            break;
+          }
-      } else { // there's just a single line left
-        if (endpoints[i]->lines.erase(Nr)) {
-          //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        }
+      }
-      if (endpoints[i]->lines.empty()) {
-        //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-        if (endpoints[i] != NULL) {
-          delete (endpoints[i]);
-          endpoints[i] = NULL;
+        }
+      }
+    }
+  }
-  if (!triangles.empty())
-    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
+}
+;
-/** Add triangle to TriangleMap of this boundary line.
- * \param *triangle to add
- */
-void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+{
-  Info FunctionInfo(__func__);
-  DoLog(0) && (Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl);
-  triangles.insert(TrianglePair(triangle->Nr, triangle));
+}
+;
-/** Checks whether we have a common endpoint with given \a *line.
- * \param *line other line to test
- * \return true - common endpoint present, false - not connected
- */
-bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
+{
-  Info FunctionInfo(__func__);
-  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
-    return true;
-  else
-    return false;
+}
+;
-/** Checks whether the adjacent triangles of a baseline are convex or not.
- * We sum the two angles of each height vector with respect to the center of the baseline.
- * If greater/equal M_PI than we are convex.
- * \param *out output stream for debugging
- * \return true - triangles are convex, false - concave or less than two triangles connected
- */
-bool BoundaryLineSet::CheckConvexityCriterion() const
+{
-  Info FunctionInfo(__func__);
-  double angle = CalculateConvexity();
-  if (angle > -MYEPSILON) {
-    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
-    return true;
-  } else {
-    DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
-    return false;
+  }
+}
-/** Calculates the angle between two triangles with respect to their normal vector.
- * We sum the two angles of each height vector with respect to the center of the baseline.
- * \return angle > 0 then convex, if < 0 then concave
- */
-double BoundaryLineSet::CalculateConvexity() const
+{
-  Info FunctionInfo(__func__);
-  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
-  // get the two triangles
-  if (triangles.size() != 2) {
-    DoeLog(0) && (eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
-    return true;
+  }
-  // check normal vectors
-  // have a normal vector on the base line pointing outwards
-  //Log() << Verbose(0) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
-  BaseLineCenter = (1./2.)*((*endpoints[0]->node->node) + (*endpoints[1]->node->node));
-  BaseLine = (*endpoints[0]->node->node) - (*endpoints[1]->node->node);
-  //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
-  BaseLineNormal.Zero();
-  NormalCheck.Zero();
-  double sign = -1.;
-  int i = 0;
-  class BoundaryPointSet *node = NULL;
-  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 += runner->second->NormalVector;
-    NormalCheck *= sign;
-    sign = -sign;
-    if (runner->second->NormalVector.NormSquared() > MYEPSILON)
-      BaseLineNormal = runner->second->NormalVector;   // yes, copy second on top of first
-    else {
-      DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
+    }
-    node = runner->second->GetThirdEndpoint(this);
-    if (node != NULL) {
-      //Log() << Verbose(0) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
-      helper[i] = (*node->node->node) - BaseLineCenter;
-      helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
-      //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
-      i++;
-    } else {
-      DoeLog(1) && (eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
-      return true;
+    }
+  }
-  //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
-  if (NormalCheck.NormSquared() < MYEPSILON) {
-    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl);
-    return true;
+  }
-  BaseLineNormal.Scale(-1.);
-  double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
-  return (angle - M_PI);
+}
-/** Checks whether point is any of the two endpoints this line contains.
- * \param *point point to test
- * \return true - point is of the line, false - is not
- */
-bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 2; i++)
-    if (point == endpoints[i])
-      return true;
-  return false;
+}
+;
-/** Returns other endpoint of the line.
- * \param *point other endpoint
- * \return NULL - if endpoint not contained in BoundaryLineSet::lines, or pointer to BoundaryPointSet otherwise
- */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
+{
-  Info FunctionInfo(__func__);
-  if (endpoints[0] == point)
-    return endpoints[1];
-  else if (endpoints[1] == point)
-    return endpoints[0];
-  else
-    return NULL;
+}
+;
-/** Returns other triangle of the line.
- * \param *point other endpoint
- * \return NULL - if triangle not contained in BoundaryLineSet::triangles, or pointer to BoundaryTriangleSet otherwise
- */
-class BoundaryTriangleSet *BoundaryLineSet::GetOtherTriangle(const BoundaryTriangleSet * const triangle) const
+{
-  Info FunctionInfo(__func__);
-  if (triangles.size() == 2) {
-    for (TriangleMap::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); ++TriangleRunner)
-      if (TriangleRunner->second != triangle)
-        return TriangleRunner->second;
+  }
-  return NULL;
+}
+;
-/** output operator for BoundaryLineSet.
- * \param &ost output stream
- * \param &a boundary line
- */
-ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
+{
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->getName() << " at " << *a.endpoints[1]->node->node << "]";
-  return ost;
+}
+;
-// ======================================== Triangles on Boundary =================================
-/** Constructor for BoundaryTriangleSet.
- */
-BoundaryTriangleSet::BoundaryTriangleSet() :
-  Nr(-1)
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++) {
-    endpoints[i] = NULL;
-    lines[i] = NULL;
+  }
+}
+;
-/** Constructor for BoundaryTriangleSet with three lines.
- * \param *line[3] lines that make up the triangle
- * \param number number of triangle
- */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
-  Nr(number)
+{
-  Info FunctionInfo(__func__);
-  // set number
-  // set lines
-  for (int i = 0; i < 3; i++) {
-    lines[i] = line[i];
-    lines[i]->AddTriangle(this);
+  }
-  // get ascending order of endpoints
-  PointMap OrderMap;
-  for (int i = 0; i < 3; i++) {
-    // for all three lines
-    for (int j = 0; j < 2; j++) { // for both endpoints
-      OrderMap.insert(pair<int, class BoundaryPointSet *> (line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
-      // and we don't care whether insertion fails
+    }
+  }
-  // set endpoints
-  int Counter = 0;
-  DoLog(0) && (Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl);
-  for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
-    endpoints[Counter] = runner->second;
-    DoLog(0) && (Log() << Verbose(0) << " " << *endpoints[Counter] << endl);
-    Counter++;
+  }
-  ASSERT(Counter >= 3,"We have a triangle with only two distinct endpoints!");
-};
-/** Destructor of BoundaryTriangleSet.
- * Removes itself from each of its lines' LineMap and removes them if necessary.
- * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
- */
-BoundaryTriangleSet::~BoundaryTriangleSet()
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++) {
-    if (lines[i] != NULL) {
-      if (lines[i]->triangles.erase(Nr)) {
-        //Log() << Verbose(0) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      }
-      if (lines[i]->triangles.empty()) {
-        //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-        delete (lines[i]);
-        lines[i] = NULL;
+      }
+    }
+  }
-  //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl;
+}
+;
-/** Calculates the normal vector for this triangle.
- * Is made unique by comparison with \a OtherVector to point in the other direction.
- * \param &OtherVector direction vector to make normal vector unique.
- */
-void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+{
-  Info FunctionInfo(__func__);
-  // get normal vector
-  NormalVector = Plane(*(endpoints[0]->node->node),
-                       *(endpoints[1]->node->node),
-                       *(endpoints[2]->node->node)).getNormal();
-  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (NormalVector.ScalarProduct(OtherVector) > 0.)
-    NormalVector.Scale(-1.);
-  DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl);
+}
+;
-/** 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 *out output stream for debugging
- * \param *MolCenter offset vector of line
- * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
- * \param *Intersection intersection on plane on return
- * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
- */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+{
-  Info FunctionInfo(__func__);
-  Vector CrossPoint;
-  Vector helper;
-  try {
-    Line centerLine = makeLineThrough(*MolCenter, *x);
-    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
-    if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
-      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
-      return true;
-    }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
-      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
-      return true;
-    }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
-      DoLog(1) && (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;
-    do {
-      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
-      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
-      CrossPoint = line1.getIntersection(line2);
-      helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
-      CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-      const double s = CrossPoint.ScalarProduct(helper)/helper.NormSquared();
-      DoLog(1) && (Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl);
-      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
-        DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
-        return false;
+      }
-      i++;
-    } while (i < 3);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
-    return true;
+  }
-  catch (MathException &excp) {
-    Log() << Verbose(1) << excp;
-    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
-    return false;
+  }
+}
+;
-/** Finds the point on the triangle to the point \a *x.
- * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point.
- * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the
- * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down.
- * 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
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl);
-  GetCenter(&Direction);
-  try {
-    Line l = makeLineThrough(*x, Direction);
-    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
+  }
-  catch (MathException &excp) {
-    (*ClosestPoint) = (*x);
+  }
-  // 2. Calculate in plane part of line (x, intersection)
-  Vector InPlane = (*x) - (*ClosestPoint); // points from plane intersection to straight-down point
-  InPlane.ProjectOntoPlane(NormalVector);
-  InPlane += *ClosestPoint;
-  DoLog(2) && (Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl);
-  DoLog(2) && (Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl);
-  DoLog(2) && (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 = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
-    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
-    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
-    CrossDirection[i] = CrossPoint[i] - InPlane;
-    CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-    const double s = CrossPoint[i].ScalarProduct(Direction)/Direction.NormSquared();
-    DoLog(2) && (Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl);
-    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
-          CrossPoint[i] += (*endpoints[i%3]->node->node);  // make cross point absolute again
-      DoLog(2) && (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) = 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) = 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) = (*endpoints[i]->node->node);
+      }
+    }
+  }
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl);
-  return ShortestDistance;
+}
+;
-/** Checks whether lines is any of the three boundary lines this triangle contains.
- * \param *line line to test
- * \return true - line is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
-    if (line == lines[i])
-      return true;
-  return false;
+}
+;
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point point to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
-    if (point == endpoints[i])
-      return true;
-  return false;
+}
+;
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point TesselPoint to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
-    if (point == endpoints[i]->node)
-      return true;
-  return false;
+}
+;
-/** Checks whether three given \a *Points coincide with triangle's endpoints.
- * \param *Points[3] pointer to BoundaryPointSet
- * \return true - is the very triangle, false - is not
- */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (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]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])
-  ));
+}
+;
-/** Checks whether three given \a *Points coincide with triangle's endpoints.
- * \param *Points[3] pointer to BoundaryPointSet
- * \return true - is the very triangle, false - is not
- */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
+{
-  Info FunctionInfo(__func__);
-  return (((endpoints[0] == T->endpoints[0]) || (endpoints[0] == T->endpoints[1]) || (endpoints[0] == T->endpoints[2])) && ((endpoints[1] == T->endpoints[0]) || (endpoints[1] == T->endpoints[1]) || (endpoints[1] == T->endpoints[2])) && ((endpoints[2] == T->endpoints[0]) || (endpoints[2] == T->endpoints[1]) || (endpoints[2] == T->endpoints[2])
-  ));
+}
+;
-/** Returns the endpoint which is not contained in the given \a *line.
- * \param *line baseline defining two endpoints
- * \return pointer third endpoint or NULL if line does not belong to triangle.
- */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
+{
-  Info FunctionInfo(__func__);
-  // sanity check
-  if (!ContainsBoundaryLine(line))
-    return NULL;
-  for (int i = 0; i < 3; i++)
-    if (!line->ContainsBoundaryPoint(endpoints[i]))
-      return endpoints[i];
-  // actually, that' impossible :)
-  return NULL;
+}
+;
-/** Returns the baseline which does not contain the given boundary point \a *point.
- * \param *point endpoint which is neither endpoint of the desired line
- * \return pointer to desired third baseline
- */
-class BoundaryLineSet *BoundaryTriangleSet::GetThirdLine(const BoundaryPointSet * const point) const
+{
-  Info FunctionInfo(__func__);
-  // sanity check
-  if (!ContainsBoundaryPoint(point))
-    return NULL;
-  for (int i = 0; i < 3; i++)
-    if (!lines[i]->ContainsBoundaryPoint(point))
-      return lines[i];
-  // actually, that' impossible :)
-  return NULL;
+}
+;
-/** Calculates the center point of the triangle.
- * Is third of the sum of all endpoints.
- * \param *center central point on return.
- */
-void BoundaryTriangleSet::GetCenter(Vector * const center) const
+{
-  Info FunctionInfo(__func__);
-  center->Zero();
-  for (int i = 0; i < 3; i++)
-    (*center) += (*endpoints[i]->node->node);
-  center->Scale(1. / 3.);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl);
+}
-/**
- * gets the Plane defined by the three triangle Basepoints
- */
-Plane BoundaryTriangleSet::getPlane() const{
-  ASSERT(endpoints[0] && endpoints[1] && endpoints[2], "Triangle not fully defined");
-  return Plane(*endpoints[0]->node->node,
-               *endpoints[1]->node->node,
-               *endpoints[2]->node->node);
+}
-Vector BoundaryTriangleSet::getEndpoint(int i) const{
-  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
-  return *endpoints[i]->node->node;
+}
-string BoundaryTriangleSet::getEndpointName(int i) const{
-  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
-  return endpoints[i]->node->getName();
+}
-/** output operator for BoundaryTriangleSet.
- * \param &ost output stream
- * \param &a boundary triangle
- */
-ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a)
+{
-  ost << "[" << a.Nr << "|" << a.getEndpointName(0) << "," << a.getEndpointName(1) << "," << a.getEndpointName(2) << "]";
-  //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
-  //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
-  return ost;
+}
+;
-// ======================================== Polygons on Boundary =================================
-/** Constructor for BoundaryPolygonSet.
- */
-BoundaryPolygonSet::BoundaryPolygonSet() :
-  Nr(-1)
+{
-  Info FunctionInfo(__func__);
+}
+;
-/** Destructor of BoundaryPolygonSet.
- * Just clears endpoints.
- * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
- */
-BoundaryPolygonSet::~BoundaryPolygonSet()
+{
-  Info FunctionInfo(__func__);
-  endpoints.clear();
-  DoLog(1) && (Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl);
+}
+;
-/** Calculates the normal vector for this triangle.
- * Is made unique by comparison with \a OtherVector to point in the other direction.
- * \param &OtherVector direction vector to make normal vector unique.
- * \return allocated vector in normal direction
- */
-Vector * BoundaryPolygonSet::GetNormalVector(const Vector &OtherVector) const
+{
-  Info FunctionInfo(__func__);
-  // get normal vector
-  Vector TemporaryNormal;
-  Vector *TotalNormal = new Vector;
-  PointSet::const_iterator Runner[3];
-  for (int i = 0; i < 3; i++) {
-    Runner[i] = endpoints.begin();
-    for (int j = 0; j < i; j++) { // go as much further
-      Runner[i]++;
-      if (Runner[i] == endpoints.end()) {
-        DoeLog(0) && (eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
-        performCriticalExit();
+      }
+    }
+  }
-  TotalNormal->Zero();
-  int counter = 0;
-  for (; Runner[2] != endpoints.end();) {
-    TemporaryNormal = Plane(*((*Runner[0])->node->node),
-                            *((*Runner[1])->node->node),
-                            *((*Runner[2])->node->node)).getNormal();
-    for (int i = 0; i < 3; i++) // increase each of them
-      Runner[i]++;
-    (*TotalNormal) += TemporaryNormal;
+  }
-  TotalNormal->Scale(1. / (double) counter);
-  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (TotalNormal->ScalarProduct(OtherVector) > 0.)
-    TotalNormal->Scale(-1.);
-  DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl);
-  return TotalNormal;
+}
+;
-/** Calculates the center point of the triangle.
- * Is third of the sum of all endpoints.
- * \param *center central point on return.
- */
-void BoundaryPolygonSet::GetCenter(Vector * const center) const
+{
-  Info FunctionInfo(__func__);
-  center->Zero();
-  int counter = 0;
-  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    (*center) += (*(*Runner)->node->node);
-    counter++;
+  }
-  center->Scale(1. / (double) counter);
-  DoLog(1) && (Log() << Verbose(1) << "Center is at " << *center << "." << endl);
+}
-/** Checks whether the polygons contains all three endpoints of the triangle.
- * \param *triangle triangle to test
- * \return true - triangle is contained polygon, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryTriangle(const BoundaryTriangleSet * const triangle) const
+{
-  Info FunctionInfo(__func__);
-  return ContainsPresentTupel(triangle->endpoints, 3);
+}
+;
-/** Checks whether the polygons contains both endpoints of the line.
- * \param *line line to test
- * \return true - line is of the triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+{
-  Info FunctionInfo(__func__);
-  return ContainsPresentTupel(line->endpoints, 2);
+}
+;
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point point to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+{
-  Info FunctionInfo(__func__);
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    DoLog(0) && (Log() << Verbose(0) << "Checking against " << **Runner << endl);
-    if (point == (*Runner)) {
-      DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
-      return true;
+    }
+  }
-  DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
-  return false;
+}
+;
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point TesselPoint to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+{
-  Info FunctionInfo(__func__);
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
-    if (point == (*Runner)->node) {
-      DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
-      return true;
+    }
-  DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
-  return false;
+}
+;
-/** Checks whether given array of \a *Points coincide with polygons's endpoints.
- * \param **Points pointer to an array of BoundaryPointSet
- * \param dim dimension of array
- * \return true - set of points is contained in polygon, false - is not
- */
-bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const
+{
-  Info FunctionInfo(__func__);
-  int counter = 0;
-  DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
-  for (int i = 0; i < dim; i++) {
-    DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl);
-    if (ContainsBoundaryPoint(Points[i])) {
-      counter++;
+    }
+  }
-  if (counter == dim)
-    return true;
-  else
-    return false;
+}
+;
-/** Checks whether given PointList coincide with polygons's endpoints.
- * \param &endpoints PointList
- * \return true - set of points is contained in polygon, false - is not
- */
-bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const
+{
-  Info FunctionInfo(__func__);
-  size_t counter = 0;
-  DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << **Runner << endl);
-    if (ContainsBoundaryPoint(*Runner))
-      counter++;
+  }
-  if (counter == endpoints.size())
-    return true;
-  else
-    return false;
+}
+;
-/** Checks whether given set of \a *Points coincide with polygons's endpoints.
- * \param *P pointer to BoundaryPolygonSet
- * \return true - is the very triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
+{
-  return ContainsPresentTupel((const PointSet) P->endpoints);
+}
+;
-/** Gathers all the endpoints' triangles in a unique set.
- * \return set of all triangles
- */
-TriangleSet * BoundaryPolygonSet::GetAllContainedTrianglesFromEndpoints() const
+{
-  Info FunctionInfo(__func__);
-  pair<TriangleSet::iterator, bool> Tester;
-  TriangleSet *triangles = new TriangleSet;
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
-    for (LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)
-      for (TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {
-        //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl;
-        if (ContainsBoundaryTriangle(Sprinter->second)) {
-          Tester = triangles->insert(Sprinter->second);
-          if (Tester.second)
-            DoLog(0) && (Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl);
+        }
+      }
-  DoLog(1) && (Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl);
-  return triangles;
+}
+;
-/** Fills the endpoints of this polygon from the triangles attached to \a *line.
- * \param *line lines with triangles attached
- * \return true - polygon contains endpoints, false - line was NULL
- */
-bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)
+{
-  Info FunctionInfo(__func__);
-  pair<PointSet::iterator, bool> Tester;
-  if (line == NULL)
-    return false;
-  DoLog(1) && (Log() << Verbose(1) << "Filling polygon from line " << *line << endl);
-  for (TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {
-    for (int i = 0; i < 3; i++) {
-      Tester = endpoints.insert((Runner->second)->endpoints[i]);
-      if (Tester.second)
-        DoLog(1) && (Log() << Verbose(1) << "  Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl);
+    }
+  }
-  return true;
+}
+;
-/** output operator for BoundaryPolygonSet.
- * \param &ost output stream
- * \param &a boundary polygon
- */
-ostream &operator <<(ostream &ost, const BoundaryPolygonSet &a)
+{
-  ost << "[" << a.Nr << "|";
-  for (PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {
-    ost << (*Runner)->node->getName();
-    Runner++;
-    if (Runner != a.endpoints.end())
-      ost << ",";
+  }
-  ost << "]";
-  return ost;
+}
+;
-// =========================================================== class TESSELPOINT ===========================================
-/** Constructor of class TesselPoint.
- */
-TesselPoint::TesselPoint()
+{
-  //Info FunctionInfo(__func__);
-  node = NULL;
-  nr = -1;
+}
+;
-/** Destructor for class TesselPoint.
- */
-TesselPoint::~TesselPoint()
+{
-  //Info FunctionInfo(__func__);
+}
+;
-/** Prints LCNode to screen.
- */
-ostream & operator <<(ostream &ost, const TesselPoint &a)
+{
-  ost << "[" << a.getName() << "|" << *a.node << "]";
-  return ost;
+}
+;
-/** Prints LCNode to screen.
- */
-ostream & TesselPoint::operator <<(ostream &ost)
+{
-  Info FunctionInfo(__func__);
-  ost << "[" << (nr) << "|" << this << "]";
-  return ost;
+}
+;
-// =========================================================== class POINTCLOUD ============================================
-/** Constructor of class PointCloud.
- */
-PointCloud::PointCloud()
+{
-  //Info FunctionInfo(__func__);
+}
+;
-/** Destructor for class PointCloud.
- */
-PointCloud::~PointCloud()
+{
-  //Info FunctionInfo(__func__);
+}
+;
-// ============================ CandidateForTesselation =============================
-/** Constructor of class CandidateForTesselation.
- */
-CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) :
-  BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
+{
-  Info FunctionInfo(__func__);
+}
+;
-/** Constructor of class CandidateForTesselation.
- */
-CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
-  BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
+{
-        Info FunctionInfo(__func__);
-  OptCenter = OptCandidateCenter;
-  OtherOptCenter = OtherOptCandidateCenter;
-};
-/** Destructor for class CandidateForTesselation.
- */
-CandidateForTesselation::~CandidateForTesselation()
+{
+}
+;
-/** Checks validity of a given sphere of a candidate line.
- * Sphere must touch all candidates and the baseline endpoints and there must be no other atoms inside.
- * \param RADIUS radius of sphere
- * \param *LC LinkedCell structure with other atoms
- * \return true - sphere is valid, false - sphere contains other points
- */
-bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell *LC) const
+{
-  Info FunctionInfo(__func__);
-  const double radiusSquared = RADIUS * RADIUS;
-  list<const Vector *> VectorList;
-  VectorList.push_back(&OptCenter);
-  //VectorList.push_back(&OtherOptCenter);  // don't check the other (wrong) center
-  if (!pointlist.empty())
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
-  else
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
-  // check baseline for OptCenter and OtherOptCenter being on sphere's surface
-  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-    for (int i = 0; i < 2; i++) {
-      const double distance = fabs((*VRunner)->DistanceSquared(*BaseLine->endpoints[i]->node->node) - radiusSquared);
-      if (distance > HULLEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
-        return false;
+      }
+    }
+  }
-  // check Candidates for OptCenter and OtherOptCenter being on sphere's surface
-  for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
-    const TesselPoint *Walker = *Runner;
-    for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-      const double distance = fabs((*VRunner)->DistanceSquared(*Walker->node) - radiusSquared);
-      if (distance > HULLEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
-        return false;
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl);
+      }
+    }
+  }
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
-  bool flag = true;
-  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-    // get all points inside the sphere
-    TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));
-    DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << (*VRunner) << ":" << endl);
-    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(*(*VRunner)) << "." << endl);
-    // remove baseline's endpoints and candidates
-    for (int i = 0; i < 2; i++) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl);
-      ListofPoints->remove(BaseLine->endpoints[i]->node);
+    }
-    for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl);
-      ListofPoints->remove(*Runner);
+    }
-    if (!ListofPoints->empty()) {
-      DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
-      flag = false;
-      DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
-      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-        DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << " at distance " << setprecision(13) << (*Runner)->node->distance(*(*VRunner)) << setprecision(6) << "." << endl);
-      // check with animate_sphere.tcl VMD script
-      if (ThirdPoint != NULL) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
-      } else {
-        DoeLog(1) && (eLog() << Verbose(1) << "Check by: ... missing third point ..." << endl);
-        DoeLog(1) && (eLog() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
+      }
+    }
-    delete (ListofPoints);
+  }
-  return flag;
+}
+;
-/** output operator for CandidateForTesselation.
- * \param &ost output stream
- * \param &a boundary line
- */
-ostream & operator <<(ostream &ost, const CandidateForTesselation &a)
+{
-  ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->getName() << "," << a.BaseLine->endpoints[1]->node->getName() << "] with ";
-  if (a.pointlist.empty())
-    ost << "no candidate.";
-  else {
-    ost << "candidate";
-    if (a.pointlist.size() != 1)
-      ost << "s ";
-    else
-      ost << " ";
-    for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)
-      ost << *(*Runner) << " ";
-    ost << " at angle " << (a.ShortestAngle) << ".";
+  }
-  return ost;
+}
+;
-// =========================================================== class TESSELATION ===========================================
 /** Constructor of class Tesselation.
 …
   int num = 0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
     (*Center) += (*GetPoint()->node);
+    (*Center) += (GetPoint()->getPosition());
     num++;
+  }
 …
       C++;
       for (; C != PointsOnBoundary.end(); C++) {
         tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
+        tmp = A->second->node->DistanceSquared(B->second->node->getPosition());
         distance = tmp * tmp;
         tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
+        tmp = A->second->node->DistanceSquared(C->second->node->getPosition());
         distance += tmp * tmp;
         tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
+        tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
         distance += tmp * tmp;
         DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
 …
     // 2. next, we have to check whether all points reside on only one side of the triangle
     // 3. construct plane vector
     PlaneVector = Plane(*A->second->node->node,
                         *baseline->second.first->second->node->node,
                         *baseline->second.second->second->node->node).getNormal();
+    PlaneVector = Plane(A->second->node->getPosition(),
+                        baseline->second.first->second->node->getPosition(),
+                        baseline->second.second->second->node->getPosition()).getNormal();
     DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
     // 4. loop over all points
 …
         continue;
       // 4a. project onto plane vector
+      TrialVector = (*checker->second->node->node);
+      TrialVector.SubtractVector(*A->second->node->node);
+      TrialVector = (checker->second->node->getPosition() - A->second->node->getPosition());
       distance = TrialVector.ScalarProduct(PlaneVector);
       if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
 …
+      }
       // 4d. Check whether the point is inside the triangle (check distance to each node
       tmp = checker->second->node->node->DistanceSquared(*A->second->node->node);
+      tmp = checker->second->node->DistanceSquared(A->second->node->getPosition());
       int innerpoint = 0;
       if ((tmp < A->second->node->node->DistanceSquared(*baseline->second.first->second->node->node)) && (tmp < A->second->node->node->DistanceSquared(*baseline->second.second->second->node->node)))
+      if ((tmp < A->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < A->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
         innerpoint++;
       tmp = checker->second->node->node->DistanceSquared(*baseline->second.first->second->node->node);
       if ((tmp < baseline->second.first->second->node->node->DistanceSquared(*A->second->node->node)) && (tmp < baseline->second.first->second->node->node->DistanceSquared(*baseline->second.second->second->node->node)))
+      tmp = checker->second->node->DistanceSquared(baseline->second.first->second->node->getPosition());
+      if ((tmp < baseline->second.first->second->node->DistanceSquared(A->second->node->getPosition())) && (tmp < baseline->second.first->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
         innerpoint++;
       tmp = checker->second->node->node->DistanceSquared(*baseline->second.second->second->node->node);
       if ((tmp < baseline->second.second->second->node->node->DistanceSquared(*baseline->second.first->second->node->node)) && (tmp < baseline->second.second->second->node->node->DistanceSquared(*A->second->node->node)))
+      tmp = checker->second->node->DistanceSquared(baseline->second.second->second->node->getPosition());
+      if ((tmp < baseline->second.second->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < baseline->second.second->second->node->DistanceSquared(A->second->node->getPosition())))
         innerpoint++;
       // 4e. If so, break 4. loop and continue with next candidate in 1. loop
 …
         // prepare some auxiliary vectors
         Vector BaseLineCenter, BaseLine;
         BaseLineCenter = 0.5 * ((*baseline->second->endpoints[0]->node->node) +
                                 (*baseline->second->endpoints[1]->node->node));
         BaseLine = (*baseline->second->endpoints[0]->node->node) - (*baseline->second->endpoints[1]->node->node);
+        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                (baseline->second->endpoints[1]->node->getPosition()));
+        BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition());
         // offset to center of triangle
 …
         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
         PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
         TempVector = CenterVector - (*baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+        TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
         //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
         if (PropagationVector.ScalarProduct(TempVector) > 0) // make sure normal propagation vector points outward from baseline
 …
             // first check direction, so that triangles don't intersect
             VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
+            VirtualNormalVector = (target->second->node->getPosition()) - BaseLineCenter;
             VirtualNormalVector.ProjectOntoPlane(NormalVector);
             TempAngle = VirtualNormalVector.Angle(PropagationVector);
 …
             // check for linear dependence
             TempVector = (*baseline->second->endpoints[0]->node->node) - (*target->second->node->node);
             helper = (*baseline->second->endpoints[1]->node->node) - (*target->second->node->node);
+            TempVector = (baseline->second->endpoints[0]->node->getPosition()) - (target->second->node->getPosition());
+            helper = (baseline->second->endpoints[1]->node->getPosition()) - (target->second->node->getPosition());
             helper.ProjectOntoPlane(TempVector);
             if (fabs(helper.NormSquared()) < MYEPSILON) {
 …
             // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
             flag = true;
             VirtualNormalVector = Plane(*(baseline->second->endpoints[0]->node->node),
                                         *(baseline->second->endpoints[1]->node->node),
                                         *(target->second->node->node)).getNormal();
             TempVector = (1./3.) * ((*baseline->second->endpoints[0]->node->node) +
                                     (*baseline->second->endpoints[1]->node->node) +
                                     (*target->second->node->node));
+            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()),
+                                        (baseline->second->endpoints[1]->node->getPosition()),
+                                        (target->second->node->getPosition())).getNormal();
+            TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                    (baseline->second->endpoints[1]->node->getPosition()) +
+                                    (target->second->node->getPosition()));
             TempVector -= (*Center);
             // make it always point outward
 …
             } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
               // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
               helper = (*target->second->node->node) - BaseLineCenter;
+              helper = (target->second->node->getPosition()) - BaseLineCenter;
               helper.ProjectOntoPlane(BaseLine);
               // ...the one with the smaller angle is the better candidate
               TempVector = (*target->second->node->node) - BaseLineCenter;
+              TempVector = (target->second->node->getPosition()) - BaseLineCenter;
               TempVector.ProjectOntoPlane(VirtualNormalVector);
               TempAngle = TempVector.Angle(helper);
               TempVector = (*winner->second->node->node) - BaseLineCenter;
+              TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
               TempVector.ProjectOntoPlane(VirtualNormalVector);
               if (TempAngle < TempVector.Angle(helper)) {
 …
             BLS[2] = LineChecker[1]->second;
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
           BTS->GetCenter(&helper);
+          BTS->GetCenter(helper);
           helper -= (*Center);
           helper *= -1;
 …
     DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
     // get the next triangle
     triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->getPosition(), BoundaryPoints);
     if (triangles != NULL)
       BTS = triangles->front();
 …
     DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
     // get the intersection point
     if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) {
+    if (BTS->GetIntersectionInsideTriangle(*Center, Walker->getPosition(), Intersection)) {
       DoLog(0) && (Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl);
       // we have the intersection, check whether in- or outside of boundary
       if ((Center->DistanceSquared(*Walker->node) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
+      if ((Center->DistanceSquared(Walker->getPosition()) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
         // inside, next!
         DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
 …
   DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
   for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
     DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
   // remove triangles's endpoints
 …
     DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
       DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
+      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
+  }
   delete (ListofPoints);
 …
         if (List != NULL) {
           for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             if ((*Runner)->node->at(map[0]) > maxCoordinate[map[0]]) {
               DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);
               maxCoordinate[map[0]] = (*Runner)->node->at(map[0]);
+            if ((*Runner)->at(map[0]) > maxCoordinate[map[0]]) {
+              DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << (*Runner)->getPosition() << "." << endl);
+              maxCoordinate[map[0]] = (*Runner)->at(map[0]);
               MaxPoint[map[0]] = (*Runner);
+            }
 …
     // construct center of circle
     CircleCenter = 0.5 * ((*BaseLine->endpoints[0]->node->node) + (*BaseLine->endpoints[1]->node->node));
+    CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
     // construct normal vector of circle
     CirclePlaneNormal = (*BaseLine->endpoints[0]->node->node) - (*BaseLine->endpoints[1]->node->node);
+    CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
     double radius = CirclePlaneNormal.NormSquared();
 …
   // construct center of circle
   CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
                         (*CandidateLine.BaseLine->endpoints[1]->node->node));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
   // construct normal vector of circle
   CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
                       (*CandidateLine.BaseLine->endpoints[1]->node->node);
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
   // calculate squared radius of circle
 …
     // construct SearchDirection and an "outward pointer"
     SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
     helper = CircleCenter - (*ThirdPoint->node->node);
+    helper = CircleCenter - (ThirdPoint->node->getPosition());
     if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
 …
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
   TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition());
   DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
 …
   // create normal vector
   BTS->GetCenter(&Center);
+  BTS->GetCenter(Center);
   Center -= CandidateLine.OptCenter;
   BTS->SphereCenter = CandidateLine.OptCenter;
 …
   // create normal vector
   BTS->GetCenter(&Center);
+  BTS->GetCenter(Center);
   Center.SubtractVector(*OptCenter);
   BTS->SphereCenter = *OptCenter;
 …
   Vector DistanceToIntersection[2], BaseLine;
   double distance[2];
   BaseLine = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
+  BaseLine = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
   for (int i = 0; i < 2; i++) {
     DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
+    DistanceToIntersection[i] = (*ClosestPoint) - (Base->endpoints[i]->node->getPosition());
     distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
+  }
 …
   // calculate volume
   volume = CalculateVolumeofGeneralTetraeder(*Base->endpoints[1]->node->node, *OtherBase->endpoints[0]->node->node, *OtherBase->endpoints[1]->node->node, *Base->endpoints[0]->node->node);
+  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->getPosition(), OtherBase->endpoints[0]->node->getPosition(), OtherBase->endpoints[1]->node->getPosition(), Base->endpoints[0]->node->getPosition());
   // delete the temporary other base line and the closest points
 …
               OrthogonalizedOben = Oben;
               aCandidate = (*a->node) - (*Candidate->node);
+              aCandidate = (a->getPosition()) - (Candidate->getPosition());
               OrthogonalizedOben.ProjectOntoPlane(aCandidate);
               OrthogonalizedOben.Normalize();
 …
               OrthogonalizedOben.Scale(scaleFactor);
               Center = 0.5 * ((*Candidate->node) + (*a->node));
+              Center = 0.5 * ((Candidate->getPosition()) + (a->getPosition()));
               Center += OrthogonalizedOben;
               AngleCheck = Center - (*a->node);
+              AngleCheck = Center - (a->getPosition());
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
 …
   // construct center of circle
   CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
                         (*CandidateLine.BaseLine->endpoints[1]->node->node));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
   // construct normal vector of circle
   CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
                       (*CandidateLine.BaseLine->endpoints[1]->node->node);
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
   RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
 …
       // get cell for the starting point
       if (LC->SetIndexToVector(&CircleCenter)) {
+      if (LC->SetIndexToVector(CircleCenter)) {
         for (int i = 0; i < NDIM; i++) // store indices of this cell
           N[i] = LC->n[i];
 …
                   // find center on the plane
                   GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
+                  GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
                   DoLog(1) && (Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl);
                   try {
                     NewNormalVector = Plane(*(CandidateLine.BaseLine->endpoints[0]->node->node),
                                             *(CandidateLine.BaseLine->endpoints[1]->node->node),
                                             *(Candidate->node)).getNormal();
+                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
+                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
+                                            (Candidate->getPosition())).getNormal();
                     DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
                     radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(NewPlaneCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter);
                     DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
                     DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
                     DoLog(1) && (Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl);
                     if (radius < RADIUS * RADIUS) {
                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(NewPlaneCenter);
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
                       if (fabs(radius - otherradius) < HULLEPSILON) {
                         // construct both new centers
 …
  * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
  */
 DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector &x, const LinkedCell* LC) const
+{
   Info FunctionInfo(__func__);
 …
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end()) {
               points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(*x), FindPoint->second));
+              points->insert(DistanceToPointPair(FindPoint->second->node->DistanceSquared(x), FindPoint->second));
               DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
+            }
 …
  * \return closest BoundaryLineSet or NULL in degenerate case.
  */
 BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector &x, const LinkedCell* LC) const
+{
   Info FunctionInfo(__func__);
 …
   // for each point, check its lines, remember closest
   DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl);
+  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << x << " ... " << endl);
   BoundaryLineSet *ClosestLine = NULL;
   double MinDistance = -1.;
 …
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       // calculate closest point on line to desired point
       helper = 0.5 * ((*(LineRunner->second)->endpoints[0]->node->node) +
                       (*(LineRunner->second)->endpoints[1]->node->node));
       Center = (*x) - helper;
       BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
                  (*(LineRunner->second)->endpoints[1]->node->node);
+      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) +
+                      ((LineRunner->second)->endpoints[1]->node->getPosition()));
+      Center = (x) - helper;
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
       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 = (*x) - (*(LineRunner->second)->endpoints[0]->node->node) - Center;
+        helper = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition()) - Center;
         const double lengthA = helper.ScalarProduct(BaseLine);
         helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
+        helper = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition()) - Center;
         const double lengthB = helper.ScalarProduct(BaseLine);
         if (lengthB * lengthA < 0) { // if have different sign
 …
  * \return BoundaryTriangleSet of nearest triangle or NULL.
  */
 TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector &x, const LinkedCell* LC) const
+{
   Info FunctionInfo(__func__);
 …
   // for each point, check its lines, remember closest
   DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl);
+  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << x << " ... " << endl);
   LineSet ClosestLines;
   double MinDistance = 1e+16;
 …
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
                  (*(LineRunner->second)->endpoints[1]->node->node);
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
       const double lengthBase = BaseLine.NormSquared();
       BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
+      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition());
       const double lengthEndA = BaseLineIntersection.NormSquared();
       BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
+      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
       const double lengthEndB = BaseLineIntersection.NormSquared();
 …
       } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
         BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
+        BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
         Center = BaseLineIntersection;
         Center.ProjectOntoPlane(BaseLine);
 …
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector &x, const LinkedCell* LC) const
+{
   Info FunctionInfo(__func__);
 …
   double MinAlignment = 2. * M_PI;
   for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
     (*Runner)->GetCenter(&Center);
     helper = (*x) - Center;
+    (*Runner)->GetCenter(Center);
+    helper = (x) - Center;
     const double Alignment = helper.Angle((*Runner)->NormalVector);
     if (Alignment < MinAlignment) {
 …
+{
   Info FunctionInfo(__func__);
   TriangleIntersectionList Intersections(&Point, this, LC);
+  TriangleIntersectionList Intersections(Point, this, LC);
   return Intersections.IsInside();
 …
+  }
   triangle->GetCenter(&Center);
+  triangle->GetCenter(Center);
   DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
   DistanceToCenter = Center - Point;
 …
     Center = Point - triangle->NormalVector; // points towards MolCenter
     DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);
     if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
+    if (triangle->GetIntersectionInsideTriangle(Center, DistanceToCenter, Intersection)) {
       DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);
       return 0.;
 …
   } else {
     // calculate smallest distance
     distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
+    distance = triangle->GetClosestPointInsideTriangle(Point, Intersection);
     DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);
 …
+{
   Info FunctionInfo(__func__);
   TriangleIntersectionList Intersections(&Point, this, LC);
+  TriangleIntersectionList Intersections(Point, this, LC);
   return Intersections.GetSmallestDistance();
 …
+{
   Info FunctionInfo(__func__);
   TriangleIntersectionList Intersections(&Point, this, LC);
+  TriangleIntersectionList Intersections(Point, this, LC);
   return Intersections.GetClosestTriangle();
 …
  * @return list of the all points linked to the provided one
  */
 TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
+{
   Info FunctionInfo(__func__);
 …
   // construct one orthogonal vector
+  if (Reference != NULL) {
+    AngleZero = (*Reference) - (*Point->node);
+    AngleZero.ProjectOntoPlane(PlaneNormal);
+  }
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
+    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
+    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
+  AngleZero = (Reference) - (Point->getPosition());
+  AngleZero.ProjectOntoPlane(PlaneNormal);
+  if ((AngleZero.NormSquared() < MYEPSILON)) {
+    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
+    AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
 …
   // go through all connected points and calculate angle
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper = (*(*listRunner)->node) - (*Point->node);
+    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
     helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
 …
  * @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
+ * @param *Reference Reference vector for zero angle or (0,0,0) 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
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
+{
   Info FunctionInfo(__func__);
 …
+  }
   DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl);
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << Reference << "." << endl);
   // calculate central point
   TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
 …
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
     helper = Plane(*((*TesselA)->node),
                    *((*TesselB)->node),
                    *((*TesselC)->node)).getNormal();
+    helper = Plane(((*TesselA)->getPosition()),
+                   ((*TesselB)->getPosition()),
+                   ((*TesselC)->getPosition())).getNormal();
     DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);
     counter++;
 …
   // construct one orthogonal vector
   if (Reference != NULL) {
     AngleZero = (*Reference) - (*Point->node);
+  if (!Reference.IsZero()) {
+    AngleZero = (Reference) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
+  }
   if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
     DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
     AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
+  if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON )) {
+    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
+    AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
 …
   pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper = (*(*listRunner)->node) - (*Point->node);
+    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
     helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
 …
   // copy old location for the volume
   OldPoint = (*point->node->node);
+  OldPoint = (point->node->getPosition());
   // get list of connected points
 …
           StartNode--;
           //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
           Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
+          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
           StartNode = MiddleNode;
           StartNode++;
 …
             StartNode = connectedPath->begin();
           //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
           Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
           OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
+          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
           OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
 …
         AddTesselationTriangle();
         // calculate volume summand as a general tetraeder
         volume += CalculateVolumeofGeneralTetraeder(*TPS[0]->node->node, *TPS[1]->node->node, *TPS[2]->node->node, OldPoint);
+        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
         // advance number
         count++;
 …
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
   class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->getPosition(), BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
 …
   class BoundaryLineSet *BestLine = NULL;
   for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
     BaseLine = (*Runner->second->endpoints[0]->node->node) -
                (*Runner->second->endpoints[1]->node->node);
     CenterToPoint = 0.5 * ((*Runner->second->endpoints[0]->node->node) +
                            (*Runner->second->endpoints[1]->node->node));
     CenterToPoint -= (*point->node);
+    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) -
+               (Runner->second->endpoints[1]->node->getPosition());
+    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) +
+                           (Runner->second->endpoints[1]->node->getPosition()));
+    CenterToPoint -= (point->getPosition());
     angle = CenterToPoint.Angle(BaseLine);
     if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {

src/tesselation.hpp

-              r8d6d31
+              r8f4df1
  * tesselation.hpp
+ *
+ *  The tesselation class is meant to contain the envelope (concave, convex or neither) of a set of Vectors. As we actually mean this stuff for atoms, we have to encapsulate it all a bit.
+ *  The tesselation class is meant to contain the envelope (concave, convex or neither) of a set of Vectors.
+ *  As we actually mean this stuff for atoms, we have to encapsulate it all a bit.
+ *
  *  Created on: Aug 3, 2009
 …
 #include <stack>
+#include "BoundaryMaps.hpp"
+#include "PointCloud.hpp"
+#include "TesselPoint.hpp"
 #include "atom_particleinfo.hpp"
 #include "Helpers/helpers.hpp"
 #include "LinearAlgebra/Vector.hpp"
 /****************************************** forward declarations *****************************/
 …
 class BoundaryLineSet;
 class BoundaryTriangleSet;
+class CandidateForTesselation;
 class LinkedCell;
-class TesselPoint;
-class PointCloud;
 class Tesselation;
 class Plane;
 …
 // ======================================================= some template functions =========================================
-#define IndexToIndex map <int, int>
-#define PointMap map < int, class BoundaryPointSet * >
-#define PointSet set < class BoundaryPointSet * >
-#define PointList list < class BoundaryPointSet * >
-#define PointPair pair < int, class BoundaryPointSet * >
-#define PointTestPair pair < PointMap::iterator, bool >
-#define CandidateList list <class CandidateForTesselation *>
-#define CandidateMap map <class BoundaryLineSet *, class CandidateForTesselation *>
-#define LineMap multimap < int, class BoundaryLineSet * >
-#define LineSet set < class BoundaryLineSet * >
-#define LineList list < class BoundaryLineSet * >
-#define LinePair pair < int, class BoundaryLineSet * >
-#define LineTestPair pair < LineMap::iterator, bool >
-#define TriangleMap map < int, class BoundaryTriangleSet * >
-#define TriangleSet set < class BoundaryTriangleSet * >
-#define TriangleList list < class BoundaryTriangleSet * >
-#define TrianglePair pair < int, class BoundaryTriangleSet * >
-#define TriangleTestPair pair < TrianglePair::iterator, bool >
-#define PolygonMap map < int, class BoundaryPolygonSet * >
-#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> >
-#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
-#define TesselPointList list <TesselPoint *>
-#define TesselPointSet set <TesselPoint *>
-#define ListOfTesselPointList list<list <TesselPoint *> *>
-enum centers {Opt, OtherOpt};
 /********************************************** declarations *******************************/
-template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
+{
-  if (endpoint1->Nr < endpoint2->Nr) {
-    endpoints[0] = endpoint1;
-    endpoints[1] = endpoint2;
-  } else {
-    endpoints[0] = endpoint2;
-    endpoints[1] = endpoint1;
+  }
-};
-// ======================================================== class BoundaryPointSet =========================================
-class BoundaryPointSet {
-  public:
-    BoundaryPointSet();
-    BoundaryPointSet(TesselPoint * const Walker);
-    ~BoundaryPointSet();
-    void AddLine(BoundaryLineSet * const line);
-    LineMap lines;
-    int LinesCount;
-    TesselPoint *node;
-    double value;
-    int Nr;
-};
-ostream & operator << (ostream &ost, const BoundaryPointSet &a);
-// ======================================================== class BoundaryLineSet ==========================================
-class BoundaryLineSet {
-  public:
-    BoundaryLineSet();
-    BoundaryLineSet(BoundaryPointSet * const Point[2], const int number);
-    BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number);
-    ~BoundaryLineSet();
-    void AddTriangle(BoundaryTriangleSet * const triangle);
-    bool IsConnectedTo(const BoundaryLineSet * const line) const;
-    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
-    bool CheckConvexityCriterion() const;
-    double CalculateConvexity() const;
-    class BoundaryPointSet *GetOtherEndpoint(const BoundaryPointSet * const point) const;
-    class BoundaryTriangleSet *GetOtherTriangle(const BoundaryTriangleSet * const triangle) const;
-    class BoundaryPointSet *endpoints[2];
-    TriangleMap triangles;
-    int Nr;
-    bool skipped;
-};
-ostream & operator << (ostream &ost, const BoundaryLineSet &a);
-// ======================================================== class BoundaryTriangleSet =======================================
-class BoundaryTriangleSet {
-  public:
-    BoundaryTriangleSet();
-    BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number);
-    ~BoundaryTriangleSet();
-    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;
-    class BoundaryLineSet *GetThirdLine(const BoundaryPointSet * const point) const;
-    bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
-    bool IsPresentTupel(const BoundaryTriangleSet * const T) const;
-    Plane getPlane() const;
-    Vector getEndpoint(int) const;
-    std::string getEndpointName(int) const;
-    class BoundaryPointSet *endpoints[3];
-    class BoundaryLineSet *lines[3];
-    Vector NormalVector;
-    Vector SphereCenter;
-    int Nr;
-  private:
-};
-ostream & operator << (ostream &ost, const BoundaryTriangleSet &a);
-// ======================================================== class BoundaryTriangleSet =======================================
-/** Set of BoundaryPointSet.
- * This is just meant as a container for a group of endpoints, extending the node, line, triangle concept. However, this has
- * only marginally something to do with the tesselation. Hence, there is no incorporation into the bookkeeping of the Tesselation
- * class (i.e. no allocation, no deletion).
- * \note we assume that the set of endpoints reside (more or less) on a plane.
- */
-class BoundaryPolygonSet {
-  public:
-    BoundaryPolygonSet();
-    ~BoundaryPolygonSet();
-    Vector * GetNormalVector(const Vector &NormalVector) const;
-    void GetCenter(Vector *center) const;
-    bool ContainsBoundaryLine(const BoundaryLineSet * const line) const;
-    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
-    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
-    bool ContainsBoundaryTriangle(const BoundaryTriangleSet * const point) const;
-    bool ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const;
-    bool ContainsPresentTupel(const BoundaryPolygonSet * const P) const;
-    bool ContainsPresentTupel(const PointSet &endpoints) const;
-    TriangleSet * GetAllContainedTrianglesFromEndpoints() const;
-    bool FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line);
-    PointSet endpoints;
-    int Nr;
-};
-ostream & operator << (ostream &ost, const BoundaryPolygonSet &a);
-// =========================================================== class TESSELPOINT ===========================================
-/** Is a single point of the set of Vectors, also a super-class to be inherited and and its functions to be implemented.
- */
-class TesselPoint : virtual public ParticleInfo {
-public:
-  TesselPoint();
-  virtual ~TesselPoint();
-  Vector *node;   // pointer to position of the dot in space
-  virtual ostream & operator << (ostream &ost);
-};
-ostream & operator << (ostream &ost, const TesselPoint &a);
-// =========================================================== class POINTCLOUD ============================================
-/** Super-class for all point clouds structures, also molecules. They have to inherit this structure and implement the virtual function to access the Vectors.
- * This basically encapsulates a list structure.
- */
-class PointCloud {
-public:
-  PointCloud();
-  virtual ~PointCloud();
-  virtual const char * const GetName() const { return "unknown"; };
-  virtual Vector *GetCenter() const { return NULL; };
-  virtual TesselPoint *GetPoint() const { return NULL; };
-  virtual int GetMaxId() const { return 0; };
-  virtual void GoToNext() const {};
-  virtual void GoToFirst() const {};
-  virtual bool IsEmpty() const { return true; };
-  virtual bool IsEnd() const { return true; };
-};
-// ======================================================== class CandidateForTesselation =========================================
-class CandidateForTesselation {
-  public :
-  CandidateForTesselation(BoundaryLineSet* currentBaseLine);
-  CandidateForTesselation(TesselPoint* candidate, BoundaryLineSet* currentBaseLine, BoundaryPointSet *point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
-  ~CandidateForTesselation();
-  bool CheckValidity(const double RADIUS, const LinkedCell *LC) const;
-  TesselPointList pointlist;
-  const BoundaryLineSet * BaseLine;
-  const BoundaryPointSet * ThirdPoint;
-  const BoundaryTriangleSet *T;
-  Vector OldCenter;
-  Vector OptCenter;
-  Vector OtherOptCenter;
-  double ShortestAngle;
-  double OtherShortestAngle;
-};
-ostream & operator <<(ostream &ost, const  CandidateForTesselation &a);
 // =========================================================== class TESSELATION ===========================================
 …
     void AddTesselationPoint(TesselPoint* Candidate, const int n);
     void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
     void AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
+    void AddTesselationLine(const Vector * OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddExistingTesselationTriangleLine(class BoundaryLineSet *FindLine, int n);
 …
     ListOfTesselPointList * GetPathsOfConnectedPoints(const TesselPoint* const Point) const;
     ListOfTesselPointList * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
     TesselPointList * GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
     TesselPointList * GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
+    TesselPointList * GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const;
+    TesselPointList * GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const;
     class BoundaryPointSet * GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
     TriangleList * FindTriangles(const TesselPoint* const Points[3]) const;
     TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
     BoundaryTriangleSet * FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const;
+    TriangleList * FindClosestTrianglesToVector(const Vector &x, const LinkedCell* LC) const;
+    BoundaryTriangleSet * FindClosestTriangleToVector(const Vector &x, const LinkedCell* LC) const;
     bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
     double GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const;
 …
     BoundaryTriangleSet * GetClosestTriangleOnSurface(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;
+    DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector &x, const LinkedCell* LC) const;
+    BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector &x, const LinkedCell* LC) const;
     // print for debugging

src/tesselationhelpers.cpp

-              r8d6d31
+              r8f4df1
 #include <fstream>
+#include "BoundaryLineSet.hpp"
+#include "BoundaryPointSet.hpp"
+#include "BoundaryPolygonSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "CandidateForTesselation.hpp"
 #include "Helpers/Info.hpp"
 #include "linkedcell.hpp"
 …
 #include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/Line.hpp"
 #include "vector_ops.hpp"
+#include "LinearAlgebra/vector_ops.hpp"
 #include "Helpers/Verbose.hpp"
 #include "LinearAlgebra/Plane.hpp"
 …
  * \param *c third point
  */
 void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c)
+void GetCenterofCircumcircle(Vector &Center, const Vector &a, const Vector &b, const Vector &c)
+{
         Info FunctionInfo(__func__);
 …
   helper[2] = SideC.NormSquared()*(SideA.NormSquared()+SideB.NormSquared() - SideC.NormSquared());
   Center->Zero();
   *Center += helper[0] * a;
   *Center += helper[1] * b;
   *Center += helper[2] * c;
   Center->Scale(1./(helper[0]+helper[1]+helper[2]));
   Log() << Verbose(1) << "INFO: Center (2nd algo) is at " << *Center << "." << endl;
+  Center.Zero();
+  Center += helper[0] * a;
+  Center += helper[1] * b;
+  Center += helper[2] * c;
+  Center.Scale(1./(helper[0]+helper[1]+helper[2]));
+  Log() << Verbose(1) << "INFO: Center (2nd algo) is at " << Center << "." << endl;
 };
 …
  * @return point which is second closest to the provided one
  */
 TesselPoint* FindSecondClosestTesselPoint(const Vector* Point, const LinkedCell* const LC)
+TesselPoint* FindSecondClosestTesselPoint(const Vector& Point, const LinkedCell* const LC)
+{
         Info FunctionInfo(__func__);
 …
         if (List != NULL) {
           for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             helper = (*Point) - (*(*Runner)->node);
+            helper = (Point) - ((*Runner)->getPosition());
             double currentNorm = helper. Norm();
             if (currentNorm < distance) {
 …
  * @return point which is closest to the provided one, NULL if none found
  */
 TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
+TesselPoint* FindClosestTesselPoint(const Vector& Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
+{
         Info FunctionInfo(__func__);
 …
         if (List != NULL) {
           for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             helper = (*Point) - (*(*Runner)->node);
+            helper = (Point) - ((*Runner)->getPosition());
             double currentNorm = helper.NormSquared();
             if (currentNorm < distance) {
 …
         Info FunctionInfo(__func__);
   // construct the plane of the two baselines (i.e. take both their directional vectors)
   Vector Baseline = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
   Vector OtherBaseline = (*OtherBase->endpoints[1]->node->node) - (*OtherBase->endpoints[0]->node->node);
+  Vector Baseline = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
+  Vector OtherBaseline = (OtherBase->endpoints[1]->node->getPosition()) - (OtherBase->endpoints[0]->node->getPosition());
   Vector Normal = Baseline;
   Normal.VectorProduct(OtherBaseline);
 …
   // project one offset point of OtherBase onto this plane (and add plane offset vector)
   Vector NewOffset = (*OtherBase->endpoints[0]->node->node) - (*Base->endpoints[0]->node->node);
+  Vector NewOffset = (OtherBase->endpoints[0]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
   NewOffset.ProjectOntoPlane(Normal);
   NewOffset += (*Base->endpoints[0]->node->node);
+  NewOffset += (Base->endpoints[0]->node->getPosition());
   Vector NewDirection = NewOffset + OtherBaseline;
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
   Line line1 = makeLineThrough(*(Base->endpoints[0]->node->node),*(Base->endpoints[1]->node->node));
+  Line line1 = makeLineThrough((Base->endpoints[0]->node->getPosition()),(Base->endpoints[1]->node->getPosition()));
   Line line2 = makeLineThrough(NewOffset, NewDirection);
   *Intersection = line1.getIntersection(line2);
   Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
+  Normal = (*Intersection) - (Base->endpoints[0]->node->getPosition());
   DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);
 …
       *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
       for (i=0;i<NDIM;i++)
         *vrmlfile << Walker->node->at(i)-center->at(i) << " ";
+        *vrmlfile << Walker->at(i)-center->at(i) << " ";
       *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
       cloud->GoToNext();
 …
       for (i=0;i<3;i++) { // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
           *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->at(j)-center->at(j) << " ";
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->at(j)-center->at(j) << " ";
         *vrmlfile << "\t";
+      }
 …
     Vector *center = cloud->GetCenter();
     // make the circumsphere's center absolute again
     Vector helper = (1./3.) * ((*Tess->LastTriangle->endpoints[0]->node->node) +
                                (*Tess->LastTriangle->endpoints[1]->node->node) +
                                (*Tess->LastTriangle->endpoints[2]->node->node));
+    Vector helper = (1./3.) * ((Tess->LastTriangle->endpoints[0]->node->getPosition()) +
+                               (Tess->LastTriangle->endpoints[1]->node->getPosition()) +
+                               (Tess->LastTriangle->endpoints[2]->node->getPosition()));
     helper -= (*center);
     // and add to file plus translucency object
 …
       *rasterfile << "2" << endl << "  ";  // 2 is sphere type
       for (int j=0;j<NDIM;j++) { // and for each node all NDIM coordinates
         const double tmp = Walker->node->at(j)-center->at(j);
+        const double tmp = Walker->at(j)-center->at(j);
         *rasterfile << ((fabs(tmp) < MYEPSILON) ? 0 : tmp) << " ";
+      }
 …
       for (i=0;i<3;i++) {  // print each node
         for (int j=0;j<NDIM;j++) { // and for each node all NDIM coordinates
           const double tmp = TriangleRunner->second->endpoints[i]->node->node->at(j)-center->at(j);
+          const double tmp = TriangleRunner->second->endpoints[i]->node->at(j)-center->at(j);
           *rasterfile << ((fabs(tmp) < MYEPSILON) ? 0 : tmp) << " ";
+        }
 …
       LookupList[Walker->nr] = Counter++;
       for (int i=0;i<NDIM;i++) {
         const double tmp = Walker->node->at(i);
+        const double tmp = Walker->at(i);
         *tecplot << ((fabs(tmp) < MYEPSILON) ? 0 : tmp) << " ";
+      }
 …
     TriangleSet *triangles = TesselStruct->GetAllTriangles(PointRunner->second);
     for (TriangleSet::iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
       totalarea += CalculateAreaofGeneralTriangle(*(*TriangleRunner)->endpoints[0]->node->node, *(*TriangleRunner)->endpoints[1]->node->node, *(*TriangleRunner)->endpoints[2]->node->node);
+      totalarea += CalculateAreaofGeneralTriangle((*TriangleRunner)->endpoints[0]->node->getPosition() , (*TriangleRunner)->endpoints[1]->node->getPosition() , (*TriangleRunner)->endpoints[2]->node->getPosition());
+    }
     ConcavityPerLine *= totalarea;
 …
       line = (*TriangleRunner)->GetThirdLine(PointRunner->second);
       triangle = line->GetOtherTriangle(*TriangleRunner);
       area = CalculateAreaofGeneralTriangle(*triangle->endpoints[0]->node->node, *triangle->endpoints[1]->node->node, *triangle->endpoints[2]->node->node);
       area += CalculateAreaofGeneralTriangle(*(*TriangleRunner)->endpoints[0]->node->node, *(*TriangleRunner)->endpoints[1]->node->node, *(*TriangleRunner)->endpoints[2]->node->node);
+      area = CalculateAreaofGeneralTriangle(triangle->endpoints[0]->node->getPosition() , triangle->endpoints[1]->node->getPosition() , triangle->endpoints[2]->node->getPosition());
+      area += CalculateAreaofGeneralTriangle((*TriangleRunner)->endpoints[0]->node->getPosition() , (*TriangleRunner)->endpoints[1]->node->getPosition() , (*TriangleRunner)->endpoints[2]->node->getPosition());
       area *= -line->CalculateConvexity();
       if (area > 0)
 …
     distance = 0.;
     for (TriangleMap::const_iterator TriangleRunner = Convex->TrianglesOnBoundary.begin(); TriangleRunner != Convex->TrianglesOnBoundary.end(); TriangleRunner++) {
       const double CurrentDistance = Convex->GetDistanceSquaredToTriangle(*PointRunner->second->node->node, TriangleRunner->second);
+      const double CurrentDistance = Convex->GetDistanceSquaredToTriangle(PointRunner->second->node->getPosition() , TriangleRunner->second);
       if (CurrentDistance < distance)
         distance = CurrentDistance;

src/tesselationhelpers.hpp

-              r8d6d31
+              r8f4df1
 class BoundaryLineSet;
 class BoundaryTriangleSet;
+class BoundaryPolygonSet;
+class CandidateForTesselation;
 class LinkedCell;
 class TesselPoint;
 …
 /********************************************** definitions *********************************/
-#define HULLEPSILON 1e-9 //!< TODO: Get rid of HULLEPSILON, points to numerical instabilities
 /********************************************** declarations *******************************/
 void GetSphere(Vector * const Center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS);
 void GetCenterOfSphere(Vector* const Center, const Vector &a, const Vector &b, const Vector &c, Vector * const NewUmkreismittelpunkt, const Vector* const Direction, const Vector* const AlternativeDirection, const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius);
 void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c);
+void GetSphere(Vector & Center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS);
+void GetCenterOfSphere(Vector& Center, const Vector &a, const Vector &b, const Vector &c, Vector &NewUmkreismittelpunkt, const Vector &Direction, const Vector* const AlternativeDirection, const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius);
+void GetCenterofCircumcircle(Vector & Center, const Vector &a, const Vector &b, const Vector &c);
 double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection);
 double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d);
 …
 bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3]);
 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation *candidate2);
 TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC);
 TesselPoint* FindSecondClosestTesselPoint(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);
 …
 void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct);
 void CalculateConstrictionPerBoundaryPoint(const Tesselation * const TesselStruct, const Tesselation * const Convex);
 double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle);
+double DistanceToTrianglePlane(const Vector &x, const BoundaryTriangleSet * const triangle);
 bool CheckListOfBaselines(const Tesselation * const TesselStruct);

src/triangleintersectionlist.cpp

-              r8d6d31
+              r8f4df1
 #include "Helpers/Info.hpp"
+#include "BoundaryMaps.hpp"
+#include "BoundaryPointSet.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
 #include "tesselation.hpp"
 #include "LinearAlgebra/Vector.hpp"
 …
+ *
  */
 TriangleIntersectionList::TriangleIntersectionList(const Vector *x, const Tesselation *surface, const LinkedCell* LC) :
+TriangleIntersectionList::TriangleIntersectionList(const Vector &x, const Tesselation *surface, const LinkedCell* LC) :
   Point(x),
   Tess(surface),
 …
+  }
   // return reference, if none can be found
   return *Point;
+  return Point;
 };
 …
   if (runner != IntersectionList.end()) {
     // if we have found one, check Scalarproduct between the vector
     Vector TestVector = (*Point) - (*(*runner).second);
+    Vector TestVector = (Point) - (*(*runner).second);
     if (fabs(TestVector.NormSquared()) < MYEPSILON) //
       return true;
 …
     // get intersection with triangle plane
     Intersection = new Vector;
     (*TriangleRunner)->GetClosestPointInsideTriangle(Point, Intersection);
     //Log() << Verbose(1) << "Intersection between " << *Point << " and " << **TriangleRunner << " is at " << *Intersection << "." << endl;
+    (*TriangleRunner)->GetClosestPointInsideTriangle(Point, *Intersection);
+    //Log() << Verbose(1) << "Intersection between " << Point << " and " << **TriangleRunner << " is at " << *Intersection << "." << endl;
     IntersectionList.insert( pair<BoundaryTriangleSet *, Vector * > (*TriangleRunner, Intersection) );
+  }
 …
   if (DistanceList.empty())
     for (TriangleVectorMap::const_iterator runner = IntersectionList.begin(); runner != IntersectionList.end(); runner++)
       DistanceList.insert( pair<double, BoundaryTriangleSet *> (Point->distance(*(*runner).second), (*runner).first) );
+      DistanceList.insert( pair<double, BoundaryTriangleSet *> (Point.distance(*(*runner).second), (*runner).first) );
   //for (DistanceTriangleMap::const_iterator runner = DistanceList.begin(); runner != DistanceList.end(); runner++)

src/triangleintersectionlist.hpp

-              r8d6d31
+              r8f4df1
+{
 public:
   TriangleIntersectionList(const Vector *x, const Tesselation *surface, const LinkedCell* LC);
+  TriangleIntersectionList(const Vector &x, const Tesselation *surface, const LinkedCell* LC);
   ~TriangleIntersectionList();
 …
   void FillDistanceList() const;
   const Vector *Point;
+  const Vector &Point;
   const Tesselation *Tess;
   const LinkedCell *Vicinity;

src/unittests/AnalysisCorrelationToPointUnitTest.cpp

-              r8d6d31
+              r8f4df1
   TestMolecule = World::getInstance().createMolecule();
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 0., 1. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 1., 1. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 1., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule->AddAtom(Walker);

src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

-              r8d6d31
+              r8f4df1
   ASSERT_DO(Assert::Throw);
+  setVerbosity(5);
   atom *Walker = NULL;
 …
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 0., 1. );
   TestSurfaceMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 1., 1. );
   TestSurfaceMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 1., 0. );
   TestSurfaceMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestSurfaceMolecule->AddAtom(Walker);
 …
   TestSurfaceMolecule = World::getInstance().createMolecule();
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
   *Walker->node = Vector(4., 0., 4. );
   TestSurfaceMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
   *Walker->node = Vector(0., 4., 4. );
   TestSurfaceMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
   *Walker->node = Vector(4., 4., 0. );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(4., 0., 4. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0., 4., 4. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(4., 4., 0. ));
   TestSurfaceMolecule->AddAtom(Walker);
   // add inner atoms
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
   *Walker->node = Vector(0.5, 0.5, 0.5 );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0.5, 0.5, 0.5 ));
   TestSurfaceMolecule->AddAtom(Walker);

src/unittests/AnalysisPairCorrelationUnitTest.cpp

-              r8d6d31
+              r8f4df1
   TestMolecule = World::getInstance().createMolecule();
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 0., 1. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 1., 1. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 1., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule->AddAtom(Walker);

src/unittests/CountBondsUnitTest.cpp

-              r8d6d31
+              r8f4df1
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(-0.2418, 0.9350, 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(-0.2418, 0.9350, 0. ));
   TestMolecule1->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(0.9658, 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0.9658, 0., 0. ));
   TestMolecule1->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = oxygen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(oxygen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule1->AddAtom(Walker);
 …
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(-0.2418, 0.9350, 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(-0.2418, 0.9350, 0. ));
   TestMolecule2->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(0.9658, 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0.9658, 0., 0. ));
   TestMolecule2->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = oxygen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(oxygen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule2->AddAtom(Walker);

src/unittests/LinkedCellUnitTest.cpp

-              r8d6d31
+              r8f4df1
         Walker = World::getInstance().createAtom();
         CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
         Walker->type = hydrogen;
         *Walker->node = Vector(x, y, z );
+        Walker->setType(hydrogen);
+        Walker->setPosition(Vector(x, y, z ));
         TestMolecule->AddAtom(Walker);
+      }
 …
   atom *newAtom = World::getInstance().createAtom();
   newAtom->setName("test");
   newAtom->x= Vector(1,1,1);
+  newAtom->setPosition(Vector(1,1,1));
   CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToNode(newAtom) );
   World::getInstance().destroyAtom(newAtom);
 …
   newAtom = World::getInstance().createAtom();
   newAtom->setName("test");
   newAtom->x = Vector(0,-1,0);
+  newAtom->setPosition(Vector(0,-1,0));
   CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToNode(newAtom) );
   World::getInstance().destroyAtom(newAtom);
 …
       for (double z=0.5;z<3;z+=1.) {
         tester = Vector(x,y,z);
         CPPUNIT_ASSERT_EQUAL( true, LC->SetIndexToVector(&tester) );
+        CPPUNIT_ASSERT_EQUAL( true, LC->SetIndexToVector(tester) );
+      }
   // check corners of each cell
 …
         tester= Vector(x,y,z);
         cout << "Tester is at " << tester << "." << endl;
         CPPUNIT_ASSERT_EQUAL( true, LC->SetIndexToVector(&tester) );
+        CPPUNIT_ASSERT_EQUAL( true, LC->SetIndexToVector(tester) );
+      }
   // check out of bounds
 …
         tester = Vector(x,y,z);
         cout << "The following test is supposed to fail and produce an ERROR." << endl;
         CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToVector(&tester) );
+        CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToVector(tester) );
+      }
   // check nonsense vectors
   tester= Vector(-423598,3245978,29349);
   cout << "The following test is supposed to fail and produce an ERROR." << endl;
   CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToVector(&tester) );
+  CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToVector(tester) );
 };
 …
   tester= Vector(0.5,0.5,0.5);
   LC->SetIndexToVector(&tester);
+  LC->SetIndexToVector(tester);
   LC->GetNeighbourBounds(lower, upper);
   for (int i=0;i<NDIM;i++)
 …
   // get all atoms
   tester= Vector(1.5,1.5,1.5);
   CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(&tester) );
+  CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(tester) );
   ListOfPoints = LC->GetallNeighbours();
   size = ListOfPoints->size();
 …
   // get all atoms in one corner
   tester= Vector(0.5, 0.5, 0.5);
   CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(&tester) );
+  CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(tester) );
   ListOfPoints = LC->GetallNeighbours();
   size=ListOfPoints->size();
   CPPUNIT_ASSERT_EQUAL( (size_t)8, size );
   for(molecule::iterator iter = TestMolecule->begin(); iter != TestMolecule->end(); ++iter){
     if (((*iter)->x[0] <2) && ((*iter)->x[1] <2) && ((*iter)->x[2] <2)) {
+    if (((*iter)->at(0) <2) && ((*iter)->at(1) <2) && ((*iter)->at(2) <2)) {
       ListOfPoints->remove(*iter);
       size--;
 …
   // get all atoms from one corner
   tester = Vector(0.5, 0.5, 0.5);
   CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(&tester) );
+  CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(tester) );
   ListOfPoints = LC->GetallNeighbours(3);
   size=ListOfPoints->size();
 …
   // get all points around central arom with radius 1.
   tester= Vector(1.5,1.5,1.5);
   CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(&tester) );
+  CPPUNIT_ASSERT_EQUAL ( true, LC->SetIndexToVector(tester) );
   ListOfPoints = LC->GetPointsInsideSphere(1., &tester);
   size = ListOfPoints->size();
   CPPUNIT_ASSERT_EQUAL( (size_t)7, size );
   for(molecule::iterator iter = TestMolecule->begin(); iter!=TestMolecule->end();++iter){
     if (((*iter)->x.DistanceSquared(tester) - 1.) < MYEPSILON ) {
+    if (((*iter)->DistanceSquared(tester) - 1.) < MYEPSILON ) {
       ListOfPoints->remove(*iter);
       size--;

src/unittests/ParserUnitTest.cpp

-              r8d6d31
+              r8f4df1
   input << Tremolo_Atomdata2;
   testParser->load(&input);
   CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->x[0] == 3.0);
+  CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->at(0) == 3.0);
   input.clear();
+}
 …
   input << Tremolo_velocity;
   testParser->load(&input);
   CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->v[0] == 3.0);
+  CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->AtomicVelocity[0] == 3.0);
   input.clear();
+}
 …
   // with the maximum number of fields and minimal information, default values are printed
   atom* newAtom = World::getInstance().createAtom();
   newAtom->type = World::getInstance().getPeriode()->FindElement(1);
+  newAtom->setType(1);
   testParser->setFieldsForSave("x=3 u=3 F stress Id neighbors=5 imprData GroupMeasureTypeNo Type extType name resName chainID resSeq occupancy tempFactor segID Charge charge GrpTypeNo torsion");
   testParser->save(&output);
 …
   atom *Walker = NULL;
   Walker = World::getInstance().createAtom();
   Walker->type = World::getInstance().getPeriode()->FindElement(8);
   Walker->x = Vector(0,0,0);
+  Walker->setType(8);
+  Walker->setPosition(Vector(0,0,0));
   Walker = World::getInstance().createAtom();
   Walker->type = World::getInstance().getPeriode()->FindElement(1);
   Walker->x = Vector(0.758602,0,0.504284);
+  Walker->setType(1);
+  Walker->setPosition(Vector(0.758602,0,0.504284));
   Walker = World::getInstance().createAtom();
   Walker->type = World::getInstance().getPeriode()->FindElement(1);
   Walker->x = Vector(0.758602,0,-0.504284);
+  Walker->setType(1);
+  Walker->setPosition(Vector(0.758602,0,-0.504284));
   CPPUNIT_ASSERT_EQUAL(3, World::getInstance().numAtoms());

src/unittests/analysisbondsunittest.cpp

-              r8d6d31
+              r8f4df1
 #include <cstring>
-#include "World.hpp"
 #include "analysis_bonds.hpp"
 #include "analysisbondsunittest.hpp"
 …
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "LinearAlgebra/Vector.hpp"
 #include "World.hpp"
 …
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(1.5, 0., 1.5 );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1.5, 0., 1.5 ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 1.5, 1.5 );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1.5, 1.5 ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(1.5, 1.5, 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1.5, 1.5, 0. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = carbon;
   *Walker->node = Vector(0.5, 0.5, 0.5 );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0.5, 0.5, 0.5 ));
   TestMolecule->AddAtom(Walker);

src/unittests/bondgraphunittest.cpp

-              r8d6d31
+              r8f4df1
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = carbon;
   *Walker->node = Vector(1., 0., 1. );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(1., 0., 1. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = carbon;
   *Walker->node = Vector(0., 1., 1. );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0., 1., 1. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = carbon;
   *Walker->node = Vector(1., 1., 0. );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(1., 1., 0. ));
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = carbon;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule->AddAtom(Walker);

src/unittests/listofbondsunittest.cpp

-              r8d6d31
+              r8f4df1
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 0., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 1., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(1., 1., 0. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
   Walker->type = hydrogen;
   *Walker->node = Vector(0., 0., 0. );
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
   TestMolecule->AddAtom(Walker);

src/unittests/tesselation_boundarytriangleunittest.cpp

-              r8d6d31
+              r8f4df1
 #include <cstring>
+#include <iostream>
 #include "defs.hpp"
+#include "tesselation.hpp"
+#include "TesselPoint.hpp"
+#include "BoundaryPointSet.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "CandidateForTesselation.hpp"
 #include "tesselation_boundarytriangleunittest.hpp"
 …
   // create nodes
   tesselpoints[0] = new TesselPoint;
   tesselpoints[0]->node = new Vector(0., 0., 0.);
+  tesselpoints[0]->setPosition(Vector(0., 0., 0.));
   tesselpoints[0]->setName("1");
   tesselpoints[0]->nr = 1;
   points[0] = new BoundaryPointSet(tesselpoints[0]);
   tesselpoints[1] = new TesselPoint;
   tesselpoints[1]->node = new Vector(0., 1., 0.);
+  tesselpoints[1]->setPosition(Vector(0., 1., 0.));
   tesselpoints[1]->setName("2");
   tesselpoints[1]->nr = 2;
   points[1] = new BoundaryPointSet(tesselpoints[1]);
   tesselpoints[2] = new TesselPoint;
   tesselpoints[2]->node = new Vector(1., 0., 0.);
+  tesselpoints[2]->setPosition(Vector(1., 0., 0.));
   tesselpoints[2]->setName("3");
   tesselpoints[2]->nr = 3;
 …
   for (int i=0;i<3;++i) {
     // TesselPoint does not delete its vector as it only got a reference
-    delete tesselpoints[i]->node;
     delete tesselpoints[i];
+  }
 …
   // simple test on y line
   Point = Vector(-1.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(-4.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( 16., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 16., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on x line
   Point = Vector(0.5,-1.,0.);
   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(0.5,-6.,0.);
   CPPUNIT_ASSERT_EQUAL( 36., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 36., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on slanted line
   Point = Vector(1.,1.,0.);
   CPPUNIT_ASSERT_EQUAL( 0.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 0.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(5.,5.,0.);
   CPPUNIT_ASSERT_EQUAL( 40.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 40.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on first node
   Point = Vector(-1.,-1.,0.);
   CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on second node
   Point = Vector(0.,2.,0.);
   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,1.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on third node
   Point = Vector(2.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(1.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // straight down/up
   Point = Vector(1./3.,1./3.,+5.);
   CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(1./3.,1./3.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(1./3.,1./3.,-5.);
   CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(1./3.,1./3.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on y line
   Point = Vector(-1.,0.5,+2.);
   CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(-1.,0.5,-3.);
   CPPUNIT_ASSERT_EQUAL( 10., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 10., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on x line
   Point = Vector(0.5,-1.,+1.);
   CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(0.5,-1.,-2.);
   CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on slanted line
   Point = Vector(1.,1.,+3.);
   CPPUNIT_ASSERT_EQUAL( 9.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 9.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
   Point = Vector(1.,1.,-4.);
   CPPUNIT_ASSERT_EQUAL( 16.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 16.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.5,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on first node
   Point = Vector(-1.,-1.,5.);
   CPPUNIT_ASSERT_EQUAL( 27., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 27., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on second node
   Point = Vector(0.,2.,5.);
   CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(0.,1.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 …
   // simple test on third node
   Point = Vector(2.,0.,5.);
   CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
   Point = Vector(1.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

src/unittests/tesselation_boundarytriangleunittest.hpp

r8d6d31	r8f4df1
15	15	#include "linkedcell.hpp"
16	16	#include "tesselation.hpp"
	17
	18	class TesselPoint;
17	19
18	20	/******************************************** Test classes ************************************/

src/unittests/tesselation_insideoutsideunittest.cpp

-              r8d6d31
+              r8f4df1
 #include <cstring>
+#include <iostream>
 #include "defs.hpp"
+#include "tesselation.hpp"
+#include "TesselPoint.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "CandidateForTesselation.hpp"
 #include "tesselation_insideoutsideunittest.hpp"
 #include "Helpers/Verbose.hpp"
 …
 CPPUNIT_TEST_SUITE_REGISTRATION( TesselationInOutsideTest );
 void TesselationInOutsideTest::setUp()
+{
 …
   class TesselPoint *Walker;
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 0., 0.);
+  Walker->setPosition(Vector(0., 0., 0.));
   Walker->setName("1");
   Walker->nr = 1;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 1., 0.);
+  Walker->setPosition(Vector(0., 1., 0.));
   Walker->setName("2");
   Walker->nr = 2;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 0., 0.);
+  Walker->setPosition(Vector(1., 0., 0.));
   Walker->setName("3");
   Walker->nr = 3;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 1., 0.);
+  Walker->setPosition(Vector(1., 1., 0.));
   Walker->setName("4");
   Walker->nr = 4;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 0., 1.);
+  Walker->setPosition(Vector(0., 0., 1.));
   Walker->setName("5");
   Walker->nr = 5;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 1., 1.);
+  Walker->setPosition(Vector(0., 1., 1.));
   Walker->setName("6");
   Walker->nr = 6;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 0., 1.);
+  Walker->setPosition(Vector(1., 0., 1.));
   Walker->setName("7");
   Walker->nr = 7;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 1., 1.);
+  Walker->setPosition(Vector(1., 1., 1.));
   Walker->setName("8");
   Walker->nr = 8;
 …
   delete(LinkedList);
   delete(TesselStruct);
+  for (LinkedCell::LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
+    delete((*Runner)->node);
+  for (LinkedCell::LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
     delete(*Runner);
+  }
   Corners.clear();
   logger::purgeInstance();

src/unittests/tesselationunittest.cpp

-              r8d6d31
+              r8f4df1
 #include "defs.hpp"
+#include "tesselation.hpp"
+#include "TesselPoint.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "CandidateForTesselation.hpp"
 #include "tesselationunittest.hpp"
 …
   class TesselPoint *Walker;
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 0., -1.);
+  Walker->setPosition(Vector(1., 0., -1.));
   Walker->setName("1");
   Walker->nr = 1;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 1., -1.);
+  Walker->setPosition(Vector(-1., 1., -1.));
   Walker->setName("2");
   Walker->nr = 2;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., -1., -1.);
+  Walker->setPosition(Vector(-1., -1., -1.));
   Walker->setName("3");
   Walker->nr = 3;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 0., 1.);
+  Walker->setPosition(Vector(-1., 0., 1.));
   Walker->setName("4");
   Walker->nr = 4;
 …
   delete(LinkedList);
   delete(TesselStruct);
+  for (LinkedCell::LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
+    delete((*Runner)->node);
+  for (LinkedCell::LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
     delete(*Runner);
+  }
   Corners.clear();
   logger::purgeInstance();

src/unittests/vectorunittest.cpp

-              r8d6d31
+              r8f4df1
 #include "Helpers/Log.hpp"
 #include "LinearAlgebra/Vector.hpp"
 #include "vector_ops.hpp"
+#include "LinearAlgebra/vector_ops.hpp"
 #include "vectorunittest.hpp"
 #include "LinearAlgebra/Plane.hpp"
 …
   errorLogger::purgeInstance();
 };
+/** UnitTest for Constructors and Vector::IsZero() and Vector::IsOne().
+ */
+void VectorTest::AssignmentTest()
+{
+  // test with zero
+  zero.at(0) = 0;
+  zero.at(1) = 0;
+  zero.at(2) = 0;
+  double zero_array[3] = {0., 0., 0.};
+  CPPUNIT_ASSERT_EQUAL( zero, Vector(0,0,0));
+  CPPUNIT_ASSERT_EQUAL( zero, Vector(0.,0.,0.));
+  CPPUNIT_ASSERT_EQUAL( zero, Vector(zero_array[0], zero_array[1], zero_array[2]));
+  CPPUNIT_ASSERT_EQUAL( zero, Vector(zero_array));
+  // test with unit
+  unit.at(0) = 1;
+  unit.at(1) = 0;
+  unit.at(2) = 0;
+  double unit_array[3] = {1., 0., 0.};
+  CPPUNIT_ASSERT_EQUAL( unit, Vector(1,0,0));
+  CPPUNIT_ASSERT_EQUAL( unit, Vector(1.,0.,0.));
+  CPPUNIT_ASSERT_EQUAL( unit, Vector(unit_array[0], unit_array[1], unit_array[2]));
+  CPPUNIT_ASSERT_EQUAL( unit, Vector(unit_array));
+  // test with two
+  two.at(0) = 2;
+  two.at(1) = 1;
+  two.at(2) = 0;
+  double two_array[3] = {2., 1., 0.};
+  CPPUNIT_ASSERT_EQUAL( two, Vector(2,1,0));
+  CPPUNIT_ASSERT_EQUAL( two, Vector(2.,1.,0.));
+  CPPUNIT_ASSERT_EQUAL( two, Vector(two_array[0], two_array[1], two_array[2]));
+  CPPUNIT_ASSERT_EQUAL( two, Vector(two_array));
+  // test with three
+  three.at(0) = 1;
+  three.at(1) = 2;
+  three.at(2) = 3;
+  double three_array[3] = {1., 2., 3.};
+  CPPUNIT_ASSERT_EQUAL( three, Vector(1,2,3));
+  CPPUNIT_ASSERT_EQUAL( three, Vector(1.,2.,3.));
+  CPPUNIT_ASSERT_EQUAL( three, Vector(three_array[0], three_array[1], three_array[2]));
+  CPPUNIT_ASSERT_EQUAL( three, Vector(three_array));
+}
 /** UnitTest for Constructors and Vector::IsZero() and Vector::IsOne().

src/unittests/vectorunittest.hpp

-              r8d6d31
+              r8f4df1
+{
     CPPUNIT_TEST_SUITE( VectorTest) ;
+    CPPUNIT_TEST ( AssignmentTest );
     CPPUNIT_TEST ( UnityTest );
     CPPUNIT_TEST ( SimpleAlgebraTest );
 …
     void tearDown();
+    void AssignmentTest();
     void UnityTest();
     void OperatorAlgebraTest();

Context Navigation

Changes in / [8d6d31:8f4df1]

Legend:

Download in other formats: