esuf_ProlateSpheroid.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluUniformFieldProlateSpheroid.cpp
// description:    source file of the function returning the Eshelby solution of an Prolate-spheroidal inclusion
//                 loaded by the uniform strain/stress field
// author(s):      Jan Novak, Tomas Janda
// place of birth: Czech Technical University in Prague, Faculty of Civil engineering
// last edit:      08. 08. 2013
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "macros.h"
#include "elementaryFunctions.h"
#include "esuf_ProlateSpheroid.h"
#include "esei_ProlateSpheroid.h"
#include "legendreIntegrals.h"
#include "types.h"

namespace mumech {

//miscelaneous macros
//-------------------
//coordinates
#define X1    x[0]
#define X2    x[1]
#define X3    x[2]
//elliptic integrals
#define __I( component )           eSInt[component]



//********************************************************************************************************
//                                            PUBLIC FUNCTIONS
//********************************************************************************************************

//********************************************************************************************************
// description: Function gives a component of the Eshelby tensor
// last edit:   08. 08. 2013
//-------------------------------------------------------------------------------------------------------
// note:        Solution adopted from Toshio Mura's book (1982) and Eshelby article (1957)
//              @sort_a[3] - sorted ellipsoidal semiaxes (a1>a2>a3)
//              @eSInt[13] - values of elliptical integrals (potentials)
//                           always for lambda = 0. Even for external fields !!!
//              @nu        - Poisson's ratio
//              @flag      - Eshelby tensor component S_ijkl
//********************************************************************************************************
double eshelbySoluUniformFieldProlateSpheroid::eshelbyTensCompUniformField( const double sort_a[3],
                                               const double eSInt[13], double nu,
                                               EshelbyTensComponent flag )
{
  double a1 = sort_a[0], a2 = sort_a[1], a3 = sort_a[2];
  double S_ijkl = 0.;

  switch( flag ){
    /* 1st row */
  case _S1111_ :
    S_ijkl = MULT * ( 3. * SQR( a1 ) * __I(_Int11_) + ( 1. - 2. * nu ) * __I(_Int1_) );
    break;
  case _S1122_ :
    S_ijkl = MULT * (      SQR( a2 ) * __I(_Int12_) - ( 1. - 2. * nu ) * __I(_Int1_) );
    break;
  case _S1133_ :
    S_ijkl = MULT * (      SQR( a3 ) * __I(_Int13_) - ( 1. - 2. * nu ) * __I(_Int1_) );
    break;
    /* 2nd row */
  case _S2211_ :
    S_ijkl = MULT * (      SQR( a1 ) * __I(_Int21_) - ( 1. - 2. * nu ) * __I(_Int2_) );
    break;
  case _S2222_ :
    S_ijkl = MULT * ( 3. * SQR( a2 ) * __I(_Int22_) + ( 1. - 2. * nu ) * __I(_Int2_) );
    break;
  case _S2233_ :
    S_ijkl = MULT * (      SQR( a3 ) * __I(_Int23_) - ( 1. - 2. * nu ) * __I(_Int2_) );
    break;
    /* 3rd row */
  case _S3311_ :
    S_ijkl = MULT * (      SQR( a1 ) * __I(_Int31_) - ( 1. - 2. * nu ) * __I(_Int3_) );
    break;
  case _S3322_ :
    S_ijkl = MULT * (      SQR( a2 ) * __I(_Int32_) - ( 1. - 2. * nu ) * __I(_Int3_) );
    break;
  case _S3333_ :
    S_ijkl = MULT * ( 3. * SQR( a3 ) * __I(_Int33_) + ( 1. - 2. * nu ) * __I(_Int3_) );
    break;
    /* !!! FEEP notation !!! e_ij = {e_11, e_22, e_33, e_12, e_32, e_31}^T */
    /* 4th row */
  case _S1212_ :
    S_ijkl = .5 * MULT * ( ( 1. - 2. * nu ) * ( __I(_Int1_) + __I(_Int2_) ) +
               ( SQR( a1 ) + SQR( a2 ) ) * __I(_Int12_) );
    break;
    /* 5th row */
  case _S2323_ :
    S_ijkl = .5 * MULT * ( ( 1. - 2. * nu ) * ( __I(_Int2_) + __I(_Int3_) ) +
               ( SQR( a2 ) + SQR( a3 ) ) * __I(_Int23_) );
    break;
    /* 6th row */
  case _S1313_ :
    S_ijkl = .5 * MULT * ( ( 1. - 2. * nu ) * ( __I(_Int1_) + __I(_Int3_) ) +
               ( SQR( a1 ) + SQR( a3 ) ) * __I(_Int13_) );
    break;
  default:
    S_ijkl = 0.;
    break;
  }//end of switch: flag

  return S_ijkl;
}//end of function: eshelbyTensCompUniformFieldProlateSpheroid( )



//********************************************************************************************************
// description: Function gives D_ijkl(x) generalized Eshelby's tensor.
//              Function requiers "constructor2"
// last edit:   08. 08. 2013
//-------------------------------------------------------------------------------------------------------
// note:        @D[36]     - Perturbation (generalized Eshelby's) tensor to be evaluated and saved
//              @S[12]     - Eshelby's tensor for internal fields (uniform fields)
//              @J[13]     - Ferers-Dysons' elliptic integrals
//              @dJi[9]    - first derivatives of Ferers-Dysons' elliptic integrals Ii
//              @dJij[27]  - first derivatives of Ferers-Dysons' elliptic integrals Iij
//              @ddJi[27]  - second derivatives of Ferers-Dysons' elliptic integrals Ii
//              @ddJij[81] - second derivatives of Ferers-Dysons' elliptic integrals Ii
//              @sort_a[3] - sorted semiaxes' dimensions
//              @x[3]      - coordinates of a point

//********************************************************************************************************
void eshelbySoluUniformFieldProlateSpheroid::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])
{
  // if ( pos == PPF_INT_POINT ) {
  //   CleanVector( D, 36 );
  //   //1st row
  //   D[_1111_] = S[__1111_]; D[_1122_] = S[__1122_]; D[_1133_] = S[__1133_];
  //   //2nd row
  //   D[_2211_] = S[__2211_]; D[_2222_] = S[__2222_]; D[_2233_] = S[__2233_];
  //   //3rd row
  //   D[_3311_] = S[__3311_]; D[_3322_] = S[__3322_]; D[_3333_] = S[__3333_];
  //   //4th row
  //   D[_1212_] = S[__1212_];
  //   //5th row
  //   D[_2323_] = S[__2323_];
  //   //6th row
  //   D[_1313_] = S[__1313_];
  // }
  // else if ( pos == PPF_EXT_POINT ) {
  
  double auxS[36];//auxiliary array to save Eshelby tensor
    //Position dependent Eshelby tensor:
    this->giveSijkl( auxS, J, sort_a, nu , false);

    // The following components or Dijkl are (probably) obtained by expressing equation (11.41) from Mura book in Maple software.

    /*1st row*/
    D[_1111_] = this->multTRN*(-(X1*(this->_1Plus2nu*J1_1 + J1_11*X1 - SQR(sort_a[0])*(5*J11_1 + X1*J11_11)))) +
      auxS[_1111_];

    D[_1122_] = this->multTRN*(this->_2nu*J1_1*X1 + SQR(X1)*(-J1_22 + SQR(sort_a[0])*J11_22) +
                   X2*(-J2_2 + SQR(sort_a[0])*J21_2)) + auxS[_1122_];

    D[_1133_] = this->multTRN*(this->_2nu*J1_1*X1 + SQR(X1)*(-J1_33 + SQR(sort_a[0])*J11_33) +
                   X3*(-J3_3 + SQR(sort_a[0])*J31_3)) + auxS[_1133_];

    D[_1112_] = this->multTRN*(X1*(-3*J1_2 - J1_21*X1 + SQR(sort_a[0])*(3*J11_2 + X1*J11_21)) +
                   2*this->_1MinNu*J2_1*X2);

    D[_1123_] = this->multTRN*(SQR(X1)*(-J1_32 + SQR(sort_a[0])*J11_32) + X2*(-J2_3 + SQR(sort_a[0])*J21_3));

    D[_1113_] = this->multTRN*(X1*(-3*J1_3 - J1_31*X1 + SQR(sort_a[0])*(3*J11_3 + X1*J11_31)) +
                   2*this->_1MinNu*J3_1*X3);

    /*2nd row*/
    D[_2211_] = this->multTRN*(this->_2nu*J2_2*X2 + X1*(-J1_1 + SQR(sort_a[1])*J12_1) +
                   SQR(X2)*(-J2_11 + SQR(sort_a[1])*J22_11)) + auxS[_2211_];

    D[_2222_] = this->multTRN*(-(X2*(this->_1Plus2nu*J2_2 + J2_22*X2 - SQR(sort_a[1])*(5*J22_2 + X2*J22_22)))) +
      auxS[_2222_];

    D[_2233_] = this->multTRN*(this->_2nu*J2_2*X2 + SQR(X2)*(-J2_33 + SQR(sort_a[1])*J22_33) +
                   X3*(-J3_3 + SQR(sort_a[1])*J32_3)) + auxS[_2233_];

    D[_2212_] = this->multTRN*(this->_1Min2nu*J1_2*X1 + SQR(sort_a[1])*X1*J12_2 +
                   X2*(-2*J2_1 - J2_21*X2 + SQR(sort_a[1])*(2*J22_1 + X2*J22_21)));

    D[_2223_] = this->multTRN*(X2*(-3*J2_3 - J2_32*X2 + SQR(sort_a[1])*(3*J22_3 + X2*J22_32)) +
                   2*this->_1MinNu*J3_2*X3);

    D[_2213_] = this->multTRN*(X1*(-J1_3 + SQR(sort_a[1])*J12_3) + SQR(X2)*(-J2_31 + SQR(sort_a[1])*J22_31));


    /*3rd row*/
    D[_3311_] = this->multTRN*(this->_2nu*J3_3*X3 + X1*(-J1_1 + SQR(sort_a[2])*J13_1) +
                   SQR(X3)*(-J3_11 + SQR(sort_a[2])*J33_11)) + auxS[_3311_];

    D[_3322_] = this->multTRN*(this->_2nu*J3_3*X3 + X2*(-J2_2 + SQR(sort_a[2])*J23_2) +
                   SQR(X3)*(-J3_22 + SQR(sort_a[2])*J33_22)) + auxS[_3322_];

    D[_3333_] = this->multTRN*(-(X3*(this->_1Plus2nu*J3_3 + J3_33*X3 - SQR(sort_a[2])*(5*J33_3 + X3*J33_33)))) +
      auxS[_3333_];

    D[_3312_] = this->multTRN*(X1*(-J1_2 + SQR(sort_a[2])*J13_2) + SQR(X3)*(-J3_21 + SQR(sort_a[2])*J33_21));

    D[_3323_] = this->multTRN*(this->_1Min2nu*J2_3*X2 + SQR(sort_a[2])*X2*J23_3 +
                   X3*(-2*J3_2 - J3_32*X3 + SQR(sort_a[2])*(2*J33_2 + X3*J33_32)));

    D[_3313_] = this->multTRN*(this->_1Min2nu*J1_3*X1 + SQR(sort_a[2])*X1*J13_3 +
                   X3*(-2*J3_1 - J3_31*X3 + SQR(sort_a[2])*(2*J33_1 + X3*J33_31)));

    /*4th row*/
    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)));

    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)) +
                   2*this->_1MinNu*J2_1*X2);

    D[_1233_] = this->multTRN*(X1*(this->_2nu*J1_2 + X2*(-J2_33 + SQR(sort_a[0])*J12_33)));

    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 +
                   SQR(sort_a[0])*X2*J12_2 + SQR(sort_a[0])*X1*X2*J12_21) + auxS[_1212_];

    D[_1223_] = this->multTRN*(X1*(-J2_3 - J2_32*X2 + SQR(sort_a[0])*(J12_3 + X2*J12_32)) +
                   this->_1MinNu*J3_1*X3);

    D[_1213_] = this->multTRN*(X2*(-J2_3 - J2_31*X1 + SQR(sort_a[0])*(J12_3 + X1*J12_31)) +
                   this->_1MinNu*J3_2*X3);

    /*5th row*/
    D[_2311_] = this->multTRN*(X2*(this->_2nu*J2_3 + X3*(-J3_11 + SQR(sort_a[1])*J23_11)));

    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)));

    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)) +
                   2*this->_1MinNu*J3_2*X3);

    D[_2312_] = this->multTRN*(this->_1MinNu*J1_3*X1 + X3*(-J3_1 - J3_21*X2 + SQR(sort_a[1])*(J23_1 + X2*J23_21)));

    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 +
                   SQR(sort_a[1])*X3*J23_3 + SQR(sort_a[1])*X2*X3*J23_32) + auxS[_2323_];

    D[_2313_] = this->multTRN*(this->_1MinNu*J1_2*X1 + X2*(-J3_1 - J3_31*X3 + SQR(sort_a[1])*(J23_1 + X3*J23_31)));

    /*6th row*/
    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)));

    D[_1322_] = this->multTRN*(X1*(this->_2nu*J1_3 + X3* (-J3_22 + SQR(sort_a[0])*J13_22)));

    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)) +
                   2*this->_1MinNu*J3_1*X3);

    D[_1312_] = this->multTRN*(this->_1MinNu*J2_3*X2 + X3*(-J3_2 - J3_21*X1 + SQR(sort_a[0])*(J13_2 + X1*J13_21)));

    D[_1323_] = this->multTRN*(X1*(-J3_2 - J3_32*X3 + SQR(sort_a[0])*(J13_2 + X3*J13_32)) +
                   this->_1MinNu*J2_1*X2);

    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 +
                   SQR(sort_a[0])*X3*J13_3 + SQR(sort_a[0])*X1*X3*J13_31) + auxS[_1313_];

    // TomTag:
    // No Mandel notation.
    // eps_star_ij = D_ijkl * eps_tau_kl
    // e.g. eps_11 = ... + D_1112 eps_tau_12 + ... + D_1121 eps_tau_21
    // Therefore thre right half of matrix D is multiplied by two.
    D[_1112_] *= 2.0;
    D[_1123_] *= 2.0;
    D[_1113_] *= 2.0;
    D[_2212_] *= 2.0;
    D[_2223_] *= 2.0;
    D[_2213_] *= 2.0;
    D[_3312_] *= 2.0;
    D[_3323_] *= 2.0;
    D[_3313_] *= 2.0;
    D[_1212_] *= 2.0;
    D[_1223_] *= 2.0;
    D[_1213_] *= 2.0;
    D[_2312_] *= 2.0;
    D[_2323_] *= 2.0;
    D[_2313_] *= 2.0;
    D[_1312_] *= 2.0;
    D[_1323_] *= 2.0;
    D[_1313_] *= 2.0;
  // }
  // else _errorr( "giveDijkl: invalid position of a given point");
  
}//end of function: giveDijkl( )


