source: src/Tesselation/BoundaryLineSet.cpp@ 693a80

/*
 * 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.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "BoundaryLineSet.hpp"

#include <iostream>

#include "BoundaryPointSet.hpp"
#include "BoundaryTriangleSet.hpp"
#include "TesselPoint.hpp"

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "tesselationhelpers.hpp"
#include "LinearAlgebra/Vector.hpp"

using namespace std;

/** Constructor of BoundaryLineSet.
 */
BoundaryLineSet::BoundaryLineSet() :
  Nr(-1)
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 2; i++)
    endpoints[i] = NULL;
}
;

/** Constructor of BoundaryLineSet with two endpoints.
 * Adds line automatically to each endpoints' LineMap
 * \param *Point[2] array of two boundary points
 * \param number number of the list
 */
BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
{
  Info FunctionInfo(__func__);
  // set number
  Nr = number;
  // set endpoints in ascending order
  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
  // add this line to the hash maps of both endpoints
  Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
  Point[1]->AddLine(this); //
  // set skipped to false
  skipped = false;
  // clear triangles list
  LOG(0, "New Line with endpoints " << *this << ".");
}
;

/** Constructor of BoundaryLineSet with two endpoints.
 * Adds line automatically to each endpoints' LineMap
 * \param *Point1 first boundary point
 * \param *Point2 second boundary point
 * \param number number of the list
 */
BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number) :
  Nr(number),
  skipped(false)
{
  Info FunctionInfo(__func__);
  // set endpoints in ascending order
  SetEndpointsOrdered(endpoints, Point1, Point2);
  // add this line to the hash maps of both endpoints
  Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
  Point2->AddLine(this); //
  // clear triangles list
  LOG(0, "New Line with endpoints " << *this << ".");
}
;

/** Destructor for BoundaryLineSet.
 * Removes itself from each endpoints' LineMap, calling RemoveTrianglePoint() when point not connected anymore.
 * \note When removing lines from a class Tesselation, use RemoveTesselationLine()
 */
BoundaryLineSet::~BoundaryLineSet()
{
  Info FunctionInfo(__func__);
  int Numbers[2];

  // get other endpoint number of finding copies of same line
  if (endpoints[1] != NULL)
    Numbers[0] = endpoints[1]->Nr;
  else
    Numbers[0] = -1;
  if (endpoints[0] != NULL)
    Numbers[1] = endpoints[0]->Nr;
  else
    Numbers[1] = -1;

  for (int i = 0; i < 2; i++) {
    if (endpoints[i] != NULL) {
      if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
        pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
        for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
          if ((*Runner).second == this) {
            //LOG(0, "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << ".");
            endpoints[i]->lines.erase(Runner);
            break;
          }
      } else { // there's just a single line left
        if (endpoints[i]->lines.erase(Nr)) {
          //LOG(0, "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << ".");
        }
      }
      if (endpoints[i]->lines.empty()) {
        //LOG(0, *endpoints[i] << " has no more lines it's attached to, erasing.");
        if (endpoints[i] != NULL) {
          delete (endpoints[i]);
          endpoints[i] = NULL;
        }
      }
    }
  }
  if (!triangles.empty())
    ELOG(2, "Memory Leak! I " << *this << " am still connected to some triangles.");
}
;

/** Add triangle to TriangleMap of this boundary line.
 * \param *triangle to add
 */
void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
{
  Info FunctionInfo(__func__);
  LOG(0, "Add " << triangle->Nr << " to line " << *this << ".");
  triangles.insert(TrianglePair(triangle->Nr, triangle));
}
;

/** Checks whether we have a common endpoint with given \a *line.
 * \param *line other line to test
 * \return true - common endpoint present, false - not connected
 */
bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
{
  Info FunctionInfo(__func__);
  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
    return true;
  else
    return false;
}
;

/** Checks whether the adjacent triangles of a baseline are convex or not.
 * We sum the two angles of each height vector with respect to the center of the baseline.
 * If greater/equal M_PI than we are convex.
 * \param *out output stream for debugging
 * \return true - triangles are convex, false - concave or less than two triangles connected
 */
