STensor33.h

// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
//
// Contributor(s):
//   Eric Bechet
//

#ifndef STENSOR33_H
#define STENSOR33_H

#include "STensor3.h"
#include "fullMatrix.h"
#include "Numeric.h"

// concrete class for general 3rd-order tensor in three-dimensional space

class STensor33 {
protected:
  // 000 001 002 010 ... 211 212 220 221 222
  double _val[27];

public:
  inline int getIndex(int i, int j, int k) const
  {
    static int _index[3][3][3] = {{{0, 1, 2}, {3, 4, 5}, {6, 7, 8}},
                                  {{9, 10, 11}, {12, 13, 14}, {15, 16, 17}},
                                  {{18, 19, 20}, {21, 22, 23}, {24, 25, 26}}};
    return _index[i][j][k];
  }
  STensor33(const STensor33 &other)
  {
    for(int i = 0; i < 27; i++) _val[i] = other._val[i];
  }
  // default constructor, null tensor
  STensor33(const double v = 0.0)
  {
    for(int i = 0; i < 3; i++)
      for(int j = 0; j < 3; j++)
        for(int k = 0; k < 3; k++) _val[getIndex(i, j, k)] = v;
  }
  inline double &operator()(int i, int j, int k)
  {
    return _val[getIndex(i, j, k)];
  }
  inline double operator()(int i, int j, int k) const
  {
    return _val[getIndex(i, j, k)];
  }
  STensor33 operator+(const STensor33 &other) const
  {
    STensor33 res(*this);
    for(int i = 0; i < 27; i++) res._val[i] += other._val[i];
    return res;
  }
  STensor33 &operator=(const STensor33 &other)
  {
    for(int i = 0; i < 27; i++) _val[i] = other._val[i];
    return *this;
  }
  STensor33 &operator+=(const STensor33 &other)
  {
    for(int i = 0; i < 27; i++) _val[i] += other._val[i];
    return *this;
  }
  STensor33 &operator-=(const STensor33 &other)
  {
    for(int i = 0; i < 27; i++) _val[i] -= other._val[i];
    return *this;
  }
  STensor33 &operator*=(const double &other)
  {
    for(int i = 0; i < 27; i++) _val[i] *= other;
    return *this;
  }
  STensor33 transpose(int n, int m) const
  {
    STensor33 ithis;
    if((n == 0 && m == 1) || (n == 1 && m == 0)) {
      for(int i = 0; i < 3; i++)
        for(int j = 0; j < 3; j++)
          for(int k = 0; k < 3; k++) ithis(i, j, k) = (*this)(j, i, k);
      return ithis;
    }
    if((n == 0 && m == 2) || (n == 2 && m == 0)) {
      for(int i = 0; i < 3; i++)
        for(int j = 0; j < 3; j++)
          for(int k = 0; k < 3; k++) ithis(i, j, k) = (*this)(k, j, i);
      return ithis;
    }
    if((n == 1 && m == 2) || (n == 2 && m == 1)) {
      for(int i = 0; i < 3; i++)
        for(int j = 0; j < 3; j++)
          for(int k = 0; k < 3; k++) ithis(i, j, k) = (*this)(i, k, j);
      return ithis;
    }
    return ithis += (*this);
  }
  /*  STensor33& operator *= (const STensor33 &other)
    {
  // to be implemented
      return *this;
    }*/
  void print(const char *) const;

  const double *data() const { return _val; }
  double *data() { return _val; }
};

// tensor product
inline void tensprod(const SVector3 &a, const STensor3 &b, STensor33 &c)
{
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) c(i, j, k) = a(i) * b(j, k);
}
inline void tensprod(const STensor3 &a, const SVector3 &b, STensor33 &c)
{
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) c(i, j, k) = a(i, j) * b(k);
}

inline double dot(const STensor33 &a, const STensor33 &b)
{
  double prod = 0;
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) prod += a(i, j, k) * b(i, j, k);
  return prod;
}

// full contracted product
inline STensor33 operator*(const STensor33 &t, double m)
{
  STensor33 val(t);
  val *= m;
  return val;
}
inline STensor33 operator*(double m, const STensor33 &t)
{
  STensor33 val(t);
  val *= m;
  return val;
}

inline STensor3 operator*(const STensor33 &t, const SVector3 &m)
{
  STensor3 val(0.);
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) val(i, j) += t(i, j, k) * m(k);
  return val;
}
inline STensor3 operator*(const SVector3 &m, const STensor33 &t)
{
  STensor3 val(0.);
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) val(j, k) += m(i) * t(i, j, k);
  return val;
}

inline SVector3 operator*(const STensor33 &t, const STensor3 &m)
{
  SVector3 val(0.);
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) val(i) += t(i, j, k) * m(k, j);
  return val;
}
inline SVector3 operator*(const STensor3 &m, const STensor33 &t)
{
  SVector3 val(0.);
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) val(k) += m(j, i) * t(i, j, k);
  return val;
}

inline double operator*(const STensor33 &m, const STensor33 &t)
{
  double val(0.);
  for(int i = 0; i < 3; i++)
    for(int j = 0; j < 3; j++)
      for(int k = 0; k < 3; k++) val += m(i, j, k) * t(k, j, i);
  return val;
}

#endif