Changeset ef0e6d for src

Timestamp:

Aug 4, 2009, 1:53:55 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

9f97c5

Parents:

62bb91

Message:

find_(non)_convex_border TesselStruct is not free'd anymore, mol->TesselStruct used, find_non_convex_border fixes but incomplete

Location:

src

Files:

• boundary.cpp (modified) (23 diffs)
• boundary.hpp (modified) (1 diff)

src/boundary.cpp

@@ -439 +439 @@
  * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
  */
-void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
+void Find_convex_border(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename)
 {
   bool BoundaryFreeFlag = false;
@@ -446 +446 @@
   cout << Verbose(1) << "Begin of find_convex_border" << endl;
 
-  if (TesselStruct != NULL) // free if allocated
-    delete(TesselStruct);
-  TesselStruct = new class Tesselation;
+  if (mol->TesselStruct != NULL) // free if allocated
+    delete(mol->TesselStruct);
+  mol->TesselStruct = new class Tesselation;
 
   // 1. Find all points on the boundary
@@ -476 +476 @@
   for (int axis = 0; axis < NDIM; axis++)
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
-        TesselStruct->AddPoint(runner->second.second);
-
-  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
+        mol->TesselStruct->AddPoint(runner->second.second);
+
+  *out << Verbose(2) << "I found " << mol->TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
   // now we have the whole set of edge points in the BoundaryList
 
@@ -489 +489 @@
 
   // 3a. guess starting triangle
-  TesselStruct->GuessStartingTriangle(out);
+  mol->TesselStruct->GuessStartingTriangle(out);
 
   // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
-  TesselStruct->TesselateOnBoundary(out, mol);
+  mol->TesselStruct->TesselateOnBoundary(out, mol);
 
   // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
-  if (!TesselStruct->InsertStraddlingPoints(out, mol))
+  if (!mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList))
     *out << Verbose(1) << "Insertion of straddling points failed!" << endl;
 
+  *out << Verbose(2) << "I created " << mol->TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << mol->TesselStruct->LinesOnBoundary.size() << " lines and " << mol->TesselStruct->PointsOnBoundary.size() << " points." << endl;
+
+  // 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());
+      write_tecplot_file(out, tecplot, mol->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());
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
+      rasterplot->close();
+      delete(rasterplot);
+    }
+  }
+
   // 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
-  if (!TesselStruct->CorrectConcaveBaselines(out))
+  if (!mol->TesselStruct->CorrectConcaveBaselines(out))
     *out << Verbose(1) << "Correction of concave baselines failed!" << endl;
 
-  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
+  // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
+//  if (!mol->TesselStruct->CorrectConcaveTesselPoints(out))
+//    *out << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
+
+  *out << Verbose(2) << "I created " << mol->TesselStruct->TrianglesOnBoundary.size() << " triangles with " << mol->TesselStruct->LinesOnBoundary.size() << " lines and " << mol->TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -510 +538 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -518 +546 @@
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      write_raster3d_file(out, rasterplot, TesselStruct, mol);
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
     }
   }
+
 
   // free reference lists
@@ -600 +629 @@
   class Tesselation *TesselStruct = NULL;
   LinkedCell LCList(mol, 10.);
-  Find_convex_border(out, mol, TesselStruct, &LCList, NULL);
+  Find_convex_border(out, mol, &LCList, NULL);
   double clustervolume;
   if (ClusterVolume == 0)
@@ -728 +757 @@
   Vector FillerDistance;
   double FillIt = false;
-  atom *Walker = NULL, *Runner = NULL;
+  atom *Walker = NULL;
   bond *Binder = NULL;
+  int i;
+  LinkedCell *LCList[List->ListOfMolecules.size()];
+
+  *out << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;
+
+  i=0;
+  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
+    *out << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
+    LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
+    if ((*ListRunner)->TesselStruct == NULL) {
+      *out << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
+      Find_non_convex_border((ofstream *)&cout, (*ListRunner), LCList[i], NULL, 5.);
+    }
+    i++;
+  }
 
   // Center filler at origin
@@ -735 +779 @@
   filler->Center.Zero();
 
+  filler->CountAtoms(out);
+  atom * CopyAtoms[filler->AtomCount];
+  int nr = 0;
+
   // calculate filler grid in [0,1]^3
   FillerDistance.Init(distance[0], distance[1], distance[2]);
   FillerDistance.InverseMatrixMultiplication(M);
