Changeset a7c344

Timestamp:

Jun 19, 2010, 4:06:59 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

c39cc4

Parents:

b32dbb (diff), 27ac00 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'StructureRefactoring' of jupiter:espack into StructureRefactoring

Location:

src

Files:

• Exceptions/MultipleSolutionsException.hpp (added)
• Exceptions/SkewException.cpp (added)
• Exceptions/SkewException.hpp (added)
• Legacy/oldmenu.cpp (modified) (3 diffs)
• Line.cpp (modified) (2 diffs)
• Line.hpp (modified) (2 diffs)
• Makefile.am (modified) (2 diffs)
• Plane.cpp (modified) (4 diffs)
• Plane.hpp (modified) (2 diffs)
• Space.cpp (modified) (1 diff)
• Space.hpp (modified) (1 diff)
• atom.cpp (modified) (1 diff)
• boundary.cpp (modified) (1 diff)
• molecule_geometry.cpp (modified) (2 diffs)
• periodentafel.cpp (modified) (6 diffs)
• tesselation.cpp (modified) (9 diffs)
• tesselation.hpp (modified) (2 diffs)
• tesselationhelpers.cpp (modified) (2 diffs)
• unittests/LineUnittest.cpp (added)
• unittests/LineUnittest.hpp (added)
• unittests/Makefile.am (modified) (4 diffs)
• unittests/PlaneUnittest.cpp (modified) (2 diffs)
• unittests/PlaneUnittest.hpp (modified) (2 diffs)
• unittests/vectorunittest.cpp (modified) (1 diff)
• unittests/vectorunittest.hpp (modified) (1 diff)
• vector.cpp (modified) (5 diffs)
• vector.hpp (modified) (5 diffs)
• vector_ops.cpp (modified) (2 diffs)
• vector_ops.hpp (modified) (1 diff)

src/Legacy/oldmenu.cpp

@@ -27 +27 @@
 #include "vector_ops.hpp"
 #include "Plane.hpp"
+#include "Line.hpp"
 
 #include "UIElements/UIFactory.hpp"
@@ -185 +186 @@
           // rotate vector around first angle
           first->x = x;
-          first->x = RotateVector(first->x,z,b - M_PI);
+          first->x = Line(zeroVec,z).rotateVector(first->x,b - M_PI);
           Log() << Verbose(0) << "Rotated vector: " << first->x << endl,
           // remove the projection onto the rotation plane of the second angle
@@ -201 +202 @@
           // rotate another vector around second angle
           n = y;
-          n = RotateVector(n,x,c - M_PI);
+          n = Line(zeroVec,x).rotateVector(n,c - M_PI);
           Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;
 

src/Line.cpp

@@ -11 +11 @@
 
 #include "vector.hpp"
-
-Line::Line(Vector &_origin, Vector &_direction) :
-  origin(new Vector(_origin)),
+#include "log.hpp"
+#include "verbose.hpp"
+#include "gslmatrix.hpp"
+#include "info.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
+#include "Exceptions/SkewException.hpp"
+#include "Plane.hpp"
+
+using namespace std;
+
+Line::Line(const Vector &_origin, const Vector &_direction) :
   direction(new Vector(_direction))
 {
   direction->Normalize();
-}
+  origin.reset(new Vector(_origin.partition(*direction).second));
+}
+
+Line::Line(const Line &src) :
+  origin(new Vector(*src.origin)),
+  direction(new Vector(*src.direction))
+{}
 
 Line::~Line()
