Changes in src/boundary.cpp [8f215d:112b09]

File:

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

src/boundary.cpp

@@ -3 +3 @@
  * Implementations and super-function for envelopes
  */
+
+#include "Helpers/MemDebug.hpp"
 
 #include "World.hpp"
@@ -139 +141 @@
 {
         Info FunctionInfo(__func__);
-  atom *Walker = NULL;
   PointMap PointsOnBoundary;
   LineMap LinesOnBoundary;
@@ -165 +166 @@
 
     // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
-    Walker = mol->start;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      ProjectedVector = Walker->x - (*MolCenter);
+    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      ProjectedVector = (*iter)->x - (*MolCenter);
       ProjectedVector.ProjectOntoPlane(AxisVector);
 
@@ -182 +181 @@
         angle = 2. * M_PI - angle;
       }
-      DoLog(1) && (Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
-      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
+    DoLog(1) && (Log() << Verbose(1) << "Inserting " << **iter << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
+      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, (*iter))));
       if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
         DoLog(2) && (Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl);
         DoLog(2) && (Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl);
-        DoLog(2) && (Log() << Verbose(2) << "New vector: " << *Walker << endl);
+        DoLog(2) && (Log() << Verbose(2) << "New vector: " << **iter << endl);
         const double ProjectedVectorNorm = ProjectedVector.NormSquared();
         if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
           BoundaryTestPair.first->second.first = ProjectedVectorNorm;
-          BoundaryTestPair.first->second.second = Walker;
+          BoundaryTestPair.first->second.second = (*iter);
           DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
         } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper = Walker->x - (*MolCenter);
+          helper = (*iter)->x;
+          helper -= *MolCenter;
           const double oldhelperNorm = helper.NormSquared();
           helper = BoundaryTestPair.first->second.second->x - (*MolCenter);
           if (helper.NormSquared() < oldhelperNorm) {
-            BoundaryTestPair.first->second.second = Walker;
+            BoundaryTestPair.first->second.second = (*iter);
             DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl);
           } else {
@@ -226 +226 @@
     do { // do as long as we still throw one out per round
       flag = false;
-      Boundaries::iterator left = BoundaryPoints[axis].end();
-      Boundaries::iterator right = BoundaryPoints[axis].end();
-      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+      Boundaries::iterator left = BoundaryPoints[axis].begin();
+      Boundaries::iterator right = BoundaryPoints[axis].begin();
+      Boundaries::iterator runner = BoundaryPoints[axis].begin();
+      bool LoopOnceDone = false;
+      while (!LoopOnceDone) {
+        runner = right;
+        right++;
         // set neighbours correctly
         if (runner == BoundaryPoints[axis].begin()) {
@@ -236 +240 @@
         }
         left--;
-        right = runner;
-        right++;
         if (right == BoundaryPoints[axis].end()) {
           right = BoundaryPoints[axis].begin();
+          LoopOnceDone = true;
         }
         // check distance
@@ -279 +282 @@
             DoLog(1) && (Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl);
             BoundaryPoints[axis].erase(runner);
+            runner = right;
             flag = true;
           }
@@ -359 +363 @@
 
   // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append("_intermed");
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(tecplot, TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append("_intermed");
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(rasterplot, TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
-    }
-  }
+  StoreTrianglesinFile(mol, TesselStruct, filename, "_intermed");
 
   // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
@@ -397 +382 @@
           TesselStruct->FlipBaseline(line);
           DoLog(1) && (Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl);
+          LineRunner = TesselStruct->LinesOnBoundary.begin(); // LineRunner may have been erase if line was deleted from LinesOnBoundary
         }
       }
@@ -409 +395 @@
 
   // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(tecplot, TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(rasterplot, TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
-    }
-  }
-
+  StoreTrianglesinFile(mol, TesselStruct, filename, "");
 
   // free reference lists
@@ -678 +646 @@
 /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
  * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
+ * TODO: Here, we need a VolumeOfGeneralEnvelope (i.e. non-convex one)
  * \param *out output stream for debugging
  * \param *configuration needed for path to store convex envelope file
@@ -695 +664 @@
   int repetition[NDIM] = { 1, 1, 1 };
   int TotalNoClusters = 1;
-  atom *Walker = NULL;
   double totalmass = 0.;
   double clustervolume = 0.;
@@ -706 +674 @@
   GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
   BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
-  LinkedCell LCList(mol, 10.);
-  FindConvexBorder(mol, TesselStruct, &LCList, NULL);
+  LinkedCell *LCList = new LinkedCell(mol, 10.);
+  FindConvexBorder(mol, TesselStruct, (const LinkedCell *&)LCList, NULL);
+  delete (LCList);
+
 
   // some preparations beforehand
