Kernel/Math/Vector.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef MECURYDPM_VECTOR_H
#define MECURYDPM_VECTOR_H
 
#include <cmath>
#include <sstream>
#include <iostream>
#include <cstdlib>
#include <array>
 
#include "GeneralDefine.h"
#include "Logger.h"
 
template<unsigned int N>
class SmallVector;
 
class Vec3D
{
public:
    
    /*
     * \todo: Make these private.
     * \todo what is the idea of this constructor?
     * These should be private so we can implement things like a cvec etc.
     * Use getters / setters.
     */
//    Vec3D(int i);
    
    // private:
    Mdouble X, Y, Z;
    
    Vec3D()
    { setZero(); }
    
    Vec3D(const SmallVector<3>& vector);
 
    Vec3D(const std::array<Mdouble, 3>& vector)
    {
        X = vector[0];
        Y = vector[1];
        Z = vector[2];
    };
 
    Vec3D(const Mdouble x, const Mdouble y, const Mdouble z)
    {
        X = x;
        Y = y;
        Z = z;
    }
 
    void setZero();
    
    void setNaN();
    
    bool isZero() const
    { return X == 0.0 && Y == 0.0 && Z == 0.0; }
    
    bool isNaN() const;
    
    Vec3D operator+(const Vec3D& a) const
    {
        return Vec3D(X + a.X, Y + a.Y, Z + a.Z);
    }
    
    inline Vec3D operator-(const Vec3D a) const
    {
        return Vec3D(X - a.X, Y - a.Y, Z - a.Z);
    };
 
    inline bool operator==(const Vec3D& a) const
    {
        return X == a.X and Y == a.Y and Z == a.Z;
    };
 
    Vec3D multiplyElementwise(const Vec3D& a) const {
        return Vec3D(X*a.X, Y*a.Y, Z*a.Z);
    }
 
    Vec3D divideElementwise(const Vec3D& a) const {
        return Vec3D(X/a.X, Y/a.Y, Z/a.Z);
    }
 
    Vec3D signedSquare() const {
        return Vec3D(fabs(X)*X, fabs(Y)*Y, fabs(Z)*Z);
    }
 
    Vec3D operator*(const Mdouble a) const
    {
        return Vec3D(X * a, Y * a, Z * a);
    }
    
    Vec3D operator/(Mdouble a) const {
        return Vec3D(X / a, Y / a, Z / a);
    }
    
    Vec3D& operator+=(const Vec3D& a)
    {
        X += a.X;
        Y += a.Y;
        Z += a.Z;
        return *this;
    }
 
    bool operator>=(const Vec3D& a) const {
        return X>=a.X && Y>=a.Y && Z>=a.Z;
    }
 
    bool operator<(const Vec3D& a) const {
        return X<a.X && Y<a.Y && Z<a.Z;
    }
 
    Vec3D& operator-=(const Vec3D& a)
    {
        X -= a.X;
        Y -= a.Y;
        Z -= a.Z;
        return *this;
    }
 
    Vec3D& operator*=(Mdouble a) {
        X *= a;
        Y *= a;
        Z *= a;
        return *this;
    }
 
    Vec3D& operator/=(const Mdouble a)
    {
        X /= a;
        Y /= a;
        Z /= a;
        return *this;
    }
 
 
    static Mdouble dot(const Vec3D& a, const Vec3D& b);
    
    inline Mdouble dot(const Vec3D& b) const { return Vec3D::dot(*this, b); }
 
    static Vec3D max(const Vec3D& a, const Vec3D& b);
    
    static Vec3D min(const Vec3D& a, const Vec3D& b);
 
    static double max(const Vec3D& a) {return std::max(std::max(a.X,a.Y),a.Z);}
 
    static double min(const Vec3D& a) {return std::min(std::min(a.X,a.Y),a.Z);}
 
    static Vec3D square(const Vec3D& a);
    
    void normalise();
 
    void setLength(Mdouble length);
    
    static Vec3D sqrt(const Vec3D& a);
    
    static Vec3D cross(const Vec3D& a, const Vec3D& b);
    
    inline Vec3D cross(const Vec3D& b) const { return Vec3D::cross(*this, b); }
    
    static Mdouble getDistance(const Vec3D& a, const Vec3D& b);
    
    static Mdouble getDistanceSquared(const Vec3D& a, const Vec3D& b) {
        const double X = a.X-b.X;
        const double Y = a.Y-b.Y;
        const double Z = a.Z-b.Z;
        return (X * X + Y * Y + Z * Z);
        //return getLengthSquared(a - b);
    }
 
    
    static Mdouble getLength(const Vec3D& a);
    
    static Mdouble getLengthSquared(const Vec3D& a)
    {
        return (a.X * a.X + a.Y * a.Y + a.Z * a.Z);
    }
    
    Mdouble getLength() const;
    
    Mdouble getLengthSquared() const;
    
    Mdouble getComponent(int index) const;
    
    void setComponent(int index, double val);
 
 
    Mdouble& operator[] (int index)
    {
        switch (index)
        {
            case 0:
                return X;
            case 1:
                return Y;
            case 2:
                return Z;
            default:
                throw std::out_of_range("Index out of bounds");
        }
    }
 
    const Mdouble& operator[] (int index) const
    {
        switch (index)
        {
            case 0:
                return X;
            case 1:
                return Y;
            case 2:
                return Z;
            default:
                throw std::out_of_range("Index out of bounds");
        }
    }
 
    inline Mdouble& x()
    { return X; }
    
    inline Mdouble x() const
    { return X; }
    
    inline Mdouble& y()
    { return Y; }
    
    inline Mdouble y() const
    { return Y; }
    
    inline Mdouble& z()
    { return Z; }
    
    inline Mdouble z() const
    { return Z; }
    
    inline void setX(Mdouble x)
    { X = x; }
    
    inline void setY(Mdouble y)
    { Y = y; }
    
    inline void setZ(Mdouble z)
    { Z = z; }
    
    inline Mdouble getX() const
    { return X; }
    
    inline Mdouble getY() const
    { return Y; }
    
    inline Mdouble getZ() const
    { return Z; }
    
    inline void set(Mdouble x, Mdouble y, Mdouble z)
    {
        X = x;
        Y = y;
        Z = z;
    }
    
    Mdouble getRadialCoordinateSquared() const;
    
    Mdouble getRadialCoordinate() const;
 
    Vec3D getSphericalCoordinates() const;
 
    Vec3D getCylindricalCoordinates() const;
    
    Vec3D getFromCylindricalCoordinates() const;
    
    Vec3D getCylindricalTensorField(const Vec3D& position) const;
    
    bool isEqualTo(const Vec3D& other, double tol) const;
    
    static Vec3D getUnitVector(const Vec3D& a);
    
    friend std::ostream& operator<<(std::ostream& os, const Vec3D& a);
    
    friend std::istream& operator>>(std::istream& is, Vec3D& a);
    
    inline friend Vec3D operator-(const Vec3D& a) {
        return Vec3D(-a.X, -a.Y, -a.Z);
    }
    
    friend Vec3D operator*(Mdouble a, const Vec3D& b) {
        return Vec3D(b.X * a, b.Y * a, b.Z * a);
    }
 
};
 
#endif