Changes in src/vector.cpp [4b94bb:112b09]

File:

• src/vector.cpp (modified) (19 diffs)

src/vector.cpp

@@ -12 +12 @@
 #include "Helpers/Assert.hpp"
 #include "Helpers/fast_functions.hpp"
-#include "Exceptions/MathException.hpp"
 
 #include <iostream>
-#include <gsl/gsl_blas.h>
-
 
 using namespace std;
@@ -27 +24 @@
 Vector::Vector()
 {
-  content = gsl_vector_calloc (NDIM);
+  x[0] = x[1] = x[2] = 0.;
 };
 
@@ -36 +33 @@
 Vector::Vector(const Vector& src)
 {
-  content = gsl_vector_alloc(NDIM);
-  gsl_vector_memcpy(content, src.content);
+  x[0] = src[0];
+  x[1] = src[1];
+  x[2] = src[2];
 }
 
@@ -44 +42 @@
 Vector::Vector(const double x1, const double x2, const double x3)
 {
-  content = gsl_vector_alloc(NDIM);
-  gsl_vector_set(content,0,x1);
-  gsl_vector_set(content,1,x2);
-  gsl_vector_set(content,2,x3);
-};
-
-Vector::Vector(gsl_vector *_content) :
-  content(_content)
-{}
+  x[0] = x1;
+  x[1] = x2;
+  x[2] = x3;
+};
 
 /**
@@ -60 +53 @@
   // check for self assignment
   if(&src!=this){
-    gsl_vector_memcpy(content, src.content);
+    x[0] = src[0];
+    x[1] = src[1];
+    x[2] = src[2];
   }
   return *this;
@@ -67 +62 @@
 /** Desctructor of class vector.
  */
-Vector::~Vector() {
-  gsl_vector_free(content);
-};
+Vector::~Vector() {};
 
 /** Calculates square of distance between this and another vector.
@@ -79 +72 @@
   double res = 0.;
   for (int i=NDIM;i--;)
-    res += (at(i)-y[i])*(at(i)-y[i]);
+    res += (x[i]-y[i])*(x[i]-y[i]);
   return (res);
 };
@@ -97 +90 @@
 }
 
+/** Calculates distance between this and another vector in a periodic cell.
+ * \param *y array to second vector
+ * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
+ * \return \f$| x - y |\f$
+ */
+double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
+{
+  double res = distance(y), tmp, matrix[NDIM*NDIM];
+    Vector Shiftedy, TranslationVector;
+    int N[NDIM];
+    matrix[0] = cell_size[0];
+    matrix[1] = cell_size[1];
+    matrix[2] = cell_size[3];
+    matrix[3] = cell_size[1];
+    matrix[4] = cell_size[2];
+    matrix[5] = cell_size[4];
+    matrix[6] = cell_size[3];
+    matrix[7] = cell_size[4];
+    matrix[8] = cell_size[5];
+    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
+    for (N[0]=-1;N[0]<=1;N[0]++)
+      for (N[1]=-1;N[1]<=1;N[1]++)
+        for (N[2]=-1;N[2]<=1;N[2]++) {
+          // create the translation vector
+          TranslationVector.Zero();
+          for (int i=NDIM;i--;)
+            TranslationVector[i] = (double)N[i];
+          TranslationVector.MatrixMultiplication(matrix);
+          // add onto the original vector to compare with
+          Shiftedy = y + TranslationVector;
+          // get distance and compare with minimum so far
+          tmp = distance(Shiftedy);
+          if (tmp < res) res = tmp;
+        }
+    return (res);
+};
+
+/** Calculates distance between this and another vector in a periodic cell.
+ * \param *y array to second vector
+ * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
+ * \return \f$| x - y |^2\f$
+ */
+double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
+{
+  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
+    Vector Shiftedy, TranslationVector;
+    int N[NDIM];
+    matrix[0] = cell_size[0];
+    matrix[1] = cell_size[1];
+    matrix[2] = cell_size[3];
+    matrix[3] = cell_size[1];
+    matrix[4] = cell_size[2];
+    matrix[5] = cell_size[4];
+    matrix[6] = cell_size[3];
+    matrix[7] = cell_size[4];
+    matrix[8] = cell_size[5];
+    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
+    for (N[0]=-1;N[0]<=1;N[0]++)
+      for (N[1]=-1;N[1]<=1;N[1]++)
+        for (N[2]=-1;N[2]<=1;N[2]++) {
+          // create the translation vector
+          TranslationVector.Zero();
+          for (int i=NDIM;i--;)
+            TranslationVector[i] = (double)N[i];
+          TranslationVector.MatrixMultiplication(matrix);
+          // add onto the original vector to compare with
+          Shiftedy = y + TranslationVector;
+          // get distance and compare with minimum so far
+          tmp = DistanceSquared(Shiftedy);
+          if (tmp < res) res = tmp;
+        }
+    return (res);
+};
+
+/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
+ * \param *out ofstream for debugging messages
+ * Tries to translate a vector into each adjacent neighbouring cell.
+ */
+void Vector::KeepPeriodic(const double * const matrix)
+{
+  //  int N[NDIM];
+  //  bool flag = false;
+    //vector Shifted, TranslationVector;
+  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
+  //  Log() << Verbose(2) << "Vector is: ";
+  //  Output(out);
+  //  Log() << Verbose(0) << endl;
+    InverseMatrixMultiplication(matrix);
+    for(int i=NDIM;i--;) { // correct periodically
+      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
+        at(i) += ceil(at(i));
+      } else {
+        at(i) -= floor(at(i));
+      }
+    }
+    MatrixMultiplication(matrix);
+  //  Log() << Verbose(2) << "New corrected vector is: ";
+  //  Output(out);
+  //  Log() << Verbose(0) << endl;
+  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
+};
+
 /** Calculates scalar product between this and another vector.
  * \param *y array to second vector
@@ -104 +199 @@
 {
   double res = 0.;
-  gsl_blas_ddot(content, y.content, &res);
+  for (int i=NDIM;i--;)
+    res += x[i]*y[i];
   return (res);
 };
@@ -118 +214 @@
 {
   Vector tmp;
-  for(int i=NDIM;i--;)
-    tmp[i] = at((i+1)%NDIM)*y[(i+2)%NDIM] - at((i+2)%NDIM)*y[(i+1)%NDIM];
+  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;
 };
@@ -212 +309 @@
 bool Vector::IsZero() const
 {
-  return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < MYEPSILON);
+  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
 };
 
@@ -241 +338 @@
   bool status = true;
   for (int i=0;i<NDIM;i++) {
-    if (fabs(at(i) - a[i]) > MYEPSILON)
+    if (fabs(x[i] - a[i]) > MYEPSILON)
       status = false;
   }
@@ -269 +366 @@
 double& Vector::operator[](size_t i){
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return *gsl_vector_ptr (content, i);
+  return x[i];
 }
 
 const double& Vector::operator[](size_t i) const{
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return *gsl_vector_ptr (content, i);
+  return x[i];
 }
 
@@ -285 +382 @@
 }
 
-gsl_vector* Vector::get(){
-  return content;
+double* Vector::get(){
+  return x;
 }
 
@@ -401 +498 @@
 {
   for (int i=NDIM;i--;)
-    at(i) *= factor[i];
-};
-
-void Vector::ScaleAll(const Vector &factor){
-  gsl_vector_mul(content, factor.content);
-}
+    x[i] *= factor[i];
+};
+
 
 
 void Vector::Scale(const double factor)
 {
-  gsl_vector_scale(content,factor);
+  for (int i=NDIM;i--;)
+    x[i] *= factor;
+};
+
+/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
+ * \param *M matrix of box
+ * \param *Minv inverse matrix
+ */
+void Vector::WrapPeriodically(const double * const M, const double * const Minv)
+{
+  MatrixMultiplication(Minv);
+  // truncate to [0,1] for each axis
+  for (int i=0;i<NDIM;i++) {
+    //x[i] += 0.5;  // set to center of box
+    while (x[i] >= 1.)
+      x[i] -= 1.;
+    while (x[i] < 0.)
+      x[i] += 1.;
+  }
+  MatrixMultiplication(M);
 };
 