-  for(int i=0;i<NDIM;i++)
+  *out << Verbose(1) << "INFO: Grid steps are ";
+  for(int i=0;i<NDIM;i++) {
     N[i] = (int) ceil(1./FillerDistance.x[i]);
+    *out << N[i];
+    if (i != NDIM-1)
+      *out<< ", ";
+    else
+      *out << "." << endl;
+  }
 
   // go over [0,1]^3 filler grid
@@ -746 +801 @@
       for (n[2] = 0; n[2] < N[2]; n[2]++) {
         // calculate position of current grid vector in untransformed box
-        CurrentPosition.Init(n[0], n[1], n[2]);
+        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
         CurrentPosition.MatrixMultiplication(M);
+        *out << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
         // Check whether point is in- or outside
         FillIt = true;
+        i=0;
         for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
-          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
+          // get linked cell list
+          if ((*ListRunner)->TesselStruct == NULL) {
+            *out << Verbose(1) << "ERROR: TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
+            FillIt = false;
+          } else
+            FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInnerPoint(out, CurrentPosition, LCList[i++]));
         }
 
@@ -761 +823 @@
           for (int i=0;i<NDIM;i++)
             FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+          *out << Verbose(3) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
 
           // go through all atoms
+          nr=0;
           Walker = filler->start;
-          while (Walker != filler->end) {
+          while (Walker->next != filler->end) {
             Walker = Walker->next;
             // copy atom ...
-            *Runner = new atom(Walker);
+            CopyAtoms[Walker->nr] = new atom(Walker);
 
             // create atomic random translation vector ...
@@ -793 +857 @@
             // ... and put at new position
             if (DoRandomRotation)
-              Runner->x.MatrixMultiplication(Rotations);
-            Runner->x.AddVector(&AtomTranslations);
-            Runner->x.AddVector(&FillerTranslations);
-            Runner->x.AddVector(&CurrentPosition);
+              CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
+            CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
+            CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
+            CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
+
             // insert into Filling
-            Filling->AddAtom(Runner);
+            *out << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
+            Filling->AddAtom(CopyAtoms[Walker->nr]);
           }
 
           // go through all bonds and add as well
           Binder = filler->first;
-          while(Binder != filler->last) {
+          while(Binder->next != filler->last) {
             Binder = Binder->next;
+            *out << 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);
           }
         } else {
@@ -811 +879 @@
         }
       }
+  *out << Verbose(0) << "End of FillBoxWithMolecule" << endl;
+
   return Filling;
 };
@@ -823 +893 @@
  * \para RADIUS radius of the virtual sphere
  */
-void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
+void Find_non_convex_border(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename, const double RADIUS)
 {
   int N = 0;
-  bool freeTess = false;
   bool freeLC = false;
   ofstream *tempstream = NULL;
@@ -833 +902 @@
 
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
-  if (Tess == NULL) {
+  if (mol->TesselStruct == NULL) {
     *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
-    Tess = new Tesselation;
-    freeTess = true;
+    mol->TesselStruct = new Tesselation;
+  } else {
+    delete(mol->TesselStruct);
+    *out << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
+    mol->TesselStruct = new Tesselation;
   }
   LineMap::iterator baseline;
@@ -849 +921 @@
   }
 
-  Tess->Find_starting_triangle(out, RADIUS, LCList);
-
-  baseline = Tess->LinesOnBoundary.begin();
-  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
+  mol->TesselStruct->Find_starting_triangle(out, RADIUS, LCList);
+
+  baseline = mol->TesselStruct->LinesOnBoundary.begin();
+  while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (flag)) {
     if (baseline->second->TrianglesCount == 1) {
-      failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
+      failflag = mol->TesselStruct->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
@@ -861 +933 @@
       //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
       if (baseline->second->TrianglesCount != 2)
-        cout << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
+        *out << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
     }
 
     N++;
     baseline++;
-    if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
-      baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+    if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (flag)) {
+      baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
       flag = false;
     }
 
     // write temporary envelope
