muMECH  1.0
esuf_Sphere.cpp
Go to the documentation of this file.
1 //********************************************************************************************************
2 // code: ### ### ##### #### ## ##
3 // ## ## ######## ## ## ## ## ##
4 // ## ## ## ## ## ### ## ######
5 // ## ## ## ## ## ## ## ## ##
6 // ##### ## ## ##### #### ## ##
7 // ##
8 //
9 // name: eshelbySoluUniformFieldSphere.cpp
10 // author(s): Jan Novak
11 // language: C, C++
12 // license: This program is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Lesser General Public License as published by
14 // the Free Software Foundation; either version 2 of the License, or
15 // (at your option) any later version.
16 //
17 // This program is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 // GNU Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public License
23 // along with this program; if not, write to the Free Software
24 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
25 //********************************************************************************************************
26 
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30 
31 #include "macros.h"
32 #include "elementaryFunctions.h"
33 #include "esuf_Sphere.h"
34 #include "esei_Sphere.h"
35 
36 
37 namespace mumech {
38 
39 //miscelaneous macros
40 #define X1 x[0]
41 #define X2 x[1]
42 #define X3 x[2]
43 
44 
45 //********************************************************************************************************
46 // PUBLIC FUNCTIONS
47 //********************************************************************************************************
48 
49 
50 //********************************************************************************************************
51 // description: Function gives a component of the Eshelby tensor
52 // last edit: 26. 08. 2009
53 //-------------------------------------------------------------------------------------------------------
54 // note: Solution adopted from Toshio Mura's book (1982)
55 // @flag - Eshelby tensor component S_ijkl
56 // @nu - Poisson's ratio
57 // @a - radius of a sphere
58 //********************************************************************************************************
59 
61  double nu,
62  double a )
63 {
64  //UNREFERENCED_PARAMETER(a);
65 
66  double S_ijkl = 0.;
67 
68  if( flag == _S1111_ || flag == _S2222_ || flag == _S3333_ ){
69  S_ijkl = ( 7. - 5. * nu ) / ( 15. * ( 1. - nu ) );
70  }
71  else if ( flag == _S1122_ || flag == _S2233_ || flag == _S3311_ ||
72  flag == _S1133_ || flag == _S2211_ || flag == _S3322_ ){
73  S_ijkl = ( 5. * nu - 1. ) / ( 15. * ( 1. - nu ) );
74  }
75  else if( flag == _S1212_ || flag == _S2323_ || flag == _S1313_ ){
76  S_ijkl = ( 4. - 5. * nu ) / ( 15. * ( 1. - nu ) );
77  }
78  else{
79  S_ijkl = 0.;
80  }
81 
82  return S_ijkl;
83 }//end of function: giveEshelbyTensCompUniformFieldSphere( )
84 
85 //********************************************************************************************************
86 // description: Function initialize an 1D array by Eshelby tensor components
87 // last edit: 07. 06. 2010
88 //-------------------------------------------------------------------------------------------------------
89 // note: Solution adopted from Toshio Mura's book (1982)
90 // @eshTens - pointer to array of Eshelby tensor components S_ijkl
91 // @a - radius of a sphere
92 //********************************************************************************************************
93 void eshelbySoluUniformFieldSphere :: eshelbyTensUniformField (double eshTens[12], const double sort_a[3],
94  const double eSInt[13] )
95 {
96  double a = sort_a[0];
97  /* 1st row */
98  eshTens[0] = this->giveEshelbyTensCompUniformField( _S1111_, nu, a );//S_1111
99  eshTens[1] = this->giveEshelbyTensCompUniformField( _S1122_, nu, a );//S_1122
100  eshTens[2] = eshTens[1];//S_1133
101  /* 2nd row */
102  eshTens[3] = eshTens[1];//S_2211
103  eshTens[4] = eshTens[0];//S_2222
104  eshTens[5] = eshTens[1];//S_2233
105  /* 3rd row */
106  eshTens[6] = eshTens[1];//S_3311
107  eshTens[7] = eshTens[1];//S_3322
108  eshTens[8] = eshTens[0];//S_3333
109  /* !!! FEEP notation !!! e_ij = {e_11, e_22, e_33, e_12, e_32, e_31}^T */
110  /* 4th row - 2* ... because of the Voigth notation! */
111  eshTens[9] = 2. * this->giveEshelbyTensCompUniformField( _S1212_, nu, a ); //2 * S_1212
112  /* 5th row - 2* ... because of the Voigth notation! */
113  eshTens[10] = eshTens[9];//2 * S_2323
114  /* 6th row - 2* ... because of the Voigth notation! */
115  eshTens[11] = eshTens[9];//2 * S_1313
116 
117  return ;
118 }//end of function: giveEshelbyTensor( )
119 
120 
123 {
124  /* !!! FEEP notation !!! e_ij = {e_11, e_22, e_33, 2e_12, 2e_32, 2e_31}^T */
125 
126  double aux = SQR( S[_11_] ) + S[_11_] * S[_12_] - 2. * SQR( S[_12_] );
127  double kii = ( S[_11_] + S[_12_] ) / aux;
128  double kij = - S[_12_] / aux;
129  double gii = 1. / S[_44_];
130 
131  SInv[_11_] = kii; SInv[_12_] = kij; SInv[_13_] = kij;
132  SInv[_21_] = kij; SInv[_22_] = kii; SInv[_23_] = kij;
133  SInv[_31_] = kij; SInv[_32_] = kij; SInv[_33_] = kii;
134 
135  SInv[_44_] = gii;
136  SInv[_55_] = gii;
137  SInv[_66_] = gii;
138 }
139 
140 //********************************************************************************************************
141 // description: Function gives D_ijkl(x) generalized Eshelby's tensor.
142 // Function requiers "constructor2"
143 // last edit: 13. 07. 2010
144 //-------------------------------------------------------------------------------------------------------
145 // note: @D[36] - Perturbation (generalized Eshelby's) tensor to be evaluated and saved
146 // @S[12] - Eshelby's tensor for internal fields (uniform fields)
147 // @J[13] - Ferers-Dysons' integrals elliptic integrals
148 // @dJi[9] - first derivatives of Ferers-Dysons' integrals elliptic integrals Ii
149 // @dJij[27] - first derivatives of Ferers-Dysons' integrals elliptic integrals Iij
150 // @ddJi[27] - second derivatives of Ferers-Dysons' integrals elliptic integrals Ii
151 // @ddJij[81] - second derivatives of Ferers-Dysons' integrals elliptic integrals Ii
152 // @aa - squared power of a_i = a = a_1 = a_2 = a_3 = sphere_radius
153 // @x[3] - coordinates of a point
154 //********************************************************************************************************
155 void eshelbySoluUniformFieldSphere::giveDijkl( double D[36], const double S[12], const double J[13], const double dJi[9],
156  const double dJij[27], const double ddJi[27], const double ddJij[81],
157  const double sort_a[3], const double x[3])
158 {
159  double aa = SQR(sort_a[0]);
160  // if ( pos == PPF_INT_POINT ) {
161  // CleanVector( D, 36 );
162  // //1st row
163  // D[_1111_] = S[__1111_]; D[_1122_] = S[__1122_]; D[_1133_] = S[__1133_];
164  // //2nd row
165  // D[_2211_] = S[__2211_]; D[_2222_] = S[__2222_]; D[_2233_] = S[__2233_];
166  // //3rd row
167  // D[_3311_] = S[__3311_]; D[_3322_] = S[__3322_]; D[_3333_] = S[__3333_];
168  // //4th row
169  // D[_1212_] = S[__1212_];
170  // //5th row
171  // D[_2323_] = S[__2323_];
172  // //6th row
173  // D[_1313_] = S[__1313_];
174  // }
175  // else if ( pos == PPF_EXT_POINT ) {
176 
177  double auxS[36];//auxiliary array to save Eshelby tensor
178  //Position dependent Eshelby tensor:
179  this->giveSijkl( auxS, J, sort_a, nu , false);
180 
181  /*1st row*/
182  D[_1111_] = this->multTRN*(-(X1*(this->_1Plus2nu*J1_1 + J1_11*X1 - aa*(5*J11_1 + X1*J11_11)))) +
183  auxS[_1111_];
184 
185  D[_1122_] = this->multTRN*(this->_2nu*J1_1*X1 + SQR(X1)*(-J1_22 + aa*J11_22) + X2*(-J2_2 + aa*J21_2)) +
186  auxS[_1122_];
187 
188  D[_1133_] = this->multTRN*(this->_2nu*J1_1*X1 + SQR(X1)*(-J1_33 + aa*J11_33) + X3*(-J3_3 + aa*J31_3)) +
189  auxS[_1133_];
190 
191  D[_1123_] = this->multTRN*(SQR(X1)*(-J1_32 + aa*J11_32) + X2*(-J2_3 + aa*J21_3));
192 
193  D[_1113_] = this->multTRN*(X1*(-3*J1_3 - J1_31*X1 + aa*(3*J11_3 + X1*J11_31)) +
194  2*this->_1MinNu*J3_1*X3);
195 
196  D[_1112_] = this->multTRN*(X1*(-3*J1_2 - J1_21*X1 + aa*(3*J11_2 + X1*J11_21)) +
197  2*this->_1MinNu*J2_1*X2);
198 
199  /*2nd row*/
200  D[_2211_] = this->multTRN*(this->_2nu*J2_2*X2 + X1*(-J1_1 + aa*J12_1) + SQR(X2)*(-J2_11 + aa*J22_11)) +
201  auxS[_2211_];
202 
203  D[_2222_] = this->multTRN*(-(X2*(this->_1Plus2nu*J2_2 + J2_22*X2 - aa*(5*J22_2 + X2*J22_22)))) +
204  auxS[_2222_];
205 
206  D[_2233_] = this->multTRN*(this->_2nu*J2_2*X2 + SQR(X2)*(-J2_33 + aa*J22_33) +
207  X3*(-J3_3 + aa*J32_3)) + auxS[_2233_];
208 
209  D[_2223_] = this->multTRN*(X2*(-3*J2_3 - J2_32*X2 + aa*(3*J22_3 + X2*J22_32)) +
210  2*this->_1MinNu*J3_2*X3);
211 
212  D[_2213_] = this->multTRN*(X1*(-J1_3 + aa*J12_3) + SQR(X2)*(-J2_31 + aa*J22_31));
213 
214  D[_2212_] = this->multTRN*(this->_1Min2nu*J1_2*X1 + aa*X1*J12_2 +
215  X2*(-2*J2_1 - J2_21*X2 + aa*(2*J22_1 + X2*J22_21)));
216 
217  /*3rd row*/
218  D[_3311_] = this->multTRN*(this->_2nu*J3_3*X3 + X1*(-J1_1 + aa*J13_1) +
219  SQR(X3)*(-J3_11 + aa*J33_11)) + auxS[_3311_];
220 
221  D[_3322_] = this->multTRN*(this->_2nu*J3_3*X3 + X2*(-J2_2 + aa*J23_2) +
222  SQR(X3)*(-J3_22 + aa*J33_22)) + auxS[_3322_];
223 
224  D[_3333_] = this->multTRN*(-(X3*(this->_1Plus2nu*J3_3 + J3_33*X3 - aa*(5*J33_3 + X3*J33_33)))) +
225  auxS[_3333_];
226 
227  D[_3323_] = this->multTRN*(this->_1Min2nu*J2_3*X2 + aa*X2*J23_3 +
228  X3*(-2*J3_2 - J3_32*X3 + aa*(2*J33_2 + X3*J33_32)));
229 
230  D[_3313_] = this->multTRN*(this->_1Min2nu*J1_3*X1 + aa*X1*J13_3 +
231  X3*(-2*J3_1 - J3_31*X3 + aa*(2*J33_1 + X3*J33_31)));
232 
233  D[_3312_] = this->multTRN*(X1*(-J1_2 + aa*J13_2) + SQR(X3)*(-J3_21 + aa*J33_21));
234 
235  /*4th row*/
236  D[_1211_] = this->multTRN*(2*J1_2*X1 + X2*(-2*J2_1 - J2_11*X1 + aa*(2*J12_1 + X1*J12_11)));
237 
238  D[_1222_] = this->multTRN*(this->_2nu*J1_2*X1 + X1*(-2*J2_2 - J2_22*X2 + aa*(2*J12_2 + X2*J12_22)) +
239  2*this->_1MinNu*J2_1*X2);
240 
241  D[_1233_] = this->multTRN*(X1*(this->_2nu*J1_2 + X2*(-J2_33 + aa*J12_33)));
242 
243  D[_1212_] = this->multTRN*(this->_1MinNu*J1_1*X1 - J2_1*X1 - nu*J2_2*X2 - J2_21*X1*X2 + aa*X1*J12_1 +
244  aa*X2*J12_2 + aa*X1*X2*J12_21) + auxS[_1212_];
245 
246  D[_1223_] = this->multTRN*(X1*(-J2_3 - J2_32*X2 + aa*(J12_3 + X2*J12_32)) +
247  this->_1MinNu*J3_1*X3);
248 
249  D[_1213_] = this->multTRN*(X2*(-J2_3 - J2_31*X1 + aa*(J12_3 + X1*J12_31)) +
250  this->_1MinNu*J3_2*X3);
251 
252  /*5th row*/
253  D[_2311_] = this->multTRN*(X2*(this->_2nu*J2_3 + X3*(-J3_11 + aa*J23_11)));
254 
255  D[_2322_] = this->multTRN*(2*J2_3*X2 + X3*(-2*J3_2 - J3_22*X2 + aa*(2*J23_2 + X2*J23_22)));
256 
257  D[_2333_] = this->multTRN*(this->_2nu*J2_3*X2 + X2*(-2*J3_3 - J3_33*X3 + aa*(2*J23_3 + X3*J23_33)) +
258  2*this->_1MinNu*J3_2*X3);
259 
260  D[_2312_] = this->multTRN*(this->_1MinNu*J1_3*X1 + X3*(-J3_1 - J3_21*X2 + aa*(J23_1 + X2*J23_21)));
261 
262  D[_2323_] = this->multTRN*(this->_1MinNu*J2_2*X2 - J3_2*X2 - nu*J3_3*X3 - J3_32*X2*X3 + aa*X2*J23_2 +
263  aa*X3*J23_3 + aa*X2*X3*J23_32) + auxS[_2323_];
264 
265  D[_2313_] = this->multTRN*(this->_1MinNu*J1_2*X1 + X2*(-J3_1 - J3_31*X3 + aa*(J23_1 + X3*J23_31)));
266 
267  /*6th row*/
268  D[_1311_] = this->multTRN*(2*J1_3*X1 + X3*(-2*J3_1 - J3_11*X1 + aa*(2*J13_1 + X1*J13_11)));
269 
270  D[_1322_] = this->multTRN*(X1*(this->_2nu*J1_3 + X3* (-J3_22 + aa*J13_22)));
271 
272  D[_1333_] = this->multTRN*(this->_2nu*J1_3*X1 + X1*(-2*J3_3 - J3_33*X3 + aa*(2*J13_3 + X3*J13_33)) +
273  2*this->_1MinNu*J3_1*X3);
274 
275  D[_1312_] = this->multTRN*(this->_1MinNu*J2_3*X2 + X3*(-J3_2 - J3_21*X1 + aa*(J13_2 + X1*J13_21)));
276 
277  D[_1323_] = this->multTRN*(X1*(-J3_2 - J3_32*X3 + aa*(J13_2 + X3*J13_32)) + this->_1MinNu*J2_1*X2);
278 
279  D[_1313_] = this->multTRN*(this->_1MinNu*J1_1*X1 - J3_1*X1 - nu*J3_3*X3 - J3_31*X1*X3 + aa*X1*J13_1 +
280  aa*X3*J13_3 + aa*X1*X3*J13_31) + auxS[_1313_];
281 
282  // TomTag:
283  // No Mandel notation.
284  // eps_star_ij = D_ijkl * eps_tau_kl
285  // e.g. eps_11 = ... + D_1112 eps_tau_12 + ... + D_1121 eps_tau_21
286  // Therefore thre right half of matrix D is multiplied by two.
287  D[_1112_] *= 2.0;
288  D[_1123_] *= 2.0;
289  D[_1113_] *= 2.0;
291  D[_2212_] *= 2.0;
292  D[_2223_] *= 2.0;
293  D[_2213_] *= 2.0;
295  D[_3312_] *= 2.0;
296  D[_3323_] *= 2.0;
297  D[_3313_] *= 2.0;
299  D[_1212_] *= 2.0;
300  D[_1223_] *= 2.0;
301  D[_1213_] *= 2.0;
303  D[_2312_] *= 2.0;
304  D[_2323_] *= 2.0;
305  D[_2313_] *= 2.0;
307  D[_1312_] *= 2.0;
308  D[_1323_] *= 2.0;
309  D[_1313_] *= 2.0;
310  // }
311  // else _errorr( "giveDijkl: invalid position of a given point");
312 
313  return ;
314 }//end of function: giveDijkl( )
315 
316 
317 //********************************************************************************************************
318 // PROTECTED FUNCTIONS
319 //********************************************************************************************************
320 
321 //********************************************************************************************************
322 // description: Function gives the displacement perturbation tensor of internal fields.
323 // Function requiers "constructor2"
324 // last edit: 15. 07. 2010
325 //-------------------------------------------------------------------------------------------------------
326 // note: @Lint[18] - Perturbation (Eshelby's-like) displacement tensor to be evaluated and saved
327 // @J[13] - Ferers-Dysons' integrals elliptic integrals
328 // @aa - squared power of a_i = a = a_1 = a_2 = a_3 = sphere_radius
329 // @x[3] - coordinates of a point
330 //********************************************************************************************************
331 void eshelbySoluUniformFieldSphere::giveLijkINT( double Lint[18], const double J[13],const double a[3], const double x[3] )
332 {
333  double aa = SQR(a[0]);
334 
335  /* 1st row */
336  Lint[_111_] = this->multTRN * ( X1 * ( this->_1Min2nu * J1 + 3. * aa * J11 ) );
337  Lint[_122_] = this->multTRN * ( X1 * ( this->_2nu * J1 - J2 + aa * J12 ) );
338  Lint[_133_] = this->multTRN * ( X1 * ( this->_2nu * J1 - J3 + aa * J13 ) );
339  Lint[_112_] = this->multTRN * ( X2 * ( this->_1Min2nu * J2 + aa * J12 ) );
340  Lint[_123_] = 0.;
341  Lint[_113_] = this->multTRN * ( X3 * ( this->_1Min2nu * J3 + aa * J13 ) );
342  /* 2nd row */
343  Lint[_211_] = this->multTRN * ( X2 * ( this->_2nu * J2 - J1 + aa * J21 ) );
344  Lint[_222_] = this->multTRN * ( X2 * ( this->_1Min2nu * J2 + 3. * aa * J22 ) );
345  Lint[_233_] = this->multTRN * ( X2 * ( this->_2nu * J2 - J3 + aa * J23 ) );
346  Lint[_212_] = this->multTRN * ( X1 * ( this->_1Min2nu * J1 + aa * J21 ) );
347  Lint[_223_] = this->multTRN * ( X3 * ( this->_1Min2nu * J3 + aa * J23 ) );
348  Lint[_213_] = 0.;
349  /* 3rd row */
350  Lint[_311_] = this->multTRN * ( X3 * ( this->_2nu * J3 - J1 + aa * J31 ) );
351  Lint[_322_] = this->multTRN * ( X3 * ( this->_2nu * J3 - J2 + aa * J32 ) );
352  Lint[_333_] = this->multTRN * ( X3 * ( this->_1Min2nu * J3 + 3. * aa * J33 ) );
353  Lint[_312_] = 0.;
354  Lint[_323_] = this->multTRN * ( X2 * ( this->_1Min2nu * J2 + aa * J32 ) );
355  Lint[_313_] = this->multTRN * ( X1 * ( this->_1Min2nu * J1 + aa * J31 ) );
356 
357  //Mandel notation adjusting
358  /* 1nd row */
359  Lint[_112_] *= 2;
360  Lint[_123_] *= 2;
361  Lint[_113_] *= 2;
362  /* 2nd row */
363  Lint[_212_] *= 2;
364  Lint[_223_] *= 2;
365  Lint[_213_] *= 2;
366  /* 3rd row */
367  Lint[_312_] *= 2;
368  Lint[_323_] *= 2;
369  Lint[_313_] *= 2;
370 
371  return ;
372 }//end of function: giveLijkINT( )
373 
374 //********************************************************************************************************
375 // description: Function gives the displacement perturbation tensor of external fields.
376 // Function requiers "constructor2"
377 // last edit: 15. 07. 2010
378 //-------------------------------------------------------------------------------------------------------
379 // note: @Lext[18] - Perturbation (generalized Eshelby's) displacement tensor to be evaluated
380 // and saved
381 // @Lint[18] - Perturbation (Eshelby's-like) displacement tensor
382 // @dJi[9] - first derivatives of Ferers-Dysons' integrals elliptic integrals Ii
383 // @dJij[27] - first derivatives of Ferers-Dysons' integrals elliptic integrals Iij
384 // @aa - squared power of a_i = a = a_1 = a_2 = a_3 = sphere_radius
385 // @x[3] - coordinates of a point
386 // @nu - Poisson's ratio of the matrix material
387 //********************************************************************************************************
388 void eshelbySoluUniformFieldSphere::giveLijkEXT( double Lext[18], const double Lint[18], const double dJi[9],
389  const double dJij[27], const double sort_a[3], const double x[3])
390 {
391  double aa = SQR(sort_a[0]);
392 
393  double xx1 = SQR( X1 ), xx2 = SQR( X2 ), xx3 = SQR( X3 ), x1x2 = X1 * X2, x2x3 = X2 * X3,
394  x1x3 = X1 * X3;
395 
396  for( int i = 0; i < 18; i++ ){
397  Lext[i] = Lint[i];
398  }
399 
400  /* 1st row */
401  Lext[_111_] += this->multTRN * xx1 * ( -J1_1 + aa * J11_1 );
402  Lext[_122_] += this->multTRN * x1x2 * ( -J2_2 + aa * J12_2 );
403  Lext[_133_] += this->multTRN * x1x3 * ( -J3_3 + aa * J13_3 );
404  Lext[_112_] += 2 * this->multTRN * xx1 * ( -J1_2 + aa * J11_2 );
405  Lext[_123_] += 2 * this->multTRN * x1x2 * ( -J2_3 + aa * J12_3 );
406  Lext[_113_] += 2 * this->multTRN * xx1 * ( -J1_3 + aa * J11_3 );
407  /* 2nd row */
408  Lext[_211_] += this->multTRN * x1x2 * ( -J1_1 + aa * J21_1 );
409  Lext[_222_] += this->multTRN * xx2 * ( -J2_2 + aa * J22_2 );
410  Lext[_233_] += this->multTRN * x2x3 * ( -J3_3 + aa * J23_3 );
411  Lext[_212_] += 2 * this->multTRN * x1x2 * ( -J1_2 + aa * J21_2 );
412  Lext[_223_] += 2 * this->multTRN * xx2 * ( -J2_3 + aa * J22_3 );
413  Lext[_213_] += 2 * this->multTRN * x1x2 * ( -J1_3 + aa * J21_3 );
414  /* 3rd row */
415  Lext[_311_] += this->multTRN * x1x3 * ( -J1_1 + aa * J31_1 );
416  Lext[_322_] += this->multTRN * x2x3 * ( -J2_2 + aa * J32_2 );
417  Lext[_333_] += this->multTRN * xx3 * ( -J3_3 + aa * J33_3 );
418  Lext[_312_] += 2 * this->multTRN * x1x3 * ( -J1_2 + aa * J31_2 );
419  Lext[_323_] += 2 * this->multTRN * x2x3 * ( -J2_3 + aa * J32_3 );
420  Lext[_313_] += 2 * this->multTRN * x1x3 * ( -J1_3 + aa * J31_3 );
421 
422  return ;
423 }//end of function: giveLijkEXT( )
424 
425 //********************************************************************************************************
426 // description: Function gives S_ijkl(x) position dependent Eshelby's tensor. (This is not
427 // the strain perturbation tensor! Just its part ). The returned members of the
428 // S_ijkl matrix are in tensorial (theoretical) notation!
429 // Function requiers "constructor2"
430 // last edit: 15. 06. 2010
431 //-------------------------------------------------------------------------------------------------------
432 // note: @S[36] - Eshelby position dependent tensor
433 // @J[13] - Ferers-Dysons' integrals elliptic integrals
434 // @aa - squared power of a_i = a = a_1 = a_2 = a_3 = sphere_radius
435 // @nu - Poisson's ratio of the matrix material
436 //********************************************************************************************************
437 void eshelbySoluUniformFieldSphere::giveSijkl( double S[36], const double J[13], const double sort_a[3], double nu , bool newFormulation)
438 {
439  double aa = SQR(sort_a[0]);
440 
441  double mult = this->multTRN;
442  // double mult = 1.;
443 
444  //1st row
445  S[_1111_] = mult * ( this->_1Min2nu * J1 + 3. * aa * J11 );
446  S[_1122_] = mult * ( this->_2nu * J1 - J2 + aa * J21 );
447  S[_1133_] = mult * ( this->_2nu * J1 - J3 + aa * J31 );
448  S[_1123_] = S[_1113_] = S[_1112_] = 0.;
449  //2nd row
450  S[_2211_] = mult * ( this->_2nu * J2 - J1 + aa * J12 );
451  S[_2222_] = mult * ( this->_1Min2nu * J2 + 3. * aa * J22 );
452  S[_2233_] = mult * ( this->_2nu * J2 - J3 + aa * J32 );
453  S[_2223_] = S[_2213_] = S[_2212_] = 0.;
454  //3rd row
455  S[_3311_] = mult * ( this->_2nu * J3 - J1 + aa * J13 );
456  S[_3322_] = mult * ( this->_2nu * J3 - J2 + aa * J23 );
457  S[_3333_] = mult * ( this->_1Min2nu * J3 + 3. * aa * J33 );
458  S[_3323_] = S[_3313_] = S[_3312_] = 0.;
459  //4th row
460  S[_1211_] = S[_1222_] = S[_1233_] = 0.;
461  S[_1212_] = mult * ( this->_1MinNu * J1 - nu * J2 + aa * J12 );
462  S[_1223_] = S[_1213_] = 0.;
463  //5th row
464  S[_2311_] = S[_2322_] = S[_2333_] = S[_2312_] = 0.;
465  S[_2323_] = mult * ( this->_1MinNu * J2 - nu * J3 + aa * J23 );
466  S[_2313_] = 0.;
467  //6th row
468  S[_1311_] = S[_1322_] = S[_1333_] = S[_1312_] = S[_1323_] = 0.;
469  S[_1313_] = mult * ( this->_1MinNu * J1 - nu * J3 + aa * J13 );
470 
471  return ;
472 }//end of function: giveSijkl( )
473 
474 } // end of namespace mumech
475 
476 /*end of file*/
#define J11_3
Definition: types.h:490
#define J1_1
Definition: types.h:451
#define J31_2
Definition: types.h:485
#define J2_3
Definition: types.h:457
#define J11_32
Definition: types.h:553
#define J12_2
Definition: types.h:478
#define _2223_
Definition: types.h:326
#define _2311_
Definition: types.h:343
#define J12_33
Definition: types.h:567
#define _1323_
Definition: types.h:354
#define J33_11
Definition: types.h:648
#define J13_2
Definition: types.h:479
#define _11_
Definition: types.h:253
#define J2_33
Definition: types.h:528
#define J12_32
Definition: types.h:566
#define J23_33
Definition: types.h:619
#define _1222_
Definition: types.h:337
#define J1_2
Definition: types.h:452
#define J3_3
Definition: types.h:461
#define J12_11
Definition: types.h:557
#define _212_
Definition: types.h:371
#define _3322_
Definition: types.h:330
#define J3_2
Definition: types.h:460
#define J33
Definition: types.h:444
#define J13_22
Definition: types.h:575
#define J2_31
Definition: types.h:526
#define _233_
Definition: types.h:370
#define J11_22
Definition: types.h:549
#define J21_1
Definition: types.h:468
#define J12_22
Definition: types.h:562
#define J31
Definition: types.h:442
#define J1_31
Definition: types.h:513
#define J12_3
Definition: types.h:491
#define J13_11
Definition: types.h:570
#define _1123_
Definition: types.h:319
#define _2212_
Definition: types.h:325
#define _44_
Definition: types.h:265
#define J2_11
Definition: types.h:518
#define J22_21
Definition: types.h:600
#define J33_22
Definition: types.h:653
double nu
nu of matrix
Definition: esuf.h:51
#define _21_
Definition: types.h:257
#define J2_21
Definition: types.h:522
#define _3313_
Definition: types.h:334
#define J22_32
Definition: types.h:605
#define J3_32
Definition: types.h:540
#define _2213_
Definition: types.h:327
#define J11_2
Definition: types.h:477
#define J23_22
Definition: types.h:614
#define _122_
Definition: types.h:362
#define J13_33
Definition: types.h:580
#define _1322_
Definition: types.h:351
#define _312_
Definition: types.h:378
#define _1212_
Definition: types.h:339
#define J11_33
Definition: types.h:554
#define _2322_
Definition: types.h:344
#define _113_
Definition: types.h:366
#define J3
Definition: types.h:432
#define J33_2
Definition: types.h:487
void giveDijkl(double D[12], const double S[], const double J[], const double dJi[], const double dJij[], const double ddJi[], const double ddJij[], const double sort_a[3], const double x[])
#define _1223_
Definition: types.h:340
The header file of usefull macros.
#define J3_31
Definition: types.h:539
Collection of the functions of basic manipulations, some of them can be used instead of their counter...
#define J2
Definition: types.h:431
#define J3_11
Definition: types.h:531
#define J1_32
Definition: types.h:514
#define _133_
Definition: types.h:363
#define J23_2
Definition: types.h:483
#define J2_1
Definition: types.h:455
#define _1113_
Definition: types.h:320
#define J33_33
Definition: types.h:658
#define _123_
Definition: types.h:365
#define _2333_
Definition: types.h:345
#define _12_
Definition: types.h:254
#define _2233_
Definition: types.h:324
#define _1213_
Definition: types.h:341
#define _2222_
Definition: types.h:323
#define _31_
Definition: types.h:261
#define J31_3
Definition: types.h:498
#define J31_1
Definition: types.h:472
#define X3
Definition: esuf_Sphere.cpp:42
#define X2
Definition: esuf_Sphere.cpp:41
#define J32
Definition: types.h:443
#define J12_21
Definition: types.h:561
#define J13
Definition: types.h:436
#define J23_21
Definition: types.h:613
#define J3_33
Definition: types.h:541
#define _1111_
Definition: types.h:315
void eshelbyTensUniformField(double eshTens[12], const double sort_a[3], const double eSInt[13])
Definition: esuf_Sphere.cpp:93
#define J11
Definition: types.h:434
#define J13_21
Definition: types.h:574
#define X1
Definition: esuf_Sphere.cpp:40
#define J12_1
Definition: types.h:465
#define _1122_
Definition: types.h:316
#define J12_31
Definition: types.h:565
#define _112_
Definition: types.h:364
#define J21_3
Definition: types.h:494
#define J1_11
Definition: types.h:505
#define J22_22
Definition: types.h:601
#define J33_1
Definition: types.h:474
#define J23_3
Definition: types.h:496
#define J13_32
Definition: types.h:579
#define _222_
Definition: types.h:369
#define J22_3
Definition: types.h:495
EshelbyTensComponent
Definition: types.h:221
#define _55_
Definition: types.h:267
#define SQR(a)
Definition: macros.h:97
#define J22
Definition: types.h:439
#define J33_21
Definition: types.h:652
#define _1211_
Definition: types.h:336
#define _3323_
Definition: types.h:333
#define J21_2
Definition: types.h:481
#define _2312_
Definition: types.h:346
#define J22_33
Definition: types.h:606
#define J22_31
Definition: types.h:604
#define _2211_
Definition: types.h:322
#define J22_1
Definition: types.h:469
#define J13_3
Definition: types.h:492
#define J1
Definition: types.h:430
#define _32_
Definition: types.h:262
#define J23_32
Definition: types.h:618
#define _1333_
Definition: types.h:352
#define _13_
Definition: types.h:255
void giveSijkl(double S[36], const double J[], const double sort_a[], double nu, bool newFormulation)
#define J33_32
Definition: types.h:657
#define J32_3
Definition: types.h:499
#define J11_31
Definition: types.h:552
#define _1312_
Definition: types.h:353
#define J11_1
Definition: types.h:464
#define _322_
Definition: types.h:376
#define _3311_
Definition: types.h:329
#define J21
Definition: types.h:438
#define J23_31
Definition: types.h:617
#define _23_
Definition: types.h:259
#define _211_
Definition: types.h:368
#define J22_11
Definition: types.h:596
void giveLijkEXT(double Lext[18], const double Lint[18], const double dJi[9], const double dJij[27], const double sort_a[3], const double x[3])
#define J2_2
Definition: types.h:456
#define J33_31
Definition: types.h:656
#define J1_33
Definition: types.h:515
#define J3_21
Definition: types.h:535
#define J2_22
Definition: types.h:523
#define J32_2
Definition: types.h:486
#define J12
Definition: types.h:435
#define _213_
Definition: types.h:373
void giveLijkINT(double Lint[18], const double J[13], const double a[3], const double x[3])
double multTRN
1./mult
Definition: esuf.h:57
#define _2313_
Definition: types.h:348
#define J11_21
Definition: types.h:548
#define _3333_
Definition: types.h:331
#define J1_21
Definition: types.h:509
#define J13_1
Definition: types.h:466
#define J2_32
Definition: types.h:527
Class eshelbySoluUniformFieldSphere; return the Eshelby solution of a spherical inclusion loaded by t...
#define J3_22
Definition: types.h:536
Class eshelbySoluEllipticIntegralsSphere.
#define J13_31
Definition: types.h:578
#define _1311_
Definition: types.h:350
#define _1133_
Definition: types.h:317
#define J3_1
Definition: types.h:459
#define J22_2
Definition: types.h:482
#define _111_
Definition: types.h:361
#define J1_3
Definition: types.h:453
#define _323_
Definition: types.h:379
#define _66_
Definition: types.h:269
#define J23_1
Definition: types.h:470
#define _311_
Definition: types.h:375
#define _22_
Definition: types.h:258
#define _313_
Definition: types.h:380
#define _33_
Definition: types.h:263
#define J33_3
Definition: types.h:500
#define _333_
Definition: types.h:377
#define _1313_
Definition: types.h:355
#define J11_11
Definition: types.h:544
#define J23
Definition: types.h:440
void giveEshelbyTensorInverse(double SInv[12], const double S[])
#define J23_11
Definition: types.h:609
#define _1233_
Definition: types.h:338
#define _223_
Definition: types.h:372
#define _2323_
Definition: types.h:347
double giveEshelbyTensCompUniformField(EshelbyTensComponent flag, double nu, double a)
Definition: esuf_Sphere.cpp:60
#define _1112_
Definition: types.h:318
#define _3312_
Definition: types.h:332
#define J1_22
Definition: types.h:510