diff --git a/Solver/STensor33.cpp b/Solver/STensor33.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..96582fe0ae878a93f54296fd1ef93b9e3b8a5cc2
--- /dev/null
+++ b/Solver/STensor33.cpp
@@ -0,0 +1,20 @@
+#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]);
+}
+
diff --git a/Solver/STensor33.h b/Solver/STensor33.h
new file mode 100644
index 0000000000000000000000000000000000000000..7795c866b6ea452e6e9936355282185d052b5c11
--- /dev/null
+++ b/Solver/STensor33.h
@@ -0,0 +1,185 @@
+#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
+
diff --git a/Solver/STensor43.cpp b/Solver/STensor43.cpp
index cb46e4ebc1bd4ed8533c438565cd3eb14083a48b..90904f64c9954a6040fa88290fdbce7e8b91b1bc 100644
--- a/Solver/STensor43.cpp
+++ b/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]);
}
+
+
+
diff --git a/Solver/STensor43.h b/Solver/STensor43.h
index 351c39657fae03aa4b94327e976f469d3069ef43..9c983f8f65f1d2ebbd5948e65d6905e6ea8ab61b 100644
--- a/Solver/STensor43.h
+++ b/Solver/STensor43.h
@@ -1,12 +1,13 @@
+
#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
+