esuf_Penny.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluUniformFieldPenny.cpp
// description:    source file of the function returning the Eshelby solution of an ellipsoidal inclusion
//                 loaded by the uniform strain/stress field
// author(s):      Jan Novak

// last edit:      11. 08. 2010
// 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_Penny.h"
#include "esei_Penny.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:   16. 10. 2009
//-------------------------------------------------------------------------------------------------------
// note:        Solution adopted from Toshio Mura's book (1982) and Eshelby article (1957)
//              @sort_a[3] - sorted flat-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 eshelbySoluUniformFieldPenny::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. * __I(_Int33_) + ( 1. - 2. * nu ) * __I(_Int3_) ); // J[13] for flat shape contains I33 * a3^2
    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: eshelbyTensCompUniformFieldPenny( )



//********************************************************************************************************
// description: Function gives D_ijkl(x) generalized Eshelby's tensor.
//              Function requiers "constructor2"
// last edit:   12. 07. 2010
//-------------------------------------------------------------------------------------------------------
// 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
//              @nu        - Poisson's ratio of the matrix material
//              @pos       - position of the pint (_INT_POINT_/_EXT_POINT_)
//********************************************************************************************************
void eshelbySoluUniformFieldPenny::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])
{
  _errorr("termit: misto teto fce byla puvodne volana rodicovska, proc?");

  // Temporary test of giveSijkl () implementation.
  const int length = 36;
  double S1[length];
  double S2[length];
  this->giveSijkl( S1, J, sort_a, nu, false );
  this->giveSijkl( S2, J, sort_a, nu, true );
  double err = 0.0;
  for (int i=0; i < length; i++)
    err += (S2[i] - S1[i]) * (S2[i] - S1[i]);
  
  
  //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 );

    for (int i = 0; i < 36; i++)
      D[i] = auxS[i];

    // 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");
  
  return ;
}//end of function: giveDijkl( )



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

//********************************************************************************************************
// description: Function gives the displacement perturbation tensor of internal fields.
//              Function requiers "constructor2"
// last edit:   15 07. 2010
//-------------------------------------------------------------------------------------------------------
// 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 eshelbySoluUniformFieldPenny::giveLijkINT( double Lint[18], const double J[13], const double sort_a[3],
                            const double x[3] )
{
  _errorr("termit: misto teto fce byla puvodne volana rodicovska, proc?");

  CleanVector(Lint, 18); // No displacements of external points for flat shape.
  return ;
}//end of function: giveLijkINT( )

//********************************************************************************************************
// description: Function gives the displacement perturbation tensor of external fields.
//              Function requiers "constructor2"
// last edit:   15. 07. 2010
//-------------------------------------------------------------------------------------------------------
// 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 eshelbySoluUniformFieldPenny::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])
{
  //UNREFERENCED_PARAMETER(nu);
  CleanVector(Lext, 18); // No displacements of external points for flat shape.
  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:   12. 07. 2010
//-------------------------------------------------------------------------------------------------------
// 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 eshelbySoluUniformFieldPenny::giveSijkl( double S[36], const double J[13], const double sort_a[3],
                          double nu, bool newFormulation )
{
  _errorr("termit: misto teto fce byla puvodne volana rodicovska, proc?");

  double mult = this->multTRN;

  if (!newFormulation)
    {
      // Original formulation of S tensor.
      // Differs from Mura book Eq.(11.16), but Eshelby tensor is OK when compared to HELP results.
      //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 );
      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. * J33 ); // J[13] for flat shape contains I33 * a3^2
      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 );
    }
  else
    {
      // New formulation (15.08.2013), corresponds to Eq.(11.16).
      //1st row
      S[_1111_] = mult * ( this->_1Min2nu * J1 + 3. * SQR( sort_a[0] ) * J11 );
      S[_1122_] = mult * ( this->_2nu * J1 - J1 + SQR( sort_a[1] ) * J12 );
      S[_1133_] = mult * ( this->_2nu * J1 - J1 + SQR( sort_a[2] ) * J13 );
      S[_1123_] = S[_1113_] = S[_1112_] = 0.;
      //2nd row
      S[_2211_] = mult * ( this->_2nu * J2 - J2 + SQR( sort_a[0] ) * J12 );
      S[_2222_] = mult * ( this->_1Min2nu * J2 + 3. * SQR( sort_a[1] ) * J22 );
      S[_2233_] = mult * ( this->_2nu * J2 - J2 + SQR( sort_a[2] ) * J23 );
      S[_2223_] = S[_2213_] = S[_2212_] = 0.;
      //3rd row
      S[_3311_] = mult * ( this->_2nu * J3 - J3 + SQR( sort_a[0] ) * J31 );
      S[_3322_] = mult * ( this->_2nu * J3 - J3 + SQR( sort_a[1] ) * J23 );
      S[_3333_] = mult * ( this->_1Min2nu * J3 + 3. * J33 ); // J[13] for flat shape contains I33 * a3^2
      S[_3323_] = S[_3313_] = S[_3312_] = 0.;
      //4th row
      S[_1211_] = S[_1222_] = S[_1233_] = 0.;
      S[_1212_] = mult * ( this->_1Min2nu / 2.0 * (J1 + J2) + ( SQR( sort_a[0] ) + SQR( sort_a[1] ) ) / 2.0 * J12 );
      S[_1223_] = S[_1213_] = 0.;
      //5th row
      S[_2311_] = S[_2322_] = S[_2333_] = S[_2312_] = 0.;
      S[_2323_] = mult * ( this->_1Min2nu / 2.0 * (J2 + J3) + ( SQR( sort_a[1] ) + SQR( sort_a[2] ) ) / 2.0 * J23 );
      S[_2313_] = 0.;
      //6th row
      S[_1311_] = S[_1322_] = S[_1333_] = S[_1312_] = S[_1323_] = 0.;
      S[_1313_] = mult * ( this->_1Min2nu / 2.0 * (J1 + J3) + ( SQR( sort_a[0] ) + SQR( sort_a[2] ) ) / 2.0 * J13 );
    }

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

} // end of namespace mumech

/*end of file*/