//********************************************************************************************************
//                                            PROTECTED FUNCTIONS
//********************************************************************************************************

//********************************************************************************************************
// description: Function gives the displacement perturbation tensor of internal fields.
//              Function requiers "constructor2"
// last edit:   08. 08. 2013
//-------------------------------------------------------------------------------------------------------
// note:        @Lint[18]  - Perturbation (Eshelby's-like) displacement tensor to be evaluated and saved
//              @J[13]     - Ferers-Dysons' elliptic integrals
//              @sort_a[3] - sorted semiaxes' dimensions
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluUniformFieldProlateSpheroid::giveLijkINT( double Lint[18], const double J[13], const double sort_a[3],
                            const double x[3] )
{
  /* 1st row */
  Lint[_111_] = this->multTRN * ( X1 * ( this->_1Min2nu * J1 + 3. * SQR( sort_a[0] ) * J11 ) );
  Lint[_122_] = this->multTRN * ( X1 * ( this->_2nu * J1 - J2 + SQR( sort_a[0] ) * J12 ) );
  Lint[_133_] = this->multTRN * ( X1 * ( this->_2nu * J1 - J3 + SQR( sort_a[0] ) * J13 ) );
  Lint[_112_] = this->multTRN * ( X2 * ( this->_1Min2nu * J2 + SQR( sort_a[0] ) * J12 ) );
  Lint[_123_] = 0.;
  Lint[_113_] = this->multTRN * ( X3 * ( this->_1Min2nu * J3 + SQR( sort_a[0] ) * J13 ) );
  /* 2nd row */
  Lint[_211_] = this->multTRN * ( X2 * ( this->_2nu * J2 - J1 + SQR( sort_a[1] ) * J21 ) );
  Lint[_222_] = this->multTRN * ( X2 * ( this->_1Min2nu * J2 + 3. * SQR( sort_a[1] ) * J22 ) );
  Lint[_233_] = this->multTRN * ( X2 * ( this->_2nu * J2 - J3 + SQR( sort_a[1] ) * J23 ) );
  Lint[_212_] = this->multTRN * ( X1 * ( this->_1Min2nu * J1 + SQR( sort_a[1] ) * J21 ) );
  Lint[_223_] = this->multTRN * ( X3 * ( this->_1Min2nu * J3 + SQR( sort_a[1] ) * J23 ) );
  Lint[_213_] = 0.;
  /* 3rd row */
  Lint[_311_] = this->multTRN * ( X3 * ( this->_2nu * J3 - J1 + SQR( sort_a[2] ) * J31 ) );
  Lint[_322_] = this->multTRN * ( X3 * ( this->_2nu * J3 - J2 + SQR( sort_a[2] ) * J32 ) );
  Lint[_333_] = this->multTRN * ( X3 * ( this->_1Min2nu * J3 + 3. * SQR( sort_a[2] ) * J33 ) );
  Lint[_312_] = 0.;
  Lint[_323_] = this->multTRN * ( X2 * ( this->_1Min2nu * J2 + SQR( sort_a[2] ) * J32 ) );
  Lint[_313_] = this->multTRN * ( X1 * ( this->_1Min2nu * J1 + SQR( sort_a[2] ) * J31 ) );

  //Mandel notation adjusting
  /* 1nd row */
  Lint[_112_] *= 2;
  Lint[_123_] *= 2;
  Lint[_113_] *= 2;
  /* 2nd row */
  Lint[_212_] *= 2;
  Lint[_223_] *= 2;
  Lint[_213_] *= 2;
  /* 3rd row */
  Lint[_312_] *= 2;
  Lint[_323_] *= 2;
  Lint[_313_] *= 2;

  return ;
}//end of function: giveLijkINT( )

