Changes in / [192f6e:4e10f5]

Files:

• src/Actions/ActionRegistry.hpp (modified) (1 diff)
• src/Actions/AnalysisAction/PairCorrelationAction.cpp (modified) (2 diffs)
• src/Actions/AtomAction/AddAction.cpp (modified) (1 diff)
• src/Actions/AtomsCalculation_impl.hpp (modified) (1 diff)
• src/Actions/MapOfActions.cpp (modified) (1 diff)
• src/Actions/MoleculeAction/BondFileAction.cpp (modified) (1 diff)
• src/Actions/MoleculeAction/FillWithMoleculeAction.cpp (modified) (1 diff)
• src/Actions/MoleculeAction/SuspendInWaterAction.cpp (modified) (1 diff)
• src/Actions/MoleculeAction/TranslateAction.cpp (modified) (1 diff)
• src/Actions/WorldAction/AddEmptyBoundaryAction.cpp (modified) (3 diffs)
• src/Actions/WorldAction/CenterInBoxAction.cpp (modified) (1 diff)
• src/Actions/WorldAction/CenterOnEdgeAction.cpp (modified) (3 diffs)
• src/Actions/WorldAction/ChangeBoxAction.cpp (modified) (2 diffs)
• src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp (modified) (1 diff)
• src/Actions/WorldAction/RepeatBoxAction.cpp (modified) (6 diffs)
• src/Actions/WorldAction/ScaleBoxAction.cpp (modified) (3 diffs)
• src/Box.cpp (added)
• src/Box.hpp (added)
• src/ConfigFileBuffer.cpp (modified) (1 diff)
• src/ConfigFileBuffer.hpp (modified) (1 diff)
• src/Exceptions/CustomException.cpp (modified) (1 diff)
• src/Exceptions/CustomException.hpp (modified) (1 diff)
• src/Exceptions/NotInvertibleException.cpp (added)
• src/Exceptions/NotInvertibleException.hpp (added)
• src/Helpers/Assert.hpp (modified) (1 diff)
• src/Helpers/MemDebug.cpp (modified) (2 diffs)
• src/Helpers/MemDebug.hpp (modified) (3 diffs)
• src/Line.cpp (modified) (1 diff)
• src/Line.hpp (modified) (1 diff)
• src/Makefile.am (modified) (8 diffs)
• src/Matrix.cpp (added)
• src/Matrix.hpp (added)
• src/MatrixContent.hpp (added)
• src/Parser/MpqcParser.hpp (modified) (1 diff)
• src/Parser/PcpParser.cpp (modified) (6 diffs)
• src/Parser/PcpParser.hpp (modified) (1 diff)
• src/Patterns/Cacheable.hpp (modified) (1 diff)
• src/Plane.hpp (modified) (1 diff)
• src/Shapes/BaseShapes.cpp (added)
• src/Shapes/BaseShapes.hpp (added)
• src/Shapes/BaseShapes_impl.hpp (added)
• src/Shapes/Shape.cpp (added)
• src/Shapes/Shape.hpp (added)
• src/Shapes/ShapeOps.cpp (added)
• src/Shapes/ShapeOps.hpp (added)
• src/Shapes/ShapeOps_impl.hpp (added)
• src/Shapes/Shape_impl.hpp (added)
• src/UIElements/CommandLineUI/CommandLineDialog.cpp (modified) (4 diffs)
• src/UIElements/CommandLineUI/CommandLineDialog.hpp (modified) (3 diffs)
• src/UIElements/Dialog.cpp (modified) (4 diffs)
• src/UIElements/Dialog.hpp (modified) (5 diffs)
• src/UIElements/QT4/QTDialog.cpp (modified) (6 diffs)
• src/UIElements/QT4/QTDialog.hpp (modified) (2 diffs)
• src/UIElements/QT4/QTMainWindow.cpp (modified) (1 diff)
• src/UIElements/TextUI/TextDialog.cpp (modified) (4 diffs)
• src/UIElements/TextUI/TextDialog.hpp (modified) (3 diffs)
• src/UIElements/TextUI/TextWindow.hpp (modified) (1 diff)
• src/UIElements/Views/StreamStringView.hpp (modified) (1 diff)
• src/VectorContent.hpp (added)
• src/VectorSet.hpp (modified) (3 diffs)
• src/World.cpp (modified) (5 diffs)
• src/World.hpp (modified) (3 diffs)
• src/analysis_bonds.cpp (modified) (1 diff)
• src/analysis_correlation.cpp (modified) (16 diffs)
• src/analysis_correlation.hpp (modified) (1 diff)
• src/analyzer.cpp (modified) (1 diff)
• src/atom.cpp (modified) (2 diffs)
• src/atom.hpp (modified) (3 diffs)
• src/atom_graphnodeinfo.cpp (modified) (1 diff)
• src/atom_particleinfo.cpp (modified) (1 diff)
• src/atom_particleinfo.hpp (modified) (1 diff)
• src/atom_trajectoryparticle.hpp (modified) (1 diff)
• src/bond.cpp (modified) (1 diff)
• src/bondgraph.cpp (modified) (1 diff)
• src/bondgraph.hpp (modified) (1 diff)
• src/boundary.cpp (modified) (8 diffs)
• src/boundary.hpp (modified) (1 diff)
• src/builder.cpp (modified) (1 diff)
• src/config.cpp (modified) (7 diffs)
• src/datacreator.cpp (modified) (1 diff)
• src/datacreator.hpp (modified) (1 diff)
• src/element.hpp (modified) (1 diff)
• src/ellipsoid.cpp (modified) (4 diffs)
• src/errorlogger.cpp (modified) (1 diff)
• src/errorlogger.hpp (modified) (1 diff)
• src/graph.cpp (modified) (1 diff)
• src/graph.hpp (modified) (2 diffs)
• src/gslvector.cpp (modified) (2 diffs)
• src/gslvector.hpp (modified) (1 diff)
• src/helpers.cpp (modified) (2 diffs)
• src/helpers.hpp (modified) (2 diffs)
• src/joiner.cpp (modified) (1 diff)
• src/linkedcell.cpp (modified) (1 diff)
• src/logger.cpp (modified) (1 diff)
• src/logger.hpp (modified) (1 diff)
• src/memoryallocator.cpp (deleted)
• src/memoryallocator.hpp (deleted)
• src/memoryusageobserver.cpp (deleted)
• src/memoryusageobserver.hpp (deleted)
• src/molecule.cpp (modified) (8 diffs)
• src/molecule.hpp (modified) (4 diffs)
• src/molecule_dynamics.cpp (modified) (1 diff)
• src/molecule_fragmentation.cpp (modified) (4 diffs)
• src/molecule_geometry.cpp (modified) (13 diffs)
• src/molecule_graph.cpp (modified) (5 diffs)
• src/molecule_pointcloud.cpp (modified) (1 diff)
• src/moleculelist.cpp (modified) (6 diffs)
• src/parser.cpp (modified) (1 diff)
• src/periodentafel.hpp (modified) (1 diff)
• src/stackclass.hpp (modified) (1 diff)
• src/tesselation.cpp (modified) (1 diff)
• src/tesselationhelpers.cpp (modified) (3 diffs)
• src/tesselationhelpers.hpp (modified) (2 diffs)
• src/triangleintersectionlist.cpp (modified) (1 diff)
• src/unittests/ActOnAllUnitTest.cpp (modified) (1 diff)
• src/unittests/LineUnittest.cpp (modified) (2 diffs)
• src/unittests/LineUnittest.hpp (modified) (2 diffs)
• src/unittests/Makefile.am (modified) (8 diffs)
• src/unittests/MatrixUnittest.cpp (added)
• src/unittests/MatrixUnittest.hpp (added)
• src/unittests/PlaneUnittest.cpp (modified) (1 diff)
• src/unittests/ShapeUnittest.cpp (added)
• src/unittests/ShapeUnittest.hpp (added)
• src/unittests/linearsystemofequationsunittest.cpp (modified) (1 diff)
• src/unittests/tesselation_insideoutsideunittest.cpp (modified) (1 diff)
• src/unittests/vectorunittest.cpp (modified) (3 diffs)
• src/unittests/vectorunittest.hpp (modified) (2 diffs)
• src/vector.cpp (modified) (14 diffs)
• src/vector.hpp (modified) (7 diffs)
• src/vector_ops.cpp (modified) (1 diff)
• src/verbose.cpp (modified) (1 diff)
• src/verbose.hpp (modified) (1 diff)
• test_all.sh (added)
• tests/regression/Domain/4/post/test.conf (modified) (1 diff)

src/Actions/ActionRegistry.hpp

@@ -9 +9 @@
 #define ACTIONREGISTRY_HPP_
 
-#include <iostream>
+#include <iosfwd>
 #include <string>
 #include <map>

src/Actions/AnalysisAction/PairCorrelationAction.cpp

@@ -12 +12 @@
 #include "boundary.hpp"
 #include "linkedcell.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "element.hpp"
@@ -67 +68 @@
   dialog = UIFactory::getInstance().makeDialog();
   if (type == "P")
-    dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
+    dialog->queryVector("position", &Point, false, MapOfActions::getInstance().getDescription("position"));
   if (type == "S")
     dialog->queryMolecule("molecule-by-id", &Boundary, MapOfActions::getInstance().getDescription("molecule-by-id"));

src/Actions/AtomAction/AddAction.cpp

@@ -41 +41 @@
 
   dialog->queryElement(NAME, &elements, MapOfActions::getInstance().getDescription(NAME));
-  dialog->queryVector("position", &position, World::getInstance().getDomain(), true, MapOfActions::getInstance().getDescription("position"));
+  dialog->queryVector("position", &position, true, MapOfActions::getInstance().getDescription("position"));
   cout << "pre-dialog" << endl;
 

src/Actions/AtomsCalculation_impl.hpp

@@ -11 +11 @@
 #include "Actions/AtomsCalculation.hpp"
 #include "Actions/Calculation_impl.hpp"
-
-#include <iostream>
 
 using namespace std;

src/Actions/MapOfActions.cpp

@@ -17 +17 @@
 #include <boost/optional.hpp>
 #include <boost/program_options.hpp>
+
+#include <iostream>
 
 #include "CommandLineParser.hpp"

src/Actions/MoleculeAction/BondFileAction.cpp

@@ -24 +24 @@
 #include "config.hpp"
 #include "defs.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"

src/Actions/MoleculeAction/FillWithMoleculeAction.cpp

@@ -62 +62 @@
 
   dialog->queryString(NAME, &filename, MapOfActions::getInstance().getDescription(NAME));
-  dialog->queryVector("distances", &distances, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("distances"));
-  dialog->queryVector("lengths", &lengths, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("lengths"));
+  dialog->queryVector("distances", &distances, false, MapOfActions::getInstance().getDescription("distances"));
+  dialog->queryVector("lengths", &lengths, false, MapOfActions::getInstance().getDescription("lengths"));
   dialog->queryBoolean("DoRotate", &DoRotate, MapOfActions::getInstance().getDescription("DoRotate"));
   dialog->queryDouble("MaxDistance", &MaxDistance, MapOfActions::getInstance().getDescription("MaxDistance"));

src/Actions/MoleculeAction/SuspendInWaterAction.cpp

@@ -22 +22 @@
 #include "boundary.hpp"
 #include "config.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "config.hpp"

src/Actions/MoleculeAction/TranslateAction.cpp

@@ -55 +55 @@
   bool periodic = false;
 
-  dialog->queryVector(NAME, &x, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &x, false, MapOfActions::getInstance().getDescription(NAME));
   dialog->queryMolecule("molecule-by-id", &mol, MapOfActions::getInstance().getDescription("molecule-by-id"));
   dialog->queryBoolean("periodic", &periodic, MapOfActions::getInstance().getDescription("periodic"));

src/Actions/WorldAction/AddEmptyBoundaryAction.cpp

@@ -41 +41 @@
   int j=0;
 
-  dialog->queryVector(NAME, &boundary, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &boundary, false, MapOfActions::getInstance().getDescription(NAME));
 
   if(dialog->display()) {
@@ -59 +59 @@
     }
     // set new box size
-    double * const cell_size = World::getInstance().getDomain();
+    double * const cell_size = new double[6];
     for (j=0;j<6;j++)
       cell_size[j] = 0.;
