Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 23b547e0a280f674d85f218c4cce67c6c3f2ee85)
+++ src/boundary.cpp	(revision 0a4f7fd0df08c7bf001c21de0e7b7462fb56da1d)
@@ -18,4 +18,5 @@
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
+#include "Plane.hpp"
 
 #include<gsl/gsl_poly.h>
@@ -80,5 +81,5 @@
               DistanceVector.SubtractVector(&Neighbour->second.second->x);
               do { // seek for neighbour pair where it flips
-                  OldComponent = DistanceVector.x[Othercomponent];
+                  OldComponent = DistanceVector[Othercomponent];
                   Neighbour++;
                   if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
@@ -88,6 +89,6 @@
                   //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
-                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
-                  DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
+                  OldComponent) - DistanceVector[Othercomponent] / fabs(
+                  DistanceVector[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
               if (runner != Neighbour) {
                   OtherNeighbour = Neighbour;
@@ -102,8 +103,8 @@
                   //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
                   // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
-                  w1 = fabs(OtherVector.x[Othercomponent]);
-                  w2 = fabs(DistanceVector.x[Othercomponent]);
-                  tmp = fabs((w1 * DistanceVector.x[component] + w2
-                      * OtherVector.x[component]) / (w1 + w2));
+                  w1 = fabs(OtherVector[Othercomponent]);
+                  w2 = fabs(DistanceVector[Othercomponent]);
+                  tmp = fabs((w1 * DistanceVector[component] + w2
+                      * OtherVector[component]) / (w1 + w2));
                   // mark if it has greater diameter
                   //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
@@ -159,7 +160,7 @@
     AngleReferenceVector.Zero();
     AngleReferenceNormalVector.Zero();
-    AxisVector.x[axis] = 1.;
-    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
-    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+    AxisVector[axis] = 1.;
+    AngleReferenceVector[(axis + 1) % NDIM] = 1.;
+    AngleReferenceNormalVector[(axis + 2) % NDIM] = 1.;
 
     Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
@@ -636,5 +637,7 @@
       const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
       const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
-      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
+      x = Plane(*(runner->second->endpoints[0]->node->node),
+                *(runner->second->endpoints[1]->node->node),
+                *(runner->second->endpoints[2]->node->node)).getNormal();
       x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
       const double h = x.Norm(); // distance of CoG to triangle
@@ -751,14 +754,14 @@
     Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
     for (int i = 0; i < NDIM; i++)
-      BoxLengths.x[i] = GreatestDiameter[i];
+      BoxLengths[i] = GreatestDiameter[i];
     mol->CenterEdge(&BoxLengths);
   } else {
-    BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
-    BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
-    BoxLengths.x[2] = minimumvolume - cellvolume;
+    BoxLengths[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
+    BoxLengths[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
+    BoxLengths[2] = minimumvolume - cellvolume;
     double x0 = 0.;
     double x1 = 0.;
     double x2 = 0.;
-    if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1], BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
+    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
       Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
     else {
@@ -769,6 +772,6 @@
     cellvolume = 1.;
     for (int i = 0; i < NDIM; i++) {
-      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
-      cellvolume *= BoxLengths.x[i];
+      BoxLengths[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+      cellvolume *= BoxLengths[i];
     }
 
@@ -780,5 +783,5 @@
   // update Box of atoms by boundary
   mol->SetBoxDimension(&BoxLengths);
-  Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
 };
 
@@ -839,5 +842,5 @@
   FillerDistance.InverseMatrixMultiplication(M);
   for(int i=0;i<NDIM;i++)
-    N[i] = (int) ceil(1./FillerDistance.x[i]);
+    N[i] = static_cast<int>(ceil(1./FillerDistance[i]));
   Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;
 
@@ -880,5 +883,5 @@
           // create molecule random translation vector ...
           for (int i=0;i<NDIM;i++)
-            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+            FillerTranslations[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
           Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
 
@@ -892,5 +895,5 @@
             // create atomic random translation vector ...
             for (int i=0;i<NDIM;i++)
-              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+              AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
 
             // ... and rotation matrix
