source: src/Tesselation/unittests/Tesselation_BoundaryTriangleUnitTest.cpp@ c41ce0

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

/*
 * Tesselation_BoundaryTriangleUnitTest.cpp
 *
 *  Created on: Jan 13, 2010
 *      Author: heber
 */

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

using namespace std;

#include <cppunit/CompilerOutputter.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
#include <cppunit/ui/text/TestRunner.h>

#include <cstring>
#include <iostream>

#include "CodePatterns/Log.hpp"
#include "Helpers/defs.hpp"
#include "Atom/TesselPoint.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/VectorSet.hpp"
#include "Tesselation/BoundaryPointSet.hpp"
#include "Tesselation/BoundaryLineSet.hpp"
#include "Tesselation/BoundaryTriangleSet.hpp"
#include "Tesselation/CandidateForTesselation.hpp"

#include "Tesselation_BoundaryTriangleUnitTest.hpp"

#ifdef HAVE_TESTRUNNER
#include "UnitTestMain.hpp"
#endif /*HAVE_TESTRUNNER*/

const double TesselationBoundaryTriangleTest::SPHERERADIUS=2.;

/********************************************** Test classes **************************************/

// Registers the fixture into the 'registry'
CPPUNIT_TEST_SUITE_REGISTRATION( TesselationBoundaryTriangleTest );


void TesselationBoundaryTriangleTest::createTriangle(const std::vector<Vector> &Vectors)
{
  CPPUNIT_ASSERT_EQUAL( (size_t)NDIM, Vectors.size() );

  // create nodes
  for (size_t count = 0; count < NDIM; ++count) {
    tesselpoints[count] = new TesselPoint;
    tesselpoints[count]->setPosition( Vectors[count] );
    tesselpoints[count]->setName(toString(count));
    tesselpoints[count]->setNr(count);
    points[count] = new BoundaryPointSet(tesselpoints[count]);
  }

  // create line
  lines[0] = new BoundaryLineSet(points[0], points[1], 0);
  lines[1] = new BoundaryLineSet(points[1], points[2], 1);
  lines[2] = new BoundaryLineSet(points[0], points[2], 2);

  // create triangle
  triangle = new BoundaryTriangleSet(lines, 0);
  Plane p(Vectors[0], Vectors[1], Vectors[2]);
  triangle->GetNormalVector(p.getNormal());
}

/** This cleanly removes a triangle created via createTriangle() from memory.
 *
 */
void TesselationBoundaryTriangleTest::removeTriangle()
{
  delete(triangle);
  for (int i=0;i<NDIM;++i) {
    // TesselPoint does not delete its vector as it only got a reference
    delete tesselpoints[i];
  }
}

void TesselationBoundaryTriangleTest::setUp()
{
  setVerbosity(5);

  // create nodes
  std::vector<Vector> Vectors;
  Vectors.push_back( Vector(0., 0., 0.) );
  Vectors.push_back( Vector(0., 1., 0.) );
  Vectors.push_back( Vector(1., 0., 0.) );

  // create triangle
  createTriangle(Vectors);
}

void TesselationBoundaryTriangleTest::tearDown()
{
  removeTriangle();
  logger::purgeInstance();
  errorLogger::purgeInstance();
}

/** UnitTest for Tesselation::IsInsideTriangle()
 */
void TesselationBoundaryTriangleTest::IsInsideTriangleTest()
{
  // inside points
  {
    // check each endnode
    for (size_t i=0; i< NDIM; ++i) {
      const Vector testPoint(triangle->endpoints[i]->node->getPosition());
      LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
      CPPUNIT_ASSERT( triangle->IsInsideTriangle( testPoint ) );
    }
  }

  {
    // check points along each BoundaryLine
    for (size_t i=0; i< NDIM; ++i) {
      Vector offset = triangle->endpoints[i%3]->node->getPosition();
      Vector direction = triangle->endpoints[(i+1)%3]->node->getPosition() - offset;
      for (double s = 0.1; s < 1.; s+= 0.1) {
        Vector testPoint = offset + s*direction;
        LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
        CPPUNIT_ASSERT( triangle->IsInsideTriangle( testPoint ) );
      }
    }
  }

  {
    // check central point
    Vector center;
    triangle->GetCenter(center);
    LOG(1, "INFO: Testing whether " << center << " is an inner point.");
    CPPUNIT_ASSERT( triangle->IsInsideTriangle( center ) );
  }

  // outside points
  {
    // check points outside (i.e. those not in xy-plane through origin)
    double n[3];
    const double boundary = 4.;
    const double step = 1.;
    // go through the cube and check each point
    for (n[0] = -boundary; n[0] <= boundary; n[0]+=step)
      for (n[1] = -boundary; n[1] <= boundary; n[1]+=step)
        for (n[2] = -boundary; n[2] <= boundary; n[2]+=step) {
          const Vector testPoint(n[0], n[1], n[2]);
          if (n[2] != 0) {
            LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
            CPPUNIT_ASSERT( !triangle->IsInsideTriangle( testPoint ) );
          }
        }
  }

  {
    // check points within the plane but outside of triangle
    double n[3];
    const double boundary = 4.;
    const double step = 1.;
    n[2] = 0;
    for (n[0] = -boundary; n[0] <= boundary; n[0]+=step)
      for (n[1] = -boundary; n[1] <= boundary; n[1]+=step) {
        const Vector testPoint(n[0], n[1], n[2]);
        if ((n[0] >=0) && (n[1] >=0) && (n[0]<=1) && (n[1]<=1)) {
          if (n[0]+n[1] <= 1) {
            LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
            CPPUNIT_ASSERT( triangle->IsInsideTriangle( testPoint) );
          } else {
            LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
            CPPUNIT_ASSERT( !triangle->IsInsideTriangle( testPoint) );
          }
        } else {
          LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
          CPPUNIT_ASSERT( !triangle->IsInsideTriangle( testPoint) );
        }
      }
  }
}

