apfVector.h

Go to the documentation of this file.

/*
 * Copyright 2011 Scientific Computation Research Center
 *
 * This work is open source software, licensed under the terms of the
 * BSD license as described in the LICENSE file in the top-level directory.
 */

#ifndef APF_VECTOR_H
#define APF_VECTOR_H

#include <cmath>
#include "apfArray.h"
#include <ostream>

namespace apf {

extern double const pi;

template <std::size_t N>
class Vector;

template <std::size_t N>
class Vector : public Array<double,N>
{
  public:
    Vector() {}
    Vector(double const* v)
    {
      for (std::size_t i=0; i < N; ++i)
        (*this)[i] = v[i];
    }
    Vector<N> operator+(Vector<N> const& b) const
    {
      Vector<N> c;
      for (std::size_t i=0; i < N; ++i)
        c[i] = (*this)[i] + b[i];
      return c;
    }
    Vector<N>& operator+=(Vector<N> const& b)
    {
      for (std::size_t i=0; i < N; ++i)
        (*this)[i] += b[i];
      return (*this);
    }
    Vector<N> operator-(Vector<N> const& b) const
    {
      Vector<N> c;
      for (std::size_t i=0; i < N; ++i)
        c[i] = (*this)[i] - b[i];
      return c;
    }
    Vector<N> operator*(double s) const
    {
      Vector<N> c;
      for (std::size_t i=0; i < N; ++i)
        c[i] = (*this)[i] * s;
      return c;
    }
    Vector<N> operator/(double s) const
    {
      Vector<N> c;
      for (std::size_t i=0; i < N; ++i)
        c[i] = (*this)[i] / s;
      return c;
    }
    double operator*(Vector<N> const& b) const
    {
      double r=0;
      for (std::size_t i=0; i < N; ++i)
        r += (*this)[i] * b[i];
      return r;
    }
    double getLength() const {return sqrt((*this)*(*this));}
    Vector<N> normalize() const {return (*this) / getLength();}
    void zero()
    {
      for (std::size_t i=0; i < N; ++i)
        (*this)[i] = 0.0;
    }
};

inline Vector<3> cross(Vector<3> const& a, Vector<3> const& b)
{
  Vector<3> r;
  r[0] = a[1]*b[2] - a[2]*b[1];
  r[1] = a[2]*b[0] - a[0]*b[2];
  r[2] = a[0]*b[1] - a[1]*b[0];
  return r;
}

template<std::size_t N>
Vector<N> project(Vector<N> const& a, Vector<N> const& b)
{
  return b*((a*b)/(b*b));
}

template<std::size_t N>
Vector<N> reject(Vector<N> const& a, Vector<N> const& b)
{
  return a - project(a, b);
}

class Vector3 : public Vector<3>
{
  public:
    Vector3() {}
    Vector3(Vector<3> const& other):
      Vector<3>(other)
    {}
    Vector3(double a, double b, double c)
    {
      (*this)[0] = a;
      (*this)[1] = b;
      (*this)[2] = c;
    }
    Vector3(double const* abc)
    {
      (*this)[0] = abc[0];
      (*this)[1] = abc[1];
      (*this)[2] = abc[2];
    }
    void toArray(double* abc) const
    {
      abc[0] = (*this)[0];
      abc[1] = (*this)[1];
      abc[2] = (*this)[2];
    }
    void fromArray(const double* abc)
    {
      (*this)[0] = abc[0];
      (*this)[1] = abc[1];
      (*this)[2] = abc[2];
    }
    double x() const {return (*this)[0];}
    double y() const {return (*this)[1];}
    double z() const {return (*this)[2];}
    double& x() {return (*this)[0];}
    double& y() {return (*this)[1];}
    double& z() {return (*this)[2];}
};

}

std::ostream& operator<<(std::ostream& s, apf::Vector3 const& v);

#endif