Tensor.hpp

Go to the documentation of this file.

// This file is part of GranOO, a workbench for DEM simulation.
//
// Author(s)    : - Damien Andre       IRCER/UNILIM, Limoges France
//                  <damien.andre@unilim.fr>
//                - Jean-luc Charles   Arts et Metiers ParisTech, CNRS, I2M, Bordeaux France
//                  <jean-luc.charles@ensam.eu>
//                - Jeremie Girardot   Arts et Metiers ParisTech, CNRS, I2M, Bordeaux France
//                  <jeremie.girardot@ensam.eu>
//                - Cedric Hubert      LAMIH/UPHF, Valenciennes France
//                  <cedric.hubert@uphf.fr>
//                - Ivan Iordanoff     Arts et Metiers ParisTech, CNRS, I2M, Bordeaux France
//                  <ivan.iordanoff@ensam.eu>
//
// Copyright (C) 2008-2019 D. Andre, JL. Charles, J. Girardot, C. Hubert, I. Iordanoff
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
#ifndef _GranOO_Geom_Tensor_hpp_
#define _GranOO_Geom_Tensor_hpp_
 
#include <cmath>
#include <iostream>
#include "GranOO3/3rdParty/Eigen/Dense"
 
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
 
#include "GranOO3/Geom/Vector.hpp"
 
 
namespace GranOO3
{
  namespace Geom
  {
 
    class Frame;
 
    std::ostream& operator<< (std::ostream& o, const Tensor& m);
    Vector operator* (const Tensor &m, const Vector &v);
    Tensor operator* (const double &d, const Tensor &m);
    Tensor operator* (const Tensor &m, const double &d);
    Tensor operator* (const Tensor &m1, const Tensor &m2);
    Tensor operator/ (const Tensor &m, const double &d);
    Tensor operator- (const Tensor &m1, const Tensor &m2);
    Tensor operator+ (const Tensor &m1, const Tensor &m2);
    bool  operator== (const Tensor &m1, const Tensor &m2);
 
 
    class Tensor
    {
    public:
      // Usefull for internal use only, NOT DOCUMENTED !
      static std::string class_ID() {return "Tensor";}
      static const int N = 9;
 
    public:
      // Constructors and destructor
      Tensor();
      Tensor(const Tensor& m);
      Tensor(const Tensor& m, const Frame& from, const Frame& to);
 
      Tensor(double, double, double,
             double, double, double,
             double, double, double);
 
      // Internal operators
      Tensor& operator=(const Tensor& m);
      void operator+= (const Tensor& m);
      void operator/= (double d);
 
      // Usefull methods
      void clear();
      void eigen_value(Vector& v, Quaternion& q) const;
      double determinant() const;
      Tensor transpose() const;
      Tensor inverse() const;
      bool is_null() const;
      bool is_nan() const;
 
      // Component
      double & xx();
      double & xy();
      double & xz();
 
      double & yx();
      double & yy();
      double & yz();
 
      double & zx();
      double & zy();
      double & zz();
 
      const double & xx() const;
      const double & xy() const;
      const double & xz() const;
 
      const double & yx() const;
      const double & yy() const;
      const double & yz() const;
 
      const double & zx() const;
      const double & zy() const;
      const double & zz() const;
 
 
      // Component (for python wrapper)
      void set_xx(const double&);
      void set_xy(const double&);
      void set_xz(const double&);
      void set_yx(const double&);
      void set_yy(const double&);
      void set_yz(const double&);
      void set_zx(const double&);
      void set_zy(const double&);
      void set_zz(const double&);
 
      // Component (for python wrapper)
      double get_xx() const;
      double get_xy() const;
      double get_xz() const;
      double get_yx() const;
      double get_yy() const;
      double get_yz() const;
      double get_zx() const;
      double get_zy() const;
      double get_zz() const;
 
      template <typename Axis> double& val();
      template <typename Axis> const double& val() const;
 
      const double & operator() (unsigned int i, unsigned int j) const;
      double& operator() (unsigned int i, unsigned int j);
 
      double trace() const;
 
 
    private:
      //BOOST SERIALIZATION
      friend class boost::serialization::access;
      template<class Archive>
      void serialize(Archive & ar, const unsigned int);
 
