Changes in / [fd4905:6c438f]

Location:

src

Files:

• Box.cpp (modified) (2 diffs)
• Box.hpp (modified) (2 diffs)
• Descriptors/AtomShapeDescriptor.cpp (modified) (2 diffs)
• Descriptors/AtomShapeDescriptor.hpp (modified) (2 diffs)
• Descriptors/AtomShapeDescriptor_impl.hpp (modified) (2 diffs)
• Exceptions/Makefile.am (modified) (2 diffs)
• Exceptions/NotOnSurfaceException.cpp (added)
• Exceptions/NotOnSurfaceException.hpp (added)
• Exceptions/ShapeException.cpp (added)
• Exceptions/ShapeException.hpp (added)
• Helpers/helpers.cpp (modified) (1 diff)
• Helpers/helpers.hpp (modified) (2 diffs)
• LinearAlgebra/Line.cpp (modified) (4 diffs)
• LinearAlgebra/Line.hpp (modified) (4 diffs)
• LinearAlgebra/LineSegment.cpp (added)
• LinearAlgebra/LineSegment.hpp (added)
• LinearAlgebra/LineSegmentSet.cpp (added)
• LinearAlgebra/LineSegmentSet.hpp (added)
• LinearAlgebra/Makefile.am (modified) (2 diffs)
• LinearAlgebra/Matrix.cpp (modified) (2 diffs)
• LinearAlgebra/Plane.cpp (modified) (4 diffs)
• LinearAlgebra/Plane.hpp (modified) (4 diffs)
• LinearAlgebra/Vector.cpp (modified) (2 diffs)
• LinearAlgebra/Vector.hpp (modified) (2 diffs)
• Makefile.am (modified) (2 diffs)
• Shapes/BaseShapes.cpp (modified) (4 diffs)
• Shapes/BaseShapes.hpp (modified) (1 diff)
• Shapes/BaseShapes_impl.hpp (modified) (2 diffs)
• Shapes/Shape.cpp (modified) (7 diffs)
• Shapes/Shape.hpp (modified) (4 diffs)
• Shapes/ShapeOps.cpp (modified) (8 diffs)
• Shapes/ShapeOps_impl.hpp (modified) (1 diff)
• Shapes/Shape_impl.hpp (modified) (7 diffs)
• unittests/BaseShapesUnittest.cpp (deleted)
• unittests/BaseShapesUnittest.hpp (deleted)
• unittests/LineUnittest.cpp (modified) (14 diffs)
• unittests/Makefile.am (modified) (4 diffs)
• unittests/MatrixUnittest.cpp (modified) (5 diffs)
• unittests/PlaneUnittest.cpp (modified) (6 diffs)
• unittests/ShapeUnittest.cpp (modified) (3 diffs)
• unittests/ShapeUnittest.hpp (modified) (2 diffs)

src/Box.cpp

@@ -29 +29 @@
 #include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/Plane.hpp"
+#include "Shapes/BaseShapes.hpp"
+#include "Shapes/ShapeOps.hpp"
 
 #include "Helpers/Assert.hpp"
@@ -218 +220 @@
 }
 
+Shape Box::getShape() const{
+  return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
+}
+
 const Box::Conditions_t Box::getConditions(){
   return conditions;

src/Box.hpp

@@ -11 +11 @@
 class Matrix;
 class Vector;
+class Shape;
 class Plane;
 
@@ -99 +100 @@
   double periodicDistance(const Vector &point1,const Vector &point2) const;
 
+  Shape getShape() const;
   const Conditions_t getConditions();
   void setCondition(int,BoundaryCondition_t);

src/Descriptors/AtomShapeDescriptor.cpp

@@ -9 +9 @@
  * AtomShapeDescriptor.cpp
  *
- *  Created on: Aug 9, 2010
- *      Author: heber
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
  */
 
@@ -24 +24 @@
 
 #include "atom.hpp"
-#include "Patterns/ObservedContainer_impl.hpp"
-#include "Shapes/Shape.hpp"
 
-using namespace std;
-
-
-AtomShapeDescriptor_impl::AtomShapeDescriptor_impl(Shape _shape) :
+AtomShapeDescriptor_impl::AtomShapeDescriptor_impl(const Shape &_shape) :
   shape(_shape)
 {}
 
-AtomShapeDescriptor_impl::~AtomShapeDescriptor_impl()
-{}
+AtomShapeDescriptor_impl::~AtomShapeDescriptor_impl(){}
 
-bool AtomShapeDescriptor_impl::predicate(std::pair<atomId_t,atom*> atom) {
+bool AtomShapeDescriptor_impl::predicate(std::pair<atomId_t,atom*> atom){
   return shape.isInside(atom.second->getPosition());
 }
 
-AtomDescriptor AtomByShape(Shape &_shape){
-  return AtomDescriptor(AtomDescriptor::impl_ptr(new AtomShapeDescriptor_impl(_shape)));
+
+AtomDescriptor AtomByShape(const Shape &shape){
+  return AtomDescriptor(AtomDescriptor::impl_ptr(new AtomShapeDescriptor_impl(shape)));
 }

src/Descriptors/AtomShapeDescriptor.hpp

@@ -2 +2 @@
  * AtomShapeDescriptor.hpp
  *
- *  Created on: Aug 9, 2010
- *      Author: heber
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
  */
 
@@ -14 +14 @@
 class Shape;
 
-AtomDescriptor AtomByShape(Shape &_shape);
+AtomDescriptor AtomByShape(const Shape &shape);
 
 #endif /* ATOMSHAPEDESCRIPTOR_HPP_ */

src/Descriptors/AtomShapeDescriptor_impl.hpp

@@ -2 +2 @@
  * AtomShapeDescriptor_impl.hpp
  *
- *  Created on: Aug 9, 2010
- *      Author: heber
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
  */
 
-#ifndef ATOMSHAPEDESCRIPTOR_IMPL_HPP
-#define ATOMSHAPEDESCRIPTOR_IMPL_HPP
+#ifndef ATOMSHAPEDESCRIPTOR_IMPL_HPP_
+#define ATOMSHAPEDESCRIPTOR_IMPL_HPP_
 
 #include "Descriptors/AtomDescriptor_impl.hpp"
+
 #include "Shapes/Shape.hpp"
 
@@ -15 +16 @@
 {
 public:
-  AtomShapeDescriptor_impl(Shape _shape);
-  virtual ~AtomShapeDescriptor_impl();
+  AtomShapeDescriptor_impl(const Shape&);
+  virtual  ~AtomShapeDescriptor_impl();
 
   bool predicate(std::pair<atomId_t,atom*> atom);
-
 private:
   Shape shape;
 };
 
-#endif //ATOMSHAPEDESCRIPTOR_IMPL_HPP
+#endif /* ATOMSHAPEDESCRIPTOR_IMPL_HPP_ */

src/Exceptions/Makefile.am

@@ -15 +15 @@
   MissingValueException.cpp \
   NotInvertibleException.cpp \
+  NotOnSurfaceException.cpp \
+  ShapeException.cpp \
   ParseError.cpp \
   SkewException.cpp \
   ZeroVectorException.cpp
+        
                                   
 EXCEPTIONHEADER = \
@@ -26 +29 @@
   MissingValueException.hpp \
   NotInvertibleException.hpp \
+  NotInvertibleException.hpp \
+  NotOnSurfaceException.hpp \
   ParseError.hpp \
   SkewException.hpp \

src/Helpers/helpers.cpp

@@ -153 +153 @@
   exit(255);
 }
+
+sign_t sign(double value){
+  if(fabs(value)<MYEPSILON){
+    return Zero;
+  }
+  if(value<0)
+    return Minus;
+  else
+    return Plus;
+}

src/Helpers/helpers.hpp

@@ -43 +43 @@
     const char *file, const int line);
   //__attribute__ ((__return__));
+
+typedef enum {
+  Minus = -1,
+  Zero = 0,
+  Plus = +1
+} sign_t;
 
 double ask_value(const char *text);
@@ -52 +58 @@
 int CompareDoubles (const void * a, const void * b);
 void performCriticalExit();
+sign_t sign(double value);
 
 /********************************************** helpful template functions *********************************/

src/LinearAlgebra/Line.cpp

@@ -50 +50 @@
 Line::~Line()
 {}
+
+Line &Line::operator=(const Line& rhs){
+  if(this!=&rhs){
+    origin.reset(new Vector(*rhs.origin));
+    direction.reset(new Vector(*rhs.direction));
+  }
+  return *this;
+}
 
 
@@ -224 +232 @@
 }
 
+Line Line::rotateLine(const Line &rhs, double alpha) const{
+  Vector lineOrigin = rotateVector(rhs.getOrigin(),alpha);
+  Vector helper = rhs.getDirection();
+  // rotate the direction without considering the ofset
+  pair<Vector,Vector> parts = helper.partition(*direction);
+  Vector lineDirection = parts.first;
+  Vector a = parts.second;
+  if(!a.IsZero()){
+    // construct a vector that is orthogonal to a and direction and has length |a|
+    Vector y = a;
+    // direction is normalized, so the result has length |a|
+    y.VectorProduct(*direction);
+
+    lineDirection += cos(alpha) * a + sin(alpha) * y;
+  }
+  return Line(lineOrigin,lineDirection);
+}
+
+Plane Line::rotatePlane(const Plane &rhs, double alpha) const{
+  vector<Vector> points = rhs.getPointsOnPlane();
+  transform(points.begin(),
+            points.end(),
+            points.begin(),
+            boost::bind(&Line::rotateVector,this,_1,alpha));
+  return Plane(points[0],points[1],points[2]);
+}
+
 Plane Line::getOrthogonalPlane(const Vector &origin) const{
   return Plane(getDirection(),origin);
@@ -247 +282 @@
 }
 