@@ -67 +67 @@
       cell_size[j] = (Max[i]-Min[i]+2.*boundary[i]);
     }
+    World::getInstance().setDomain(cell_size);
+    delete[] cell_size;
     // translate all atoms, such that Min is aty (0,0,0)
     AtomRunner = AllAtoms.begin();

src/Actions/WorldAction/CenterInBoxAction.cpp

@@ -34 +34 @@
   Dialog *dialog = UIFactory::getInstance().makeDialog();
 
-  double * cell_size = World::getInstance().getDomain();
+  Box& cell_size = World::getInstance().getDomain();
   dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
 
   if(dialog->display()) {
-    World::getInstance().setDomain(cell_size);
     // center
     vector<molecule *> AllMolecules = World::getInstance().getAllMolecules();

src/Actions/WorldAction/CenterOnEdgeAction.cpp

@@ -13 +13 @@
 #include "vector.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
 
 #include <iostream>
@@ -37 +38 @@
   Vector Min;
   Vector Max;
-  int j=0;
 
   dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
@@ -57 +57 @@
     }
     // set new box size
-    double * const cell_size = World::getInstance().getDomain();
-    for (j=0;j<6;j++)
-      cell_size[j] = 0.;
-    j=-1;
+    Matrix domain;
     for (int i=0;i<NDIM;i++) {
-      j += i+1;
-      cell_size[j] = (Max[i]-Min[i]);
+      double tmp = Max[i]-Min[i];
+      tmp = fabs(tmp)>=1. ? tmp : 1.0;
+      domain.at(i,i) = tmp;
     }
-    World::getInstance().setDomain(cell_size);
+    World::getInstance().setDomain(domain);
     // translate all atoms, such that Min is aty (0,0,0)
     for (vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner)

src/Actions/WorldAction/ChangeBoxAction.cpp

@@ -12 +12 @@
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
+#include "Matrix.hpp"
 
 #include <iostream>
@@ -34 +36 @@
   Dialog *dialog = UIFactory::getInstance().makeDialog();
 
-  double * cell_size = World::getInstance().getDomain();
+  Box& cell_size = World::getInstance().getDomain();
   dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
 
   if(dialog->display()) {
-    DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size[0] << "," << cell_size[1] << "," << cell_size[2] << "," << cell_size[3] << "," << cell_size[4] << "," << cell_size[5] << "," << endl);
-    World::getInstance().setDomain(cell_size);
+    DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size.getM() << endl);
     delete dialog;
     return Action::success;

src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp

@@ -41 +41 @@
 
   dialog->queryDouble(NAME, &radius, MapOfActions::getInstance().getDescription(NAME));
-  dialog->queryVector("position", &point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
+  dialog->queryVector("position", &point, false, MapOfActions::getInstance().getDescription("position"));
 
   if(dialog->display()) {

src/Actions/WorldAction/RepeatBoxAction.cpp

@@ -13 +13 @@
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 
 #include <iostream>
@@ -46 +48 @@
   MoleculeListClass *molecules = World::getInstance().getMolecules();
 
-  dialog->queryVector(NAME, &Repeater, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &Repeater, false, MapOfActions::getInstance().getDescription(NAME));
   //dialog->queryMolecule("molecule-by-id", &mol,MapOfActions::getInstance().getDescription("molecule-by-id"));
   vector<molecule *> AllMolecules;
@@ -59 +61 @@
   if(dialog->display()) {
     (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
-    double * const cell_size = World::getInstance().getDomain();
-    double *M = ReturnFullMatrixforSymmetric(cell_size);
+    Matrix M = World::getInstance().getDomain().getM();
+    Matrix newM = M;
     Vector x,y;
     int n[NDIM];
+    Matrix repMat;
     for (int axis = 0; axis < NDIM; axis++) {
       Repeater[axis] = floor(Repeater[axis]);
@@ -69 +72 @@
         Repeater[axis] = 1;
       }
-      cell_size[ ((axis == 1) ? 2 : (axis == 2) ? 5 : 0) ] *= Repeater[axis];
+      repMat.at(axis,axis) = Repeater[axis];
     }
+    newM *= repMat;
+    World::getInstance().setDomain(newM);
 
-    World::getInstance().setDomain(cell_size);
     molecule *newmol = NULL;
     Vector ** vectors = NULL;
@@ -99 +103 @@
                 DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
               x = y;
-              x.MatrixMultiplication(M);
+              x *= M;
               newmol = World::getInstance().createMolecule();
               molecules->insert(newmol);
@@ -119 +123 @@
       }
     }
-    delete(M);
     delete dialog;
     return Action::success;

src/Actions/WorldAction/ScaleBoxAction.cpp

@@ -14 +14 @@
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
+#include "Matrix.hpp"
 
 #include <iostream>
@@ -37 +39 @@
   Vector Scaler;
   double x[NDIM];
-  int j=0;
 
-  dialog->queryVector(NAME, &Scaler, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &Scaler, false, MapOfActions::getInstance().getDescription(NAME));
 
   if(dialog->display()) {
@@ -49 +50 @@
       (*AtomRunner)->x.ScaleAll(x);
     }
-    j = -1;
-    double * const cell_size = World::getInstance().getDomain();
+
+    Matrix M = World::getInstance().getDomain().getM();
+    Matrix scale;
+
     for (int i=0;i<NDIM;i++) {
-      j += i+1;
-      cell_size[j]*=x[i];
+      scale.at(i,i) = x[i];
     }
-    World::getInstance().setDomain(cell_size);
+    M *= scale;
+    World::getInstance().setDomain(M);
 
     delete dialog;

src/ConfigFileBuffer.cpp

@@ -9 +9 @@
 #include "helpers.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "World.hpp"

src/ConfigFileBuffer.hpp

@@ -17 +17 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 
 /****************************************** forward declarations *****************************/

src/Exceptions/CustomException.cpp

@@ -9 +9 @@
 
 #include "CustomException.hpp"
+#include <iostream>
 
 using namespace std;

src/Exceptions/CustomException.hpp

@@ -10 +10 @@
 
 #include <exception>
-#include <iostream>
+#include <iosfwd>
+#include <string>
 
 /**

src/Helpers/Assert.hpp

@@ -11 +11 @@
 #include<sstream>
 #include<string>
-#include<iostream>
+#include<iosfwd>
 #include<vector>
 #include<map>

src/Helpers/MemDebug.cpp

@@ -6 +6 @@
  */
 
-#ifndef NDEBUG
-#ifndef NO_MEMDEBUG
+// NDEBUG implies NO_MEMDEBUG
+#ifdef NDEBUG
+# ifndef NO_MEMDEBUG
+#   define NO_MEMDEBUG
+# endif
+#endif
+
+// NO_MEMDEBUG and MEMDEBUG are mutually exclusive, but at least one must be set
+#ifdef NO_MEMDEBUG
+# ifdef MEMDEBUG
+#   undef MEMDEBUG
+# endif
+#else
+# ifndef MEMDEBUG
+#   define MEMDEBUG
+# endif
+#endif
+
+#ifdef MEMDEBUG
 
 #include <iostream>
@@ -489 +506 @@
 }
 #endif
-#endif

src/Helpers/MemDebug.hpp

@@ -21 +21 @@
  * your sourcefiles.
  */
-#ifndef NDEBUG
-#ifndef NO_MEMDEBUG
 
-#ifndef MEMDEBUG
-#define MEMDEBUG
+// Set all flags in a way that makes sense
+
+// NDEBUG implies NO_MEMDEBUG
+#ifdef NDEBUG
+# ifndef NO_MEMDEBUG
+#   define NO_MEMDEBUG
+# endif
 #endif
+
+// NO_MEMDEBUG and MEMDEBUG are mutually exclusive, but at least one must be set
+#ifdef NO_MEMDEBUG
+# ifdef MEMDEBUG
+#   undef MEMDEBUG
+# endif
+#else
+# ifndef MEMDEBUG
+#   define MEMDEBUG
+# endif
+#endif
+
+#ifdef MEMDEBUG
 
 #include <new>
@@ -83 +99 @@
 #endif
 
-#endif
-#endif
-
-
-#ifdef NDEBUG
-#undef MEMDEBUG
-#endif
-
-#ifndef MEMDEBUG
+#else
 // memory debugging was disabled
 
@@ -104 +112 @@
 
 #endif
+
+
 #endif /* MEMDEBUG_HPP_ */

src/Line.cpp

@@ -215 +215 @@
 }
 
