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

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>
@@ -301 +296 @@
  * \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
@@ -306 +302 @@
  * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
  */
-void FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
+void FindConvexBorder(const molecule* mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
 {
         Info FunctionInfo(__func__);
@@ -317 +313 @@
 
   // 1. Find all points on the boundary
-  if (BoundaryPoints == NULL) {
-      BoundaryFreeFlag = true;
-      BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
+  if (BoundaryPts == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
   } else {
-      DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
+    BoundaryPoints = BoundaryPts;
+    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
@@ -343 +339 @@
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
         if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
-          DoeLog(2) && (eLog()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl);
+          DoLog(2) && (Log()<< 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);
@@ -677 +673 @@
 
   IsAngstroem = configuration->GetIsAngstroem();
+  BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
   GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
-  BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
   LinkedCell *LCList = new LinkedCell(mol, 10.);
-  FindConvexBorder(mol, TesselStruct, (const LinkedCell *&)LCList, NULL);
+  FindConvexBorder(mol, BoundaryPoints, TesselStruct, (const LinkedCell *&)LCList, NULL);
   delete (LCList);
+  delete[] BoundaryPoints;
 
 
@@ -689 +686 @@
   else
     clustervolume = ClusterVolume;
+
+  delete TesselStruct;
 
   for (int i = 0; i < NDIM; i++)
@@ -741 +740 @@
     mol->CenterInBox();
   }
+  delete GreatestDiameter;
   // update Box of atoms by boundary
   mol->SetBoxDimension(&BoxLengths);
@@ -769 +769 @@
   int N[NDIM];
   int n[NDIM];
-  const Matrix &M = World::getInstance().getDomain().getM();
-  Matrix Rotations;
-  const Matrix &MInverse = World::getInstance().getDomain().getMinv();
+  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
+  double Rotations[NDIM*NDIM];
+  double *MInverse = InverseMatrix(M);
   Vector AtomTranslations;
   Vector FillerTranslations;
@@ -804 +804 @@
 
   // calculate filler grid in [0,1]^3
-  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
+  FillerDistance = Vector(distance[0], distance[1], distance[2]);
+  FillerDistance.InverseMatrixMultiplication(M);
   for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance[i]);
@@ -817 +818 @@
       for (n[2] = 0; n[2] < N[2]; n[2]++) {
         // calculate position of current grid vector in untransformed box
-        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
+        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
+        CurrentPosition.MatrixMultiplication(M);
         // create molecule random translation vector ...
         for (int i=0;i<NDIM;i++)
@@ -838 +840 @@
             }
 
-            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])                                        );
+            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])                                     ;
           }
 
@@ -852 +854 @@
           Inserter = (*iter)->x;
           if (DoRandomRotation)
-            Inserter *= Rotations;
+            Inserter.MatrixMultiplication(Rotations);
           Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
 
           // check whether inserter is inside box
-          Inserter *= MInverse;
+          Inserter.MatrixMultiplication(MInverse);
           FillIt = true;
           for (int i=0;i<NDIM;i++)
             FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
-          Inserter *= M;
+          Inserter.MatrixMultiplication(M);
 
           // Check whether point is in- or outside
@@ -895 +897 @@
             }
       }
+  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
+    delete LCList[*ListRunner];
+    delete TesselStruct[(*ListRunner)];
+  }
+  delete[](M);
+  delete[](MInverse);
 
   return Filling;