/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * CandidateForTesselation.cpp * * Created on: Jul 29, 2010 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "Helpers/MemDebug.hpp" #include "CandidateForTesselation.hpp" #include #include #include "BoundaryLineSet.hpp" #include "BoundaryPointSet.hpp" #include "TesselPoint.hpp" #include "Helpers/Assert.hpp" #include "Helpers/Info.hpp" //#include "Line.hpp" #include "linkedcell.hpp" #include "Helpers/Log.hpp" //#include "Plane.hpp" #include "LinearAlgebra/Vector.hpp" #include "Helpers/Verbose.hpp" using namespace std; const double CandidateForTesselation::HULLEPSILON = 1e-9; /** Constructor of class CandidateForTesselation. */ CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) : BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI) { Info FunctionInfo(__func__); } ; /** Constructor of class CandidateForTesselation. */ CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, const Vector &OptCandidateCenter, const Vector &OtherOptCandidateCenter) : BaseLine(line), ThirdPoint(point), T(NULL), OptCenter(OptCandidateCenter), OtherOptCenter(OtherOptCandidateCenter), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI) { Info FunctionInfo(__func__); }; /** Destructor for class CandidateForTesselation. */ CandidateForTesselation::~CandidateForTesselation() { } ; /** Checks validity of a given sphere of a candidate line. * Sphere must touch all candidates and the baseline endpoints and there must be no other atoms inside. * \param RADIUS radius of sphere * \param *LC LinkedCell structure with other atoms * \return true - sphere is valid, false - sphere contains other points */ bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell *LC) const { Info FunctionInfo(__func__); const double radiusSquared = RADIUS * RADIUS; list VectorList; VectorList.push_back(&OptCenter); //VectorList.push_back(&OtherOptCenter); // don't check the other (wrong) center if (!pointlist.empty()) DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl); else DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl); // check baseline for OptCenter and OtherOptCenter being on sphere's surface for (list::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { for (int i = 0; i < 2; i++) { const double distance = fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->getPosition()) - radiusSquared); if (distance > HULLEPSILON) { DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl); return false; } } } // check Candidates for OptCenter and OtherOptCenter being on sphere's surface for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) { const TesselPoint *Walker = *Runner; for (list::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { const double distance = fabs((*VRunner)->DistanceSquared(Walker->getPosition()) - radiusSquared); if (distance > HULLEPSILON) { DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl); return false; } else { DoLog(1) && (Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl); } } } DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl); bool flag = true; for (list::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { // get all points inside the sphere TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner)); DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << (*VRunner) << ":" << endl); for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) DoLog(1) && (Log() << Verbose(1) << " " << *(*Runner) << " with distance " << (*Runner)->distance(*(*VRunner)) << "." << endl); // remove baseline's endpoints and candidates for (int i = 0; i < 2; i++) { DoLog(1) && (Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl); ListofPoints->remove(BaseLine->endpoints[i]->node); } for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) { DoLog(1) && (Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl); ListofPoints->remove(*Runner); } if (!ListofPoints->empty()) { DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl); flag = false; DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl); for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) DoeLog(1) && (eLog() << Verbose(1) << " " << *(*Runner) << " at distance " << setprecision(13) << (*Runner)->distance(*(*VRunner)) << setprecision(6) << "." << endl); // check with animate_sphere.tcl VMD script if (ThirdPoint != NULL) { 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); } else { DoeLog(1) && (eLog() << Verbose(1) << "Check by: ... missing third point ..." << endl); 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); } } delete (ListofPoints); } return flag; } ; /** output operator for CandidateForTesselation. * \param &ost output stream * \param &a boundary line */ ostream & operator <<(ostream &ost, const CandidateForTesselation &a) { ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->getName() << "," << a.BaseLine->endpoints[1]->node->getName() << "] with "; if (a.pointlist.empty()) ost << "no candidate."; else { ost << "candidate"; if (a.pointlist.size() != 1) ost << "s "; else ost << " "; for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++) ost << *(*Runner) << " "; ost << " at angle " << (a.ShortestAngle) << "."; } return ost; } ;