+std::vector<Vector> Line::getSphereIntersections() const{
+  std::vector<Vector> res;
+
+  // line is kept in normalized form, so we can skip a lot of calculations
+  double discriminant = 1-origin->NormSquared();
+  // we might have 2, 1 or 0 solutions, depending on discriminant
+  if(discriminant>=0){
+    if(discriminant==0){
+      res.push_back(*origin);
+    }
+    else{
+      Vector helper = sqrt(discriminant)*(*direction);
+      res.push_back(*origin+helper);
+      res.push_back(*origin-helper);
+    }
+  }
+  return res;
+}
+
 Line makeLineThrough(const Vector &x1, const Vector &x2){
   if(x1==x2){

src/Line.hpp

@@ -38 +38 @@
   Plane getOrthogonalPlane(const Vector &origin) const;
 
+  std::vector<Vector> getSphereIntersections() const;
+
 private:
   std::auto_ptr<Vector> origin;

src/Makefile.am

@@ -38 +38 @@
   vector.cpp
                            
-LINALGHEADER = gslmatrix.hpp \
+LINALGHEADER = \
+  gslmatrix.hpp \
   gslvector.hpp \
   linearsystemofequations.hpp \
@@ -75 +76 @@
   Actions/MethodAction.hpp \
   Actions/Process.hpp
-  
+
+EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
+                                  Exceptions/LinearDependenceException.cpp \
+                                  Exceptions/MathException.cpp \
+                                  Exceptions/NotInvertibleException.cpp \
+                                  Exceptions/SkewException.cpp \
+                                  Exceptions/ZeroVectorException.cpp
+                                  
+EXCEPTIONHEADER = Exceptions/CustomException.hpp \
+                                  Exceptions/LinearDependenceException.hpp \
+                                  Exceptions/MathException.hpp \
+                                  Exceptions/NotInvertibleException.hpp \
+                                  Exceptions/SkewException.hpp \
+                                  Exceptions/ZeroVectorException.hpp
 
 PARSERSOURCE = \
@@ -101 +115 @@
   Patterns/Observer.hpp \
   Patterns/Singleton.hpp
+  
+SHAPESOURCE = \
+  Shapes/BaseShapes.cpp \
+  Shapes/Shape.cpp \
+  Shapes/ShapeOps.cpp
+SHAPEHEADER = \
+  Shapes/BaseShapes.hpp \
+  Shapes/Shape.hpp \
+  Shapes/ShapeOps.hpp
+  
 
 QTUIMOC_HEADER = UIElements/QT4/QTDialog.hpp \
@@ -134 +158 @@
   Descriptors/MoleculeNameDescriptor.hpp \
   Descriptors/MoleculePtrDescriptor.hpp
-                                   
-EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
-                                  Exceptions/LinearDependenceException.cpp \
-                                  Exceptions/MathException.cpp \
-                                  Exceptions/SkewException.cpp \
-                                  Exceptions/ZeroVectorException.cpp
-                                  
-EXCEPTIONHEADER = Exceptions/CustomException.hpp \
-                                  Exceptions/LinearDependenceException.hpp \
-                                  Exceptions/MathException.hpp \
-                                  Exceptions/SkewException.hpp \
-                                  Exceptions/ZeroVectorException.hpp
                                   
 QTUISOURCE = ${QTUIMOC_TARGETS} \
@@ -165 +177 @@
   ${ACTIONSSOURCE} \
   ${ATOMSOURCE} \
+  ${EXCEPTIONSOURCE} \
   ${PATTERNSOURCE} \
   ${PARSERSOURCE} \
+  ${SHAPESOURCE} \
   ${DESCRIPTORSOURCE} \
   ${HELPERSOURCE} \
-  ${EXCEPTIONSOURCE} \
   bond.cpp \
   bondgraph.cpp \
   boundary.cpp \
+  Box.cpp \
   CommandLineParser.cpp \
   config.cpp \
@@ -188 +202 @@
   log.cpp \
   logger.cpp \
+  Matrix.cpp \
   moleculelist.cpp \
   molecule.cpp \
@@ -212 +227 @@
   ${ACTIONSHEADER} \
   ${ATOMHEADER} \
+  ${EXCEPTIONHEADER} \
   ${PARSERHEADER} \
   ${PATTERNHEADER} \
+  ${SHAPEHEADER} \
   ${DESCRIPTORHEADER} \
-  ${EXCEPTIONHEADER} \
   bond.hpp \
   bondgraph.hpp \
   boundary.hpp \
+  Box.hpp \
   CommandLineParser.hpp \
   config.hpp \
@@ -236 +253 @@
   log.hpp \
   logger.hpp \
+  Matrix.hpp \
   molecule.hpp \
   molecule_template.hpp \

src/Parser/MpqcParser.hpp

@@ -11 +11 @@
 #include "FormatParser.hpp"
 
-#include <iostream>
+#include <iosfwd>
 
 /**

src/Parser/PcpParser.cpp

@@ -7 +7 @@
 
 #include <iostream>
+#include <iomanip>
 
 #include "atom.hpp"
@@ -20 +21 @@
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 
 /** Constructor of PcpParser.
@@ -210 +213 @@
   // Unit cell and magnetic field
   ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
-  double * const cell_size = World::getInstance().getDomain();
+  double *cell_size = new double[6];
   cell_size[0] = BoxLength[0];
   cell_size[1] = BoxLength[3];
@@ -217 +220 @@
   cell_size[4] = BoxLength[7];
   cell_size[5] = BoxLength[8];
+  World::getInstance().setDomain(cell_size);
+  delete[] cell_size;
   //if (1) fprintf(stderr,"\n");
 
@@ -327 +332 @@
 void PcpParser::save(std::ostream* file)
 {
-  const double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   class ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
   if (!file->fail()) {
@@ -412 +417 @@
     *file << endl;
     *file << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
-    *file << cell_size[0] << "\t" << endl;
-    *file << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
-    *file << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
+    *file << domain.at(0,0) << "\t" << endl;
+    *file << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
+    *file << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     // FIXME
     *file << endl;

src/Parser/PcpParser.hpp

@@ -9 +9 @@
 #define PCPPARSER_HPP_
 
-#include <iostream>
+#include <iosfwd>
 #include "Parser/FormatParser.hpp"
 

src/Patterns/Cacheable.hpp

@@ -12 +12 @@
 #include <boost/function.hpp>
 #include <boost/shared_ptr.hpp>
-#include <iostream>
 
 #include "Helpers/Assert.hpp"

src/Plane.hpp

@@ -11 +11 @@
 #include <memory>
 #include <vector>
-#include <iostream>
+#include <iosfwd>
 #include "Space.hpp"
 #include "Exceptions/LinearDependenceException.hpp"

src/UIElements/CommandLineUI/CommandLineDialog.cpp

@@ -27 +27 @@
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 
 #include "atom.hpp"
@@ -77 +78 @@
 }
 
-void CommandLineDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string _description) {
-  registerQuery(new VectorCommandLineQuery(title,target,cellSize,check, _description));
-}
-
-void CommandLineDialog::queryBox(const char* title, double ** const cellSize, string _description) {
+void CommandLineDialog::queryVector(const char* title, Vector *target, bool check, string _description) {
+  registerQuery(new VectorCommandLineQuery(title,target,check, _description));
+}
+
+void CommandLineDialog::queryBox(const char* title, Box* cellSize, string _description) {
   registerQuery(new BoxCommandLineQuery(title,cellSize,_description));
 }
@@ -222 +223 @@
 }
 
-CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, const double *const _cellSize, bool _check, string _description) :
-    Dialog::VectorQuery(title,_target,_cellSize,_check, _description)
+CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, bool _check, string _description) :
+    Dialog::VectorQuery(title,_target,_check, _description)
 {}
 
@@ -244 +245 @@
 
 
-CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, double ** const _cellSize, string _description) :
+CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, Box* _cellSize, string _description) :
     Dialog::BoxQuery(title,_cellSize, _description)
 {}

src/UIElements/CommandLineUI/CommandLineDialog.hpp

@@ -35 +35 @@
   virtual void queryAtom(const char*,atom**, std::string = "");
   virtual void queryMolecule(const char*,molecule**,std::string = "");
-  virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
-  virtual void queryBox(const char*,double ** const, std::string = "");
+  virtual void queryVector(const char*,Vector *,bool, std::string = "");
+  virtual void queryBox(const char*,Box *, std::string = "");
   virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
 
@@ -99 +99 @@
   class VectorCommandLineQuery : public Dialog::VectorQuery {
   public:
-    VectorCommandLineQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
+    VectorCommandLineQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
     virtual ~VectorCommandLineQuery();
     virtual bool handle();
@@ -106 +106 @@
   class BoxCommandLineQuery : public Dialog::BoxQuery {
   public:
-    BoxCommandLineQuery(std::string title,double ** const _cellSize, std::string _description = "");
+    BoxCommandLineQuery(std::string title,Box* _cellSize, std::string _description = "");
     virtual ~BoxCommandLineQuery();
     virtual bool handle();

src/UIElements/Dialog.cpp

@@ -10 +10 @@
 #include "Dialog.hpp"
 
+#include "verbose.hpp"
 #include "atom.hpp"
 #include "element.hpp"
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 
 using namespace std;
@@ -185 +188 @@
 // Vector Queries
 
-Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description) :
+Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,bool _check, std::string _description) :
   Query(title, _description),
-  cellSize(_cellSize),
   check(_check),
   target(_target)
@@ -205 +207 @@
 // Box Queries
 
-Dialog::BoxQuery::BoxQuery(std::string title, double ** const _cellSize, std::string _description) :
+Dialog::BoxQuery::BoxQuery(std::string title, Box* _cellSize, std::string _description) :
   Query(title, _description),
   target(_cellSize)
@@ -218 +220 @@
 
 void Dialog::BoxQuery::setResult() {
-  for (int i=0;i<6;i++) {
-    (*target)[i] = tmp[i];
-  }
+  (*target)= ReturnFullMatrixforSymmetric(tmp);
 }
 

src/UIElements/Dialog.hpp

@@ -14 +14 @@
 
 class atom;
+class Box;
 class element;
 class molecule;
@@ -43 +44 @@
   virtual void queryAtom(const char*,atom**,std::string = "")=0;
   virtual void queryMolecule(const char*,molecule**, std::string = "")=0;
-  virtual void queryVector(const char*,Vector *,const double *const,bool, std::string = "")=0;
-  virtual void queryBox(const char*,double ** const, std::string = "")=0;
+  virtual void queryVector(const char*,Vector *,bool, std::string = "")=0;
+  virtual void queryBox(const char*,Box*, std::string = "")=0;
   virtual void queryElement(const char*, std::vector<element *> *, std::string = "")=0;
 
@@ -176 +177 @@
   class VectorQuery : public Query {
   public:
-      VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
+      VectorQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
       virtual ~VectorQuery();
       virtual bool handle()=0;
@@ -182 +183 @@
     protected:
       Vector *tmp;
-      const double *const cellSize;
       bool check;
     private:
@@ -190 +190 @@
   class BoxQuery : public Query {
   public:
-      BoxQuery(std::string title,double ** const _cellSize, std::string _description = "");
+      BoxQuery(std::string title,Box *_cellSize, std::string _description = "");
       virtual ~BoxQuery();
       virtual bool handle()=0;
       virtual void setResult();
     protected:
-      double *tmp;
+      double* tmp;
     private:
-      double **target;
+      Box* target;
   };
 

src/UIElements/QT4/QTDialog.cpp

@@ -29 +29 @@
 #include "molecule.hpp"
 #include "Descriptors/MoleculeIdDescriptor.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 
 
@@ -93 +95 @@
 }
 
-void QTDialog::queryBox(char const*, double**, string){
+void QTDialog::queryBox(char const*, Box*, string){
   // TODO
   ASSERT(false, "Not implemented yet");
@@ -123 +125 @@
 }
 
-void QTDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check,string) {
-  registerQuery(new VectorQTQuery(title,target,cellSize,check,inputLayout,this));
+void QTDialog::queryVector(const char* title, Vector *target, bool check,string) {
+  registerQuery(new VectorQTQuery(title,target,check,inputLayout,this));
 }
 
@@ -291 +293 @@
 }
 
-QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
-    Dialog::VectorQuery(title,_target,_cellSize,_check),
-    parent(_parent)
-{
-  // About the j: I don't know why it was done this way, but it was used this way in Vector::AskPosition, so I reused it
-  int j = -1;
+QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
+    Dialog::VectorQuery(title,_target,_check),
+    parent(_parent)
+{
+  const Matrix& M = World::getInstance().getDomain().getM();
   const char *coords[3] = {"x:","y:","z:"};
   mainLayout= new QHBoxLayout();
@@ -304 +305 @@
   mainLayout->addLayout(subLayout);
   for(int i=0; i<3; i++) {
-    j+=i+1;
     coordLayout[i] = new QHBoxLayout();
     subLayout->addLayout(coordLayout[i]);
@@ -310 +310 @@
     coordLayout[i]->addWidget(coordLabel[i]);
     coordInput[i] = new QDoubleSpinBox();
-    coordInput[i]->setRange(0,cellSize[j]);
+    coordInput[i]->setRange(0,M.at(i,i));
     coordInput[i]->setDecimals(3);
     coordLayout[i]->addWidget(coordInput[i]);

src/UIElements/QT4/QTDialog.hpp

@@ -43 +43 @@
   virtual void queryAtom(const char*,atom**,std::string = "");
   virtual void queryMolecule(const char*,molecule**,std::string = "");
-  virtual void queryVector(const char*,Vector *,const double *const,bool,std::string = "");
-  virtual void queryBox(const char*,double ** const, std::string = "");
+  virtual void queryVector(const char*,Vector *,bool,std::string = "");
+  virtual void queryBox(const char*,Box*, std::string = "");
   virtual void queryElement(const char*,std::vector<element *> *_target,std::string = "");
 
@@ -124 +124 @@
     class VectorQTQuery : public Dialog::VectorQuery {
     public:
-      VectorQTQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check,QBoxLayout *,QTDialog *);
+      VectorQTQuery(std::string title,Vector *_target,bool _check,QBoxLayout *,QTDialog *);
       virtual ~VectorQTQuery();
       virtual bool handle();

src/UIElements/QT4/QTMainWindow.cpp

@@ -21 +21 @@
 #include "atom.hpp"
 #include "molecule.hpp"
+#include "verbose.hpp"
 #include "Actions/Action.hpp"
 #include "Actions/ActionRegistry.hpp"

src/UIElements/TextUI/TextDialog.cpp

@@ -25 +25 @@
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 
 using namespace std;
@@ -70 +72 @@
 }
 
-void TextDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string description) {
-  registerQuery(new VectorTextQuery(title,target,cellSize,check,description));
-}
-
-void TextDialog::queryBox(const char* title,double ** const cellSize, string description) {
+void TextDialog::queryVector(const char* title, Vector *target, bool check, string description) {
+  registerQuery(new VectorTextQuery(title,target,check,description));
+}
+
+void TextDialog::queryBox(const char* title,Box* cellSize, string description) {
   registerQuery(new BoxTextQuery(title,cellSize,description));
 }
@@ -270 +272 @@
 }
 
-TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check, std::string _description) :
-    Dialog::VectorQuery(title,_target,_cellSize,_check,_description)
+TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, bool _check, std::string _description) :
+    Dialog::VectorQuery(title,_target,_check,_description)
 {}
 
@@ -280 +282 @@
   Log() << Verbose(0) << getTitle();
 
+  const Matrix &M = World::getInstance().getDomain().getM();
   char coords[3] = {'x','y','z'};
-  int j = -1;
   for (int i=0;i<3;i++) {
-    j += i+1;
     do {
-      Log() << Verbose(0) << coords[i] << "[0.." << cellSize[j] << "]: ";
+      Log() << Verbose(0) << coords[i] << "[0.." << M.at(i,i) << "]: ";
       cin >> (*tmp)[i];
-    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= cellSize[j])) && (check));
+    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= M.at(i,i))) && (check));
   }
   return true;
 }
 
-TextDialog::BoxTextQuery::BoxTextQuery(std::string title, double ** const _cellSize, std::string _description) :
+TextDialog::BoxTextQuery::BoxTextQuery(std::string title, Box* _cellSize, std::string _description) :
     Dialog::BoxQuery(title,_cellSize,_description)
 {}

src/UIElements/TextUI/TextDialog.hpp

@@ -32 +32 @@
   virtual void queryAtom(const char*,atom**,std::string = "");
   virtual void queryMolecule(const char*,molecule**,std::string = "");
-  virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
-  virtual void queryBox(const char*,double ** const, std::string = "");
+  virtual void queryVector(const char*,Vector *,bool, std::string = "");
+  virtual void queryBox(const char*,Box*, std::string = "");
   virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
 
@@ -96 +96 @@
   class VectorTextQuery : public Dialog::VectorQuery {
   public:
-    VectorTextQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = NULL);
+    VectorTextQuery(std::string title,Vector *_target,bool _check, std::string _description = NULL);
     virtual ~VectorTextQuery();
     virtual bool handle();
@@ -103 +103 @@
   class BoxTextQuery : public Dialog::BoxQuery {
   public:
-    BoxTextQuery(std::string title,double ** const _cellSize, std::string _description = NULL);
+    BoxTextQuery(std::string title,Box* _cellSize, std::string _description = NULL);
     virtual ~BoxTextQuery();
     virtual bool handle();

src/UIElements/TextUI/TextWindow.hpp

@@ -11 +11 @@
 #include "MainWindow.hpp"
 
+#include <string>
 #include <set>
 

src/UIElements/Views/StreamStringView.hpp

@@ -10 +10 @@
 
 #include <boost/function.hpp>
-#include <iostream>
+#include <iosfwd>
 
 #include "Views/StringView.hpp"

src/VectorSet.hpp

@@ -12 +12 @@
 #include <functional>
 #include <algorithm>
+#include <limits>
 
 /**
@@ -19 +20 @@
  */
 
-class Vector;
+#include "vector.hpp"
+#include <list>
 
 // this tests, whether we actually have a Vector
@@ -48 +50 @@
   void translate(const Vector &translater){
     // this is needed to allow template lookup
-    transform(this->begin(),this->end(),this->begin(),bind1st(plus<Vector>(),translater));
+    transform(this->begin(),this->end(),this->begin(),std::bind1st(std::plus<Vector>(),translater));
+  }
+
+  double minDistSquared(const Vector &point){
+    if(!this->size())
+      return std::numeric_limits<double>::infinity();
+    std::list<double> helper;
+    helper.resize(this->size());
+    transform(this->begin(),this->end(),
+              helper.begin(),
+              std::bind2nd(std::mem_fun_ref(&Vector::DistanceSquared),point));
+    return *min_element(helper.begin(),helper.end());
   }
 };
 
+// allows simpler definition of VectorSets
+#define VECTORSET(container_type) VectorSet<container_type<Vector> >
+
 #endif /* VECTORSET_HPP_ */

src/World.cpp

@@ -22 +22 @@
 #include "Actions/ManipulateAtomsProcess.hpp"
 #include "Helpers/Assert.hpp"
+#include "Box.hpp"
+#include "Matrix.hpp"
 
 #include "Patterns/Singleton_impl.hpp"
@@ -74 +76 @@
 // system
 
-double * World::getDomain() {
-  return cell_size;
+Box& World::getDomain() {
+  return *cell_size;
+}
+
+void World::setDomain(const Matrix &mat){
+  *cell_size = mat;
 }
 
 void World::setDomain(double * matrix)
 {
-  OBSERVE;
-  cell_size[0] = matrix[0];
-  cell_size[1] = matrix[1];
-  cell_size[2] = matrix[2];
-  cell_size[3] = matrix[3];
-  cell_size[4] = matrix[4];
-  cell_size[5] = matrix[5];
+  Matrix M = ReturnFullMatrixforSymmetric(matrix);
+  cell_size->setM(M);
 }
 
@@ -139 +140 @@
   molecules.erase(id);
 }
-
-double *World::cell_size = NULL;
 
 atom *World::createAtom(){
@@ -303 +302 @@
     molecules_deprecated(new MoleculeListClass(this))
 {
-  cell_size = new double[6];
-  cell_size[0] = 20.;
-  cell_size[1] = 0.;
-  cell_size[2] = 20.;
-  cell_size[3] = 0.;
-  cell_size[4] = 0.;
-  cell_size[5] = 20.;
+  cell_size = new Box;
+  Matrix domain;
+  domain.at(0,0) = 20;
+  domain.at(1,1) = 20;
+  domain.at(2,2) = 20;
+  cell_size->setM(domain);
   defaultName = "none";
   molecules_deprecated->signOn(this);
@@ -317 +315 @@
 {
   molecules_deprecated->signOff(this);
-  delete[] cell_size;
+  delete cell_size;
   delete molecules_deprecated;
   delete periode;

src/World.hpp

@@ -34 +34 @@
 class AtomDescriptor_impl;
 template<typename T> class AtomsCalculation;
+class Box;
 class config;
 class ManipulateAtomsProcess;
+class Matrix;
 class molecule;
 class MoleculeDescriptor;
@@ -125 +127 @@
    * get the domain size as a symmetric matrix (6 components)
    */
-  double * getDomain();
+  Box& getDomain();
+
+  /**
+   * Set the domain size from a matrix object
+   *
+   * Matrix needs to be symmetric
+   */
+  void setDomain(const Matrix &mat);
 
   /**
@@ -257 +266 @@
   periodentafel *periode;
   config *configuration;
-  static double *cell_size;
+  Box *cell_size;
   std::string defaultName;
   class ThermoStatContainer *Thermostats;

src/analysis_bonds.cpp

@@ -13 +13 @@
 #include "element.hpp"
 #include "info.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"

src/analysis_correlation.cpp

@@ -9 +9 @@
 
 #include <iostream>
+#include <iomanip>
 
 #include "analysis_correlation.hpp"
@@ -19 +20 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 
 
@@ -34 +37 @@
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
-  double *domain = World::getInstance().getDomain();
+  Box &domain = World::getInstance().getDomain();
 
   if (molecules->ListOfMolecules.empty()) {
@@ -75 +78 @@
                 for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                   if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
-                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
+                    distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
                     //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
@@ -135 +138 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        periodicX = *(*iter)->node;
-        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+        periodicX = FullInverseMatrix * (*(*iter)->node); // 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]++)
           for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
             for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-              checkX = Vector(n[0], n[1], n[2]) + periodicX;
-              checkX.MatrixMultiplication(FullMatrix);
+              checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
               for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
                 if ((*MolOtherWalker)->ActiveFlag) {
@@ -157 +158 @@
                       for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                         if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
-                          periodicOtherX = *(*runner)->node;
-                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+                          periodicOtherX = FullInverseMatrix * (*(*runner)->node); // 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]++)
                             for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
                               for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
-                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
-                                checkOtherX.MatrixMultiplication(FullMatrix);
+                                checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
                                 distance = checkX.distance(checkOtherX);
                                 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
@@ -176 +175 @@
             }
       }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
     }
   }
@@ -195 +192 @@
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
-  double *cell_size = World::getInstance().getDomain();
+  Box &domain = World::getInstance().getDomain();
 
   if (molecules->ListOfMolecules.empty()) {
@@ -211 +208 @@
         for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
           if ((*type == NULL) || ((*iter)->type == *type)) {
-            distance = (*iter)->node->PeriodicDistance(*point, cell_size);
+            distance = domain.periodicDistance(*(*iter)->node,*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) ) );
@@ -246 +243 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
@@ -253 +250 @@
         for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
           if ((*type == NULL) || ((*iter)->type == *type)) {
-            periodicX = *(*iter)->node;
-            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+            periodicX = FullInverseMatrix * (*(*iter)->node); // 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]++)
               for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
                 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
-                  checkX.MatrixMultiplication(FullMatrix);
+                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                   distance = checkX.distance(*point);
                   DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
@@ -267 +262 @@
           }
       }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
     }
 
@@ -352 +345 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
@@ -359 +352 @@
         for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
           if ((*type == NULL) || ((*iter)->type == *type)) {
-            periodicX = *(*iter)->node;
-            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
             // go through every range in xyz and get distance
             ShortestDistance = -1.;
@@ -366 +358 @@
               for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
                 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
-                  checkX.MatrixMultiplication(FullMatrix);
+                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                   TriangleIntersectionList Intersections(&checkX,Surface,LC);
                   distance = Intersections.GetSmallestDistance();
@@ -381 +372 @@
           }
       }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
     }
 

src/analysis_correlation.hpp

@@ -26 +26 @@
 
 #include "atom.hpp"
+#include "verbose.hpp"
 
 /****************************************** forward declarations *****************************/

src/analyzer.cpp

@@ -14 +14 @@
 #include "datacreator.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
 #include "periodentafel.hpp"
+#include "verbose.hpp"
 
 // include config.h

src/atom.cpp

@@ -12 +12 @@
 #include "element.hpp"
 #include "lists.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
 #include "vector.hpp"
 #include "World.hpp"
 #include "molecule.hpp"
+#include "Shapes/Shape.hpp"
+
+#include <iomanip>
 
 /************************************* Functions for class atom *************************************/
@@ -112 +114 @@
  * \return true - is inside, false - is not
  */
-bool atom::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
-{
-  return (node->IsInParallelepiped(offset, parallelepiped));
+bool atom::IsInShape(const Shape& shape) const
+{
+  return shape.isInside(*node);
 };
 

src/atom.hpp

@@ -18 +18 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 #include <list>
 #include <vector>
@@ -35 +35 @@
 class World;
 class molecule;
+class Shape;
 
 /********************************************** declarations *******************************/
@@ -66 +67 @@
   double DistanceToVector(const Vector &origin) const;
   double DistanceSquaredToVector(const Vector &origin) const;
-  bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
+  bool IsInShape(const Shape&) const;
 
   // getter and setter

src/atom_graphnodeinfo.cpp

@@ -9 +9 @@
 
 #include "atom_graphnodeinfo.hpp"
-#include "memoryallocator.hpp"
 
 /** Constructor of class GraphNodeInfo.
  */
-GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
+GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(0), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(0) {};
 
 /** Destructor of class GraphNodeInfo.

src/atom_particleinfo.cpp

@@ -9 +9 @@
 
 #include "atom_particleinfo.hpp"
-#include "memoryallocator.hpp"
+
+#include <iostream>
 
 /** Constructor of ParticleInfo.

src/atom_particleinfo.hpp

@@ -18 +18 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
+#include <string>
 
 /****************************************** forward declarations *****************************/

src/atom_trajectoryparticle.hpp

@@ -20 +20 @@
 #include <fstream>
 
+#include <gsl/gsl_inline.h>
 #include <gsl/gsl_randist.h>
 

src/bond.cpp

@@ -7 +7 @@
 #include "Helpers/MemDebug.hpp"
 
+#include "verbose.hpp"
 #include "atom.hpp"
 #include "bond.hpp"

src/bondgraph.cpp

@@ -15 +15 @@
 #include "element.hpp"
 #include "info.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"

src/bondgraph.hpp

@@ -18 +18 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 
 /*********************************************** defines ***********************************/

src/boundary.cpp

@@ -15 +15 @@
 #include "info.hpp"
 #include "linkedcell.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "tesselation.hpp"
@@ -22 +22 @@
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+
+#include <iostream>
+#include <iomanip>
 
 #include<gsl/gsl_poly.h>
@@ -769 +774 @@
   int N[NDIM];
   int n[NDIM];
-  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-  double Rotations[NDIM*NDIM];
-  double *MInverse = InverseMatrix(M);
+  const Matrix &M = World::getInstance().getDomain().getM();
+  Matrix Rotations;
+  const Matrix &MInverse = World::getInstance().getDomain().getMinv();
   Vector AtomTranslations;
   Vector FillerTranslations;
@@ -804 +809 @@
 
   // calculate filler grid in [0,1]^3
-  FillerDistance = Vector(distance[0], distance[1], distance[2]);
-  FillerDistance.InverseMatrixMultiplication(M);
+  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
   for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance[i]);
@@ -818 +822 @@
       for (n[2] = 0; n[2] < N[2]; n[2]++) {
         // calculate position of current grid vector in untransformed box
-        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
-        CurrentPosition.MatrixMultiplication(M);
+        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
         // create molecule random translation vector ...
         for (int i=0;i<NDIM;i++)
@@ -840 +843 @@
             }
 
-            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
-            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
-            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
-            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
-            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
-            Rotations[7] =               sin(phi[1])                                                ;
-            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
-            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
-            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+            Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
+            Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
+            Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
+            Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
+            Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
+            Rotations.set(1,2,              sin(phi[1])                                                    );
+            Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
+            Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
+            Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
           }
 
@@ -854 +857 @@
           Inserter = (*iter)->x;
           if (DoRandomRotation)
-            Inserter.MatrixMultiplication(Rotations);
+            Inserter *= Rotations;
           Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
 
           // check whether inserter is inside box
-          Inserter.MatrixMultiplication(MInverse);
+          Inserter *= MInverse;
           FillIt = true;
           for (int i=0;i<NDIM;i++)
             FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
-          Inserter.MatrixMultiplication(M);
+          Inserter *= M;
 
           // Check whether point is in- or outside
@@ -901 +904 @@
     delete TesselStruct[(*ListRunner)];
   }
-  delete[](M);
-  delete[](MInverse);
 
   return Filling;

src/boundary.hpp

@@ -12 +12 @@
 
 #include <fstream>
-#include <iostream>
+#include <iosfwd>
 
 // STL headers

src/builder.cpp

@@ -60 +60 @@
 #include "config.hpp"
 #include "log.hpp"
-#include "memoryusageobserver.hpp"
 #include "molecule.hpp"
 #include "periodentafel.hpp"

src/config.cpp

@@ -19 +19 @@
 #include "info.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "periodentafel.hpp"
 #include "ThermoStatContainer.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 
 /************************************* Functions for class config ***************************/
@@ -679 +681 @@
   // Unit cell and magnetic field
   ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
-  double * const cell_size = World::getInstance().getDomain();
+  double * cell_size = new double[6];
   cell_size[0] = BoxLength[0];
   cell_size[1] = BoxLength[3];
@@ -686 +688 @@
   cell_size[4] = BoxLength[7];
   cell_size[5] = BoxLength[8];
+  World::getInstance().setDomain(cell_size);
+  delete cell_size;
   //if (1) fprintf(stderr,"\n");
 
@@ -883 +887 @@
 
   ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
-  double * const cell_size = World::getInstance().getDomain();
+  double * cell_size = new double[6];
   cell_size[0] = BoxLength[0];
   cell_size[1] = BoxLength[3];
@@ -890 +894 @@
   cell_size[4] = BoxLength[7];
   cell_size[5] = BoxLength[8];
+  World::getInstance().setDomain(cell_size);
+  delete[] cell_size;
   if (1) fprintf(stderr,"\n");
   config::DoPerturbation = 0;
@@ -1027 +1033 @@
   // bring MaxTypes up to date
   mol->CountElements();
-  const double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   ofstream * const output = new ofstream(filename, ios::out);
   if (output != NULL) {
@@ -1098 +1104 @@
     *output << endl;
     *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
-    *output << cell_size[0] << "\t" << endl;
-    *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
-    *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
+    *output << domain.at(0,0) << "\t" << endl;
+    *output << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
+    *output << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     // FIXME
     *output << endl;

src/datacreator.cpp

@@ -12 +12 @@
 #include "helpers.hpp"
 #include "parser.hpp"
+#include "verbose.hpp"
+
+#include <iomanip>
 
 //=========================== FUNCTIONS============================

src/datacreator.hpp

@@ -10 +10 @@
 using namespace std;
 
-#include <iostream>
+#include <iosfwd>
 
 /****************************************** forward declarations *****************************/

src/element.hpp

@@ -16 +16 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 #include <string>
 

src/ellipsoid.cpp

@@ -21 +21 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 
@@ -34 +35 @@
   Vector helper, RefPoint;
   double distance = -1.;
-  double Matrix[NDIM*NDIM];
+  Matrix Matrix;
   double InverseLength[3];
   double psi,theta,phi; // euler angles in ZX'Z'' convention
@@ -51 +52 @@
   theta = EllipsoidAngle[1];
   phi = EllipsoidAngle[2];
-  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
-  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
-  Matrix[2] = sin(psi)*sin(theta);
-  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
-  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
-  Matrix[5] = -cos(psi)*sin(theta);
-  Matrix[6] = sin(theta)*sin(phi);
-  Matrix[7] = sin(theta)*cos(phi);
-  Matrix[8] = cos(theta);
-  helper.MatrixMultiplication(Matrix);
+  Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
+  Matrix.set(2,0, sin(psi)*sin(theta));
+  Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
+  Matrix.set(2,1, -cos(psi)*sin(theta));
+  Matrix.set(0,2, sin(theta)*sin(phi));
+  Matrix.set(1,2, sin(theta)*cos(phi));
+  Matrix.set(2,2, cos(theta));
+  helper *= Matrix;
   helper.ScaleAll(InverseLength);
   //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
@@ -73 +74 @@
   phi = -EllipsoidAngle[2];
   helper.ScaleAll(EllipsoidLength);
-  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
-  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
-  Matrix[2] = sin(psi)*sin(theta);
-  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
-  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
-  Matrix[5] = -cos(psi)*sin(theta);
-  Matrix[6] = sin(theta)*sin(phi);
-  Matrix[7] = sin(theta)*cos(phi);
-  Matrix[8] = cos(theta);
-  helper.MatrixMultiplication(Matrix);
+  Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
+  Matrix.set(2,0, sin(psi)*sin(theta));
+  Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
+  Matrix.set(2,1, -cos(psi)*sin(theta));
+  Matrix.set(0,2, sin(theta)*sin(phi));
+  Matrix.set(1,2, sin(theta)*cos(phi));
+  Matrix.set(2,2, cos(theta));
+  helper *= Matrix;
   //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;
 

src/errorlogger.cpp

@@ -9 +9 @@
 
 #include <fstream>
+#include <iostream>
 #include "errorlogger.hpp"
 #include "verbose.hpp"

src/errorlogger.hpp

@@ -9 +9 @@
 #define ERRORLOGGER_HPP_
 
-#include <iostream>
+#include <iosfwd>
 
 #include "Patterns/Singleton.hpp"

src/graph.cpp

@@ -13 +13 @@
 #include "config.hpp"
 #include "graph.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"

src/graph.hpp

@@ -27 +27 @@
 class molecule;
 
-class Graph;
 class SubGraph;
 class Node;
@@ -34 +33 @@
 /********************************************** definitions *********************************/
 
-#define NodeMap pair < int, class Node* >
-#define EdgeMap multimap < class Node*, class Edge* >
+typedef std::pair < int, class Node* > NodeMap;
+typedef std::multimap < class Node*, class Edge* > EdgeMap;
 
-#define KeyStack deque<int>
-#define KeySet set<int>
-#define NumberValuePair pair<int, double>
-#define Graph map <KeySet, NumberValuePair, KeyCompare >
-#define GraphPair pair <KeySet, NumberValuePair >
-#define KeySetTestPair pair<KeySet::iterator, bool>
-#define GraphTestPair pair<Graph::iterator, bool>
+typedef std::deque<int> KeyStack;
+typedef std::set<int> KeySet;
+typedef std::pair<int, double> NumberValuePair;
 
-
-/******************************** Some small functions and/or structures **********************************/
-
+// needed for definition of Graph and GraphTestPair
 struct KeyCompare
 {
   bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
 };
+
+typedef std::map <KeySet, NumberValuePair, KeyCompare > Graph;
+typedef std::pair <KeySet, NumberValuePair > GraphPair;
+typedef std::pair<KeySet::iterator, bool> KeySetTestPair;
+typedef std::pair<Graph::iterator, bool> GraphTestPair;
+
+
+/******************************** Some small functions and/or structures **********************************/
 
 //bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)

src/gslvector.cpp

@@ -10 +10 @@
 #include <cassert>
 #include <cmath>
+#include <iostream>
 
 #include "gslvector.hpp"
 #include "defs.hpp"
 #include "vector.hpp"
+#include "VectorContent.hpp"
 
 /** Constructor of class GSLVector.
@@ -69 +71 @@
  */
 void GSLVector::SetFromVector(Vector &v){
-  gsl_vector_memcpy (vector, v.get());
+  gsl_vector_memcpy (vector, v.get()->content);
 }
 

src/gslvector.hpp

@@ -18 +18 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 #include <gsl/gsl_vector.h>
 

src/helpers.cpp

@@ -8 +8 @@
 #include "helpers.hpp"
 #include "Helpers/fast_functions.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryusageobserver.hpp"
+
+#include <iostream>
 
 /********************************************** helpful functions *********************************/
@@ -117 +119 @@
 };
 
-/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
- * \param *symm 6-dim array of unique symmetric matrix components
- * \return allocated NDIM*NDIM array with the symmetric matrix
- */
-double * ReturnFullMatrixforSymmetric(const double * const symm)
-{
-  double *matrix = new double[NDIM * NDIM];
-  matrix[0] = symm[0];
-  matrix[1] = symm[1];
-  matrix[2] = symm[3];
-  matrix[3] = symm[1];
-  matrix[4] = symm[2];
-  matrix[5] = symm[4];
-  matrix[6] = symm[3];
-  matrix[7] = symm[4];
-  matrix[8] = symm[5];
-  return matrix;
-};
-
-/** Calculate the inverse of a 3x3 matrix.
- * \param *matrix NDIM_NDIM array
- */
-double * InverseMatrix( const double * const A)
-{
-  double *B = new double[NDIM * NDIM];
-  double detA = RDET3(A);
-  double detAReci;
-
-  for (int i=0;i<NDIM*NDIM;++i)
-    B[i] = 0.;
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-  }
-  return B;
-};
-
-
-
 /** Comparison function for GSL heapsort on distances in two molecules.
  * \param *a

src/helpers.hpp

@@ -20 +20 @@
 #include "defs.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 
 /********************************************** definitions *********************************/
@@ -52 +51 @@
 bool IsValidNumber( const char *string);
 int CompareDoubles (const void * a, const void * b);
-double * ReturnFullMatrixforSymmetric(const double * const cell_size);
-double * InverseMatrix(const double * const A);
 void performCriticalExit();
 

src/joiner.cpp

@@ -14 +14 @@
 #include "datacreator.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
 #include "periodentafel.hpp"
+#include "verbose.hpp"
 
 //============================== MAIN =============================

src/linkedcell.cpp

@@ -10 +10 @@
 #include "helpers.hpp"
 #include "linkedcell.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"

src/logger.cpp

@@ -9 +9 @@
 
 #include <fstream>
+#include <iostream>
 #include "logger.hpp"
 #include "verbose.hpp"

src/logger.hpp

@@ -9 +9 @@
 #define LOGGER_HPP_
 
-#include <iostream>
+#include <iosfwd>
 
 #include "Patterns/Singleton.hpp"

src/molecule.cpp

@@ -5 +5 @@
  */
 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
 #include "Helpers/MemDebug.hpp"
 
 #include <cstring>
 #include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_heapsort.h>
 
 #include "World.hpp"
@@ -22 +30 @@
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
+
 #include "periodentafel.hpp"
 #include "stackclass.hpp"
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 #include "Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
@@ -284 +294 @@
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
-  double *matrix = NULL;
+  const Matrix &matrix =  World::getInstance().getDomain().getM();
   bond *Binder = NULL;
-  double * const cell_size = World::getInstance().getDomain();
 
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -307 +316 @@
       } // (signs are correct, was tested!)
     }
