includeLinks/Vector3D.hpp

Go to the documentation of this file.

/* ***************************************************************** 
    MESQUITE -- The Mesh Quality Improvement Toolkit

    Copyright 2004 Sandia Corporation and Argonne National
    Laboratory.  Under the terms of Contract DE-AC04-94AL85000 
    with Sandia Corporation, the U.S. Government retains certain 
    rights in this software.

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License 
    (lgpl.txt) along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
    diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov, 
    pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov      
   
  ***************************************************************** */
#ifndef MESQUITE_VECTOR3D_HPP
#define MESQUITE_VECTOR3D_HPP

#include "Mesquite.hpp"

#ifdef MSQ_USE_OLD_IO_HEADERS
   class ostream;
#else
#  include <iosfwd>
#endif

#ifdef MSQ_USE_OLD_C_HEADERS
#  include <assert.h>
#  include <string.h>
#else
#  include <cassert>
#  include <cstring>
#endif

namespace Mesquite
{
   class Matrix3D;
   class MsqError;
   
  class Vector3D
  {
  public:
    // Constructors
    Vector3D();
    Vector3D(const double& x, const double& y, const double& z);
    Vector3D(const double xyz[3]);
    Vector3D(const Vector3D& to_copy);

    // *** virtual destructor *** Do not use for Vector3D, we need to keep
    // the deallocation of those objects very fast
    
    // Functions to get the coordinates
    double x() const;
    double y() const;
    double z() const;
    void get_coordinates(double& x, double& y, double& z) const;
    void get_coordinates(double xyz[3]) const;
    const double& operator[](size_t index) const; // 0-based
    
    // Functions to set the coordinates.
    void x(const double x);
    void y(const double y);
    void z(const double z);
    void set(const double x, const double y, const double z);
    void set(const double xyz[3]);
    void set(const Vector3D& to_copy);
    // Subscripts on non-consts both get and set coords
    double& operator[](size_t index); // 0-based
    Vector3D& operator=(const Vector3D &to_copy);
    Vector3D& operator=(const double &to_copy);
    
    // Functions to modify existing coordinates
    Vector3D operator-() const;  //- unary negation.
    Vector3D& operator*=(const double scalar);
    Vector3D& operator/=(const double scalar);
    Vector3D& operator*=(const Vector3D &rhs); //- cross product
    Vector3D& operator+=(const Vector3D &rhs);
    Vector3D& operator-=(const Vector3D &rhs);
    
    // Binary operators (like a+b).
    friend const Vector3D operator+(const Vector3D &lhs,
                              const Vector3D &rhs);
    friend const Vector3D operator-(const Vector3D &lhs,
                              const Vector3D &rhs);
    friend const Vector3D operator*(const Vector3D &lhs,
                              const double scalar); 
    friend const Vector3D operator*(const double scalar,
                              const Vector3D &rhs); 
    friend const Vector3D operator/(const Vector3D &lhs,
                              const double scalar); //- lhs / scalar
    friend double operator%(const Vector3D &v1,
                            const Vector3D &v2); 
    friend double inner(const Vector3D v1[],
                        const Vector3D v2[], int n); 
    friend double operator%(const double scalar,
                            const Vector3D &v2); 
    friend double operator%(const Vector3D &v1,
                            const double scalar); 
    friend const Vector3D operator*(const Vector3D &v1, 
                              const Vector3D &v2); 
 
    friend void eqAx(Vector3D& v, const Matrix3D& A, const Vector3D& x);
    friend void plusEqAx(Vector3D& v, const Matrix3D& A, const Vector3D& x);
    friend void plusEqTransAx(Vector3D& v, const Matrix3D& A, const Vector3D& x);
    