+LinePoint Line::getLinePoint(const Vector &point) const{
+  ASSERT(isContained(point),"Line point queried for point not on line");
+  Vector helper = point - (*origin);
+  double param = helper.ScalarProduct(*direction);
+  return LinePoint(*this,param);
+}
+
+LinePoint Line::posEndpoint() const{
+  return LinePoint(*this, numeric_limits<double>::infinity());
+}
+LinePoint Line::negEndpoint() const{
+  return LinePoint(*this,-numeric_limits<double>::infinity());
+}
+
+bool operator==(const Line &x,const Line &y){
+  return *x.origin == *y.origin && *x.direction == *y.direction;
+}
+
 Line makeLineThrough(const Vector &x1, const Vector &x2){
   if(x1==x2){
@@ -252 +305 @@
   }
   return Line(x1,x1-x2);
+}
+
+/******************************** Points on the line ********************/
+
+LinePoint::LinePoint(const LinePoint &src) :
+  line(src.line),param(src.param)
+{}
+
+LinePoint::LinePoint(const Line &_line, double _param) :
+  line(_line),param(_param)
+{}
+
+LinePoint& LinePoint::operator=(const LinePoint &src){
+  line=src.line;
+  param=src.param;
+  return *this;
+}
+
+Vector LinePoint::getPoint() const{
+  ASSERT(!isInfinite(),"getPoint() on infinite LinePoint called");
+  return (*line.origin)+param*(*line.direction);
+}
+
+Line LinePoint::getLine() const{
+  return line;
+}
+
+bool LinePoint::isInfinite() const{
+  return isPosInfinity() || isNegInfinity();
+}
+bool LinePoint::isPosInfinity() const{
+  return param == numeric_limits<double>::infinity();
+}
+bool LinePoint::isNegInfinity() const{
+  return param ==-numeric_limits<double>::infinity();
+}
+
+bool operator==(const LinePoint &x, const LinePoint &y){
+  ASSERT(x.line==y.line,"Operation on two points of different lines");
+  return x.param == y.param;
+
+}
+bool operator<(const LinePoint &x, const LinePoint &y){
+  ASSERT(x.line==y.line,"Operation on two points of different lines");
+  return x.param<y.param;
 }
 

src/LinearAlgebra/Line.hpp

@@ -17 +17 @@
 class Vector;
 class Plane;
+class LinePoint;
 
 class Line : public Space
 {
+  friend bool operator==(const Line&,const Line&);
+  friend class LinePoint;
 public:
   Line(const Vector &_origin, const Vector &_direction);
   Line(const Line& _src);
   virtual ~Line();
+
+  Line &operator=(const Line& rhs);
 
   virtual double distance(const Vector &point) const;
@@ -36 +41 @@
 
   Vector rotateVector(const Vector &rhs, double alpha) const;
+  Line rotateLine(const Line &rhs, double alpha) const;
+  Plane rotatePlane(const Plane &rhs, double alpha) const;
 
   Plane getOrthogonalPlane(const Vector &origin) const;
@@ -41 +48 @@
   std::vector<Vector> getSphereIntersections() const;
 
+  LinePoint getLinePoint(const Vector&) const;
+  LinePoint posEndpoint() const;
+  LinePoint negEndpoint() const;
+
+
+
 private:
   std::auto_ptr<Vector> origin;
   std::auto_ptr<Vector> direction;
 };
+
+bool operator==(const Line&,const Line&);
 
 std::ostream & operator << (std::ostream& ost, const Line &m);
@@ -53 +68 @@
 Line makeLineThrough(const Vector &x1, const Vector &x2);
 
+/**
+ * Class for representing points on a line
+ * These objects allow comparison of points on the same line as well as specifying the
+ * infinite "endpoints" of a line.
+ */
+class LinePoint{
+  friend class Line;
+  friend bool operator==(const LinePoint&, const LinePoint&);
+  friend bool operator<(const LinePoint&, const LinePoint&);
+public:
+  LinePoint(const LinePoint&);
+  LinePoint& operator=(const LinePoint&);
+  Vector getPoint() const;
+  Line getLine() const;
+  bool isInfinite() const;
+  bool isPosInfinity() const;
+  bool isNegInfinity() const;
+
+private:
+  LinePoint(const Line&,double);
+  Line line;
+  double param;
+};
+
+bool operator==(const LinePoint&, const LinePoint&);
+bool operator<(const LinePoint&, const LinePoint&);
+
+inline bool operator!= (const LinePoint& x, const LinePoint& y) { return !(x==y); }
+inline bool operator>  (const LinePoint& x, const LinePoint& y) { return y<x; }
+inline bool operator<= (const LinePoint& x, const LinePoint& y) { return !(y<x); }
+inline bool operator>= (const LinePoint& x, const LinePoint& y) { return !(x<y); }
+
+
 #endif /* LINE_HPP_ */

src/LinearAlgebra/Makefile.am

@@ -14 +14 @@
   gslvector.cpp \
   Line.cpp \
+  LineSegment.cpp \
+  LineSegmentSet.cpp \
   linearsystemofequations.cpp \
   Matrix.cpp \
@@ -27 +29 @@
   gslvector.hpp \
   Line.hpp \
+  LineSegment.hpp \
+  LineSegmentSet.hpp \
   linearsystemofequations.hpp \
   Matrix.hpp \

src/LinearAlgebra/Matrix.cpp

@@ -230 +230 @@
        - at(2,1)*at(1,2)*at(0,0)
        - at(2,2)*at(1,0)*at(0,1);
-}
-
-Matrix Matrix::transpose() const
-{
-  Matrix copy(*this);
-  copy.transpose();
-  return copy;
-}
-
-
-Matrix &Matrix::transpose()
-{
-  double tmp;
-  for (int i=0;i<NDIM;i++)
-    for (int j=i+1;j<NDIM;j++) {
-      tmp = at(j,i);
-      at(j,i) = at(i,j);
-      at(i,j) = tmp;
-    }
-  return *this;
 }
 
@@ -273 +253 @@
 }
 
+Matrix Matrix::transpose() const{
+  MatrixContent *newContent = new MatrixContent();
+  newContent->content = gsl_matrix_alloc(NDIM, NDIM);
+  gsl_matrix_transpose_memcpy(newContent->content, content->content);
+  Matrix res = Matrix(newContent);
+  return res;
+}
+
+Matrix &Matrix::transpose()
+{
+  double tmp;
+  for (int i=0;i<NDIM;i++)
+    for (int j=i+1;j<NDIM;j++) {
+      tmp = at(j,i);
+      at(j,i) = at(i,j);
+      at(i,j) = tmp;
+    }
+  return *this;
+}
+
+
 Vector Matrix::transformToEigenbasis()
 {

src/LinearAlgebra/Plane.cpp

@@ -27 +27 @@
 #include "Helpers/Verbose.hpp"
 #include "Helpers/Assert.hpp"
+#include "Helpers/helpers.hpp"
 #include <cmath>
 #include "LinearAlgebra/Line.hpp"
@@ -118 +119 @@
 Plane::~Plane()
 {}
+
+Plane &Plane::operator=(const Plane &rhs){
+  if(&rhs!=this){
+    normalVector.reset(new Vector(*rhs.normalVector));
+    offset = rhs.offset;
+  }
+  return *this;
+}
 
 
@@ -200 +209 @@
 }
 
+bool Plane::onSameSide(const Vector &point1,const Vector &point2) const{
+  return sign(point1.ScalarProduct(*normalVector)-offset) ==
+         sign(point2.ScalarProduct(*normalVector)-offset);
+}
+
 /************ Methods inherited from Space ****************/
 
@@ -219 +233 @@
 // Operators
 
+bool operator==(const Plane &x,const Plane &y){
+  return *x.normalVector == *y.normalVector && x.offset == y.offset;
+}
+
 ostream &operator << (ostream &ost,const Plane &p){
   ost << "<" << p.getNormal() << ";x> - " << p.getOffset() << "=0";

src/LinearAlgebra/Plane.hpp

@@ -21 +21 @@
 class Plane : public Space
 {
+  friend bool operator==(const Plane&,const Plane&);
   typedef std::auto_ptr<Vector> vec_ptr;
 public:
@@ -29 +30 @@
   Plane(const Plane& plane);
   virtual ~Plane();
+
+  Plane &operator=(const Plane&);
 
   // Accessor Functions
@@ -63 +66 @@
   Line getOrthogonalLine(const Vector &origin) const;
 
+  /**
+   * Test if two points are on the same side of the plane
+   */
+  bool onSameSide(const Vector&,const Vector&) const;
+
   /****** Methods inherited from Space ***********/
 
@@ -73 +81 @@
 };
 
+bool operator==(const Plane&,const Plane&);
+
 std::ostream &operator<< (std::ostream &ost,const Plane& p);
 

src/LinearAlgebra/Vector.cpp

@@ -237 +237 @@
   (*this) *= 1/factor;
 };