@@ -24 +38 @@
 
 double Line::distance(const Vector &point) const{
-  return fabs(point.ScalarProduct(*direction) - origin->ScalarProduct(*direction));
+  // get any vector from line to point
+  Vector helper = point - *origin;
+  // partition this vector along direction
+  // the residue points from the line to the point
+  return helper.partition(*direction).second.Norm();
 }
 
 Vector Line::getClosestPoint(const Vector &point) const{
-  double factor = point.ScalarProduct(*direction) - origin->ScalarProduct(*direction);
-  return (point - (factor * (*direction)));
-}
+  // get any vector from line to point
+  Vector helper = point - *origin;
+  // partition this vector along direction
+  // add only the part along the direction
+  return *origin + helper.partition(*direction).first;
+}
+
+Vector Line::getDirection() const{
+  return *direction;
+}
+
+Vector Line::getOrigin() const{
+  return *origin;
+}
+
+vector<Vector> Line::getPointsOnLine() const{
+  vector<Vector> res;
+  res.reserve(2);
+  res.push_back(*origin);
+  res.push_back(*origin+*direction);
+  return res;
+}
+
+/** Calculates the intersection of the two lines that are both on the same plane.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
+ * \param *out output stream for debugging
+ * \param *Line1a first vector of first line
+ * \param *Line1b second vector of first line
+ * \param *Line2a first vector of second line
+ * \param *Line2b second vector of second line
+ * \return true - \a this will contain the intersection on return, false - lines are parallel
+ */
+Vector Line::getIntersection(const Line& otherLine) const{
+  Info FunctionInfo(__func__);
+
+  pointset line1Points = getPointsOnLine();
+
+  Vector Line1a = line1Points[0];
+  Vector Line1b = line1Points[1];
+
+  pointset line2Points = otherLine.getPointsOnLine();
+
+  Vector Line2a = line2Points[0];
+  Vector Line2b = line2Points[1];
+
+  Vector res;
+
+  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
+
+  M->SetAll(1.);
+  for (int i=0;i<3;i++) {
+    M->Set(0, i, Line1a[i]);
+    M->Set(1, i, Line1b[i]);
+    M->Set(2, i, Line2a[i]);
+    M->Set(3, i, Line2b[i]);
+  }
+
+  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  //for (int i=0;i<4;i++) {
+  //  for (int j=0;j<4;j++)
+  //    cout << "\t" << M->Get(i,j);
+  //  cout << endl;
+  //}
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
+    throw SkewException(__FILE__,__LINE__);
+  }
+
+  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
+
+
+  // constuct a,b,c
+  Vector a = Line1b - Line1a;
+  Vector b = Line2b - Line2a;
+  Vector c = Line2a - Line1a;
+  Vector d = Line2b - Line1b;
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
+   res.Zero();
+   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  // check for parallelity
+  Vector parallel;
+  double factor = 0.;
+  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+    parallel = Line1a - Line2a;
+    factor = parallel.ScalarProduct(a)/a.Norm();
+    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+      res = Line2a;
+      Log() << Verbose(1) << "Lines conincide." << endl;
+      return res;
+    } else {
+      parallel = Line1a - Line2b;
+      factor = parallel.ScalarProduct(a)/a.Norm();
+      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+        res = Line2b;
+        Log() << Verbose(1) << "Lines conincide." << endl;
+        return res;
+      }
+    }
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    res.Zero();
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1 = c;
+  temp1.VectorProduct(b);
+  temp2 = a;
+  temp2.VectorProduct(b);
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+
+  // construct intersection
+  res = a;
+  res.Scale(s);
+  res += Line1a;
+  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
+
+  return res;
+}
+
+/** Rotates the vector by an angle of \a alpha around this line.
+ * \param rhs Vector to rotate
+ * \param alpha rotation angle in radian
+ */
+Vector Line::rotateVector(const Vector &rhs, double alpha) const{
+  Vector helper = rhs;
+
+  // translate the coordinate system so that the line goes through (0,0,0)
+  helper -= *origin;
+
+  // partition the vector into a part that gets rotated and a part that lies along the line
+  pair<Vector,Vector> parts = helper.partition(*direction);
+
+  // we just keep anything that is along the axis
+  Vector res = parts.first;
+
+  // the rest has to be rotated
+  Vector a = parts.second;
+  // we only have to do the rest, if we actually could partition the vector
+  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);
+
+    res += cos(alpha) * a + sin(alpha) * y;
+  }
+
+  // translate the coordinate system back
+  res += *origin;
+  return res;
+}
+
+Plane Line::getOrthogonalPlane(const Vector &origin) const{
+  return Plane(getDirection(),origin);
+}
+
+Line makeLineThrough(const Vector &x1, const Vector &x2){
+  if(x1==x2){
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+  return Line(x1,x1-x2);
+}

src/Line.hpp

@@ -12 +12 @@
 
 #include <memory>
+#include <vector>
 
 class Vector;
+class Plane;
 
 class Line : public Space
 {
 public:
-  Line(Vector &_origin, Vector &_direction);
+  Line(const Vector &_origin, const Vector &_direction);
+  Line(const Line& _src);
   virtual ~Line();
 
-  virtual double distance(const Vector &point) const=0;
-  virtual Vector getClosestPoint(const Vector &point) const=0;
+  virtual double distance(const Vector &point) const;
+  virtual Vector getClosestPoint(const Vector &point) const;
+
+  Vector getDirection() const;
+  Vector getOrigin() const;
+
+  std::vector<Vector> getPointsOnLine() const;
+
+  Vector getIntersection(const Line& otherLine) const;
+
+  Vector rotateVector(const Vector &rhs, double alpha) const;
+
+  Plane getOrthogonalPlane(const Vector &origin) const;
 
 private:
@@ -29 +43 @@
 };
 
+/**
+ * Named constructor to make a line through two points
+ */
+Line makeLineThrough(const Vector &x1, const Vector &x2);
+
 #endif /* LINE_HPP_ */

src/Makefile.am

@@ -88 +88 @@
                                   Exceptions/LinearDependenceException.cpp \
                                   Exceptions/MathException.cpp \
+                                  Exceptions/SkewException.cpp \
                                   Exceptions/ZeroVectorException.cpp
                                   
@@ -93 +94 @@
                                   Exceptions/LinearDependenceException.hpp \
                                   Exceptions/MathException.hpp \
+                                  Exceptions/SkewException.hpp \
                                   Exceptions/ZeroVectorException.hpp
 

src/Plane.cpp

@@ -14 +14 @@
 #include "Helpers/Assert.hpp"
 #include <cmath>
+#include "Line.hpp"
+#include "Exceptions/MultipleSolutionsException.hpp"
 
 /**
@@ -42 +44 @@
 /**
  * Constructs a plane from two direction vectors and a offset.
- * If no offset is given a plane through origin is assumed
  */
 Plane::Plane(const Vector &y1, const Vector &y2, double _offset) throw(ZeroVectorException,LinearDependenceException) :
@@ -113 +114 @@
 }
 