/** UnitTest for Tesselation::IsInsideTriangle()
 *
 * We test for some specific points that occured in larger examples of the code.
 *
 */
void TesselationBoundaryTriangleTest::IsInsideTriangle_specificTest()
{
  {
    removeTriangle();
    // test is from --create-micelle 200 --radius 30. --position "0,0,0" of sles.data
    // failure is: Intersection (23.1644,24.1867,65.1272) is not inside triangle [659|Na2451,O3652,Na3762].
    const Vector testPoint(1.57318,1.57612,10.9874);
    const Vector testIntersection(23.1644,24.1867,65.1272);
    std::vector<Vector> vectors;
    vectors.push_back( Vector(23.0563,30.4673,73.7555) );
    vectors.push_back( Vector(25.473,25.1512,68.5467) );
    vectors.push_back( Vector(23.1644,24.1867,65.1272) );
    createTriangle(vectors);
    CPPUNIT_ASSERT( triangle->IsInsideTriangle( testIntersection ) );
  }
  {
    removeTriangle();
    // test is from --create-micelle 200 --radius 30. --position "0,0,0" of sles.data
    // failure is: Intersection (20.6787,70.655,71.5657) is not inside triangle [622|Na1197,Na2166,O3366].
    // fix is lower LINALG_MYEPSILON (not std::numeric_limits<doubble>*100 but *1e-4)
    const Vector testPoint(1.57318,14.185,61.2155);
    const Vector testIntersection(20.67867516517798,70.65496977054023,71.56572984946152);
    std::vector<Vector> vectors;
    vectors.push_back( Vector(22.9592,68.7791,77.5907) );
    vectors.push_back( Vector(18.4729,72.0386,68.08839999999999) );
    vectors.push_back( Vector(20.3834,71.0154,70.1443) );
    createTriangle(vectors);
    CPPUNIT_ASSERT( triangle->IsInsideTriangle( testIntersection ) );
  }
  {
    removeTriangle();
    // test is from --create-micelle 200 --radius 30. --position "0,0,0" of sles.data
    // failure is:Intersection (27.56537519896,13.40256646925,6.672946688134) is not inside triangle [702|Na5016,O6388,Na6498].
    // note that the Intersection cannot possibly lie be within the triangle!
    // we test now against correct intersection
    const Vector testIntersection(14.6872,36.204,39.8043);
    std::vector<Vector> vectors;
    vectors.push_back( Vector(10.7513,43.4247,48.4127) );
    vectors.push_back( Vector(13.7119,37.0827,47.4203) );
    vectors.push_back( Vector(14.6872,36.204,39.8043) );
    createTriangle(vectors);
    CPPUNIT_ASSERT( triangle->IsInsideTriangle( testIntersection ) );
  }
}

/** UnitTest for Tesselation::GetClosestPointInsideTriangle()
 *
 * We check whether this function always returns a intersection inside the
 * triangle.
 *
 */
void TesselationBoundaryTriangleTest::GetClosestPointInsideTriangleTest()
{
  Vector TestIntersection;

  {
    // march through a cube mesh on triangle in xy plane
    double n[3];
    const double boundary = 4.;
    const double step = 1.;
    // go through the cube and check each point
    for (n[0] = -boundary; n[0] <= boundary; n[0]+=step)
      for (n[1] = -boundary; n[1] <= boundary; n[1]+=step)
        for (n[2] = -boundary; n[2] <= boundary; n[2]+=step) {
          const Vector testPoint(n[0], n[1], n[2]);
          triangle->GetClosestPointInsideTriangle(testPoint, TestIntersection);
          CPPUNIT_ASSERT( triangle->IsInsideTriangle( TestIntersection ));
        }
  }

  removeTriangle();
  // create better triangle;
  VECTORSET(std::vector) Vectors;
  Vectors.push_back( Vector(0., 0., 0.) );
  Vectors.push_back( Vector(0., 1., 0.) );
  Vectors.push_back( Vector(1., 0., 0.) );
  RealSpaceMatrix M;
  M.setRotation(M_PI/3., M_PI/4., 2.*M_PI/3.);
  Vectors.transform(M);
  createTriangle(Vectors);

  {
    // march through a cube mesh on rotated triangle
    double n[3];
    const double boundary = 4.;
    const double step = 1.;
    // go through the cube and check each point
    for (n[0] = -boundary; n[0] <= boundary; n[0]+=step)
      for (n[1] = -boundary; n[1] <= boundary; n[1]+=step)
        for (n[2] = -boundary; n[2] <= boundary; n[2]+=step) {
          const Vector testPoint(n[0], n[1], n[2]);
          triangle->GetClosestPointInsideTriangle(testPoint, TestIntersection);
          CPPUNIT_ASSERT( triangle->IsInsideTriangle( TestIntersection ));
        }
  }
}

/** UnitTest for Tesselation::GetClosestPointInsideTriangle()
 *
 * We test for some specific points that occured in larger examples of the code.
 *
 */
void TesselationBoundaryTriangleTest::GetClosestPointInsideTriangle_specificTest()
{
  {
    removeTriangle();
    // test is from --create-micelle 200 --radius 30. --position "0,0,0" of sles.data
    // failure is:Intersection (27.56537519896,13.40256646925,6.672946688134) is not inside triangle [702|Na5016,O6388,Na6498].
    // note that the Intersection cannot possibly lie be within the triangle
    // however, it is on its plane (off only by 2.7e-12)
    // testPoint2 is corrected version
    const Vector testPoint(27.56537519896,13.40256646925,6.672946688134);
    const Vector testPoint2(14.6872,36.204,39.8043);
    Vector testIntersection;
    std::vector<Vector> vectors;
    vectors.push_back( Vector(10.7513,43.4247,48.4127) );
    vectors.push_back( Vector(13.7119,37.0827,47.4203) );
    vectors.push_back( Vector(14.6872,36.204,39.8043) );
    createTriangle(vectors);
    triangle->GetClosestPointInsideTriangle(testPoint, testIntersection);
    CPPUNIT_ASSERT ( testPoint != testIntersection );
    CPPUNIT_ASSERT ( testPoint2 == testIntersection );
  }
}

/** UnitTest for Tesselation::IsInnerPoint()
 */
void TesselationBoundaryTriangleTest::GetClosestPointOnPlaneTest()
{
  Vector TestIntersection;
  Vector Point;

  // simple test on y line
  Point = Vector(-1.,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(-4.,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( 16., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on x line
  Point = Vector(0.5,-1.,0.);
  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(0.5,-6.,0.);
  CPPUNIT_ASSERT_EQUAL( 36., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on slanted line
  Point = Vector(1.,1.,0.);
  CPPUNIT_ASSERT_EQUAL( 0.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(5.,5.,0.);
  CPPUNIT_ASSERT_EQUAL( 40.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on first node
  Point = Vector(-1.,-1.,0.);
  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on second node
  Point = Vector(0.,2.,0.);
  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,1.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on third node
  Point = Vector(2.,0.,0.);
  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(1.,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
};

/** UnitTest for Tesselation::IsInnerPoint()
 */
void TesselationBoundaryTriangleTest::GetClosestPointOffPlaneTest()
{
  Vector TestIntersection;
  Vector Point;

  // straight down/up
  Point = Vector(1./3.,1./3.,+5.);
  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(1./3.,1./3.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(1./3.,1./3.,-5.);
  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(1./3.,1./3.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on y line
  Point = Vector(-1.,0.5,+2.);
  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(-1.,0.5,-3.);
  CPPUNIT_ASSERT_EQUAL( 10., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on x line
  Point = Vector(0.5,-1.,+1.);
  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(0.5,-1.,-2.);
  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on slanted line
  Point = Vector(1.,1.,+3.);
  CPPUNIT_ASSERT_EQUAL( 9.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
  Point = Vector(1.,1.,-4.);
  CPPUNIT_ASSERT_EQUAL( 16.5, triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.5,0.5,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on first node
  Point = Vector(-1.,-1.,5.);
  CPPUNIT_ASSERT_EQUAL( 27., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on second node
  Point = Vector(0.,2.,5.);
  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(0.,1.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );

  // simple test on third node
  Point = Vector(2.,0.,5.);
  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(Point, TestIntersection) );
  Point = Vector(1.,0.,0.);
  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
};