+
+Vector Vector::getNormalized() const{
+  Vector res= *this;
+  res.Normalize();
+  return res;
+}
 
 /** Zeros all components of this vector.
@@ -578 +584 @@
 // some comonly used vectors
 const Vector zeroVec(0,0,0);
-const Vector e1(1,0,0);
-const Vector e2(0,1,0);
-const Vector e3(0,0,1);
+const Vector unitVec[NDIM]={Vector(1,0,0),Vector(0,1,0),Vector(0,0,1)};

src/LinearAlgebra/Vector.hpp

@@ -83 +83 @@
   double NormSquared() const;
   void Normalize();
+  Vector getNormalized() const;
   void Zero();
   void One(const double one);
@@ -109 +110 @@
 // some commonly used and fixed vectors
 const extern Vector zeroVec;
-const extern Vector e1;
-const extern Vector e2;
-const extern Vector e3;
+const extern Vector unitVec[NDIM];
 
 ostream & operator << (ostream& ost, const Vector &m);

src/Makefile.am

@@ -58 +58 @@
   Actions/MethodAction.hpp \
   Actions/Process.hpp
+
 
 PATTERNSOURCE = \
@@ -227 +228 @@
   lists.hpp \
   molecule.hpp \
-  molecule_template.hpp \
   parser.hpp \
   periodentafel.hpp \

src/Shapes/BaseShapes.cpp

@@ -22 +22 @@
 #include "Shapes/BaseShapes.hpp"
 #include "Shapes/BaseShapes_impl.hpp"
-
-#include <cmath>
+#include "Shapes/ShapeOps.hpp"
+
+#include "defs.hpp"
 
 #include "Helpers/Assert.hpp"
 #include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/Plane.hpp"
+#include "LinearAlgebra/LineSegment.hpp"
+#include "LinearAlgebra/LineSegmentSet.hpp"
+
+#include <cmath>
+#include <algorithm>
 
 bool Sphere_impl::isInside(const Vector &point){
@@ -32 +40 @@
 }
 
+bool Sphere_impl::isOnSurface(const Vector &point){
+  return fabs(point.NormSquared()-1)<MYEPSILON;
+}
+
+Vector Sphere_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  return point;
+}
+
+LineSegmentSet Sphere_impl::getLineIntersections(const Line &line){
+  LineSegmentSet res(line);
+  std::vector<Vector> intersections = line.getSphereIntersections();
+  if(intersections.size()==2){
+    res.insert(LineSegment(intersections[0],intersections[1]));
+  }
+  return res;
+}
+
+string Sphere_impl::toString(){
+  return "Sphere()";
+}
 
 /**
@@ -66 +97 @@
 }
 
+Shape Sphere(const Vector &center,double radius){
+  return translate(resize(Sphere(),radius),center);
+}
+
+Shape Ellipsoid(const Vector &center, const Vector &radius){
+  return translate(stretch(Sphere(),radius),center);
+}
+
 bool Cuboid_impl::isInside(const Vector &point){
   return (point[0]>=0 && point[0]<=1) && (point[1]>=0 && point[1]<=1) && (point[2]>=0 && point[2]<=1);
+}
+
+bool Cuboid_impl::isOnSurface(const Vector &point){
+  bool retVal = isInside(point);
+  // test all borders of the cuboid
+  // double fabs
+  retVal = retVal &&
+           (((fabs(point[0]-1.)  < MYEPSILON) || (fabs(point[0])  < MYEPSILON)) ||
+            ((fabs(point[1]-1.)  < MYEPSILON) || (fabs(point[1])  < MYEPSILON)) ||
+            ((fabs(point[2]-1.)  < MYEPSILON) || (fabs(point[2])  < MYEPSILON)));
+  return retVal;
+}
+
+Vector Cuboid_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res;
+  // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
+  for(int i=NDIM;i--;){
+    if(fabs(fabs(point[i])-1)<MYEPSILON){
+      // add the scaled (-1/+1) Vector to the set of surface vectors
+      res[i] = point[i];
+    }
+  }
+  ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");
+
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line){
+  LineSegmentSet res(line);
+  // get the intersection on each of the six faces
+  vector<Vector> intersections;
+  intersections.resize(2);
+  int c=0;
+  int x[2]={-1,+1};
+  for(int i=NDIM;i--;){
+    for(int p=0;p<2;++p){
+      if(c==2) goto end; // I know this sucks, but breaking two loops is stupid
+      Vector base;
+      base[i]=x[p];
+      // base now points to the surface and is normal to it at the same time
+      Plane p(base,base);
+      Vector intersection = p.GetIntersection(line);
+      if(isInside(intersection)){
+        // if we have a point on the edge it might already be contained
+        if(c==1 && intersections[0]==intersection)
+          continue;
+        intersections[c++]=intersection;
+      }
+    }
+  }
+  end:
+  if(c==2){
+    res.insert(LineSegment(intersections[0],intersections[1]));
+  }
+  return res;
+}
+
+string Cuboid_impl::toString(){
+  return "Cuboid()";
 }
 
@@ -83 +185 @@
   return Shape(impl);
 }
+
+Shape Cuboid(const Vector &corner1, const Vector &corner2){
+  // make sure the two edges are upper left front and lower right back
+  Vector sortedC1;
+  Vector sortedC2;
+  for(int i=NDIM;i--;){
+    sortedC1[i] = min(corner1[i],corner2[i]);
+    sortedC2[i] = max(corner1[i],corner2[i]);
+    ASSERT(corner1[i]!=corner2[i],"Given points for cuboid edges did not define a valid space");
+  }
+  // get the middle point
+  Vector middle = (1./2.)*(sortedC1+sortedC2);
+  Vector factors = sortedC2-middle;
+  return translate(stretch(Cuboid(),factors),middle);
+}

src/Shapes/BaseShapes.hpp

@@ -12 +12 @@
 
 Shape Sphere();
+Shape Sphere(const Vector &center,double radius);
+Shape Ellipsoid(const Vector &center, const Vector &radius);
 Shape Cuboid();
+Shape Cuboid(const Vector &corner1, const Vector &corner2);
 
 #endif /* BASESHAPES_HPP_ */

src/Shapes/BaseShapes_impl.hpp

@@ -9 +9 @@
 #define BASESHAPES_IMPL_HPP_
 
+class Vector;
+class LineSegmentSet;
+
 #include "Shapes/Shape_impl.hpp"
 
 class Sphere_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 };
@@ -18 +25 @@
 class Cuboid_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 };

src/Shapes/Shape.cpp

@@ -20 +20 @@
 #include "Helpers/MemDebug.hpp"
 
+#include "Helpers/Assert.hpp"
+#include "LinearAlgebra/Vector.hpp"
+
 #include "Shape.hpp"
 #include "Shape_impl.hpp"
 
-
-#include "Helpers/Assert.hpp"
-#include "LinearAlgebra/Vector.hpp"
+#include <string>
+
 
 Shape::Shape(const Shape& src) :
@@ -35 +37 @@
 bool Shape::isInside(const Vector &point) const{
   return impl->isInside(point);
+}
+
+bool Shape::isOnSurface(const Vector &point) const{
+  return impl->isOnSurface(point);
+}
+
+Vector Shape::getNormal(const Vector &point) const throw (NotOnSurfaceException){
+  return impl->getNormal(point);
+}
+
+LineSegmentSet Shape::getLineIntersections(const Line &line){
+  return impl->getLineIntersections(line);
 }
 
@@ -52 +66 @@
 }
 
+std::string Shape::toString() const{
+  return impl->toString();
+}
+
 Shape::impl_ptr Shape::getImpl() const{
   return impl;
@@ -85 +103 @@
 bool AndShape_impl::isInside(const Vector &point){
   return lhs->isInside(point) && rhs->isInside(point);
+}
+
+bool AndShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector AndShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet AndShape_impl::getLineIntersections(const Line &line){
+  return intersect(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
+}
+
+string AndShape_impl::toString(){
+  return string("(") + lhs->toString() + string("&&") + rhs->toString() + string(")");
 }
 
@@ -122 +191 @@
 }
 
+bool OrShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && !rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && !lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector OrShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet OrShape_impl::getLineIntersections(const Line &line){
+  return merge(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
+}
+
+string OrShape_impl::toString(){
+  return string("(") + lhs->toString() + string("||") + rhs->toString() + string(")");
+}
+
 std::vector<Vector> OrShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
   std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
@@ -156 +276 @@
 }
 
+bool NotShape_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(point);
+}
+
+Vector NotShape_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  return -1*arg->getNormal(point);
+}
+
+LineSegmentSet NotShape_impl::getLineIntersections(const Line &line){
+  return invert(arg->getLineIntersections(line));
+}
+
+string NotShape_impl::toString(){
+  return string("!") + arg->toString();
+}
+
 std::vector<Vector> NotShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
   // surfaces are the same, only normal direction is different
@@ -165 +301 @@
   return Shape(newImpl);
 }
+
+/**************** global operations *********************************/
+ostream &operator<<(ostream &ost,const Shape &shape){
+  ost << shape.toString();
+  return ost;
+}

src/Shapes/Shape.hpp

@@ -10 +10 @@
 
 #include <boost/shared_ptr.hpp>
+#include <iosfwd>
+
+#include "Exceptions/NotOnSurfaceException.hpp"
 
 #include <vector>
@@ -15 +18 @@
 class Vector;
 class Shape_impl;
+class LineSegmentSet;
+class Line;
 
 class Shape
@@ -27 +32 @@
 
   bool isInside(const Vector &point) const;
+  bool isOnSurface(const Vector &point) const;
+  Vector getNormal(const Vector &point) const throw(NotOnSurfaceException);
+
+  LineSegmentSet getLineIntersections(const Line&);
   std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 
   Shape &operator=(const Shape& rhs);
 
+  std::string toString() const;
 protected:
   impl_ptr getImpl() const;
@@ -45 +55 @@
 Shape operator!(const Shape&);
 
+std::ostream &operator<<(std::ostream&,const Shape&);
+
 #endif /* SHAPE_HPP_ */

src/Shapes/ShapeOps.cpp

@@ -23 +23 @@
 #include "Shapes/ShapeOps_impl.hpp"
 
+#include "LinearAlgebra/Vector.hpp"
 #include "Helpers/Assert.hpp"
 
