Context Navigation

-              r47191b
+              r45ef76
 #include "vector.hpp"
+#include "log.hpp"
+#include "verbose.hpp"
+#include "gslmatrix.hpp"
+#include "info.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
+#include "Exceptions/SkewException.hpp"
+Line::Line(Vector &_origin, Vector &_direction) :
+  origin(new Vector(_origin)),
+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()
 …
 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;
+}
+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;
+}
+Line makeLineThrough(const Vector &x1, const Vector &x2){
+  if(x1==x2){
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+  return Line(x1,x1-x2);
+}

src/Line.hpp

-              r47191b
+              r45ef76
 #include <memory>
+#include <vector>
 class Vector;
 …
+{
 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;
 private:
 …
 };
+/**
+ * Named constructor to make a line through two points
+ */
+Line makeLineThrough(const Vector &x1, const Vector &x2);
 #endif /* LINE_HPP_ */

src/Plane.cpp

-              r47191b
+              r45ef76
+}
 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);
 …
  * \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 Vector &Origin, const Vector &LineVector) const
+{
   Info FunctionInfo(__func__);

src/Plane.hpp

-              r47191b
+              r45ef76
    * 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 Vector &Origin, const Vector &LineVector) const;
   /****** Methods inherited from Space ***********/

src/unittests/Makefile.am

-              r47191b
+              r45ef76
   InfoUnitTest \
   LinearSystemOfEquationsUnitTest \
+  LineUnittest \
   LinkedCellUnitTest \
   ListOfBondsUnitTest \
 …
   infounittest.cpp \
   linearsystemofequationsunittest.cpp \
+  LineUnittest.cpp \
   LinkedCellUnitTest.cpp \
   listofbondsunittest.cpp \
 …
   infounittest.hpp \
   linearsystemofequationsunittest.hpp \
+  LineUnittest.hpp \
   LinkedCellUnitTest.hpp \
   listofbondsunittest.hpp \
 …
 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/vectorunittest.cpp

-              r47191b
+              r45ef76
   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 );
 };

src/vector.cpp

-              r47191b
+              r45ef76
 };
+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.
  * \param *matrix NDIM_NDIM array
 …
   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

-              r47191b
+              r45ef76
 #include <memory>
+#include <vector>
 #include "defs.hpp"
 …
 /********************************************** declarations *******************************/
+class Vector;
+typedef std::vector<Vector> pointset;
 /** Single vector.
 …
   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

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 45ef76 for src

Legend:

Download in other formats: