source: ThirdParty/LinearAlgebra/src/unittests/LineUnitTest.cpp@ e7ff10

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

/*
 * LineUnittest.cpp
 *
 *  Created on: May 27, 2010
 *      Author: crueger
 */

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


#include "Exceptions.hpp"
#include "Vector.hpp"

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

#include <iostream>
#include <cmath>

#include "LineUnitTest.hpp"

#include "RealSpaceMatrix.hpp"

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

CPPUNIT_TEST_SUITE_REGISTRATION( LineUnittest );

void LineUnittest::setUp(){
  // three lines along the axes
  la1 = new Line(zeroVec,unitVec[0]);
  la2 = new Line(zeroVec,unitVec[1]);
  la3 = new Line(zeroVec,unitVec[2]);

  // the lines along the planes defined by two coordinate axes
  lp1 = new Line(unitVec[0],unitVec[0]-unitVec[1]);
  lp2 = new Line(unitVec[1],unitVec[1]-unitVec[2]);
  lp3 = new Line(unitVec[2],unitVec[2]-unitVec[0]);
}
void LineUnittest::tearDown(){
  delete la1;
  delete la2;
  delete la3;

  delete lp1;
  delete lp2;
  delete lp3;
}

void LineUnittest::constructionErrorTest(){
  // test some constructions

  // direction+origin should never fail
  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,unitVec[0]));
  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,unitVec[1]));
  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,unitVec[2]));

  // two points fails if both points are the same
  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[0],unitVec[1]));
  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[1],unitVec[2]));
  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[2],unitVec[0]));
  // for zerovectors
  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[0],zeroVec));
  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[1],zeroVec));
  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[2],zeroVec));
  // now we pass two times the same point
  CPPUNIT_ASSERT_THROW(makeLineThrough(zeroVec,zeroVec),LinearDependenceException);
  CPPUNIT_ASSERT_THROW(makeLineThrough(unitVec[0],unitVec[0]),LinearDependenceException);
  CPPUNIT_ASSERT_THROW(makeLineThrough(unitVec[1],unitVec[1]),LinearDependenceException);
  CPPUNIT_ASSERT_THROW(makeLineThrough(unitVec[2],unitVec[2]),LinearDependenceException);

}

bool testDirection(const Vector &dir1,const Vector &dir2){
  return (dir1==dir2) || (dir1==-1*dir2);
}

void LineUnittest::constructionResultTest(){
  // test all directions
  CPPUNIT_ASSERT(testDirection(la1->getDirection(),unitVec[0]));
  CPPUNIT_ASSERT(testDirection(la2->getDirection(),unitVec[1]));
  CPPUNIT_ASSERT(testDirection(la3->getDirection(),unitVec[2]));

  // test origins
  CPPUNIT_ASSERT_EQUAL(la1->getOrigin(),zeroVec);
  CPPUNIT_ASSERT_EQUAL(la2->getOrigin(),zeroVec);
  CPPUNIT_ASSERT_EQUAL(la2->getOrigin(),zeroVec);

  // test if desired points are on the lines
  CPPUNIT_ASSERT(la1->isContained(zeroVec));
  CPPUNIT_ASSERT(la2->isContained(zeroVec));
  CPPUNIT_ASSERT(la3->isContained(zeroVec));

  CPPUNIT_ASSERT(la1->isContained(unitVec[0]));
  CPPUNIT_ASSERT(la2->isContained(unitVec[1]));
  CPPUNIT_ASSERT(la3->isContained(unitVec[2]));

  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]));
  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]));
  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]));

  CPPUNIT_ASSERT(lp1->isContained(unitVec[1]));
  CPPUNIT_ASSERT(lp2->isContained(unitVec[2]));
  CPPUNIT_ASSERT(lp3->isContained(unitVec[0]));
}