+/*************** Base case ***********************/
+
+ShapeOpsBase_impl::ShapeOpsBase_impl(const Shape::impl_ptr &_arg) :
+  arg(_arg){}
+
+ShapeOpsBase_impl::~ShapeOpsBase_impl(){}
+
+bool ShapeOpsBase_impl::isInside(const Vector &point){
+  return arg->isInside(translateIn(point));
+}
+
+bool ShapeOpsBase_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(translateIn(point));
+}
+
+Vector ShapeOpsBase_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector helper = translateIn(point);
+  if(!arg->isOnSurface(helper)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res = translateOutNormal(arg->getNormal(helper));
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet ShapeOpsBase_impl::getLineIntersections(const Line &line){
+  Line newLine(translateIn(line.getOrigin()),translateIn(line.getDirection()));
+  LineSegmentSet res(line);
+  LineSegmentSet helper = getArg()->getLineIntersections(newLine);
+  for(LineSegmentSet::iterator iter = helper.begin();iter!=helper.end();++iter){
+    LinePoint lpBegin = iter->getBegin();
+    LinePoint lpEnd = iter->getBegin();
+    // translate both linepoints
+    lpBegin = lpBegin.isNegInfinity()?
+                line.negEndpoint():
+                line.getLinePoint(translateOutPos(lpBegin.getPoint()));
+    lpEnd = lpEnd.isPosInfinity()?
+              line.posEndpoint():
+              line.getLinePoint(translateOutPos(lpEnd.getPoint()));
+    res.insert(LineSegment(lpBegin,lpEnd));
+  }
+  return res;
+}
+
+std::vector<Vector> ShapeOpsBase_impl::getHomogeneousPointsOnSurface(const size_t N) const {
+  return getArg()->getHomogeneousPointsOnSurface(N);;
+}
+
+Shape::impl_ptr ShapeOpsBase_impl::getArg() const{
+  return arg;
+}
+
 /********************* Resize ********************/
 
 Resize_impl::Resize_impl(const Shape::impl_ptr &_arg,double _size) :
-  arg(_arg), size(_size)
+  ShapeOpsBase_impl(_arg), size(_size)
 {
   ASSERT(size>0,"Cannot resize a Shape to size zero or below");
@@ -36 +89 @@
 
 bool Resize_impl::isInside(const Vector& point){
-  return arg->isInside((1/size) * point);
+  return getArg()->isInside((1/size) * point);
+}
+
+Vector Resize_impl::translateIn(const Vector& point){
+  return (1/size) * point;
+}
+
+Vector Resize_impl::translateOutPos(const Vector& point){
+  return size * point;
+}
+
+Vector Resize_impl::translateOutNormal(const Vector& point){
+  return point;
+}
+
+string Resize_impl::toString(){
+  stringstream sstr;
+  sstr << "resize(" << getArg()->toString() << "," << size << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Resize_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     *iter *= size;
@@ -56 +127 @@
 
 Translate_impl::Translate_impl(const Shape::impl_ptr &_arg, const Vector &_offset) :
-  arg(_arg),offset(_offset)
+  ShapeOpsBase_impl(_arg),offset(_offset)
 {}
 
@@ -62 +133 @@
 
 bool Translate_impl::isInside(const Vector& point){
-  return arg->isInside(point-offset);
+  return getArg()->isInside(point-offset);
+}
+
+Vector Translate_impl::translateIn(const Vector& point){
+  return point-offset;
+}
+
+Vector Translate_impl::translateOutPos(const Vector& point){
+  return point+offset;
+}
+
+Vector Translate_impl::translateOutNormal(const Vector& point){
+  return point;
+}
+
+string Translate_impl::toString(){
+  stringstream sstr;
+  sstr << "translate(" << getArg()->toString() << "," << offset << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Translate_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     *iter += offset;
@@ -81 +170 @@
 
 Stretch_impl::Stretch_impl(const Shape::impl_ptr &_arg, const Vector &_factors) :
-  arg(_arg),factors(_factors)
+  ShapeOpsBase_impl(_arg),factors(_factors)
 {
   ASSERT(factors[0]>0,"cannot stretch a shape by a negative amount");
@@ -96 +185 @@
   Vector helper=point;
   helper.ScaleAll(reciFactors);
-  return arg->isInside(helper);
+  return getArg()->isInside(helper);
+}
+
+Vector Stretch_impl::translateIn(const Vector& point){
+  Vector helper=point;
+  helper.ScaleAll(reciFactors);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutPos(const Vector& point){
+  Vector helper=point;
+  helper.ScaleAll(factors);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutNormal(const Vector& point){
+  Vector helper=point;
+  // the normalFactors are derived from appearances of the factors
+  // with in the vectorproduct
+  Vector normalFactors;
+  normalFactors[0]=factors[1]*factors[2];
+  normalFactors[1]=factors[0]*factors[2];
+  normalFactors[2]=factors[0]*factors[1];
+  helper.ScaleAll(normalFactors);
+  return helper;
+}
+
+string Stretch_impl::toString(){
+  stringstream sstr;
+  sstr << "stretch(" << getArg()->toString() << "," << factors << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Stretch_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     (*iter).ScaleAll(reciFactors);
@@ -115 +234 @@
 
 Transform_impl::Transform_impl(const Shape::impl_ptr &_arg, const Matrix &_transformation) :
-  arg(_arg),transformation(_transformation)
+  ShapeOpsBase_impl(_arg),transformation(_transformation)
 {
   transformationInv = transformation.invert();
@@ -123 +242 @@
 
 bool Transform_impl::isInside(const Vector& point){
-  return arg->isInside(transformationInv * point);
+  return getArg()->isInside(transformationInv * point);
+}
+
+Vector Transform_impl::translateIn(const Vector& point){
+  return transformationInv * point;
+}
+
+Vector Transform_impl::translateOutPos(const Vector& point){
+  return transformation * point;
+}
+
+Vector Transform_impl::translateOutNormal(const Vector& point){
+  Matrix mat = transformation.invert().transpose();
+  return mat * point;
+}
+
+string Transform_impl::toString(){
+  stringstream sstr;
+  sstr << "transform(" << getArg()->toString() << "," << transformation << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Transform_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     *iter = transformation * (*iter);

src/Shapes/ShapeOps_impl.hpp

@@ -13 +13 @@
 #include "LinearAlgebra/Matrix.hpp"
 
-class Resize_impl :  public Shape_impl
+#include <vector>
+
+class LineSegment;
+
+class ShapeOpsBase_impl : public Shape_impl{
+public:
+  ShapeOpsBase_impl(const Shape::impl_ptr&);
+  virtual ~ShapeOpsBase_impl();
+  virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw (NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
+protected:
+  virtual Vector translateIn(const Vector &point)=0;
+  virtual Vector translateOutPos(const Vector &point)=0;
+  virtual Vector translateOutNormal(const Vector &point)=0;
+  Shape::impl_ptr getArg() const;
+private:
+  Shape::impl_ptr arg;
+};
+
+class Resize_impl :  public ShapeOpsBase_impl
 {
 public:
   Resize_impl(const Shape::impl_ptr&,double);
   virtual ~Resize_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   double size;
 };
 
-class Translate_impl :  public Shape_impl
+class Translate_impl :  public ShapeOpsBase_impl
 {
 public:
   Translate_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Translate_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   Vector offset;
 };
 
-class Stretch_impl : public Shape_impl
+class Stretch_impl : public ShapeOpsBase_impl
 {
 public:
   Stretch_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Stretch_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   Vector factors;
   Vector reciFactors;
 };
 
-class Transform_impl : public Shape_impl
+class Transform_impl : public ShapeOpsBase_impl
 {
 public:
   Transform_impl(const Shape::impl_ptr&, const Matrix&);
   virtual ~Transform_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   Matrix transformation;
   Matrix transformationInv;

src/Shapes/Shape_impl.hpp

@@ -12 +12 @@
 
 #include "Shapes/Shape.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/LineSegment.hpp"
+#include "LinearAlgebra/LineSegmentSet.hpp"
+#include "LinearAlgebra/Vector.hpp"
 
-class Vector;
 
 class Shape_impl {
@@ -20 +23 @@
   virtual ~Shape_impl(){};
   virtual bool isInside(const Vector &point)=0;
+  virtual bool isOnSurface(const Vector &point)=0;
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException)=0;
+  virtual LineSegmentSet getLineIntersections(const Line&)=0;
+  virtual std::string toString()=0;
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const=0;
 };
@@ -27 +34 @@
   virtual bool isInside(const Vector &point){
     return true;
+  }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  virtual LineSegmentSet getLineIntersections(const Line &line){
+    LineSegmentSet res(line);
+    res.insert(LineSegment(line.negEndpoint(),line.posEndpoint()));
+    return res;
+  }
+  virtual std::string toString(){
+    return "Everywhere()";
   }
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const {
@@ -37 +58 @@
   virtual bool isInside(const Vector &point){
     return false;
+  }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  virtual LineSegmentSet getLineIntersections(const Line &line){
+    return LineSegmentSet(line);
+  }
+  virtual std::string toString(){
+    return "Nowhere()";
   }
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const {
@@ -49 +82 @@
   virtual ~AndShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
@@ -60 +97 @@
   virtual ~OrShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
@@ -71 +112 @@
   virtual ~NotShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:

src/unittests/LineUnittest.cpp

@@ -41 +41 @@
 void LineUnittest::setUp(){
   // three lines along the axes
-  la1 = new Line(zeroVec,e1);
-  la2 = new Line(zeroVec,e2);
-  la3 = new Line(zeroVec,e3);
+  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(e1,e1-e2);
-  lp2 = new Line(e2,e2-e3);
-  lp3 = new Line(e3,e3-e1);
+  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(){
@@ -64 +64 @@
 
   // direction+origin should never fail
-  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,e1));
-  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,e2));
-  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,e3));
+  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(e1,e2));
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e2,e3));
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e3,e1));
+  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(e1,zeroVec));
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e2,zeroVec));
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e3,zeroVec));
+  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(e1,e1),LinearDependenceException);
-  CPPUNIT_ASSERT_THROW(makeLineThrough(e2,e2),LinearDependenceException);
-  CPPUNIT_ASSERT_THROW(makeLineThrough(e3,e3),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);
 
 }
