Vector.h

#ifndef UTIL_VECTOR_H
#define UTIL_VECTOR_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <util/global.h>
#include <util/space/Dimension.h>

#ifdef UTIL_MPI
#include <util/mpi/MpiTraits.h>
#endif

#include <iostream>
#include <cmath>

namespace Util
{

   class Vector
   {

   public:



      Vector();

      Vector(const Vector& v);

      explicit Vector(double scalar);

      explicit Vector(const double* v);

      Vector(double x, double y, double z=0.0);


      Vector& zero();

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);



      Vector& operator = (const Vector& v);

      Vector& operator = (const double* v);



      void operator += (const Vector& dv);

      void operator -= (const Vector& dv);

      void operator *= (double s);

      void operator /= (double s);



      const double& operator [] (int i) const;

      double& operator [] (int i);



      double square() const;

      double abs() const;

      double dot(const Vector& v) const;

      double projection(const Vector& p) const;



      Vector& add(const Vector& v1, const Vector& v2);

      Vector& subtract(const Vector& v1, const Vector& v2);

      Vector& multiply(const Vector& v, double s);

      Vector& divide(const Vector& v, double s);

      Vector& cross(const Vector& v1, const Vector& v2);

      Vector& versor(const Vector& v);

      Vector& parallel(const Vector& v, const Vector& p);

      Vector& transverse(const Vector& v, const Vector& p);

      int minId(const Vector& v);

      int maxId(const Vector& v);



      static void initStatic();

      static const Vector Zero;

      #ifdef UTIL_MPI

      static void commitMpiType();
      #endif


   private:

      static const int Width = 24;

      static const int Precision = 16;

      double elem_[Dimension];

   //friends:

      friend bool operator == (const Vector& v1, const Vector& v2);

      friend bool operator == (const Vector& v1, const double* v2);

      friend 
      std::istream& operator >> (std::istream& in, Vector &vector);

