vector3.h

Go to the documentation of this file.

#ifndef VECTOR3_H
#define VECTOR3_H
 
#include <cmath>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <array>
 
#ifdef _MCD_CHECK
#endif
 
namespace ModuleBase
{
    // Small epsilon value for numerical comparisons
    constexpr double epsilon = 1e-10;
 
template <class T> class Vector3
{
  public:
    T x;
    T y;
    T z;
 
    Vector3(const T &x1 = 0, const T &y1 = 0, const T &z1 = 0)
        : x(x1), y(y1), z(z1)
    {}
    Vector3(const Vector3<T> &v)
        : x(v.x), y(v.y), z(v.z)
    {} 
    explicit Vector3(const std::array<T,3> &v)
        : x(v[0]), y(v[1]), z(v[2])
    {}
 
    template <typename U>
    explicit Vector3(const Vector3<U>& other)
        : x(static_cast<T>(other.x)), y(static_cast<T>(other.y)), z(static_cast<T>(other.z))
    {}
 
    Vector3(Vector3<T> &&v) noexcept
        : x(v.x), y(v.y), z(v.z)
    {}
 
    void set(const T &x1, const T &y1, const T &z1)
    {
        x = x1;
        y = y1;
        z = z1;
    }
 
    Vector3<T> &operator=(const Vector3<T> &u)
    {
        x = u.x;
        y = u.y;
        z = u.z;
        return *this;
    }
 
    Vector3<T> &operator=(const T &u)
    {
        x = u;
        y = u;
        z = u;
        return *this;
    }
 
    Vector3<T> &operator=(Vector3<T> &&u) noexcept
    {
        x = u.x;
        y = u.y;
        z = u.z;
        return *this;
    }
 
    Vector3<T> &operator+=(const Vector3<T> &u)
    {
        x += u.x;
        y += u.y;
        z += u.z;
        return *this;
    }
 
    Vector3<T> &operator-=(const Vector3<T> &u)
    {
        x -= u.x;
        y -= u.y;
        z -= u.z;
        return *this;
    }
 
    Vector3<T> &operator*=(const T &s)
    {
        x *= s;
        y *= s;
        z *= s;
        return *this;
    }
 
    Vector3<T> &operator/=(const T &s)
    {
        assert(s != 0); // Avoid division by zero
        x /= s;
        y /= s;
        z /= s;
        return *this;
    }
 
    Vector3<T> operator-() const
    {
        return Vector3<T>(-x, -y, -z);
    }
 
    T operator[](int index) const
    {
        assert(index >= 0 && index < 3);
        //return (&x)[index]; // this is undefind behavior and breaks with icpx
        T const* ptr[3] = {&x, &y, &z};
        return *ptr[index];
    }
 
    T &operator[](int index)
    {
        assert(index >= 0 && index < 3);
        //return (&x)[index]; // this is undefind behavior and breaks with icpx
        T* ptr[3] = {&x, &y, &z};
        return *ptr[index];
    }
 
    inline T norm2(void) const
    {
        return x * x + y * y + z * z;
    }
 
    inline T norm(void) const
    {
        return sqrt(norm2());
    }
 
    Vector3<T> &normalize(void)
    {
        const T m = norm();
        if (m > epsilon) // Avoid division by zero
        {
            x /= m;
            y /= m;
            z /= m;
        }
        return *this;
    } // Peize Lin update return 2019-09-08
 
    Vector3<T> &reverse(void)
    {
        x = -x;
        y = -y;
        z = -z;
        return *this;
    } // Peize Lin update return 2019-09-08
 
    void print(const int precision = 5) const;
};
 
template <class T> inline Vector3<T> operator+(const Vector3<T> &u, const Vector3<T> &v)
{
    return Vector3<T>(u.x + v.x, u.y + v.y, u.z + v.z);
}
 
template <class T> inline Vector3<T> operator-(const Vector3<T> &u, const Vector3<T> &v)
{
    return Vector3<T>(u.x - v.x, u.y - v.y, u.z - v.z);
}
 
template <class T> inline T operator*(const Vector3<T> &u, const Vector3<T> &v)
{
    return (u.x * v.x + u.y * v.y + u.z * v.z);
}
 
template <class T> inline Vector3<T> operator*(const T &s, const Vector3<T> &u)
{
    return Vector3<T>(u.x * s, u.y * s, u.z * s);
}
 
template <class T> inline Vector3<T> operator*(const Vector3<T> &u, const T &s)
{
    return Vector3<T>(u.x * s, u.y * s, u.z * s);
} // mohan add 2009-5-10
 
template <class T> inline Vector3<T> operator/(const Vector3<T> &u, const T &s)
{
    return Vector3<T>(u.x / s, u.y / s, u.z / s);
}
 
template <class T> inline Vector3<T> operator/(const T &s, const Vector3<T> &u)
{
    return Vector3<T>(s/u.x, s/u.y, s/u.z);
}
 
template <class T> inline T dot(const Vector3<T> &u, const Vector3<T> &v)
{
    return (u.x * v.x + u.y * v.y + u.z * v.z);
}
 
template <class T> inline Vector3<T> operator^(const Vector3<T> &u, const Vector3<T> &v)
{
    return Vector3<T>(u.y * v.z - u.z * v.y, -u.x * v.z + u.z * v.x, u.x * v.y - u.y * v.x);
}
 
template <class T> inline Vector3<T> cross(const Vector3<T> &u, const Vector3<T> &v)
{
    return Vector3<T>(u.y * v.z - u.z * v.y, -u.x * v.z + u.z * v.x, u.x * v.y - u.y * v.x);
}
// s = u.(v x w)
// template <class T> T TripleScalarProduct(Vector3<T> u, Vector3<T> v, Vector3<T> w)
//{
//  return T((u.x * (v.y * w.z - v.z * w.y)) +
//           (u.y * (-v.x * w.z + v.z * w.x)) +
//           (u.z * (v.x * w.y - v.y * w.x)));
// }
 
// Overload the < operator for sorting
template <class T> bool operator<(const Vector3<T> &u, const Vector3<T> &v)
{
    if (u.x < v.x)
    {
        return true;
    }
    if (u.x > v.x)
    {
        return false;
    }
    if (u.y < v.y)
    {
        return true;
    }
    if (u.y > v.y)
    {
        return false;
    }
    if (u.z < v.z)
    {
        return true;
    }
    return false;
}
 
// whether m1 != m2
template <class T> inline bool operator!=(const Vector3<T> &u, const Vector3<T> &v)
{
    return !(u == v);
}
// whether u == v
template <class T> inline bool operator==(const Vector3<T> &u, const Vector3<T> &v)
{
    if (u.x == v.x && u.y == v.y && u.z == v.z)
    {
        return true;
    }
    return false;
}
 
template <class T> void Vector3<T>::print(const int precision) const
{
    // Ensure precision is non-negative
    int valid_precision = precision > 0 ? precision : 5;
    std::cout.precision(valid_precision);
    std::cout << "(" << std::setw(10) << x << "," << std::setw(10) << y << "," << std::setw(10) << z << ")"
              << std::endl;
    return;
}
 
template <class T> static std::ostream &operator<<(std::ostream &os, const Vector3<T> &u)
{
    os << "(" << std::setw(10) << u.x << "," << std::setw(10) << u.y << "," << std::setw(10) << u.z << ")";
    return os;
}
 
} // namespace ModuleBase
 
#endif