SymTensor.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_SymTensor_hpp_
#define _GranOO_Geom_SymTensor_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"
#include "GranOO3/Geom/Tensor.hpp"
 
 
namespace GranOO3
{
  namespace Geom
  {
 
    class Frame;
    class SymTensor;
 
    std::ostream& operator<< (std::ostream& o, const SymTensor& m);
    Vector operator* (const SymTensor &m, const Vector &v);
    SymTensor operator* (const double &d, const SymTensor &m);
    SymTensor operator* (const SymTensor &m, const double &d);
    double operator* (const SymTensor &m1, const SymTensor &m2); // dot product
    Tensor operator* (const SymTensor &m1, const Tensor &m2);
    Tensor operator* (const Tensor &m1, const SymTensor &m2);
    SymTensor operator/ (const SymTensor &m, const double &d);
    SymTensor operator- (const SymTensor &m1, const SymTensor &m2);
    SymTensor operator+ (const SymTensor &m1, const SymTensor &m2);
    bool      operator==(const SymTensor &m1, const SymTensor &m2);
 
    // A tensor that uses the Voigt's notation for symetric tensor    //
    // The tensor is stored as a 6 components vector in the following //
    // order : xx, yy, zz, xy, xz, yz                                 //
    class SymTensor
    {
    public:
      // Usefull for internal use only, NOT DOCUMENTED !
      static std::string class_ID() {return "SymTensor";}
      static const int N = 6;
      static const SymTensor id;
 
    public:
      // Constructors and destructor
      SymTensor();
      SymTensor(const SymTensor& m);
      SymTensor(const SymTensor& m, const Frame& from, const Frame& to);
 
      SymTensor(double, double, double,
        double, double, double);
 
      // Internal operators
      SymTensor& operator=(const SymTensor& m);
      SymTensor& operator=(const Tensor& m);
      void operator+= (const SymTensor& m);
      void operator/= (double d);
 
      // Usefull methods
      void clear();
      void eigen_value(Vector& v, Quaternion& q) const;
      bool is_null() const;
      bool is_nan() const;
 
      // Invariants
      double trace() const;
      double invariant1() const;
      double invariant2() const;
      double invariant3() 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);
 
      const double & operator() (unsigned int i) const;
      double& operator() (unsigned int i);
 
    private:
      //BOOST SERIALIZATION
      friend class boost::serialization::access;
      template<class Archive>
      void serialize(Archive & ar, const unsigned int);
 
    public:
      Eigen::Matrix<double, 6, 1> coord;
    };
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
 
    template<class Archive>
    void SymTensor::serialize(Archive & ar, const unsigned int ) {
      ar & xx(); ar & yy(); ar & zz();
      ar & xy(); ar & xz(); ar & yz();
    }
 
    inline
    SymTensor::SymTensor(const SymTensor& M)
      : coord(M.coord) {
    }
 
    inline
    SymTensor::SymTensor(double XX, double YY, double ZZ,
             double XY, double XZ, double YZ)
      : coord(){
      coord(0) = XX;
      coord(1) = YY;
      coord(2) = ZZ;
      coord(3) = XY;
      coord(4) = XZ;
      coord(5) = YZ;
    }
 
    inline
    SymTensor::SymTensor()
      : coord() {
    }
 
    // Operator
    inline SymTensor
    operator/(const SymTensor &m, const double &d) {
      return SymTensor(m.xx()/d, m.yy()/d, m.zz()/d,
               m.xy()/d, m.xz()/d, m.yz()/d);
    }
    
    inline SymTensor &
    SymTensor::operator=(const SymTensor& M) {
      coord = M.coord;
      return *this;
    }
 
 
    inline SymTensor &
    SymTensor::operator=(const Tensor& M) {
      xx() = M.xx();
      yy() = M.yy();
      zz() = M.zz();
 
      xy() = M.xy();
      xz() = M.xz();
      yz() = M.yz();
 
      return *this;
    }
 
    
    inline void
    SymTensor::clear() {
      xx() = xy() = xz() = 0.;
      yx() = yy() = yz() = 0.;
      zx() = zy() = zz() = 0.;
    }
 
    inline bool
    SymTensor::is_null() const {
      return (xy() == 0. && xz() == 0. && xz() == 0. &&
              yy() == 0. && yz() == 0. && yz() == 0. &&
              zy() == 0. && zz() == 0. && zz() == 0);
    }
 