-Vector Plane::getOffsetVector() {
+Vector Plane::getOffsetVector() const {
   return getOffset()*getNormal();
 }
 
-vector<Vector> Plane::getPointsOnPlane(){
+vector<Vector> Plane::getPointsOnPlane() const{
   std::vector<Vector> res;
   res.reserve(3);
@@ -147 +148 @@
  * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
  */
-Vector Plane::GetIntersection(const Vector &Origin, const Vector &LineVector)
+Vector Plane::GetIntersection(const Line& line) const
 {
   Info FunctionInfo(__func__);
   Vector res;
 
-  // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
-  Vector Direction = LineVector - Origin;
-  Direction.Normalize();
-  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
-  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
-  double factor1 = Direction.ScalarProduct(*normalVector.get());
-  if (fabs(factor1) < MYEPSILON) { // Uniqueness: line parallel to plane?
-    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  double factor2 = Origin.ScalarProduct(*normalVector.get());
-  if (fabs(factor2-offset) < MYEPSILON) { // Origin is in-plane
-    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
-    res = Origin;
-    return res;
-  }
-
+  double factor1 = getNormal().ScalarProduct(line.getDirection());
+  if(fabs(factor1)<MYEPSILON){
+    // the plane is parallel... under all circumstances this is bad luck
+    // we no have either no or infinite solutions
+    if(isContained(line.getOrigin())){
+      throw MultipleSolutionsException<Vector>(__FILE__,__LINE__,line.getOrigin());
+    }
+    else{
+      throw LinearDependenceException(__FILE__,__LINE__);
+    }
+  }
+
+  double factor2 = getNormal().ScalarProduct(line.getOrigin());
   double scaleFactor = (offset-factor2)/factor1;
 
-  //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
-  Direction.Scale(scaleFactor);
-  res = Origin + Direction;
-  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
-
-  // test whether resulting vector really is on plane
-  ASSERT(fabs(res.ScalarProduct(*normalVector) - offset) < MYEPSILON,
-         "Calculated line-Plane intersection does not lie on plane.");
+  res = line.getOrigin() + scaleFactor * line.getDirection();
+
+  // tests to make sure the resulting vector really is on plane and line
+  ASSERT(isContained(res),"Calculated line-Plane intersection does not lie on plane.");
+  ASSERT(line.isContained(res),"Calculated line-Plane intersection does not lie on line.");
   return res;
 };
+
+Vector Plane::mirrorVector(const Vector &rhs) const {
+  Vector helper = getVectorToPoint(rhs);
+  // substract twice the Vector to the plane
+  return rhs+2*helper;
+}
+
+Line Plane::getOrthogonalLine(const Vector &origin) const{
+  return Line(origin,getNormal());
+}
 
 /************ Methods inherited from Space ****************/

src/Plane.hpp

@@ -17 +17 @@
 
 class Vector;
+class Line;
 
 class Plane : public Space
@@ -42 +43 @@
    * Same as getOffset()*getNormal();
    */
-  Vector getOffsetVector();
+  Vector getOffsetVector() const;
 
   /**
    * returns three seperate points on this plane
    */
-  std::vector<Vector> getPointsOnPlane();
+  std::vector<Vector> getPointsOnPlane() const;
 
   // some calculations
-  Vector GetIntersection(const Vector &Origin, const Vector &LineVector);
+  Vector GetIntersection(const Line &Line) const;
+
+  Vector mirrorVector(const Vector &rhs) const;
+
+  /**
+   * get a Line that is orthogonal to this plane, going through a chosen
+   * point.
+   *
+   * The point does not have to lie on the plane itself.
+   */
+  Line getOrthogonalLine(const Vector &origin) const;
 
   /****** Methods inherited from Space ***********/

src/Space.cpp

@@ -19 +19 @@
 }
 
+Vector Space::getVectorToPoint(const Vector &origin) const{
+  Vector support = getClosestPoint(origin);
+  return support-origin;
+}
+
 bool Space::isContained(const Vector &point) const{
   return (distance(point)) < MYEPSILON;

src/Space.hpp

@@ -17 +17 @@
   virtual ~Space();
 
+  /**
+   * Calculates the distance between a Space and a Vector.
+   */
   virtual double distance(const Vector &point) const=0;
+
+  /**
+   * get the closest point inside the space to another point
+   */
   virtual Vector getClosestPoint(const Vector &point) const=0;
+
+  /**
+   * get the shortest Vector from a point to a Space.
+   *
+   * The Vector always points from the given Vector to the point in space
+   * returned by Plane::getClosestPoint().
+   */
+  virtual Vector getVectorToPoint(const Vector &point) const;
+
+  /**
+   * Test wether a point is contained in the space.
+   *
+   * returns true, when the point lies inside and false
+   * otherwise.
+   */
   virtual bool isContained(const Vector &point) const;
+
+  /**
+   * Tests if this space contains the center of the coordinate system.
+   */
   virtual bool hasZero() const;
 

src/atom.cpp

@@ -68 +68 @@
 atom *atom::GetTrueFather()
 {
-  atom *walker = this;
-  do {
-    if (walker == walker->father) // top most father is the one that points on itself
-      break;
-    walker = walker->father;
-  } while (walker != NULL);
-  return walker;
+  if(father == this){ // top most father is the one that points on itself
+    return this;
+  }
+  else if(!father) {
+    return 0;
+  }
+  else {
+    return father->GetTrueFather();
+  }
 };
 

src/boundary.cpp

@@ -620 +620 @@
   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
     { // go through every triangle, calculate volume of its pyramid with CoG as peak
-      x = (*runner->second->endpoints[0]->node->node) - (*runner->second->endpoints[1]->node->node);
-      y = (*runner->second->endpoints[0]->node->node) - (*runner->second->endpoints[2]->node->node);
-      const double a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(*runner->second->endpoints[1]->node->node));
-      const double b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(*runner->second->endpoints[2]->node->node));
-      const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(*runner->second->endpoints[1]->node->node));
+      x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1);
+      y = runner->second->getEndpoint(0) - runner->second->getEndpoint(2);
+      const double a = x.Norm();
+      const double b = y.Norm();
+      const double c = runner->second->getEndpoint(2).distance(runner->second->getEndpoint(1));
       const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
-      x = Plane(*(runner->second->endpoints[0]->node->node),
-                *(runner->second->endpoints[1]->node->node),
-                *(runner->second->endpoints[2]->node->node)).getNormal();
-      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(x));
+      x = runner->second->getPlane().getNormal();
+      x.Scale(runner->second->getEndpoint(1).ScalarProduct(x));
       const double h = x.Norm(); // distance of CoG to triangle
       const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)