-    matrix = ReturnFullMatrixforSymmetric(cell_size);
-    Orthovector1.MatrixMultiplication(matrix);
+    Orthovector1 *= matrix;
     InBondvector -= Orthovector1; // subtract just the additional translation
-    delete[](matrix);
     bondlength = InBondvector.Norm();
 //    Log() << Verbose(4) << "Corrected InBondvector is now: ";
@@ -541 +548 @@
       break;
   }
-  delete[](matrix);
 
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
@@ -660 +666 @@
  * @param three vectors forming the matrix that defines the shape of the parallelpiped
  */
-molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
+molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
   molecule *copy = World::getInstance().createMolecule();
 
-  ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
+  BOOST_FOREACH(atom *iter,atoms){
+    if(iter->IsInShape(region)){
+      copy->AddCopyAtom(iter);
+    }
+  }
 
   //TODO: copy->BuildInducedSubgraph(this);
@@ -750 +760 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
-  double * const cell_size = World::getInstance().getDomain();
-  cell_size[0] = dim->at(0);
-  cell_size[1] = 0.;
-  cell_size[2] = dim->at(1);
-  cell_size[3] = 0.;
-  cell_size[4] = 0.;
-  cell_size[5] = dim->at(2);
+  Matrix domain;
+  for(int i =0; i<NDIM;++i)
+    domain.at(i,i) = dim->at(i);
+  World::getInstance().setDomain(domain);
 };
 
