esuf.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluUniformField.cpp
// author(s):      Jan Novak
// 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 "esuf.h"
#include "esei.h"
#include "eshelbySoluLambda.h"
#include "transformTensors.h"
#include "elementaryFunctions.h"
#include "matrix.h"
#include "inclusion.h"
#include "problem.h"
#include "polynomialRootSolution.h"


namespace mumech {

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


eshelbySoluUniformField :: eshelbySoluUniformField (const InclusionRecord3D *i)
{
  I = i;
  
  this->nu        = I->P->matrix->nu();
  this->_2nu      = 2. * nu;
  this->_1Plus2nu = 1. + this->_2nu;
  this->_1MinNu   = 1. - nu;
  this->mult      = 8. * PI * this->_1MinNu;
  this->MULT      = 1. / ( 8. * PI * ( 1. - nu ) );
  this->multTRN   = 1./this->mult;
  this->_1Min2nu  = 1. - this->_2nu;
  this->_2nuMin1  =  this->_2nu - 1.;
  this->_3Min4nu  = 3. - 2. * this->_2nu;
}


void eshelbySoluUniformField :: giveEshelbyDisplacementUniformField (double displacement[3],
                                     const double pertDispTens[18],
                                     const double unifStrain[6])
{
  //u_1
  displacement[0] = pertDispTens[0] * unifStrain[0] + pertDispTens[1] * unifStrain[1] +
    pertDispTens[2] * unifStrain[2] + pertDispTens[3] * unifStrain[3] + pertDispTens[4] * unifStrain[4] +
    pertDispTens[5] * unifStrain[5];
  //u_2
  displacement[1] = pertDispTens[6] * unifStrain[0] + pertDispTens[7] * unifStrain[1] +
    pertDispTens[8] * unifStrain[2] + pertDispTens[9] * unifStrain[3] + pertDispTens[10] * unifStrain[4] +
    pertDispTens[11] * unifStrain[5];
  //u_3
  displacement[2] = pertDispTens[12] * unifStrain[0] + pertDispTens[13] * unifStrain[1] +
    pertDispTens[14] * unifStrain[2] + pertDispTens[15] * unifStrain[3] + pertDispTens[16] * unifStrain[4] +
    pertDispTens[17] * unifStrain[5];
  
  return ;
}//end of function: giveEshelbyDisplacementUniformField( )

void eshelbySoluUniformField :: giveEshelbyStrainUniformField (double strain[6], const double pertTens[36],
                                   const double unifStrain[6])
{
  //e_11
  strain[0] = pertTens[0] * unifStrain[0] + pertTens[1] * unifStrain[1] + pertTens[2] * unifStrain[2] +
              pertTens[3] * unifStrain[3] + pertTens[4] * unifStrain[4] + pertTens[5] * unifStrain[5];
  //e_22
  strain[1] = pertTens[6] * unifStrain[0] + pertTens[7] * unifStrain[1] + pertTens[8] * unifStrain[2] +
              pertTens[9] * unifStrain[3] + pertTens[10] * unifStrain[4] + pertTens[11] * unifStrain[5];
  //e_33
  strain[2] = pertTens[12] * unifStrain[0] + pertTens[13] * unifStrain[1] + pertTens[14] * unifStrain[2] +
              pertTens[15] * unifStrain[3] + pertTens[16] * unifStrain[4] + pertTens[17] * unifStrain[5];
  /* !!! FEEP notation !!! e_ij = {e_11, e_22, e_33, e_12, e_32, e_31}^T */
  //2*e_12
  strain[3] = pertTens[18] * unifStrain[0] + pertTens[19] * unifStrain[1] + pertTens[20] * unifStrain[2] +
              pertTens[21] * unifStrain[3] + pertTens[22] * unifStrain[4] + pertTens[23] * unifStrain[5];
  //2*e_32
  strain[4] = pertTens[24] * unifStrain[0] + pertTens[25] * unifStrain[1] + pertTens[26] * unifStrain[2] +
              pertTens[27] * unifStrain[3] + pertTens[28] * unifStrain[4] + pertTens[29] * unifStrain[5];
  //2*e_31
  strain[5] = pertTens[30] * unifStrain[0] + pertTens[31] * unifStrain[1] + pertTens[32] * unifStrain[2] +
              pertTens[33] * unifStrain[3] + pertTens[34] * unifStrain[4] + pertTens[35] * unifStrain[5];
  
  return ;
}//end of function: giveEshelbyStrainUniformField( )

//
void eshelbySoluUniformField :: giveEshelbyStrainOfOnePoint (Point *point)
{
  //GLOBAL->LOCAL transform of the point coordinates
  I->transformCoords_G2L(point->x, point->loc_x);
  
  // lambda parameter
  point->la = eshelbySoluLambda::giveLambda( I->a, point->loc_x, I->shape );
  
  // Ferers-Dyson's elliptic integrals
  I->ellInt->giveEllipticIntegrals (point->eInt, point->la, false);
  // Derivatives of Ferers-Dyson's elliptic integrals
  I->ellInt->giveDerivativesOfEllipticIntegrals (point, false);
  
  // Eshelby tensor
  this->giveEshelbyTensor( point->S, point->eInt);
  // stress/strain perturbation tensor
  this->giveStrainPerturbTensor( point);
  
  // strain perturbation
  this->giveEshelbyStrainUniformField (point->strain[0], point->D, I->locEigStrain);
  //LOCAL->GLOBAL transformation of resulting fields
  I->rotateStrain_L2G(point->strain[0], point->strain[0]);
  
  return ;
}//end of function: giveEshelbyStrainOfOnePoint( )


//
void eshelbySoluUniformField :: giveEshelbyFieldsOfOnePoint (Point *point, int lc, int nlc,
                                 bool disp, bool strn)
{
  //GLOBAL->LOCAL transform of the point coordinates
  I->transformCoords_G2L(point->x, point->loc_x);
  
  // lambda parameter
  point->la = eshelbySoluLambda::giveLambda( I->a, point->loc_x, I->shape );
  
  // Ferers-Dyson's elliptic integrals
  I->ellInt->giveEllipticIntegrals (point->eInt, point->la, false);
  // Derivatives of Ferers-Dyson's elliptic integrals
  I->ellInt->giveDerivativesOfEllipticIntegrals (point, false);
  
  // CALCULATION OF STRAIN PERTURBATION FIELD
  if (strn) {
    // Eshelby tensor
    this->giveEshelbyTensor( point->S, point->eInt);
    // stress/strain perturbation tensor
    this->giveStrainPerturbTensor( point);
    
    for (int i=0; i<nlc; i++) {
      // LOCAL FIELDS
      this->giveEshelbyStrainUniformField (point->strain[i], point->D, I->locEigStrain_LC[i+lc]);
      // GLOBAL FIELDS
      I->rotateStrain_L2G(point->strain[i], point->strain[i]);
    }
  }
  
  // CALCULATION OF DISPLACEMENT PERTURBATION FIELD
  if (disp) {
    //displacement perturbation tensor
    this->giveDisplacementPerturbTensor_EXTpoint(point);
    
    for (int i=0; i<nlc; i++) {
      // LOCAL FIELDS
      this->giveEshelbyDisplacementUniformField (point->displacement[i], point->L, I->locEigStrain_LC[i+lc]);
      // GLOBAL FIELDS
      I->rotateDisplc_L2G(point->displacement[i], point->displacement[i]);
    }
  }
  
  return ;
}//end of function: giveEshelbyFieldsOfOnePoint( )


//
void eshelbySoluUniformField :: giveEshelbyDisplacementOfOnePoint (double **globPert_displc, const double *coords, int lc, int nlc)
{
  double loc_coords[3];
  I->transformCoords_G2L(coords, loc_coords);

  double L[18];
  this->giveDisplacementPerturbTensor_INTpoint (L, loc_coords);

  for (int i=0; i<nlc; i++) {
    this->giveEshelbyDisplacementUniformField (globPert_displc[i], L, I->locEigStrain_LC[i+lc]);
    I->rotateDisplc_L2G (globPert_displc[i], globPert_displc[i]);
  }
}


//
void eshelbySoluUniformField :: giveEshelbyTensor (double S[12], const double eInt[13])
{
  this->eshelbyTensUniformField (S, I->a, eInt);
}


//********************************************************************************************************
//                                            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 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 eshelbySoluUniformField :: 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: eshelbyTensCompUniformFieldEllipsoid( )

//********************************************************************************************************
// description: Function initialize an 1D array by Eshelby tensor components
// last edit:   15. 02. 2010
//-------------------------------------------------------------------------------------------------------
// note:        Solution adopted from Toshio Mura's book (1982) and Eshelby article  (1957)
//              @eshTens - pointer to array of Eshelby tensor components S_ijkl
//              @sort_a[3] - sorted ellipsoidal semiaxes (a1>a2>a3)
//              @eInt[13] - values of elliptical integrals (potentials)
//********************************************************************************************************
void eshelbySoluUniformField :: eshelbyTensUniformField (double eshTens[12],
                             const double sort_a[3], const double eInt[13])
{
  /* 1st row */
  eshTens[0] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S1111_ );
  eshTens[1] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S1122_ );
  eshTens[2] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S1133_ );
  /* 2nd row */
  eshTens[3] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S2211_ );
  eshTens[4] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S2222_ );
  eshTens[5] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S2233_ );
  /* 3rd row */
  eshTens[6] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S3311_ );
  eshTens[7] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S3322_ );
  eshTens[8] = this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S3333_ );

  /* !!! FEEP - Mandel notation !!! e_ij = {e_11, e_22, e_33, 2e_12, 2e_32, 2e_31}^T */
  /* 4th row */
  eshTens[9] = 2. * this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S1212_ );

  /* 5th row */
  eshTens[10] = 2. * this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S2323_ );

  /* 6th row */
  eshTens[11] = 2. * this->eshelbyTensCompUniformField( sort_a, eInt, nu, _S1313_ );
  
  return ;
}//end of function: eshelbyTensUniformFieldEllipsoid( )