src/molecule_geometry.cpp

@@ -16 +16 @@
 #include "molecule.hpp"
 #include "World.hpp"
+#include "Plane.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -251 +252 @@
 void molecule::Mirror(const Vector *n)
 {
-  ActOnAllVectors( &Vector::Mirror, *n );
+  Plane p(*n,0);
+  // TODO: replace with simpler construct (e.g. Boost::foreach)
+  // once the structure of the atom list is fully reworked
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (*Walker->node) = p.mirrorVector(*Walker->node);
+  }
 };
 

src/periodentafel.cpp

@@ -212 +212 @@
 
   // fill elements DB
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDELEMENTSDB, MAXSTRINGSIZE-strlen(filename));
+  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDELEMENTSDB);
   infile.open(filename);
   if (infile != NULL) {
@@ -258 +256 @@
 
   // fill valence DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
+  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDVALENCEDB);
   infile.open(filename);
   if (infile != NULL) {
@@ -277 +273 @@
 
   // fill valence DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
+  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDORBITALDB);
   infile.open(filename);
   if (infile != NULL) {
@@ -296 +290 @@
 
   // fill H-BondDistance DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
+  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDHBONDDISTANCEDB);
   infile.open(filename);
   if (infile != NULL) {
@@ -318 +310 @@
 
   // fill H-BondAngle DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
+  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDHBONDANGLEDB);
   infile.open(filename);
   if (infile != NULL) {
@@ -363 +353 @@
   char filename[MAXSTRINGSIZE];
 
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDELEMENTSDB, MAXSTRINGSIZE-strlen(filename));
+  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDELEMENTSDB);
   f.open(filename);
   if (f != NULL) {

src/tesselation.cpp

@@ -17 +17 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -469 +470 @@
 
   try {
-    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
-  }
-  catch (LinearDependenceException &excp) {
-    Log() << Verbose(1) << excp;
-    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
-    return false;
-  }
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
-
-  if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
-    return true;
-  }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
-    return true;
-  }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
-    return true;
-  }
-  // Calculate cross point between one baseline and the line from the third endpoint to intersection
-  int i = 0;
-  do {
-    try {
-      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
-                                                    *(endpoints[(i+1)%3]->node->node),
-                                                    *(endpoints[(i+2)%3]->node->node),
-                                                    *Intersection);
+    Line centerLine = makeLineThrough(*MolCenter, *x);
+    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);
+
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
+
+    if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
+      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
+      return true;
+    }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
+      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
+      return true;
+    }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
+      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
+      return true;
+    }
+    // Calculate cross point between one baseline and the line from the third endpoint to intersection
+    int i = 0;
+    do {
+      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
+      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
+      CrossPoint = line1.getIntersection(line2);
       helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
       CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
@@ -505 +499 @@
       if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
         DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
-        i=4;
-        break;
+        return false;
       }
       i++;