@@ -849 +856 @@
 bool molecule::CheckBounds(const Vector *x) const
 {
-  double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   bool result = true;
-  int j =-1;
   for (int i=0;i<NDIM;i++) {
-    j += i+1;
-    result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
+    result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
   }
   //return result;

src/molecule.hpp

@@ -7 +7 @@
 #define MOLECULES_HPP_
 
-using namespace std;
-
 /*********************************************** includes ***********************************/
 
-// GSL headers
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_heapsort.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_randist.h>
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 
 //// STL headers
@@ -31 +24 @@
 #include "types.hpp"
 #include "graph.hpp"
-#include "stackclass.hpp"
 #include "tesselation.hpp"
 #include "Patterns/Observer.hpp"
@@ -53 +45 @@
 class periodentafel;
 class Vector;
+class Shape;
+template <class> class StackClass;
 
 /******************************** Some definitions for easier reading **********************************/
@@ -316 +310 @@
 
   molecule *CopyMolecule();
-  molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;
+  molecule* CopyMoleculeFromSubRegion(const Shape&) const;
 
   /// Fragment molecule by two different approaches:

src/molecule_dynamics.cpp

@@ -13 +13 @@
 #include "element.hpp"
 #include "info.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "parser.hpp"
 #include "Plane.hpp"
 #include "ThermoStatContainer.hpp"
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
 
 /************************************* Functions for class molecule *********************************/

src/molecule_fragmentation.cpp

@@ -17 +17 @@
 #include "helpers.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "periodentafel.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include "stackclass.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -1716 +1719 @@
   atom *Walker = NULL;
   atom *OtherWalker = NULL;
-  double * const cell_size = World::getInstance().getDomain();
-  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
+  Matrix matrix = World::getInstance().getDomain().getM();
   enum Shading *ColorList = NULL;
   double tmp;
@@ -1757 +1759 @@
           Translationvector[i] = (tmp < 0) ? +1. : -1.;
       }