    inline bool
    SymTensor::is_nan() const {
      if (xx()!=xx())
        return true;
      if (yy()!=yy())
        return true;
      if (zz()!=zz())
        return true;
 
      if (xy()!=xy())
        return true;
      if (xz()!=xz())
        return true;
      if (yz()!=yz())
        return true;
 
      return false;
    }
 
    inline double & SymTensor::xx() {
      return coord(0);
    }
 
    inline double & SymTensor::xy() {
      return coord(3);
    }
 
    inline double & SymTensor::xz() {
      return coord(4);
    }
 
    inline double & SymTensor::yx() {
      return coord(3);
    }
 
    inline double & SymTensor::yy() {
      return coord(1);
    }
 
    inline double & SymTensor::yz() {
      return coord(5);
    }
 
    inline double & SymTensor::zx() {
      return coord(4);
    }
 
    inline double & SymTensor::zy() {
      return coord(5);
    }
 
    inline double & SymTensor::zz() {
      return coord(2);
    }
 
    inline const double & SymTensor::xx() const {
      return coord(0);
    }
 
    inline const double & SymTensor::xy() const {
      return coord(3);
    }
 
    inline const double & SymTensor::xz() const {
      return coord(4);
    }
 
    inline const double & SymTensor::yx() const {
      return coord(3);
    }
 
    inline const double & SymTensor::yy() const {
      return coord(1);
    }
 
    inline const double & SymTensor::yz() const {
      return coord(5);
    }
 
    inline const double & SymTensor::zx() const {
      return coord(4);
    }
 
    inline const double & SymTensor::zy() const {
      return coord(5);
    }
 
    inline const double & SymTensor::zz() const {
      return coord(2);
    }
 
    template<> inline double& SymTensor::val<XX>() {
      return xx();
    }
 
    template<> inline double& SymTensor::val<XY>() {
      return xy();
    }
 
    template<> inline double& SymTensor::val<XZ>() {
      return xz();
    }
 
    template<> inline double& SymTensor::val<YX>() {
      return yx();
    }
 
    template<> inline double& SymTensor::val<YY>() {
      return yy();
    }
 
    template<> inline double& SymTensor::val<YZ>() {
      return yz();
    }
 
    template<> inline double& SymTensor::val<ZX>() {
      return zx();
    }
 
    template<> inline double& SymTensor::val<ZY>() {
      return zy();
    }
 
    template<> inline double& SymTensor::val<ZZ>() {
      return zz();
    }
 
    template<> inline const double& SymTensor::val<XX>() const {
      return xx();
    }
 
    template<> inline const double& SymTensor::val<XY>() const {
      return xy();
    }
 
    template<> inline const double& SymTensor::val<XZ>() const {
      return xz();
    }
 
    template<> inline const double& SymTensor::val<YX>() const {
      return yx();
    }
 
    template<> inline const double& SymTensor::val<YY>() const {
      return yy();
    }
 
    template<> inline const double& SymTensor::val<YZ>() const {
      return yz();
    }
 
    template<> inline const double& SymTensor::val<ZX>() const {
      return zx();
    }
 
    template<> inline const double& SymTensor::val<ZY>() const {
      return zy();
    }
 
    template<> inline const double& SymTensor::val<ZZ>() const {
      return zz();
    }
 
    inline void SymTensor::set_xx(const double& val) {
      xx() = val;
    }
 
    inline void SymTensor::set_xy(const double& val) {
      xy() = val;
    }
 
    inline void SymTensor::set_xz(const double& val) {
      xz() = val;
    }
 
    inline void SymTensor::set_yx(const double& val) {
      yx() = val;
    }
 
    inline void SymTensor::set_yy(const double& val) {
      yy() = val;
    }
 
    inline void SymTensor::set_yz(const double& val) {
      yz() = val;
    }
 
    inline void SymTensor::set_zx(const double& val) {
      zx() = val;
    }
 
    inline void SymTensor::set_zy(const double& val) {
      zy() = val;
    }
 
    inline void SymTensor::set_zz(const double& val) {
      zz() = val;
    }
 
    inline double SymTensor::get_xx() const {
      return xx();
    }
 
    inline double SymTensor::get_xy() const {
      return xy();
    }
 
