Changes in src/boundary.cpp [bdc91e:986ed3]

File:

• src/boundary.cpp (modified) (15 diffs)

src/boundary.cpp

@@ -15 +15 @@
 #include "info.hpp"
 #include "linkedcell.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "tesselation.hpp"
@@ -22 +22 @@
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+
+#include <iostream>
+#include <iomanip>
 
 #include<gsl/gsl_poly.h>
@@ -296 +301 @@
  * \param *out output stream for debugging
  * \param *mol molecule structure with Atom's and Bond's.
- * \param *BoundaryPts set of boundary points to use or NULL
  * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
  * \param *LCList atoms in LinkedCell list
@@ -302 +306 @@
  * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
  */
-void FindConvexBorder(const molecule* mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
+void FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
 {
         Info FunctionInfo(__func__);
@@ -313 +317 @@
 
   // 1. Find all points on the boundary
-  if (BoundaryPts == NULL) {
-    BoundaryFreeFlag = true;
-    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
+  if (BoundaryPoints == NULL) {
+      BoundaryFreeFlag = true;
+      BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
   } else {
-    BoundaryPoints = BoundaryPts;
-    DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
+      DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
   }
 
 // printing all inserted for debugging
-  for (int axis=0; axis < NDIM; axis++) {
-    DoLog(1) && (Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl);
-    int i=0;
-    for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-      if (runner != BoundaryPoints[axis].begin())
-        DoLog(0) && (Log() << Verbose(0) << ", " << i << ": " << *runner->second.second);
-      else
-        DoLog(0) && (Log() << Verbose(0) << i << ": " << *runner->second.second);
-      i++;
-    }
-    DoLog(0) && (Log() << Verbose(0) << endl);
-  }
+  for (int axis=0; axis < NDIM; axis++)
+    {
+      DoLog(1) && (Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl);
+      int i=0;
+      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        if (runner != BoundaryPoints[axis].begin())
+          DoLog(0) && (Log() << Verbose(0) << ", " << i << ": " << *runner->second.second);
+        else
+          DoLog(0) && (Log() << Verbose(0) << i << ": " << *runner->second.second);
+        i++;
+      }
+      DoLog(0) && (Log() << Verbose(0) << endl);
+    }
 
   // 2. fill the boundary point list
@@ -339 +343 @@
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
         if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
-          DoLog(2) && (Log()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present." << endl);
+          DoeLog(2) && (eLog()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl);
 
   DoLog(0) && (Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl);
@@ -673 +677 @@
 
   IsAngstroem = configuration->GetIsAngstroem();
+  GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
   BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
-  GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
   LinkedCell *LCList = new LinkedCell(mol, 10.);
-  FindConvexBorder(mol, BoundaryPoints, TesselStruct, (const LinkedCell *&)LCList, NULL);
+  FindConvexBorder(mol, TesselStruct, (const LinkedCell *&)LCList, NULL);
   delete (LCList);
-  delete[] BoundaryPoints;
 
 
@@ -686 +689 @@
   else
     clustervolume = ClusterVolume;
-
-  delete TesselStruct;
 
   for (int i = 0; i < NDIM; i++)
@@ -740 +741 @@
     mol->CenterInBox();
   }
-  delete GreatestDiameter;
   // update Box of atoms by boundary
   mol->SetBoxDimension(&BoxLengths);
@@ -769 +769 @@
   int N[NDIM];
   int n[NDIM];
-  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-  double Rotations[NDIM*NDIM];
-  double *MInverse = InverseMatrix(M);
+  const Matrix &M = World::getInstance().getDomain().getM();
+  Matrix Rotations;
+  const Matrix &MInverse = World::getInstance().getDomain().getMinv();
   Vector AtomTranslations;
   Vector FillerTranslations;
@@ -804 +804 @@
 
   // calculate filler grid in [0,1]^3
-  FillerDistance = Vector(distance[0], distance[1], distance[2]);
-  FillerDistance.InverseMatrixMultiplication(M);
+  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
   for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance[i]);
@@ -818 +817 @@
       for (n[2] = 0; n[2] < N[2]; n[2]++) {
         // calculate position of current grid vector in untransformed box
-        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
-        CurrentPosition.MatrixMultiplication(M);
+        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
         // create molecule random translation vector ...
         for (int i=0;i<NDIM;i++)
@@ -840 +838 @@
             }
 
-            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
-            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
-            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
-            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
-            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
-            Rotations[7] =               sin(phi[1])                                                ;
-            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
-            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
-            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+            Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
+            Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
+            Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
+            Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
+            Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
+            Rotations.set(1,2,              sin(phi[1])                                                    );
+            Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
+            Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
+            Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
           }
 
@@ -854 +852 @@
           Inserter = (*iter)->x;
           if (DoRandomRotation)
-            Inserter.MatrixMultiplication(Rotations);
+            Inserter *= Rotations;
           Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
 
           // check whether inserter is inside box
-          Inserter.MatrixMultiplication(MInverse);
+          Inserter *= MInverse;
           FillIt = true;
           for (int i=0;i<NDIM;i++)
             FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
-          Inserter.MatrixMultiplication(M);
+          Inserter *= M;
 
           // Check whether point is in- or outside
@@ -897 +895 @@
             }
       }
-  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
-    delete LCList[*ListRunner];
-    delete TesselStruct[(*ListRunner)];
-  }
-  delete[](M);
-  delete[](MInverse);
 
   return Filling;