source: src/Tesselation/CandidateForTesselation.cpp@ d94e7e

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * CandidateForTesselation.cpp
 *
 *  Created on: Jul 29, 2010
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "CandidateForTesselation.hpp"

#include <iostream>
#include <iomanip>

#include "BoundaryLineSet.hpp"
#include "BoundaryPointSet.hpp"
#include "Atom/TesselPoint.hpp"

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinkedCell/linkedcell.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_deprecated structure with other atoms
 * \return true - sphere is valid, false - sphere contains other points
 */
bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell_deprecated *LC) const
{
  //Info FunctionInfo(__func__);

  const double radiusSquared = RADIUS * RADIUS;
  list<const Vector *> VectorList;
  VectorList.push_back(&OptCenter);
  //VectorList.push_back(&OtherOptCenter);  // don't check the other (wrong) center

  if (!pointlist.empty())
    LOG(3, "DEBUG: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");
  else
    LOG(3, "DEBUG: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");
  // check baseline for OptCenter and OtherOptCenter being on sphere's surface
  for (list<const Vector *>::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) {
        ELOG(1, "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << ".");
        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 Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
      const double distance = fabs((*VRunner)->DistanceSquared(Walker->getPosition()) - radiusSquared);
      if (distance > HULLEPSILON) {
        ELOG(1, "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << ".");
        return false;
      } else {
        LOG(3, "DEBUG: Candidate " << *Walker << " is inside by " << distance << ".");
      }
    }
  }

  LOG(2, "DEBUG: Checking whether sphere contains no others points ...");
  bool flag = true;
  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
    // get all points inside the sphere
    TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));

    {
      LOG(3, "DEBUG: The following atoms are inside sphere at " << (*VRunner) << ":");
      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
        LOG(3, "DEBUG:  " << *(*Runner) << " with distance " << (*Runner)->distance(*(*VRunner)) << ".");
    }

    // remove baseline's endpoints and candidates
    for (int i = 0; i < 2; i++) {
      LOG(3, "DEBUG: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << ".");
      ListofPoints->remove(BaseLine->endpoints[i]->node);
    }
    for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
      LOG(3, "DEBUG: removing candidate tesselpoint " << *(*Runner) << ".");
      ListofPoints->remove(*Runner);
    }
    if (!ListofPoints->empty()) {
      ELOG(1, "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere.");
      flag = false;
      ELOG(1, "External atoms inside of sphere at " << *(*VRunner) << ":");
      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
        ELOG(1, "  " << *(*Runner) << " at distance " << setprecision(13) << (*Runner)->distance(*(*VRunner)) << setprecision(6) << ".");

      // check with animate_sphere.tcl VMD script
      if (ThirdPoint != NULL) {
        ELOG(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));
      } else {
        ELOG(1, "Check by: ... missing third point ...");
        ELOG(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));
      }
    }
    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;
}
;