//********************************************************************************************************
// description: Function gives the displacement perturbation tensor of external fields.
//              Function requiers "constructor2"
// last edit:   08. 08. 2013
//-------------------------------------------------------------------------------------------------------
// note:        @Lext[18]  - Perturbation (generalized Eshelby's) displacement tensor to be evaluated
//                           and saved
//              @Lint[18]  - Perturbation (Eshelby's-like) displacement tensor
//              @dJi[9]    - first derivatives of Ferers-Dysons' elliptic integrals Ii
//              @dJij[27]  - first derivatives of Ferers-Dysons' elliptic integrals Iij
//              @sort_a[3] - sorted semiaxes' dimensions
//              @x[3]      - coordinates of a point
//              @nu        - Poisson's ratio of the matrix material
//********************************************************************************************************
void eshelbySoluUniformFieldProlateSpheroid::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])
{
  double xx1 = SQR( X1 ), xx2 = SQR( X2 ), xx3 = SQR( X3 ), x1x2 = X1 * X2, x2x3 = X2 * X3,
    x1x3 = X1 * X3;

  for( int i = 0; i < 18; i++ ){
    Lext[i] = Lint[i];
  }

  /* 1st row */
  Lext[_111_] += this->multTRN * xx1 * ( -J1_1 + SQR( sort_a[0] ) * J11_1 );
  Lext[_122_] += this->multTRN * x1x2 * ( -J2_2 + SQR( sort_a[0] ) * J12_2 );
  Lext[_133_] += this->multTRN * x1x3 * ( -J3_3 + SQR( sort_a[0] ) * J13_3 );
  Lext[_112_] += 2 * this->multTRN * xx1 * ( -J1_2 + SQR( sort_a[0] ) * J11_2 );
  Lext[_123_] += 2 * this->multTRN * x1x2 * ( -J2_3 + SQR( sort_a[0] ) * J12_3 );
  Lext[_113_] += 2 * this->multTRN * xx1 * ( -J1_3 + SQR( sort_a[0] ) * J11_3 );
  /* 2nd row */
  Lext[_211_] += this->multTRN * x1x2 * ( -J1_1 + SQR( sort_a[1] ) * J21_1 );
  Lext[_222_] += this->multTRN * xx2 * ( -J2_2 + SQR( sort_a[1] ) * J22_2 );
  Lext[_233_] += this->multTRN * x2x3 * ( -J3_3 + SQR( sort_a[1] ) * J23_3 );
  Lext[_212_] += 2 * this->multTRN * x1x2 * ( -J1_2 + SQR( sort_a[1] ) * J21_2 );
  Lext[_223_] += 2 * this->multTRN * xx2 * ( -J2_3 + SQR( sort_a[1] ) * J22_3 );
  Lext[_213_] += 2 * this->multTRN * x1x2 * ( -J1_3 + SQR( sort_a[1] ) * J21_3 );
  /* 3rd row */
  Lext[_311_] += this->multTRN * x1x3 * ( -J1_1 + SQR( sort_a[2] ) * J31_1 );
  Lext[_322_] += this->multTRN * x2x3 * ( -J2_2 + SQR( sort_a[2] ) * J32_2 );
  Lext[_333_] += this->multTRN * xx3 * ( -J3_3 + SQR( sort_a[2] ) * J33_3 );
  Lext[_312_] += 2 * this->multTRN * x1x3 * ( -J1_2 + SQR( sort_a[2] ) * J31_2 );
  Lext[_323_] += 2 * this->multTRN * x2x3 * ( -J2_3 + SQR( sort_a[2] ) * J32_3 );
  Lext[_313_] += 2 * this->multTRN * x1x3 * ( -J1_3 + SQR( sort_a[2] ) * J31_3 );

  return ;
}//end of function: giveLijkEXT( )