-      Translationvector.MatrixMultiplication(matrix);
+      Translationvector *= matrix;
       //Log() << Verbose(3) << "Translation vector is ";
       Log() << Verbose(0) << Translationvector <<  endl;
@@ -1788 +1790 @@
   delete(AtomStack);
   delete[](ColorList);
-  delete[](matrix);
   DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
 };

src/molecule_geometry.cpp

@@ -5 +5 @@
  *      Author: heber
  */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 
 #include "Helpers/MemDebug.hpp"
@@ -14 +18 @@
 #include "helpers.hpp"
 #include "leastsquaremin.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 #include <boost/foreach.hpp>
+
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_multimin.h>
 
 
@@ -33 +42 @@
   const Vector *Center = DetermineCenterOfAll();
   const Vector *CenterBox = DetermineCenterOfBox();
-  double * const cell_size = World::getInstance().getDomain();
-  double *M = ReturnFullMatrixforSymmetric(cell_size);
-  double *Minv = InverseMatrix(M);
+  Box &domain = World::getInstance().getDomain();
 
   // go through all atoms
   ActOnAllVectors( &Vector::SubtractVector, *Center);
   ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
-  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
-
-  delete[](M);
-  delete[](Minv);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter->node = domain.WrapPeriodically(*iter->node);
+  }
+
   delete(Center);
   delete(CenterBox);
@@ -56 +63 @@
 {
   bool status = true;
-  double * const cell_size = World::getInstance().getDomain();
-  double *M = ReturnFullMatrixforSymmetric(cell_size);
-  double *Minv = InverseMatrix(M);
+  Box &domain = World::getInstance().getDomain();
 
   // go through all atoms
-  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
-
-  delete[](M);
-  delete[](Minv);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter->node = domain.WrapPeriodically(*iter->node);
+  }
+
   return status;
 };
@@ -153 +158 @@
 {
   Vector *a = new Vector(0.5,0.5,0.5);
-
-  const double *cell_size = World::getInstance().getDomain();
-  double *M = ReturnFullMatrixforSymmetric(cell_size);
-  a->MatrixMultiplication(M);
-  delete[](M);
-
+  const Matrix &M = World::getInstance().getDomain().getM();
+  (*a) *= M;
   return a;
 };
@@ -244 +245 @@
 void molecule::TranslatePeriodically(const Vector *trans)
 {
-  double * const cell_size = World::getInstance().getDomain();
-  double *M = ReturnFullMatrixforSymmetric(cell_size);
-  double *Minv = InverseMatrix(M);
+  Box &domain = World::getInstance().getDomain();
 
   // go through all atoms
   ActOnAllVectors( &Vector::AddVector, *trans);
-  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
-
-  delete[](M);
-  delete[](Minv);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter->node = domain.WrapPeriodically(*iter->node);
+  }
+
 };
 
@@ -264 +263 @@
   OBSERVE;
   Plane p(*n,0);
-  BOOST_FOREACH( atom* iter, atoms ){
+  BOOST_FOREACH(atom* iter, atoms ){
     (*iter->node) = p.mirrorVector(*iter->node);
   }
@@ -274 +273 @@
 void molecule::DeterminePeriodicCenter(Vector &center)
 {
-  double * const cell_size = World::getInstance().getDomain();
-  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
-  double *inversematrix = InverseMatrix(matrix);
+  const Matrix &matrix = World::getInstance().getDomain().getM();
+  const Matrix &inversematrix = World::getInstance().getDomain().getM();
   double tmp;
   bool flag;
@@ -288 +286 @@
       if ((*iter)->type->Z != 1) {
 #endif
-        Testvector = (*iter)->x;
-        Testvector.MatrixMultiplication(inversematrix);
+        Testvector = inversematrix * (*iter)->x;
         Translationvector.Zero();
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
@@ -306 +303 @@
         }
         Testvector += Translationvector;
-        Testvector.MatrixMultiplication(matrix);
+        Testvector *= matrix;
         Center += Testvector;
         Log() << Verbose(1) << "vector is: " << Testvector << endl;
@@ -313 +310 @@
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
           if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
-            Testvector = (*Runner)->GetOtherAtom((*iter))->x;
-            Testvector.MatrixMultiplication(inversematrix);
+            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
             Testvector += Translationvector;
-            Testvector.MatrixMultiplication(matrix);
+            Testvector *= matrix;
             Center += Testvector;
             Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
@@ -325 +321 @@
     }
   } while (!flag);
-  delete[](matrix);
-  delete[](inversematrix);
 
   Center.Scale(1./static_cast<double>(getAtomCount()));
@@ -388 +382 @@
     DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
     // the eigenvectors specify the transformation matrix
-    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
+    Matrix M = Matrix(evec->data);
+    BOOST_FOREACH(atom* iter, atoms){
+      (*iter->node) *= M;
+    }
     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
 

src/molecule_graph.cpp

@@ -5 +5 @@
  *      Author: heber
  */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 
 #include "Helpers/MemDebug.hpp"
@@ -18 +22 @@
 #include "linkedcell.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
 #include "Helpers/fast_functions.hpp"
 #include "Helpers/Assert.hpp"
-
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include "stackclass.hpp"
 
 struct BFSAccounting
@@ -121 +127 @@
   LinkedCell *LC = NULL;
   bool free_BG = false;
-  double * const cell_size = World::getInstance().getDomain();
+  Box &domain = World::getInstance().getDomain();
 
   if (BG == NULL) {
@@ -178 +184 @@
                           //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 = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);
+                          const double distance = domain.periodicDistanceSquared(OtherWalker->x,Walker->x);
                           const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
 //                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
@@ -579 +585 @@
     // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
     DoLog(0) && (Log() << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << DFS.CurrentGraphNr << "." << endl);
-    LeafWalker->Leaf->Output((ofstream *)&cout);
+    LeafWalker->Leaf->Output((ofstream *)&(Log() << Verbose(0)));
     DoLog(0) && (Log() << Verbose(0) << endl);
 

src/molecule_pointcloud.cpp

@@ -11 +11 @@
 #include "config.hpp"
 #include "info.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 

src/moleculelist.cpp

@@ -5 +5 @@
  */
 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
 #include "Helpers/MemDebug.hpp"
 
 #include <cstring>
+
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_heapsort.h>
 
 #include "World.hpp"
@@ -19 +26 @@
 #include "linkedcell.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
 #include "periodentafel.hpp"
 #include "Helpers/Assert.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include "stackclass.hpp"
 
 #include "Helpers/Assert.hpp"
@@ -631 +641 @@
   int FragmentCounter = 0;
   ofstream output;
-  double cell_size_backup[6];
-  double * const cell_size = World::getInstance().getDomain();
-
-  // backup cell_size
-  for (int i=0;i<6;i++)
-    cell_size_backup[i] = cell_size[i];
+  Matrix cell_size = World::getInstance().getDomain().getM();
+  Matrix cell_size_backup = cell_size;
+
   // store the fragments as config and as xyz
   for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
@@ -674 +681 @@
     (*ListRunner)->CenterEdge(&BoxDimension);
     (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
-    int j = -1;
     for (int k = 0; k < NDIM; k++) {
-      j += k + 1;
       BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
-      cell_size[j] = BoxDimension[k] * 2.;
-    }
+      cell_size.at(k,k) = BoxDimension[k] * 2.;
+    }
+    World::getInstance().setDomain(cell_size);
     (*ListRunner)->Translate(&BoxDimension);
 
@@ -725 +731 @@
 
   // restore cell_size
-  for (int i=0;i<6;i++)
-    cell_size[i] = cell_size_backup[i];
+  World::getInstance().setDomain(cell_size_backup);
 
   return result;
@@ -887 +892 @@
   // center at set box dimensions
   mol->CenterEdge(&center);
-  World::getInstance().getDomain()[0] = center[0];
-  World::getInstance().getDomain()[1] = 0;
-  World::getInstance().getDomain()[2] = center[1];
-  World::getInstance().getDomain()[3] = 0;
-  World::getInstance().getDomain()[4] = 0;
-  World::getInstance().getDomain()[5] = center[2];
+  Matrix domain;
+  for(int i =0;i<NDIM;++i)
+    domain.at(i,i) = center[i];
+  World::getInstance().setDomain(domain);
   insert(mol);
 }

src/parser.cpp

@@ -12 +12 @@
 
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
+#include "verbose.hpp"
 
 // include config.h

src/periodentafel.hpp

@@ -9 +9 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 #include <map>
-#include <iterator>
 
 #include "unittests/periodentafelTest.hpp"

src/stackclass.hpp

@@ -12 +12 @@
 
 #include "verbose.hpp"
-#include "memoryallocator.hpp"
+#include "log.hpp"
 
 /****************************************** forward declarations *****************************/

src/tesselation.cpp

@@ -9 +9 @@
 
 #include <fstream>
+#include <iomanip>
 
 #include "helpers.hpp"

src/tesselationhelpers.cpp

@@ -21 +21 @@
 #include "verbose.hpp"
 #include "Plane.hpp"
