Skip to content
Snippets Groups Projects
Commit d52f3c8d authored by Gauthier Becker's avatar Gauthier Becker
Browse files

Fix bug with STensor43 (index 77 appears twice and no index 78...)

parent 1ad6a125
No related branches found
No related tags found
No related merge requests found
...@@ -10,14 +10,14 @@ ...@@ -10,14 +10,14 @@
class STensor43 { class STensor43 {
protected: protected:
// 0000 0001 0002 0010 ... 2211 2212 2220 2221 2222 // 0000 0001 0002 0010 ... 2211 2212 2220 2221 2222
double _val[81]; double _val[81];
public: public:
inline int getIndex(int i, int j, int k, int l) const inline int getIndex(int i, int j, int k, int l) const
{ {
static int _index[3][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}}}, static int _index[3][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}}},
{{{27,28,29},{30,31,32},{33,34,35}},{{36,37,38},{39,40,41},{42,43,44}},{{45,46,47},{48,49,50},{51,52,53}}}, {{{27,28,29},{30,31,32},{33,34,35}},{{36,37,38},{39,40,41},{42,43,44}},{{45,46,47},{48,49,50},{51,52,53}}},
{{{54,55,56},{57,58,59},{60,61,62}},{{63,64,65},{66,67,68},{69,70,71}},{{72,73,74},{75,76,77},{77,79,80}}}}; {{{54,55,56},{57,58,59},{60,61,62}},{{63,64,65},{66,67,68},{69,70,71}},{{72,73,74},{75,76,77},{78,79,80}}}};
return _index[i][j][k][l]; return _index[i][j][k][l];
} }
STensor43(const STensor43& other) STensor43(const STensor43& other)
...@@ -33,7 +33,7 @@ class STensor43 { ...@@ -33,7 +33,7 @@ class STensor43 {
for (int l = 0; l < 3; l++) for (int l = 0; l < 3; l++)
if ((i==k)&&(j==l)) if ((i==k)&&(j==l))
_val[getIndex(i, j, k, l)]=v; _val[getIndex(i, j, k, l)]=v;
else else
_val[getIndex(i, j, k, l)]=0.0; _val[getIndex(i, j, k, l)]=0.0;
} }
inline double &operator()(int i, int j,int k, int l) inline double &operator()(int i, int j,int k, int l)
...@@ -43,7 +43,7 @@ class STensor43 { ...@@ -43,7 +43,7 @@ class STensor43 {
inline double operator()(int i, int j, int k, int l) const inline double operator()(int i, int j, int k, int l) const
{ {
return _val[getIndex(i, j, k ,l)]; return _val[getIndex(i, j, k ,l)];
} }
STensor43 operator + (const STensor43 &other) const STensor43 operator + (const STensor43 &other) const
{ {
STensor43 res(*this); STensor43 res(*this);
...@@ -55,12 +55,12 @@ class STensor43 { ...@@ -55,12 +55,12 @@ class STensor43 {
for (int i = 0; i < 81; i++) _val[i] += other._val[i]; for (int i = 0; i < 81; i++) _val[i] += other._val[i];
return *this; return *this;
} }
STensor43& operator *= (const double &other) STensor43& operator *= (const double &other)
{ {
for (int i = 0; i < 81; i++) _val[i] *= other; for (int i = 0; i < 81; i++) _val[i] *= other;
return *this; return *this;
} }
/* STensor43& operator *= (const STensor43 &other) /* STensor43& operator *= (const STensor43 &other)
{ {
// to be implemented // to be implemented
return *this; return *this;
...@@ -70,7 +70,7 @@ class STensor43 { ...@@ -70,7 +70,7 @@ class STensor43 {
// tensor product // tensor product
inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c) inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c)
{ {
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
for (int k = 0; k < 3; k++) for (int k = 0; k < 3; k++)
...@@ -79,8 +79,8 @@ inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c) ...@@ -79,8 +79,8 @@ inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c)
} }
inline double dot(const STensor43 &a, const STensor43 &b) inline double dot(const STensor43 &a, const STensor43 &b)
{ {
double prod=0; double prod=0;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
for (int k = 0; k < 3; k++) for (int k = 0; k < 3; k++)
...@@ -90,21 +90,21 @@ inline double dot(const STensor43 &a, const STensor43 &b) ...@@ -90,21 +90,21 @@ inline double dot(const STensor43 &a, const STensor43 &b)
} }
inline STensor43 operator*(const STensor43 &t, double m) inline STensor43 operator*(const STensor43 &t, double m)
{ {
STensor43 val(t); STensor43 val(t);
val *= m; val *= m;
return val; return val;
} }
inline STensor43 operator*(double m,const STensor43 &t) inline STensor43 operator*(double m,const STensor43 &t)
{ {
STensor43 val(t); STensor43 val(t);
val *= m; val *= m;
return val; return val;
} }
inline STensor3 operator*(const STensor43 &t, const STensor3 &m) inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
{ {
STensor3 val(0.); STensor3 val(0.);
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
...@@ -115,7 +115,7 @@ inline STensor3 operator*(const STensor43 &t, const STensor3 &m) ...@@ -115,7 +115,7 @@ inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
} }
inline STensor3 operator*( const STensor3 &m , const STensor43 &t) inline STensor3 operator*( const STensor3 &m , const STensor43 &t)
{ {
STensor3 val(0.); STensor3 val(0.);
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment