No commit message

No commit message

Solver/STensor33.cpp

+#include "STensor33.h"
+
+void STensor33::print (const char *s) const
+{
+  char format[512];
+  const char l[256]="%12.5E %12.5E %12.5E \n";
+  sprintf (format, " tensor3 %s : \n %s %s %s \n %s %s %s \n %s %s %s \n",s, l,l,l, l,l,l, l,l,l);
+  printf(format,s,_val[ 0],_val[ 1],_val[ 2],
+                  _val[ 3],_val[ 4],_val[ 5],
+                  _val[ 6],_val[ 7],_val[ 8],
+         
+                  _val[ 9],_val[10],_val[11],
+                  _val[12],_val[13],_val[14],
+                  _val[15],_val[16],_val[17],
+         
+                  _val[18],_val[19],_val[20],
+                  _val[21],_val[22],_val[23],
+                  _val[24],_val[25],_val[26]);
+}
+

Solver/STensor33.h

+#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 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;
+};
+
+// 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
+

Solver/STensor43.cpp

@@ -3,9 +3,21 @@
 
 void STensor43::print (const char *s) const
 {
-  printf("STensor43::print to be implemented \n");
-/*  printf(" tensor4 %s : \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n",s,
-         (*this)(0,0),(*this)(0,1),(*this)(0,2),
-         (*this)(1,0),(*this)(1,1),(*this)(1,2),
-         (*this)(2,0),(*this)(2,1),(*this)(2,2));*/
+  char format[2048];
+  const char l[256]="%12.5E %12.5E %12.5E  %12.5E %12.5E %12.5E  %12.5E %12.5E %12.5E \n";
+  sprintf (format, " tensor4 %s : \n %s %s %s \n %s %s %s \n %s %s %s \n",s, l,l,l, l,l,l, l,l,l);
+  printf(format,s,_val[ 0],_val[ 1],_val[ 2], _val[ 3],_val[ 4],_val[ 5], _val[ 6],_val[ 7],_val[ 8],
+                  _val[ 9],_val[10],_val[11], _val[12],_val[13],_val[14], _val[15],_val[16],_val[17],
+                  _val[18],_val[19],_val[20], _val[21],_val[22],_val[23], _val[24],_val[25],_val[26],
+         
+                  _val[27],_val[28],_val[29], _val[30],_val[31],_val[32], _val[33],_val[34],_val[35],
+                  _val[36],_val[37],_val[38], _val[39],_val[40],_val[41], _val[42],_val[43],_val[44],
+                  _val[45],_val[46],_val[47], _val[48],_val[49],_val[50], _val[51],_val[52],_val[53],
+         
+                  _val[54],_val[55],_val[56], _val[57],_val[58],_val[59], _val[60],_val[61],_val[62],
+                  _val[63],_val[64],_val[65], _val[66],_val[67],_val[68], _val[69],_val[70],_val[71],
+                  _val[72],_val[73],_val[74], _val[75],_val[76],_val[77], _val[78],_val[79],_val[80]);
 }
+
+
+

Solver/STensor43.h

+
 #ifndef _STENSOR43_H_
 #define _STENSOR43_H_
 
-#include "STensor3.h"
+#include "STensor33.h"
 #include "fullMatrix.h"
 #include "Numeric.h"
 
 
-// concrete class for general 3x3 matrix
+// concrete class for general 3rd-order tensor in three-dimensional space
 
 class STensor43 {
  protected:
@@ -31,13 +32,10 @@ class STensor43 {
       for (int j = 0; j < 3; j++)
         for (int k = 0; k < 3; k++)
           for (int l = 0; l < 3; l++)
-         {
-	     _val[getIndex(i, j, k, l)]= 0.;
             if ((i==k)&&(j==l))
-              _val[getIndex(i, j, k, l)]+=0.5*v;
-            if ((i==l)&&(j==k))
-              _val[getIndex(i, j, k, l)]+=0.5*v;
-	 }
+              _val[getIndex(i, j, k, l)]=v;
+            else
+              _val[getIndex(i, j, k, l)]=0.0;
   }
   inline double &operator()(int i, int j,int k, int l)
   {
@@ -63,14 +61,7 @@ class STensor43 {
     for (int i = 0; i < 81; i++) _val[i] *= other;
     return *this;
   }
-/*  STensor43& operator *= (const STensor43 &other)
-  {
-// to be implemented
-    return *this;
-  }*/
-  void print(const char *) const;
-
- STensor43 transpose (int n, int m) const
+  STensor43 transpose (int n, int m) const
   {
     STensor43 ithis;
     if ((n==0 && m==1) || (n==1 && m==0))
@@ -129,7 +120,12 @@ class STensor43 {
     }
     return ithis+=(*this);
   }
-
+/*  STensor43& operator *= (const STensor43 &other)
+  {
+// to be implemented
+    return *this;
+  }*/
+  void print(const char *) const;
 };
 
 // tensor product
@@ -142,6 +138,23 @@ inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c)
             c(i,j,k,l)=a(i,j)*b(k,l);
 }
 