-
-double DetGet(gsl_matrix * const A, const int inPlace)
-{
-        Info FunctionInfo(__func__);
-  /*
-  inPlace = 1 => A is replaced with the LU decomposed copy.
-  inPlace = 0 => A is retained, and a copy is used for LU.
-  */
-
-  double det;
-  int signum;
-  gsl_permutation *p = gsl_permutation_alloc(A->size1);
-  gsl_matrix *tmpA=0;
-
-  if (inPlace)
-  tmpA = A;
-  else {
-  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
-  gsl_matrix_memcpy(tmpA , A);
-  }
-
-
-  gsl_linalg_LU_decomp(tmpA , p , &signum);
-  det = gsl_linalg_LU_det(tmpA , signum);
-  gsl_permutation_free(p);
-  if (! inPlace)
-  gsl_matrix_free(tmpA);
-
-  return det;
-};
+#include "Matrix.hpp"
 
 void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
 {
         Info FunctionInfo(__func__);
-  gsl_matrix *A = gsl_matrix_calloc(3,3);
+  Matrix mat;
   double m11, m12, m13, m14;
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a[i]);
-    gsl_matrix_set(A, i, 1, b[i]);
-    gsl_matrix_set(A, i, 2, c[i]);
-  }
-  m11 = DetGet(A, 1);
+    mat.set(i, 0, a[i]);
+    mat.set(i, 1, b[i]);
+    mat.set(i, 2, c[i]);
+  }
+  m11 = mat.determinant();
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
-    gsl_matrix_set(A, i, 1, b[i]);
-    gsl_matrix_set(A, i, 2, c[i]);
-  }
-  m12 = DetGet(A, 1);
+    mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    mat.set(i, 1, b[i]);
+    mat.set(i, 2, c[i]);
+  }
+  m12 = mat.determinant();
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
-    gsl_matrix_set(A, i, 1, a[i]);
-    gsl_matrix_set(A, i, 2, c[i]);
-  }
-  m13 = DetGet(A, 1);
+    mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    mat.set(i, 1, a[i]);
+    mat.set(i, 2, c[i]);
+  }
+  m13 = mat.determinant();
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
-    gsl_matrix_set(A, i, 1, a[i]);
-    gsl_matrix_set(A, i, 2, b[i]);
-  }
-  m14 = DetGet(A, 1);
+    mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    mat.set(i, 1, a[i]);
+    mat.set(i, 2, b[i]);
+  }
+  m14 = mat.determinant();
 
   if (fabs(m11) < MYEPSILON)
@@ -95 +66 @@
   if (fabs(a.distance(*center) - RADIUS) > MYEPSILON)
     DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.distance(*center) - RADIUS) << " from a than RADIUS." << endl);
-
-  gsl_matrix_free(A);
 };
 
@@ -248 +217 @@
   Vector x3;
   Vector x4;
-};
-
-/**
- * Intersection calculation function.
- *
- * @param x to find the result for
- * @param function parameter
- */
-double MinIntersectDistance(const gsl_vector * x, void *params)
-{
-        Info FunctionInfo(__func__);
-  double retval = 0;
-  struct Intersection *I = (struct Intersection *)params;
-  Vector intersection;
-  for (int i=0;i<NDIM;i++)
-    intersection[i] = gsl_vector_get(x, i);
-
-  Vector SideA = I->x1 -I->x2 ;
-  Vector HeightA = intersection - I->x1;
-  HeightA.ProjectOntoPlane(SideA);
-
-  Vector SideB = I->x3 - I->x4;
-  Vector HeightB = intersection - I->x3;
-  HeightB.ProjectOntoPlane(SideB);
-
-  retval = HeightA.ScalarProduct(HeightA) + HeightB.ScalarProduct(HeightB);
-  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
-
-  return retval;
-};
-
-
-/**
- * Calculates whether there is an intersection between two lines. The first line
- * always goes through point 1 and point 2 and the second line is given by the
- * connection between point 4 and point 5.
- *
- * @param point 1 of line 1
- * @param point 2 of line 1
- * @param point 1 of line 2
- * @param point 2 of line 2
- *
- * @return true if there is an intersection between the given lines, false otherwise
- */
-bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
-{
-        Info FunctionInfo(__func__);
-  bool result;
-
-  struct Intersection par;
-    par.x1 = point1;
-    par.x2 = point2;
-    par.x3 = point3;
-    par.x4 = point4;
-
-    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
-    gsl_multimin_fminimizer *s = NULL;
-    gsl_vector *ss, *x;
-    gsl_multimin_function minexFunction;
-
-    size_t iter = 0;
-    int status;
-    double size;
-
-    /* Starting point */
-    x = gsl_vector_alloc(NDIM);
-    gsl_vector_set(x, 0, point1[0]);
-    gsl_vector_set(x, 1, point1[1]);
-    gsl_vector_set(x, 2, point1[2]);
-
-    /* Set initial step sizes to 1 */
-    ss = gsl_vector_alloc(NDIM);
-    gsl_vector_set_all(ss, 1.0);
-
-    /* Initialize method and iterate */
-    minexFunction.n = NDIM;
-    minexFunction.f = &MinIntersectDistance;
-    minexFunction.params = (void *)&par;
-
-    s = gsl_multimin_fminimizer_alloc(T, NDIM);
-    gsl_multimin_fminimizer_set(s, &minexFunction, x, ss);
-
-    do {
-        iter++;
-        status = gsl_multimin_fminimizer_iterate(s);
-
-        if (status) {
-          break;
-        }
-
-        size = gsl_multimin_fminimizer_size(s);
-        status = gsl_multimin_test_size(size, 1e-2);
-
-        if (status == GSL_SUCCESS) {
-          DoLog(1) && (Log() << Verbose(1) << "converged to minimum" <<  endl);
-        }
-    } while (status == GSL_CONTINUE && iter < 100);
-
-    // check whether intersection is in between or not
-  Vector intersection;
-  double t1, t2;
-  for (int i = 0; i < NDIM; i++) {
-    intersection[i] = gsl_vector_get(s->x, i);
-  }
-
-  Vector SideA = par.x2 - par.x1;
-  Vector HeightA = intersection - par.x1;
-
-  t1 = HeightA.ScalarProduct(SideA)/SideA.ScalarProduct(SideA);
-
-  Vector SideB = par.x4 - par.x3;
-  Vector HeightB = intersection - par.x3;
-
-  t2 = HeightB.ScalarProduct(SideB)/SideB.ScalarProduct(SideB);
-
-  Log() << Verbose(1) << "Intersection " << intersection << " is at "
-    << t1 << " for (" << point1 << "," << point2 << ") and at "
-    << t2 << " for (" << point3 << "," << point4 << "): ";
-
-  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
-    DoLog(1) && (Log() << Verbose(1) << "true intersection." << endl);
-    result = true;
-  } else {
-    DoLog(1) && (Log() << Verbose(1) << "intersection out of region of interest." << endl);
-    result = false;
-  }
-
-  // free minimizer stuff
-    gsl_vector_free(x);
-    gsl_vector_free(ss);
-    gsl_multimin_fminimizer_free(s);
-
-  return result;
 };
 

src/tesselationhelpers.hpp

@@ -20 +20 @@
 #endif
 
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_permutation.h>
-#include <gsl/gsl_vector.h>
-
-#include <iostream>
+#include <iosfwd>
 
 #include "defs.hpp"
@@ -47 +41 @@
 /********************************************** declarations *******************************/
 
-double DetGet(gsl_matrix * const A, const int inPlace);
 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);
 double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection);
-double MinIntersectDistance(const gsl_vector * x, void *params);
-bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4);
 double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d);
 double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C);

src/triangleintersectionlist.cpp

@@ -18 +18 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "verbose.hpp"
 
 /** Constructor for class TriangleIntersectionList.

src/unittests/ActOnAllUnitTest.cpp

@@ -14 +14 @@
 #include "../test/ActOnAlltest.hpp"
 #include "ActOnAllUnitTest.hpp"
-#include "memoryallocator.hpp"
 #include "vector.hpp"
 

src/unittests/LineUnittest.cpp

@@ -17 +17 @@
 
 #include <iostream>
+#include <cmath>
 
 using namespace std;
@@ -352 +353 @@
   CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
 }
+
+void LineUnittest::sphereIntersectionTest(){
+  {
+    std::vector<Vector> res = la1->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT(testDirection(res[0],e1));
+    CPPUNIT_ASSERT(testDirection(res[1],e1));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+
+  {
+    std::vector<Vector> res = la2->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT(testDirection(res[0],e2));
+    CPPUNIT_ASSERT(testDirection(res[1],e2));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+
+  {
+    std::vector<Vector> res = la3->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT(testDirection(res[0],e3));
+    CPPUNIT_ASSERT(testDirection(res[1],e3));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+
+  {
+    std::vector<Vector> res = lp1->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT((res[0]==e1) || (res[0]==e2));
+    CPPUNIT_ASSERT((res[1]==e1) || (res[1]==e2));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+
+  {
+    std::vector<Vector> res = lp2->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT((res[0]==e2) || (res[0]==e3));
+    CPPUNIT_ASSERT((res[1]==e2) || (res[1]==e3));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+
+  {
+    std::vector<Vector> res = lp3->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT((res[0]==e3) || (res[0]==e1));
+    CPPUNIT_ASSERT((res[1]==e3) || (res[1]==e1));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+}

src/unittests/LineUnittest.hpp

@@ -22 +22 @@
   CPPUNIT_TEST ( intersectionTest );
   CPPUNIT_TEST ( rotationTest );
+  CPPUNIT_TEST ( sphereIntersectionTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -33 +34 @@
   void intersectionTest();
   void rotationTest();
+  void sphereIntersectionTest();
 
 private:

src/unittests/Makefile.am

@@ -31 +31 @@
   LogUnitTest \
   manipulateAtomsTest \
-  MemoryUsageObserverUnitTest \
-  MemoryAllocatorUnitTest \
+  MatrixUnittest \
   MoleculeDescriptorTest \
   ObserverTest \
@@ -38 +37 @@
   periodentafelTest \
   PlaneUnittest \
+  ShapeUnittest \
   SingletonTest \
   StackClassUnitTest \
@@ -75 +75 @@
   listofbondsunittest.cpp \
   logunittest.cpp \
+  MatrixUnittest.cpp \
   manipulateAtomsTest.cpp \
-  memoryallocatorunittest.cpp  \
-  memoryusageobserverunittest.cpp \
   MoleculeDescriptorTest.cpp \
   ObserverTest.cpp \
@@ -83 +82 @@
   periodentafelTest.cpp \
   PlaneUnittest.cpp \
+  ShapeUnittest.cpp \
   SingletonTest.cpp \
   stackclassunittest.cpp \
@@ -112 +112 @@
   logunittest.hpp \
   manipulateAtomsTest.hpp \
-  memoryallocatorunittest.hpp  \
-  memoryusageobserverunittest.hpp \
+  MatrixUnittest.hpp \
   MoleculeDescriptorTest.hpp \
   periodentafelTest.hpp \
@@ -189 +188 @@
 manipulateAtomsTest_LDADD = ${ALLLIBS}
 
-MemoryAllocatorUnitTest_SOURCES = UnitTestMain.cpp ../memoryallocator.hpp ../memoryallocator.cpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp memoryallocatorunittest.cpp memoryallocatorunittest.hpp 
-MemoryAllocatorUnitTest_LDADD = ${ALLLIBS}
-
-MemoryUsageObserverUnitTest_SOURCES = UnitTestMain.cpp ../memoryallocator.hpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp memoryusageobserverunittest.cpp memoryusageobserverunittest.hpp 
-MemoryUsageObserverUnitTest_LDADD = ${ALLLIBS}
+MatrixUnittest_SOURCES = UnitTestMain.cpp MatrixUnittest.cpp MatrixUnittest.hpp
+MatrixUnittest_LDADD = ${ALLLIBS}
 
 MoleculeDescriptorTest_SOURCES = UnitTestMain.cpp MoleculeDescriptorTest.cpp MoleculeDescriptorTest.hpp
@@ -210 +206 @@
 PlaneUnittest_LDADD = ${ALLLIBS}
 
+ShapeUnittest_SOURCES = UnitTestMain.cpp ShapeUnittest.cpp ShapeUnittest.hpp
+ShapeUnittest_LDADD = ${ALLLIBS}
+
 SingletonTest_SOURCES = UnitTestMain.cpp SingletonTest.cpp SingletonTest.hpp
 SingletonTest_LDADD = ${ALLLIBS} $(BOOST_LIB) ${BOOST_THREAD_LIB}
@@ -225 +224 @@
 Tesselation_InOutsideUnitTest_LDADD = ${ALLLIBS}
 
-TestRunner_SOURCES = TestRunnerMain.cpp ../memoryallocator.hpp ../memoryallocator.cpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp $(TESTSOURCES) $(TESTHEADERS)
+TestRunner_SOURCES = TestRunnerMain.cpp $(TESTSOURCES) $(TESTHEADERS)
 TestRunner_LDADD = ${ALLLIBS}
 

src/unittests/PlaneUnittest.cpp

@@ -11 +11 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+
+#include <cmath>
 
 #ifdef HAVE_TESTRUNNER

src/unittests/linearsystemofequationsunittest.cpp

@@ -12 +12 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+#include <cmath>
 
 #include "linearsystemofequationsunittest.hpp"

src/unittests/tesselation_insideoutsideunittest.cpp

@@ -17 +17 @@
 #include "tesselation.hpp"
 #include "tesselation_insideoutsideunittest.hpp"
+#include "verbose.hpp"
 
 #ifdef HAVE_TESTRUNNER

src/unittests/vectorunittest.cpp

@@ -15 +15 @@
 #include "defs.hpp"
 #include "log.hpp"
-#include "memoryusageobserver.hpp"
 #include "vector.hpp"
 #include "vector_ops.hpp"
@@ -21 +20 @@
 #include "Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
+#include "Matrix.hpp"
 
 #ifdef HAVE_TESTRUNNER
@@ -214 +214 @@
   CPPUNIT_ASSERT(testVector.ScalarProduct(three) < MYEPSILON);
 }
-
-
-/**
- * UnitTest for Vector::IsInParallelepiped().
- */
-void VectorTest::IsInParallelepipedTest()
-{
-  double parallelepiped[NDIM*NDIM];
-  parallelepiped[0] = 1;
-  parallelepiped[1] = 0;
-  parallelepiped[2] = 0;
-  parallelepiped[3] = 0;
-  parallelepiped[4] = 1;
-  parallelepiped[5] = 0;
-  parallelepiped[6] = 0;
-  parallelepiped[7] = 0;
-  parallelepiped[8] = 1;
-
-  fixture = zero;
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(2.5,2.5,2.5);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(1.,1.,1.);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(3.5,3.5,3.5);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(2.,2.,2.);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(2.,3.,2.);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(-2.,2.,-1.);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-}
-