@@ -90 +90 @@
 void LineUnittest::constructionResultTest(){
   // test all directions
-  CPPUNIT_ASSERT(testDirection(la1->getDirection(),e1));
-  CPPUNIT_ASSERT(testDirection(la2->getDirection(),e2));
-  CPPUNIT_ASSERT(testDirection(la3->getDirection(),e3));
+  CPPUNIT_ASSERT(testDirection(la1->getDirection(),unitVec[0]));
+  CPPUNIT_ASSERT(testDirection(la2->getDirection(),unitVec[1]));
+  CPPUNIT_ASSERT(testDirection(la3->getDirection(),unitVec[2]));
 
   // test origins
@@ -104 +104 @@
   CPPUNIT_ASSERT(la3->isContained(zeroVec));
 
-  CPPUNIT_ASSERT(la1->isContained(e1));
-  CPPUNIT_ASSERT(la2->isContained(e2));
-  CPPUNIT_ASSERT(la3->isContained(e3));
-
-  CPPUNIT_ASSERT(lp1->isContained(e1));
-  CPPUNIT_ASSERT(lp2->isContained(e2));
-  CPPUNIT_ASSERT(lp3->isContained(e3));
-
-  CPPUNIT_ASSERT(lp1->isContained(e2));
-  CPPUNIT_ASSERT(lp2->isContained(e3));
-  CPPUNIT_ASSERT(lp3->isContained(e1));
+  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]));
 }
 
@@ -124 +124 @@
 
   // multiples of the second support vector
-  CPPUNIT_ASSERT(la1->isContained(e1));
-  CPPUNIT_ASSERT(la2->isContained(e2));
-  CPPUNIT_ASSERT(la3->isContained(e3));
-
-  CPPUNIT_ASSERT(la1->isContained(2*e1));
-  CPPUNIT_ASSERT(la2->isContained(2*e2));
-  CPPUNIT_ASSERT(la3->isContained(2*e3));
-
-  CPPUNIT_ASSERT(la1->isContained(3*e1));
-  CPPUNIT_ASSERT(la2->isContained(3*e2));
-  CPPUNIT_ASSERT(la3->isContained(3*e3));
+  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*e1));
-  CPPUNIT_ASSERT(la2->isContained(-1*e2));
-  CPPUNIT_ASSERT(la3->isContained(-1*e3));
-
-  CPPUNIT_ASSERT(la1->isContained(-2*e1));
-  CPPUNIT_ASSERT(la2->isContained(-2*e2));
-  CPPUNIT_ASSERT(la3->isContained(-2*e3));
+  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(e2));
-  CPPUNIT_ASSERT(!la2->isContained(e3));
-  CPPUNIT_ASSERT(!la3->isContained(e1));
-
-  CPPUNIT_ASSERT(!la1->isContained(2*e2));
-  CPPUNIT_ASSERT(!la2->isContained(2*e3));
-  CPPUNIT_ASSERT(!la3->isContained(2*e1));
-
-  CPPUNIT_ASSERT(!la1->isContained(-1*e2));
-  CPPUNIT_ASSERT(!la2->isContained(-1*e3));
-  CPPUNIT_ASSERT(!la3->isContained(-1*e1));
+  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(e1));
-  CPPUNIT_ASSERT(lp2->isContained(e2));
-  CPPUNIT_ASSERT(lp3->isContained(e3));
-
-  CPPUNIT_ASSERT(lp1->isContained(e2));
-  CPPUNIT_ASSERT(lp2->isContained(e3));
-  CPPUNIT_ASSERT(lp3->isContained(e1));
-
-  CPPUNIT_ASSERT(lp1->isContained(e1+2*(e1-e2)));
-  CPPUNIT_ASSERT(lp2->isContained(e2+2*(e2-e3)));
-  CPPUNIT_ASSERT(lp3->isContained(e3+2*(e3-e1)));
-
-  CPPUNIT_ASSERT(lp1->isContained(e1-2*(e1-e2)));
-  CPPUNIT_ASSERT(lp2->isContained(e2-2*(e2-e3)));
-  CPPUNIT_ASSERT(lp3->isContained(e3-2*(e3-e1)));
+  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])));
 }
 
@@ -189 +189 @@
   // axes and plane lines
   fixture = la1->getIntersection(*lp1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
   fixture = la2->getIntersection(*lp2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
   fixture = la3->getIntersection(*lp3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   fixture = la1->getIntersection(*lp3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
   fixture = la2->getIntersection(*lp1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
   fixture = la3->getIntersection(*lp2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   // two plane lines
   fixture = lp1->getIntersection(*lp2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
   fixture = lp2->getIntersection(*lp3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
   fixture = lp3->getIntersection(*lp1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
 
   // When we have two times the same line, we check if the point is on the line
@@ -254 +254 @@
 
   // rotate vectors on the axis around their lines
-  fixture = la1->rotateVector(e1,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = la2->rotateVector(e2,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = la3->rotateVector(e3,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-
-  fixture = la1->rotateVector(e1,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = la2->rotateVector(e2,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = la3->rotateVector(e3,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  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*e1,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,2*e1);
-  fixture = la2->rotateVector(2*e2,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,2*e2);
-  fixture = la3->rotateVector(2*e3,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,2*e3);
-
-  fixture = la1->rotateVector(2*e1,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,2*e1);
-  fixture = la2->rotateVector(2*e2,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,2*e2);
-  fixture = la3->rotateVector(2*e3,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,2*e3);
+  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*e1,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
-  fixture = la2->rotateVector(-1*e2,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
-  fixture = la3->rotateVector(-1*e3,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
-
-  fixture = la1->rotateVector(-1*e1,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
-  fixture = la2->rotateVector(-1*e2,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
-  fixture = la3->rotateVector(-1*e3,2.);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+  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
-  // e2 around e1
-  fixture = la1->rotateVector(e2,0);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = la1->rotateVector(e2,1./2.*M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
-  fixture = la1->rotateVector(e2,M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
-  fixture = la1->rotateVector(e2,2*M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-
-  // e3 around e2
-  fixture = la2->rotateVector(e3,0);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-  fixture = la2->rotateVector(e3,1./2.*M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
-  fixture = la2->rotateVector(e3,M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
-  fixture = la2->rotateVector(e3,2*M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-
-  // e1 around e3
-  fixture = la3->rotateVector(e1,0);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = la3->rotateVector(e1,1./2.*M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
-  fixture = la3->rotateVector(e1,M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
-  fixture = la3->rotateVector(e1,2*M_PI);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  // 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]);
 
 
@@ -335 +335 @@
 
   // Vectors on the line
-  fixture = lp1->rotateVector(e1,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = lp1->rotateVector(e2,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-
-  fixture = lp2->rotateVector(e2,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = lp2->rotateVector(e3,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-
-  fixture = lp3->rotateVector(e3,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-  fixture = lp3->rotateVector(e1,1.);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  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
@@ -370 +370 @@
     std::vector<Vector> res = la1->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT(testDirection(res[0],e1));
-    CPPUNIT_ASSERT(testDirection(res[1],e1));
+    CPPUNIT_ASSERT(testDirection(res[0],unitVec[0]));
+    CPPUNIT_ASSERT(testDirection(res[1],unitVec[0]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -378 +378 @@
     std::vector<Vector> res = la2->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT(testDirection(res[0],e2));
-    CPPUNIT_ASSERT(testDirection(res[1],e2));
+    CPPUNIT_ASSERT(testDirection(res[0],unitVec[1]));
+    CPPUNIT_ASSERT(testDirection(res[1],unitVec[1]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -386 +386 @@
     std::vector<Vector> res = la3->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT(testDirection(res[0],e3));
-    CPPUNIT_ASSERT(testDirection(res[1],e3));
+    CPPUNIT_ASSERT(testDirection(res[0],unitVec[2]));
+    CPPUNIT_ASSERT(testDirection(res[1],unitVec[2]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -394 +394 @@
     std::vector<Vector> res = lp1->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT((res[0]==e1) || (res[0]==e2));
-    CPPUNIT_ASSERT((res[1]==e1) || (res[1]==e2));
+    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]);
   }
@@ -402 +402 @@
     std::vector<Vector> res = lp2->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT((res[0]==e2) || (res[0]==e3));
-    CPPUNIT_ASSERT((res[1]==e2) || (res[1]==e3));
+    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]);
   }
@@ -410 +410 @@
     std::vector<Vector> res = lp3->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT((res[0]==e3) || (res[0]==e1));
-    CPPUNIT_ASSERT((res[1]==e3) || (res[1]==e1));
-    CPPUNIT_ASSERT(res[0]!=res[1]);
-  }
-}
+    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]);
+  }
+}

src/unittests/Makefile.am

@@ -16 +16 @@
   atomsCalculationTest \
   AtomDescriptorTest \
-  BaseShapesUnitTest \
   BondGraphUnitTest \
   BoxUnittest \
@@ -75 +74 @@
   AtomDescriptorTest.cpp \
   atomsCalculationTest.cpp \
-  BaseShapesUnittest.cpp \
   bondgraphunittest.cpp \
   BoxUnittest.cpp \
@@ -114 +112 @@
   AtomDescriptorTest.hpp \
   atomsCalculationTest.hpp \
-  BaseShapesUnittest.hpp \
   bondgraphunittest.hpp \
   BoxUnittest.hpp \
@@ -168 +165 @@
 AtomDescriptorTest_LDADD = ${ALLLIBS}
 
-BaseShapesUnitTest_SOURCES = UnitTestMain.cpp BaseShapesUnittest.cpp BaseShapesUnittest.hpp
-BaseShapesUnitTest_LDADD = ${ALLLIBS}
-
 BondGraphUnitTest_SOURCES = UnitTestMain.cpp bondgraphunittest.cpp bondgraphunittest.hpp
 BondGraphUnitTest_LDADD = ${ALLLIBS}

src/unittests/MatrixUnittest.cpp

@@ -25 +25 @@
 
 #include "MatrixUnittest.hpp"
+#include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/Matrix.hpp"
-#include "LinearAlgebra/Vector.hpp"
 #include "Exceptions/NotInvertibleException.hpp"
 
@@ -53 +53 @@
   }
   perm1 = new Matrix();
-  perm1->column(0) = e1;
-  perm1->column(1) = e3;
-  perm1->column(2) = e2;
+  perm1->column(0) = unitVec[0];
+  perm1->column(1) = unitVec[2];
+  perm1->column(2) = unitVec[1];
 
 
   perm2 = new Matrix();
-  perm2->column(0) = e2;
-  perm2->column(1) = e1;
-  perm2->column(2) = e3;
+  perm2->column(0) = unitVec[1];
+  perm2->column(1) = unitVec[0];
+  perm2->column(2) = unitVec[2];
 
   perm3 = new Matrix();
-  perm3->column(0) = e2;
-  perm3->column(1) = e3;
-  perm3->column(2) = e1;
+  perm3->column(0) = unitVec[1];
+  perm3->column(1) = unitVec[2];
+  perm3->column(2) = unitVec[0];
 
   perm4 = new Matrix();
-  perm4->column(0) = e3;
-  perm4->column(1) = e2;
-  perm4->column(2) = e1;
+  perm4->column(0) = unitVec[2];
+  perm4->column(1) = unitVec[1];
+  perm4->column(2) = unitVec[0];
 
   perm5 = new Matrix();
-  perm5->column(0) = e3;
-  perm5->column(1) = e1;
-  perm5->column(2) = e2;
+  perm5->column(0) = unitVec[2];
+  perm5->column(1) = unitVec[0];
+  perm5->column(2) = unitVec[1];
 
 }
@@ -122 +122 @@
 
   mat.one();
-  CPPUNIT_ASSERT_EQUAL(mat.row(0),e1);
-  CPPUNIT_ASSERT_EQUAL(mat.row(1),e2);
-  CPPUNIT_ASSERT_EQUAL(mat.row(2),e3);
-  CPPUNIT_ASSERT_EQUAL(mat.column(0),e1);
-  CPPUNIT_ASSERT_EQUAL(mat.column(1),e2);
-  CPPUNIT_ASSERT_EQUAL(mat.column(2),e3);
+  CPPUNIT_ASSERT_EQUAL(mat.row(0),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(mat.row(1),unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(mat.row(2),unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(mat.column(0),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(mat.column(1),unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(mat.column(2),unitVec[2]);
 
   Vector t1=Vector(1.,1.,1.);
@@ -232 +232 @@
 
   res = (*diagonal)*(*perm1);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),e1);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),2*e2);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
   res = (*diagonal)*(*perm2);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),2*e2);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),e1);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),3*e3);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),3*unitVec[2]);
   res = (*diagonal)*(*perm3);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),2*e2);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),e1);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
   res = (*diagonal)*(*perm4);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),2*e2);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),e1);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
   res = (*diagonal)*(*perm5);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),e1);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),2*e2);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
 }
 
@@ -310 +310 @@
 
 void MatrixUnittest::VecMultTest(){
-  CPPUNIT_ASSERT_EQUAL((*zero)*e1,zeroVec);
-  CPPUNIT_ASSERT_EQUAL((*zero)*e2,zeroVec);
-  CPPUNIT_ASSERT_EQUAL((*zero)*e3,zeroVec);
+  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[0],zeroVec);
+  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[1],zeroVec);
+  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[2],zeroVec);
   CPPUNIT_ASSERT_EQUAL((*zero)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*one)*e1,e1);