    // Comparison functions
    friend bool operator==(const Vector3D &lhs, const Vector3D &rhs);
    friend bool operator!=(const Vector3D &lhs, const Vector3D &rhs);
    static double distance_between(const Vector3D& p1,
                                   const Vector3D& p2);
    int within_tolerance_box(const Vector3D &compare_to,
                             double tolerance) const;
    //- Compare two Vector3Ds to see if they are spatially equal.  
    // Return TRUE if difference in x, y, and z are all within tolerance.
    // Essentially checks to see if 'this' lies within a box centered on
    // 'compare_to' with sides of length ('tolerance' * 2).

    // Length functions
    inline double length_squared() const;
    inline double length() const;
    friend  double length(Vector3D* const v,int n); 
    friend  double Linf(Vector3D* const v,int n); 

    inline void set_length(const double new_length);
    inline void normalize();
    
    // Utility functions.  All angle functions work in radians.
    static double interior_angle(const Vector3D &a,
                                 const Vector3D &b,
                                 MsqError& err);
    //- Interior angle: acos((a%b)/(|a||b|))
    static Vector3D interpolate(const double param, const Vector3D &p1,
                                const Vector3D &p2);
    //- Interpolate between two points. Returns (1-param)*v1 + param*v2.

    const double* to_array() const 
      { return mCoords; }

  protected:
    double mCoords[3];
  };

  // Constructors
  inline Vector3D::Vector3D() 
  {
    mCoords[0] = 0;
    mCoords[1] = 0;
    mCoords[2] = 0;
  }
  inline Vector3D::Vector3D(const double &x,
                            const double &y,
                            const double &z) 
  {
    mCoords[0] = x;
    mCoords[1] = y;
    mCoords[2] = z;
  }
  inline Vector3D::Vector3D(const double xyz[3]) 
  { msq_stdc::memcpy(mCoords, xyz, 3*sizeof(double)); }
  inline Vector3D::Vector3D(const Vector3D& to_copy) 
  { msq_stdc::memcpy(mCoords, to_copy.mCoords, 3*sizeof(double)); }
  
  // Functions to get coordinates
  inline double Vector3D::x() const
  { return mCoords[0]; }
  inline double Vector3D::y() const
  { return mCoords[1]; }
  inline double Vector3D::z() const
  { return mCoords[2]; }
  inline void Vector3D::get_coordinates(double &x, double &y, double &z) const
  {
    x = mCoords[0]; 
    y = mCoords[1]; 
    z = mCoords[2];
  }
  inline void Vector3D::get_coordinates(double xyz[3]) const
  { msq_stdc::memcpy(xyz, mCoords, 3*sizeof(double)); }
  inline const double& Vector3D::operator[](size_t index) const
  {
    return mCoords[index];
  }
  
  // Functions to set coordinates
  inline void Vector3D::x( const double x )
  { mCoords[0] = x; }
  inline void Vector3D::y( const double y )
  { mCoords[1] = y; }
  inline void Vector3D::z( const double z )
  { mCoords[2] = z; }
  inline void Vector3D::set(const double x,
                            const double y,
                            const double z)
  {
    mCoords[0] = x;
    mCoords[1] = y;
    mCoords[2] = z;
  }
  inline void Vector3D::set(const double xyz[3])
  { msq_stdc::memcpy(mCoords, xyz, 3*sizeof(double)); }
  inline void Vector3D::set(const Vector3D& to_copy)
  { msq_stdc::memcpy(mCoords, to_copy.mCoords, 3*sizeof(double)); }
  inline double& Vector3D::operator[](size_t index)
  { return mCoords[index]; }
   
  inline Vector3D& Vector3D::operator=(const Vector3D &to_copy)  
  {
    mCoords[0] = to_copy.mCoords[0];
    mCoords[1] = to_copy.mCoords[1];
    mCoords[2] = to_copy.mCoords[2];
//    memcpy(mCoords, to_copy.mCoords, 3*sizeof(double));
    return *this;
  }

  inline Vector3D& Vector3D::operator=(const double &to_copy)  
  {
    mCoords[0] = to_copy;
    mCoords[1] = to_copy;
    mCoords[2] = to_copy;
    return *this;
  }
  
