CFDEMcoupling  2.4
 All Classes
mathExtra.H
1 /* ----------------------------------------------------------------------
2  CFDEMcoupling - Open Source CFD-DEM coupling
3 
4  CFDEMcoupling is part of the CFDEMproject
5  www.cfdem.com
6  Christoph Goniva, christoph.goniva@cfdem.com
7  Copyright 2009-2012 JKU Linz
8  Copyright 2012- DCS Computing GmbH, Linz
9 -------------------------------------------------------------------------------
10 License
11  This file is part of CFDEMcoupling.
12 
13  CFDEMcoupling is free software; you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by the
15  Free Software Foundation; either version 3 of the License, or (at your
16  option) any later version.
17 
18  CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with CFDEMcoupling; if not, write to the Free Software Foundation,
25  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
26 
27 Description
28  This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
29  and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
30 
31  Copyright of this contribution:
32  Copyright 2014- TU Graz, IPPT
33 ------------------------------------------------------------------------- */
34 
35 #ifndef CFDEM_MATH_EXTRA_H
36 #define CFDEM_MATH_EXTRA_H
37 
38 #include <vector>
39 //#include "math.h"
40 #include "stdio.h"
41 #include "string.h"
42 #include "error.h"
43 #include "ctype.h"
44 
45 #define TOLERANCE_ORTHO 1e-10
46 
47 namespace MathExtra
48 {
49 
50 // inline void outerProduct(double *vec1, double *vec2, double **m);
51 // inline double spheroidGeometry(int index, double bi, double ai);
52 
53 
54 
55 //--------------------------------------------------------------------
56 // Outer Product of two vectors
57 inline void outerProduct(double *vec1, double *vec2, double **m)
58 {
59  int i, j;
60  //debug output
61  for( i = 0; i < 3; ++i )
62 // printf("OUTER PRODUCT: Input: vec1 element %d = %g", i, vec1[i]);
63  for( i = 0; i < 3; ++i )
64 // printf("OUTER PRODUCT: Input: vec2 element %d=%g", i, vec2[i]);
65 
66  //calculation
67  for( i = 0; i < 3; ++i )
68  for( j = 0; j < 3; ++j )
69  {
70  m[i][j] = vec1[i] * vec2[j];
71  printf("OUTER PRODUCT: Result: m[%d][%d]=%g", i, j, m[i][j]);
72  }
73 
74 
75 }
76 //--------------------------------------------------------------------
77 // Compute the major, minor axis and eccentricity parameters of a prolate spheroid
78 inline bool spheroidGeometry(double radius, double aspectRatio, //inputs
79  double& ai, double& bi, //outputs
80  double& ei, double& Le //outputs
81  ) //
82 {
83  //INPUT
84  // radius ...volume-equivalent radius of the spheroid
85  // aspectRatio ...major/minor aspect ratio
86 
87  //OUTPUT
88  // ai ...
89  // bi ...
90  // ei ...
91  // Le ...
92 
93  if(radius<=0.0) //avoid troubles in case radius is 0 or negative
94  return false;
95 
96  ai = radius * std::pow(aspectRatio*aspectRatio,0.33333333333333333333333);
97  bi = ai / aspectRatio;
98  ei = std::sqrt(
99  1.0
100  - 1.0 / (aspectRatio*aspectRatio)
101  );
102  Le = std::log(
103  (1.0+ei)
104  /(1.0-ei)
105  );
106 
107  return true;
108 }
109 
110 //--------------------------------------------------------------------
111 // Compute the major, minor axis and eccentricity parameters of a prolate spheroid
112 inline double Pi()
113 {
114  return 3.1415926535897932384626433832795;
115 }
116 
117 
118 //--------------------------------------------------------------------
119 // Compute the major, minor axis and eccentricity parameters of a prolate spheroid
120 inline bool spheroidGeometry2(double radius, double aspectRatio, //inputs
121  double& ai, double& bi, //outputs
122  double& XAe, double& YAe, //outputs
123  double& XCe, double& YCe, //outputs
124  double& YHe
125  ) //
126 {
127  //INPUT
128  // radius ...volume-equivalent radius of the spheroid
129  // aspectRatio ...major/minor aspect ratio
130 
131  //OUTPUT
132  // XAe ...Eccentricity dependet parameter
133  // YAe ...Eccentricity dependet parameter
134  // XCe ...Eccentricity dependet parameter
135  // XCe ...Eccentricity dependet parameter
136  // YCe ...Eccentricity dependet parameter
137  // YHe ...Eccentricity dependet parameter
138 
139 
140  double ei(0.0), Le(0.0);
141  bool result =
142  spheroidGeometry(radius, aspectRatio, //inputs
143  ai, bi, //outputs
144  ei, Le //outputs
145  );
146  if(!result)
147  return false;
148 
149  XAe= 2.6666666666666666666666667
150  *ei*ei*ei
151  /(-2.0*ei+(1.0+ei*ei)*Le);
152  YAe= 5.333333333333333333333333333
153  *ei*ei*ei
154  /(2.0*ei+(3*ei*ei-1.0)*Le);
155  XCe= 1.333333333333333333333333333
156  *ei*ei*ei
157  *(1.0-ei*ei)
158  /(2.0*ei-(1.0-ei*ei)*Le);
159  YCe= 1.3333333333333333333333
160  *ei*ei*ei
161  *(2.0-ei*ei)
162  /(-2.0*ei+(1.0+ei*ei)*Le);
163  YHe= 1.3333333333333333333333
164  *ei*ei*ei*ei*ei
165  /(-2.0*ei+(1.0+ei*ei)*Le);
166 
167  return true;
168 
169 }
170 
171 //--------------------------------------------------------------------
172 // zeroize a 3x3x3 tensor
173 inline void zeroize333(double tensor[3][3][3] )
174 {
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;
179 }
180 
181 //--------------------------------------------------------------------
182 // zeroize a 3x3 tensor
183 inline void zeroize33(double tensor[3][3] )
184 {
185  for(int iX=0; iX<3; iX++)
186  for(int iY=0; iY<3; iY++)
187  tensor[iX][iY] = 0.0;
188 }
189 
190 //--------------------------------------------------------------------
191 // multiply a 3x3x3 tensor with a scalar
192 inline void multiply333(double scalar, double tensor[3][3][3] )
193 {
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;
198 }
199 //--------------------------------------------------------------------
200 // Compute a dot and dyadic product of with a vector
201 inline void permutationTensor(double tensor[3][3][3] )
202 {
203  zeroize333(tensor);
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;
210 }
211 
212 
213 
214 //--------------------------------------------------------------------
215 // Compute a dot product of the permutation tensor and
216 // then a dyadic product of with a vector
217 inline bool permutationDotDyadic(
218  double vector[3],
219  double tensor[3][3][3]
220  )
221 {
222  //Generate permutation tensor
223  double permutationT[3][3][3];
224  permutationTensor(permutationT);
225 
226  //Step 1: compute dot prodcut of permutation tensor
227  double permutationDotProd[3][3];
228  zeroize33(permutationDotProd);
229 
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]
234  * vector[iZ];
235 
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]
240  * vector[iZ];
241  return true;
242 
243 }
244 
245 //--------------------------------------------------------------------
246 // Compute a dot and dyadic product of with a vector
247 inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3],
248  double tensor33[3][3],
249  double result[3]
250  )
251 {
252  result[0]=0.0;result[1]=0.0;result[2]=0.0;
253 
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];
258 
259  return true;
260 }
261 
262 
263 }; //end of namespace
264 
265 #endif