-  CPPUNIT_ASSERT_EQUAL((*one)*e2,e2);
-  CPPUNIT_ASSERT_EQUAL((*one)*e3,e3);
+  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[0],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[1],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[2],unitVec[2]);
   CPPUNIT_ASSERT_EQUAL((*one)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*diagonal)*e1,e1);
-  CPPUNIT_ASSERT_EQUAL((*diagonal)*e2,2*e2);
-  CPPUNIT_ASSERT_EQUAL((*diagonal)*e3,3*e3);
+  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[0],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[1],2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[2],3*unitVec[2]);
   CPPUNIT_ASSERT_EQUAL((*diagonal)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm1)*e1,e1);
-  CPPUNIT_ASSERT_EQUAL((*perm1)*e2,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm1)*e3,e2);
+  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[0],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[1],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[2],unitVec[1]);
   CPPUNIT_ASSERT_EQUAL((*perm1)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm2)*e1,e2);
-  CPPUNIT_ASSERT_EQUAL((*perm2)*e2,e1);
-  CPPUNIT_ASSERT_EQUAL((*perm2)*e3,e3);
+  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[0],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[1],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[2],unitVec[2]);
   CPPUNIT_ASSERT_EQUAL((*perm2)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm3)*e1,e2);
-  CPPUNIT_ASSERT_EQUAL((*perm3)*e2,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm3)*e3,e1);
+  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[0],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[1],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[2],unitVec[0]);
   CPPUNIT_ASSERT_EQUAL((*perm3)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm4)*e1,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm4)*e2,e2);
-  CPPUNIT_ASSERT_EQUAL((*perm4)*e3,e1);
+  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[0],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[1],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[2],unitVec[0]);
   CPPUNIT_ASSERT_EQUAL((*perm4)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm5)*e1,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm5)*e2,e1);
-  CPPUNIT_ASSERT_EQUAL((*perm5)*e3,e2);
+  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[0],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[1],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[2],unitVec[1]);
   CPPUNIT_ASSERT_EQUAL((*perm5)*zeroVec,zeroVec);
 

src/unittests/PlaneUnittest.cpp

@@ -36 +36 @@
 
 void PlaneUnittest::setUp(){
-  p1 = new Plane(e1,e2,e3);
-  p2 = new Plane(e1,e2,zeroVec);
-  p3 = new Plane(e1,zeroVec,e3);
-  p4 = new Plane(zeroVec,e2,e3);
+  p1 = new Plane(unitVec[0],unitVec[1],unitVec[2]);
+  p2 = new Plane(unitVec[0],unitVec[1],zeroVec);
+  p3 = new Plane(unitVec[0],zeroVec,unitVec[2]);
+  p4 = new Plane(zeroVec,unitVec[1],unitVec[2]);
 }
 
@@ -54 +54 @@
 
   // three points
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,e2,e3));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],unitVec[1],unitVec[2]));
   // when only two points are differnt this gives an error
-  CPPUNIT_ASSERT_THROW(Plane(e1,e2,e2),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[1],unitVec[1]),LinearDependenceException);
   // same with only one point
-  CPPUNIT_ASSERT_THROW(Plane(e1,e1,e1),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[0],unitVec[0]),LinearDependenceException);
 
   // use two vector giving two directions
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,e2,0));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],unitVec[1],0));
   // and again this is actually only one vector
-  CPPUNIT_ASSERT_THROW(Plane(e1,e1,0),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[0],0),LinearDependenceException);
   // Zero vector does not give a good direction
-  CPPUNIT_ASSERT_THROW(Plane(e1,zeroVec,0),ZeroVectorException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],zeroVec,0),ZeroVectorException);
 
   // use a normalvector and an scalar offset
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,0));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],0));
   // The zero vector is no good as a normalvector
   CPPUNIT_ASSERT_THROW(Plane(zeroVec,0),ZeroVectorException);
 
   // use a normalvector and an offset vector
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,zeroVec));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],zeroVec));
   // and the bad zeroVector again
   CPPUNIT_ASSERT_THROW(Plane(zeroVec,zeroVec),ZeroVectorException);
