Changeset bdb143 for src


Ignore:
Timestamp:
Jun 28, 2010, 3:38:05 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:
13892b, 51a5c8
Parents:
b2531f (diff), c39cc4 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Frederik Heber <heber@…> (06/28/10 14:53:18)
git-committer:
Frederik Heber <heber@…> (06/28/10 15:38:05)
Message:

Merge branch 'StructureRefactoring' into stable

Conflicts:

molecuilder/src/analysis_bonds.cpp
molecuilder/src/analysis_bonds.hpp
molecuilder/src/builder.cpp
molecuilder/src/unittests/Makefile.am
molecuilder/src/unittests/TestRunnerMain.cpp

  • TESTFIX: Tesselations/heptan - Due to different implementations of GetCircumCenter() sequence of triangles nodes has changed in degenerate case
  • TESTFIX: Tesselation/1-3 - Due to convexity criterion fourth argument (i.e. the value at the node) has changed.

Signed-off-by: Frederik Heber <heber@…>

Location:
src
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • src/Legacy/oldmenu.cpp

    rb2531f rbdb143  
    1010
    1111#include "Legacy/oldmenu.hpp"
     12#include "analysis_bonds.hpp"
    1213#include "analysis_correlation.hpp"
    1314#include "World.hpp"
    … …  
    509510  Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;
    510511  Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
     512  Log() << Verbose(0) << " h - count the number of hydrogen bonds" << endl;
    511513  Log() << Verbose(0) << "all else - go back" << endl;
    512514  Log() << Verbose(0) << "===============================================" << endl;
    … …  
    597599        output->close();
    598600        delete(output);
     601      }
     602      break;
     603    case 'h':
     604      {
     605        int Z1;
     606        cout << "Please enter first interface element: ";
     607        cin >> Z1;
     608        const element * InterfaceElement = World::getInstance().getPeriode()->FindElement(Z1);
     609        int Z2;
     610        cout << "Please enter second interface element: ";
     611        cin >> Z2;
     612        const element * InterfaceElement2 = World::getInstance().getPeriode()->FindElement(Z2);
     613        cout << endl << "There are " << CountHydrogenBridgeBonds(World::getInstance().getMolecules(), InterfaceElement, InterfaceElement2) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << " and " << (InterfaceElement2 != 0 ? InterfaceElement2->name : "None") << "." << endl;
    599614      }
    600615      break;

  • src/analysis_bonds.cpp

    rb2531f rbdb143  
    121121 * \param *molecules molecules to count bonds
    122122 * \param *InterfaceElement or NULL
    123  */
    124 int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL)
     123 * \param *Interface2Element or NULL
     124 */
     125int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL, const element * Interface2Element = NULL)
    125126{
    126127  int count = 0;
    … …  
    128129  double Otherangle = 0.;
    129130  bool InterfaceFlag = false;
     131  bool Interface2Flag = false;
    130132  bool OtherHydrogenFlag = true;
    131133  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); ++MolWalker) {
    … …  
    145147              OtherHydrogens = 0;
    146148              InterfaceFlag = (InterfaceElement == NULL);
     149              Interface2Flag = (Interface2Element == NULL);
    147150              for (BondList::const_iterator BondRunner = (*Runner)->ListOfBonds.begin(); BondRunner != (*Runner)->ListOfBonds.end(); BondRunner++) {
    148151                atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
    … …  
    155158                }
    156159                InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement);
     160                Interface2Flag = Interface2Flag || (OtherAtom->type == Interface2Element);
    157161              }
    158162              DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl);
    … …  
    168172                  break;
    169173              }
    170               if (InterfaceFlag && OtherHydrogenFlag) {
     174              if (InterfaceFlag && Interface2Flag && OtherHydrogenFlag) {
    171175                // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
    172176                for (BondList::const_iterator BondRunner = (*Walker)->ListOfBonds.begin(); BondRunner != (*Walker)->ListOfBonds.end(); BondRunner++) {

  • src/analysis_bonds.hpp

    rb2531f rbdb143  
    3333void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, const element *type1, const element *type2, double &Min, double &Mean, double &Max);
    3434
    35 int CountHydrogenBridgeBonds(MoleculeListClass * const molecules, const element * InterfaceElement);
     35int CountHydrogenBridgeBonds(MoleculeListClass * const molecules, const element * InterfaceElement, const element * Interface2Element);
    3636int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second);
    3737int CountBondsOfThree(MoleculeListClass * const molecules, const element * const first, const element * const second, const element * const third);

  • src/boundary.cpp

    rb2531f rbdb143  
    10041004//  // Purges surplus triangles.
    10051005//  TesselStruct->RemoveDegeneratedTriangles();
    1006 
    1007   // check envelope for consistency
    1008   status = CheckListOfBaselines(TesselStruct);
     1006//
     1007//  // check envelope for consistency
     1008//  status = CheckListOfBaselines(TesselStruct);
    10091009
    10101010  cout << "before correction" << endl;

  • src/builder.cpp

    rb2531f rbdb143  
    6363#include "molecule.hpp"
    6464#include "periodentafel.hpp"
     65#include "tesselationhelpers.hpp"
    6566#include "UIElements/UIFactory.hpp"
    6667#include "UIElements/TextUI/TextUIFactory.hpp"

  • src/tesselation.cpp

    rb2531f rbdb143  
    231231{
    232232  Info FunctionInfo(__func__);
     233  double angle = CalculateConvexity();
     234  if (angle > -MYEPSILON) {
     235    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
     236    return true;
     237  } else {
     238    DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
     239    return false;
     240  }
     241}
     242
     243
     244/** Calculates the angle between two triangles with respect to their normal vector.
     245 * We sum the two angles of each height vector with respect to the center of the baseline.
     246 * \return angle > 0 then convex, if < 0 then concave
     247 */
     248double BoundaryLineSet::CalculateConvexity() const
     249{
     250  Info FunctionInfo(__func__);
    233251  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
    234252  // get the two triangles
    … …  
    279297  BaseLineNormal.Scale(-1.);
    280298  double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
    281   if ((angle - M_PI) > -MYEPSILON) {
    282     DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
    283     return true;
    284   } else {
    285     DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
    286     return false;
    287   }
     299  return (angle - M_PI);
    288300}
    289301
    … …  
    304316/** Returns other endpoint of the line.
    305317 * \param *point other endpoint
    306  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
     318 * \return NULL - if endpoint not contained in BoundaryLineSet::lines, or pointer to BoundaryPointSet otherwise
    307319 */
    308320class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
    … …  
    315327  else
    316328    return NULL;
     329}
     330;
     331
     332/** Returns other triangle of the line.
     333 * \param *point other endpoint
     334 * \return NULL - if triangle not contained in BoundaryLineSet::triangles, or pointer to BoundaryTriangleSet otherwise
     335 */
     336class BoundaryTriangleSet *BoundaryLineSet::GetOtherTriangle(const BoundaryTriangleSet * const triangle) const
     337{
     338  Info FunctionInfo(__func__);
     339  if (triangles.size() == 2) {
     340    for (TriangleMap::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); ++TriangleRunner)
     341      if (TriangleRunner->second != triangle)
     342        return TriangleRunner->second;
     343  }
     344  return NULL;
    317345}
    318346;
    … …  
    661689;
    662690
     691/** Returns the baseline which does not contain the given boundary point \a *point.
     692 * \param *point endpoint which is neither endpoint of the desired line
     693 * \return pointer to desired third baseline
     694 */
     695class BoundaryLineSet *BoundaryTriangleSet::GetThirdLine(const BoundaryPointSet * const point) const
     696{
     697  Info FunctionInfo(__func__);
     698  // sanity check
     699  if (!ContainsBoundaryPoint(point))
     700    return NULL;
     701  for (int i = 0; i < 3; i++)
     702    if (!lines[i]->ContainsBoundaryPoint(point))
     703      return lines[i];
     704  // actually, that' impossible :)
     705  return NULL;
     706}
     707;
     708
    663709/** Calculates the center point of the triangle.
    664710 * Is third of the sum of all endpoints.
    … …  
    11101156    TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));
    11111157
    1112     DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl);
     1158    DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << (*VRunner) << ":" << endl);
    11131159    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    1114       DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(OtherOptCenter) << "." << endl);
     1160      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(*(*VRunner)) << "." << endl);
    11151161
    11161162    // remove baseline's endpoints and candidates
    … …  
    11281174      DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
    11291175      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    1130         DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << endl);
     1176        DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << " at distance " << setprecision(13) << (*Runner)->node->distance(*(*VRunner)) << setprecision(6) << "." << endl);
     1177
     1178      // check with animate_sphere.tcl VMD script
     1179      if (ThirdPoint != NULL) {
     1180        DoeLog(1) && (eLog() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
     1181      } else {
     1182        DoeLog(1) && (eLog() << Verbose(1) << "Check by: ... missing third point ..." << endl);
     1183        DoeLog(1) && (eLog() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
     1184      }
    11311185    }
    11321186    delete (ListofPoints);
    11331187
    1134     // check with animate_sphere.tcl VMD script
    1135     if (ThirdPoint != NULL) {
    1136       DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
    1137     } else {
    1138       DoLog(1) && (Log() << Verbose(1) << "Check by: ... missing third point ..." << endl);
    1139       DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
    1140     }
    11411188  }
    11421189  return flag;
    … …  
    32903337                        }
    32913338                      } else {
    3292                         DoLog(1) && (Log() << Verbose(1) << "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius) << endl);
     3339                        DoeLog(0) && (eLog() << Verbose(1) << "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius) << endl);
    32933340                      }
    32943341                    } else {
    … …  
    45904637
    45914638  DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl);
    4592   IndexToIndex::iterator it;
    4593   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
     4639  for (IndexToIndex::iterator it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
    45944640    DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
    45954641
    … …  
    46094655  int count = 0;
    46104656
    4611   for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner) {
     4657  // iterate over all degenerated triangles
     4658  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); !DegeneratedTriangles->empty(); TriangleKeyRunner = DegeneratedTriangles->begin()) {
     4659    DoLog(0) && (Log() << Verbose(0) << "Checking presence of triangles " << TriangleKeyRunner->first << " and " << TriangleKeyRunner->second << "." << endl);
     4660    // both ways are stored in the map, only use one
     4661    if (TriangleKeyRunner->first > TriangleKeyRunner->second)
     4662      continue;
     4663
     4664    // determine from the keys in the map the two _present_ triangles
    46124665    finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
    46134666    if (finder != TrianglesOnBoundary.end())
    46144667      triangle = finder->second;
    46154668    else
    4616       break;
     4669      continue;
    46174670    finder = TrianglesOnBoundary.find(TriangleKeyRunner->second);
    46184671    if (finder != TrianglesOnBoundary.end())
    46194672      partnerTriangle = finder->second;
    46204673    else
    4621       break;
    4622 
     4674      continue;
     4675
     4676    // determine which lines are shared by the two triangles
    46234677    bool trianglesShareLine = false;
    46244678    for (int i = 0; i < 3; ++i)

  • src/tesselation.hpp

    rb2531f rbdb143  
    138138    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
    139139    bool CheckConvexityCriterion() const;
     140    double CalculateConvexity() const;
    140141    class BoundaryPointSet *GetOtherEndpoint(const BoundaryPointSet * const point) const;
     142    class BoundaryTriangleSet *GetOtherTriangle(const BoundaryTriangleSet * const triangle) const;
    141143
    142144    class BoundaryPointSet *endpoints[2];
    … …  
    164166    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
    165167    class BoundaryPointSet *GetThirdEndpoint(const BoundaryLineSet * const line) const;
     168    class BoundaryLineSet *GetThirdLine(const BoundaryPointSet * const point) const;
    166169    bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
    167170    bool IsPresentTupel(const BoundaryTriangleSet * const T) const;

  • src/tesselationhelpers.cpp

    rb2531f rbdb143  
    1212#include "info.hpp"
    1313#include "linkedcell.hpp"
     14#include "linearsystemofequations.hpp"
    1415#include "log.hpp"
    1516#include "tesselation.hpp"
    … …  
    185186  beta = M_PI - SideC.Angle(SideA);
    186187  gamma = M_PI - SideA.Angle(SideB);
    187   //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
     188  Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
    188189  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
    189190    DoeLog(2) && (eLog()<< Verbose(2) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl);
    … …  
    198199  (*Center) += helper;
    199200  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
     201  Log() << Verbose(1) << "INFO: Center (1st algo) is at " << *Center << "." << endl;
     202
     203//  LinearSystemOfEquations LSofEq(NDIM,NDIM);
     204//  double *matrix = new double[NDIM*NDIM];
     205//  matrix[0] = 0.;
     206//  matrix[1] = a.DistanceSquared(b);
     207//  matrix[2] = a.DistanceSquared(c);
     208//  matrix[3] = a.DistanceSquared(b);
     209//  matrix[4] = 0.;
     210//  matrix[5] = b.DistanceSquared(c);
     211//  matrix[6] = a.DistanceSquared(c);
     212//  matrix[7] = b.DistanceSquared(c);
     213//  matrix[8] = 0.;
     214//  cout << "Matrix is: ";
     215//  for (int i=0;i<NDIM*NDIM;i++)
     216//    cout << matrix[i] << "\t";
     217//  cout << endl;
     218//  LSofEq.SetA(matrix);
     219//  delete[](matrix);
     220//  LSofEq.Setb(new Vector(1.,1.,1.));
     221//  LSofEq.SetSymmetric(true);
     222//  helper.Zero();
     223//  if (!LSofEq.GetSolutionAsVector(helper)) {
     224//    DoLog(0) && (eLog()<< Verbose(0) << "Could not solve the linear system in GetCenterofCircumCircle()!" << endl);
     225//  }
     226//  cout << "Solution is " << helper << endl;
     227  // is equivalent to the three lines below
     228  helper[0] = SideA.NormSquared()*(SideB.NormSquared()+SideC.NormSquared() - SideA.NormSquared());
     229  helper[1] = SideB.NormSquared()*(SideC.NormSquared()+SideA.NormSquared() - SideB.NormSquared());
     230  helper[2] = SideC.NormSquared()*(SideA.NormSquared()+SideB.NormSquared() - SideC.NormSquared());
     231
     232  Center->Zero();
     233  *Center += helper[0] * a;
     234  *Center += helper[1] * b;
     235  *Center += helper[2] * c;
     236  Center->Scale(1./(helper[0]+helper[1]+helper[2]));
     237  Log() << Verbose(1) << "INFO: Center (2nd algo) is at " << *Center << "." << endl;
    200238};
    201239
    … …  
    419457/** Calculates the volume of a general tetraeder.
    420458 * \param *a first vector
    421  * \param *a first vector
    422  * \param *a first vector
    423  * \param *a first vector
     459 * \param *b second vector
     460 * \param *c third vector
     461 * \param *d fourth vector
    424462 * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
    425463 */
    … …  
    439477  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
    440478  return volume;
     479};
     480
     481/** Calculates the area of a general triangle.
     482 * We use the Heron's formula of area, [Bronstein, S. 138]
     483 * \param &A first vector
     484 * \param &B second vector
     485 * \param &C third vector
     486 * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
     487 */
     488double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C)
     489{
     490  Info FunctionInfo(__func__);
     491
     492  const double sidea = B.distance(C);
     493  const double sideb = A.distance(C);
     494  const double sidec = A.distance(B);
     495  const double s = (sidea+sideb+sidec)/2.;
     496
     497  const double area = sqrt(s*(s-sidea)*(s-sideb)*(s-sidec));
     498  return area;
    441499};
    442500
    … …  
    882940  class BoundaryPointSet *point = NULL;
    883941  class BoundaryLineSet *line = NULL;
    884 
    885   // calculate remaining concavity
     942  class BoundaryTriangleSet *triangle = NULL;
     943  double ConcavityPerLine = 0.;
     944  double ConcavityPerTriangle = 0.;
     945  double area = 0.;
     946  double totalarea = 0.;
     947
    886948  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
    887949    point = PointRunner->second;
    888950    DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
    889     point->value = 0;
     951
     952    // calculate mean concavity over all connected line
     953    ConcavityPerLine = 0.;
    890954    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    891955      line = LineRunner->second;
    892956      //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    893       if (!line->CheckConvexityCriterion())
    894         point->value += 1;
    895     }
    896   }
    897 };
    898 
     957      ConcavityPerLine -= line->CalculateConvexity();
     958    }
     959    ConcavityPerLine /= point->lines.size();
     960
     961    // weigh with total area of the surrounding triangles
     962    totalarea  = 0.;
     963    TriangleSet *triangles = TesselStruct->GetAllTriangles(PointRunner->second);
     964    for (TriangleSet::iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
     965      totalarea += CalculateAreaofGeneralTriangle(*(*TriangleRunner)->endpoints[0]->node->node, *(*TriangleRunner)->endpoints[1]->node->node, *(*TriangleRunner)->endpoints[2]->node->node);
     966    }
     967    ConcavityPerLine *= totalarea;
     968
     969    // calculate mean concavity over all attached triangles
     970    ConcavityPerTriangle = 0.;
     971    for (TriangleSet::const_iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
     972      line = (*TriangleRunner)->GetThirdLine(PointRunner->second);
     973      triangle = line->GetOtherTriangle(*TriangleRunner);
     974      area = CalculateAreaofGeneralTriangle(*triangle->endpoints[0]->node->node, *triangle->endpoints[1]->node->node, *triangle->endpoints[2]->node->node);
     975      area += CalculateAreaofGeneralTriangle(*(*TriangleRunner)->endpoints[0]->node->node, *(*TriangleRunner)->endpoints[1]->node->node, *(*TriangleRunner)->endpoints[2]->node->node);
     976      area *= -line->CalculateConvexity();
     977      if (area > 0)
     978        ConcavityPerTriangle += area;
     979//      else
     980//        ConcavityPerTriangle -= area;
     981    }
     982    ConcavityPerTriangle /= triangles->size()/totalarea;
     983    delete(triangles);
     984
     985    // add up
     986    point->value = ConcavityPerLine + ConcavityPerTriangle;
     987  }
     988};
     989
     990
     991
     992/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
     993 * Sets BoundaryPointSet::value equal to the nearest distance to convex envelope.
     994 * \param *out output stream for debugging
     995 * \param *TesselStruct pointer to Tesselation structure
     996 * \param *Convex pointer to convex Tesselation structure as reference
     997 */
     998void CalculateConstrictionPerBoundaryPoint(const Tesselation * const TesselStruct, const Tesselation * const Convex)
     999{
     1000  Info FunctionInfo(__func__);
     1001  double distance = 0.;
     1002
     1003  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     1004    DoeLog(1) && (eLog() << Verbose(1) << "INFO: Current point is " << * PointRunner->second << "." << endl);
     1005
     1006    distance = 0.;
     1007    for (TriangleMap::const_iterator TriangleRunner = Convex->TrianglesOnBoundary.begin(); TriangleRunner != Convex->TrianglesOnBoundary.end(); TriangleRunner++) {
     1008      const double CurrentDistance = Convex->GetDistanceSquaredToTriangle(*PointRunner->second->node->node, TriangleRunner->second);
     1009      if (CurrentDistance < distance)
     1010        distance = CurrentDistance;
     1011    }
     1012
     1013    PointRunner->second->value = distance;
     1014  }
     1015};
    8991016
    9001017/** Checks whether each BoundaryLineSet in the Tesselation has two triangles.

  • src/tesselationhelpers.hpp

    rb2531f rbdb143  
    4343/********************************************** definitions *********************************/
    4444
    45 #define HULLEPSILON 1e-10
     45#define HULLEPSILON 1e-9 //!< TODO: Get rid of HULLEPSILON, points to numerical instabilities
    4646
    4747/********************************************** declarations *******************************/
    … …  
    5555bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4);
    5656double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d);
     57double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C);
    5758double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector);
    5859
    … …  
    6869void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud);
    6970void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct);
     71void CalculateConstrictionPerBoundaryPoint(const Tesselation * const TesselStruct, const Tesselation * const Convex);
    7072double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle);
    7173

  • src/unittests/CountBondsUnitTest.cpp

    rb2531f rbdb143  
    154154  Translator  = Vector(3,0,0);
    155155  TestMolecule2->Translate(&Translator);
    156   CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL) );
    157   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, oxygen) );
     156  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
     157  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, oxygen, NULL) );
    158158  //OutputTestMolecule(TestMolecule2, "testmolecule2-1.xyz");
    159159  Translator = Vector(-3,0,0);
    … …  
    163163  Translator = Vector(0,3,0);
    164164  TestMolecule2->Translate(&Translator);
    165   CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL) );
     165  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    166166  //OutputTestMolecule(TestMolecule2, "testmolecule2-2.xyz");
    167167  Translator = Vector(0,-3,0);
    … …  
    172172  TestMolecule2->Scale((const double ** const)&mirror);
    173173  TestMolecule2->Translate(&Translator);
    174   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     174  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    175175  //OutputTestMolecule(TestMolecule2, "testmolecule2-3.xyz");
    176176  Translator = Vector(0,3,0);
    … …  
    181181  Translator = Vector(2,1,0);
    182182  TestMolecule2->Translate(&Translator);
    183   CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL) );
     183  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    184184  //OutputTestMolecule(TestMolecule2, "testmolecule2-4.xyz");
    185185  Translator = Vector(-2,-1,0);
    … …  
    189189  Translator = Vector(0,0,3);
    190190  TestMolecule2->Translate(&Translator);
    191   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     191  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    192192  //OutputTestMolecule(TestMolecule2, "testmolecule2-5.xyz");
    193193  Translator = Vector(0,0,-3);
    … …  
    198198  TestMolecule2->Scale((const double ** const)&mirror);
    199199  TestMolecule2->Translate(&Translator);
    200   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     200  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    201201  //OutputTestMolecule(TestMolecule2, "testmolecule2-6.xyz");
    202202  Translator = Vector(3,0,0);
    … …  
    208208  TestMolecule2->Scale((const double ** const)&mirror);
    209209  TestMolecule2->Translate(&Translator);
    210   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     210  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    211211  //OutputTestMolecule(TestMolecule2, "testmolecule2-7.xyz");
    212212  Translator = Vector(-3,0,0);
    … …  
    219219  TestMolecule2->Translate(&Translator);
    220220  //OutputTestMolecule(TestMolecule2, "testmolecule2-8.xyz");
    221   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     221  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    222222  Translator = Vector(0,-3,0);
    223223  TestMolecule2->Translate(&Translator);

  • src/unittests/TestRunnerMain.cpp

    rb2531f rbdb143  
    1717#include "../../../TestRunnerClient.h"
    1818#include "../../../TestRunnerClient.cpp"
     19
    1920#else
    20 #include "UnitTestMain.cpp"
    21 #endif
    2221
     22#include "unittests/UnitTestMain.cpp"
     23
     24#endif /* HAVE_ECUT */

  • src/unittests/UnitTestMain.cpp

    rb2531f rbdb143  
    66 */
    77
    8 #include <cppunit/CompilerOutputter.h>
    9 #include <cppunit/extensions/TestFactoryRegistry.h>
    10 #include <cppunit/ui/text/TestRunner.h>
    11 
    128// include config.h
    139#ifdef HAVE_CONFIG_H
    1410#include <config.h>
    1511#endif
     12
     13#include <cppunit/CompilerOutputter.h>
     14#include <cppunit/extensions/TestFactoryRegistry.h>
     15#include <cppunit/ui/text/TestRunner.h>
    1616
    1717/********************************************** Main routine **************************************/
    … …  
    3535  return wasSucessful ? 0 : 1;
    3636};
     37
Note: See TracChangeset for help on using the changeset viewer.