35 #ifndef CFDEM_MATH_EXTRA_H
36 #define CFDEM_MATH_EXTRA_H
45 #define TOLERANCE_ORTHO 1e-10
57 inline void outerProduct(
double *vec1,
double *vec2,
double **m)
61 for( i = 0; i < 3; ++i )
63 for( i = 0; i < 3; ++i )
67 for( i = 0; i < 3; ++i )
68 for( j = 0; j < 3; ++j )
70 m[i][j] = vec1[i] * vec2[j];
71 printf(
"OUTER PRODUCT: Result: m[%d][%d]=%g", i, j, m[i][j]);
78 inline bool spheroidGeometry(
double radius,
double aspectRatio,
79 double& ai,
double& bi,
80 double& ei,
double& Le
96 ai = radius * std::pow(aspectRatio*aspectRatio,0.33333333333333333333333);
97 bi = ai / aspectRatio;
100 - 1.0 / (aspectRatio*aspectRatio)
114 return 3.1415926535897932384626433832795;
120 inline bool spheroidGeometry2(
double radius,
double aspectRatio,
121 double& ai,
double& bi,
122 double& XAe,
double& YAe,
123 double& XCe,
double& YCe,
140 double ei(0.0), Le(0.0);
142 spheroidGeometry(radius, aspectRatio,
149 XAe= 2.6666666666666666666666667
151 /(-2.0*ei+(1.0+ei*ei)*Le);
152 YAe= 5.333333333333333333333333333
154 /(2.0*ei+(3*ei*ei-1.0)*Le);
155 XCe= 1.333333333333333333333333333
158 /(2.0*ei-(1.0-ei*ei)*Le);
159 YCe= 1.3333333333333333333333
162 /(-2.0*ei+(1.0+ei*ei)*Le);
163 YHe= 1.3333333333333333333333
165 /(-2.0*ei+(1.0+ei*ei)*Le);
173 inline void zeroize333(
double tensor[3][3][3] )
175 for(
int iX=0; iX<3; iX++)
176 for(
int iY=0; iY<3; iY++)
177 for(
int iZ=0; iZ<3; iZ++)
178 tensor[iX][iY][iZ] = 0.0;
183 inline void zeroize33(
double tensor[3][3] )
185 for(
int iX=0; iX<3; iX++)
186 for(
int iY=0; iY<3; iY++)
187 tensor[iX][iY] = 0.0;
192 inline void multiply333(
double scalar,
double tensor[3][3][3] )
194 for(
int iX=0; iX<3; iX++)
195 for(
int iY=0; iY<3; iY++)
196 for(
int iZ=0; iZ<3; iZ++)
197 tensor[iX][iY][iZ] *= scalar;
201 inline void permutationTensor(
double tensor[3][3][3] )
204 tensor[0][1][2] = 1.0;
205 tensor[1][2][0] = 1.0;
206 tensor[2][0][1] = 1.0;
207 tensor[0][2][1] =-1.0;
208 tensor[2][1][0] =-1.0;
209 tensor[1][0][2] =-1.0;
217 inline bool permutationDotDyadic(
219 double tensor[3][3][3]
223 double permutationT[3][3][3];
224 permutationTensor(permutationT);
227 double permutationDotProd[3][3];
228 zeroize33(permutationDotProd);
230 for(
int iX=0; iX<3; iX++)
231 for(
int iY=0; iY<3; iY++)
232 for(
int iZ=0; iZ<3; iZ++)
233 permutationDotProd[iX][iY] += permutationT[iX][iY][iZ]
236 for(
int iX=0; iX<3; iX++)
237 for(
int iY=0; iY<3; iY++)
238 for(
int iZ=0; iZ<3; iZ++)
239 tensor[iX][iY][iZ] = permutationDotProd[iX][iY]
247 inline bool doubleDotTensor333Tensor33(
double tensor333[3][3][3],
248 double tensor33[3][3],
252 result[0]=0.0;result[1]=0.0;result[2]=0.0;
254 for(
int iX=0; iX<3; iX++)
255 for(
int iY=0; iY<3; iY++)
256 for(
int iZ=0; iZ<3; iZ++)
257 result[iX] += tensor333[iX][iY][iZ] * tensor33[iY][iZ];