+inline void tensprod(const SVector3 &a, const STensor33 &b, STensor43 &c)
+{
+    for (int i = 0; i < 3; i++)
+      for (int j = 0; j < 3; j++)
+        for (int k = 0; k < 3; k++)
+          for (int l = 0; l < 3; l++)
+            c(i,j,k,l)=a(i)*b(j,k,l);
+}
+inline void tensprod(const STensor33 &a, const SVector3 &b, STensor43 &c)
+{
+    for (int i = 0; i < 3; i++)
+      for (int j = 0; j < 3; j++)
+        for (int k = 0; k < 3; k++)
+          for (int l = 0; l < 3; l++)
+            c(i,j,k,l)=a(i,j,k)*b(l);
+}
+
 inline double dot(const STensor43 &a, const STensor43 &b)
 {
   double prod=0;
@@ -153,13 +166,13 @@ inline double dot(const STensor43 &a, const STensor43 &b)
   return prod;
 }
 
+// full contracted product
 inline STensor43 operator*(const STensor43 &t, double m)
 {
   STensor43 val(t);
   val *= m;
   return val;
 }
-
 inline STensor43 operator*(double m,const STensor43 &t)
 {
   STensor43 val(t);
@@ -167,6 +180,27 @@ inline STensor43 operator*(double m,const STensor43 &t)
   return val;
 }
 
+inline STensor33 operator*(const STensor43 &t, const SVector3 &m)
+{
+  STensor33 val(0.);
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; j++)
+      for (int k = 0; k < 3; k++)
+        for (int l = 0; l < 3; l++)
+          val(i,j,k)+=t(i,j,k,l)*m(l);
+  return val;
+}
+inline STensor33 operator*( const SVector3 &m , const STensor43 &t)
+{
+  STensor33 val(0.);
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; j++)
+      for (int k = 0; k < 3; k++)
+        for (int l = 0; l < 3; l++)
+          val(j,k,l)+=m(i)*t(i,j,k,l);
+  return val;
+}
+
 inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
 {
   STensor3 val(0.);
@@ -174,10 +208,9 @@ inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
     for (int j = 0; j < 3; j++)
       for (int k = 0; k < 3; k++)
         for (int l = 0; l < 3; l++)
-          val(i,j)+=t(i,j,k,l)*m(k,l);
+          val(i,j)+=t(i,j,k,l)*m(l,k);
   return val;
 }
-
 inline STensor3 operator*( const STensor3 &m , const STensor43 &t)
 {
   STensor3 val(0.);
@@ -185,10 +218,44 @@ inline STensor3 operator*( const STensor3 &m , const STensor43 &t)
     for (int j = 0; j < 3; j++)
       for (int k = 0; k < 3; k++)
         for (int l = 0; l < 3; l++)
-          val(k,l)+=t(i,j,k,l)*m(i,j);
+          val(k,l)+=m(j,i)*t(i,j,k,l);
   return val;
 }
 
+inline SVector3 operator*(const STensor43 &t, const STensor33 &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++)
+        for (int l = 0; l < 3; l++)
+          val(i)+=t(i,j,k,l)*m(l,k,j);
+  return val;
+}
+inline SVector3 operator*( const STensor33 &m , const STensor43 &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++)
+        for (int l = 0; l < 3; l++)
+          val(l)+=m(k,j,i)*t(i,j,k,l);
+  return val;
+}
+
+inline double operator*( const STensor43 &m , const STensor43 &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++)
+        for (int l = 0; l < 3; l++)
+          val+=m(i,j,k,l)*t(l,k,j,i);
+  return val;
+}
+
+
 
 
 #endif
+