    public:
      Eigen::Matrix<double, 3, 3> coord;
    };
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
 
    template<class Archive>
    void Tensor::serialize(Archive & ar, const unsigned int ) {
      ar & xx(); ar & yx(); ar & zx();
      ar & xy(); ar & yy(); ar & zy();
      ar & xz(); ar & yz(); ar & zz();
    }
 
    inline
    Tensor::Tensor(const Tensor& M)
      : coord(M.coord) {
    }
 
    inline
    Tensor::Tensor(double XX, double XY, double XZ,
                   double YX, double YY, double YZ,
                   double ZX, double ZY, double ZZ)
      : coord(){
      xx() = XX; xy() = XY;  xz() = XZ;
      yx() = YX; yy() = YY;  yz() = YZ;
      zx() = ZX; zy() = ZY;  zz() = ZZ;
    }
 
    inline
    Tensor::Tensor()
      : coord() {
    }
 
    // Operator
    inline Tensor
    operator/(const Tensor &m, const double &d) {
      return Tensor(m.xx()/d, m.xy()/d, m.xz()/d,
                    m.yx()/d, m.yy()/d, m.yz()/d,
                    m.zx()/d, m.zy()/d, m.zz()/d);
    }
 
    inline Tensor &
    Tensor::operator=(const Tensor& M) {
      coord = M.coord;
      return *this;
    }
 
    inline double
    Tensor::determinant() const {
      return (xx()*(yy()*zz()-yz()*zy()) - xy()*(yx()*zz()-yz()*zx()) + xz()*(yx()*zy()-yy()*zx()));
    }
 
    inline Tensor
    Tensor::transpose() const {
      return Tensor(xx(), yx(), zx(), xy(), yy(), zy(), xz(), yz(), zz());
    }
 
    inline Tensor
    Tensor::inverse() const {
      const double det = determinant();
      return Tensor(yy()*zz() - zy()*yz(),
                    zx()*yz() - yx()*zz(),
                    yx()*zy() - zx()*yy(),
                    zy()*xz() - xy()*zz(),
                    xx()*zz() - zx()*xz(),
                    zx()*xy() - xx()*zy(),
                    xy()*yz() - yy()*xz(),
                    yx()*xz() - xx()*yz(),
                    xx()*yy() - yx()*xy()
                    ).transpose()/det;
    }
 
    inline void
    Tensor::clear() {
      xx() = xy() = xz() = 0.;
      yx() = yy() = yz() = 0.;
      zx() = zy() = zz() = 0.;
    }
 
    inline bool
    Tensor::is_null() const {
      return (xy() == 0. && xz() == 0. && xz() == 0. &&
              yy() == 0. && yz() == 0. && yz() == 0. &&
              zy() == 0. && zz() == 0. && zz() == 0);
    }
 
    inline bool
    Tensor::is_nan() const {
      if (xx()!=xx())
        return true;
      if (xy()!=xy())
        return true;
      if (xz()!=xz())
        return true;
 
      if (yx()!=yx())
        return true;
      if (yy()!=yy())
        return true;
      if (yz()!=yz())
        return true;
 
      if (zx()!=zx())
        return true;
      if (zy()!=zy())
        return true;
      if (zz()!=zz())
        return true;
 
      return false;
    }
 
    inline double & Tensor::xx() {
      return coord(0,0);
    }
 
    inline double & Tensor::xy() {
      return coord(0,1);
    }
 
    inline double & Tensor::xz() {
      return coord(0,2);
    }
 
    inline double & Tensor::yx() {
      return coord(1,0);
    }
 
    inline double & Tensor::yy() {
      return coord(1,1);
    }
 
    inline double & Tensor::yz() {
      return coord(1,2);
    }
 
    inline double & Tensor::zx() {
      return coord(2,0);
    }
 
    inline double & Tensor::zy() {
      return coord(2,1);
    }
 
    inline double & Tensor::zz() {
      return coord(2,2);
    }
 
    inline const double & Tensor::xx() const {
      return coord(0,0);
    }
 
    inline const double & Tensor::xy() const {
      return coord(0,1);
    }
 
    inline const double & Tensor::xz() const {
      return coord(0,2);
    }
 
    inline const double & Tensor::yx() const {
      return coord(1,0);
    }
 