//
void eshelbySoluUniformField :: giveEshelbyTensorInverse (double SInv[12], const double S[12])
{
  giveInverseMatrix6x6to12(SInv,S);
}

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


void eshelbySoluUniformField :: giveStrainPerturbTensor (Point *point)
{
  this->giveDijkl (point->D, point->S, point->eInt, point->dJi, point->dJij, point->ddJi, point->ddJij, I->a, point->loc_x );
}

//********************************************************************************************************
// 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 eshelbySoluUniformField :: giveSijkl( double S[36], const double J[13], const double sort_a[3], double nu, bool newFormulation )
{
  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. * 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 );
    }
  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. * 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->_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 );
    }

#ifdef TOM_DEBUG
  std::cout << "Vypis slozek tenzoru S. (giveSijkl())" << std::endl;
  for (int i = 0; i < 36; i++)
    std::cout << "S[" << i << "] = " << S[i] << std::endl;
  std::cout << std::endl;
#endif

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


//********************************************************************************************************
// 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!
//              See:
// last edit:   20. 11. 2014
//-------------------------------------------------------------------------------------------------------
// note:        @S[36]        - Eshelby position dependent tensor
//              @sort_a[3]    - sorted semiaxes' dimensions
//              @stiffMat[36] - Full stiffness matrix of the matrix material
//              @M_partition  - number of points in Gauss quadrature
//              @N_partition  - number of points in Gauss quadrature
//********************************************************************************************************
void  eshelbySoluUniformField :: giveSijkl( double S[36], const double sort_a[3], const double stiffMat[36], int M_partition, int N_partition )
{

  if( S == NULL )  _errorr( "Field for Eshelby tensor must be allocated" );
  if( sort_a == NULL )  _errorr( "Field for sorted semiaxes must be allocated" );
  if( stiffMat == NULL )  _errorr( "Field for stiffness matrix of the matrix must be allocated" );


  double xiBar[3], zeta[3], K[3][3], N[3][3];
  double epsilon[3][3][3] = {};
  double STen[3][3][3][3] = {};
  double G[3][3][3][3] = {};
  int idx[3][3] = {};
  double mltIdx[6][6] = {};
  double D;

  double *x1 = new double[N_partition];
  double *w1 = new double[N_partition];
  polynomialRootSolution::GauLegF(x1, w1, N_partition, 0.0, 2.0*PI);
  
  double *x2 = new double[M_partition];
  double *w2 = new double[M_partition];
  polynomialRootSolution::GauLegF(x2, w2, M_partition, -1.0, 1.0);

  //initialization of permutation tensor
  epsilon[0][1][2] = epsilon[1][2][0] = epsilon[2][0][1] = 1.;
  epsilon[0][2][1] = epsilon[2][1][0] = epsilon[1][0][2] = -1;

  // idx initialization for 4order tensor transformation
  idx[0][0] = 0; idx[0][1] = 3; idx[0][2] = 5;
  idx[1][0] = 3; idx[1][1] = 1; idx[1][2] = 4;
  idx[2][0] = 5; idx[2][1] = 4; idx[2][2] = 2;

  // value which for multiplication of stiffness matrix element to get tensorial value
  mltIdx[0][0] = 1.; mltIdx[0][1] = 1.; mltIdx[0][2] = 1.; mltIdx[0][3] = .5; mltIdx[0][4] = .5; mltIdx[0][5] = .5;
  mltIdx[1][0] = 1.; mltIdx[1][1] = 1.; mltIdx[1][2] = 1.; mltIdx[1][3] = .5; mltIdx[1][4] = .5; mltIdx[1][5] = .5;
  mltIdx[2][0] = 1.; mltIdx[2][1] = 1.; mltIdx[2][2] = 1.; mltIdx[2][3] = .5; mltIdx[2][4] = .5; mltIdx[2][5] = .5;
  mltIdx[3][0] = 1.; mltIdx[3][1] = 1.; mltIdx[3][2] = 1.; mltIdx[3][3] = .5; mltIdx[3][4] = .5; mltIdx[3][5] = .5;
  mltIdx[4][0] = 1.; mltIdx[4][1] = 1.; mltIdx[4][2] = 1.; mltIdx[4][3] = .5; mltIdx[4][4] = .5; mltIdx[4][5] = .5;
  mltIdx[5][0] = 1.; mltIdx[5][1] = 1.; mltIdx[5][2] = 1.; mltIdx[5][3] = .5; mltIdx[5][4] = .5; mltIdx[5][5] = .5;

  int ii, jj;
  double dum;
  for (int p = 0; p < M_partition; p++) {
    for (int q = 0; q < N_partition; q++) {
      zeta[0] = sqrt(1.0 - x2[p] * x2[p]) * cos(x1[q]);
      zeta[1] = sqrt(1.0 - x2[p] * x2[p]) * sin(x1[q]);
      zeta[2] = x2[p];

      xiBar[0] = zeta[0] / sort_a[0];
      xiBar[1] = zeta[1] / sort_a[1];
      xiBar[2] = zeta[2] / sort_a[2];

      for (int i = 0; i < 3; i++)
        for (int k = 0; k < 3; k++) {
          dum = 0.0;
          for (int j = 0; j < 3; j++) {
            for (int l = 0; l < 3; l++) {
              ii = idx[i][j];
              jj = idx[k][l];
              dum += mltIdx[ii][jj]*stiffMat[ii*6+jj] * xiBar[j] * xiBar[l];
            }
          }
          K[i][k] = dum;
        }

      D = 0.0;
      for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++) {
          dum = 0.0;
          for (int k = 0; k < 3; k++) {
            D += epsilon[i][j][k] * K[i][0] * K[j][1] * K[k][2];
            for (int l = 0; l < 3; l++) {
              for (int m = 0; m < 3; m++) {
                for (int n = 0; n < 3; n++) {
                  dum += 0.5 * epsilon[i][k][l] * epsilon[j][m][n] * K[k][m] * K[l][n];
                }
              }
            }
          }
          N[i][j] = dum;
        }


       for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++)
           for (int k = 0; k < 3; k++)
             for (int l = 0; l < 3; l++) {
               G[i][j][k][l] = xiBar[k] * xiBar[l] * N[i][j] / D;
             }



       for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++) {
           N[i][j] = K[i][j] = 0.0; // cleaning

           for (int k = 0; k < 3; k++)
             for (int l = 0; l < 3; l++)
               for (int m = 0; m < 3; m++)
                 for (int n = 0; n < 3; n++) {
                   ii = idx[m][n];
                   jj = idx[k][l];
                   STen[i][j][k][l] += 1.0 / (8.0 * PI) * mltIdx[ii][jj]*stiffMat[ii*6+jj] * (G[i][m][j][n] + G[j][m][i][n]) * w1[q] * w2[p];
                 }

         }
    }
  }

  S[0]  = STen[0][0][0][0];
  S[1]  = STen[0][0][1][1];
  S[2]  = STen[0][0][2][2];
  S[3]  = (STen[0][0][0][1] + STen[0][0][1][0]);
  S[4]  = (STen[0][0][1][2] + STen[0][0][2][1]);
  S[5]  = (STen[0][0][0][2] + STen[0][0][2][0]);

  S[6]  = STen[1][1][0][0];
  S[7]  = STen[1][1][1][1];
  S[8]  = STen[1][1][2][2];
  S[9]  = (STen[1][1][0][1] + STen[1][1][1][0]);
  S[10] = (STen[1][1][1][2] + STen[1][1][2][1]);
  S[11] = (STen[1][1][0][2] + STen[1][1][2][0]);

  S[12] = STen[2][2][0][0];
  S[13] = STen[2][2][1][1];
  S[14] = STen[2][2][2][2];
  S[15] = (STen[2][2][0][1] + STen[2][2][1][0]);
  S[16] = (STen[2][2][1][2] + STen[2][2][2][1]);
  S[17] = (STen[2][2][0][2] + STen[2][2][2][0]);

  S[18] = STen[1][2][0][0];
  S[19] = STen[1][2][1][1];
  S[20] = STen[1][2][2][2];
  S[21] = (STen[0][1][0][1] + STen[0][1][1][0]);
  S[22] = (STen[0][1][1][2] + STen[0][1][2][1]);
  S[23] = (STen[0][1][0][2] + STen[0][1][2][0]);

  S[24] = STen[0][2][0][0];
  S[25] = STen[0][2][1][1];
  S[26] = STen[0][2][2][2];
  S[27] = (STen[1][2][0][1] + STen[1][2][1][0]);
  S[28] = (STen[1][2][1][2] + STen[1][2][2][1]);
  S[29] = (STen[1][2][0][2] + STen[1][2][2][0]);

  S[30] = STen[0][1][0][0];
  S[31] = STen[0][1][1][1];
  S[32] = STen[0][1][2][2];
  S[33] = (STen[0][2][0][1] + STen[0][2][1][0]);
  S[34] = (STen[0][2][1][2] + STen[0][2][2][1]);
  S[35] = (STen[0][2][0][2] + STen[0][2][2][0]);

  delete [] x1;
  delete [] w1;
  delete [] x2;
  delete [] w2;
}