bool BoundaryLineSet::CheckConvexityCriterion() const
{
  Info FunctionInfo(__func__);
  double angle = CalculateConvexity();
  if (angle > -MYEPSILON) {
    LOG(0, "ACCEPT: Angle is greater than pi: convex.");
    return true;
  } else {
    LOG(0, "REJECT: Angle is less than pi: concave.");
    return false;
  }
}


/** Calculates the angle between two triangles with respect to their normal vector.
 * We sum the two angles of each height vector with respect to the center of the baseline.
 * \return angle > 0 then convex, if < 0 then concave
 */
double BoundaryLineSet::CalculateConvexity() const
{
  Info FunctionInfo(__func__);
  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
  // get the two triangles
  if (triangles.size() != 2) {
    ELOG(0, "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!");
    return true;
  }
  // check normal vectors
  // have a normal vector on the base line pointing outwards
  //LOG(0, "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << ".");
  BaseLineCenter = (1./2.)*((endpoints[0]->node->getPosition()) + (endpoints[1]->node->getPosition()));
  BaseLine = (endpoints[0]->node->getPosition()) - (endpoints[1]->node->getPosition());

  //LOG(0, "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << ".");

  BaseLineNormal.Zero();
  NormalCheck.Zero();
  double sign = -1.;
  int i = 0;
  class BoundaryPointSet *node = NULL;
  for (TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
    //LOG(0, "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << ".");
    NormalCheck += runner->second->NormalVector;
    NormalCheck *= sign;
    sign = -sign;
    if (runner->second->NormalVector.NormSquared() > MYEPSILON)
      BaseLineNormal = runner->second->NormalVector;   // yes, copy second on top of first
    else {
      ELOG(0, "Triangle " << *runner->second << " has zero normal vector!");
    }
    node = runner->second->GetThirdEndpoint(this);
    if (node != NULL) {
      //LOG(0, "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << ".");
      helper[i] = (node->node->getPosition()) - BaseLineCenter;
      helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
      //LOG(0, "INFO: Height vector with respect to baseline is " << helper[i] << ".");
      i++;
    } else {
      ELOG(1, "I cannot find third node in triangle, something's wrong.");
      return true;
    }
  }
  //LOG(0, "INFO: BaselineNormal is " << BaseLineNormal << ".");
  if (NormalCheck.NormSquared() < MYEPSILON) {
    LOG(0, "ACCEPT: Normalvectors of both triangles are the same: convex.");
    return true;
  }
  BaseLineNormal.Scale(-1.);
  double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
  return (angle - M_PI);
}

/** Checks whether point is any of the two endpoints this line contains.
 * \param *point point to test
 * \return true - point is of the line, false - is not
 */
bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 2; i++)
    if (point == endpoints[i])
      return true;
  return false;
}
;

/** Returns other endpoint of the line.
 * \param *point other endpoint
 * \return NULL - if endpoint not contained in BoundaryLineSet::lines, or pointer to BoundaryPointSet otherwise
 */
class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
{
  Info FunctionInfo(__func__);
  if (endpoints[0] == point)
    return endpoints[1];
  else if (endpoints[1] == point)
    return endpoints[0];
  else
    return NULL;
}
;

/** Returns other triangle of the line.
 * \param *point other endpoint
 * \return NULL - if triangle not contained in BoundaryLineSet::triangles, or pointer to BoundaryTriangleSet otherwise
 */
class BoundaryTriangleSet *BoundaryLineSet::GetOtherTriangle(const BoundaryTriangleSet * const triangle) const
{
  Info FunctionInfo(__func__);
  if (triangles.size() == 2) {
    for (TriangleMap::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); ++TriangleRunner)
      if (TriangleRunner->second != triangle)
        return TriangleRunner->second;
  }
  return NULL;
}
;

/** output operator for BoundaryLineSet.
 * \param &ost output stream
 * \param &a boundary line
 */
ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
{
  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << " at " << a.endpoints[0]->node->getPosition() << "," << a.endpoints[1]->node->getName() << " at " << a.endpoints[1]->node->getPosition() << "]";
  return ost;
}
;