void LineUnittest::isContainedTest(){
  // Zerovector on the axes lines
  CPPUNIT_ASSERT(la1->isContained(zeroVec));
  CPPUNIT_ASSERT(la2->isContained(zeroVec));
  CPPUNIT_ASSERT(la3->isContained(zeroVec));

  // multiples of the second support vector
  CPPUNIT_ASSERT(la1->isContained(unitVec[0]));
  CPPUNIT_ASSERT(la2->isContained(unitVec[1]));
  CPPUNIT_ASSERT(la3->isContained(unitVec[2]));

  CPPUNIT_ASSERT(la1->isContained(2*unitVec[0]));
  CPPUNIT_ASSERT(la2->isContained(2*unitVec[1]));
  CPPUNIT_ASSERT(la3->isContained(2*unitVec[2]));

  CPPUNIT_ASSERT(la1->isContained(3*unitVec[0]));
  CPPUNIT_ASSERT(la2->isContained(3*unitVec[1]));
  CPPUNIT_ASSERT(la3->isContained(3*unitVec[2]));

  // negative multiples
  CPPUNIT_ASSERT(la1->isContained(-1*unitVec[0]));
  CPPUNIT_ASSERT(la2->isContained(-1*unitVec[1]));
  CPPUNIT_ASSERT(la3->isContained(-1*unitVec[2]));

  CPPUNIT_ASSERT(la1->isContained(-2*unitVec[0]));
  CPPUNIT_ASSERT(la2->isContained(-2*unitVec[1]));
  CPPUNIT_ASSERT(la3->isContained(-2*unitVec[2]));

  // points that should not be on the lines
  CPPUNIT_ASSERT(!la1->isContained(unitVec[1]));
  CPPUNIT_ASSERT(!la2->isContained(unitVec[2]));
  CPPUNIT_ASSERT(!la3->isContained(unitVec[0]));

  CPPUNIT_ASSERT(!la1->isContained(2*unitVec[1]));
  CPPUNIT_ASSERT(!la2->isContained(2*unitVec[2]));
  CPPUNIT_ASSERT(!la3->isContained(2*unitVec[0]));

  CPPUNIT_ASSERT(!la1->isContained(-1*unitVec[1]));
  CPPUNIT_ASSERT(!la2->isContained(-1*unitVec[2]));
  CPPUNIT_ASSERT(!la3->isContained(-1*unitVec[0]));

  // For the plane lines
  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]));
  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]));
  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]));

  CPPUNIT_ASSERT(lp1->isContained(unitVec[1]));
  CPPUNIT_ASSERT(lp2->isContained(unitVec[2]));
  CPPUNIT_ASSERT(lp3->isContained(unitVec[0]));

  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]+2*(unitVec[0]-unitVec[1])));
  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]+2*(unitVec[1]-unitVec[2])));
  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]+2*(unitVec[2]-unitVec[0])));

  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]-2*(unitVec[0]-unitVec[1])));
  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]-2*(unitVec[1]-unitVec[2])));
  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]-2*(unitVec[2]-unitVec[0])));
}

void LineUnittest::intersectionTest(){
  Vector fixture;

  // intersection of the axis lines
  fixture = la1->getIntersection(*la2);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = la2->getIntersection(*la3);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = la3->getIntersection(*la1);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);

  // axes and plane lines
  fixture = la1->getIntersection(*lp1);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = la2->getIntersection(*lp2);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = la3->getIntersection(*lp3);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  fixture = la1->getIntersection(*lp3);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = la2->getIntersection(*lp1);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = la3->getIntersection(*lp2);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  // two plane lines
  fixture = lp1->getIntersection(*lp2);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = lp2->getIntersection(*lp3);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
  fixture = lp3->getIntersection(*lp1);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);

  // When we have two times the same line, we check if the point is on the line
  fixture = la1->getIntersection(*la1);
  CPPUNIT_ASSERT(la1->isContained(fixture));
  fixture = la2->getIntersection(*la2);
  CPPUNIT_ASSERT(la2->isContained(fixture));
  fixture = la3->getIntersection(*la3);
  CPPUNIT_ASSERT(la3->isContained(fixture));

  fixture = lp1->getIntersection(*lp1);
  CPPUNIT_ASSERT(lp1->isContained(fixture));
  fixture = lp2->getIntersection(*lp2);
  CPPUNIT_ASSERT(lp2->isContained(fixture));
  fixture = lp3->getIntersection(*lp3);
  CPPUNIT_ASSERT(lp3->isContained(fixture));

  // lines that are askew should produce an Error
  CPPUNIT_ASSERT_THROW(lp1->getIntersection(*la3),SkewException);
  CPPUNIT_ASSERT_THROW(lp2->getIntersection(*la1),SkewException);
  CPPUNIT_ASSERT_THROW(lp3->getIntersection(*la2),SkewException);

  CPPUNIT_ASSERT_THROW(la1->getIntersection(*lp2),SkewException);
  CPPUNIT_ASSERT_THROW(la2->getIntersection(*lp3),SkewException);
  CPPUNIT_ASSERT_THROW(la3->getIntersection(*lp1),SkewException);
}