      friend 
      std::ostream& operator << (std::ostream& out, const Vector &vector);

   };

   bool operator == (const Vector& v1, const Vector& v2);

   bool operator == (const Vector& v1, const double* v2);

   bool operator == (const double* v1, const Vector& v2);

   bool operator != (const Vector& v1, const Vector& v2);

   bool operator != (const Vector& v1, const double* v2);

   bool operator != (const double* v1, const Vector& v2);

   std::istream& operator >> (std::istream& in, Vector &vector);

   std::ostream& operator << (std::ostream& out, const Vector &vector);

   #ifdef UTIL_MPI

   template <>
   class MpiTraits<Vector>
   {
   public:
      static MPI::Datatype type;   
      static bool hasType;         
   };
   #endif

   // Inline methods

   /*
   * Default constructor
   */
   inline
   Vector::Vector()
   {}

   /*
   * Copy constructor
   */
   inline
   Vector::Vector(const Vector& v)
   {
      elem_[0] = v.elem_[0];
      elem_[1] = v.elem_[1];
      elem_[2] = v.elem_[2];
   }

   /*
   * Constructor, initialize all elements to a scalar value.
   */
   inline
   Vector::Vector(double scalar)
   {
      elem_[0] = scalar;
      elem_[1] = scalar;
      elem_[2] = scalar;
   }

   /*
   * Construct Vector from C double[3] array.
   */
   inline
   Vector::Vector(const double* v)
   {
      elem_[0] = v[0];
      elem_[1] = v[1];
      elem_[2] = v[2];
   }

   /*
   * Construct Vector from its coordinates.
   */
   inline
   Vector::Vector(double x, double y, double z)
   {
      elem_[0] = x;
      elem_[1] = y;
      elem_[2] = z;
   }

   /*
   * Set all elements of a 3D vector to zero.
   */
   inline
   Vector& Vector::zero()
   {
      elem_[0] = 0.0;
      elem_[1] = 0.0;
      elem_[2] = 0.0;
      return *this;
   }

   /*
   * Assignment.
   */
   inline
   Vector& Vector::operator = (const Vector& v)
   {
      elem_[0] = v.elem_[0];
      elem_[1] = v.elem_[1];
      elem_[2] = v.elem_[2];
      return *this;
   }

   /*
   * Assignment from C double[] array.
   */
   inline
   Vector& Vector::operator = (const double* v)
   {
      elem_[0] = v[0];
      elem_[1] = v[1];
      elem_[2] = v[2];
      return *this;
   }

   /*
   * Add vector dv to this vector.
   */
   inline
   void Vector::operator += (const Vector& dv)
   {
      elem_[0] += dv.elem_[0];
      elem_[1] += dv.elem_[1];
      elem_[2] += dv.elem_[2];
   }

   /*
   * Subtract vector dv from this vector.
   */
   inline
   void Vector::operator -= (const Vector& dv)
   {
      elem_[0] -= dv.elem_[0];
      elem_[1] -= dv.elem_[1];
      elem_[2] -= dv.elem_[2];
   }

   /*
   * Multiply this vector by scalar s.
   */
   inline
   void Vector::operator *= (double s)
   {
      elem_[0] *= s;
      elem_[1] *= s;
      elem_[2] *= s;
   }

   /*
   * Divide this vector by scalar s.
   */
   inline
   void Vector::operator /= (double s)
   {
      elem_[0] /= s;
      elem_[1] /= s;
      elem_[2] /= s;
   }

   /*
   * Return one Cartesian element by value.
   */
   inline
   const double& Vector::operator [] (int i) const
   { 
      assert(i >=0); 
      assert(i < Dimension); 
      return elem_[i]; 
   }

   /*
   * Return a reference to one element of the vector.
   */
   inline
   double& Vector::operator [] (int i)
   {
      assert(i >=0); 
      assert(i < Dimension); 
      return elem_[i]; 
   }

   /*
   * Return square magnitude of this vector.
   */
   inline
   double Vector::square() const
   {
      return (elem_[0]*elem_[0] + elem_[1]*elem_[1] + elem_[2]*elem_[2]);
   }

   /*
   * Return absolute magnitude of this vector.
   */
   inline
   double Vector::abs() const
   {  return sqrt(square()); }

   /*
   * Return dot product of this vector and vector v.
   */
   inline
   double Vector::dot(const Vector& v) const
   {
      return elem_[0]*v.elem_[0] + elem_[1]*v.elem_[1] + elem_[2]*v.elem_[2];
   }

   /*
   * Return projection of this vector along vector p.
   */
   inline
   double Vector::projection(const Vector& p) const
   {
      double abs_p = p.abs();
      if (abs_p > 1.0E-8) {
         return dot(p)/abs_p;
      } else {
         return 0.0;
      }
   }

   /*
   * Add vectors v1 and v2.
   *
   * Upon return, *this = v1 + v2.
   */
   inline
   Vector& Vector::add(const Vector& v1, const Vector& v2)
   {
      elem_[0] = v1.elem_[0] + v2.elem_[0];
      elem_[1] = v1.elem_[1] + v2.elem_[1];
      elem_[2] = v1.elem_[2] + v2.elem_[2];
      return *this;
   }

   /*
   * Subtract vector v2 from v1.
   *
   * Upon return, *this = v1 - v2.
   * \return modified invoking vector
   */
   inline
   Vector& Vector::subtract(const Vector& v1, const Vector& v2)
   {
      elem_[0] = v1.elem_[0] - v2.elem_[0];
      elem_[1] = v1.elem_[1] - v2.elem_[1];
      elem_[2] = v1.elem_[2] - v2.elem_[2];
      return *this;
   }

   /*
   * Multiply a vector v by a scalar s.
   *
   * Upon return, *this = v*s.
   */
   inline
   Vector& Vector::multiply(const Vector& v, double s)
   {
      elem_[0] = v.elem_[0]*s;
      elem_[1] = v.elem_[1]*s;
      elem_[2] = v.elem_[2]*s;
      return *this;
   }

   /*
   * Divide vector v by scalar s.
   *
   * Upon return, *this = v/s;
   */
   inline
   Vector& Vector::divide(const Vector& v, double s)
   {
      elem_[0] = v.elem_[0]/s;
      elem_[1] = v.elem_[1]/s;
      elem_[2] = v.elem_[2]/s;
      return *this;
   }

   /*
   * Calculate cross product of vectors v1 and v2.
   *
   * Upon return, *this = v1 x v2.
   */
   inline
   Vector& Vector::cross(const Vector& v1, const Vector& v2)
   {
      elem_[0] = v1.elem_[1]*v2.elem_[2] - v1.elem_[2]*v2.elem_[1];
      elem_[1] = v1.elem_[2]*v2.elem_[0] - v1.elem_[0]*v2.elem_[2];
      elem_[2] = v1.elem_[0]*v2.elem_[1] - v1.elem_[1]*v2.elem_[0];
      return *this;
   }

   /*
   * Calculate unit vector parallel to v: this = v/|v|
   */
   inline
   Vector& Vector::versor(const Vector& v)
   {
      double vabs = v.abs();
      if (vabs > 1.0E-8) {
         vabs = 1.0/vabs;
         elem_[0] = v.elem_[0]*vabs;
         elem_[1] = v.elem_[1]*vabs;
         elem_[2] = v.elem_[2]*vabs;
      } else {
         elem_[0] = 1.0;
         elem_[1] = 0.0;
         elem_[2] = 0.0;
      };
      return *this;
   }

   /*
   * Calculate component of vector v parallel to vector p.
   *
   * Upon return, the invoking vector is equal to the vector projection of
   * vector v along a direction parallel to vector p.
   */
   inline
   Vector& Vector::parallel(const Vector& v, const Vector& p)
   {
      double pp, fac;
      pp = p.square();
      if (pp > 1.0E-8) {
         fac = v.dot(p)/pp;
      } else {
         fac = 0.0;
      }
      elem_[0] = p.elem_[0]*fac;
      elem_[1] = p.elem_[1]*fac;
      elem_[2] = p.elem_[2]*fac;
      return *this;
   }

   /*
   * Calculate component of vector v transverse to vector p.
   *
   * Upon return, the invoking vector is equal to the vector
   * projection of vector v perpendicular to vector p.
   */
   inline
   Vector& Vector::transverse(const Vector& v, const Vector& p)
   {
      double pp, fac;
      pp = p.square();
      if (pp > 1.0E-8) {
         fac = v.dot(p)/pp;
      } else {
         fac = 0.0;
      }
      elem_[0] = v.elem_[0] - p.elem_[0]*fac;
      elem_[1] = v.elem_[1] - p.elem_[1]*fac;
      elem_[2] = v.elem_[2] - p.elem_[2]*fac;
      return *this;
   }

   /*
   * Compute the index corresponding to minimum element in a vector.
   */
   inline
   int minId(const Vector& v)
   {
      int id = 0;
      for (int i = 1; i < DimensionSq; i++) {
         if (v[i] < v[id])
            id = i;
      }
      return id;
   }

   /*
   * Compute the index corresponding to maximum element in a vector.
   */
   inline
   int maxId(const Vector& v)
   {
      int id = 0;
      for (int i = 1; i < DimensionSq; i++) {
         if (v[i] > v[id])
            id = i;
      }
      return id;
   }

   /*
   * Serialize to/from an archive.
   */
   template <class Archive>
   inline void Vector::serialize(Archive& ar, const unsigned int version)
   {
      ar & elem_[0];
      ar & elem_[1];
      ar & elem_[2];
   }

}
#endif