src/unittests/vectorunittest.hpp

@@ -27 +27 @@
     CPPUNIT_TEST ( ProjectionTest );
     CPPUNIT_TEST ( NormalsTest );
-    CPPUNIT_TEST ( IsInParallelepipedTest );
     CPPUNIT_TEST_SUITE_END();
 
@@ -45 +44 @@
     void LineIntersectionTest();
     void VectorRotationTest();
-    void IsInParallelepipedTest();
 
 private:

src/vector.cpp

@@ -8 +8 @@
 
 #include "vector.hpp"
+#include "VectorContent.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
 #include "Helpers/Assert.hpp"
 #include "Helpers/fast_functions.hpp"
+#include "Exceptions/MathException.hpp"
 
 #include <iostream>
+#include <gsl/gsl_blas.h>
+
 
 using namespace std;
@@ -24 +28 @@
 Vector::Vector()
 {
-  content = gsl_vector_calloc (NDIM);
+  content = new VectorContent();
 };
 
@@ -33 +37 @@
 Vector::Vector(const Vector& src)
 {
-  content = gsl_vector_alloc(NDIM);
-  gsl_vector_set(content,0,src[0]);
-  gsl_vector_set(content,1,src[1]);
-  gsl_vector_set(content,2,src[2]);
+  content = new VectorContent();
+  gsl_vector_memcpy(content->content, src.content->content);
 }
 
@@ -43 +45 @@
 Vector::Vector(const double x1, const double x2, const double x3)
 {
-  content = gsl_vector_alloc(NDIM);
-  gsl_vector_set(content,0,x1);
-  gsl_vector_set(content,1,x2);
-  gsl_vector_set(content,2,x3);
-};
+  content = new VectorContent();
+  gsl_vector_set(content->content,0,x1);
+  gsl_vector_set(content->content,1,x2);
+  gsl_vector_set(content->content,2,x3);
+};
+
+Vector::Vector(VectorContent *_content) :
+  content(_content)
+{}
 
 /**
@@ -55 +61 @@
   // check for self assignment
   if(&src!=this){
-    gsl_vector_set(content,0,src[0]);
-    gsl_vector_set(content,1,src[1]);
-    gsl_vector_set(content,2,src[2]);
+    gsl_vector_memcpy(content->content, src.content->content);
   }
   return *this;
@@ -65 +69 @@
  */
 Vector::~Vector() {
-  gsl_vector_free(content);
+  delete content;
 };
 
@@ -94 +98 @@
 }
 
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |\f$
- */
-double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
-{
-  double res = distance(y), tmp, matrix[NDIM*NDIM];
-    Vector Shiftedy, TranslationVector;
-    int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
-    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-    for (N[0]=-1;N[0]<=1;N[0]++)
-      for (N[1]=-1;N[1]<=1;N[1]++)
-        for (N[2]=-1;N[2]<=1;N[2]++) {
-          // create the translation vector
-          TranslationVector.Zero();
-          for (int i=NDIM;i--;)
-            TranslationVector[i] = (double)N[i];
-          TranslationVector.MatrixMultiplication(matrix);
-          // add onto the original vector to compare with
-          Shiftedy = y + TranslationVector;
-          // get distance and compare with minimum so far
-          tmp = distance(Shiftedy);
-          if (tmp < res) res = tmp;
-        }
-    return (res);
-};
-
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |^2\f$
- */
-double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
-{
-  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
-    Vector Shiftedy, TranslationVector;
-    int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
-    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-    for (N[0]=-1;N[0]<=1;N[0]++)
-      for (N[1]=-1;N[1]<=1;N[1]++)
-        for (N[2]=-1;N[2]<=1;N[2]++) {
-          // create the translation vector
-          TranslationVector.Zero();
-          for (int i=NDIM;i--;)
-            TranslationVector[i] = (double)N[i];
-          TranslationVector.MatrixMultiplication(matrix);
-          // add onto the original vector to compare with
-          Shiftedy = y + TranslationVector;
-          // get distance and compare with minimum so far
-          tmp = DistanceSquared(Shiftedy);
-          if (tmp < res) res = tmp;
-        }
-    return (res);
-};
-
-/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
- * \param *out ofstream for debugging messages
- * Tries to translate a vector into each adjacent neighbouring cell.
- */
-void Vector::KeepPeriodic(const double * const matrix)
-{
-  //  int N[NDIM];
-  //  bool flag = false;
-    //vector Shifted, TranslationVector;
-  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
-  //  Log() << Verbose(2) << "Vector is: ";
-  //  Output(out);
-  //  Log() << Verbose(0) << endl;
-    InverseMatrixMultiplication(matrix);
-    for(int i=NDIM;i--;) { // correct periodically
-      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
-        at(i) += ceil(at(i));
-      } else {
-        at(i) -= floor(at(i));
-      }
-    }
-    MatrixMultiplication(matrix);
-  //  Log() << Verbose(2) << "New corrected vector is: ";
-  //  Output(out);
-  //  Log() << Verbose(0) << endl;
-  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
-};
-
 /** Calculates scalar product between this and another vector.
  * \param *y array to second vector
@@ -203 +105 @@
 {
   double res = 0.;
-  for (int i=NDIM;i--;)
-    res += at(i)*y[i];
+  gsl_blas_ddot(content->content, y.content->content, &res);
   return (res);
 };
@@ -369 +270 @@
 double& Vector::operator[](size_t i){
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return *gsl_vector_ptr (content, i);
+  return *gsl_vector_ptr (content->content, i);
 }
 
 const double& Vector::operator[](size_t i) const{
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return *gsl_vector_ptr (content, i);
+  return *gsl_vector_ptr (content->content, i);
 }
 
@@ -385 +286 @@
 }
 
-gsl_vector* Vector::get(){
+VectorContent* Vector::get(){
   return content;
 }
@@ -504 +405 @@
 };
 
+void Vector::ScaleAll(const Vector &factor){
+  gsl_vector_mul(content->content, factor.content->content);
+}
 
 
 void Vector::Scale(const double factor)
 {
-  for (int i=NDIM;i--;)
-    at(i) *= factor;
-};
-
-/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
- * \param *M matrix of box
- * \param *Minv inverse matrix
- */
-void Vector::WrapPeriodically(const double * const M, const double * const Minv)
-{
-  MatrixMultiplication(Minv);
-  // truncate to [0,1] for each axis
-  for (int i=0;i<NDIM;i++) {
-    //at(i) += 0.5;  // set to center of box
-    while (at(i) >= 1.)
-      at(i) -= 1.;
-    while (at(i) < 0.)
-      at(i) += 1.;
-  }
-  MatrixMultiplication(M);
+  gsl_vector_scale(content->content,factor);
 };
 
@@ -546 +431 @@
   return make_pair(res,helper);
 }
-
-/** Do a matrix multiplication.
- * \param *matrix NDIM_NDIM array
- */
-void Vector::MatrixMultiplication(const double * const M)
-{
-  Vector tmp;
-  // do the matrix multiplication
-  for(int i=NDIM;i--;)
-    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
-
-  (*this) = tmp;
-};
-
-/** Do a matrix multiplication with the \a *A' inverse.
- * \param *matrix NDIM_NDIM array
- */
-bool Vector::InverseMatrixMultiplication(const double * const A)
-{
-  double B[NDIM*NDIM];
-  double detA = RDET3(A);
-  double detAReci;
-
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-
-    MatrixMultiplication(B);
-
-    return true;
-  } else {
-    return false;
-  }
-};
-
 
 /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
@@ -679 +520 @@
 void Vector::AddVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    at(i) += y[i];
+  gsl_vector_add(content->content, y.content->content);
 }
 
@@ -688 +528 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    at(i) -= y[i];
-}
-
-/**
- * Checks whether this vector is within the parallelepiped defined by the given three vectors and
- * their offset.
- *
- * @param offest for the origin of the parallelepiped
- * @param three vectors forming the matrix that defines the shape of the parallelpiped
- */
-bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
-{
-  Vector a = (*this)-offset;
-  a.InverseMatrixMultiplication(parallelepiped);
-  bool isInside = true;
-
-  for (int i=NDIM;i--;)
-    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
-
-  return isInside;
+  gsl_vector_sub(content->content, y.content->content);
 }
 

src/vector.hpp

@@ -11 +11 @@
 #endif
 
-#include <iostream>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_multimin.h>
+#include <iosfwd>
 
 #include <memory>
@@ -24 +22 @@
 
 class Vector;
+class Matrix;
+struct VectorContent;
 
 typedef std::vector<Vector> pointset;
@@ -31 +31 @@
  */
 class Vector : public Space{
+  friend Vector operator*(const Matrix&,const Vector&);
+  friend class Matrix;
 public:
-
   Vector();
   Vector(const double x1, const double x2, const double x3);
@@ -42 +43 @@
   double DistanceSquared(const Vector &y) const;
   double DistanceToSpace(const Space& space) const;
-  double PeriodicDistance(const Vector &y, const double * const cell_size) const;
-  double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
   double ScalarProduct(const Vector &y) const;
   double Angle(const Vector &y) const;
@@ -58 +57 @@
   Vector Projection(const Vector &y) const;
   void ScaleAll(const double *factor);
+  void ScaleAll(const Vector &factor);
   void Scale(const double factor);
-  void MatrixMultiplication(const double * const M);
-  bool InverseMatrixMultiplication(const double * const M);
-  void KeepPeriodic(const double * const matrix);
   bool GetOneNormalVector(const Vector &x1);
   bool MakeNormalTo(const Vector &y1);
-  bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
-  void WrapPeriodically(const double * const M, const double * const Minv);
   std::pair<Vector,Vector> partition(const Vector&) const;
   std::pair<pointset,Vector> partition(const pointset&) const;
@@ -79 +74 @@
 
   // Access to internal structure
-  gsl_vector* get();
+  VectorContent* get();
 
   // Methods that are derived directly from other methods
@@ -104 +99 @@
 
 private:
-  gsl_vector *content;
+  Vector(VectorContent *);
+  VectorContent *content;
 
 };

src/vector_ops.cpp

@@ -23 +23 @@
 #include <gsl/gsl_permutation.h>
 #include <gsl/gsl_vector.h>
+#include <gsl/gsl_multimin.h>
 
 /**

src/verbose.cpp

@@ -5 +5 @@
 #include "info.hpp"
 #include "verbose.hpp"
+#include <iostream>
 
 /** Prints the tabs according to verbosity stored in the temporary constructed class.

src/verbose.hpp

@@ -18 +18 @@
 #endif
 
-#include <iostream>
+#include <iosfwd>
 
 /************************************* Class Verbose & Binary *******************************/

tests/regression/Domain/4/post/test.conf

@@ -47 +47 @@
 BoxLength                       # (Length of a unit cell)
 1       
-0       0       
+0       1       
 0       0       2