//********************************************************************************************************
// description: Function gives D_ijkl(x) generalized Eshelby's tensor.
//              Function requiers "constructor2"
//              D_ijkl is adjusted to mandel notation. The resulting D tensor thus differs from the values in Mura book.
// 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
//              @pos       - position of the pint (_INT_POINT_/_EXT_POINT_)
//********************************************************************************************************
void eshelbySoluUniformField::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])
{
  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.

  // if and else variants are identical with differen arrangement
  if (1) {
    // 11 diagonala
    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[_2222_] = this->multTRN*(-X2*(this->_1Plus2nu*J2_2 + J2_22*X2 - SQR(sort_a[1])*(5*J22_2 + X2*J22_22))) + auxS[_2222_];
    D[_3333_] = this->multTRN*(-X3*(this->_1Plus2nu*J3_3 + J3_33*X3 - SQR(sort_a[2])*(5*J33_3 + X3*J33_33))) + auxS[_3333_];
  
    // 11
    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[_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[_2211_] = this->multTRN*( this->_2nu*J2_2*X2 + SQR(X2)*(-J2_11 + SQR(sort_a[1])*J22_11) + X1 *(-J1_1 + SQR(sort_a[1])*J12_1) ) + auxS[_2211_];
    D[_3311_] = this->multTRN*( this->_2nu*J3_3*X3 + SQR(X3)*(-J3_11 + SQR(sort_a[2])*J33_11) + X1 *(-J1_1 + SQR(sort_a[2])*J13_1) ) + auxS[_3311_];
    D[_3322_] = this->multTRN*( this->_2nu*J3_3*X3 + SQR(X3)*(-J3_22 + SQR(sort_a[2])*J33_22) + X2 *(-J2_2 + SQR(sort_a[2])*J23_2) ) + auxS[_3322_];
    
    // 12
    D[_1123_] = this->multTRN*( X2 *(-J2_3  + SQR(sort_a[0])*J21_3 ) + SQR(X1)*(-J1_32 + SQR(sort_a[0])*J11_32));
    D[_2213_] = this->multTRN*( X1 *(-J1_3  + SQR(sort_a[1])*J12_3 ) + SQR(X2)*(-J2_31 + SQR(sort_a[1])*J22_31));
    D[_3312_] = this->multTRN*( X1 *(-J1_2  + SQR(sort_a[2])*J13_2 ) + SQR(X3)*(-J3_21 + SQR(sort_a[2])*J33_21));
    
    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[_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);
    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[_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[_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)));
    
    // 21
    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[_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[_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[_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[_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[_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[_1233_] = this->multTRN*(X1*(this->_2nu*J1_2 + X2*(-J2_33 + SQR(sort_a[0])*J12_33)));
    D[_2311_] = this->multTRN*(X2*(this->_2nu*J2_3 + X3*(-J3_11 + SQR(sort_a[1])*J23_11)));
    D[_1322_] = this->multTRN*(X1*(this->_2nu*J1_3 + X3*(-J3_22 + SQR(sort_a[0])*J13_22)));
  
    // 22 diagonala
    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[_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[_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_];
    
    // 22 zbytek
    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[_2313_] = this->multTRN*( this->_1MinNu*J1_2*X1 + X2*(-J3_1 - J3_31*X3 + SQR(sort_a[1])*(J23_1 + X3*J23_31)));
    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*( this->_1MinNu*J2_1*X2 + X1*(-J3_2 - J3_32*X3 + SQR(sort_a[0])*(J13_2 + X3*J13_32)));
    D[_1223_] = this->multTRN*( this->_1MinNu*J3_1*X3 + X1*(-J2_3 - J2_32*X2 + SQR(sort_a[0])*(J12_3 + X2*J12_32)));
    D[_1213_] = this->multTRN*( this->_1MinNu*J3_2*X3 + X2*(-J2_3 - J2_31*X1 + SQR(sort_a[0])*(J12_3 + X1*J12_31)));
  }
  else {
    /*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;
    
  //D[_1112_] *= SQRT_2;
  //D[_1123_] *= SQRT_2;
  //D[_1113_] *= SQRT_2;
  //D[_2212_] *= SQRT_2;
  //D[_2223_] *= SQRT_2;
  //D[_2213_] *= SQRT_2;
  //D[_3312_] *= SQRT_2;
  //D[_3323_] *= SQRT_2;
  //D[_3313_] *= SQRT_2;
  //D[_1211_] *= SQRT_2;
  //D[_1222_] *= SQRT_2;
  //D[_1233_] *= SQRT_2;
  //D[_1212_] *= 2.;
  //D[_1223_] *= 2.;
  //D[_1213_] *= 2.;
  //D[_2311_] *= SQRT_2;
  //D[_2322_] *= SQRT_2;
  //D[_2333_] *= SQRT_2;
  //D[_2312_] *= 2.;
  //D[_2323_] *= 2.;
  //D[_2313_] *= 2.;
  //D[_1311_] *= SQRT_2;
  //D[_1322_] *= SQRT_2;
  //D[_1333_] *= SQRT_2;
  //D[_1312_] *= 2.;
  //D[_1323_] *= 2.;
  //D[_1313_] *= 2.;
    
#ifdef TOM_DEBUG
  // 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]);
    
  std::cout << "Vypis slozek tenzoru S. (ve funkci giveDijkl())" << std::endl;
  for (int i = 0; i < 12; i++)
    std::cout << "S[" << i << "] = " << S[i] << std::endl;
  std::cout << std::endl;
  std::cout << "Vypis slozek tenzoru D. (giveDijkl())" << std::endl;
  for (int i = 0; i < 36; i++)
    std::cout << "D[" << i << "] = " << D[i] << std::endl;
  std::cout << std::endl;
#endif
  
  return ;
}//end of function: giveDijkl( )


void eshelbySoluUniformField :: giveDisplacementPerturbTensor_EXTpoint (Point *point)
{
  this->giveLijkINT (point->L, point->eInt, I->a, point->loc_x );
  this->giveLijkEXT (point->L, point->L, point->dJi, point->dJij, I->a, point->loc_x);
}

void eshelbySoluUniformField :: giveDisplacementPerturbTensor_INTpoint (double L[18], const double x[3])
{
  this->giveLijkINT (L, I->eInt, I->a, x);
}


void eshelbySoluUniformField :: 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 );
  double xx2 = SQR( X2 );
  double xx3 = SQR( X3 );
  double x1x2 = X1 * X2;
  double x2x3 = X2 * X3;
  double x1x3 = X1 * X3;
  
  CopyVector(Lint, Lext, 18);
  
  /* 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( )

void eshelbySoluUniformField :: 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( )

} // end of namespace mumech

/*end of file*/