Changes in / [4e10f5:192f6e]


Ignore:
Files:
4 added
22 deleted
108 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/ActionRegistry.hpp

    r4e10f5 r192f6e  
    99#define ACTIONREGISTRY_HPP_
    1010
    11 #include <iosfwd>
     11#include <iostream>
    1212#include <string>
    1313#include <map>

  • src/Actions/AnalysisAction/PairCorrelationAction.cpp

    r4e10f5 r192f6e  
    1212#include "boundary.hpp"
    1313#include "linkedcell.hpp"
    14 #include "verbose.hpp"
    1514#include "log.hpp"
    1615#include "element.hpp"
    … …  
    6867  dialog = UIFactory::getInstance().makeDialog();
    6968  if (type == "P")
    70     dialog->queryVector("position", &Point, false, MapOfActions::getInstance().getDescription("position"));
     69    dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
    7170  if (type == "S")
    7271    dialog->queryMolecule("molecule-by-id", &Boundary, MapOfActions::getInstance().getDescription("molecule-by-id"));

  • src/Actions/AtomAction/AddAction.cpp

    r4e10f5 r192f6e  
    4141
    4242  dialog->queryElement(NAME, &elements, MapOfActions::getInstance().getDescription(NAME));
    43   dialog->queryVector("position", &position, true, MapOfActions::getInstance().getDescription("position"));
     43  dialog->queryVector("position", &position, World::getInstance().getDomain(), true, MapOfActions::getInstance().getDescription("position"));
    4444  cout << "pre-dialog" << endl;
    4545

  • src/Actions/AtomsCalculation_impl.hpp

    r4e10f5 r192f6e  
    1111#include "Actions/AtomsCalculation.hpp"
    1212#include "Actions/Calculation_impl.hpp"
     13
     14#include <iostream>
    1315
    1416using namespace std;

  • src/Actions/MapOfActions.cpp

    r4e10f5 r192f6e  
    1717#include <boost/optional.hpp>
    1818#include <boost/program_options.hpp>
    19 
    20 #include <iostream>
    2119
    2220#include "CommandLineParser.hpp"

  • src/Actions/MoleculeAction/BondFileAction.cpp

    r4e10f5 r192f6e  
    2424#include "config.hpp"
    2525#include "defs.hpp"
    26 #include "verbose.hpp"
    2726#include "log.hpp"
    2827#include "molecule.hpp"

  • src/Actions/MoleculeAction/FillWithMoleculeAction.cpp

    r4e10f5 r192f6e  
    6262
    6363  dialog->queryString(NAME, &filename, MapOfActions::getInstance().getDescription(NAME));
    64   dialog->queryVector("distances", &distances, false, MapOfActions::getInstance().getDescription("distances"));
    65   dialog->queryVector("lengths", &lengths, false, MapOfActions::getInstance().getDescription("lengths"));
     64  dialog->queryVector("distances", &distances, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("distances"));
     65  dialog->queryVector("lengths", &lengths, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("lengths"));
    6666  dialog->queryBoolean("DoRotate", &DoRotate, MapOfActions::getInstance().getDescription("DoRotate"));
    6767  dialog->queryDouble("MaxDistance", &MaxDistance, MapOfActions::getInstance().getDescription("MaxDistance"));

  • src/Actions/MoleculeAction/SuspendInWaterAction.cpp

    r4e10f5 r192f6e  
    2222#include "boundary.hpp"
    2323#include "config.hpp"
    24 #include "verbose.hpp"
    2524#include "log.hpp"
    2625#include "config.hpp"

  • src/Actions/MoleculeAction/TranslateAction.cpp

    r4e10f5 r192f6e  
    5555  bool periodic = false;
    5656
    57   dialog->queryVector(NAME, &x, false, MapOfActions::getInstance().getDescription(NAME));
     57  dialog->queryVector(NAME, &x, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    5858  dialog->queryMolecule("molecule-by-id", &mol, MapOfActions::getInstance().getDescription("molecule-by-id"));
    5959  dialog->queryBoolean("periodic", &periodic, MapOfActions::getInstance().getDescription("periodic"));

  • src/Actions/WorldAction/AddEmptyBoundaryAction.cpp

    r4e10f5 r192f6e  
    4141  int j=0;
    4242
    43   dialog->queryVector(NAME, &boundary, false, MapOfActions::getInstance().getDescription(NAME));
     43  dialog->queryVector(NAME, &boundary, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    4444
    4545  if(dialog->display()) {
    … …  
    5959    }
    6060    // set new box size
    61     double * const cell_size = new double[6];
     61    double * const cell_size = World::getInstance().getDomain();
    6262    for (j=0;j<6;j++)
    6363      cell_size[j] = 0.;
    … …  
    6767      cell_size[j] = (Max[i]-Min[i]+2.*boundary[i]);
    6868    }
    69     World::getInstance().setDomain(cell_size);
    70     delete[] cell_size;
    7169    // translate all atoms, such that Min is aty (0,0,0)
    7270    AtomRunner = AllAtoms.begin();

  • src/Actions/WorldAction/CenterInBoxAction.cpp

    r4e10f5 r192f6e  
    3434  Dialog *dialog = UIFactory::getInstance().makeDialog();
    3535
    36   Box& cell_size = World::getInstance().getDomain();
     36  double * cell_size = World::getInstance().getDomain();
    3737  dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
    3838
    3939  if(dialog->display()) {
     40    World::getInstance().setDomain(cell_size);
    4041    // center
    4142    vector<molecule *> AllMolecules = World::getInstance().getAllMolecules();

  • src/Actions/WorldAction/CenterOnEdgeAction.cpp

    r4e10f5 r192f6e  
    1313#include "vector.hpp"
    1414#include "World.hpp"
    15 #include "Matrix.hpp"
    1615
    1716#include <iostream>
    … …  
    3837  Vector Min;
    3938  Vector Max;
     39  int j=0;
    4040
    4141  dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
    … …  
    5757    }
    5858    // set new box size
    59     Matrix domain;
     59    double * const cell_size = World::getInstance().getDomain();
     60    for (j=0;j<6;j++)
     61      cell_size[j] = 0.;
     62    j=-1;
    6063    for (int i=0;i<NDIM;i++) {
    61       double tmp = Max[i]-Min[i];
    62       tmp = fabs(tmp)>=1. ? tmp : 1.0;
    63       domain.at(i,i) = tmp;
     64      j += i+1;
     65      cell_size[j] = (Max[i]-Min[i]);
    6466    }
    65     World::getInstance().setDomain(domain);
     67    World::getInstance().setDomain(cell_size);
    6668    // translate all atoms, such that Min is aty (0,0,0)
    6769    for (vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner)

  • src/Actions/WorldAction/ChangeBoxAction.cpp

    r4e10f5 r192f6e  
    1212#include "verbose.hpp"
    1313#include "World.hpp"
    14 #include "Box.hpp"
    15 #include "Matrix.hpp"
    1614
    1715#include <iostream>
    … …  
    3634  Dialog *dialog = UIFactory::getInstance().makeDialog();
    3735
    38   Box& cell_size = World::getInstance().getDomain();
     36  double * cell_size = World::getInstance().getDomain();
    3937  dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
    4038
    4139  if(dialog->display()) {
    42     DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size.getM() << endl);
     40    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);
     41    World::getInstance().setDomain(cell_size);
    4342    delete dialog;
    4443    return Action::success;

  • src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp

    r4e10f5 r192f6e  
    4141
    4242  dialog->queryDouble(NAME, &radius, MapOfActions::getInstance().getDescription(NAME));
    43   dialog->queryVector("position", &point, false, MapOfActions::getInstance().getDescription("position"));
     43  dialog->queryVector("position", &point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
    4444
    4545  if(dialog->display()) {

  • src/Actions/WorldAction/RepeatBoxAction.cpp

    r4e10f5 r192f6e  
    1313#include "molecule.hpp"
    1414#include "vector.hpp"
    15 #include "Matrix.hpp"
    1615#include "verbose.hpp"
    1716#include "World.hpp"
    18 #include "Box.hpp"
    1917
    2018#include <iostream>
    … …  
    4846  MoleculeListClass *molecules = World::getInstance().getMolecules();
    4947
    50   dialog->queryVector(NAME, &Repeater, false, MapOfActions::getInstance().getDescription(NAME));
     48  dialog->queryVector(NAME, &Repeater, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    5149  //dialog->queryMolecule("molecule-by-id", &mol,MapOfActions::getInstance().getDescription("molecule-by-id"));
    5250  vector<molecule *> AllMolecules;
    … …  
    6159  if(dialog->display()) {
    6260    (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
    63     Matrix M = World::getInstance().getDomain().getM();
    64     Matrix newM = M;
     61    double * const cell_size = World::getInstance().getDomain();
     62    double *M = ReturnFullMatrixforSymmetric(cell_size);
    6563    Vector x,y;
    6664    int n[NDIM];
    67     Matrix repMat;
    6865    for (int axis = 0; axis < NDIM; axis++) {
    6966      Repeater[axis] = floor(Repeater[axis]);
    … …  
    7269        Repeater[axis] = 1;
    7370      }
    74       repMat.at(axis,axis) = Repeater[axis];
     71      cell_size[ ((axis == 1) ? 2 : (axis == 2) ? 5 : 0) ] *= Repeater[axis];
    7572    }
    76     newM *= repMat;
    77     World::getInstance().setDomain(newM);
    7873
     74    World::getInstance().setDomain(cell_size);
    7975    molecule *newmol = NULL;
    8076    Vector ** vectors = NULL;
    … …  
    10399                DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
    104100              x = y;
    105               x *= M;
     101              x.MatrixMultiplication(M);
    106102              newmol = World::getInstance().createMolecule();
    107103              molecules->insert(newmol);
    … …  
    123119      }
    124120    }
     121    delete(M);
    125122    delete dialog;
    126123    return Action::success;

  • src/Actions/WorldAction/ScaleBoxAction.cpp

    r4e10f5 r192f6e  
    1414#include "verbose.hpp"
    1515#include "World.hpp"
    16 #include "Box.hpp"
    17 #include "Matrix.hpp"
    1816
    1917#include <iostream>
    … …  
    3937  Vector Scaler;
    4038  double x[NDIM];
     39  int j=0;
    4140
    42   dialog->queryVector(NAME, &Scaler, false, MapOfActions::getInstance().getDescription(NAME));
     41  dialog->queryVector(NAME, &Scaler, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    4342
    4443  if(dialog->display()) {
    … …  
    5049      (*AtomRunner)->x.ScaleAll(x);
    5150    }
    52 
    53     Matrix M = World::getInstance().getDomain().getM();
    54     Matrix scale;
    55 
     51    j = -1;
     52    double * const cell_size = World::getInstance().getDomain();
    5653    for (int i=0;i<NDIM;i++) {
    57       scale.at(i,i) = x[i];
     54      j += i+1;
     55      cell_size[j]*=x[i];
    5856    }
    59     M *= scale;
    60     World::getInstance().setDomain(M);
     57    World::getInstance().setDomain(cell_size);
    6158
    6259    delete dialog;

  • src/ConfigFileBuffer.cpp

    r4e10f5 r192f6e  
    99#include "helpers.hpp"
    1010#include "lists.hpp"
    11 #include "verbose.hpp"
    1211#include "log.hpp"
    1312#include "World.hpp"

  • src/ConfigFileBuffer.hpp

    r4e10f5 r192f6e  
    1717#endif
    1818
    19 #include <iosfwd>
     19#include <iostream>
    2020
    2121/****************************************** forward declarations *****************************/

  • src/Exceptions/CustomException.cpp

    r4e10f5 r192f6e  
    99
    1010#include "CustomException.hpp"
    11 #include <iostream>
    1211
    1312using namespace std;

  • src/Exceptions/CustomException.hpp

    r4e10f5 r192f6e  
    1010
    1111#include <exception>
    12 #include <iosfwd>
    13 #include <string>
     12#include <iostream>
    1413
    1514/**

  • src/Helpers/Assert.hpp

    r4e10f5 r192f6e  
    1111#include<sstream>
    1212#include<string>
    13 #include<iosfwd>
     13#include<iostream>
    1414#include<vector>
    1515#include<map>

  • src/Helpers/MemDebug.cpp

    r4e10f5 r192f6e  
    66 */
    77
    8 // NDEBUG implies NO_MEMDEBUG
    9 #ifdef NDEBUG
    10 # ifndef NO_MEMDEBUG
    11 #   define NO_MEMDEBUG
    12 # endif
    13 #endif
    14 
    15 // NO_MEMDEBUG and MEMDEBUG are mutually exclusive, but at least one must be set
    16 #ifdef NO_MEMDEBUG
    17 # ifdef MEMDEBUG
    18 #   undef MEMDEBUG
    19 # endif
    20 #else
    21 # ifndef MEMDEBUG
    22 #   define MEMDEBUG
    23 # endif
    24 #endif
    25 
    26 #ifdef MEMDEBUG
     8#ifndef NDEBUG
     9#ifndef NO_MEMDEBUG
    2710
    2811#include <iostream>
    … …  
    506489}
    507490#endif
     491#endif

  • src/Helpers/MemDebug.hpp

    r4e10f5 r192f6e  
    2121 * your sourcefiles.
    2222 */
     23#ifndef NDEBUG
     24#ifndef NO_MEMDEBUG
    2325
    24 // Set all flags in a way that makes sense
    25 
    26 // NDEBUG implies NO_MEMDEBUG
    27 #ifdef NDEBUG
    28 # ifndef NO_MEMDEBUG
    29 #   define NO_MEMDEBUG
    30 # endif
     26#ifndef MEMDEBUG
     27#define MEMDEBUG
    3128#endif
    32 
    33 // NO_MEMDEBUG and MEMDEBUG are mutually exclusive, but at least one must be set
    34 #ifdef NO_MEMDEBUG
    35 # ifdef MEMDEBUG
    36 #   undef MEMDEBUG
    37 # endif
    38 #else
    39 # ifndef MEMDEBUG
    40 #   define MEMDEBUG
    41 # endif
    42 #endif
    43 
    44 #ifdef MEMDEBUG
    4529
    4630#include <new>
    … …  
    9983#endif
    10084
    101 #else
     85#endif
     86#endif
     87
     88
     89#ifdef NDEBUG
     90#undef MEMDEBUG
     91#endif
     92
     93#ifndef MEMDEBUG
    10294// memory debugging was disabled
    10395
    … …  
    112104
    113105#endif
    114 
    115 
    116106#endif /* MEMDEBUG_HPP_ */

  • src/Line.cpp

    r4e10f5 r192f6e  
    215215}
    216216
    217 std::vector<Vector> Line::getSphereIntersections() const{
    218   std::vector<Vector> res;
    219 
    220   // line is kept in normalized form, so we can skip a lot of calculations
    221   double discriminant = 1-origin->NormSquared();
    222   // we might have 2, 1 or 0 solutions, depending on discriminant
    223   if(discriminant>=0){
    224     if(discriminant==0){
    225       res.push_back(*origin);
    226     }
    227     else{
    228       Vector helper = sqrt(discriminant)*(*direction);
    229       res.push_back(*origin+helper);
    230       res.push_back(*origin-helper);
    231     }
    232   }
    233   return res;
    234 }
    235 
    236217Line makeLineThrough(const Vector &x1, const Vector &x2){
    237218  if(x1==x2){

  • src/Line.hpp

    r4e10f5 r192f6e  
    3838  Plane getOrthogonalPlane(const Vector &origin) const;
    3939
    40   std::vector<Vector> getSphereIntersections() const;
    41 
    4240private:
    4341  std::auto_ptr<Vector> origin;

  • src/Makefile.am

    r4e10f5 r192f6e  
    3838  vector.cpp
    3939                           
    40 LINALGHEADER = \
    41   gslmatrix.hpp \
     40LINALGHEADER = gslmatrix.hpp \
    4241  gslvector.hpp \
    4342  linearsystemofequations.hpp \
    … …  
    7675  Actions/MethodAction.hpp \
    7776  Actions/Process.hpp
    78 
    79 EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
    80                                   Exceptions/LinearDependenceException.cpp \
    81                                   Exceptions/MathException.cpp \
    82                                   Exceptions/NotInvertibleException.cpp \
    83                                   Exceptions/SkewException.cpp \
    84                                   Exceptions/ZeroVectorException.cpp
    85                                  
    86 EXCEPTIONHEADER = Exceptions/CustomException.hpp \
    87                                   Exceptions/LinearDependenceException.hpp \
    88                                   Exceptions/MathException.hpp \
    89                                   Exceptions/NotInvertibleException.hpp \
    90                                   Exceptions/SkewException.hpp \
    91                                   Exceptions/ZeroVectorException.hpp
     77 
    9278
    9379PARSERSOURCE = \
    … …  
    115101  Patterns/Observer.hpp \
    116102  Patterns/Singleton.hpp
    117  
    118 SHAPESOURCE = \
    119   Shapes/BaseShapes.cpp \
    120   Shapes/Shape.cpp \
    121   Shapes/ShapeOps.cpp
    122 SHAPEHEADER = \
    123   Shapes/BaseShapes.hpp \
    124   Shapes/Shape.hpp \
    125   Shapes/ShapeOps.hpp
    126  
    127103
    128104QTUIMOC_HEADER = UIElements/QT4/QTDialog.hpp \
    … …  
    158134  Descriptors/MoleculeNameDescriptor.hpp \
    159135  Descriptors/MoleculePtrDescriptor.hpp
     136                                   
     137EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
     138                                  Exceptions/LinearDependenceException.cpp \
     139                                  Exceptions/MathException.cpp \
     140                                  Exceptions/SkewException.cpp \
     141                                  Exceptions/ZeroVectorException.cpp
     142                                 
     143EXCEPTIONHEADER = Exceptions/CustomException.hpp \
     144                                  Exceptions/LinearDependenceException.hpp \
     145                                  Exceptions/MathException.hpp \
     146                                  Exceptions/SkewException.hpp \
     147                                  Exceptions/ZeroVectorException.hpp
    160148                                 
    161149QTUISOURCE = ${QTUIMOC_TARGETS} \
    … …  
    177165  ${ACTIONSSOURCE} \
    178166  ${ATOMSOURCE} \
    179   ${EXCEPTIONSOURCE} \
    180167  ${PATTERNSOURCE} \
    181168  ${PARSERSOURCE} \
    182   ${SHAPESOURCE} \
    183169  ${DESCRIPTORSOURCE} \
    184170  ${HELPERSOURCE} \
     171  ${EXCEPTIONSOURCE} \
    185172  bond.cpp \
    186173  bondgraph.cpp \
    187174  boundary.cpp \
    188   Box.cpp \
    189175  CommandLineParser.cpp \
    190176  config.cpp \
    … …  
    202188  log.cpp \
    203189  logger.cpp \
    204   Matrix.cpp \
    205190  moleculelist.cpp \
    206191  molecule.cpp \
    … …  
    227212  ${ACTIONSHEADER} \
    228213  ${ATOMHEADER} \
    229   ${EXCEPTIONHEADER} \
    230214  ${PARSERHEADER} \
    231215  ${PATTERNHEADER} \
    232   ${SHAPEHEADER} \
    233216  ${DESCRIPTORHEADER} \
     217  ${EXCEPTIONHEADER} \
    234218  bond.hpp \
    235219  bondgraph.hpp \
    236220  boundary.hpp \
    237   Box.hpp \
    238221  CommandLineParser.hpp \
    239222  config.hpp \
    … …  
    253236  log.hpp \
    254237  logger.hpp \
    255   Matrix.hpp \
    256238  molecule.hpp \
    257239  molecule_template.hpp \

  • src/Parser/MpqcParser.hpp

    r4e10f5 r192f6e  
    1111#include "FormatParser.hpp"
    1212
    13 #include <iosfwd>
     13#include <iostream>
    1414
    1515/**

  • src/Parser/PcpParser.cpp

    r4e10f5 r192f6e  
    77
    88#include <iostream>
    9 #include <iomanip>
    109
    1110#include "atom.hpp"
    … …  
    2120#include "verbose.hpp"
    2221#include "World.hpp"
    23 #include "Matrix.hpp"
    24 #include "Box.hpp"
    2522
    2623/** Constructor of PcpParser.
    … …  
    213210  // Unit cell and magnetic field
    214211  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    215   double *cell_size = new double[6];
     212  double * const cell_size = World::getInstance().getDomain();
    216213  cell_size[0] = BoxLength[0];
    217214  cell_size[1] = BoxLength[3];
    … …  
    220217  cell_size[4] = BoxLength[7];
    221218  cell_size[5] = BoxLength[8];
    222   World::getInstance().setDomain(cell_size);
    223   delete[] cell_size;
    224219  //if (1) fprintf(stderr,"\n");
    225220
    … …  
    332327void PcpParser::save(std::ostream* file)
    333328{
    334   const Matrix &domain = World::getInstance().getDomain().getM();
     329  const double * const cell_size = World::getInstance().getDomain();
    335330  class ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
    336331  if (!file->fail()) {
    … …  
    417412    *file << endl;
    418413    *file << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    419     *file << domain.at(0,0) << "\t" << endl;
    420     *file << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
    421     *file << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     414    *file << cell_size[0] << "\t" << endl;
     415    *file << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     416    *file << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
    422417    // FIXME
    423418    *file << endl;

  • src/Parser/PcpParser.hpp

    r4e10f5 r192f6e  
    99#define PCPPARSER_HPP_
    1010
    11 #include <iosfwd>
     11#include <iostream>
    1212#include "Parser/FormatParser.hpp"
    1313

  • src/Patterns/Cacheable.hpp

    r4e10f5 r192f6e  
    1212#include <boost/function.hpp>
    1313#include <boost/shared_ptr.hpp>
     14#include <iostream>
    1415
    1516#include "Helpers/Assert.hpp"

  • src/Plane.hpp

    r4e10f5 r192f6e  
    1111#include <memory>
    1212#include <vector>
    13 #include <iosfwd>
     13#include <iostream>
    1414#include "Space.hpp"
    1515#include "Exceptions/LinearDependenceException.hpp"

  • src/UIElements/CommandLineUI/CommandLineDialog.cpp

    r4e10f5 r192f6e  
    2727#include "verbose.hpp"
    2828#include "World.hpp"
    29 #include "Box.hpp"
    3029
    3130#include "atom.hpp"
    … …  
    7877}
    7978
    80 void CommandLineDialog::queryVector(const char* title, Vector *target, bool check, string _description) {
    81   registerQuery(new VectorCommandLineQuery(title,target,check, _description));
    82 }
    83 
    84 void CommandLineDialog::queryBox(const char* title, Box* cellSize, string _description) {
     79void CommandLineDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string _description) {
     80  registerQuery(new VectorCommandLineQuery(title,target,cellSize,check, _description));
     81}
     82
     83void CommandLineDialog::queryBox(const char* title, double ** const cellSize, string _description) {
    8584  registerQuery(new BoxCommandLineQuery(title,cellSize,_description));
    8685}
    … …  
    223222}
    224223
    225 CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, bool _check, string _description) :
    226     Dialog::VectorQuery(title,_target,_check, _description)
     224CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, const double *const _cellSize, bool _check, string _description) :
     225    Dialog::VectorQuery(title,_target,_cellSize,_check, _description)
    227226{}
    228227
    … …  
    245244
    246245
    247 CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, Box* _cellSize, string _description) :
     246CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, double ** const _cellSize, string _description) :
    248247    Dialog::BoxQuery(title,_cellSize, _description)
    249248{}

  • src/UIElements/CommandLineUI/CommandLineDialog.hpp

    r4e10f5 r192f6e  
    3535  virtual void queryAtom(const char*,atom**, std::string = "");
    3636  virtual void queryMolecule(const char*,molecule**,std::string = "");
    37   virtual void queryVector(const char*,Vector *,bool, std::string = "");
    38   virtual void queryBox(const char*,Box *, std::string = "");
     37  virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
     38  virtual void queryBox(const char*,double ** const, std::string = "");
    3939  virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
    4040
    … …  
    9999  class VectorCommandLineQuery : public Dialog::VectorQuery {
    100100  public:
    101     VectorCommandLineQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
     101    VectorCommandLineQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
    102102    virtual ~VectorCommandLineQuery();
    103103    virtual bool handle();
    … …  
    106106  class BoxCommandLineQuery : public Dialog::BoxQuery {
    107107  public:
    108     BoxCommandLineQuery(std::string title,Box* _cellSize, std::string _description = "");
     108    BoxCommandLineQuery(std::string title,double ** const _cellSize, std::string _description = "");
    109109    virtual ~BoxCommandLineQuery();
    110110    virtual bool handle();

  • src/UIElements/Dialog.cpp

    r4e10f5 r192f6e  
    1010#include "Dialog.hpp"
    1111
    12 #include "verbose.hpp"
    1312#include "atom.hpp"
    1413#include "element.hpp"
    1514#include "molecule.hpp"
    1615#include "vector.hpp"
    17 #include "Matrix.hpp"
    18 #include "Box.hpp"
    1916
    2017using namespace std;
    … …  
    188185// Vector Queries
    189186
    190 Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,bool _check, std::string _description) :
     187Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description) :
    191188  Query(title, _description),
     189  cellSize(_cellSize),
    192190  check(_check),
    193191  target(_target)
    … …  
    207205// Box Queries
    208206
    209 Dialog::BoxQuery::BoxQuery(std::string title, Box* _cellSize, std::string _description) :
     207Dialog::BoxQuery::BoxQuery(std::string title, double ** const _cellSize, std::string _description) :
    210208  Query(title, _description),
    211209  target(_cellSize)
    … …  
    220218
    221219void Dialog::BoxQuery::setResult() {
    222   (*target)= ReturnFullMatrixforSymmetric(tmp);
     220  for (int i=0;i<6;i++) {
     221    (*target)[i] = tmp[i];
     222  }
    223223}
    224224

  • src/UIElements/Dialog.hpp

    r4e10f5 r192f6e  
    1414
    1515class atom;
    16 class Box;
    1716class element;
    1817class molecule;
    … …  
    4443  virtual void queryAtom(const char*,atom**,std::string = "")=0;
    4544  virtual void queryMolecule(const char*,molecule**, std::string = "")=0;
    46   virtual void queryVector(const char*,Vector *,bool, std::string = "")=0;
    47   virtual void queryBox(const char*,Box*, std::string = "")=0;
     45  virtual void queryVector(const char*,Vector *,const double *const,bool, std::string = "")=0;
     46  virtual void queryBox(const char*,double ** const, std::string = "")=0;
    4847  virtual void queryElement(const char*, std::vector<element *> *, std::string = "")=0;
    4948
    … …  
    177176  class VectorQuery : public Query {
    178177  public:
    179       VectorQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
     178      VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
    180179      virtual ~VectorQuery();
    181180      virtual bool handle()=0;
    … …  
    183182    protected:
    184183      Vector *tmp;
     184      const double *const cellSize;
    185185      bool check;
    186186    private:
    … …  
    190190  class BoxQuery : public Query {
    191191  public:
    192       BoxQuery(std::string title,Box *_cellSize, std::string _description = "");
     192      BoxQuery(std::string title,double ** const _cellSize, std::string _description = "");
    193193      virtual ~BoxQuery();
    194194      virtual bool handle()=0;
    195195      virtual void setResult();
    196196    protected:
    197       double* tmp;
     197      double *tmp;
    198198    private:
    199       Box* target;
     199      double **target;
    200200  };
    201201

  • src/UIElements/QT4/QTDialog.cpp

    r4e10f5 r192f6e  
    2929#include "molecule.hpp"
    3030#include "Descriptors/MoleculeIdDescriptor.hpp"
    31 #include "Matrix.hpp"
    32 #include "Box.hpp"
    3331
    3432
    … …  
    9593}
    9694
    97 void QTDialog::queryBox(char const*, Box*, string){
     95void QTDialog::queryBox(char const*, double**, string){
    9896  // TODO
    9997  ASSERT(false, "Not implemented yet");
    … …  
    125123}
    126124
    127 void QTDialog::queryVector(const char* title, Vector *target, bool check,string) {
    128   registerQuery(new VectorQTQuery(title,target,check,inputLayout,this));
     125void QTDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check,string) {
     126  registerQuery(new VectorQTQuery(title,target,cellSize,check,inputLayout,this));
    129127}
    130128
    … …  
    293291}
    294292
    295 QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
    296     Dialog::VectorQuery(title,_target,_check),
    297     parent(_parent)
    298 {
    299   const Matrix& M = World::getInstance().getDomain().getM();
     293QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
     294    Dialog::VectorQuery(title,_target,_cellSize,_check),
     295    parent(_parent)
     296{
     297  // 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
     298  int j = -1;
    300299  const char *coords[3] = {"x:","y:","z:"};
    301300  mainLayout= new QHBoxLayout();
    … …  
    305304  mainLayout->addLayout(subLayout);
    306305  for(int i=0; i<3; i++) {
     306    j+=i+1;
    307307    coordLayout[i] = new QHBoxLayout();
    308308    subLayout->addLayout(coordLayout[i]);
    … …  
    310310    coordLayout[i]->addWidget(coordLabel[i]);
    311311    coordInput[i] = new QDoubleSpinBox();
    312     coordInput[i]->setRange(0,M.at(i,i));
     312    coordInput[i]->setRange(0,cellSize[j]);
    313313    coordInput[i]->setDecimals(3);
    314314    coordLayout[i]->addWidget(coordInput[i]);

  • src/UIElements/QT4/QTDialog.hpp

    r4e10f5 r192f6e  
    4343  virtual void queryAtom(const char*,atom**,std::string = "");
    4444  virtual void queryMolecule(const char*,molecule**,std::string = "");
    45   virtual void queryVector(const char*,Vector *,bool,std::string = "");
    46   virtual void queryBox(const char*,Box*, std::string = "");
     45  virtual void queryVector(const char*,Vector *,const double *const,bool,std::string = "");
     46  virtual void queryBox(const char*,double ** const, std::string = "");
    4747  virtual void queryElement(const char*,std::vector<element *> *_target,std::string = "");
    4848
    … …  
    124124    class VectorQTQuery : public Dialog::VectorQuery {
    125125    public:
    126       VectorQTQuery(std::string title,Vector *_target,bool _check,QBoxLayout *,QTDialog *);
     126      VectorQTQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check,QBoxLayout *,QTDialog *);
    127127      virtual ~VectorQTQuery();
    128128      virtual bool handle();

  • src/UIElements/QT4/QTMainWindow.cpp

    r4e10f5 r192f6e  
    2121#include "atom.hpp"
    2222#include "molecule.hpp"
    23 #include "verbose.hpp"
    2423#include "Actions/Action.hpp"
    2524#include "Actions/ActionRegistry.hpp"

  • src/UIElements/TextUI/TextDialog.cpp

    r4e10f5 r192f6e  
    2525#include "molecule.hpp"
    2626#include "vector.hpp"
    27 #include "Matrix.hpp"
    28 #include "Box.hpp"
    2927
    3028using namespace std;
    … …  
    7270}
    7371
    74 void TextDialog::queryVector(const char* title, Vector *target, bool check, string description) {
    75   registerQuery(new VectorTextQuery(title,target,check,description));
    76 }
    77 
    78 void TextDialog::queryBox(const char* title,Box* cellSize, string description) {
     72void TextDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string description) {
     73  registerQuery(new VectorTextQuery(title,target,cellSize,check,description));
     74}
     75
     76void TextDialog::queryBox(const char* title,double ** const cellSize, string description) {
    7977  registerQuery(new BoxTextQuery(title,cellSize,description));
    8078}
    … …  
    272270}
    273271
    274 TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, bool _check, std::string _description) :
    275     Dialog::VectorQuery(title,_target,_check,_description)
     272TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check, std::string _description) :
     273    Dialog::VectorQuery(title,_target,_cellSize,_check,_description)
    276274{}
    277275
    … …  
    282280  Log() << Verbose(0) << getTitle();
    283281
    284   const Matrix &M = World::getInstance().getDomain().getM();
    285282  char coords[3] = {'x','y','z'};
     283  int j = -1;
    286284  for (int i=0;i<3;i++) {
     285    j += i+1;
    287286    do {
    288       Log() << Verbose(0) << coords[i] << "[0.." << M.at(i,i) << "]: ";
     287      Log() << Verbose(0) << coords[i] << "[0.." << cellSize[j] << "]: ";
    289288      cin >> (*tmp)[i];
    290     } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= M.at(i,i))) && (check));
     289    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= cellSize[j])) && (check));
    291290  }
    292291  return true;
    293292}
    294293
    295 TextDialog::BoxTextQuery::BoxTextQuery(std::string title, Box* _cellSize, std::string _description) :
     294TextDialog::BoxTextQuery::BoxTextQuery(std::string title, double ** const _cellSize, std::string _description) :
    296295    Dialog::BoxQuery(title,_cellSize,_description)
    297296{}

  • src/UIElements/TextUI/TextDialog.hpp

    r4e10f5 r192f6e  
    3232  virtual void queryAtom(const char*,atom**,std::string = "");
    3333  virtual void queryMolecule(const char*,molecule**,std::string = "");
    34   virtual void queryVector(const char*,Vector *,bool, std::string = "");
    35   virtual void queryBox(const char*,Box*, std::string = "");
     34  virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
     35  virtual void queryBox(const char*,double ** const, std::string = "");
    3636  virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
    3737
    … …  
    9696  class VectorTextQuery : public Dialog::VectorQuery {
    9797  public:
    98     VectorTextQuery(std::string title,Vector *_target,bool _check, std::string _description = NULL);
     98    VectorTextQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = NULL);
    9999    virtual ~VectorTextQuery();
    100100    virtual bool handle();
    … …  
    103103  class BoxTextQuery : public Dialog::BoxQuery {
    104104  public:
    105     BoxTextQuery(std::string title,Box* _cellSize, std::string _description = NULL);
     105    BoxTextQuery(std::string title,double ** const _cellSize, std::string _description = NULL);
    106106    virtual ~BoxTextQuery();
    107107    virtual bool handle();

  • src/UIElements/TextUI/TextWindow.hpp

    r4e10f5 r192f6e  
    1111#include "MainWindow.hpp"
    1212
    13 #include <string>
    1413#include <set>
    1514

  • src/UIElements/Views/StreamStringView.hpp

    r4e10f5 r192f6e  
    1010
    1111#include <boost/function.hpp>
    12 #include <iosfwd>
     12#include <iostream>
    1313
    1414#include "Views/StringView.hpp"

  • src/VectorSet.hpp

    r4e10f5 r192f6e  
    1212#include <functional>
    1313#include <algorithm>
    14 #include <limits>
    1514
    1615/**
    … …  
    2019 */
    2120
    22 #include "vector.hpp"
    23 #include <list>
     21class Vector;
    2422
    2523// this tests, whether we actually have a Vector
    … …  
    5048  void translate(const Vector &translater){
    5149    // this is needed to allow template lookup
    52     transform(this->begin(),this->end(),this->begin(),std::bind1st(std::plus<Vector>(),translater));
    53   }
    54 
    55   double minDistSquared(const Vector &point){
    56     if(!this->size())
    57       return std::numeric_limits<double>::infinity();
    58     std::list<double> helper;
    59     helper.resize(this->size());
    60     transform(this->begin(),this->end(),
    61               helper.begin(),
    62               std::bind2nd(std::mem_fun_ref(&Vector::DistanceSquared),point));
    63     return *min_element(helper.begin(),helper.end());
     50    transform(this->begin(),this->end(),this->begin(),bind1st(plus<Vector>(),translater));
    6451  }
    6552};
    6653
    67 // allows simpler definition of VectorSets
    68 #define VECTORSET(container_type) VectorSet<container_type<Vector> >
    69 
    7054#endif /* VECTORSET_HPP_ */

  • src/World.cpp

    r4e10f5 r192f6e  
    2222#include "Actions/ManipulateAtomsProcess.hpp"
    2323#include "Helpers/Assert.hpp"
    24 #include "Box.hpp"
    25 #include "Matrix.hpp"
    2624
    2725#include "Patterns/Singleton_impl.hpp"
    … …  
    7674// system
    7775
    78 Box& World::getDomain() {
    79   return *cell_size;
    80 }
    81 
    82 void World::setDomain(const Matrix &mat){
    83   *cell_size = mat;
     76double * World::getDomain() {
     77  return cell_size;
    8478}
    8579
    8680void World::setDomain(double * matrix)
    8781{
    88   Matrix M = ReturnFullMatrixforSymmetric(matrix);
    89   cell_size->setM(M);
     82  OBSERVE;
     83  cell_size[0] = matrix[0];
     84  cell_size[1] = matrix[1];
     85  cell_size[2] = matrix[2];
     86  cell_size[3] = matrix[3];
     87  cell_size[4] = matrix[4];
     88  cell_size[5] = matrix[5];
    9089}
    9190
    … …  
    140139  molecules.erase(id);
    141140}
     141
     142double *World::cell_size = NULL;
    142143
    143144atom *World::createAtom(){
    … …  
    302303    molecules_deprecated(new MoleculeListClass(this))
    303304{
    304   cell_size = new Box;
    305   Matrix domain;
    306   domain.at(0,0) = 20;
    307   domain.at(1,1) = 20;
    308   domain.at(2,2) = 20;
    309   cell_size->setM(domain);
     305  cell_size = new double[6];
     306  cell_size[0] = 20.;
     307  cell_size[1] = 0.;
     308  cell_size[2] = 20.;
     309  cell_size[3] = 0.;
     310  cell_size[4] = 0.;
     311  cell_size[5] = 20.;
    310312  defaultName = "none";
    311313  molecules_deprecated->signOn(this);
    … …  
    315317{
    316318  molecules_deprecated->signOff(this);
    317   delete cell_size;
     319  delete[] cell_size;
    318320  delete molecules_deprecated;
    319321  delete periode;

  • src/World.hpp

    r4e10f5 r192f6e  
    3434class AtomDescriptor_impl;
    3535template<typename T> class AtomsCalculation;
    36 class Box;
    3736class config;
    3837class ManipulateAtomsProcess;
    39 class Matrix;
    4038class molecule;
    4139class MoleculeDescriptor;
    … …  
    127125   * get the domain size as a symmetric matrix (6 components)
    128126   */
    129   Box& getDomain();
    130 
    131   /**
    132    * Set the domain size from a matrix object
    133    *
    134    * Matrix needs to be symmetric
    135    */
    136   void setDomain(const Matrix &mat);
     127  double * getDomain();
    137128
    138129  /**
    … …  
    266257  periodentafel *periode;
    267258  config *configuration;
    268   Box *cell_size;
     259  static double *cell_size;
    269260  std::string defaultName;
    270261  class ThermoStatContainer *Thermostats;

  • src/analysis_bonds.cpp

    r4e10f5 r192f6e  
    1313#include "element.hpp"
    1414#include "info.hpp"
    15 #include "verbose.hpp"
    1615#include "log.hpp"
    1716#include "molecule.hpp"

  • src/analysis_correlation.cpp

    r4e10f5 r192f6e  
    99
    1010#include <iostream>
    11 #include <iomanip>
    1211
    1312#include "analysis_correlation.hpp"
    … …  
    2019#include "triangleintersectionlist.hpp"
    2120#include "vector.hpp"
    22 #include "Matrix.hpp"
    2321#include "verbose.hpp"
    2422#include "World.hpp"
    25 #include "Box.hpp"
    2623
    2724
    … …  
    3734  PairCorrelationMap *outmap = NULL;
    3835  double distance = 0.;
    39   Box &domain = World::getInstance().getDomain();
     36  double *domain = World::getInstance().getDomain();
    4037
    4138  if (molecules->ListOfMolecules.empty()) {
    … …  
    7875                for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
    7976                  if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    80                     distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
     77                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
    8178                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    8279                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    … …  
    138135  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    139136    if ((*MolWalker)->ActiveFlag) {
    140       Matrix FullMatrix = World::getInstance().getDomain().getM();
    141       Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     137      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     138      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    142139      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    143140      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    144141      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    145142        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    146         periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     143        periodicX = *(*iter)->node;
     144        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    147145        // go through every range in xyz and get distance
    148146        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    149147          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    150148            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    151               checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     149              checkX = Vector(n[0], n[1], n[2]) + periodicX;
     150              checkX.MatrixMultiplication(FullMatrix);
    152151              for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
    153152                if ((*MolOtherWalker)->ActiveFlag) {
    … …  
    158157                      for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
    159158                        if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    160                           periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
     159                          periodicOtherX = *(*runner)->node;
     160                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    161161                          // go through every range in xyz and get distance
    162162                          for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
    163163                            for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
    164164                              for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
    165                                 checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
     165                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
     166                                checkOtherX.MatrixMultiplication(FullMatrix);
    166167                                distance = checkX.distance(checkOtherX);
    167168                                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    … …  
    175176            }
    176177      }
     178      delete[](FullMatrix);
     179      delete[](FullInverseMatrix);
    177180    }
    178181  }
    … …  
    192195  CorrelationToPointMap *outmap = NULL;
    193196  double distance = 0.;
    194   Box &domain = World::getInstance().getDomain();
     197  double *cell_size = World::getInstance().getDomain();
    195198
    196199  if (molecules->ListOfMolecules.empty()) {
    … …  
    208211        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    209212          if ((*type == NULL) || ((*iter)->type == *type)) {
    210             distance = domain.periodicDistance(*(*iter)->node,*point);
     213            distance = (*iter)->node->PeriodicDistance(*point, cell_size);
    211214            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    212215            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
    … …  
    243246  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    244247    if ((*MolWalker)->ActiveFlag) {
    245       Matrix FullMatrix = World::getInstance().getDomain().getM();
    246       Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     248      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     249      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    247250      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    248251      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    250253        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    251254          if ((*type == NULL) || ((*iter)->type == *type)) {
    252             periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     255            periodicX = *(*iter)->node;
     256            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    253257            // go through every range in xyz and get distance
    254258            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    255259              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    256260                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    257                   checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     261                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     262                  checkX.MatrixMultiplication(FullMatrix);
    258263                  distance = checkX.distance(*point);
    259264                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    … …  
    262267          }
    263268      }
     269      delete[](FullMatrix);
     270      delete[](FullInverseMatrix);
    264271    }
    265272
    … …  
    345352  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    346353    if ((*MolWalker)->ActiveFlag) {
    347       Matrix FullMatrix = World::getInstance().getDomain().getM();
    348       Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     354      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     355      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    349356      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    350357      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    352359        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    353360          if ((*type == NULL) || ((*iter)->type == *type)) {
    354             periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     361            periodicX = *(*iter)->node;
     362            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    355363            // go through every range in xyz and get distance
    356364            ShortestDistance = -1.;
    … …  
    358366              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    359367                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    360                   checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     368                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     369                  checkX.MatrixMultiplication(FullMatrix);
    361370                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
    362371                  distance = Intersections.GetSmallestDistance();
    … …  
    372381          }
    373382      }
     383      delete[](FullMatrix);
     384      delete[](FullInverseMatrix);
    374385    }
    375386

  • src/analysis_correlation.hpp

    r4e10f5 r192f6e  
    2626
    2727#include "atom.hpp"
    28 #include "verbose.hpp"
    2928
    3029/****************************************** forward declarations *****************************/

  • src/analyzer.cpp

    r4e10f5 r192f6e  
    1414#include "datacreator.hpp"
    1515#include "helpers.hpp"
     16#include "memoryallocator.hpp"
    1617#include "parser.hpp"
    1718#include "periodentafel.hpp"
    18 #include "verbose.hpp"
    1919
    2020// include config.h

  • src/atom.cpp

    r4e10f5 r192f6e  
    1212#include "element.hpp"
    1313#include "lists.hpp"
     14#include "memoryallocator.hpp"
    1415#include "parser.hpp"
    1516#include "vector.hpp"
    1617#include "World.hpp"
    1718#include "molecule.hpp"
    18 #include "Shapes/Shape.hpp"
    19 
    20 #include <iomanip>
    2119
    2220/************************************* Functions for class atom *************************************/
    … …  
    114112 * \return true - is inside, false - is not
    115113 */
    116 bool atom::IsInShape(const Shape& shape) const
    117 {
    118   return shape.isInside(*node);
     114bool atom::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
     115{
     116  return (node->IsInParallelepiped(offset, parallelepiped));
    119117};
    120118

  • src/atom.hpp

    r4e10f5 r192f6e  
    1818#endif
    1919
    20 #include <iosfwd>
     20#include <iostream>
    2121#include <list>
    2222#include <vector>
    … …  
    3535class World;
    3636class molecule;
    37 class Shape;
    3837
    3938/********************************************** declarations *******************************/
    … …  
    6766  double DistanceToVector(const Vector &origin) const;
    6867  double DistanceSquaredToVector(const Vector &origin) const;
    69   bool IsInShape(const Shape&) const;
     68  bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
    7069
    7170  // getter and setter

  • src/atom_graphnodeinfo.cpp

    r4e10f5 r192f6e  
    99
    1010#include "atom_graphnodeinfo.hpp"
     11#include "memoryallocator.hpp"
    1112
    1213/** Constructor of class GraphNodeInfo.
    1314 */
    14 GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(0), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(0) {};
     15GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
    1516
    1617/** Destructor of class GraphNodeInfo.

  • src/atom_particleinfo.cpp

    r4e10f5 r192f6e  
    99
    1010#include "atom_particleinfo.hpp"
    11 
    12 #include <iostream>
     11#include "memoryallocator.hpp"
    1312
    1413/** Constructor of ParticleInfo.

  • src/atom_particleinfo.hpp

    r4e10f5 r192f6e  
    1818#endif
    1919
    20 #include <iosfwd>
    21 #include <string>
     20#include <iostream>
    2221
    2322/****************************************** forward declarations *****************************/

  • src/atom_trajectoryparticle.hpp

    r4e10f5 r192f6e  
    2020#include <fstream>
    2121
    22 #include <gsl/gsl_inline.h>
    2322#include <gsl/gsl_randist.h>
    2423

  • src/bond.cpp

    r4e10f5 r192f6e  
    77#include "Helpers/MemDebug.hpp"
    88
    9 #include "verbose.hpp"
    109#include "atom.hpp"
    1110#include "bond.hpp"

  • src/bondgraph.cpp

    r4e10f5 r192f6e  
    1515#include "element.hpp"
    1616#include "info.hpp"
    17 #include "verbose.hpp"
    1817#include "log.hpp"
    1918#include "molecule.hpp"

  • src/bondgraph.hpp

    r4e10f5 r192f6e  
    1818#endif
    1919
    20 #include <iosfwd>
     20#include <iostream>
    2121
    2222/*********************************************** defines ***********************************/

  • src/boundary.cpp

    r4e10f5 r192f6e  
    1515#include "info.hpp"
    1616#include "linkedcell.hpp"
    17 #include "verbose.hpp"
    1817#include "log.hpp"
     18#include "memoryallocator.hpp"
    1919#include "molecule.hpp"
    2020#include "tesselation.hpp"
    … …  
    2222#include "World.hpp"
    2323#include "Plane.hpp"
    24 #include "Matrix.hpp"
    25 #include "Box.hpp"
    26 
    27 #include <iostream>
    28 #include <iomanip>
    2924
    3025#include<gsl/gsl_poly.h>
    … …  
    774769  int N[NDIM];
    775770  int n[NDIM];
    776   const Matrix &M = World::getInstance().getDomain().getM();
    777   Matrix Rotations;
    778   const Matrix &MInverse = World::getInstance().getDomain().getMinv();
     771  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     772  double Rotations[NDIM*NDIM];
     773  double *MInverse = InverseMatrix(M);
    779774  Vector AtomTranslations;
    780775  Vector FillerTranslations;
    … …  
    809804
    810805  // calculate filler grid in [0,1]^3
    811   FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
     806  FillerDistance = Vector(distance[0], distance[1], distance[2]);
     807  FillerDistance.InverseMatrixMultiplication(M);
    812808  for(int i=0;i<NDIM;i++)
    813809    N[i] = (int) ceil(1./FillerDistance[i]);
    … …  
    822818      for (n[2] = 0; n[2] < N[2]; n[2]++) {
    823819        // calculate position of current grid vector in untransformed box
    824         CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     820        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     821        CurrentPosition.MatrixMultiplication(M);
    825822        // create molecule random translation vector ...
    826823        for (int i=0;i<NDIM;i++)
    … …  
    843840            }
    844841
    845             Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
    846             Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
    847             Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
    848             Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
    849             Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
    850             Rotations.set(1,2,              sin(phi[1])                                                    );
    851             Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
    852             Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
    853             Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
     842            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
     843            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
     844            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
     845            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
     846            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
     847            Rotations[7] =               sin(phi[1])                                                ;
     848            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
     849            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
     850            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
    854851          }
    855852
    … …  
    857854          Inserter = (*iter)->x;
    858855          if (DoRandomRotation)
    859             Inserter *= Rotations;
     856            Inserter.MatrixMultiplication(Rotations);
    860857          Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
    861858
    862859          // check whether inserter is inside box
    863           Inserter *= MInverse;
     860          Inserter.MatrixMultiplication(MInverse);
    864861          FillIt = true;
    865862          for (int i=0;i<NDIM;i++)
    866863            FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
    867           Inserter *= M;
     864          Inserter.MatrixMultiplication(M);
    868865
    869866          // Check whether point is in- or outside
    … …  
    904901    delete TesselStruct[(*ListRunner)];
    905902  }
     903  delete[](M);
     904  delete[](MInverse);
    906905
    907906  return Filling;

  • src/boundary.hpp

    r4e10f5 r192f6e  
    1212
    1313#include <fstream>
    14 #include <iosfwd>
     14#include <iostream>
    1515
    1616// STL headers

  • src/builder.cpp

    r4e10f5 r192f6e  
    6060#include "config.hpp"
    6161#include "log.hpp"
     62#include "memoryusageobserver.hpp"
    6263#include "molecule.hpp"
    6364#include "periodentafel.hpp"

  • src/config.cpp

    r4e10f5 r192f6e  
    1919#include "info.hpp"
    2020#include "lists.hpp"
    21 #include "verbose.hpp"
    2221#include "log.hpp"
    2322#include "molecule.hpp"
     23#include "memoryallocator.hpp"
    2424#include "molecule.hpp"
    2525#include "periodentafel.hpp"
    2626#include "ThermoStatContainer.hpp"
    2727#include "World.hpp"
    28 #include "Matrix.hpp"
    29 #include "Box.hpp"
    3028
    3129/************************************* Functions for class config ***************************/
    … …  
    681679  // Unit cell and magnetic field
    682680  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    683   double * cell_size = new double[6];
     681  double * const cell_size = World::getInstance().getDomain();
    684682  cell_size[0] = BoxLength[0];
    685683  cell_size[1] = BoxLength[3];
    … …  
    688686  cell_size[4] = BoxLength[7];
    689687  cell_size[5] = BoxLength[8];
    690   World::getInstance().setDomain(cell_size);
    691   delete cell_size;
    692688  //if (1) fprintf(stderr,"\n");
    693689
    … …  
    887883
    888884  ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    889   double * cell_size = new double[6];
     885  double * const cell_size = World::getInstance().getDomain();
    890886  cell_size[0] = BoxLength[0];
    891887  cell_size[1] = BoxLength[3];
    … …  
    894890  cell_size[4] = BoxLength[7];
    895891  cell_size[5] = BoxLength[8];
    896   World::getInstance().setDomain(cell_size);
    897   delete[] cell_size;
    898892  if (1) fprintf(stderr,"\n");
    899893  config::DoPerturbation = 0;
    … …  
    10331027  // bring MaxTypes up to date
    10341028  mol->CountElements();
    1035   const Matrix &domain = World::getInstance().getDomain().getM();
     1029  const double * const cell_size = World::getInstance().getDomain();
    10361030  ofstream * const output = new ofstream(filename, ios::out);
    10371031  if (output != NULL) {
    … …  
    11041098    *output << endl;
    11051099    *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    1106     *output << domain.at(0,0) << "\t" << endl;
    1107     *output << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
    1108     *output << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     1100    *output << cell_size[0] << "\t" << endl;
     1101    *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     1102    *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
    11091103    // FIXME
    11101104    *output << endl;

  • src/datacreator.cpp

    r4e10f5 r192f6e  
    1212#include "helpers.hpp"
    1313#include "parser.hpp"
    14 #include "verbose.hpp"
    15 
    16 #include <iomanip>
    1714
    1815//=========================== FUNCTIONS============================

  • src/datacreator.hpp

    r4e10f5 r192f6e  
    1010using namespace std;
    1111
    12 #include <iosfwd>
     12#include <iostream>
    1313
    1414/****************************************** forward declarations *****************************/

  • src/element.hpp

    r4e10f5 r192f6e  
    1616#endif
    1717
    18 #include <iosfwd>
     18#include <iostream>
    1919#include <string>
    2020

  • src/ellipsoid.cpp

    r4e10f5 r192f6e  
    2121#include "tesselation.hpp"
    2222#include "vector.hpp"
    23 #include "Matrix.hpp"
    2423#include "verbose.hpp"
    2524
    … …  
    3534  Vector helper, RefPoint;
    3635  double distance = -1.;
    37   Matrix Matrix;
     36  double Matrix[NDIM*NDIM];
    3837  double InverseLength[3];
    3938  double psi,theta,phi; // euler angles in ZX'Z'' convention
    … …  
    5251  theta = EllipsoidAngle[1];
    5352  phi = EllipsoidAngle[2];
    54   Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
    55   Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
    56   Matrix.set(2,0, sin(psi)*sin(theta));
    57   Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
    58   Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
    59   Matrix.set(2,1, -cos(psi)*sin(theta));
    60   Matrix.set(0,2, sin(theta)*sin(phi));
    61   Matrix.set(1,2, sin(theta)*cos(phi));
    62   Matrix.set(2,2, cos(theta));
    63   helper *= Matrix;
     53  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
     54  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
     55  Matrix[2] = sin(psi)*sin(theta);
     56  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
     57  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
     58  Matrix[5] = -cos(psi)*sin(theta);
     59  Matrix[6] = sin(theta)*sin(phi);
     60  Matrix[7] = sin(theta)*cos(phi);
     61  Matrix[8] = cos(theta);
     62  helper.MatrixMultiplication(Matrix);
    6463  helper.ScaleAll(InverseLength);
    6564  //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
    … …  
    7473  phi = -EllipsoidAngle[2];
    7574  helper.ScaleAll(EllipsoidLength);
    76   Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
    77   Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
    78   Matrix.set(2,0, sin(psi)*sin(theta));
    79   Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
    80   Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
    81   Matrix.set(2,1, -cos(psi)*sin(theta));
    82   Matrix.set(0,2, sin(theta)*sin(phi));
    83   Matrix.set(1,2, sin(theta)*cos(phi));
    84   Matrix.set(2,2, cos(theta));
    85   helper *= Matrix;
     75  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
     76  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
     77  Matrix[2] = sin(psi)*sin(theta);
     78  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
     79  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
     80  Matrix[5] = -cos(psi)*sin(theta);
     81  Matrix[6] = sin(theta)*sin(phi);
     82  Matrix[7] = sin(theta)*cos(phi);
     83  Matrix[8] = cos(theta);
     84  helper.MatrixMultiplication(Matrix);
    8685  //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;
    8786

  • src/errorlogger.cpp

    r4e10f5 r192f6e  
    99
    1010#include <fstream>
    11 #include <iostream>
    1211#include "errorlogger.hpp"
    1312#include "verbose.hpp"

  • src/errorlogger.hpp

    r4e10f5 r192f6e  
    99#define ERRORLOGGER_HPP_
    1010
    11 #include <iosfwd>
     11#include <iostream>
    1212
    1313#include "Patterns/Singleton.hpp"

  • src/graph.cpp

    r4e10f5 r192f6e  
    1313#include "config.hpp"
    1414#include "graph.hpp"
    15 #include "verbose.hpp"
    1615#include "log.hpp"
    1716#include "molecule.hpp"

  • src/graph.hpp

    r4e10f5 r192f6e  
    2727class molecule;
    2828
     29class Graph;
    2930class SubGraph;
    3031class Node;
    … …  
    3334/********************************************** definitions *********************************/
    3435
    35 typedef std::pair < int, class Node* > NodeMap;
    36 typedef std::multimap < class Node*, class Edge* > EdgeMap;
     36#define NodeMap pair < int, class Node* >
     37#define EdgeMap multimap < class Node*, class Edge* >
    3738
    38 typedef std::deque<int> KeyStack;
    39 typedef std::set<int> KeySet;
    40 typedef std::pair<int, double> NumberValuePair;
     39#define KeyStack deque<int>
     40#define KeySet set<int>
     41#define NumberValuePair pair<int, double>
     42#define Graph map <KeySet, NumberValuePair, KeyCompare >
     43#define GraphPair pair <KeySet, NumberValuePair >
     44#define KeySetTestPair pair<KeySet::iterator, bool>
     45#define GraphTestPair pair<Graph::iterator, bool>
    4146
    42 // needed for definition of Graph and GraphTestPair
     47
     48/******************************** Some small functions and/or structures **********************************/
     49
    4350struct KeyCompare
    4451{
    4552  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
    4653};
    47 
    48 typedef std::map <KeySet, NumberValuePair, KeyCompare > Graph;
    49 typedef std::pair <KeySet, NumberValuePair > GraphPair;
    50 typedef std::pair<KeySet::iterator, bool> KeySetTestPair;
    51 typedef std::pair<Graph::iterator, bool> GraphTestPair;
    52 
    53 
    54 /******************************** Some small functions and/or structures **********************************/
    5554
    5655//bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)

  • src/gslvector.cpp

    r4e10f5 r192f6e  
    1010#include <cassert>
    1111#include <cmath>
    12 #include <iostream>
    1312
    1413#include "gslvector.hpp"
    1514#include "defs.hpp"
    1615#include "vector.hpp"
    17 #include "VectorContent.hpp"
    1816
    1917/** Constructor of class GSLVector.
    … …  
    7169 */
    7270void GSLVector::SetFromVector(Vector &v){
    73   gsl_vector_memcpy (vector, v.get()->content);
     71  gsl_vector_memcpy (vector, v.get());
    7472}
    7573

  • src/gslvector.hpp

    r4e10f5 r192f6e  
    1818#endif
    1919
    20 #include <iosfwd>
     20#include <iostream>
    2121#include <gsl/gsl_vector.h>
    2222

  • src/helpers.cpp

    r4e10f5 r192f6e  
    88#include "helpers.hpp"
    99#include "Helpers/fast_functions.hpp"
    10 #include "verbose.hpp"
    1110#include "log.hpp"
    12 
    13 #include <iostream>
     11#include "memoryusageobserver.hpp"
    1412
    1513/********************************************** helpful functions *********************************/
    … …  
    119117};
    120118
     119/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
     120 * \param *symm 6-dim array of unique symmetric matrix components
     121 * \return allocated NDIM*NDIM array with the symmetric matrix
     122 */
     123double * ReturnFullMatrixforSymmetric(const double * const symm)
     124{
     125  double *matrix = new double[NDIM * NDIM];
     126  matrix[0] = symm[0];
     127  matrix[1] = symm[1];
     128  matrix[2] = symm[3];
     129  matrix[3] = symm[1];
     130  matrix[4] = symm[2];
     131  matrix[5] = symm[4];
     132  matrix[6] = symm[3];
     133  matrix[7] = symm[4];
     134  matrix[8] = symm[5];
     135  return matrix;
     136};
     137
     138/** Calculate the inverse of a 3x3 matrix.
     139 * \param *matrix NDIM_NDIM array
     140 */
     141double * InverseMatrix( const double * const A)
     142{
     143  double *B = new double[NDIM * NDIM];
     144  double detA = RDET3(A);
     145  double detAReci;
     146
     147  for (int i=0;i<NDIM*NDIM;++i)
     148    B[i] = 0.;
     149  // calculate the inverse B
     150  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
     151    detAReci = 1./detA;
     152    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
     153    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
     154    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
     155    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
     156    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
     157    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
     158    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
     159    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     160    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
     161  }
     162  return B;
     163};
     164
     165
     166
    121167/** Comparison function for GSL heapsort on distances in two molecules.
    122168 * \param *a

  • src/helpers.hpp

    r4e10f5 r192f6e  
    2020#include "defs.hpp"
    2121#include "log.hpp"
     22#include "memoryallocator.hpp"
    2223
    2324/********************************************** definitions *********************************/
    … …  
    5152bool IsValidNumber( const char *string);
    5253int CompareDoubles (const void * a, const void * b);
     54double * ReturnFullMatrixforSymmetric(const double * const cell_size);
     55double * InverseMatrix(const double * const A);
    5356void performCriticalExit();
    5457

  • src/joiner.cpp

    r4e10f5 r192f6e  
    1414#include "datacreator.hpp"
    1515#include "helpers.hpp"
     16#include "memoryallocator.hpp"
    1617#include "parser.hpp"
    1718#include "periodentafel.hpp"
    18 #include "verbose.hpp"
    1919
    2020//============================== MAIN =============================

  • src/linkedcell.cpp

    r4e10f5 r192f6e  
    1010#include "helpers.hpp"
    1111#include "linkedcell.hpp"
    12 #include "verbose.hpp"
    1312#include "log.hpp"
    1413#include "molecule.hpp"

  • src/logger.cpp

    r4e10f5 r192f6e  
    99
    1010#include <fstream>
    11 #include <iostream>
    1211#include "logger.hpp"
    1312#include "verbose.hpp"

  • src/logger.hpp

    r4e10f5 r192f6e  
    99#define LOGGER_HPP_
    1010
    11 #include <iosfwd>
     11#include <iostream>
    1212
    1313#include "Patterns/Singleton.hpp"

  • src/molecule.cpp

    r4e10f5 r192f6e  
    55 */
    66
    7 #ifdef HAVE_CONFIG_H
    8 #include <config.h>
    9 #endif
    10 
    117#include "Helpers/MemDebug.hpp"
    128
    139#include <cstring>
    1410#include <boost/bind.hpp>
    15 #include <boost/foreach.hpp>
    16 
    17 #include <gsl/gsl_inline.h>
    18 #include <gsl/gsl_heapsort.h>
    1911
    2012#include "World.hpp"
    … …  
    3022#include "log.hpp"
    3123#include "molecule.hpp"
    32 
     24#include "memoryallocator.hpp"
    3325#include "periodentafel.hpp"
    3426#include "stackclass.hpp"
    3527#include "tesselation.hpp"
    3628#include "vector.hpp"
    37 #include "Matrix.hpp"
    3829#include "World.hpp"
    39 #include "Box.hpp"
    4030#include "Plane.hpp"
    4131#include "Exceptions/LinearDependenceException.hpp"
    … …  
    294284  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
    295285  Vector InBondvector;    // vector in direction of *Bond
    296   const Matrix &matrix =  World::getInstance().getDomain().getM();
     286  double *matrix = NULL;
    297287  bond *Binder = NULL;
     288  double * const cell_size = World::getInstance().getDomain();
    298289
    299290//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
    … …  
    316307      } // (signs are correct, was tested!)
    317308    }
    318     Orthovector1 *= matrix;
     309    matrix = ReturnFullMatrixforSymmetric(cell_size);
     310    Orthovector1.MatrixMultiplication(matrix);
    319311    InBondvector -= Orthovector1; // subtract just the additional translation
     312    delete[](matrix);
    320313    bondlength = InBondvector.Norm();
    321314//    Log() << Verbose(4) << "Corrected InBondvector is now: ";
    … …  
    548541      break;
    549542  }
     543  delete[](matrix);
    550544
    551545//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
    … …  
    666660 * @param three vectors forming the matrix that defines the shape of the parallelpiped
    667661 */
    668 molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
     662molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
    669663  molecule *copy = World::getInstance().createMolecule();
    670664
    671   BOOST_FOREACH(atom *iter,atoms){
    672     if(iter->IsInShape(region)){
    673       copy->AddCopyAtom(iter);
    674     }
    675   }
     665  ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
    676666
    677667  //TODO: copy->BuildInducedSubgraph(this);
    … …  
    760750void molecule::SetBoxDimension(Vector *dim)
    761751{
    762   Matrix domain;
    763   for(int i =0; i<NDIM;++i)
    764     domain.at(i,i) = dim->at(i);
    765   World::getInstance().setDomain(domain);
     752  double * const cell_size = World::getInstance().getDomain();
     753  cell_size[0] = dim->at(0);
     754  cell_size[1] = 0.;
     755  cell_size[2] = dim->at(1);
     756  cell_size[3] = 0.;
     757  cell_size[4] = 0.;
     758  cell_size[5] = dim->at(2);
    766759};
    767760
    … …  
    856849bool molecule::CheckBounds(const Vector *x) const
    857850{
    858   const Matrix &domain = World::getInstance().getDomain().getM();
     851  double * const cell_size = World::getInstance().getDomain();
    859852  bool result = true;
     853  int j =-1;
    860854  for (int i=0;i<NDIM;i++) {
    861     result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
     855    j += i+1;
     856    result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
    862857  }
    863858  //return result;

  • src/molecule.hpp

    r4e10f5 r192f6e  
    77#define MOLECULES_HPP_
    88
     9using namespace std;
     10
    911/*********************************************** includes ***********************************/
    1012
    11 #ifdef HAVE_CONFIG_H
    12 #include <config.h>
    13 #endif
     13// GSL headers
     14#include <gsl/gsl_eigen.h>
     15#include <gsl/gsl_heapsort.h>
     16#include <gsl/gsl_linalg.h>
     17#include <gsl/gsl_matrix.h>
     18#include <gsl/gsl_multimin.h>
     19#include <gsl/gsl_vector.h>
     20#include <gsl/gsl_randist.h>
    1421
    1522//// STL headers
    … …  
    2431#include "types.hpp"
    2532#include "graph.hpp"
     33#include "stackclass.hpp"
    2634#include "tesselation.hpp"
    2735#include "Patterns/Observer.hpp"
    … …  
    4553class periodentafel;
    4654class Vector;
    47 class Shape;
    48 template <class> class StackClass;
    4955
    5056/******************************** Some definitions for easier reading **********************************/
    … …  
    310316
    311317  molecule *CopyMolecule();
    312   molecule* CopyMoleculeFromSubRegion(const Shape&) const;
     318  molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;
    313319
    314320  /// Fragment molecule by two different approaches:

  • src/molecule_dynamics.cpp

    r4e10f5 r192f6e  
    1313#include "element.hpp"
    1414#include "info.hpp"
    15 #include "verbose.hpp"
    1615#include "log.hpp"
     16#include "memoryallocator.hpp"
    1717#include "molecule.hpp"
    1818#include "parser.hpp"
    1919#include "Plane.hpp"
    2020#include "ThermoStatContainer.hpp"
    21 
    22 #include <gsl/gsl_matrix.h>
    23 #include <gsl/gsl_vector.h>
    24 #include <gsl/gsl_linalg.h>
    2521
    2622/************************************* Functions for class molecule *********************************/

  • src/molecule_fragmentation.cpp

    r4e10f5 r192f6e  
    1717#include "helpers.hpp"
    1818#include "lists.hpp"
    19 #include "verbose.hpp"
    2019#include "log.hpp"
     20#include "memoryallocator.hpp"
    2121#include "molecule.hpp"
    2222#include "periodentafel.hpp"
    2323#include "World.hpp"
    24 #include "Matrix.hpp"
    25 #include "Box.hpp"
    26 #include "stackclass.hpp"
    2724
    2825/************************************* Functions for class molecule *********************************/
    … …  
    17191716  atom *Walker = NULL;
    17201717  atom *OtherWalker = NULL;
    1721   Matrix matrix = World::getInstance().getDomain().getM();
     1718  double * const cell_size = World::getInstance().getDomain();
     1719  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    17221720  enum Shading *ColorList = NULL;
    17231721  double tmp;
    … …  
    17591757          Translationvector[i] = (tmp < 0) ? +1. : -1.;
    17601758      }
    1761       Translationvector *= matrix;
     1759      Translationvector.MatrixMultiplication(matrix);
    17621760      //Log() << Verbose(3) << "Translation vector is ";
    17631761      Log() << Verbose(0) << Translationvector <<  endl;
    … …  
    17901788  delete(AtomStack);
    17911789  delete[](ColorList);
     1790  delete[](matrix);
    17921791  DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
    17931792};

  • src/molecule_geometry.cpp

    r4e10f5 r192f6e  
    55 *      Author: heber
    66 */
    7 
    8 #ifdef HAVE_CONFIG_H
    9 #include <config.h>
    10 #endif
    117
    128#include "Helpers/MemDebug.hpp"
    … …  
    1814#include "helpers.hpp"
    1915#include "leastsquaremin.hpp"
    20 #include "verbose.hpp"
    2116#include "log.hpp"
     17#include "memoryallocator.hpp"
    2218#include "molecule.hpp"
    2319#include "World.hpp"
    2420#include "Plane.hpp"
    25 #include "Matrix.hpp"
    26 #include "Box.hpp"
    2721#include <boost/foreach.hpp>
    28 
    29 #include <gsl/gsl_eigen.h>
    30 #include <gsl/gsl_multimin.h>
    3122
    3223
    … …  
    4233  const Vector *Center = DetermineCenterOfAll();
    4334  const Vector *CenterBox = DetermineCenterOfBox();
    44   Box &domain = World::getInstance().getDomain();
     35  double * const cell_size = World::getInstance().getDomain();
     36  double *M = ReturnFullMatrixforSymmetric(cell_size);
     37  double *Minv = InverseMatrix(M);
    4538
    4639  // go through all atoms
    4740  ActOnAllVectors( &Vector::SubtractVector, *Center);
    4841  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
    49   BOOST_FOREACH(atom* iter, atoms){
    50     *iter->node = domain.WrapPeriodically(*iter->node);
    51   }
    52 
     42  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
     43
     44  delete[](M);
     45  delete[](Minv);
    5346  delete(Center);
    5447  delete(CenterBox);
    … …  
    6356{
    6457  bool status = true;
    65   Box &domain = World::getInstance().getDomain();
     58  double * const cell_size = World::getInstance().getDomain();
     59  double *M = ReturnFullMatrixforSymmetric(cell_size);
     60  double *Minv = InverseMatrix(M);
    6661
    6762  // go through all atoms
    68   BOOST_FOREACH(atom* iter, atoms){
    69     *iter->node = domain.WrapPeriodically(*iter->node);
    70   }
    71 
     63  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
     64
     65  delete[](M);
     66  delete[](Minv);
    7267  return status;
    7368};
    … …  
    158153{
    159154  Vector *a = new Vector(0.5,0.5,0.5);
    160   const Matrix &M = World::getInstance().getDomain().getM();
    161   (*a) *= M;
     155
     156  const double *cell_size = World::getInstance().getDomain();
     157  double *M = ReturnFullMatrixforSymmetric(cell_size);
     158  a->MatrixMultiplication(M);
     159  delete[](M);
     160
    162161  return a;
    163162};
    … …  
    245244void molecule::TranslatePeriodically(const Vector *trans)
    246245{
    247   Box &domain = World::getInstance().getDomain();
     246  double * const cell_size = World::getInstance().getDomain();
     247  double *M = ReturnFullMatrixforSymmetric(cell_size);
     248  double *Minv = InverseMatrix(M);
    248249
    249250  // go through all atoms
    250251  ActOnAllVectors( &Vector::AddVector, *trans);
    251   BOOST_FOREACH(atom* iter, atoms){
    252     *iter->node = domain.WrapPeriodically(*iter->node);
    253   }
    254 
     252  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
     253
     254  delete[](M);
     255  delete[](Minv);
    255256};
    256257
    … …  
    263264  OBSERVE;
    264265  Plane p(*n,0);
    265   BOOST_FOREACH(atom* iter, atoms ){
     266  BOOST_FOREACH( atom* iter, atoms ){
    266267    (*iter->node) = p.mirrorVector(*iter->node);
    267268  }
    … …  
    273274void molecule::DeterminePeriodicCenter(Vector &center)
    274275{
    275   const Matrix &matrix = World::getInstance().getDomain().getM();
    276   const Matrix &inversematrix = World::getInstance().getDomain().getM();
     276  double * const cell_size = World::getInstance().getDomain();
     277  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     278  double *inversematrix = InverseMatrix(matrix);
    277279  double tmp;
    278280  bool flag;
    … …  
    286288      if ((*iter)->type->Z != 1) {
    287289#endif
    288         Testvector = inversematrix * (*iter)->x;
     290        Testvector = (*iter)->x;
     291        Testvector.MatrixMultiplication(inversematrix);
    289292        Translationvector.Zero();
    290293        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    … …  
    303306        }
    304307        Testvector += Translationvector;
    305         Testvector *= matrix;
     308        Testvector.MatrixMultiplication(matrix);
    306309        Center += Testvector;
    307310        Log() << Verbose(1) << "vector is: " << Testvector << endl;
    … …  
    310313        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    311314          if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
    312             Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
     315            Testvector = (*Runner)->GetOtherAtom((*iter))->x;
     316            Testvector.MatrixMultiplication(inversematrix);
    313317            Testvector += Translationvector;
    314             Testvector *= matrix;
     318            Testvector.MatrixMultiplication(matrix);
    315319            Center += Testvector;
    316320            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
    … …  
    321325    }
    322326  } while (!flag);
     327  delete[](matrix);
     328  delete[](inversematrix);
    323329
    324330  Center.Scale(1./static_cast<double>(getAtomCount()));
    … …  
    382388    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    383389    // the eigenvectors specify the transformation matrix
    384     Matrix M = Matrix(evec->data);
    385     BOOST_FOREACH(atom* iter, atoms){
    386       (*iter->node) *= M;
    387     }
     390    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
    388391    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    389392

  • src/molecule_graph.cpp

    r4e10f5 r192f6e  
    55 *      Author: heber
    66 */
    7 
    8 #ifdef HAVE_CONFIG_H
    9 #include <config.h>
    10 #endif
    117
    128#include "Helpers/MemDebug.hpp"
    … …  
    2218#include "linkedcell.hpp"
    2319#include "lists.hpp"
    24 #include "verbose.hpp"
    2520#include "log.hpp"
     21#include "memoryallocator.hpp"
    2622#include "molecule.hpp"
    2723#include "World.hpp"
    2824#include "Helpers/fast_functions.hpp"
    2925#include "Helpers/Assert.hpp"
    30 #include "Matrix.hpp"
    31 #include "Box.hpp"
    32 #include "stackclass.hpp"
     26
    3327
    3428struct BFSAccounting
    … …  
    127121  LinkedCell *LC = NULL;
    128122  bool free_BG = false;
    129   Box &domain = World::getInstance().getDomain();
     123  double * const cell_size = World::getInstance().getDomain();
    130124
    131125  if (BG == NULL) {
    … …  
    184178                          //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
    185179                          (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
    186                           const double distance = domain.periodicDistanceSquared(OtherWalker->x,Walker->x);
     180                          const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);
    187181                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
    188182//                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
    … …  
    585579    // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
    586580    DoLog(0) && (Log() << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << DFS.CurrentGraphNr << "." << endl);
    587     LeafWalker->Leaf->Output((ofstream *)&(Log() << Verbose(0)));
     581    LeafWalker->Leaf->Output((ofstream *)&cout);
    588582    DoLog(0) && (Log() << Verbose(0) << endl);
    589583

  • src/molecule_pointcloud.cpp

    r4e10f5 r192f6e  
    1111#include "config.hpp"
    1212#include "info.hpp"
     13#include "memoryallocator.hpp"
    1314#include "molecule.hpp"
    1415

  • src/moleculelist.cpp

    r4e10f5 r192f6e  
    55 */
    66
    7 #ifdef HAVE_CONFIG_H
    8 #include <config.h>
    9 #endif
    10 
    117#include "Helpers/MemDebug.hpp"
    128
    139#include <cstring>
    14 
    15 #include <gsl/gsl_inline.h>
    16 #include <gsl/gsl_heapsort.h>
    1710
    1811#include "World.hpp"
    … …  
    2619#include "linkedcell.hpp"
    2720#include "lists.hpp"
    28 #include "verbose.hpp"
    2921#include "log.hpp"
    3022#include "molecule.hpp"
     23#include "memoryallocator.hpp"
    3124#include "periodentafel.hpp"
    3225#include "Helpers/Assert.hpp"
    33 #include "Matrix.hpp"
    34 #include "Box.hpp"
    35 #include "stackclass.hpp"
    3626
    3727#include "Helpers/Assert.hpp"
    … …  
    641631  int FragmentCounter = 0;
    642632  ofstream output;
    643   Matrix cell_size = World::getInstance().getDomain().getM();
    644   Matrix cell_size_backup = cell_size;
    645 
     633  double cell_size_backup[6];
     634  double * const cell_size = World::getInstance().getDomain();
     635
     636  // backup cell_size
     637  for (int i=0;i<6;i++)
     638    cell_size_backup[i] = cell_size[i];
    646639  // store the fragments as config and as xyz
    647640  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    … …  
    681674    (*ListRunner)->CenterEdge(&BoxDimension);
    682675    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
     676    int j = -1;
    683677    for (int k = 0; k < NDIM; k++) {
     678      j += k + 1;
    684679      BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
    685       cell_size.at(k,k) = BoxDimension[k] * 2.;
    686     }
    687     World::getInstance().setDomain(cell_size);
     680      cell_size[j] = BoxDimension[k] * 2.;
     681    }
    688682    (*ListRunner)->Translate(&BoxDimension);
    689683
    … …  
    731725
    732726  // restore cell_size
    733   World::getInstance().setDomain(cell_size_backup);
     727  for (int i=0;i<6;i++)
     728    cell_size[i] = cell_size_backup[i];
    734729
    735730  return result;
    … …  
    892887  // center at set box dimensions
    893888  mol->CenterEdge(&center);
    894   Matrix domain;
    895   for(int i =0;i<NDIM;++i)
    896     domain.at(i,i) = center[i];
    897   World::getInstance().setDomain(domain);
     889  World::getInstance().getDomain()[0] = center[0];
     890  World::getInstance().getDomain()[1] = 0;
     891  World::getInstance().getDomain()[2] = center[1];
     892  World::getInstance().getDomain()[3] = 0;
     893  World::getInstance().getDomain()[4] = 0;
     894  World::getInstance().getDomain()[5] = center[2];
    898895  insert(mol);
    899896}

  • src/parser.cpp

    r4e10f5 r192f6e  
    1212
    1313#include "helpers.hpp"
     14#include "memoryallocator.hpp"
    1415#include "parser.hpp"
    15 #include "verbose.hpp"
    1616
    1717// include config.h

  • src/periodentafel.hpp

    r4e10f5 r192f6e  
    99#endif
    1010
    11 #include <iosfwd>
     11#include <iostream>
    1212#include <map>
     13#include <iterator>
    1314
    1415#include "unittests/periodentafelTest.hpp"

  • src/stackclass.hpp

    r4e10f5 r192f6e  
    1212
    1313#include "verbose.hpp"
    14 #include "log.hpp"
     14#include "memoryallocator.hpp"
    1515
    1616/****************************************** forward declarations *****************************/

  • src/tesselation.cpp

    r4e10f5 r192f6e  
    99
    1010#include <fstream>
    11 #include <iomanip>
    1211
    1312#include "helpers.hpp"

  • src/tesselationhelpers.cpp

    r4e10f5 r192f6e  
    2121#include "verbose.hpp"
    2222#include "Plane.hpp"
    23 #include "Matrix.hpp"
     23
     24double DetGet(gsl_matrix * const A, const int inPlace)
     25{
     26        Info FunctionInfo(__func__);
     27  /*
     28  inPlace = 1 => A is replaced with the LU decomposed copy.
     29  inPlace = 0 => A is retained, and a copy is used for LU.
     30  */
     31
     32  double det;
     33  int signum;
     34  gsl_permutation *p = gsl_permutation_alloc(A->size1);
     35  gsl_matrix *tmpA=0;
     36
     37  if (inPlace)
     38  tmpA = A;
     39  else {
     40  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
     41  gsl_matrix_memcpy(tmpA , A);
     42  }
     43
     44
     45  gsl_linalg_LU_decomp(tmpA , p , &signum);
     46  det = gsl_linalg_LU_det(tmpA , signum);
     47  gsl_permutation_free(p);
     48  if (! inPlace)
     49  gsl_matrix_free(tmpA);
     50
     51  return det;
     52};
    2453
    2554void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
    2655{
    2756        Info FunctionInfo(__func__);
    28   Matrix mat;
     57  gsl_matrix *A = gsl_matrix_calloc(3,3);
    2958  double m11, m12, m13, m14;
    3059
    3160  for(int i=0;i<3;i++) {
    32     mat.set(i, 0, a[i]);
    33     mat.set(i, 1, b[i]);
    34     mat.set(i, 2, c[i]);
    35   }
    36   m11 = mat.determinant();
     61    gsl_matrix_set(A, i, 0, a[i]);
     62    gsl_matrix_set(A, i, 1, b[i]);
     63    gsl_matrix_set(A, i, 2, c[i]);
     64  }
     65  m11 = DetGet(A, 1);
    3766
    3867  for(int i=0;i<3;i++) {
    39     mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
    40     mat.set(i, 1, b[i]);
    41     mat.set(i, 2, c[i]);
    42   }
    43   m12 = mat.determinant();
     68    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     69    gsl_matrix_set(A, i, 1, b[i]);
     70    gsl_matrix_set(A, i, 2, c[i]);
     71  }
     72  m12 = DetGet(A, 1);
    4473
    4574  for(int i=0;i<3;i++) {
    46     mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
    47     mat.set(i, 1, a[i]);
    48     mat.set(i, 2, c[i]);
    49   }
    50   m13 = mat.determinant();
     75    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     76    gsl_matrix_set(A, i, 1, a[i]);
     77    gsl_matrix_set(A, i, 2, c[i]);
     78  }
     79  m13 = DetGet(A, 1);
    5180
    5281  for(int i=0;i<3;i++) {
    53     mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
    54     mat.set(i, 1, a[i]);
    55     mat.set(i, 2, b[i]);
    56   }
    57   m14 = mat.determinant();
     82    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     83    gsl_matrix_set(A, i, 1, a[i]);
     84    gsl_matrix_set(A, i, 2, b[i]);
     85  }
     86  m14 = DetGet(A, 1);
    5887
    5988  if (fabs(m11) < MYEPSILON)
    … …  
    6695  if (fabs(a.distance(*center) - RADIUS) > MYEPSILON)
    6796    DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.distance(*center) - RADIUS) << " from a than RADIUS." << endl);
     97
     98  gsl_matrix_free(A);
    6899};
    69100
    … …  
    217248  Vector x3;
    218249  Vector x4;
     250};
     251
     252/**
     253 * Intersection calculation function.
     254 *
     255 * @param x to find the result for
     256 * @param function parameter
     257 */
     258double MinIntersectDistance(const gsl_vector * x, void *params)
     259{
     260        Info FunctionInfo(__func__);
     261  double retval = 0;
     262  struct Intersection *I = (struct Intersection *)params;
     263  Vector intersection;
     264  for (int i=0;i<NDIM;i++)
     265    intersection[i] = gsl_vector_get(x, i);
     266
     267  Vector SideA = I->x1 -I->x2 ;
     268  Vector HeightA = intersection - I->x1;
     269  HeightA.ProjectOntoPlane(SideA);
     270
     271  Vector SideB = I->x3 - I->x4;
     272  Vector HeightB = intersection - I->x3;
     273  HeightB.ProjectOntoPlane(SideB);
     274
     275  retval = HeightA.ScalarProduct(HeightA) + HeightB.ScalarProduct(HeightB);
     276  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
     277
     278  return retval;
     279};
     280
     281
     282/**
     283 * Calculates whether there is an intersection between two lines. The first line
     284 * always goes through point 1 and point 2 and the second line is given by the
     285 * connection between point 4 and point 5.
     286 *
     287 * @param point 1 of line 1
     288 * @param point 2 of line 1
     289 * @param point 1 of line 2
     290 * @param point 2 of line 2
     291 *
     292 * @return true if there is an intersection between the given lines, false otherwise
     293 */
     294bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
     295{
     296        Info FunctionInfo(__func__);
     297  bool result;
     298
     299  struct Intersection par;
     300    par.x1 = point1;
     301    par.x2 = point2;
     302    par.x3 = point3;
     303    par.x4 = point4;
     304
     305    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
     306    gsl_multimin_fminimizer *s = NULL;
     307    gsl_vector *ss, *x;
     308    gsl_multimin_function minexFunction;
     309
     310    size_t iter = 0;
     311    int status;
     312    double size;
     313
     314    /* Starting point */
     315    x = gsl_vector_alloc(NDIM);
     316    gsl_vector_set(x, 0, point1[0]);
     317    gsl_vector_set(x, 1, point1[1]);
     318    gsl_vector_set(x, 2, point1[2]);
     319
     320    /* Set initial step sizes to 1 */
     321    ss = gsl_vector_alloc(NDIM);
     322    gsl_vector_set_all(ss, 1.0);
     323
     324    /* Initialize method and iterate */
     325    minexFunction.n = NDIM;
     326    minexFunction.f = &MinIntersectDistance;
     327    minexFunction.params = (void *)&par;
     328
     329    s = gsl_multimin_fminimizer_alloc(T, NDIM);
     330    gsl_multimin_fminimizer_set(s, &minexFunction, x, ss);
     331
     332    do {
     333        iter++;
     334        status = gsl_multimin_fminimizer_iterate(s);
     335
     336        if (status) {
     337          break;
     338        }
     339
     340        size = gsl_multimin_fminimizer_size(s);
     341        status = gsl_multimin_test_size(size, 1e-2);
     342
     343        if (status == GSL_SUCCESS) {
     344          DoLog(1) && (Log() << Verbose(1) << "converged to minimum" <<  endl);
     345        }
     346    } while (status == GSL_CONTINUE && iter < 100);
     347
     348    // check whether intersection is in between or not
     349  Vector intersection;
     350  double t1, t2;
     351  for (int i = 0; i < NDIM; i++) {
     352    intersection[i] = gsl_vector_get(s->x, i);
     353  }
     354
     355  Vector SideA = par.x2 - par.x1;
     356  Vector HeightA = intersection - par.x1;
     357
     358  t1 = HeightA.ScalarProduct(SideA)/SideA.ScalarProduct(SideA);
     359
     360  Vector SideB = par.x4 - par.x3;
     361  Vector HeightB = intersection - par.x3;
     362
     363  t2 = HeightB.ScalarProduct(SideB)/SideB.ScalarProduct(SideB);
     364
     365  Log() << Verbose(1) << "Intersection " << intersection << " is at "
     366    << t1 << " for (" << point1 << "," << point2 << ") and at "
     367    << t2 << " for (" << point3 << "," << point4 << "): ";
     368
     369  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
     370    DoLog(1) && (Log() << Verbose(1) << "true intersection." << endl);
     371    result = true;
     372  } else {
     373    DoLog(1) && (Log() << Verbose(1) << "intersection out of region of interest." << endl);
     374    result = false;
     375  }
     376
     377  // free minimizer stuff
     378    gsl_vector_free(x);
     379    gsl_vector_free(ss);
     380    gsl_multimin_fminimizer_free(s);
     381
     382  return result;
    219383};
    220384

  • src/tesselationhelpers.hpp

    r4e10f5 r192f6e  
    2020#endif
    2121
    22 #include <iosfwd>
     22#include <gsl/gsl_linalg.h>
     23#include <gsl/gsl_matrix.h>
     24#include <gsl/gsl_multimin.h>
     25#include <gsl/gsl_permutation.h>
     26#include <gsl/gsl_vector.h>
     27
     28#include <iostream>
    2329
    2430#include "defs.hpp"
    … …  
    4147/********************************************** declarations *******************************/
    4248
     49double DetGet(gsl_matrix * const A, const int inPlace);
    4350void GetSphere(Vector * const Center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS);
    4451void 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);
    4552void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c);
    4653double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection);
     54double MinIntersectDistance(const gsl_vector * x, void *params);
     55bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4);
    4756double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d);
    4857double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C);

  • src/triangleintersectionlist.cpp

    r4e10f5 r192f6e  
    1818#include "tesselation.hpp"
    1919#include "vector.hpp"
    20 #include "verbose.hpp"
    2120
    2221/** Constructor for class TriangleIntersectionList.

  • src/unittests/ActOnAllUnitTest.cpp

    r4e10f5 r192f6e  
    1414#include "../test/ActOnAlltest.hpp"
    1515#include "ActOnAllUnitTest.hpp"
     16#include "memoryallocator.hpp"
    1617#include "vector.hpp"
    1718

  • src/unittests/LineUnittest.cpp

    r4e10f5 r192f6e  
    1717
    1818#include <iostream>
    19 #include <cmath>
    2019
    2120using namespace std;
    … …  
    353352  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
    354353}
    355 
    356 void LineUnittest::sphereIntersectionTest(){
    357   {
    358     std::vector<Vector> res = la1->getSphereIntersections();
    359     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    360     CPPUNIT_ASSERT(testDirection(res[0],e1));
    361     CPPUNIT_ASSERT(testDirection(res[1],e1));
    362     CPPUNIT_ASSERT(res[0]!=res[1]);
    363   }
    364 
    365   {
    366     std::vector<Vector> res = la2->getSphereIntersections();
    367     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    368     CPPUNIT_ASSERT(testDirection(res[0],e2));
    369     CPPUNIT_ASSERT(testDirection(res[1],e2));
    370     CPPUNIT_ASSERT(res[0]!=res[1]);
    371   }
    372 
    373   {
    374     std::vector<Vector> res = la3->getSphereIntersections();
    375     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    376     CPPUNIT_ASSERT(testDirection(res[0],e3));
    377     CPPUNIT_ASSERT(testDirection(res[1],e3));
    378     CPPUNIT_ASSERT(res[0]!=res[1]);
    379   }
    380 
    381   {
    382     std::vector<Vector> res = lp1->getSphereIntersections();
    383     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    384     CPPUNIT_ASSERT((res[0]==e1) || (res[0]==e2));
    385     CPPUNIT_ASSERT((res[1]==e1) || (res[1]==e2));
    386     CPPUNIT_ASSERT(res[0]!=res[1]);
    387   }
    388 
    389   {
    390     std::vector<Vector> res = lp2->getSphereIntersections();
    391     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    392     CPPUNIT_ASSERT((res[0]==e2) || (res[0]==e3));
    393     CPPUNIT_ASSERT((res[1]==e2) || (res[1]==e3));
    394     CPPUNIT_ASSERT(res[0]!=res[1]);
    395   }
    396 
    397   {
    398     std::vector<Vector> res = lp3->getSphereIntersections();
    399     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    400     CPPUNIT_ASSERT((res[0]==e3) || (res[0]==e1));
    401     CPPUNIT_ASSERT((res[1]==e3) || (res[1]==e1));
    402     CPPUNIT_ASSERT(res[0]!=res[1]);
    403   }
    404 }

  • src/unittests/LineUnittest.hpp

    r4e10f5 r192f6e  
    2222  CPPUNIT_TEST ( intersectionTest );
    2323  CPPUNIT_TEST ( rotationTest );
    24   CPPUNIT_TEST ( sphereIntersectionTest );
    2524  CPPUNIT_TEST_SUITE_END();
    2625
    … …  
    3433  void intersectionTest();
    3534  void rotationTest();
    36   void sphereIntersectionTest();
    3735
    3836private:

  • src/unittests/Makefile.am

    r4e10f5 r192f6e  
    3131  LogUnitTest \
    3232  manipulateAtomsTest \
    33   MatrixUnittest \
     33  MemoryUsageObserverUnitTest \
     34  MemoryAllocatorUnitTest \
    3435  MoleculeDescriptorTest \
    3536  ObserverTest \
    … …  
    3738  periodentafelTest \
    3839  PlaneUnittest \
    39   ShapeUnittest \
    4040  SingletonTest \
    4141  StackClassUnitTest \
    … …  
    7575  listofbondsunittest.cpp \
    7676  logunittest.cpp \
    77   MatrixUnittest.cpp \
    7877  manipulateAtomsTest.cpp \
     78  memoryallocatorunittest.cpp  \
     79  memoryusageobserverunittest.cpp \
    7980  MoleculeDescriptorTest.cpp \
    8081  ObserverTest.cpp \
    … …  
    8283  periodentafelTest.cpp \
    8384  PlaneUnittest.cpp \
    84   ShapeUnittest.cpp \
    8585  SingletonTest.cpp \
    8686  stackclassunittest.cpp \
    … …  
    112112  logunittest.hpp \
    113113  manipulateAtomsTest.hpp \
    114   MatrixUnittest.hpp \
     114  memoryallocatorunittest.hpp  \
     115  memoryusageobserverunittest.hpp \
    115116  MoleculeDescriptorTest.hpp \
    116117  periodentafelTest.hpp \
    … …  
    188189manipulateAtomsTest_LDADD = ${ALLLIBS}
    189190
    190 MatrixUnittest_SOURCES = UnitTestMain.cpp MatrixUnittest.cpp MatrixUnittest.hpp
    191 MatrixUnittest_LDADD = ${ALLLIBS}
     191MemoryAllocatorUnitTest_SOURCES = UnitTestMain.cpp ../memoryallocator.hpp ../memoryallocator.cpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp memoryallocatorunittest.cpp memoryallocatorunittest.hpp
     192MemoryAllocatorUnitTest_LDADD = ${ALLLIBS}
     193
     194MemoryUsageObserverUnitTest_SOURCES = UnitTestMain.cpp ../memoryallocator.hpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp memoryusageobserverunittest.cpp memoryusageobserverunittest.hpp
     195MemoryUsageObserverUnitTest_LDADD = ${ALLLIBS}
    192196
    193197MoleculeDescriptorTest_SOURCES = UnitTestMain.cpp MoleculeDescriptorTest.cpp MoleculeDescriptorTest.hpp
    … …  
    206210PlaneUnittest_LDADD = ${ALLLIBS}
    207211
    208 ShapeUnittest_SOURCES = UnitTestMain.cpp ShapeUnittest.cpp ShapeUnittest.hpp
    209 ShapeUnittest_LDADD = ${ALLLIBS}
    210 
    211212SingletonTest_SOURCES = UnitTestMain.cpp SingletonTest.cpp SingletonTest.hpp
    212213SingletonTest_LDADD = ${ALLLIBS} $(BOOST_LIB) ${BOOST_THREAD_LIB}
    … …  
    224225Tesselation_InOutsideUnitTest_LDADD = ${ALLLIBS}
    225226
    226 TestRunner_SOURCES = TestRunnerMain.cpp $(TESTSOURCES) $(TESTHEADERS)
     227TestRunner_SOURCES = TestRunnerMain.cpp ../memoryallocator.hpp ../memoryallocator.cpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp $(TESTSOURCES) $(TESTHEADERS)
    227228TestRunner_LDADD = ${ALLLIBS}
    228229

  • src/unittests/PlaneUnittest.cpp

    r4e10f5 r192f6e  
    1111#include <cppunit/extensions/TestFactoryRegistry.h>
    1212#include <cppunit/ui/text/TestRunner.h>
    13 
    14 #include <cmath>
    1513
    1614#ifdef HAVE_TESTRUNNER

  • src/unittests/linearsystemofequationsunittest.cpp

    r4e10f5 r192f6e  
    1212#include <cppunit/extensions/TestFactoryRegistry.h>
    1313#include <cppunit/ui/text/TestRunner.h>
    14 #include <cmath>
    1514
    1615#include "linearsystemofequationsunittest.hpp"

  • src/unittests/tesselation_insideoutsideunittest.cpp

    r4e10f5 r192f6e  
    1717#include "tesselation.hpp"
    1818#include "tesselation_insideoutsideunittest.hpp"
    19 #include "verbose.hpp"
    2019
    2120#ifdef HAVE_TESTRUNNER

  • src/unittests/vectorunittest.cpp

    r4e10f5 r192f6e  
    1515#include "defs.hpp"
    1616#include "log.hpp"
     17#include "memoryusageobserver.hpp"
    1718#include "vector.hpp"
    1819#include "vector_ops.hpp"
    … …  
    2021#include "Plane.hpp"
    2122#include "Exceptions/LinearDependenceException.hpp"
    22 #include "Matrix.hpp"
    2323
    2424#ifdef HAVE_TESTRUNNER
    … …  
    214214  CPPUNIT_ASSERT(testVector.ScalarProduct(three) < MYEPSILON);
    215215}
     216
     217
     218/**
     219 * UnitTest for Vector::IsInParallelepiped().
     220 */
     221void VectorTest::IsInParallelepipedTest()
     222{
     223  double parallelepiped[NDIM*NDIM];
     224  parallelepiped[0] = 1;
     225  parallelepiped[1] = 0;
     226  parallelepiped[2] = 0;
     227  parallelepiped[3] = 0;
     228  parallelepiped[4] = 1;
     229  parallelepiped[5] = 0;
     230  parallelepiped[6] = 0;
     231  parallelepiped[7] = 0;
     232  parallelepiped[8] = 1;
     233
     234  fixture = zero;
     235  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     236  fixture = Vector(2.5,2.5,2.5);
     237  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     238  fixture = Vector(1.,1.,1.);
     239  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     240  fixture = Vector(3.5,3.5,3.5);
     241  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     242  fixture = Vector(2.,2.,2.);
     243  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     244  fixture = Vector(2.,3.,2.);
     245  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     246  fixture = Vector(-2.,2.,-1.);
     247  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     248}
     249

  • src/unittests/vectorunittest.hpp

    r4e10f5 r192f6e  
    2727    CPPUNIT_TEST ( ProjectionTest );
    2828    CPPUNIT_TEST ( NormalsTest );
     29    CPPUNIT_TEST ( IsInParallelepipedTest );
    2930    CPPUNIT_TEST_SUITE_END();
    3031
    … …  
    4445    void LineIntersectionTest();
    4546    void VectorRotationTest();
     47    void IsInParallelepipedTest();
    4648
    4749private:

  • src/vector.cpp

    r4e10f5 r192f6e  
    88
    99#include "vector.hpp"
    10 #include "VectorContent.hpp"
    1110#include "verbose.hpp"
    1211#include "World.hpp"
    1312#include "Helpers/Assert.hpp"
    1413#include "Helpers/fast_functions.hpp"
    15 #include "Exceptions/MathException.hpp"
    1614
    1715#include <iostream>
    18 #include <gsl/gsl_blas.h>
    19 
    2016
    2117using namespace std;
    … …  
    2824Vector::Vector()
    2925{
    30   content = new VectorContent();
     26  content = gsl_vector_calloc (NDIM);
    3127};
    3228
    … …  
    3733Vector::Vector(const Vector& src)
    3834{
    39   content = new VectorContent();
    40   gsl_vector_memcpy(content->content, src.content->content);
     35  content = gsl_vector_alloc(NDIM);
     36  gsl_vector_set(content,0,src[0]);
     37  gsl_vector_set(content,1,src[1]);
     38  gsl_vector_set(content,2,src[2]);
    4139}
    4240
    … …  
    4543Vector::Vector(const double x1, const double x2, const double x3)
    4644{
    47   content = new VectorContent();
    48   gsl_vector_set(content->content,0,x1);
    49   gsl_vector_set(content->content,1,x2);
    50   gsl_vector_set(content->content,2,x3);
    51 };
    52 
    53 Vector::Vector(VectorContent *_content) :
    54   content(_content)
    55 {}
     45  content = gsl_vector_alloc(NDIM);
     46  gsl_vector_set(content,0,x1);
     47  gsl_vector_set(content,1,x2);
     48  gsl_vector_set(content,2,x3);
     49};
    5650
    5751/**
    … …  
    6155  // check for self assignment
    6256  if(&src!=this){
    63     gsl_vector_memcpy(content->content, src.content->content);
     57    gsl_vector_set(content,0,src[0]);
     58    gsl_vector_set(content,1,src[1]);
     59    gsl_vector_set(content,2,src[2]);
    6460  }
    6561  return *this;
    … …  
    6965 */
    7066Vector::~Vector() {
    71   delete content;
     67  gsl_vector_free(content);
    7268};
    7369
    … …  
    9894}
    9995
     96/** Calculates distance between this and another vector in a periodic cell.
     97 * \param *y array to second vector
     98 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
     99 * \return \f$| x - y |\f$
     100 */
     101double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
     102{
     103  double res = distance(y), tmp, matrix[NDIM*NDIM];
     104    Vector Shiftedy, TranslationVector;
     105    int N[NDIM];
     106    matrix[0] = cell_size[0];
     107    matrix[1] = cell_size[1];
     108    matrix[2] = cell_size[3];
     109    matrix[3] = cell_size[1];
     110    matrix[4] = cell_size[2];
     111    matrix[5] = cell_size[4];
     112    matrix[6] = cell_size[3];
     113    matrix[7] = cell_size[4];
     114    matrix[8] = cell_size[5];
     115    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     116    for (N[0]=-1;N[0]<=1;N[0]++)
     117      for (N[1]=-1;N[1]<=1;N[1]++)
     118        for (N[2]=-1;N[2]<=1;N[2]++) {
     119          // create the translation vector
     120          TranslationVector.Zero();
     121          for (int i=NDIM;i--;)
     122            TranslationVector[i] = (double)N[i];
     123          TranslationVector.MatrixMultiplication(matrix);
     124          // add onto the original vector to compare with
     125          Shiftedy = y + TranslationVector;
     126          // get distance and compare with minimum so far
     127          tmp = distance(Shiftedy);
     128          if (tmp < res) res = tmp;
     129        }
     130    return (res);
     131};
     132
     133/** Calculates distance between this and another vector in a periodic cell.
     134 * \param *y array to second vector
     135 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
     136 * \return \f$| x - y |^2\f$
     137 */
     138double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
     139{
     140  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
     141    Vector Shiftedy, TranslationVector;
     142    int N[NDIM];
     143    matrix[0] = cell_size[0];
     144    matrix[1] = cell_size[1];
     145    matrix[2] = cell_size[3];
     146    matrix[3] = cell_size[1];
     147    matrix[4] = cell_size[2];
     148    matrix[5] = cell_size[4];
     149    matrix[6] = cell_size[3];
     150    matrix[7] = cell_size[4];
     151    matrix[8] = cell_size[5];
     152    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     153    for (N[0]=-1;N[0]<=1;N[0]++)
     154      for (N[1]=-1;N[1]<=1;N[1]++)
     155        for (N[2]=-1;N[2]<=1;N[2]++) {
     156          // create the translation vector
     157          TranslationVector.Zero();
     158          for (int i=NDIM;i--;)
     159            TranslationVector[i] = (double)N[i];
     160          TranslationVector.MatrixMultiplication(matrix);
     161          // add onto the original vector to compare with
     162          Shiftedy = y + TranslationVector;
     163          // get distance and compare with minimum so far
     164          tmp = DistanceSquared(Shiftedy);
     165          if (tmp < res) res = tmp;
     166        }
     167    return (res);
     168};
     169
     170/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
     171 * \param *out ofstream for debugging messages
     172 * Tries to translate a vector into each adjacent neighbouring cell.
     173 */
     174void Vector::KeepPeriodic(const double * const matrix)
     175{
     176  //  int N[NDIM];
     177  //  bool flag = false;
     178    //vector Shifted, TranslationVector;
     179  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
     180  //  Log() << Verbose(2) << "Vector is: ";
     181  //  Output(out);
     182  //  Log() << Verbose(0) << endl;
     183    InverseMatrixMultiplication(matrix);
     184    for(int i=NDIM;i--;) { // correct periodically
     185      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
     186        at(i) += ceil(at(i));
     187      } else {
     188        at(i) -= floor(at(i));
     189      }
     190    }
     191    MatrixMultiplication(matrix);
     192  //  Log() << Verbose(2) << "New corrected vector is: ";
     193  //  Output(out);
     194  //  Log() << Verbose(0) << endl;
     195  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
     196};
     197
    100198/** Calculates scalar product between this and another vector.
    101199 * \param *y array to second vector
    … …  
    105203{
    106204  double res = 0.;
    107   gsl_blas_ddot(content->content, y.content->content, &res);
     205  for (int i=NDIM;i--;)
     206    res += at(i)*y[i];
    108207  return (res);
    109208};
    … …  
    270369double& Vector::operator[](size_t i){
    271370  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
    272   return *gsl_vector_ptr (content->content, i);
     371  return *gsl_vector_ptr (content, i);
    273372}
    274373
    275374const double& Vector::operator[](size_t i) const{
    276375  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
    277   return *gsl_vector_ptr (content->content, i);
     376  return *gsl_vector_ptr (content, i);
    278377}
    279378
    … …  
    286385}
    287386
    288 VectorContent* Vector::get(){
     387gsl_vector* Vector::get(){
    289388  return content;
    290389}
    … …  
    405504};
    406505
    407 void Vector::ScaleAll(const Vector &factor){
    408   gsl_vector_mul(content->content, factor.content->content);
    409 }
    410506
    411507
    412508void Vector::Scale(const double factor)
    413509{
    414   gsl_vector_scale(content->content,factor);
     510  for (int i=NDIM;i--;)
     511    at(i) *= factor;
     512};
     513
     514/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
     515 * \param *M matrix of box
     516 * \param *Minv inverse matrix
     517 */
     518void Vector::WrapPeriodically(const double * const M, const double * const Minv)
     519{
     520  MatrixMultiplication(Minv);
     521  // truncate to [0,1] for each axis
     522  for (int i=0;i<NDIM;i++) {
     523    //at(i) += 0.5;  // set to center of box
     524    while (at(i) >= 1.)
     525      at(i) -= 1.;
     526    while (at(i) < 0.)
     527      at(i) += 1.;
     528  }
     529  MatrixMultiplication(M);
    415530};
    416531
    … …  
    431546  return make_pair(res,helper);
    432547}
     548
     549/** Do a matrix multiplication.
     550 * \param *matrix NDIM_NDIM array
     551 */
     552void Vector::MatrixMultiplication(const double * const M)
     553{
     554  Vector tmp;
     555  // do the matrix multiplication
     556  for(int i=NDIM;i--;)
     557    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
     558
     559  (*this) = tmp;
     560};
     561
     562/** Do a matrix multiplication with the \a *A' inverse.
     563 * \param *matrix NDIM_NDIM array
     564 */
     565bool Vector::InverseMatrixMultiplication(const double * const A)
     566{
     567  double B[NDIM*NDIM];
     568  double detA = RDET3(A);
     569  double detAReci;
     570
     571  // calculate the inverse B
     572  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
     573    detAReci = 1./detA;
     574    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
     575    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
     576    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
     577    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
     578    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
     579    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
     580    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
     581    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     582    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
     583
     584    MatrixMultiplication(B);
     585
     586    return true;
     587  } else {
     588    return false;
     589  }
     590};
     591
    433592
    434593/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
    … …  
    520679void Vector::AddVector(const Vector &y)
    521680{
    522   gsl_vector_add(content->content, y.content->content);
     681  for(int i=NDIM;i--;)
     682    at(i) += y[i];
    523683}
    524684
    … …  
    528688void Vector::SubtractVector(const Vector &y)
    529689{
    530   gsl_vector_sub(content->content, y.content->content);
     690  for(int i=NDIM;i--;)
     691    at(i) -= y[i];
     692}
     693
     694/**
     695 * Checks whether this vector is within the parallelepiped defined by the given three vectors and
     696 * their offset.
     697 *
     698 * @param offest for the origin of the parallelepiped
     699 * @param three vectors forming the matrix that defines the shape of the parallelpiped
     700 */
     701bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
     702{
     703  Vector a = (*this)-offset;
     704  a.InverseMatrixMultiplication(parallelepiped);
     705  bool isInside = true;
     706
     707  for (int i=NDIM;i--;)
     708    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
     709
     710  return isInside;
    531711}
    532712

  • src/vector.hpp

    r4e10f5 r192f6e  
    1111#endif
    1212
    13 #include <iosfwd>
     13#include <iostream>
     14#include <gsl/gsl_vector.h>
     15#include <gsl/gsl_multimin.h>
    1416
    1517#include <memory>
    … …  
    2224
    2325class Vector;
    24 class Matrix;
    25 struct VectorContent;
    2626
    2727typedef std::vector<Vector> pointset;
    … …  
    3131 */
    3232class Vector : public Space{
    33   friend Vector operator*(const Matrix&,const Vector&);
    34   friend class Matrix;
    3533public:
     34
    3635  Vector();
    3736  Vector(const double x1, const double x2, const double x3);
    … …  
    4342  double DistanceSquared(const Vector &y) const;
    4443  double DistanceToSpace(const Space& space) const;
     44  double PeriodicDistance(const Vector &y, const double * const cell_size) const;
     45  double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
    4546  double ScalarProduct(const Vector &y) const;
    4647  double Angle(const Vector &y) const;
    … …  
    5758  Vector Projection(const Vector &y) const;
    5859  void ScaleAll(const double *factor);
    59   void ScaleAll(const Vector &factor);
    6060  void Scale(const double factor);
     61  void MatrixMultiplication(const double * const M);
     62  bool InverseMatrixMultiplication(const double * const M);
     63  void KeepPeriodic(const double * const matrix);
    6164  bool GetOneNormalVector(const Vector &x1);
    6265  bool MakeNormalTo(const Vector &y1);
     66  bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
     67  void WrapPeriodically(const double * const M, const double * const Minv);
    6368  std::pair<Vector,Vector> partition(const Vector&) const;
    6469  std::pair<pointset,Vector> partition(const pointset&) const;
    … …  
    7479
    7580  // Access to internal structure
    76   VectorContent* get();
     81  gsl_vector* get();
    7782
    7883  // Methods that are derived directly from other methods
    … …  
    99104
    100105private:
    101   Vector(VectorContent *);
    102   VectorContent *content;
     106  gsl_vector *content;
    103107
    104108};

  • src/vector_ops.cpp

    r4e10f5 r192f6e  
    2323#include <gsl/gsl_permutation.h>
    2424#include <gsl/gsl_vector.h>
    25 #include <gsl/gsl_multimin.h>
    2625
    2726/**

  • src/verbose.cpp

    r4e10f5 r192f6e  
    55#include "info.hpp"
    66#include "verbose.hpp"
    7 #include <iostream>
    87
    98/** Prints the tabs according to verbosity stored in the temporary constructed class.

  • src/verbose.hpp

    r4e10f5 r192f6e  
    1818#endif
    1919
    20 #include <iosfwd>
     20#include <iostream>
    2121
    2222/************************************* Class Verbose & Binary *******************************/

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

    r4e10f5 r192f6e  
    4747BoxLength                       # (Length of a unit cell)
    48481       
    49 0       1       
     490       0       
    50500       0       2       
    5151
Note: See TracChangeset for help on using the changeset viewer.