source: src/Tesselation/unittests/Tesselation_InsideOutsideUnitTest.cpp@ b9d387

/*
 * 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_InsideOutsideUnitTest.cpp
 *
 *  Created on: Dec 28, 2009
 *      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 "Atom/TesselPoint.hpp"
#include "CodePatterns/Log.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "Tesselation/BoundaryLineSet.hpp"
#include "Tesselation/BoundaryTriangleSet.hpp"
#include "Tesselation/CandidateForTesselation.hpp"
#include "Tesselation/tesselation.hpp"

#include "Tesselation_InsideOutsideUnitTest.hpp"

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

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

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

/** Creates TesselPoints out of given Vector's.
 *
 * @param Vectors given vector with Vector's as positions for TesselPoint's
 */
void TesselationInOutsideTest::prepareCorners(const std::vector<Vector> &Vectors)
{
  class TesselPoint *Walker;
  size_t index = 0;
  for (std::vector<Vector>::const_iterator iter = Vectors.begin();
      iter != Vectors.end(); ++iter) {
    Walker = new TesselPoint;
    Walker->setPosition( *iter );
    Walker->setNr(index++);
    Walker->setName(toString(index)); // yes, name is one higher than index
    Corners.push_back(Walker);
  }
}

/** Prepares Vector's in such a way that the mesh forms a cube.
 *
 * @param vectors vector of Vector's to work on
 * @param factor factor to expand mesh
 * @param offset offset to translate mesh
 * @param return vector for concatenation
 */
std::vector<Vector> TesselationInOutsideTest::translateAndexpand(
    VECTORSET(std::vector) Vectors, const double factor, const Vector &offset) const
{
  RealSpaceMatrix M;
  M.setIdentity();
  M *= factor;
  Vectors.transform(M);
  Vectors.translate(offset);
  return Vectors;
}

/** Prepares Vector's in such a way that the mesh forms a cube.
 *
 */
std::vector<Vector> TesselationInOutsideTest::setupCube() const
{
  std::vector<Vector> returnVectors;
  returnVectors.push_back( Vector(0., 0., 0.) );
  returnVectors.push_back( Vector(0., 1., 0.) );
  returnVectors.push_back( Vector(1., 0., 0.) );
  returnVectors.push_back( Vector(1., 1., 0.) );
  returnVectors.push_back( Vector(0., 0., 1.) );
  returnVectors.push_back( Vector(0., 1., 1.) );
  returnVectors.push_back( Vector(1., 0., 1.) );
  returnVectors.push_back( Vector(1., 1., 1.) );
  return returnVectors;
}

/** Prepares Vector's in such a way that the mesh forms a sphere.
 *
 */
std::vector<Vector> TesselationInOutsideTest::setupSphere() const
{
  std::vector<Vector> returnVectors;
  // start with a cube
  returnVectors.push_back( Vector(0., 0., 0.) );
  returnVectors.push_back( Vector(0., 1., 0.) );
  returnVectors.push_back( Vector(1., 0., 0.) );
  returnVectors.push_back( Vector(1., 1., 0.) );
  returnVectors.push_back( Vector(0., 0., 1.) );
  returnVectors.push_back( Vector(0., 1., 1.) );
  returnVectors.push_back( Vector(1., 0., 1.) );
  returnVectors.push_back( Vector(1., 1., 1.) );
  // then add a point slightly above each face
  returnVectors.push_back( Vector(0.5, 0.5, -0.4142136) );
  returnVectors.push_back( Vector(0.5, 0.5, 1.4142136) );
  returnVectors.push_back( Vector(0.5, -0.4142136, 0.5) );
  returnVectors.push_back( Vector(0.5, 1.4142136, 0.5) );
  returnVectors.push_back( Vector(-0.4142136, 0.5, 0.5) );
  returnVectors.push_back( Vector(1.4142136, 0.5, 0.5) );
  return returnVectors;
}

/** Prepares Vector's in such a way that the mesh forms a nonconvex form.
 *
 */
std::vector<Vector> TesselationInOutsideTest::setupNonConvex() const
{
  std::vector<Vector> returnVectors;
  // make an along the x-axis elongated cuboid
  returnVectors.push_back( Vector(0., 0., 0.) );
  returnVectors.push_back( Vector(0., 1., 0.) );
  returnVectors.push_back( Vector(2., 0., 0.) );
  returnVectors.push_back( Vector(2., 1., 0.) );
  returnVectors.push_back( Vector(0., 0., 1.) );
  returnVectors.push_back( Vector(0., 1., 1.) );
  returnVectors.push_back( Vector(2., 0., 1.) );
  returnVectors.push_back( Vector(2., 1., 1.) );
  // add two lowered points in the middle of the elongation
  returnVectors.push_back( Vector(1., 0.5, 0.) );
  returnVectors.push_back( Vector(1., 0.5, 1.) );
  returnVectors.push_back( Vector(1., 0., 0.) );
  returnVectors.push_back( Vector(1., 0., 1.) );
  return returnVectors;
}

/** Creates the tesselation out of current setup in \a Corners.
 *
 */
void TesselationInOutsideTest::prepareTesselation(const double SPHERERADIUS)
{
  // create LinkedCell
  PointCloudAdaptor< TesselPointSTLList > cloud(&Corners, "TesselPointSTLList");
  LinkedList = new LinkedCell_deprecated(cloud, 2.*SPHERERADIUS);

  // create tesselation
  TesselStruct = new Tesselation;
  (*TesselStruct)(cloud, SPHERERADIUS);
}

/** Removes any currently present tesselation.
 *
 */
void TesselationInOutsideTest::removeTesselation()
{
  delete LinkedList;
  delete TesselStruct;
  for (TesselPointSTLList::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
    delete *Runner;

  LinkedList = NULL;
  TesselStruct = NULL;
}

void TesselationInOutsideTest::setUp()
{
  setVerbosity(6);
}

void TesselationInOutsideTest::tearDown()
{
  removeTesselation();
  Corners.clear();
  logger::purgeInstance();
  errorLogger::purgeInstance();
}

/** UnitTest for Tesselation::IsInnerPoint() with a cube mesh
 */
void TesselationInOutsideTest::IsInnerPointCubeTest()
{
  {
    const double SPHERERADIUS = 2.;
    prepareCorners( translateAndexpand( setupCube(), 1., Vector(0.,0.,0.)) );
    prepareTesselation(SPHERERADIUS);
    CPPUNIT_ASSERT( TesselStruct != NULL );
    PointCloudAdaptor< TesselPointSTLList > cloud(&Corners, "TesselPointSTLList");
    //TesselStruct->Output("cube", cloud);
    double n[3];
    const double boundary = 2.;
    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[0] >= 0.) && (n[1] >= 0.) && (n[2] >= 0.)) && ((n[0] <= 1.) && (n[1] <= 1.) && (n[2] <= 1.)))
            CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
          else
            CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
        }
  }
}

/** UnitTest for Tesselation::IsInnerPoint() with a sphere mesh
 */
void TesselationInOutsideTest::IsInnerPointSphereTest()
{
  {
    // with radius 2, we see the opposite point which causes FindStartingTriangle() to fail
    const double SPHERERADIUS = 1.;
    const Vector offset(0.5,0.5,0.5);
    const double factor = 1.;
    prepareCorners( translateAndexpand( setupSphere(), factor, offset) );
    prepareTesselation(SPHERERADIUS);
    CPPUNIT_ASSERT( TesselStruct != NULL );
    double n[3];
    const double boundary = 2.;
    const double step = 1.;
    const std::vector<bool> test(1, true);
    const Vector SphereCenter( Vector(0.5,0.5,0.5) + offset );
    const double SphereRadius = (Vector(0.,0.,0.)+offset - SphereCenter).Norm();
    // 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]);
          // radius is actually 2. but we simply avoid the corners in this test
          if ( testPoint.DistanceSquared(SphereCenter) <= factor*SphereRadius ) {
            LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
            CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
          } else {
            LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
            CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
          }
        }
  }
}

/** UnitTest for Tesselation::IsInnerPoint() with a non-convex mesh
 */
void TesselationInOutsideTest::IsInnerPointNonConvexTest()
{
  {
    const double SPHERERADIUS = 1.;
    prepareCorners( translateAndexpand( setupNonConvex(), 1., Vector(0.,0.,0.)) );
    prepareTesselation(SPHERERADIUS);
    CPPUNIT_ASSERT( TesselStruct != NULL );
    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[0] >= 0.) && (n[1] >= 0.) && (n[2] >= 0.)) && ((n[0] <= 2.) && (n[1] <= .5) && (n[2] <= 1.))) {
            // (y) lower half of elongated block
            LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
            CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
          } else if ( ((n[0] >= 0.) && (n[1] >= 0.5) && (n[2] >= 0.)) && ((n[0] <= 2.) && (n[1] <= 1.) && (n[2] <= 1.))) {
            // (y) upper half of elongated block
            if ( ((n[0] >= 0.) && (n[1] >= 0.5) && (n[2] >= 0.)) && ((n[0] <= 1.) && (n[1] <= 1.) && (n[2] <= 1.))) {
              // (x) left side
              if (n[0]+n[1] <= 1.) {
                LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
                CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
              } else {
                LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
                CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
              }
            } else {
              // (x) right side)
              if ((1.-n[0])+n[1] <= 1.) {
                LOG(1, "INFO: Testing whether " << testPoint << " is an inner point.");
                CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
              } else {
                LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
                CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
              }
            }
          } else {
            LOG(1, "INFO: Testing whether " << testPoint << " is not an inner point.");
            CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(testPoint, LinkedList) );
          }
        }
  }
}