void LineUnittest::rotationTest(){
  Vector fixture;

  // rotate zero Vector along the axes lines by various degrees
  fixture = la1->rotateVector(zeroVec,1.);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = la2->rotateVector(zeroVec,1.);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = la3->rotateVector(zeroVec,1.);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);

  fixture = la1->rotateVector(zeroVec,2.);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = la2->rotateVector(zeroVec,2.);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = la3->rotateVector(zeroVec,2.);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);

  // rotate vectors on the axis around their lines
  fixture = la1->rotateVector(unitVec[0],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = la2->rotateVector(unitVec[1],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = la3->rotateVector(unitVec[2],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  fixture = la1->rotateVector(unitVec[0],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = la2->rotateVector(unitVec[1],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = la3->rotateVector(unitVec[2],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  // more vectors on the axis
  fixture = la1->rotateVector(2*unitVec[0],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[0]);
  fixture = la2->rotateVector(2*unitVec[1],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[1]);
  fixture = la3->rotateVector(2*unitVec[2],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[2]);

  fixture = la1->rotateVector(2*unitVec[0],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[0]);
  fixture = la2->rotateVector(2*unitVec[1],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[1]);
  fixture = la3->rotateVector(2*unitVec[2],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[2]);

  // negative factors
  fixture = la1->rotateVector(-1*unitVec[0],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
  fixture = la2->rotateVector(-1*unitVec[1],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = la3->rotateVector(-1*unitVec[2],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);

  fixture = la1->rotateVector(-1*unitVec[0],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
  fixture = la2->rotateVector(-1*unitVec[1],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = la3->rotateVector(-1*unitVec[2],2.);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);


  // now the real rotations
  // unitVec[1] around unitVec[0]
  fixture = la1->rotateVector(unitVec[1],0);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = la1->rotateVector(unitVec[1],1./2.*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
  fixture = la1->rotateVector(unitVec[1],M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = la1->rotateVector(unitVec[1],2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);

  // unitVec[2] around unitVec[1]
  fixture = la2->rotateVector(unitVec[2],0);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
  fixture = la2->rotateVector(unitVec[2],1./2.*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
  fixture = la2->rotateVector(unitVec[2],M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
  fixture = la2->rotateVector(unitVec[2],2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  // unitVec[0] around unitVec[2]
  fixture = la3->rotateVector(unitVec[0],0);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = la3->rotateVector(unitVec[0],1./2.*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = la3->rotateVector(unitVec[0],M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
  fixture = la3->rotateVector(unitVec[0],2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);


  // and some rotation around the plane lines

  // Vectors on the line
  fixture = lp1->rotateVector(unitVec[0],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = lp1->rotateVector(unitVec[1],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);

  fixture = lp2->rotateVector(unitVec[1],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = lp2->rotateVector(unitVec[2],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  fixture = lp3->rotateVector(unitVec[2],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
  fixture = lp3->rotateVector(unitVec[0],1.);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);

  // the real stuff
  fixture = lp1->rotateVector(zeroVec,M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,Vector(1,1,0));
  fixture = lp2->rotateVector(zeroVec,M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,Vector(0,1,1));
  fixture = lp3->rotateVector(zeroVec,M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,Vector(1,0,1));

  fixture = lp1->rotateVector(zeroVec,2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = lp2->rotateVector(zeroVec,2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
  fixture = lp3->rotateVector(zeroVec,2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
}

void LineUnittest::getRotationMatrixTest()
{
  RealSpaceMatrix M;
  Vector fixture;

  // check getRotationMatrix;
  M = la1->getRotationMatrix(1.);
  fixture = M * unitVec[0];
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  M = la2->getRotationMatrix(1.);
  fixture = M * unitVec[1];
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  M = la3->getRotationMatrix(1.);
  fixture = M * unitVec[2];
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  // unitVec[1] around unitVec[0]
  fixture = la1->rotateVector(unitVec[1],0);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
  fixture = la1->rotateVector(unitVec[1],1./2.*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
  fixture = la1->rotateVector(unitVec[1],M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = la1->rotateVector(unitVec[1],2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);

  // unitVec[2] around unitVec[1]
  fixture = la2->rotateVector(unitVec[2],0);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
  fixture = la2->rotateVector(unitVec[2],1./2.*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
  fixture = la2->rotateVector(unitVec[2],M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
  fixture = la2->rotateVector(unitVec[2],2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);

  // unitVec[0] around unitVec[2]
  fixture = la3->rotateVector(unitVec[0],0);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
  fixture = la3->rotateVector(unitVec[0],1./2.*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
  fixture = la3->rotateVector(unitVec[0],M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
  fixture = la3->rotateVector(unitVec[0],2*M_PI);
  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
}

void LineUnittest::sphereIntersectionTest(){
  {
    std::vector<Vector> res = la1->getSphereIntersections();
    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    CPPUNIT_ASSERT(testDirection(res[0],unitVec[0]));
    CPPUNIT_ASSERT(testDirection(res[1],unitVec[0]));
    CPPUNIT_ASSERT(res[0]!=res[1]);
  }

  {
    std::vector<Vector> res = la2->getSphereIntersections();
    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    CPPUNIT_ASSERT(testDirection(res[0],unitVec[1]));
    CPPUNIT_ASSERT(testDirection(res[1],unitVec[1]));
    CPPUNIT_ASSERT(res[0]!=res[1]);
  }

  {
    std::vector<Vector> res = la3->getSphereIntersections();
    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    CPPUNIT_ASSERT(testDirection(res[0],unitVec[2]));
    CPPUNIT_ASSERT(testDirection(res[1],unitVec[2]));
    CPPUNIT_ASSERT(res[0]!=res[1]);
  }

  {
    std::vector<Vector> res = lp1->getSphereIntersections();
    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    CPPUNIT_ASSERT((res[0]==unitVec[0]) || (res[0]==unitVec[1]));
    CPPUNIT_ASSERT((res[1]==unitVec[0]) || (res[1]==unitVec[1]));
    CPPUNIT_ASSERT(res[0]!=res[1]);
  }

  {
    std::vector<Vector> res = lp2->getSphereIntersections();
    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    CPPUNIT_ASSERT((res[0]==unitVec[1]) || (res[0]==unitVec[2]));
    CPPUNIT_ASSERT((res[1]==unitVec[1]) || (res[1]==unitVec[2]));
    CPPUNIT_ASSERT(res[0]!=res[1]);
  }

  {
    std::vector<Vector> res = lp3->getSphereIntersections();
    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    CPPUNIT_ASSERT((res[0]==unitVec[2]) || (res[0]==unitVec[0]));
    CPPUNIT_ASSERT((res[1]==unitVec[2]) || (res[1]==unitVec[0]));
    CPPUNIT_ASSERT(res[0]!=res[1]);
  }
}