    inline double SymTensor::get_xz() const {
      return xz();
    }
 
    inline double SymTensor::get_yx() const {
      return yx();
    }
 
    inline double SymTensor::get_yy() const {
      return yy();
    }
 
    inline double SymTensor::get_yz() const {
      return yz();
    }
 
    inline double SymTensor::get_zx() const {
      return zx();
    }
 
    inline double SymTensor::get_zy() const {
      return zy();
    }
 
    inline double SymTensor::get_zz() const {
      return zz();
    }
 
    inline const double&
    SymTensor::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&
    SymTensor::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 const double&
    SymTensor::operator() (unsigned int i) const {
      SafeModeAssert(i < 5, "The component must be in [0, 5] range");
      return coord(i);
    }
 
    inline double&
    SymTensor::operator() (unsigned int i) {
      SafeModeAssert(i < 5, "The component must be in [0, 2] range");
      return coord(i);
    }
 
 
    inline double SymTensor::trace() const {
      return xx() + yy() + zz();
    }
 
    
    inline  double SymTensor::invariant1() const {
      return xx() + yy() + zz();
    }
 
    inline double SymTensor::invariant2() const {
      return xx()*yy() + yy()*zz() + xx()*zz() - xy()*xy() - yz()*yz() - xz()*xz();
    }
 
    inline double SymTensor::invariant3() const {
      return xx()*yy()*zz() + 2*xy()*yz()*xz() - xy()*xy()*zz() - yz()*yz()*xx() - xz()*xz()*yy();
    }
 
 
    //
    //EXTERNALS OPERATORS
    //
    inline
    std::ostream&
    operator<< (std::ostream& o, const SymTensor& M) {
      return o << M.xx() << '\t' << M.yy() << '\t' << M.zz() << '\n'
               << M.xy() << '\t' << M.xz() << '\t' << M.yz() << '\n';
    }
 
    inline Vector
    operator*(const SymTensor &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 SymTensor
    operator*(const double &d, const SymTensor &m) {
      return SymTensor(m.xx()*d, m.yy()*d, m.zz()*d,
               m.xy()*d, m.xz()*d, m.yz()*d);
    }
 
    inline SymTensor
    operator*(const SymTensor &m, const double &d) {
      return d*m;
    }
 
    inline double
    operator* (const SymTensor &m1, const SymTensor &m2) {
      return (m1.xx()*m2.xx() + m1.yy()*m2.yy() + m1.zz()*m2.zz() +
          m1.xy()*m2.xy() + m1.xz()*m2.xz() + m1.yz()*m2.yz());
    }
 
    inline Tensor
    operator* (const SymTensor &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 SymTensor &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 SymTensor
    operator+ (const SymTensor &m1, const SymTensor &m2) {
      return SymTensor(m1.xx()+m2.xx(), m1.yy()+m2.yy(), m1.zz()+m2.zz(),
               m1.xy()+m2.xy(), m1.xz()+m2.xz(), m1.yz()+m2.yz());
    }
 
    inline SymTensor
    operator- (const SymTensor &m1, const SymTensor &m2) {
      return SymTensor(m1.xx()-m2.xx(), m1.yy()-m2.yy(), m1.zz()-m2.zz(),
               m1.xy()-m2.xy(), m1.xz()-m2.xz(), m1.yz()-m2.yz());
    }
 
    inline bool
    operator== (const SymTensor &m1, const SymTensor &m2) {
      SymTensor t(m1.xx()-m2.xx(), m1.yy()-m2.yy(), m1.zz()-m2.zz(),
          m1.xy()-m2.xy(), m1.xz()-m2.xz(), m1.yz()-m2.yz());
      return t.is_null();
    }
 
    inline void
    SymTensor::operator+= (const SymTensor &m1) {
      SymTensor& m = *this;
      m.xx() += m1.xx();
      m.yy() += m1.yy();
      m.zz() += m1.zz();
      m.xy() += m1.xy();
      m.xz() += m1.xz();
      m.yz() += m1.yz();
    }
 
    inline void
    SymTensor::operator/= (double d) {
      SymTensor& m = *this;
      m.xx() /= d;
      m.yy() /= d;
      m.zz() /= d;
      m.xy() /= d;
      m.xz() /= d;
      m.yz() /= d;
    }
 
#endif
  }
}
 
 
 
#endif