Skip to content
Snippets Groups Projects
Select Git revision
  • 733882475c5fd9179dcab3a274f9c6a41cfca65b
  • master default protected
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

FuncGradDisc.h

Blame
  • STensor33.h 4.90 KiB
    // 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