@@ -87 +87 @@
   {
     // construct with three points on plane
-    Plane p1(e1,e2,zeroVec);
-    CPPUNIT_ASSERT(testNormal(e3,p1.getNormal()));
+    Plane p1(unitVec[0],unitVec[1],zeroVec);
+    CPPUNIT_ASSERT(testNormal(unitVec[2],p1.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());
 
-    Plane p2(e1,e3,zeroVec);
-    CPPUNIT_ASSERT(testNormal(e2,p2.getNormal()));
+    Plane p2(unitVec[0],unitVec[2],zeroVec);
+    CPPUNIT_ASSERT(testNormal(unitVec[1],p2.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());
 
-    Plane p3(e2,e3,zeroVec);
-    CPPUNIT_ASSERT(testNormal(e1,p3.getNormal()));
+    Plane p3(unitVec[1],unitVec[2],zeroVec);
+    CPPUNIT_ASSERT(testNormal(unitVec[0],p3.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
   }
   {
     // construct with two directions + offset
-    Plane p1(e1,e2,0);
-    CPPUNIT_ASSERT(testNormal(e3,p1.getNormal()));
+    Plane p1(unitVec[0],unitVec[1],0);
+    CPPUNIT_ASSERT(testNormal(unitVec[2],p1.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());
 
-    Plane p2(e1,e3,0);
-    CPPUNIT_ASSERT(testNormal(e2,p2.getNormal()));
+    Plane p2(unitVec[0],unitVec[2],0);
+    CPPUNIT_ASSERT(testNormal(unitVec[1],p2.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());
 
-    Plane p3(e2,e3,0);
-    CPPUNIT_ASSERT(testNormal(e1,p3.getNormal()));
+    Plane p3(unitVec[1],unitVec[2],0);
+    CPPUNIT_ASSERT(testNormal(unitVec[0],p3.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
   }
@@ -157 +157 @@
 void PlaneUnittest::operationsTest(){
   {
-    Vector t = (1./3.)*(e1+e2+e3);
+    Vector t = (1./3.)*(unitVec[0]+unitVec[1]+unitVec[2]);
     CPPUNIT_ASSERT(fabs(p1->distance(zeroVec)-t.Norm()) < MYEPSILON);
     CPPUNIT_ASSERT_EQUAL(t,p1->getClosestPoint(zeroVec));
   }
 
-  CPPUNIT_ASSERT(fabs(p2->distance(e3)-1) < MYEPSILON);
-  CPPUNIT_ASSERT_EQUAL(zeroVec,p2->getClosestPoint(e3));
-  CPPUNIT_ASSERT(fabs(p3->distance(e2)-1) < MYEPSILON);
-  CPPUNIT_ASSERT_EQUAL(zeroVec,p3->getClosestPoint(e2));
-  CPPUNIT_ASSERT(fabs(p4->distance(e1)-1) < MYEPSILON);
-  CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(e1));
+  CPPUNIT_ASSERT(fabs(p2->distance(unitVec[2])-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p2->getClosestPoint(unitVec[2]));
+  CPPUNIT_ASSERT(fabs(p3->distance(unitVec[1])-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p3->getClosestPoint(unitVec[1]));
+  CPPUNIT_ASSERT(fabs(p4->distance(unitVec[0])-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(unitVec[0]));
 }
 
@@ -174 +174 @@
 
   // some Vectors that lie on the planes
-  fixture = p1->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = p1->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = p1->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = p1->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = p1->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = p1->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   fixture = p2->mirrorVector(zeroVec);
   CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
-  fixture = p2->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = p2->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = p2->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = p2->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
 
   fixture = p3->mirrorVector(zeroVec);
   CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
-  fixture = p3->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = p3->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = p3->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = p3->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   fixture = p4->mirrorVector(zeroVec);
   CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
-  fixture = p4->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = p4->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = p4->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = p4->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   // some Vectors outside of the planes
   {
-    Vector t = (2./3.)*(e1+e2+e3);
+    Vector t = (2./3.)*(unitVec[0]+unitVec[1]+unitVec[2]);
     fixture = p1->mirrorVector(zeroVec);
     CPPUNIT_ASSERT_EQUAL(fixture,t);
   }
 
-  fixture = p2->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
-  fixture = p3->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
-  fixture = p4->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = p2->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
+  fixture = p3->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
+  fixture = p4->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
 }
 
@@ -221 +221 @@
   // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
   Line l1 = makeLineThrough(zeroVec,Vector(2,1,0));
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(e1, zeroVec).GetIntersection(l1) );
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unitVec[0], zeroVec).GetIntersection(l1) );
   CPPUNIT_ASSERT_EQUAL( zeroVec, fixture );
 
   // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
-  Line l2 = makeLineThrough(e1,Vector(0,1,1));
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(e2, Vector(2,1,0)).GetIntersection(l2) );
+  Line l2 = makeLineThrough(unitVec[0],Vector(0,1,1));
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unitVec[1], Vector(2,1,0)).GetIntersection(l2) );
   CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
 }

src/unittests/ShapeUnittest.cpp

@@ -40 +40 @@
 void ShapeUnittest::setUp()
 {
-  v000 =  0*e1+0*e2+0*e3;
-  v100 =  1*e1+0*e2+0*e3;
-  v200 = -1*e1+0*e2+0*e3;
-  v010 =  0*e1+1*e2+0*e3;
-  v110 =  1*e1+1*e2+0*e3;
-  v210 = -1*e1+1*e2+0*e3;
-  v020 =  0*e1-1*e2+0*e3;
-  v120 =  1*e1-1*e2+0*e3;
-  v220 = -1*e1-1*e2+0*e3;
-  v001 =  0*e1+0*e2+1*e3;
-  v101 =  1*e1+0*e2+1*e3;
-  v201 = -1*e1+0*e2+1*e3;
-  v011 =  0*e1+1*e2+1*e3;
-  v111 =  1*e1+1*e2+1*e3;
-  v211 = -1*e1+1*e2+1*e3;
-  v021 =  0*e1-1*e2+1*e3;
-  v121 =  1*e1-1*e2+1*e3;
-  v221 = -1*e1-1*e2+1*e3;
-  v002 =  0*e1+0*e2-1*e3;
-  v102 =  1*e1+0*e2-1*e3;
-  v202 = -1*e1+0*e2-1*e3;
-  v012 =  0*e1+1*e2-1*e3;
-  v112 =  1*e1+1*e2-1*e3;
-  v212 = -1*e1+1*e2-1*e3;
-  v022 =  0*e1-1*e2-1*e3;
-  v122 =  1*e1-1*e2-1*e3;
-  v222 = -1*e1-1*e2-1*e3;
+  v000 =  0*unitVec[0]+0*unitVec[1]+0*unitVec[2];
+  v100 =  1*unitVec[0]+0*unitVec[1]+0*unitVec[2];
+  v200 = -1*unitVec[0]+0*unitVec[1]+0*unitVec[2];
+  v010 =  0*unitVec[0]+1*unitVec[1]+0*unitVec[2];
+  v110 =  1*unitVec[0]+1*unitVec[1]+0*unitVec[2];
+  v210 = -1*unitVec[0]+1*unitVec[1]+0*unitVec[2];
+  v020 =  0*unitVec[0]-1*unitVec[1]+0*unitVec[2];
+  v120 =  1*unitVec[0]-1*unitVec[1]+0*unitVec[2];
+  v220 = -1*unitVec[0]-1*unitVec[1]+0*unitVec[2];
+  v001 =  0*unitVec[0]+0*unitVec[1]+1*unitVec[2];
+  v101 =  1*unitVec[0]+0*unitVec[1]+1*unitVec[2];
+  v201 = -1*unitVec[0]+0*unitVec[1]+1*unitVec[2];
+  v011 =  0*unitVec[0]+1*unitVec[1]+1*unitVec[2];
+  v111 =  1*unitVec[0]+1*unitVec[1]+1*unitVec[2];
+  v211 = -1*unitVec[0]+1*unitVec[1]+1*unitVec[2];
+  v021 =  0*unitVec[0]-1*unitVec[1]+1*unitVec[2];
+  v121 =  1*unitVec[0]-1*unitVec[1]+1*unitVec[2];
+  v221 = -1*unitVec[0]-1*unitVec[1]+1*unitVec[2];
+  v002 =  0*unitVec[0]+0*unitVec[1]-1*unitVec[2];
+  v102 =  1*unitVec[0]+0*unitVec[1]-1*unitVec[2];
+  v202 = -1*unitVec[0]+0*unitVec[1]-1*unitVec[2];
+  v012 =  0*unitVec[0]+1*unitVec[1]-1*unitVec[2];
+  v112 =  1*unitVec[0]+1*unitVec[1]-1*unitVec[2];
+  v212 = -1*unitVec[0]+1*unitVec[1]-1*unitVec[2];
+  v022 =  0*unitVec[0]-1*unitVec[1]-1*unitVec[2];
+  v122 =  1*unitVec[0]-1*unitVec[1]-1*unitVec[2];
+  v222 = -1*unitVec[0]-1*unitVec[1]-1*unitVec[2];
 }
 
@@ -129 +129 @@
   CPPUNIT_ASSERT(Everywhere().isInside(v122));
   CPPUNIT_ASSERT(Everywhere().isInside(v222));
+
+  Shape s1 = Sphere();
+
+  CPPUNIT_ASSERT(s1.isInside(v000));
+  CPPUNIT_ASSERT(s1.isInside(v100));
+  CPPUNIT_ASSERT(s1.isInside(v200));
+  CPPUNIT_ASSERT(s1.isInside(v010));
+  CPPUNIT_ASSERT(!s1.isInside(v110));
+  CPPUNIT_ASSERT(!s1.isInside(v210));
+  CPPUNIT_ASSERT(s1.isInside(v020));
+  CPPUNIT_ASSERT(!s1.isInside(v120));
+  CPPUNIT_ASSERT(!s1.isInside(v220));
+  CPPUNIT_ASSERT(s1.isInside(v001));
+  CPPUNIT_ASSERT(!s1.isInside(v101));
+  CPPUNIT_ASSERT(!s1.isInside(v201));
+  CPPUNIT_ASSERT(!s1.isInside(v011));
+  CPPUNIT_ASSERT(!s1.isInside(v111));
+  CPPUNIT_ASSERT(!s1.isInside(v211));
+  CPPUNIT_ASSERT(!s1.isInside(v021));
+  CPPUNIT_ASSERT(!s1.isInside(v121));
+  CPPUNIT_ASSERT(!s1.isInside(v221));
+  CPPUNIT_ASSERT(s1.isInside(v002));
+  CPPUNIT_ASSERT(!s1.isInside(v102));
+  CPPUNIT_ASSERT(!s1.isInside(v202));
+  CPPUNIT_ASSERT(!s1.isInside(v012));
+  CPPUNIT_ASSERT(!s1.isInside(v112));
+  CPPUNIT_ASSERT(!s1.isInside(v212));
+  CPPUNIT_ASSERT(!s1.isInside(v022));
+  CPPUNIT_ASSERT(!s1.isInside(v122));
+  CPPUNIT_ASSERT(!s1.isInside(v222));
+
+  Shape s2 = Cuboid();
+
+  CPPUNIT_ASSERT(s2.isInside(v000));
+  CPPUNIT_ASSERT(s2.isInside(v100));
+  CPPUNIT_ASSERT(!s2.isInside(v200));
+  CPPUNIT_ASSERT(s2.isInside(v010));
+  CPPUNIT_ASSERT(s2.isInside(v110));
+  CPPUNIT_ASSERT(!s2.isInside(v210));
+  CPPUNIT_ASSERT(!s2.isInside(v020));
+  CPPUNIT_ASSERT(!s2.isInside(v120));
+  CPPUNIT_ASSERT(!s2.isInside(v220));
+  CPPUNIT_ASSERT(s2.isInside(v001));
+  CPPUNIT_ASSERT(s2.isInside(v101));
+  CPPUNIT_ASSERT(!s2.isInside(v201));
+  CPPUNIT_ASSERT(s2.isInside(v011));
+  CPPUNIT_ASSERT(s2.isInside(v111));
+  CPPUNIT_ASSERT(!s2.isInside(v211));
+  CPPUNIT_ASSERT(!s2.isInside(v021));
+  CPPUNIT_ASSERT(!s2.isInside(v121));
+  CPPUNIT_ASSERT(!s2.isInside(v221));
+  CPPUNIT_ASSERT(!s2.isInside(v002));
+  CPPUNIT_ASSERT(!s2.isInside(v102));
+  CPPUNIT_ASSERT(!s2.isInside(v202));
+  CPPUNIT_ASSERT(!s2.isInside(v012));
+  CPPUNIT_ASSERT(!s2.isInside(v112));
+  CPPUNIT_ASSERT(!s2.isInside(v212));
+  CPPUNIT_ASSERT(!s2.isInside(v022));
+  CPPUNIT_ASSERT(!s2.isInside(v122));
+  CPPUNIT_ASSERT(!s2.isInside(v222));
+}
+
+void ShapeUnittest::surfaceTest(){
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v000));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v100));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v200));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v010));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v110));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v210));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v020));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v120));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v220));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v001));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v101));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v201));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v011));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v111));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v211));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v021));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v121));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v221));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v002));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v102));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v202));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v012));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v112));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v212));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v022));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v122));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v222));
+
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v000));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v100));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v200));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v010));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v110));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v210));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v020));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v120));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v220));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v001));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v101));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v201));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v011));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v111));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v211));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v021));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v121));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v221));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v002));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v102));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v202));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v012));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v112));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v212));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v022));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v122));
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v222));
+
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v000));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v100));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v200));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v010));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v110));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v210));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v020));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v120));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v220));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v001));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v101));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v201));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v011));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v111));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v211));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v021));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v121));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v221));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v002));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v102));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v202));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v012));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v112));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v212));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v022));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v122));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v222));
 }
 
