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