@@ -431 +544 @@
 }
 
+/** Do a matrix multiplication.
+ * \param *matrix NDIM_NDIM array
+ */
+void Vector::MatrixMultiplication(const double * const M)
+{
+  // do the matrix multiplication
+  at(0) = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
+  at(1) = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
+  at(2) = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
+};
+
+/** Do a matrix multiplication with the \a *A' inverse.
+ * \param *matrix NDIM_NDIM array
+ */
+bool Vector::InverseMatrixMultiplication(const double * const A)
+{
+  double B[NDIM*NDIM];
+  double detA = RDET3(A);
+  double detAReci;
+
+  // calculate the inverse B
+  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
+    detAReci = 1./detA;
+    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
+    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
+    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
+    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
+    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
+    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
+    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
+    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
+    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+
+    // do the matrix multiplication
+    at(0) = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
+    at(1) = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
+    at(2) = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
+
+    return true;
+  } else {
+    return false;
+  }
+};
+
+
 /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
  * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
@@ -458 +616 @@
   SubtractVector(x1);
   for (int i=NDIM;i--;)
-    result = result || (fabs(at(i)) > MYEPSILON);
+    result = result || (fabs(x[i]) > MYEPSILON);
 
   return result;
@@ -519 +677 @@
 void Vector::AddVector(const Vector &y)
 {
-  gsl_vector_add(content, y.content);
+  for(int i=NDIM;i--;)
+    x[i] += y[i];
 }
 
@@ -527 +686 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  gsl_vector_sub(content, y.content);
+  for(int i=NDIM;i--;)
+    x[i] -= y[i];
+}
+
+/**
+ * Checks whether this vector is within the parallelepiped defined by the given three vectors and
+ * their offset.
+ *
+ * @param offest for the origin of the parallelepiped
+ * @param three vectors forming the matrix that defines the shape of the parallelpiped
+ */
+bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
+{
+  Vector a = (*this)-offset;
+  a.InverseMatrixMultiplication(parallelepiped);
+  bool isInside = true;
+
+  for (int i=NDIM;i--;)
+    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
+
+  return isInside;
 }