esuf_Cylinder.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluUniformFieldCylinder.cpp
// author(s):      Jan Novak, Tomas Janda
// 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_Cylinder.h"
#include "esei_Cylinder.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 circular cylinder semiaxes (a2>a3, a1 -> infinity (not used))
//              @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 eshelbySoluUniformFieldCylinder::eshelbyTensCompUniformField( const double sort_a[3],
                                         const double eSInt[13], double nu,
                                         EshelbyTensComponent flag )
{
  double a2 = sort_a[1];
  double S_ijkl = 0.;

  // The following expressions holds for elliptic cylinder
  // I1 = 0.0
  // I11 * a1^2 = 0.0
  // I12 * a1^2 = I2
  // I13 * a1^2 = I3 (= I2)

  switch( flag ){

    /* 1st row */
  case _S1111_ :
    S_ijkl = MULT * ( 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( a2 ) * __I(_Int13_) - ( 1. - 2. * nu ) * __I(_Int1_) );
    break;
    /* 2nd row */
  case _S2211_ :
    S_ijkl = MULT * (      __I(_Int2_) - ( 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( a2 ) * __I(_Int23_) - ( 1. - 2. * nu ) * __I(_Int2_) );
    break;
    /* 3rd row */
  case _S3311_ :
    S_ijkl = MULT * (      __I(_Int3_) - ( 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 * ( 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( a2 ) * __I(_Int12_) + __I(_Int2_));
    break;
    /* 5th row */
  case _S2323_ :
    S_ijkl = .5 * MULT * ( ( 1. - 2. * nu ) * ( __I(_Int2_) + __I(_Int3_) ) +
               ( SQR( a2 ) + SQR( a2 ) ) * __I(_Int23_) );
    break;
    /* 6th row */
  case _S1313_ :
    S_ijkl = .5 * MULT * ( ( 1. - 2. * nu ) * ( __I(_Int1_) + __I(_Int3_) ) + SQR( a2 ) * __I(_Int13_) + __I(_Int3_));
    break;
  default:
    S_ijkl = 0.;
    break;
  }//end of switch: flag

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

//********************************************************************************************************
// 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 flat-ellipsoidal semiaxes (a1>a2>>a3)
//              @eInt[13] - values of elliptical integrals (potentials)
//              @nu - Poisson's ratio
//********************************************************************************************************
void eshelbySoluUniformFieldCylinder :: eshelbyTensUniformField (double eshTens[12],
                                                                 const double sort_a[3], const double eInt[13])
{
  double S[36];
  this->giveSijkl(S, eInt, sort_a, nu, false);

  /* 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_ );

  // TomTag
  // Convert S[36] to eshTens[12]
  eshTens[0] = S[0]; eshTens[1] = S[1]; eshTens[2] = S[2];
  eshTens[3] = S[6]; eshTens[4] = S[7]; eshTens[5] = S[8];
  eshTens[6] = S[12]; eshTens[7] = S[13]; eshTens[8] = S[14];
  eshTens[9] = 2.0 * S[21]; eshTens[10] = 2.0 * S[28]; eshTens[11] = 2.0 * S[35];
return ;
}//end of function: eshelbyTensUniformFieldCylinder( )



//********************************************************************************************************
// 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
//********************************************************************************************************
void eshelbySoluUniformFieldCylinder::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.

  // The following expressions holds for elliptic cylinder
  // I1 = 0.0
  // I11 * a1^2 = 0.0
  // I12 * a1^2 = I2
  // I13 * a1^2 = I3 (= I2)

    /*1st row*/
    D[_1111_] = auxS[_1111_];

    D[_1122_] = auxS[_1122_];

    D[_1133_] = auxS[_1133_];

    D[_1112_] = this->multTRN*(2*this->_1MinNu*J2_1*X2);

    D[_1123_] = 0.0;

    D[_1113_] = this->multTRN*(2*this->_1MinNu*J3_1*X3);

    /*2nd row*/
    D[_2211_] = this->multTRN*(this->_2nu*J2_2*X2 + 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*(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*(SQR(X2)*(-J2_31 + SQR(sort_a[1])*J22_31));


    /*3rd row*/
    D[_3311_] = this->multTRN*(this->_2nu*J3_3*X3 + 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*(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*(X3*(-2*J3_1 - J3_31*X3 + SQR(sort_a[2])*(2*J33_1 + X3*J33_31)));

    /*4th row*/
    D[_1211_] = 0.0;

    D[_1222_] = this->multTRN*(2*this->_1MinNu*J2_1*X2);

    D[_1233_] = 0.0;
      
    D[_1212_] = this->multTRN*(this->_1MinNu*J2_2*X2) + auxS[_1212_];

    D[_1223_] = this->multTRN*(this->_1MinNu*J3_1*X3);

    D[_1213_] = this->multTRN*(this->_1MinNu*J3_2*X3);

    /*5th row*/
    D[_2311_] = this->multTRN*(this->_2nu*J2_3*X2);
      
    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*(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*(X2*(-J3_1 - J3_31*X3 + SQR(sort_a[1])*(J23_1 + X3*J23_31)));

    /*6th row*/
    D[_1311_] = 0.0;

    D[_1322_] = 0.0;

    D[_1333_] = this->multTRN*(2*this->_1MinNu*J3_1*X3);

    D[_1312_] = this->multTRN*(this->_1MinNu*J2_3*X2);

    D[_1323_] = this->multTRN*(this->_1MinNu*J2_1*X2);

    D[_1313_] = this->multTRN*(this->_1MinNu*J3_3*X3) + 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");
  
  
  #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( )


//********************************************************************************************************
//                                            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 eshelbySoluUniformFieldCylinder::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 eshelbySoluUniformFieldCylinder::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 eshelbySoluUniformFieldCylinder::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 + 0.0 );  // I11 * a1^2 = 0.0 for cylinder
      S[_1122_] = mult * ( this->_2nu * J1 - J2 + J2 );      // I12 * a1^2 = I2 for cylinder
      S[_1133_] = mult * ( this->_2nu * J1 - J3 + J3 );      // I13 * a1^2 = I3 (=I2) for cylinder
      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 * J2 );
      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 * J3 );
  }
  else
  {
      // New formulation (15.08.2013), corresponds to Eq.(11.16).
      // ToDo: transrorm following expressions for a1 = infinity

      //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 + J1 );
      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 + J2 );
      S[_2223_] = S[_2213_] = S[_2212_] = 0.;
      //3rd row
      S[_3311_] = mult * ( this->_2nu * J3 - J1 + J1 ); // I13 * a3^2 = I1 for cylinder
      S[_3322_] = mult * ( this->_2nu * J3 - J2 + J2 ); // I23 * a3^2 = I2 for cylinder
      S[_3333_] = mult * ( this->_1Min2nu * J3 + 0.0 ); // I33 * a3^2 = 0.0 for cylinder
      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) + J2 / 2.0 );
      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) + J1 / 2.0 );
  }

  #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( )

} // end of namespace mumech

/*end of file*/