    inline const double & Tensor::yy() const {
      return coord(1,1);
    }
 
    inline const double & Tensor::yz() const {
      return coord(1,2);
    }
 
    inline const double & Tensor::zx() const {
      return coord(2,0);
    }
 
    inline const double & Tensor::zy() const {
      return coord(2,1);
    }
 
    inline const double & Tensor::zz() const {
      return coord(2,2);
    }
 
    template<> inline double& Tensor::val<XX>() {
      return xx();
    }
 
    template<> inline double& Tensor::val<XY>() {
      return xy();
    }
 
    template<> inline double& Tensor::val<XZ>() {
      return xz();
    }
 
    template<> inline double& Tensor::val<YX>() {
      return yx();
    }
 
    template<> inline double& Tensor::val<YY>() {
      return yy();
    }
 
    template<> inline double& Tensor::val<YZ>() {
      return yz();
    }
 
    template<> inline double& Tensor::val<ZX>() {
      return zx();
    }
 
    template<> inline double& Tensor::val<ZY>() {
      return zy();
    }
 
    template<> inline double& Tensor::val<ZZ>() {
      return zz();
    }
 
    template<> inline const double& Tensor::val<XX>() const {
      return xx();
    }
 
    template<> inline const double& Tensor::val<XY>() const {
      return xy();
    }
 
    template<> inline const double& Tensor::val<XZ>() const {
      return xz();
    }
 
    template<> inline const double& Tensor::val<YX>() const {
      return yx();
    }
 
    template<> inline const double& Tensor::val<YY>() const {
      return yy();
    }
 
    template<> inline const double& Tensor::val<YZ>() const {
      return yz();
    }
 
    template<> inline const double& Tensor::val<ZX>() const {
      return zx();
    }
 
    template<> inline const double& Tensor::val<ZY>() const {
      return zy();
    }
 
    template<> inline const double& Tensor::val<ZZ>() const {
      return zz();
    }
 
    inline void Tensor::set_xx(const double& val) {
      xx() = val;
    }
 
    inline void Tensor::set_xy(const double& val) {
      xy() = val;
    }
 
    inline void Tensor::set_xz(const double& val) {
      xz() = val;
    }
 
    inline void Tensor::set_yx(const double& val) {
      yx() = val;
    }
 
    inline void Tensor::set_yy(const double& val) {
      yy() = val;
    }
 
    inline void Tensor::set_yz(const double& val) {
      yz() = val;
    }
 
    inline void Tensor::set_zx(const double& val) {
      zx() = val;
    }
 
    inline void Tensor::set_zy(const double& val) {
      zy() = val;
    }
 
    inline void Tensor::set_zz(const double& val) {
      zz() = val;
    }
 
    inline double Tensor::get_xx() const {
      return xx();
    }
 
    inline double Tensor::get_xy() const {
      return xy();
    }
 
    inline double Tensor::get_xz() const {
      return xz();
    }
 
    inline double Tensor::get_yx() const {
      return yx();
    }
 
    inline double Tensor::get_yy() const {
      return yy();
    }
 
    inline double Tensor::get_yz() const {
      return yz();
    }
 
    inline double Tensor::get_zx() const {
      return zx();
    }
 
    inline double Tensor::get_zy() const {
      return zy();
    }
 
    inline double Tensor::get_zz() const {
      return zz();
    }
 
    inline const double&
    Tensor::operator() (unsigned int i, unsigned int j) const {
      SafeModeAssert(i < 3, "The component must be in [0, 2] range");
      SafeModeAssert(j < 3, "The component must be in [0, 2] range");
      return coord(i,j);
    }
 
    inline double&
    Tensor::operator() (unsigned int i, unsigned int j) {
      SafeModeAssert(i < 3, "The component must be in [0, 2] range");
      SafeModeAssert(j < 3, "The component must be in [0, 2] range");
      return coord(i,j);
    }
 
    inline double Tensor::trace() const {
      return xx() + yy() + zz();
    }
 
 
    //
    //EXTERNALS OPERATORS
    //
    inline
    std::ostream&
    operator<< (std::ostream& o, const Tensor& M) {
      return o << M.xx() << '\t' << M.xy() << '\t' << M.xz() << '\n'
               << M.yx() << '\t' << M.yy() << '\t' << M.yz() << '\n'
               << M.zx() << '\t' << M.zy() << '\t' << M.zz() << '\n';
    }
 