-    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
+    if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+      class BoundaryTriangleSet *triangle = (mol->TesselStruct->TrianglesOnBoundary.end()--)->second;
       sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
       if (DoTecplotOutput) {
@@ -881 +953 @@
         NameofTempFile.erase(npos, 1);
         NameofTempFile.append(TecplotSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
         tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
+        write_tecplot_file(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
         tempstream->close();
         tempstream->flush();
@@ -895 +967 @@
         NameofTempFile.erase(npos, 1);
         NameofTempFile.append(Raster3DSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
         tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, Tess, mol);
+        write_raster3d_file(out, tempstream, mol->TesselStruct, mol);
 //        // include the current position of the virtual sphere in the temporary raster3d file
 //        // make the circumsphere's center absolute again
@@ -930 +1002 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      write_tecplot_file(out, tecplot, Tess, mol, -1);
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, -1);
       tecplot->close();
       delete(tecplot);
@@ -938 +1010 @@
       OutputName.append(Raster3DSuffix);
       ofstream *raster = new ofstream(OutputName.c_str());
-      write_raster3d_file(out, raster, Tess, mol);
+      write_raster3d_file(out, raster, mol->TesselStruct, mol);
       raster->close();
       delete(raster);
@@ -948 +1020 @@
   cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl;
   int counter = 0;
-  for (testline = Tess->LinesOnBoundary.begin(); testline != Tess->LinesOnBoundary.end(); testline++) {
+  for (testline = mol->TesselStruct->LinesOnBoundary.begin(); testline != mol->TesselStruct->LinesOnBoundary.end(); testline++) {
     if (testline->second->TrianglesCount != 2) {
       cout << Verbose(2) << *testline->second << "\t" << testline->second->TrianglesCount << endl;
@@ -955 +1027 @@
   }
   if (counter == 0)
-    cout << Verbose(2) << "None." << endl;
-
-  // Tests the IsInnerAtom() function.
-  Vector x (0, 0, 0);
-  cout << Verbose(0) << "Point to check: " << x << endl;
-  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
-    << "for vector " << x << "." << endl;
-  TesselPoint* a = Tess->PointsOnBoundary.begin()->second->node;
-  cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
-  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerPoint(a, Tess, LCList)
-    << "for atom " << a << " (on boundary)." << endl;
-  LinkedNodes *List = NULL;
-  for (int i=0;i<NDIM;i++) { // each axis
-    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
-    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
-      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
-        List = LCList->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
-              a = *Runner;
-              i=3;
-              for (int j=0;j<NDIM;j++)
-                LCList->n[j] = LCList->N[j];
-              break;
-            }
-          }
-        }
-      }
-  }
-  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(a, Tess, LCList)
-    << "for atom " << a << " (inside)." << endl;
-
-
-  if (freeTess)
-    delete(Tess);
+    *out << Verbose(2) << "None." << endl;
+
+//  // Tests the IsInnerAtom() function.
+//  Vector x (0, 0, 0);
+//  *out << Verbose(0) << "Point to check: " << x << endl;
+//  *out << Verbose(0) << "Check: IsInnerPoint() returns " << mol->TesselStruct->IsInnerPoint(out, x, LCList)
+//    << "for vector " << x << "." << endl;
+//  TesselPoint* a = mol->TesselStruct->PointsOnBoundary.begin()->second->node;
+//  *out << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
+//  *out << Verbose(0) << "Check: IsInnerAtom() returns " << mol->TesselStruct->IsInnerPoint(out, a, LCList)
+//    << "for atom " << a << " (on boundary)." << endl;
+//  LinkedNodes *List = NULL;
+//  for (int i=0;i<NDIM;i++) { // each axis
+//    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
+//    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
+//      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
+//        List = LCList->GetCurrentCell();
+//        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+//        if (List != NULL) {
+//          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
+//            if (mol->TesselStruct->PointsOnBoundary.find((*Runner)->nr) == mol->TesselStruct->PointsOnBoundary.end()) {
+//              a = *Runner;
+//              i=3;
+//              for (int j=0;j<NDIM;j++)
+//                LCList->n[j] = LCList->N[j];
+//              break;
+//            }
+//          }
+//        }
+//      }
+//  }
+//  *out << Verbose(0) << "Check: IsInnerPoint() returns " << mol->TesselStruct->IsInnerPoint(out, a, LCList)
+//    << "for atom " << a << " (inside)." << endl;
+
+
   if (freeLC)
     delete(LCList);

src/boundary.hpp

@@ -36 +36 @@
 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
 molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
-void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename);
-void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
+void Find_convex_border(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename);
+void Find_non_convex_border(ofstream *out, molecule* mol, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
 Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);