Changeset fcc9b0 for src


Ignore:
Timestamp:
Jul 27, 2009, 2:53:57 PM (16 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:
a567c3
Parents:
3c4f04
Message:

FIX: DetermineCenterOfAll() did point in wrong direction, convex envelope working again, but algorithm is not stable

  • molecule::DetermineCenterOfAll() center was the sum of all scaled by -1/Number instead of 1/Number. This was corrected for in the code elsewhere. Now, it returns the correct center and where it's elsewhere called we subtract instead of add
  • Tesselation::TesselateOnBoundary() does now something meaningful. The NormalVector for the starting triangle could not be constructed, as we lacked the initial normal vector. It is however easy to construct from the center of all and the center of the starting triangle. However, the algorithm is faulty for 1_2_dimethylethane.
  • GetBoundaryPoints() - as we do not translate to center of gravity anymore, we have to subtract the center of all. Otherwise the radial method for determining boundary points does not work.

Is not yet working.

Location:
src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    r3c4f04 rfcc9b0  
    302302  LineMap LinesOnBoundary;
    303303  TriangleMap TrianglesOnBoundary;
     304  Vector *MolCenter = mol->DetermineCenterOfAll(out);
     305  Vector helper;
    304306
    305307  *out << Verbose(1) << "Finding all boundary points." << endl;
    … …  
    309311  double radius, angle;
    310312  // 3a. Go through every axis
    311   for (int axis = 0; axis < NDIM; axis++)
    312     {
    313       AxisVector.Zero();
    314       AngleReferenceVector.Zero();
    315       AngleReferenceNormalVector.Zero();
    316       AxisVector.x[axis] = 1.;
    317       AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
    318       AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
    319       //    *out << Verbose(1) << "Axisvector is ";
    320       //    AxisVector.Output(out);
    321       //    *out << " and AngleReferenceVector is ";
    322       //    AngleReferenceVector.Output(out);
    323       //    *out << "." << endl;
    324       //    *out << " and AngleReferenceNormalVector is ";
    325       //    AngleReferenceNormalVector.Output(out);
    326       //    *out << "." << endl;
    327       // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    328       Walker = mol->start;
    329       while (Walker->next != mol->end)
     313  for (int axis = 0; axis < NDIM; axis++) {
     314    AxisVector.Zero();
     315    AngleReferenceVector.Zero();
     316    AngleReferenceNormalVector.Zero();
     317    AxisVector.x[axis] = 1.;
     318    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
     319    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
     320
     321    //*out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
     322
     323    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
     324    Walker = mol->start;
     325    while (Walker->next != mol->end) {
     326      Walker = Walker->next;
     327      Vector ProjectedVector;
     328      ProjectedVector.CopyVector(&Walker->x);
     329      ProjectedVector.SubtractVector(MolCenter);
     330      ProjectedVector.ProjectOntoPlane(&AxisVector);
     331
     332      // correct for negative side
     333      radius = ProjectedVector.Norm();
     334      if (fabs(radius) > MYEPSILON)
     335        angle = ProjectedVector.Angle(&AngleReferenceVector);
     336      else
     337        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
     338
     339      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
     340      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
     341        angle = 2. * M_PI - angle;
     342      }
     343      //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
     344      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
     345          DistancePair (radius, Walker)));
     346      if (BoundaryTestPair.second) { // successfully inserted
     347      } else { // same point exists, check first r, then distance of original vectors to center of gravity
     348        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
     349        *out << Verbose(2) << "Present vector: " << BoundaryTestPair.first->second.second->x << endl;
     350        *out << Verbose(2) << "New vector: " << Walker << endl;
     351        double tmp = ProjectedVector.Norm();
     352        if (tmp > BoundaryTestPair.first->second.first) {
     353          BoundaryTestPair.first->second.first = tmp;
     354          BoundaryTestPair.first->second.second = Walker;
     355          *out << Verbose(2) << "Keeping new vector." << endl;
     356        } else if (tmp == BoundaryTestPair.first->second.first) {
     357          helper.CopyVector(&Walker->x);
     358          helper.SubtractVector(MolCenter);
     359          if (BoundaryTestPair.first->second.second->x.ScalarProduct(&BoundaryTestPair.first->second.second->x) < helper.ScalarProduct(&helper)) { // Norm() does a sqrt, which makes it a lot slower
     360            BoundaryTestPair.first->second.second = Walker;
     361            *out << Verbose(2) << "Keeping new vector." << endl;
     362          } else {
     363            *out << Verbose(2) << "Keeping present vector." << endl;
     364          }
     365        } else {
     366          *out << Verbose(2) << "Keeping present vector." << endl;
     367        }
     368      }
     369    }
     370    // printing all inserted for debugging
     371    //    {
     372    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     373    //      int i=0;
     374    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
     375    //        if (runner != BoundaryPoints[axis].begin())
     376    //          *out << ", " << i << ": " << *runner->second.second;
     377    //        else
     378    //          *out << i << ": " << *runner->second.second;
     379    //        i++;
     380    //      }
     381    //      *out << endl;
     382    //    }
     383    // 3c. throw out points whose distance is less than the mean of left and right neighbours
     384    bool flag = false;
     385    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
     386    do { // do as long as we still throw one out per round
     387      flag = false;
     388      Boundaries::iterator left = BoundaryPoints[axis].end();
     389      Boundaries::iterator right = BoundaryPoints[axis].end();
     390      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
     391        // set neighbours correctly
     392        if (runner == BoundaryPoints[axis].begin()) {
     393          left = BoundaryPoints[axis].end();
     394        } else {
     395          left = runner;
     396        }
     397        left--;
     398        right = runner;
     399        right++;
     400        if (right == BoundaryPoints[axis].end()) {
     401          right = BoundaryPoints[axis].begin();
     402        }
     403        // check distance
     404
     405        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
    330406        {
    331           Walker = Walker->next;
    332           Vector ProjectedVector;
    333           ProjectedVector.CopyVector(&Walker->x);
    334           ProjectedVector.ProjectOntoPlane(&AxisVector);
    335           // correct for negative side
    336           //if (Projection(y) < 0)
    337           //angle = 2.*M_PI - angle;
    338           radius = ProjectedVector.Norm();
    339           if (fabs(radius) > MYEPSILON)
    340             angle = ProjectedVector.Angle(&AngleReferenceVector);
    341           else
    342             angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
    343 
    344           //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
    345           if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
    346             {
    347               angle = 2. * M_PI - angle;
    348             }
    349           //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
    350           //ProjectedVector.Output(out);
    351           //*out << endl;
    352           BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
    353               DistancePair (radius, Walker)));
    354           if (BoundaryTestPair.second)
    355             { // successfully inserted
    356             }
    357           else
    358             { // same point exists, check first r, then distance of original vectors to center of gravity
    359               *out << Verbose(2)
    360                   << "Encountered two vectors whose projection onto axis "
    361                   << axis << " is equal: " << endl;
    362               *out << Verbose(2) << "Present vector: ";
    363               BoundaryTestPair.first->second.second->x.Output(out);
    364               *out << endl;
    365               *out << Verbose(2) << "New vector: ";
    366               Walker->x.Output(out);
    367               *out << endl;
    368               double tmp = ProjectedVector.Norm();
    369               if (tmp > BoundaryTestPair.first->second.first)
    370                 {
    371                   BoundaryTestPair.first->second.first = tmp;
    372                   BoundaryTestPair.first->second.second = Walker;
    373                   *out << Verbose(2) << "Keeping new vector." << endl;
    374                 }
    375               else if (tmp == BoundaryTestPair.first->second.first)
    376                 {
    377                   if (BoundaryTestPair.first->second.second->x.ScalarProduct(
    378                       &BoundaryTestPair.first->second.second->x)
    379                       < Walker->x.ScalarProduct(&Walker->x))
    380                     { // Norm() does a sqrt, which makes it a lot slower
    381                       BoundaryTestPair.first->second.second = Walker;
    382                       *out << Verbose(2) << "Keeping new vector." << endl;
    383                     }
    384                   else
    385                     {
    386                       *out << Verbose(2) << "Keeping present vector." << endl;
    387                     }
    388                 }
    389               else
    390                 {
    391                   *out << Verbose(2) << "Keeping present vector." << endl;
    392                 }
    393             }
     407          Vector SideA, SideB, SideC, SideH;
     408          SideA.CopyVector(&left->second.second->x);
     409          SideA.SubtractVector(MolCenter);
     410          SideA.ProjectOntoPlane(&AxisVector);
     411          //          *out << "SideA: ";
     412          //          SideA.Output(out);
     413          //          *out << endl;
     414
     415          SideB.CopyVector(&right->second.second->x);
     416          SideB.SubtractVector(MolCenter);
     417          SideB.ProjectOntoPlane(&AxisVector);
     418          //          *out << "SideB: ";
     419          //          SideB.Output(out);
     420          //          *out << endl;
     421
     422          SideC.CopyVector(&left->second.second->x);
     423          SideC.SubtractVector(&right->second.second->x);
     424          SideC.ProjectOntoPlane(&AxisVector);
     425          //          *out << "SideC: ";
     426          //          SideC.Output(out);
     427          //          *out << endl;
     428
     429          SideH.CopyVector(&runner->second.second->x);
     430          SideH.SubtractVector(MolCenter);
     431          SideH.ProjectOntoPlane(&AxisVector);
     432          //          *out << "SideH: ";
     433          //          SideH.Output(out);
     434          //          *out << endl;
     435
     436          // calculate each length
     437          double a = SideA.Norm();
     438          //double b = SideB.Norm();
     439          //double c = SideC.Norm();
     440          double h = SideH.Norm();
     441          // calculate the angles
     442          double alpha = SideA.Angle(&SideH);
     443          double beta = SideA.Angle(&SideC);
     444          double gamma = SideB.Angle(&SideH);
     445          double delta = SideC.Angle(&SideH);
     446          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
     447          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
     448          //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
     449          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && (h < MinDistance)) {
     450            // throw out point
     451            //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
     452            BoundaryPoints[axis].erase(runner);
     453            flag = true;
     454          }
    394455        }
    395       // printing all inserted for debugging
    396       //    {
    397       //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
    398       //      int i=0;
    399       //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    400       //        if (runner != BoundaryPoints[axis].begin())
    401       //          *out << ", " << i << ": " << *runner->second.second;
    402       //        else
    403       //          *out << i << ": " << *runner->second.second;
    404       //        i++;
    405       //      }
    406       //      *out << endl;
    407       //    }
    408       // 3c. throw out points whose distance is less than the mean of left and right neighbours
    409       bool flag = false;
    410       do
    411         { // do as long as we still throw one out per round
    412           *out << Verbose(1)
    413               << "Looking for candidates to kick out by convex condition ... "
    414               << endl;
    415           flag = false;
    416           Boundaries::iterator left = BoundaryPoints[axis].end();
    417           Boundaries::iterator right = BoundaryPoints[axis].end();
    418           for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
    419               != BoundaryPoints[axis].end(); runner++)
    420             {
    421               // set neighbours correctly
    422               if (runner == BoundaryPoints[axis].begin())
    423                 {
    424                   left = BoundaryPoints[axis].end();
    425                 }
    426               else
    427                 {
    428                   left = runner;
    429                 }
    430               left--;
    431               right = runner;
    432               right++;
    433               if (right == BoundaryPoints[axis].end())
    434                 {
    435                   right = BoundaryPoints[axis].begin();
    436                 }
    437               // check distance
    438 
    439               // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
    440                 {
    441                   Vector SideA, SideB, SideC, SideH;
    442                   SideA.CopyVector(&left->second.second->x);
    443                   SideA.ProjectOntoPlane(&AxisVector);
    444                   //          *out << "SideA: ";
    445                   //          SideA.Output(out);
    446                   //          *out << endl;
    447 
    448                   SideB.CopyVector(&right->second.second->x);
    449                   SideB.ProjectOntoPlane(&AxisVector);
    450                   //          *out << "SideB: ";
    451                   //          SideB.Output(out);
    452                   //          *out << endl;
    453 
    454                   SideC.CopyVector(&left->second.second->x);
    455                   SideC.SubtractVector(&right->second.second->x);
    456                   SideC.ProjectOntoPlane(&AxisVector);
    457                   //          *out << "SideC: ";
    458                   //          SideC.Output(out);
    459                   //          *out << endl;
    460 
    461                   SideH.CopyVector(&runner->second.second->x);
    462                   SideH.ProjectOntoPlane(&AxisVector);
    463                   //          *out << "SideH: ";
    464                   //          SideH.Output(out);
    465                   //          *out << endl;
    466 
    467                   // calculate each length
    468                   double a = SideA.Norm();
    469                   //double b = SideB.Norm();
    470                   //double c = SideC.Norm();
    471                   double h = SideH.Norm();
    472                   // calculate the angles
    473                   double alpha = SideA.Angle(&SideH);
    474                   double beta = SideA.Angle(&SideC);
    475                   double gamma = SideB.Angle(&SideH);
    476                   double delta = SideC.Angle(&SideH);
    477                   double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
    478                       < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
    479                   //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
    480                   //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
    481                   if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
    482                       < MYEPSILON) && (h < MinDistance))
    483                     {
    484                       // throw out point
    485                       //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
    486                       BoundaryPoints[axis].erase(runner);
    487                       flag = true;
    488                     }
    489                 }
    490             }
    491         }
    492       while (flag);
    493     }
     456      }
     457    } while (flag);
     458  }
     459  delete(MolCenter);
    494460  return BoundaryPoints;
    495461}
    … …  
    620586      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
    621587      for (i=0;i<NDIM;i++)
    622         *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
     588        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
    623589      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    624590    }
    … …  
    629595      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
    630596      for (i=0;i<NDIM;i++)
    631         *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
     597        *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
    632598      *vrmlfile << "\t0.03\t";
    633599      for (i=0;i<NDIM;i++)
    634         *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
     600        *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
    635601      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
    636602    }
    … …  
    641607      for (i=0;i<3;i++) { // print each node
    642608        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    643           *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
     609          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
    644610        *vrmlfile << "\t";
    645611      }
    … …  
    674640      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
    675641      for (i=0;i<NDIM;i++)
    676         *rasterfile << Walker->x.x[i]+center->x[i] << " ";
     642        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
    677643      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    678644    }
    … …  
    683649      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
    684650      for (i=0;i<NDIM;i++)
    685         *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
     651        *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
    686652      *rasterfile << "\t0.03\t";
    687653      for (i=0;i<NDIM;i++)
    688         *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
     654        *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
    689655      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
    690656    }
    … …  
    696662      for (i=0;i<3;i++) {  // print each node
    697663        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    698           *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
     664          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
    699665        *rasterfile << "\t";
    700666      }
    … …  
    714680 * \param N arbitrary number to differentiate various zones in the tecplot format
    715681 */
    716 void
    717 write_tecplot_file(ofstream *out, ofstream *tecplot,
    718     class Tesselation *TesselStruct, class molecule *mol, int N)
     682void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
    719683{
    720684  if (tecplot != NULL)
    … …  
    770734void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
    771735{
    772   atom *Walker = NULL;
    773736  bool BoundaryFreeFlag = false;
    774737  Boundaries *BoundaryPoints = NULL;
     738
     739  cout << Verbose(1) << "Begin of find_convex_border" << endl;
    775740
    776741  if (TesselStruct != NULL) // free if allocated
    … …  
    778743  TesselStruct = new class Tesselation;
    779744
    780   // 1. calculate center of gravity
    781   *out << endl;
    782   Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
    783 
    784   // 2. translate all points into CoG
    785   *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
    786   Walker = mol->start;
    787   while (Walker->next != mol->end) {
    788     Walker = Walker->next;
    789     Walker->x.Translate(CenterOfGravity);
    790   }
    791 
    792   // 3. Find all points on the boundary
     745  // 1. Find all points on the boundary
    793746  if (BoundaryPoints == NULL) {
    794747      BoundaryFreeFlag = true;
    … …  
    798751  }
    799752
    800   // 4. fill the boundary point list
     753// printing all inserted for debugging
     754  for (int axis=0; axis < NDIM; axis++)
     755    {
     756      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     757      int i=0;
     758      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
     759        if (runner != BoundaryPoints[axis].begin())
     760          *out << ", " << i << ": " << *runner->second.second;
     761        else
     762          *out << i << ": " << *runner->second.second;
     763        i++;
     764      }
     765      *out << endl;
     766    }
     767
     768  // 2. fill the boundary point list
    801769  for (int axis = 0; axis < NDIM; axis++)
    802770    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
    803771        TesselStruct->AddPoint(runner->second.second);
    804772
    805   *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount
    806       << " points on the convex boundary." << endl;
     773  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
    807774  // now we have the whole set of edge points in the BoundaryList
    808775
    … …  
    814781  //  *out << endl;
    815782
    816   // 5a. guess starting triangle
     783  // 3a. guess starting triangle
    817784  TesselStruct->GuessStartingTriangle(out);
    818785
    819   // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
     786  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
    820787  TesselStruct->TesselateOnBoundary(out, mol);
    821788
    822   *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
    823       << " triangles with " << TesselStruct->LinesOnBoundaryCount
    824       << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
    825       << endl;
    826 
    827   // 6. translate all points back from CoG
    828   *out << Verbose(1) << "Translating system back from Center of Gravity."
    829       << endl;
    830   CenterOfGravity->Scale(-1);
    831   Walker = mol->start;
    832   while (Walker->next != mol->end)
    833     {
    834       Walker = Walker->next;
    835       Walker->x.Translate(CenterOfGravity);
    836     }
    837 
    838   // 7. Store triangles in tecplot file
     789  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
     790
     791  // 4. Store triangles in tecplot file
    839792  if (filename != NULL) {
    840793    string OutputName(filename);
    … …  
    854807    delete[] (BoundaryPoints);
    855808
     809  cout << Verbose(1) << "End of find_convex_border" << endl;
    856810};
    857811
    … …  
    864818 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
    865819 */
    866 double
    867 VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
     820double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
    868821{
    869822  bool IsAngstroem = configuration->GetIsAngstroem();
    … …  
    12401193 * \param *mol the cluster as a molecule structure
    12411194 */
    1242 void
    1243 Tesselation::TesselateOnBoundary(ofstream *out, molecule *mol)
     1195void Tesselation::TesselateOnBoundary(ofstream *out, molecule *mol)
    12441196{
    12451197  bool flag;
    … …  
    12471199  class BoundaryPointSet *peak = NULL;
    12481200  double SmallestAngle, TempAngle;
    1249   Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
    1250       PropagationVector;
     1201  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, PropagationVector, *MolCenter = NULL;
    12511202  LineMap::iterator LineChecker[2];
    1252   do
    1253     {
    1254       flag = false;
    1255       for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
    1256           != LinesOnBoundary.end(); baseline++)
    1257         if (baseline->second->TrianglesCount == 1)
    1258           {
    1259             *out << Verbose(2) << "Current baseline is between "
    1260                 << *(baseline->second) << "." << endl;
    1261             // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
    1262             SmallestAngle = M_PI;
    1263             BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
    1264             // get peak point with respect to this base line's only triangle
    1265             for (int i = 0; i < 3; i++)
    1266               if ((BTS->endpoints[i] != baseline->second->endpoints[0])
    1267                   && (BTS->endpoints[i] != baseline->second->endpoints[1]))
    1268                 peak = BTS->endpoints[i];
    1269             *out << Verbose(3) << " and has peak " << *peak << "." << endl;
    1270             // normal vector of triangle
    1271             BTS->GetNormalVector(NormalVector);
    1272             *out << Verbose(4) << "NormalVector of base triangle is ";
    1273             NormalVector.Output(out);
    1274             *out << endl;
    1275             // offset to center of triangle
    1276             CenterVector.Zero();
    1277             for (int i = 0; i < 3; i++)
    1278               CenterVector.AddVector(&BTS->endpoints[i]->node->x);
    1279             CenterVector.Scale(1. / 3.);
    1280             *out << Verbose(4) << "CenterVector of base triangle is ";
    1281             CenterVector.Output(out);
    1282             *out << endl;
    1283             // vector in propagation direction (out of triangle)
    1284             // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
    1285             TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
    1286             TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
    1287             PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
    1288             TempVector.CopyVector(&CenterVector);
    1289             TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
    1290             //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
    1291             if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
    1292               PropagationVector.Scale(-1.);
    1293             *out << Verbose(4) << "PropagationVector of base triangle is ";
    1294             PropagationVector.Output(out);
    1295             *out << endl;
    1296             winner = PointsOnBoundary.end();
    1297             for (PointMap::iterator target = PointsOnBoundary.begin(); target
    1298                 != PointsOnBoundary.end(); target++)
    1299               if ((target->second != baseline->second->endpoints[0])
    1300                   && (target->second != baseline->second->endpoints[1]))
    1301                 { // don't take the same endpoints
    1302                   *out << Verbose(3) << "Target point is " << *(target->second)
    1303                       << ":";
    1304                   bool continueflag = true;
    1305 
    1306                   VirtualNormalVector.CopyVector(
    1307                       &baseline->second->endpoints[0]->node->x);
    1308                   VirtualNormalVector.AddVector(
    1309                       &baseline->second->endpoints[0]->node->x);
    1310                   VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
    1311                   VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
    1312                   TempAngle = VirtualNormalVector.Angle(&PropagationVector);
    1313                   continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
    1314                   if (!continueflag)
    1315                     {
    1316                       *out << Verbose(4)
    1317                           << "Angle between propagation direction and base line to "
    1318                           << *(target->second) << " is " << TempAngle
    1319                           << ", bad direction!" << endl;
    1320                       continue;
    1321                     }
    1322                   else
    1323                     *out << Verbose(4)
    1324                         << "Angle between propagation direction and base line to "
    1325                         << *(target->second) << " is " << TempAngle
    1326                         << ", good direction!" << endl;
    1327                   LineChecker[0] = baseline->second->endpoints[0]->lines.find(
    1328                       target->first);
    1329                   LineChecker[1] = baseline->second->endpoints[1]->lines.find(
    1330                       target->first);
    1331                   //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
    1332                   //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
    1333                   //            else
    1334                   //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
    1335                   //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
    1336                   //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
    1337                   //            else
    1338                   //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
    1339                   // check first endpoint (if any connecting line goes to target or at least not more than 1)
    1340                   continueflag = continueflag && (((LineChecker[0]
    1341                       == baseline->second->endpoints[0]->lines.end())
    1342                       || (LineChecker[0]->second->TrianglesCount == 1)));
    1343                   if (!continueflag)
    1344                     {
    1345                       *out << Verbose(4) << *(baseline->second->endpoints[0])
    1346                           << " has line " << *(LineChecker[0]->second)
    1347                           << " to " << *(target->second)
    1348                           << " as endpoint with "
    1349                           << LineChecker[0]->second->TrianglesCount
    1350                           << " triangles." << endl;
    1351                       continue;
    1352                     }
    1353                   // check second endpoint (if any connecting line goes to target or at least not more than 1)
    1354                   continueflag = continueflag && (((LineChecker[1]
    1355                       == baseline->second->endpoints[1]->lines.end())
    1356                       || (LineChecker[1]->second->TrianglesCount == 1)));
    1357                   if (!continueflag)
    1358                     {
    1359                       *out << Verbose(4) << *(baseline->second->endpoints[1])
    1360                           << " has line " << *(LineChecker[1]->second)
    1361                           << " to " << *(target->second)
    1362                           << " as endpoint with "
    1363                           << LineChecker[1]->second->TrianglesCount
    1364                           << " triangles." << endl;
    1365                       continue;
    1366                     }
    1367                   // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
    1368                   continueflag = continueflag && (!(((LineChecker[0]
    1369                       != baseline->second->endpoints[0]->lines.end())
    1370                       && (LineChecker[1]
    1371                           != baseline->second->endpoints[1]->lines.end())
    1372                       && (GetCommonEndpoint(LineChecker[0]->second,
    1373                           LineChecker[1]->second) == peak))));
    1374                   if (!continueflag)
    1375                     {
    1376                       *out << Verbose(4) << "Current target is peak!" << endl;
    1377                       continue;
    1378                     }
    1379                   // in case NOT both were found
    1380                   if (continueflag)
    1381                     { // create virtually this triangle, get its normal vector, calculate angle
    1382                       flag = true;
    1383                       VirtualNormalVector.MakeNormalVector(
    1384                           &baseline->second->endpoints[0]->node->x,
    1385                           &baseline->second->endpoints[1]->node->x,
    1386                           &target->second->node->x);
    1387                       // make it always point inward
    1388                       if (baseline->second->endpoints[0]->node->x.Projection(
    1389                           &VirtualNormalVector) > 0)
    1390                         VirtualNormalVector.Scale(-1.);
    1391                       // calculate angle
    1392                       TempAngle = NormalVector.Angle(&VirtualNormalVector);
    1393                       *out << Verbose(4) << "NormalVector is ";
    1394                       VirtualNormalVector.Output(out);
    1395                       *out << " and the angle is " << TempAngle << "." << endl;
    1396                       if (SmallestAngle > TempAngle)
    1397                         { // set to new possible winner
    1398                           SmallestAngle = TempAngle;
    1399                           winner = target;
    1400                         }
     1203
     1204  MolCenter = mol->DetermineCenterOfAll(out);
     1205  do {
     1206    flag = false;
     1207    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
     1208      if (baseline->second->TrianglesCount == 1) {
     1209        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
     1210        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
     1211        SmallestAngle = M_PI;
     1212        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
     1213        // get peak point with respect to this base line's only triangle
     1214        for (int i = 0; i < 3; i++)
     1215          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
     1216            peak = BTS->endpoints[i];
     1217        *out << Verbose(3) << " and has peak " << *peak << "." << endl;
     1218        // offset to center of triangle
     1219        CenterVector.Zero();
     1220        for (int i = 0; i < 3; i++)
     1221          CenterVector.AddVector(&BTS->endpoints[i]->node->x);
     1222        CenterVector.Scale(1. / 3.);
     1223        *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;
     1224        NormalVector.CopyVector(MolCenter);
     1225        NormalVector.SubtractVector(&CenterVector);
     1226        // normal vector of triangle
     1227        BTS->GetNormalVector(NormalVector);
     1228        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
     1229        // vector in propagation direction (out of triangle)
     1230        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
     1231        TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
     1232        TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
     1233        PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
     1234        TempVector.CopyVector(&CenterVector);
     1235        TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
     1236        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
     1237        if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
     1238          PropagationVector.Scale(-1.);
     1239        *out << Verbose(4) << "PropagationVector of base triangle is ";
     1240        PropagationVector.Output(out);
     1241        *out << endl;
     1242        winner = PointsOnBoundary.end();
     1243        for (PointMap::iterator target = PointsOnBoundary.begin(); target
     1244            != PointsOnBoundary.end(); target++)
     1245          if ((target->second != baseline->second->endpoints[0])
     1246              && (target->second != baseline->second->endpoints[1]))
     1247            { // don't take the same endpoints
     1248              *out << Verbose(3) << "Target point is " << *(target->second)
     1249                  << ":";
     1250              bool continueflag = true;
     1251
     1252              VirtualNormalVector.CopyVector(
     1253                  &baseline->second->endpoints[0]->node->x);
     1254              VirtualNormalVector.AddVector(
     1255                  &baseline->second->endpoints[0]->node->x);
     1256              VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
     1257              VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
     1258              TempAngle = VirtualNormalVector.Angle(&PropagationVector);
     1259              continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
     1260              if (!continueflag)
     1261                {
     1262                  *out << Verbose(4)
     1263                      << "Angle between propagation direction and base line to "
     1264                      << *(target->second) << " is " << TempAngle
     1265                      << ", bad direction!" << endl;
     1266                  continue;
     1267                }
     1268              else
     1269                *out << Verbose(4)
     1270                    << "Angle between propagation direction and base line to "
     1271                    << *(target->second) << " is " << TempAngle
     1272                    << ", good direction!" << endl;
     1273              LineChecker[0] = baseline->second->endpoints[0]->lines.find(
     1274                  target->first);
     1275              LineChecker[1] = baseline->second->endpoints[1]->lines.find(
     1276                  target->first);
     1277              //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
     1278              //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
     1279              //            else
     1280              //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
     1281              //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
     1282              //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
     1283              //            else
     1284              //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
     1285              // check first endpoint (if any connecting line goes to target or at least not more than 1)
     1286              continueflag = continueflag && (((LineChecker[0]
     1287                  == baseline->second->endpoints[0]->lines.end())
     1288                  || (LineChecker[0]->second->TrianglesCount == 1)));
     1289              if (!continueflag)
     1290                {
     1291                  *out << Verbose(4) << *(baseline->second->endpoints[0])
     1292                      << " has line " << *(LineChecker[0]->second)
     1293                      << " to " << *(target->second)
     1294                      << " as endpoint with "
     1295                      << LineChecker[0]->second->TrianglesCount
     1296                      << " triangles." << endl;
     1297                  continue;
     1298                }
     1299              // check second endpoint (if any connecting line goes to target or at least not more than 1)
     1300              continueflag = continueflag && (((LineChecker[1]
     1301                  == baseline->second->endpoints[1]->lines.end())
     1302                  || (LineChecker[1]->second->TrianglesCount == 1)));
     1303              if (!continueflag)
     1304                {
     1305                  *out << Verbose(4) << *(baseline->second->endpoints[1])
     1306                      << " has line " << *(LineChecker[1]->second)
     1307                      << " to " << *(target->second)
     1308                      << " as endpoint with "
     1309                      << LineChecker[1]->second->TrianglesCount
     1310                      << " triangles." << endl;
     1311                  continue;
     1312                }
     1313              // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
     1314              continueflag = continueflag && (!(((LineChecker[0]
     1315                  != baseline->second->endpoints[0]->lines.end())
     1316                  && (LineChecker[1]
     1317                      != baseline->second->endpoints[1]->lines.end())
     1318                  && (GetCommonEndpoint(LineChecker[0]->second,
     1319                      LineChecker[1]->second) == peak))));
     1320              if (!continueflag)
     1321                {
     1322                  *out << Verbose(4) << "Current target is peak!" << endl;
     1323                  continue;
     1324                }
     1325              // in case NOT both were found
     1326              if (continueflag)
     1327                { // create virtually this triangle, get its normal vector, calculate angle
     1328                  flag = true;
     1329                  VirtualNormalVector.MakeNormalVector(
     1330                      &baseline->second->endpoints[0]->node->x,
     1331                      &baseline->second->endpoints[1]->node->x,
     1332                      &target->second->node->x);
     1333                  // make it always point inward
     1334                  if (baseline->second->endpoints[0]->node->x.Projection(
     1335                      &VirtualNormalVector) > 0)
     1336                    VirtualNormalVector.Scale(-1.);
     1337                  // calculate angle
     1338                  TempAngle = NormalVector.Angle(&VirtualNormalVector);
     1339                  *out << Verbose(4) << "NormalVector is ";
     1340                  VirtualNormalVector.Output(out);
     1341                  *out << " and the angle is " << TempAngle << "." << endl;
     1342                  if (SmallestAngle > TempAngle)
     1343                    { // set to new possible winner
     1344                      SmallestAngle = TempAngle;
     1345                      winner = target;
    14011346                    }
    14021347                }
    1403             // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
    1404             if (winner != PointsOnBoundary.end())
    1405               {
    1406                 *out << Verbose(2) << "Winning target point is "
    1407                     << *(winner->second) << " with angle " << SmallestAngle
    1408                     << "." << endl;
    1409                 // create the lins of not yet present
    1410                 BLS[0] = baseline->second;
    1411                 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
    1412                 LineChecker[0] = baseline->second->endpoints[0]->lines.find(
    1413                     winner->first);
    1414                 LineChecker[1] = baseline->second->endpoints[1]->lines.find(
    1415                     winner->first);
    1416                 if (LineChecker[0]
    1417                     == baseline->second->endpoints[0]->lines.end())
    1418                   { // create
    1419                     BPS[0] = baseline->second->endpoints[0];
    1420                     BPS[1] = winner->second;
    1421                     BLS[1] = new class BoundaryLineSet(BPS,
    1422                         LinesOnBoundaryCount);
    1423                     LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
    1424                         BLS[1]));
    1425                     LinesOnBoundaryCount++;
    1426                   }
    1427                 else
    1428                   BLS[1] = LineChecker[0]->second;
    1429                 if (LineChecker[1]
    1430                     == baseline->second->endpoints[1]->lines.end())
    1431                   { // create
    1432                     BPS[0] = baseline->second->endpoints[1];
    1433                     BPS[1] = winner->second;
    1434                     BLS[2] = new class BoundaryLineSet(BPS,
    1435                         LinesOnBoundaryCount);
    1436                     LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
    1437                         BLS[2]));
    1438                     LinesOnBoundaryCount++;
    1439                   }
    1440                 else
    1441                   BLS[2] = LineChecker[1]->second;
    1442                 BTS = new class BoundaryTriangleSet(BLS,
    1443                     TrianglesOnBoundaryCount);
    1444                 TrianglesOnBoundary.insert(TrianglePair(
    1445                     TrianglesOnBoundaryCount, BTS));
    1446                 TrianglesOnBoundaryCount++;
     1348            }
     1349        // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
     1350        if (winner != PointsOnBoundary.end())
     1351          {
     1352            *out << Verbose(2) << "Winning target point is "
     1353                << *(winner->second) << " with angle " << SmallestAngle
     1354                << "." << endl;
     1355            // create the lins of not yet present
     1356            BLS[0] = baseline->second;
     1357            // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
     1358            LineChecker[0] = baseline->second->endpoints[0]->lines.find(
     1359                winner->first);
     1360            LineChecker[1] = baseline->second->endpoints[1]->lines.find(
     1361                winner->first);
     1362            if (LineChecker[0]
     1363                == baseline->second->endpoints[0]->lines.end())
     1364              { // create
     1365                BPS[0] = baseline->second->endpoints[0];
     1366                BPS[1] = winner->second;
     1367                BLS[1] = new class BoundaryLineSet(BPS,
     1368                    LinesOnBoundaryCount);
     1369                LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
     1370                    BLS[1]));
     1371                LinesOnBoundaryCount++;
    14471372              }
    14481373            else
    1449               {
    1450                 *out << Verbose(1)
    1451                     << "I could not determine a winner for this baseline "
    1452                     << *(baseline->second) << "." << endl;
     1374              BLS[1] = LineChecker[0]->second;
     1375            if (LineChecker[1]
     1376                == baseline->second->endpoints[1]->lines.end())
     1377              { // create
     1378                BPS[0] = baseline->second->endpoints[1];
     1379                BPS[1] = winner->second;
     1380                BLS[2] = new class BoundaryLineSet(BPS,
     1381                    LinesOnBoundaryCount);
     1382                LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
     1383                    BLS[2]));
     1384                LinesOnBoundaryCount++;
    14531385              }
    1454 
    1455             // 5d. If the set of lines is not yet empty, go to 5. and continue
     1386            else
     1387              BLS[2] = LineChecker[1]->second;
     1388            BTS = new class BoundaryTriangleSet(BLS,
     1389                TrianglesOnBoundaryCount);
     1390            TrianglesOnBoundary.insert(TrianglePair(
     1391                TrianglesOnBoundaryCount, BTS));
     1392            TrianglesOnBoundaryCount++;
    14561393          }
    14571394        else
    1458           *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
    1459               << " has a triangle count of "
    1460               << baseline->second->TrianglesCount << "." << endl;
    1461     }
    1462   while (flag);
    1463 
    1464 }
    1465 ;
     1395          {
     1396            *out << Verbose(1)
     1397                << "I could not determine a winner for this baseline "
     1398                << *(baseline->second) << "." << endl;
     1399          }
     1400
     1401        // 5d. If the set of lines is not yet empty, go to 5. and continue
     1402      }
     1403      else
     1404        *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
     1405            << " has a triangle count of "
     1406            << baseline->second->TrianglesCount << "." << endl;
     1407  } while (flag);
     1408  delete(MolCenter);
     1409};
    14661410
    14671411/** Adds an atom to the tesselation::PointsOnBoundary list.
    … …  
    27232667        (*it)->OptCenter.AddVector(&helper);
    27242668        Vector *center = mol->DetermineCenterOfAll(out);
    2725         (*it)->OptCenter.AddVector(center);
     2669        (*it)->OptCenter.SubtractVector(center);
    27262670        delete(center);
    27272671        // and add to file plus translucency object

  • src/molecules.cpp

    r3c4f04 rfcc9b0  
    739739};
    740740
    741 /** Returns vector pointing to center of gravity.
     741/** Returns vector pointing to center of all atoms.
    742742 * \param *out output stream for debugging
    743  * \return pointer to center of gravity vector
     743 * \return pointer to center of all vector
    744744 */
    745745Vector * molecule::DetermineCenterOfAll(ofstream *out)
    … …  
    759759      a->AddVector(&tmp);
    760760    }
    761     a->Scale(-1./Num); // divide through total mass (and sign for direction)
     761    a->Scale(1./Num); // divide through total mass (and sign for direction)
    762762  }
    763763  //cout << Verbose(1) << "Resulting center of gravity: ";
Note: See TracChangeset for help on using the changeset viewer.