-    } catch (LinearDependenceException &excp){
-      break;
-    }
-  } while (i < 3);
-  if (i == 3) {
+    } while (i < 3);
     DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
     return true;
-  } else {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl);
+  }
+  catch (MathException &excp) {
+    Log() << Verbose(1) << excp;
+    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
     return false;
   }
@@ -544 +535 @@
   GetCenter(&Direction);
   try {
-    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
-  }
-  catch (LinearDependenceException &excp) {
+    Line l = makeLineThrough(*x, Direction);
+    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
+  }
+  catch (MathException &excp) {
     (*ClosestPoint) = (*x);
   }
@@ -569 +561 @@
     Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
+    Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
+    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
     CrossDirection[i] = CrossPoint[i] - InPlane;
     CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
@@ -740 +733 @@
 }
 
+Vector BoundaryTriangleSet::getEndpoint(int i) const{
+  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
+
+  return *endpoints[i]->node->node;
+}
+
+string BoundaryTriangleSet::getEndpointName(int i) const{
+  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
+
+  return endpoints[i]->node->getName();
+}
+
 /** output operator for BoundaryTriangleSet.
  * \param &ost output stream
@@ -746 +751 @@
 ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << "," << a.endpoints[1]->node->getName() << "," << a.endpoints[2]->node->getName() << "]";
+  ost << "[" << a.Nr << "|" << a.getEndpointName(0) << "," << a.getEndpointName(1) << "," << a.getEndpointName(2) << "]";
   //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
   //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
@@ -1514 +1519 @@
         CenterVector.Zero();
         for (int i = 0; i < 3; i++)
-          CenterVector += (*BTS->endpoints[i]->node->node);
+          CenterVector += BTS->getEndpoint(i);
         CenterVector.Scale(1. / 3.);
         DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);
@@ -4868 +4873 @@
   if (LastTriangle != NULL) {
     stringstream sstr;
-    sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->endpoints[0]->node->getName() << "_" << LastTriangle->endpoints[1]->node->getName() << "_" << LastTriangle->endpoints[2]->node->getName();
+    sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
     NumberName = sstr.str();
     if (DoTecplotOutput) {

src/tesselation.hpp

@@ -171 +171 @@
 
     Plane getPlane() const;
+    Vector getEndpoint(int) const;
+    std::string getEndpointName(int) const;
 
     class BoundaryPointSet *endpoints[3];
@@ -177 +179 @@
     Vector SphereCenter;
     int Nr;
+
+  private:
+
 };
 

src/tesselationhelpers.cpp

@@ -15 +15 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
+#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -724 +725 @@
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
-  *Intersection = GetIntersectionOfTwoLinesOnPlane(*(Base->endpoints[0]->node->node),
-                                                   *(Base->endpoints[1]->node->node),
-                                                     NewOffset, NewDirection);
+  Line line1 = makeLineThrough(*(Base->endpoints[0]->node->node),*(Base->endpoints[1]->node->node));
+  Line line2 = makeLineThrough(NewOffset, NewDirection);
+  *Intersection = line1.getIntersection(line2);
   Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
   DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);

src/unittests/Makefile.am

@@ -23 +23 @@
   InfoUnitTest \
   LinearSystemOfEquationsUnitTest \
+  LineUnittest \
   LinkedCellUnitTest \
   ListOfBondsUnitTest \
@@ -64 +65 @@
   infounittest.cpp \
   linearsystemofequationsunittest.cpp \
+  LineUnittest.cpp \
   LinkedCellUnitTest.cpp \
   listofbondsunittest.cpp \
@@ -97 +99 @@
   infounittest.hpp \
   linearsystemofequationsunittest.hpp \
+  LineUnittest.hpp \
   LinkedCellUnitTest.hpp \
   listofbondsunittest.hpp \
@@ -153 +156 @@
 LinearSystemOfEquationsUnitTest_LDADD = ${ALLLIBS}
 
+LineUnittest_SOURCES = UnitTestMain.cpp LineUnittest.cpp LineUnittest.hpp
+LineUnittest_LDADD = ${ALLLIBS}
+
 LinkedCellUnitTest_SOURCES = UnitTestMain.cpp LinkedCellUnitTest.cpp LinkedCellUnitTest.hpp 
 LinkedCellUnitTest_LDADD = ${ALLLIBS}

src/unittests/PlaneUnittest.cpp

@@ -17 +17 @@
 
 #include "vector.hpp"
+#include "Line.hpp"
 
 CPPUNIT_TEST_SUITE_REGISTRATION( PlaneUnittest );
@@ -153 +154 @@
   CPPUNIT_ASSERT(fabs(p4->distance(e1)-1) < MYEPSILON);
   CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(e1));
-
-
-}
+}
+
+void PlaneUnittest::mirrorTest(){
+  Vector fixture;
+
+  // 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 = 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 = 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 = 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);
+
+  // some Vectors outside of the planes
+  {
+    Vector t = (2./3.)*(e1+e2+e3);
+    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);
+}
+
+void PlaneUnittest::LineIntersectionTest(){
+  Vector fixture;
+  // 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_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) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
+}

src/unittests/PlaneUnittest.hpp

@@ -20 +20 @@
   CPPUNIT_TEST ( pointsTest );
   CPPUNIT_TEST ( operationsTest );
+  CPPUNIT_TEST ( mirrorTest );
+  CPPUNIT_TEST ( LineIntersectionTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -30 +32 @@
   void pointsTest();
   void operationsTest();
+  void mirrorTest();
+  void LineIntersectionTest();
 
 private:

src/unittests/vectorunittest.cpp

@@ -216 +216 @@
 }
 
-/** UnitTest for line intersections.
- */
-void VectorTest::LineIntersectionTest()
-{
-  // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unit, zero).GetIntersection(zero, two) );
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(otherunit, two).GetIntersection( unit, notunit) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
-
-  // four vectors equal to zero
-  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(zero, zero, zero, zero), LinearDependenceException);
-  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // four vectors equal to unit
-  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, unit, unit, unit), LinearDependenceException);
-  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // two equal lines
-  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, two));
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // line from (1,0,0) to (2,1,0) cuts line from (1,0,0) to (0,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, otherunit) );
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // line from (1,0,0) to (0,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, zero, zero, two) );
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // line from (1,0,0) to (2,1,0) cuts line from (0,0,0) to (0,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, zero, otherunit) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0., -1., 0.), fixture );
-};
-
-/** UnitTest for vector rotations.
- */
-void VectorTest::VectorRotationTest()
-{
-  fixture = Vector(-1.,0.,0.);
-
-  // zero vector does not change
-  fixture = RotateVector(zero,unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  fixture = RotateVector(zero, two, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero,  fixture);
-
-  // vector on axis does not change
-  fixture = RotateVector(unit,unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // rotations
-  fixture = RotateVector(otherunit, unit, M_PI);
-  CPPUNIT_ASSERT_EQUAL( Vector(0.,-1.,0.), fixture );
-
-  fixture = RotateVector(otherunit, unit, 2. * M_PI);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-
-  fixture = RotateVector(otherunit,unit, 0);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-
-  fixture = RotateVector(Vector(0.,0.,1.), notunit, M_PI);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-}
 
 /**

src/unittests/vectorunittest.hpp

@@ -27 +27 @@
     CPPUNIT_TEST ( ProjectionTest );
     CPPUNIT_TEST ( NormalsTest );
-    CPPUNIT_TEST ( LineIntersectionTest );
-    CPPUNIT_TEST ( VectorRotationTest );
     CPPUNIT_TEST ( IsInParallelepipedTest );
     CPPUNIT_TEST_SUITE_END();

src/vector.cpp

@@ -213 +213 @@
 {
   Vector tmp;
-  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
-  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
-  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
+  tmp[0] = x[1]* y[2] - x[2]* y[1];
+  tmp[1] = x[2]* y[0] - x[0]* y[2];
+  tmp[2] = x[0]* y[1] - x[1]* y[0];
   (*this) = tmp;
 };
@@ -232 +232 @@
   *this -= tmp;
 };
-
-/** Calculates the minimum distance vector of this vector to the plane.
- * \param *out output stream for debugging
- * \param *PlaneNormal normal of plane
- * \param *PlaneOffset offset of plane
- * \return distance to plane
- * \return distance vector onto to plane
- */
-Vector Vector::GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
-{
-  Vector temp = (*this) - PlaneOffset;
-  temp.MakeNormalTo(PlaneNormal);
-  temp.Scale(-1.);
-  // then add connecting vector from plane to point
-  temp += (*this)-PlaneOffset;
-  double sign = temp.ScalarProduct(PlaneNormal);
-  if (fabs(sign) > MYEPSILON)
-    sign /= fabs(sign);
-  else
-    sign = 0.;
-
-  temp.Normalize();
-  temp.Scale(sign);
-  return temp;
-};
-
 
 /** Calculates the minimum distance of this vector to the plane.
@@ -551 +525 @@
   MatrixMultiplication(M);
 };
+
+std::pair<Vector,Vector> Vector::partition(const Vector &rhs) const{
+  double factor = ScalarProduct(rhs)/rhs.NormSquared();
+  Vector res= factor * rhs;
+  return make_pair(res,(*this)-res);
+}
+
+std::pair<pointset,Vector> Vector::partition(const pointset &points) const{
+  Vector helper = *this;
+  pointset res;
+  for(pointset::const_iterator iter=points.begin();iter!=points.end();++iter){
+    pair<Vector,Vector> currPart = helper.partition(*iter);
+    res.push_back(currPart.first);
+    helper = currPart.second;
+  }
+  return make_pair(res,helper);
+}
 
 /** Do a matrix multiplication.
@@ -611 +602 @@
 };
 
-/** Mirrors atom against a given plane.
- * \param n[] normal vector of mirror plane.
- */
-void Vector::Mirror(const Vector &n)
-{
-  double projection;
-  projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
-  // withdraw projected vector twice from original one
-  for (int i=NDIM;i--;)
-    at(i) -= 2.*projection*n[i];
-};
-
 /** Calculates orthonormal vector to one given vectors.
  * Just subtracts the projection onto the given vector from this vector.
@@ -633 +612 @@
   bool result = false;
   double factor = y1.ScalarProduct(*this)/y1.NormSquared();
-  Vector x1;
-  x1 = factor * y1;
+  Vector x1 = factor * y1;
   SubtractVector(x1);
   for (int i=NDIM;i--;)

src/vector.hpp

@@ -16 +16 @@
 
 #include <memory>
+#include <vector>
 
 #include "defs.hpp"
@@ -21 +22 @@
 
 /********************************************** declarations *******************************/
+
+class Vector;
+
+typedef std::vector<Vector> pointset;
 
 /** Single vector.
@@ -39 +44 @@
 
   double DistanceSquared(const Vector &y) const;
-  Vector GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
   double DistanceToSpace(const Space& space) const;
   double PeriodicDistance(const Vector &y, const double * const cell_size) const;
@@ -56 +60 @@
   void ProjectIt(const Vector &y);
   Vector Projection(const Vector &y) const;
-  void Mirror(const Vector &x);
   void ScaleAll(const double *factor);
   void Scale(const double factor);
@@ -66 +69 @@
   bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
   void WrapPeriodically(const double * const M, const double * const Minv);
+  std::pair<Vector,Vector> partition(const Vector&) const;
+  std::pair<pointset,Vector> partition(const pointset&) const;
 
   // Accessors ussually come in pairs... and sometimes even more than that

src/vector_ops.cpp

@@ -15 +15 @@
 #include "Helpers/fast_functions.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
+#include "Exceptions/SkewException.hpp"
 
 #include <gsl/gsl_linalg.h>
@@ -110 +111 @@
   return true;
 };
-
-/** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
- * \param *axis rotation axis
- * \param alpha rotation angle in radian
- */
-Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha)
-{
-  Vector a,y;
-  Vector res;
-  // normalise this vector with respect to axis
-  a = vec;
-  a.ProjectOntoPlane(axis);
-  // construct normal vector
-  try {
-    y = Plane(axis,a,0).getNormal();
-  }
-  catch (MathException &excp) {
-    // The normal vector cannot be created if there is linar dependency.
-    // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
-    return vec;
-  }
-  y.Scale(vec.Norm());
-  // scale normal vector by sine and this vector by cosine
-  y.Scale(sin(alpha));
-  a.Scale(cos(alpha));
-  res = vec.Projection(axis);
-  // add scaled normal vector onto this vector
-  res += y;
-  // add part in axis direction
-  res += a;
-  return res;
-};
-
-/** Calculates the intersection of the two lines that are both on the same plane.
- * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
- * \param *out output stream for debugging
- * \param *Line1a first vector of first line
- * \param *Line1b second vector of first line
- * \param *Line2a first vector of second line
- * \param *Line2b second vector of second line
- * \return true - \a this will contain the intersection on return, false - lines are parallel
- */
-Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b)
-{
-  Info FunctionInfo(__func__);
-
-  Vector res;
-
-  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
-
-  M->SetAll(1.);
-  for (int i=0;i<3;i++) {
-    M->Set(0, i, Line1a[i]);
-    M->Set(1, i, Line1b[i]);
-    M->Set(2, i, Line2a[i]);
-    M->Set(3, i, Line2b[i]);
-  }
-
-  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
-  //for (int i=0;i<4;i++) {
-  //  for (int j=0;j<4;j++)
-  //    cout << "\t" << M->Get(i,j);
-  //  cout << endl;
-  //}
-  if (fabs(M->Determinant()) > MYEPSILON) {
-    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
-
-
-  // constuct a,b,c
-  Vector a = Line1b - Line1a;
-  Vector b = Line2b - Line2a;
-  Vector c = Line2a - Line1a;
-  Vector d = Line2b - Line1b;
-  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
-  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
-   res.Zero();
-   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
-   throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // check for parallelity
-  Vector parallel;
-  double factor = 0.;
-  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
-    parallel = Line1a - Line2a;
-    factor = parallel.ScalarProduct(a)/a.Norm();
-    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-      res = Line2a;
-      Log() << Verbose(1) << "Lines conincide." << endl;
-      return res;
-    } else {
-      parallel = Line1a - Line2b;
-      factor = parallel.ScalarProduct(a)/a.Norm();
-      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-        res = Line2b;
-        Log() << Verbose(1) << "Lines conincide." << endl;
-        return res;
-      }
-    }
-    Log() << Verbose(1) << "Lines are parallel." << endl;
-    res.Zero();
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // obtain s
-  double s;
-  Vector temp1, temp2;
-  temp1 = c;
-  temp1.VectorProduct(b);
-  temp2 = a;
-  temp2.VectorProduct(b);
-  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
-  if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
-  else
-    s = 0.;
-  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
-
-  // construct intersection
-  res = a;
-  res.Scale(s);
-  res += Line1a;
-  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
-
-  return res;
-};

src/vector_ops.hpp

@@ -10 +10 @@
 
 bool LSQdistance(Vector &res,const Vector **vectors, int num);
-Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha);
-Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b);
 
 #endif /* VECTOR_OPS_HPP_ */