  // Functions that modify existing coordinates
  inline Vector3D Vector3D::operator-() const
  {
    return Vector3D(-mCoords[0], -mCoords[1], -mCoords[2]);
  }
  inline Vector3D& Vector3D::operator*=(const double scalar)
  {
    mCoords[0] *= scalar;
    mCoords[1] *= scalar;
    mCoords[2] *= scalar;
    return *this;
  }
  inline Vector3D& Vector3D::operator/=(const double scalar)
  {
    mCoords[0] /= scalar;
    mCoords[1] /= scalar;
    mCoords[2] /= scalar;
    return *this;
  }
  inline Vector3D& Vector3D::operator*=(const Vector3D &rhs)
  {
    double new_coords[3] = 
      {mCoords[1]*rhs.mCoords[2] - mCoords[2]*rhs.mCoords[1],
       mCoords[2]*rhs.mCoords[0] - mCoords[0]*rhs.mCoords[2],
       mCoords[0]*rhs.mCoords[1] - mCoords[1]*rhs.mCoords[0]
      };
    msq_stdc::memcpy(mCoords, new_coords, 3*sizeof(double));
    return *this;
  }
  inline Vector3D& Vector3D::operator+=(const Vector3D &rhs)
  {
    mCoords[0] += rhs.mCoords[0];
    mCoords[1] += rhs.mCoords[1];
    mCoords[2] += rhs.mCoords[2];
    return *this;
  }
  inline Vector3D& Vector3D::operator-=(const Vector3D &rhs)
  {
    mCoords[0] -= rhs.mCoords[0];
    mCoords[1] -= rhs.mCoords[1];
    mCoords[2] -= rhs.mCoords[2];
    return *this;
  }

  // Binary operators
  inline const Vector3D operator+(const Vector3D &lhs,
                            const Vector3D &rhs)
  {
    return Vector3D(lhs.x() + rhs.x(),
                    lhs.y() + rhs.y(),
                    lhs.z() + rhs.z());
  }
  inline const Vector3D operator-(const Vector3D &lhs,
                            const Vector3D &rhs)
  {
    return Vector3D(lhs.x() - rhs.x(),
                    lhs.y() - rhs.y(),
                    lhs.z() - rhs.z());
  }
  inline const Vector3D operator*(const Vector3D &lhs,
                            const double scalar)
  {
    return Vector3D(lhs.x() * scalar,
                    lhs.y() * scalar,
                    lhs.z() * scalar);
  }
  inline const Vector3D operator*(const double scalar,
                            const Vector3D &rhs)
  {
    return Vector3D(rhs.x() * scalar,
                    rhs.y() * scalar,
                    rhs.z() * scalar);
  }
  inline const Vector3D operator/(const Vector3D &lhs,
                            const double scalar)
  {
    assert (scalar != 0);
    return Vector3D(lhs.x() / scalar,
                    lhs.y() / scalar,
                    lhs.z() / scalar);
  }
  inline double operator%(const Vector3D &lhs,
                          const Vector3D &rhs) // Dot Product
  {
    return( lhs.mCoords[0] * rhs.mCoords[0] +
            lhs.mCoords[1] * rhs.mCoords[1] +
            lhs.mCoords[2] * rhs.mCoords[2] );
  }

  inline double inner(const Vector3D lhs[],
                      const Vector3D rhs[], int n)
  {
    int i;
    double dot=0;
    for (i=0; i<n; ++i) 
      dot+= lhs[i].mCoords[0] * rhs[i].mCoords[0] +
            lhs[i].mCoords[1] * rhs[i].mCoords[1] +
            lhs[i].mCoords[2] * rhs[i].mCoords[2];
    return dot;
  }