//********************************************************************************************************
// description: Function gives S_ijkl(x) position dependent Eshelby's tensor. (This is not
//              the strain perturbation tensor! Just its part ). The returned members of the
//              S_ijkl matrix are in tensorial (theoretical) notation!
//              Function requiers "constructor2"
// last edit:   08. 08. 2013
//-------------------------------------------------------------------------------------------------------
// note:        @S[36]     - Eshelby position dependent tensor
//              @J[13]     - Ferers-Dysons' elliptic integrals
//              @sort_a[3] - sorted semiaxes' dimensions
//              @nu        - Poisson's ratio of the matrix material
//********************************************************************************************************
void eshelbySoluUniformFieldProlateSpheroid::giveSijkl( double S[36], const double J[13], const double sort_a[3],
                                                  double nu , bool newFormulation)
{
  double mult = this->multTRN;

  // Mura book, eq (11.16)
  // Tensor symmetry: Sijkl = Sjikl = Sijlk

  //1st row
  S[_1111_] = mult * ( this->_1Min2nu * J1 + 3. * SQR( sort_a[0] ) * J11 );
  S[_1122_] = mult * ( this->_2nu * J1 - J2 + SQR( sort_a[0] ) * J21 ); // ?? J21 == J12 (give_Ii_and_Iij)
  S[_1133_] = mult * ( this->_2nu * J1 - J3 + SQR( sort_a[0] ) * J31 ); // ??
  S[_1123_] = S[_1113_] = S[_1112_] = 0.;
  //2nd row
  S[_2211_] = mult * ( this->_2nu * J2 - J1 + SQR( sort_a[1] ) * J12 );
  S[_2222_] = mult * ( this->_1Min2nu * J2 + 3. * SQR( sort_a[1] ) * J22 );
  S[_2233_] = mult * ( this->_2nu * J2 - J3 + SQR( sort_a[1] ) * J32 );
  S[_2223_] = S[_2213_] = S[_2212_] = 0.;
  //3rd row
  S[_3311_] = mult * ( this->_2nu * J3 - J1 + SQR( sort_a[2] ) * J13 );
  S[_3322_] = mult * ( this->_2nu * J3 - J2 + SQR( sort_a[2] ) * J23 );
  S[_3333_] = mult * ( this->_1Min2nu * J3 + 3. * SQR( sort_a[2] ) * J33 );
  S[_3323_] = S[_3313_] = S[_3312_] = 0.;
  //4th row
  S[_1211_] = S[_1222_] = S[_1233_] = 0.;
  S[_1212_] = mult * ( this->_1MinNu * J1 - nu * J2 + SQR( sort_a[0] ) * J12 );
  S[_1223_] = S[_1213_] = 0.;
    //5th row
  S[_2311_] = S[_2322_] = S[_2333_] = S[_2312_] = 0.;
  S[_2323_] = mult * ( this->_1MinNu * J2 - nu * J3 + SQR( sort_a[1] ) * J23 );
  S[_2313_] = 0.;
  //6th row
  S[_1311_] = S[_1322_] = S[_1333_] = S[_1312_] = S[_1323_] = 0.;
  S[_1313_] = mult * ( this->_1MinNu * J1 - nu * J3 + SQR( sort_a[0] ) * J13 ); // BUG?
  //S[_1313_] = mult * ( this->_1MinNu * J1 - nu * J3 + SQR( sort_a[2] ) * J13 );

  // Mura book: Sijkl = Sjikl = Sijlk

  return ;
}//end of function: giveSijkl( )

} // end of namespace mumech

/*end of file*/