Changes in / [27ac00:a7c344]


Ignore:
Location:
src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • src/Legacy/oldmenu.cpp

    r27ac00 ra7c344  
    88
    99#include "Legacy/oldmenu.hpp"
     10#include "analysis_bonds.hpp"
    1011#include "analysis_correlation.hpp"
    1112#include "World.hpp"
    … …  
    502503  Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;
    503504  Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
     505  Log() << Verbose(0) << " h - count the number of hydrogen bonds" << endl;
    504506  Log() << Verbose(0) << "all else - go back" << endl;
    505507  Log() << Verbose(0) << "===============================================" << endl;
    … …  
    592594        output->close();
    593595        delete(output);
     596      }
     597      break;
     598    case 'h':
     599      {
     600        int Z1;
     601        cout << "Please enter first interface element: ";
     602        cin >> Z1;
     603        const element * InterfaceElement = World::getInstance().getPeriode()->FindElement(Z1);
     604        int Z2;
     605        cout << "Please enter second interface element: ";
     606        cin >> Z2;
     607        const element * InterfaceElement2 = World::getInstance().getPeriode()->FindElement(Z2);
     608        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;
    594609      }
    595610      break;

  • src/analysis_bonds.cpp

    r27ac00 ra7c344  
    123123 * \param *molecules molecules to count bonds
    124124 * \param *InterfaceElement or NULL
    125  */
    126 int CountHydrogenBridgeBonds(MoleculeListClass *molecules, element * InterfaceElement = NULL)
     125 * \param *Interface2Element or NULL
     126 */
     127int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL, const element * Interface2Element = NULL)
    127128{
    128129  atom *Walker = NULL;
    … …  
    132133  double Otherangle = 0.;
    133134  bool InterfaceFlag = false;
     135  bool Interface2Flag = false;
    134136  bool OtherHydrogenFlag = true;
    135137  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    … …  
    151153              OtherHydrogens = 0;
    152154              InterfaceFlag = (InterfaceElement == NULL);
     155              Interface2Flag = (Interface2Element == NULL);
    153156              for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); BondRunner != Runner->ListOfBonds.end(); BondRunner++) {
    154157                atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Runner);
    … …  
    161164                }
    162165                InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement);
     166                Interface2Flag = Interface2Flag || (OtherAtom->type == Interface2Element);
    163167              }
    164168              DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl);
    … …  
    174178                  break;
    175179              }
    176               if (InterfaceFlag && OtherHydrogenFlag) {
     180              if (InterfaceFlag && Interface2Flag && OtherHydrogenFlag) {
    177181                // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
    178182                for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {

  • src/analysis_bonds.hpp

    r27ac00 ra7c344  
    3333void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, element *type1, element *type2, double &Min, double &Mean, double &Max);
    3434
    35 int CountHydrogenBridgeBonds(MoleculeListClass * const molecules, 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

    r27ac00 ra7c344  
    10391039//  // Purges surplus triangles.
    10401040//  TesselStruct->RemoveDegeneratedTriangles();
    1041 
    1042   // check envelope for consistency
    1043   status = CheckListOfBaselines(TesselStruct);
     1041//
     1042//  // check envelope for consistency
     1043//  status = CheckListOfBaselines(TesselStruct);
    10441044
    10451045  // store before correction

  • src/builder.cpp

    r27ac00 ra7c344  
    7070#include "molecule.hpp"
    7171#include "periodentafel.hpp"
     72#include "tesselationhelpers.hpp"
    7273#include "UIElements/UIFactory.hpp"
    7374#include "UIElements/MainWindow.hpp"
    … …  
    15191520        switch(argv[argptr-1][1]) {
    15201521          case 'h':
    1521           case 'H':
    15221522          case '?':
    15231523            DoLog(0) && (Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl);
    … …  
    15381538            DoLog(0) && (Log() << Verbose(0) << "\t-FF <MaxDistance> <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl);
    15391539            DoLog(0) && (Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl);
    1540             DoLog(0) && (Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl);
     1540            DoLog(0) && (Log() << Verbose(0) << "\t-h/-?\tGive this help screen." << endl);
     1541            DoLog(0) && (Log() << Verbose(0) << "\t-H\tCount Hydrogen bridge bonds." << endl);
    15411542            DoLog(0) && (Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl);
    15421543            DoLog(0) && (Log() << Verbose(0) << "\t-j\t<path> Store all bonds to file." << endl);
    … …  
    18191820              }
    18201821              break;
     1822            case 'H':
     1823              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-')) {
     1824                ExitFlag = 255;
     1825                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for calculating hydrogen bridge bonds: -H <Z1> <Z2>" << endl);
     1826                performCriticalExit();
     1827              } else {
     1828                if (ExitFlag == 0) ExitFlag = 1;
     1829                const element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
     1830                const element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+1]));
     1831                cout << endl << "There are " << CountHydrogenBridgeBonds(molecules, elemental, elemental2) << " hydrogen bridges with connections to " << (elemental != 0 ? elemental->name : "None") << " and " << (elemental2 != 0 ? elemental2->name : "None") << "." << endl;
     1832                argptr+=1;
     1833              }
     1834              break;
    18211835            case 'C':
    18221836              {
    … …  
    20892103              } else {
    20902104                class Tesselation *T = NULL;
     2105                class Tesselation *Convex = NULL;
    20912106                const LinkedCell *LCList = NULL;
     2107                const LinkedCell *LCListConvex = NULL;
    20922108                molecule * Boundary = NULL;
    20932109                //string filename(argv[argptr+1]);
    … …  
    21092125                if (!FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]))
    21102126                  ExitFlag = 255;
     2127                const double ConvexRadius = 20.;
     2128                LCListConvex = new LinkedCell(Boundary, 2.*ConvexRadius);
     2129//                setVerbosity(3);
     2130                if (!FindNonConvexBorder(Boundary, Convex, LCListConvex, ConvexRadius, "ConvexEnvelope"))
     2131                  ExitFlag = 255;
     2132                CalculateConstrictionPerBoundaryPoint(T, Convex);
     2133                StoreTrianglesinFile(mol, (const Tesselation *&)T, argv[argptr+1], "");
    21112134                //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
    21122135                end = clock();

  • src/tesselation.cpp

    r27ac00 ra7c344  
    230230{
    231231  Info FunctionInfo(__func__);
     232  double angle = CalculateConvexity();
     233  if (angle > -MYEPSILON) {
     234    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
     235    return true;
     236  } else {
     237    DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
     238    return false;
     239  }
     240}
     241
     242
     243/** Calculates the angle between two triangles with respect to their normal vector.
     244 * We sum the two angles of each height vector with respect to the center of the baseline.
     245 * \return angle > 0 then convex, if < 0 then concave
     246 */
     247double BoundaryLineSet::CalculateConvexity() const
     248{
     249  Info FunctionInfo(__func__);
    232250  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
    233251  // get the two triangles
    … …  
    278296  BaseLineNormal.Scale(-1.);
    279297  double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
    280   if ((angle - M_PI) > -MYEPSILON) {
    281     DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
    282     return true;
    283   } else {
    284     DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
    285     return false;
    286   }
     298  return (angle - M_PI);
    287299}
    288300
    … …  
    303315/** Returns other endpoint of the line.
    304316 * \param *point other endpoint
    305  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
     317 * \return NULL - if endpoint not contained in BoundaryLineSet::lines, or pointer to BoundaryPointSet otherwise
    306318 */
    307319class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
    … …  
    314326  else
    315327    return NULL;
     328}
     329;
     330
     331/** Returns other triangle of the line.
     332 * \param *point other endpoint
     333 * \return NULL - if triangle not contained in BoundaryLineSet::triangles, or pointer to BoundaryTriangleSet otherwise
     334 */
     335class BoundaryTriangleSet *BoundaryLineSet::GetOtherTriangle(const BoundaryTriangleSet * const triangle) const
     336{
     337  Info FunctionInfo(__func__);
     338  if (triangles.size() == 2) {
     339    for (TriangleMap::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); ++TriangleRunner)
     340      if (TriangleRunner->second != triangle)
     341        return TriangleRunner->second;
     342  }
     343  return NULL;
    316344}
    317345;
    … …  
    662690;
    663691
     692/** Returns the baseline which does not contain the given boundary point \a *point.
     693 * \param *point endpoint which is neither endpoint of the desired line
     694 * \return pointer to desired third baseline
     695 */
     696class BoundaryLineSet *BoundaryTriangleSet::GetThirdLine(const BoundaryPointSet * const point) const
     697{
     698  Info FunctionInfo(__func__);
     699  // sanity check
     700  if (!ContainsBoundaryPoint(point))
     701    return NULL;
     702  for (int i = 0; i < 3; i++)
     703    if (!lines[i]->ContainsBoundaryPoint(point))
     704      return lines[i];
     705  // actually, that' impossible :)
     706  return NULL;
     707}
     708;
     709
    664710/** Calculates the center point of the triangle.
    665711 * Is third of the sum of all endpoints.
    … …  
    11111157    TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));
    11121158
    1113     DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl);
     1159    DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << (*VRunner) << ":" << endl);
    11141160    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    1115       DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(OtherOptCenter) << "." << endl);
     1161      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(*(*VRunner)) << "." << endl);
    11161162
    11171163    // remove baseline's endpoints and candidates
    … …  
    11291175      DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
    11301176      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    1131         DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << endl);
     1177        DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << " at distance " << setprecision(13) << (*Runner)->node->distance(*(*VRunner)) << setprecision(6) << "." << endl);
     1178
     1179      // check with animate_sphere.tcl VMD script
     1180      if (ThirdPoint != NULL) {
     1181        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);
     1182      } else {
     1183        DoeLog(1) && (eLog() << Verbose(1) << "Check by: ... missing third point ..." << endl);
     1184        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);
     1185      }
    11321186    }
    11331187    delete (ListofPoints);
    11341188
    1135     // check with animate_sphere.tcl VMD script
    1136     if (ThirdPoint != NULL) {
    1137       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);
    1138     } else {
    1139       DoLog(1) && (Log() << Verbose(1) << "Check by: ... missing third point ..." << endl);
    1140       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);
    1141     }
    11421189  }
    11431190  return flag;
    … …  
    33183365                        }
    33193366                      } else {
    3320                         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);
     3367                        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);
    33213368                      }
    33223369                    } else {
    … …  
    46184665
    46194666  DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl);
    4620   IndexToIndex::iterator it;
    4621   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
     4667  for (IndexToIndex::iterator it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
    46224668    DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
    46234669
    … …  
    46374683  int count = 0;
    46384684
    4639   for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner) {
     4685  // iterate over all degenerated triangles
     4686  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); !DegeneratedTriangles->empty(); TriangleKeyRunner = DegeneratedTriangles->begin()) {
     4687    DoLog(0) && (Log() << Verbose(0) << "Checking presence of triangles " << TriangleKeyRunner->first << " and " << TriangleKeyRunner->second << "." << endl);
     4688    // both ways are stored in the map, only use one
     4689    if (TriangleKeyRunner->first > TriangleKeyRunner->second)
     4690      continue;
     4691
     4692    // determine from the keys in the map the two _present_ triangles
    46404693    finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
    46414694    if (finder != TrianglesOnBoundary.end())
    46424695      triangle = finder->second;
    46434696    else
    4644       break;
     4697      continue;
    46454698    finder = TrianglesOnBoundary.find(TriangleKeyRunner->second);
    46464699    if (finder != TrianglesOnBoundary.end())
    46474700      partnerTriangle = finder->second;
    46484701    else
    4649       break;
    4650 
     4702      continue;
     4703
     4704    // determine which lines are shared by the two triangles
    46514705    bool trianglesShareLine = false;
    46524706    for (int i = 0; i < 3; ++i)

  • src/tesselation.hpp

    r27ac00 ra7c344  
    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

    r27ac00 ra7c344  
    1010#include "info.hpp"
    1111#include "linkedcell.hpp"
     12#include "linearsystemofequations.hpp"
    1213#include "log.hpp"
    1314#include "tesselation.hpp"
    … …  
    183184  beta = M_PI - SideC.Angle(SideA);
    184185  gamma = M_PI - SideA.Angle(SideB);
    185   //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
     186  Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
    186187  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
    187188    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);
    … …  
    196197  (*Center) += helper;
    197198  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
     199  Log() << Verbose(1) << "INFO: Center (1st algo) is at " << *Center << "." << endl;
     200
     201//  LinearSystemOfEquations LSofEq(NDIM,NDIM);
     202//  double *matrix = new double[NDIM*NDIM];
     203//  matrix[0] = 0.;
     204//  matrix[1] = a.DistanceSquared(b);
     205//  matrix[2] = a.DistanceSquared(c);
     206//  matrix[3] = a.DistanceSquared(b);
     207//  matrix[4] = 0.;
     208//  matrix[5] = b.DistanceSquared(c);
     209//  matrix[6] = a.DistanceSquared(c);
     210//  matrix[7] = b.DistanceSquared(c);
     211//  matrix[8] = 0.;
     212//  cout << "Matrix is: ";
     213//  for (int i=0;i<NDIM*NDIM;i++)
     214//    cout << matrix[i] << "\t";
     215//  cout << endl;
     216//  LSofEq.SetA(matrix);
     217//  delete[](matrix);
     218//  LSofEq.Setb(new Vector(1.,1.,1.));
     219//  LSofEq.SetSymmetric(true);
     220//  helper.Zero();
     221//  if (!LSofEq.GetSolutionAsVector(helper)) {
     222//    DoLog(0) && (eLog()<< Verbose(0) << "Could not solve the linear system in GetCenterofCircumCircle()!" << endl);
     223//  }
     224//  cout << "Solution is " << helper << endl;
     225  // is equivalent to the three lines below
     226  helper[0] = SideA.NormSquared()*(SideB.NormSquared()+SideC.NormSquared() - SideA.NormSquared());
     227  helper[1] = SideB.NormSquared()*(SideC.NormSquared()+SideA.NormSquared() - SideB.NormSquared());
     228  helper[2] = SideC.NormSquared()*(SideA.NormSquared()+SideB.NormSquared() - SideC.NormSquared());
     229
     230  Center->Zero();
     231  *Center += helper[0] * a;
     232  *Center += helper[1] * b;
     233  *Center += helper[2] * c;
     234  Center->Scale(1./(helper[0]+helper[1]+helper[2]));
     235  Log() << Verbose(1) << "INFO: Center (2nd algo) is at " << *Center << "." << endl;
    198236};
    199237
    … …  
    417455/** Calculates the volume of a general tetraeder.
    418456 * \param *a first vector
    419  * \param *a first vector
    420  * \param *a first vector
    421  * \param *a first vector
     457 * \param *b second vector
     458 * \param *c third vector
     459 * \param *d fourth vector
    422460 * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
    423461 */
    … …  
    437475  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
    438476  return volume;
     477};
     478
     479/** Calculates the area of a general triangle.
     480 * We use the Heron's formula of area, [Bronstein, S. 138]
     481 * \param &A first vector
     482 * \param &B second vector
     483 * \param &C third vector
     484 * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
     485 */
     486double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C)
     487{
     488  Info FunctionInfo(__func__);
     489
     490  const double sidea = B.distance(C);
     491  const double sideb = A.distance(C);
     492  const double sidec = A.distance(B);
     493  const double s = (sidea+sideb+sidec)/2.;
     494
     495  const double area = sqrt(s*(s-sidea)*(s-sideb)*(s-sidec));
     496  return area;
    439497};
    440498
    … …  
    871929  class BoundaryPointSet *point = NULL;
    872930  class BoundaryLineSet *line = NULL;
    873 
    874   // calculate remaining concavity
     931  class BoundaryTriangleSet *triangle = NULL;
     932  double ConcavityPerLine = 0.;
     933  double ConcavityPerTriangle = 0.;
     934  double area = 0.;
     935  double totalarea = 0.;
     936
    875937  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
    876938    point = PointRunner->second;
    877939    DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
    878     point->value = 0;
     940
     941    // calculate mean concavity over all connected line
     942    ConcavityPerLine = 0.;
    879943    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    880944      line = LineRunner->second;
    881945      //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    882       if (!line->CheckConvexityCriterion())
    883         point->value += 1;
    884     }
    885   }
    886 };
    887 
     946      ConcavityPerLine -= line->CalculateConvexity();
     947    }
     948    ConcavityPerLine /= point->lines.size();
     949
     950    // weigh with total area of the surrounding triangles
     951    totalarea  = 0.;
     952    TriangleSet *triangles = TesselStruct->GetAllTriangles(PointRunner->second);
     953    for (TriangleSet::iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
     954      totalarea += CalculateAreaofGeneralTriangle(*(*TriangleRunner)->endpoints[0]->node->node, *(*TriangleRunner)->endpoints[1]->node->node, *(*TriangleRunner)->endpoints[2]->node->node);
     955    }
     956    ConcavityPerLine *= totalarea;
     957
     958    // calculate mean concavity over all attached triangles
     959    ConcavityPerTriangle = 0.;
     960    for (TriangleSet::const_iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
     961      line = (*TriangleRunner)->GetThirdLine(PointRunner->second);
     962      triangle = line->GetOtherTriangle(*TriangleRunner);
     963      area = CalculateAreaofGeneralTriangle(*triangle->endpoints[0]->node->node, *triangle->endpoints[1]->node->node, *triangle->endpoints[2]->node->node);
     964      area += CalculateAreaofGeneralTriangle(*(*TriangleRunner)->endpoints[0]->node->node, *(*TriangleRunner)->endpoints[1]->node->node, *(*TriangleRunner)->endpoints[2]->node->node);
     965      area *= -line->CalculateConvexity();
     966      if (area > 0)
     967        ConcavityPerTriangle += area;
     968//      else
     969//        ConcavityPerTriangle -= area;
     970    }
     971    ConcavityPerTriangle /= triangles->size()/totalarea;
     972    delete(triangles);
     973
     974    // add up
     975    point->value = ConcavityPerLine + ConcavityPerTriangle;
     976  }
     977};
     978
     979
     980
     981/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
     982 * Sets BoundaryPointSet::value equal to the nearest distance to convex envelope.
     983 * \param *out output stream for debugging
     984 * \param *TesselStruct pointer to Tesselation structure
     985 * \param *Convex pointer to convex Tesselation structure as reference
     986 */
     987void CalculateConstrictionPerBoundaryPoint(const Tesselation * const TesselStruct, const Tesselation * const Convex)
     988{
     989  Info FunctionInfo(__func__);
     990  double distance = 0.;
     991
     992  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     993    DoeLog(1) && (eLog() << Verbose(1) << "INFO: Current point is " << * PointRunner->second << "." << endl);
     994
     995    distance = 0.;
     996    for (TriangleMap::const_iterator TriangleRunner = Convex->TrianglesOnBoundary.begin(); TriangleRunner != Convex->TrianglesOnBoundary.end(); TriangleRunner++) {
     997      const double CurrentDistance = Convex->GetDistanceSquaredToTriangle(*PointRunner->second->node->node, TriangleRunner->second);
     998      if (CurrentDistance < distance)
     999        distance = CurrentDistance;
     1000    }
     1001
     1002    PointRunner->second->value = distance;
     1003  }
     1004};
    8881005
    8891006/** Checks whether each BoundaryLineSet in the Tesselation has two triangles.

  • src/tesselationhelpers.hpp

    r27ac00 ra7c344  
    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

    r27ac00 ra7c344  
    167167  Translator  = Vector(3,0,0);
    168168  TestMolecule2->Translate(&Translator);
    169   CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL) );
    170   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, oxygen) );
     169  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
     170  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, oxygen, NULL) );
    171171  //OutputTestMolecule(TestMolecule2, "testmolecule2-1.xyz");
    172172  Translator = Vector(-3,0,0);
    … …  
    176176  Translator = Vector(0,3,0);
    177177  TestMolecule2->Translate(&Translator);
    178   CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL) );
     178  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    179179  //OutputTestMolecule(TestMolecule2, "testmolecule2-2.xyz");
    180180  Translator = Vector(0,-3,0);
    … …  
    185185  TestMolecule2->Scale((const double ** const)&mirror);
    186186  TestMolecule2->Translate(&Translator);
    187   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     187  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    188188  //OutputTestMolecule(TestMolecule2, "testmolecule2-3.xyz");
    189189  Translator = Vector(0,3,0);
    … …  
    194194  Translator = Vector(2,1,0);
    195195  TestMolecule2->Translate(&Translator);
    196   CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL) );
     196  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    197197  //OutputTestMolecule(TestMolecule2, "testmolecule2-4.xyz");
    198198  Translator = Vector(-2,-1,0);
    … …  
    202202  Translator = Vector(0,0,3);
    203203  TestMolecule2->Translate(&Translator);
    204   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     204  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    205205  //OutputTestMolecule(TestMolecule2, "testmolecule2-5.xyz");
    206206  Translator = Vector(0,0,-3);
    … …  
    211211  TestMolecule2->Scale((const double ** const)&mirror);
    212212  TestMolecule2->Translate(&Translator);
    213   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     213  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    214214  //OutputTestMolecule(TestMolecule2, "testmolecule2-6.xyz");
    215215  Translator = Vector(3,0,0);
    … …  
    221221  TestMolecule2->Scale((const double ** const)&mirror);
    222222  TestMolecule2->Translate(&Translator);
    223   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     223  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    224224  //OutputTestMolecule(TestMolecule2, "testmolecule2-7.xyz");
    225225  Translator = Vector(-3,0,0);
    … …  
    232232  TestMolecule2->Translate(&Translator);
    233233  //OutputTestMolecule(TestMolecule2, "testmolecule2-8.xyz");
    234   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL) );
     234  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    235235  Translator = Vector(0,-3,0);
    236236  TestMolecule2->Translate(&Translator);
Note: See TracChangeset for help on using the changeset viewer.