  inline double operator%(const double scalar,
                          const Vector3D &rhs) // Dot Product
  {
    return( scalar * (rhs.mCoords[0] + rhs.mCoords[1] + rhs.mCoords[2]) );
  }
  inline double operator%(const Vector3D &lhs,
                          const double scalar) // Dot Product
  {
    return( scalar * (lhs.mCoords[0] + lhs.mCoords[1] + lhs.mCoords[2]) );
  }
  inline const Vector3D operator*(const Vector3D &lhs, 
                            const Vector3D &rhs) // Cross Product
  {
    return Vector3D(lhs.mCoords[1]*rhs.mCoords[2]-lhs.mCoords[2]*rhs.mCoords[1],
                    lhs.mCoords[2]*rhs.mCoords[0]-lhs.mCoords[0]*rhs.mCoords[2],
                    lhs.mCoords[0]*rhs.mCoords[1]-lhs.mCoords[1]*rhs.mCoords[0]);
  }
  
  // output operator
  msq_stdio::ostream& operator<<(msq_stdio::ostream &s, const Mesquite::Vector3D &v);
  
  inline double Vector3D::distance_between(const Vector3D &p1,
                                           const Vector3D &p2)
  {
    Vector3D v = p2 - p1;
    return v.length();
  }
  inline int Vector3D::within_tolerance_box(const Vector3D &compare_to,
                                            double tolerance) const
  {
    return ((msq_stdc::fabs(this->mCoords[0] - compare_to.mCoords[0]) < tolerance) &&
            (msq_stdc::fabs(this->mCoords[1] - compare_to.mCoords[1]) < tolerance) &&
            (msq_stdc::fabs(this->mCoords[2] - compare_to.mCoords[2]) < tolerance));
  }
  
  // Length functions
  inline double Vector3D::length_squared() const
  {
    return (mCoords[0]*mCoords[0] +
            mCoords[1]*mCoords[1] +
            mCoords[2]*mCoords[2]);
  }
  inline double Vector3D::length() const
  {
    return msq_stdc::sqrt(mCoords[0]*mCoords[0] +
                          mCoords[1]*mCoords[1] +
                          mCoords[2]*mCoords[2]);
  }
  inline double length(Vector3D* const v,int n) // norm for an array of Vector3Ds
  {
    double l=0;
    for (int j=0; j<n; ++j)
      l += v[j].mCoords[0]*v[j].mCoords[0] +
           v[j].mCoords[1]*v[j].mCoords[1] +
           v[j].mCoords[2]*v[j].mCoords[2];
    return msq_stdc::sqrt(l);
  }
  inline double Linf(Vector3D* const v,int n) // max entry for an array of Vector3Ds
  {
    double max=0;  
    //loop over the length of the array
    for(int i=0;i<n;++i){
      if ( max < msq_stdc::fabs(v[i][0]) )   max=msq_stdc::fabs(v[i][0]) ;
      if ( max < msq_stdc::fabs(v[i][1]) )   max=msq_stdc::fabs(v[i][1]) ;
      if ( max < msq_stdc::fabs(v[i][2]) )   max=msq_stdc::fabs(v[i][2]) ;
    }
    //return the value of the largest entry in the array
    return max;
  }
  inline void Vector3D::set_length(const double new_length)
  {
    double factor = new_length / length();
    *this *= factor;
  }
  inline void Vector3D::normalize()
  { set_length(1.0); }

  // Utility functions.
  inline Vector3D Vector3D::interpolate(const double param,
                                        const Vector3D &p1,
                                        const Vector3D &p2)
  {
    return (1-param)*p1 + param*p2;
  }

  inline bool operator==( const Vector3D& v1, const Vector3D&v2 )
    { return v1.mCoords[0] == v2.mCoords[0] &&
             v1.mCoords[1] == v2.mCoords[1] &&
             v1.mCoords[2] == v2.mCoords[2]; }

  inline bool operator!=( const Vector3D& v1, const Vector3D&v2 )
    { return v1.mCoords[0] != v2.mCoords[0] ||
             v1.mCoords[1] != v2.mCoords[1] ||
             v1.mCoords[2] != v2.mCoords[2]; }

} // namespace Mesquite

#endif