Vector.h

Go to the documentation of this file.

/*
 This file is part of the Ristra wonton project.
 Please see the license file at the root of this repository, or at:
 https://github.com/laristra/wonton/blob/master/LICENSE
*/

// Copyright 2016 Los Alamos National Laboratory                             //
//                                                                           //
// Original Author: Paul Henning, from the SCF library.                      //
//                                                                           //
// Modified for GK by: Brian Jean                                            //
//                     TSM, ASCI Problem Setup Team                          //
//                     Applied Physics Division                              //
//                     Los Alamos National Laboratory                        //
//                     505.665.6374                                          //
//                     baj@lanl.gov                                          //
//                                                                           //
// Modified for Wonton by: Rao Garimella, rao@lanl.gov                       //

#ifndef WONTON_VECTOR_H_
#define WONTON_VECTOR_H_

#include <cmath>
#include <cassert>
#include <iostream>
#include <vector>
#include <limits>
#include "wonton/support/wonton.h"

namespace Wonton {

template <int D> class Vector {
 private:
  double m_comp[D] {};

 public:

  WONTON_INLINE
  Vector() {
    for (int i = 0; i < D; i++)
      m_comp[i] = 0.0;
  }

  WONTON_INLINE
  Vector(const double& xm_comp) {
    for (int i = 0; i < D; i++)
      m_comp[i] = xm_comp;
  }

  WONTON_INLINE
  Vector(const double& xm_comp, const double& ym_comp) {
    assert(D == 2);
    m_comp[0] = xm_comp;
    m_comp[1] = ym_comp;
  }

  WONTON_INLINE
  Vector(const double& xm_comp, const double& ym_comp, const double& zm_comp) {
    assert(D == 3);
    m_comp[0] = xm_comp;
    m_comp[1] = ym_comp;
    m_comp[2] = zm_comp;
  }

  inline Vector(std::vector<double> const& invec) {
    assert(D == invec.size());
    for (int i = 0; i < D; i++) m_comp[i] = invec[i];
  }

  WONTON_INLINE
  const double& operator[](const int& i) const {
    return m_comp[i];
  }

  WONTON_INLINE
  double& operator[](const int& i) {
    return m_comp[i];
  }

  WONTON_INLINE
  Vector operator-() const {
    Vector v;
    for (int i = 0; i < D; i++) v.m_comp[i] = -m_comp[i];
    return v;
  }

  WONTON_INLINE
  Vector& operator+=(const Vector<D>& rhs) {
    for (int i = 0; i < D; i++) m_comp[i] += rhs.m_comp[i];
    return *this;
  }

  WONTON_INLINE
  Vector& operator-=(const Vector<D>& rhs) {
    for (int i = 0; i < D; i++) m_comp[i] -= rhs.m_comp[i];
    return *this;
  }

  WONTON_INLINE
  Vector& operator*=(const double& s) {
    for (int i = 0; i < D; i++) m_comp[i] *= s;
    return *this;
  }

  WONTON_INLINE
  Vector& operator/=(const double& s) {
    for (int i = 0; i < D; i++) m_comp[i] /= s;
    return *this;
  }

  WONTON_INLINE
  double norm(bool doSqrt = true) const {
    double result = 0.0;
    for (int i = 0; i < D; i++) result += (m_comp[i] * m_comp[i]);
    return doSqrt ? sqrt(result) : result;
  }

  WONTON_INLINE
  double one_norm() const {
    double result = 0.0;
    for (int i = 0; i < D; i++) result += std::fabs(m_comp[i]);
    return result;
  }

  WONTON_INLINE
  double max_norm() const {
    double result = std::fabs(m_comp[0]);
    for (int i = 0; i < D - 1; i++) {
      double abs_val = std::fabs(m_comp[i + 1]);
      if (result < abs_val)
        result = abs_val;
    }
    return result;
  }

  WONTON_INLINE
  void normalize() {
    double s = norm();
    *this /= s;
  }

  WONTON_INLINE
  void zero() {
    for (int i = 0; i < D; i++) m_comp[i] = 0;
  }

  WONTON_INLINE
  bool is_zero(double dst_tol) const {
    return (norm() < dst_tol);
  }

  WONTON_INLINE
  void fill(double val) {
    for (int i = 0; i < D; i++) m_comp[i] = val;
  }

  WONTON_INLINE
  void axis(int nonZero) {
    zero();
    m_comp[nonZero] = 1;
  }

  std::istream& readFromStream(std::istream& is) {
    for (int i = 0; i < D; i++)
      is >> m_comp[i];
    return is;
  }

  std::ostream& writeToStream(std::ostream& os) const {
    for (int i = 0; i < D; i++) {
      if (i > 0) os << ' ';
      os << m_comp[i];
    }
    return os;
  }
};

typedef Vector<3> Vector3;

typedef Vector<2> Vector2;

template<int D>
WONTON_INLINE
double dot(const Vector<D>& a, const Vector<D>& b) {
  double r = 0.0;
  for (int i = 0; i < D; i++) r += a[i] * b[i];
  return r;
}

template<int D>
WONTON_INLINE
Vector<D> operator+(const Vector<D>& a, const Vector<D>& b) {
  return Vector<D>(a) += b;
}

template<int D>
WONTON_INLINE
Vector<D> operator-(const Vector<D>& a, const Vector<D>& b) {
  return Vector<D>(a) -= b;
}

template<int D>
WONTON_INLINE
Vector<D> operator*(const Vector<D>& a, const double& s) {
  return Vector<D>(a) *= s;
}

template<int D>
WONTON_INLINE
Vector<D> operator*(const double& s, const Vector<D>& a) {
  return Vector<D>(a) *= s;
}

template<int D>
WONTON_INLINE
Vector<D> operator/(const Vector<D>& a, const double& s) {
  return Vector<D>(a) /= s;
}

template<int D>
inline
std::ostream& operator<<(std::ostream& os, const Vector<D>& v) {
  return v.writeToStream(os);
}

template<int D> inline std::istream&
operator>>(std::istream& is, Vector<D>& v) {
  return v.readFromStream(is);
}

WONTON_INLINE
double cross(const Vector<2>& a, const Vector<2>& b) {
  return (a[0] * b[1] - a[1] * b[0]);
}

WONTON_INLINE
Vector<3> cross(const Vector<3>& a, const Vector<3>& 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<int D>
WONTON_INLINE
double MaxComponent(const Vector<D>& v, int& icomp) {
  double max = v[0];
  icomp = 0;
  for (int i = 1; i < D; i++)
    if (max < v[i]) {
      max = v[i];
      icomp = i;
    }
  return max;
}

}  // namespace Wonton


#endif  // WONTON_VECTOR_H_