@@ -377 +523 @@
   }
 
+  {
+    Shape s1 = Sphere() && Cuboid(); // Sphere in upper place
+
+    CPPUNIT_ASSERT(s1.isInside(v000));
+    CPPUNIT_ASSERT(s1.isInside(v100));
+    CPPUNIT_ASSERT(!s1.isInside(v200));
+    CPPUNIT_ASSERT(s1.isInside(v010));
+    CPPUNIT_ASSERT(!s1.isInside(v110));
+    CPPUNIT_ASSERT(!s1.isInside(v210));
+    CPPUNIT_ASSERT(!s1.isInside(v020));
+    CPPUNIT_ASSERT(!s1.isInside(v120));
+    CPPUNIT_ASSERT(!s1.isInside(v220));
+    CPPUNIT_ASSERT(s1.isInside(v001));
+    CPPUNIT_ASSERT(!s1.isInside(v101));
+    CPPUNIT_ASSERT(!s1.isInside(v201));
+    CPPUNIT_ASSERT(!s1.isInside(v011));
+    CPPUNIT_ASSERT(!s1.isInside(v111));
+    CPPUNIT_ASSERT(!s1.isInside(v211));
+    CPPUNIT_ASSERT(!s1.isInside(v021));
+    CPPUNIT_ASSERT(!s1.isInside(v121));
+    CPPUNIT_ASSERT(!s1.isInside(v221));
+    CPPUNIT_ASSERT(!s1.isInside(v002));
+    CPPUNIT_ASSERT(!s1.isInside(v102));
+    CPPUNIT_ASSERT(!s1.isInside(v202));
+    CPPUNIT_ASSERT(!s1.isInside(v012));
+    CPPUNIT_ASSERT(!s1.isInside(v112));
+    CPPUNIT_ASSERT(!s1.isInside(v212));
+    CPPUNIT_ASSERT(!s1.isInside(v022));
+    CPPUNIT_ASSERT(!s1.isInside(v122));
+    CPPUNIT_ASSERT(!s1.isInside(v222));
+
+    CPPUNIT_ASSERT(s1.isOnSurface(v000));
+    CPPUNIT_ASSERT(s1.isOnSurface(v100));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v200));
+    CPPUNIT_ASSERT(s1.isOnSurface(v010));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v110));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v210));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v020));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v120));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v220));
+    CPPUNIT_ASSERT(s1.isOnSurface(v001));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v101));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v201));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v011));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v111));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v211));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v021));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v121));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v221));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v002));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v102));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v202));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v012));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v112));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v212));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v022));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v122));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v222));
+
+
+    Shape s2 = Sphere() || Cuboid(); // Should be same as Cuboid
+
+    CPPUNIT_ASSERT(s2.isInside(v000));
+    CPPUNIT_ASSERT(s2.isInside(v100));
+    CPPUNIT_ASSERT(s2.isInside(v200));
+    CPPUNIT_ASSERT(s2.isInside(v010));
+    CPPUNIT_ASSERT(s2.isInside(v110));
+    CPPUNIT_ASSERT(!s2.isInside(v210));
+    CPPUNIT_ASSERT(s2.isInside(v020));
+    CPPUNIT_ASSERT(!s2.isInside(v120));
+    CPPUNIT_ASSERT(!s2.isInside(v220));
+    CPPUNIT_ASSERT(s2.isInside(v001));
+    CPPUNIT_ASSERT(s2.isInside(v101));
+    CPPUNIT_ASSERT(!s2.isInside(v201));
+    CPPUNIT_ASSERT(s2.isInside(v011));
+    CPPUNIT_ASSERT(s2.isInside(v111));
+    CPPUNIT_ASSERT(!s2.isInside(v211));
+    CPPUNIT_ASSERT(!s2.isInside(v021));
+    CPPUNIT_ASSERT(!s2.isInside(v121));
+    CPPUNIT_ASSERT(!s2.isInside(v221));
+    CPPUNIT_ASSERT(s2.isInside(v002));
+    CPPUNIT_ASSERT(!s2.isInside(v102));
+    CPPUNIT_ASSERT(!s2.isInside(v202));
+    CPPUNIT_ASSERT(!s2.isInside(v012));
+    CPPUNIT_ASSERT(!s2.isInside(v112));
+    CPPUNIT_ASSERT(!s2.isInside(v212));
+    CPPUNIT_ASSERT(!s2.isInside(v022));
+    CPPUNIT_ASSERT(!s2.isInside(v122));
+    CPPUNIT_ASSERT(!s2.isInside(v222));
+
+    CPPUNIT_ASSERT(!s2.isOnSurface(v000));
+    CPPUNIT_ASSERT(s2.isOnSurface(v100));
+    CPPUNIT_ASSERT(s2.isOnSurface(v200));
+    CPPUNIT_ASSERT(s2.isOnSurface(v010));
+    CPPUNIT_ASSERT(s2.isOnSurface(v110));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v210));
+    CPPUNIT_ASSERT(s2.isOnSurface(v020));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v120));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v220));
+    CPPUNIT_ASSERT(s2.isOnSurface(v001));
+    CPPUNIT_ASSERT(s2.isOnSurface(v101));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v201));
+    CPPUNIT_ASSERT(s2.isOnSurface(v011));
+    CPPUNIT_ASSERT(s2.isOnSurface(v111));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v211));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v021));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v121));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v221));
+    CPPUNIT_ASSERT(s2.isOnSurface(v002));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v102));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v202));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v012));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v112));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v212));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v022));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v122));
+    CPPUNIT_ASSERT(!s2.isOnSurface(v222));
+  }
+
 }
+
+void ShapeUnittest::PointsOnSurfaceTest(){
+  Shape s = Sphere();
+  const size_t N = 200;
+  std::vector<Vector> PointsOnSurface = s.getHomogeneousPointsOnSurface(N);
+  for (std::vector<Vector>::const_iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
+    CPPUNIT_ASSERT(fabs(1. - (*iter).NormSquared()) < MYEPSILON);
+  }
+  CPPUNIT_ASSERT_EQUAL(N, PointsOnSurface.size());
+}
+

src/unittests/ShapeUnittest.hpp

@@ -20 +20 @@
   CPPUNIT_TEST_SUITE( ShapeUnittest) ;
   CPPUNIT_TEST ( baseShapesTest );
+  CPPUNIT_TEST ( surfaceTest );
   CPPUNIT_TEST ( assignmentTest );
   CPPUNIT_TEST ( operatorTest );
+  CPPUNIT_TEST ( PointsOnSurfaceTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -29 +31 @@
 
   void baseShapesTest();
+  void surfaceTest();
   void assignmentTest();
   void operatorTest();
+  void PointsOnSurfaceTest();
 
   // a lot of vectors.