    inline Vector
    operator*(const Tensor &m, const Vector &v) {
      return Vector(m.xx()*v.x() + m.xy()*v.y() + m.xz()*v.z(),
                    m.yx()*v.x() + m.yy()*v.y() + m.yz()*v.z(),
                    m.zx()*v.x() + m.zy()*v.y() + m.zz()*v.z());
 
    }
 
    inline Tensor
    operator*(const double &d, const Tensor &m) {
      return Tensor(m.xx()*d, m.xy()*d, m.xz()*d,
                    m.yx()*d, m.yy()*d, m.yz()*d,
                    m.zx()*d, m.zy()*d, m.zz()*d);
    }
 
    inline Tensor
    operator*(const Tensor &m, const double &d) {
      return d*m;
    }
 
    inline Tensor
    operator* (const Tensor &m1, const Tensor &m2) {
      return Tensor(m1.xx()*m2.xx() + m1.xy()*m2.yx() + m1.xz()*m2.zx(),
                    m1.xx()*m2.xy() + m1.xy()*m2.yy() + m1.xz()*m2.zy(),
                    m1.xx()*m2.xz() + m1.xy()*m2.yz() + m1.xz()*m2.zz(),
 
                    m1.yx()*m2.xx() + m1.yy()*m2.yx() + m1.yz()*m2.zx(),
                    m1.yx()*m2.xy() + m1.yy()*m2.yy() + m1.yz()*m2.zy(),
                    m1.yx()*m2.xz() + m1.yy()*m2.yz() + m1.yz()*m2.zz(),
 
                    m1.zx()*m2.xx() + m1.zy()*m2.yx() + m1.zz()*m2.zx(),
                    m1.zx()*m2.xy() + m1.zy()*m2.yy() + m1.zz()*m2.zy(),
                    m1.zx()*m2.xz() + m1.zy()*m2.yz() + m1.zz()*m2.zz());
    }
 
    inline Tensor
    operator+ (const Tensor &m1, const Tensor &m2) {
      return Tensor(m1.xx()+m2.xx(), m1.xy()+m2.xy(), m1.xz()+m2.xz(),
                    m1.yx()+m2.yx(), m1.yy()+m2.yy(), m1.yz()+m2.yz(),
                    m1.zx()+m2.zx(), m1.zy()+m2.zy(), m1.zz()+m2.zz());
    }
 
    inline Tensor
    operator- (const Tensor &m1, const Tensor &m2) {
      return Tensor(m1.xx()-m2.xx(), m1.xy()-m2.xy(), m1.xz()-m2.xz(),
                    m1.yx()-m2.yx(), m1.yy()-m2.yy(), m1.yz()-m2.yz(),
                    m1.zx()-m2.zx(), m1.zy()-m2.zy(), m1.zz()-m2.zz());
    }
 
    inline bool
    operator== (const Tensor &m1, const Tensor &m2) {
       Tensor t(m1.xx()-m2.xx(), m1.xy()-m2.xy(), m1.xz()-m2.xz(),
        m1.yx()-m2.yx(), m1.yy()-m2.yy(), m1.yz()-m2.yz(),
        m1.zx()-m2.zx(), m1.zy()-m2.zy(), m1.zz()-m2.zz());
       return t.is_null();
    }
 
    inline void
    Tensor::operator+= (const Tensor &m1) {
      Tensor& m = *this;
      m.xx() += m1.xx();
      m.yy() += m1.yy();
      m.zz() += m1.zz();
      m.xy() += m1.xy();
      m.yx() += m1.yx();
      m.xz() += m1.xz();
      m.zx() += m1.zx();
      m.yz() += m1.yz();
      m.zy() += m1.zy();
    }
 
    inline void
    Tensor::operator/= (double d) {
      Tensor& m = *this;
      m.xx() /= d;
      m.yy() /= d;
      m.zz() /= d;
      m.xy() /= d;
      m.yx() /= d;
      m.xz() /= d;
      m.zx() /= d;
      m.yz() /= d;
      m.zy() /= d;
    }
 
#endif
  }
}
 
 
 
#endif