@@ -719 +689 @@
 
   // sum up the atomic masses
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      totalmass += Walker->type->mass;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      totalmass += (*iter)->type->mass;
   }
   DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl);
@@ -804 +772 @@
   Vector Inserter;
   double FillIt = false;
-  atom *Walker = NULL;
   bond *Binder = NULL;
   double phi[NDIM];
@@ -811 +778 @@
 
   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
-    if ((*ListRunner)->AtomCount > 0) {
+    if ((*ListRunner)->getAtomCount() > 0) {
       DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
       LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
@@ -822 +789 @@
   filler->CenterEdge(&Inserter);
   filler->Center.Zero();
+  const int FillerCount = filler->getAtomCount();
   DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
-  Binder = filler->first;
-  while(Binder->next != filler->last) {
-    Binder = Binder->next;
-    DoLog(2) && (Log() << Verbose(2) << "  " << *Binder << endl);
-  }
-
-  filler->CountAtoms();
-  atom * CopyAtoms[filler->AtomCount];
+  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
+    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+      if ((*BondRunner)->leftatom == *AtomRunner)
+        DoLog(2) && (Log() << Verbose(2) << "  " << *(*BondRunner) << endl);
+
+  atom * CopyAtoms[FillerCount];
 
   // calculate filler grid in [0,1]^3
@@ -855 +821 @@
 
         // go through all atoms
-        for (int i=0;i<filler->AtomCount;i++)
+        for (int i=0;i<FillerCount;i++)
           CopyAtoms[i] = NULL;
-        Walker = filler->start;
-        while (Walker->next != filler->end) {
-          Walker = Walker->next;
+        for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){
 
           // create atomic random translation vector ...
@@ -883 +847 @@
 
           // ... and put at new position
-          Inserter = Walker->x;
+          Inserter = (*iter)->x;
           if (DoRandomRotation)
             Inserter.MatrixMultiplication(Rotations);
@@ -907 +871 @@
             DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl);
             // copy atom ...
-            CopyAtoms[Walker->nr] = Walker->clone();
-            CopyAtoms[Walker->nr]->x = Inserter;
-            Filling->AddAtom(CopyAtoms[Walker->nr]);
-            DoLog(4) && (Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl);
+            CopyAtoms[(*iter)->nr] = (*iter)->clone();
+            CopyAtoms[(*iter)->nr]->x = Inserter;
+            Filling->AddAtom(CopyAtoms[(*iter)->nr]);
+            DoLog(4) && (Log() << Verbose(4) << "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->nr]->x) << "." << endl);
           } else {
             DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl);
-            CopyAtoms[Walker->nr] = NULL;
+            CopyAtoms[(*iter)->nr] = NULL;
             continue;
           }
         }
         // go through all bonds and add as well
-        Binder = filler->first;
-        while(Binder->next != filler->last) {
-          Binder = Binder->next;
-          if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
-            Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
-            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
-          }
-        }
+        for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
+          for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+            if ((*BondRunner)->leftatom == *AtomRunner) {
+              Binder = (*BondRunner);
+              if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
+                Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
+                Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
+              }
+            }
       }
-  Free(&M);
-  Free(&MInverse);
+  delete[](M);
+  delete[](MInverse);
 
   return Filling;
@@ -952 +917 @@
   bool TesselationFailFlag = false;
 
+  mol->getAtomCount();
+
   if (TesselStruct == NULL) {
     DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl);
@@ -1023 +990 @@
 //  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
 //  //->InsertStraddlingPoints(mol, LCList);
-//  mol->GoToFirst();
+//  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
 //  class TesselPoint *Runner = NULL;
-//  while (!mol->IsEnd()) {
-//    Runner = mol->GetPoint();
+//    Runner = *iter;
 //    Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
 //    if (!->IsInnerPoint(Runner, LCList)) {
@@ -1034 +1000 @@
 //      Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
 //    }
-//    mol->GoToNext();
 //  }
 
@@ -1043 +1008 @@
   status = CheckListOfBaselines(TesselStruct);
 
+  cout << "before correction" << endl;
+
   // store before correction
-  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
+  StoreTrianglesinFile(mol, TesselStruct, filename, "");
 
 //  // correct degenerated polygons
@@ -1054 +1021 @@
   // write final envelope
   CalculateConcavityPerBoundaryPoint(TesselStruct);
-  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
+  cout << "after correction" << endl;
+  StoreTrianglesinFile(mol, TesselStruct, filename, "");
 
   if (freeLC)