eshelbySoluLambda.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluLambda.cpp
// description:    source file of the functions calculating the values of 'lambda' parameter and
//                 its derivatives
// author(s):      Jan Novak

// last edit:      02. 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cmath>

#include "macros.h"
#include "polynomialRootSolution.h"
#include "eshelbySoluLambda.h"

namespace mumech {

//********************************************************************************************************
// description: constructor
// last edit:   07. 09. 2009
//********************************************************************************************************
eshelbySoluLambda::eshelbySoluLambda( ){ }

//********************************************************************************************************
// description: destructor
// last edit:   07. 09. 2009
//********************************************************************************************************
eshelbySoluLambda::~eshelbySoluLambda( ){ }



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


//********************************************************************************************************
// description: Function gives the lambda value
// last edit:   12. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 72-73
//              @a[3]  - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3]  - coordinates of a point
//              @shape - inclusion shape
//********************************************************************************************************
double eshelbySoluLambda :: giveLambda (const double a[3], const double x[3], InclusionGeometry shape)
{
  return eshelbySoluLambda :: giveLambda(a, x[0], x[1], x[2], shape);
}

// Overloaded giveLambda method
double eshelbySoluLambda :: giveLambda (const double a[3], const double x1, const double x2, const double x3, InclusionGeometry shape)
{
  double lambda;
  if (shape == IS_ELLIPSOID  || shape == IS_OBLATE_SPHEROID || shape == IS_PROLATE_SPHEROID) {  // || shape == IS_ELLIPTIC_CYLINDER || shape == IS_CYLINDER){
//     if( ( SQR( x[0]/a[0]) + SQR( x[1]/a[1]) +  SQR( x[2]/a[2]) ) <= ( 1. + MACHINE_EPS ) ){
//       lambda = 0.;//patch: 23. 09. 2009
//     }
//     else{
      double a1 = a[0], a2 = a[1], a3 = a[2],
    aa = SQR( a1 ) + SQR( a2 ) + SQR( a3 ) - SQR( x1 ) - SQR( x2 ) - SQR( x3 ),
    b =  SQR( a1 ) * SQR( a2 ) + SQR( a1 ) * SQR( a3 ) + SQR( a2 ) * SQR( a3 ) -
    SQR( a2 ) * SQR( x1 ) - SQR( a3 ) * SQR( x1 ) - SQR( a1 ) * SQR( x2 ) -
    SQR( a3 ) * SQR( x2 ) - SQR( a1 ) * SQR( x3 ) - SQR( a2 ) * SQR( x3 ),
    c =  SQR( a1 ) * SQR( a2 ) * SQR( a3 ) - SQR( a2 ) * SQR( a3 ) * SQR( x1 ) -
    SQR( a1 ) * SQR( a3 ) * SQR( x2 ) - SQR( a1 ) * SQR( a2 ) * SQR( x3 );
    std::complex<double> roots[3];

      polynomialRootSolution giveRoots;

      //calculation of generally complex roots
      giveRoots.GetCubicPolyRoots( aa, b, c, roots );
      //elimination of complex roots
      roots[1].real(( ABS( imag(roots[1]) ) > MACHINE_EPS )? -INFTY : real(roots[1]));
      roots[2].real(( ABS( imag(roots[2]) ) > MACHINE_EPS )? -INFTY : real(roots[2]));

      lambda = ( MAX( MAX( real(roots[0]), real(roots[1]) ), real(roots[2]) ) );
      //    }
  }
  else if( shape == IS_SPHERE ){
    //this is most likely not necessary to evaluate in the case of spherical inclusion
    lambda = SQR( x1 ) + SQR( x2 ) + SQR( x3 ) - SQR( a[0] );
  }
  else{
    _errorr( "giveLambda: unknown inclusion shape\n");
  }

  return lambda;
}//end of function: giveLambda( )

//********************************************************************************************************
// description: Function gives generally complex lambda value
// last edit:   02. 08. 2010
//-------------------------------------------------------------------------------------------------------
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 72-73
//              @a[3]  - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3]  - coordinates of a point
//              @shape - inclusion shape
//********************************************************************************************************
std::complex<double> eshelbySoluLambda::giveLambda( const double a[3], std::complex<double> x[3], InclusionGeometry shape )
{
  std::complex<double> lambda;

  if (shape == IS_ELLIPSOID){
    double a1 = a[0], a2 = a[1], a3 = a[2];
    std::complex<double> x1 = x[0], x2 = x[1], x3 = x[2],
      aa = SQR( a1 ) + SQR( a2 ) + SQR( a3 ) - SQR( x1 ) - SQR( x2 ) - SQR( x3 ),
      b =  SQR( a1 ) * SQR( a2 ) + SQR( a1 ) * SQR( a3 ) + SQR( a2 ) * SQR( a3 ) -
      SQR( a2 ) * SQR( x1 ) - SQR( a3 ) * SQR( x1 ) - SQR( a1 ) * SQR( x2 ) -
      SQR( a3 ) * SQR( x2 ) - SQR( a1 ) * SQR( x3 ) - SQR( a2 ) * SQR( x3 ),
      c =  SQR( a1 ) * SQR( a2 ) * SQR( a3 ) - SQR( a2 ) * SQR( a3 ) * SQR( x1 ) -
      SQR( a1 ) * SQR( a3 ) * SQR( x2 ) - SQR( a1 ) * SQR( a2 ) * SQR( x3 ),
      roots[3];
    polynomialRootSolution giveRoots;

    //calculation of generally complex roots
    giveRoots.GetCubicPolyRoots( aa, b, c, roots );
    //seeking for maximum root
    lambda = ( CReMAX( CReMAX( roots[0], roots[1] ), roots[2] ) );
  }
  else if( shape == IS_SPHERE ){
    //this is most likely not necessary to evaluate in the case of spherical inclusion
    lambda = SQR( x[0] ) + SQR( x[1] ) + SQR( x[2] ) - SQR( a[0] );
  }
  else{
    _errorr( "giveLambda: unknown inclusion shape\n");
  }

  return lambda;
}//end of function: giveLambda( )


//********************************************************************************************************
// description: Function gives the first or seccond derivative of Lambda
// last edit:   07. 09. 2009
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 73
//              @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3] - coordinates of a point
//              @lambda - lambda parameter calculated according 'giveLambda()' function
//              @direction - derivative direction
//********************************************************************************************************
double eshelbySoluLambda::giveLambdaDerivative( const double a[3], const double x[3], double lambda,
                        derivativeDirection direction )
{
  double derivative = 0.;

  switch( direction ){
  case _x1_ :
    derivative = this->giveLambdaFirstDerivative( a, x, lambda, _x1_ );
    break;
  case _x1x1_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x1x1_ );
    break;
  case _x1x2_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x1x2_ );
    break;
  case _x1x3_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x1x3_ );
    break;
  case _x2_ :
    derivative = this->giveLambdaFirstDerivative( a, x, lambda, _x2_ );
    break;
  case _x2x1_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x2x1_ );
    break;
  case _x2x2_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x2x2_ );
    break;
  case _x2x3_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x2x3_ );
    break;
  case _x3_ :
    derivative = this->giveLambdaFirstDerivative( a, x, lambda, _x3_ );
    break;
  case _x3x1_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x3x1_ );
    break;
  case _x3x2_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x3x2_ );
    break;
  case _x3x3_ :
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x3x3_ );
    break;
  default:
    _errorr( "giveLambdaDerivative: invalid derivative direction ");
    break;
  }//end of switch: direction

  return( derivative );
}//end of function: giveLambdaDerivative()

//********************************************************************************************************
// description: Function gives the first or seccond derivative of Lambda
// last edit:   07. 09. 2009
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 73
//              @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3] - coordinates of a point
//              @lambda - lambda parameter calculated according 'giveLambda()' function
//              @direction_1 - direction of the first derivative: _x1_, _x2_, _x3_
//              @direction_2 - direction of the seccond derivative: _x1_, _x2_, _x3_ or must be '_empty_'
//                             to return the first derivative according 'direction_1' only
//********************************************************************************************************
double eshelbySoluLambda::giveLambdaDerivative( const double a[3], const double x[3], double lambda,
                        derivativeDirection direction_1,
                        derivativeDirection direction_2 )
{
  double derivative = 0.;

  if( direction_1 ==  _x1_ && direction_2 ==  _empty_ ){
    derivative = this->giveLambdaFirstDerivative( a, x, lambda, _x1_ );
  }
  else if( direction_1 ==  _x1_ && direction_2 ==  _x1_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x1x1_ );
  }
  else if( direction_1 ==  _x1_ && direction_2 ==  _x2_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x1x2_ );
  }
  else if( direction_1 ==  _x1_ && direction_2 ==  _x3_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x1x3_ );
  }
  else if( direction_1 ==  _x2_ && direction_2 ==  _empty_ ){
    derivative = this->giveLambdaFirstDerivative( a, x, lambda, _x2_ );
  }
  else if( direction_1 ==  _x2_ && direction_2 ==  _x1_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x2x1_ );
  }
  else if( direction_1 ==  _x2_ && direction_2 ==  _x2_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x2x2_ );
  }
  else if( direction_1 ==  _x2_ && direction_2 ==  _x3_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x2x3_ );
  }
  else if( direction_1 ==  _x3_ && direction_2 ==  _empty_ ){
    derivative = this->giveLambdaFirstDerivative( a, x, lambda, _x3_ );
  }
  else if( direction_1 ==  _x3_ && direction_2 ==  _x1_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x3x1_ );
  }
  else if( direction_1 ==  _x3_ && direction_2 ==  _x2_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x3x2_ );
  }
  else if( direction_1 ==  _x3_ && direction_2 ==  _x3_ ){
    derivative = this->giveLambdaSeccondDerivative( a, x, lambda, _x3x3_ );
  }
  else{
    _errorr( "giveLambdaDerivative: invalid derivative direction ");
  }

  return( derivative );
}//end of function: giveLambdaDerivative()

//********************************************************************************************************
// description: Function gives the first complex derivative of Lambda
// last edit:   02. 08. 2010
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 73
//              @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3] - coordinates of a point
//              @lambda - lambda parameter calculated according 'giveLambda()' function
//              @direction - derivative direction: _x1_, _x2_, _x3_
//********************************************************************************************************
std::complex<double> eshelbySoluLambda::c_dLambda( const double a[3], std::complex<double> x[3], std::complex<double> lambda,
                         derivativeDirection direction )
{
  int dir = direction - 1;

  return( 2. * x[dir] / ( ( SQR( a[dir] ) + lambda ) * ( SQR( x[0] )/SQR( SQR( a[0] ) + lambda ) +
                             SQR( x[1] )/SQR( SQR( a[1] ) + lambda ) +
                             SQR( x[2] )/SQR( SQR( a[2] ) + lambda ) ) ) );
}//end of function: c_dLambda()

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


//********************************************************************************************************
// description: Function gives the first derivative of Lambda
// last edit:   07. 09. 2009
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 73
//              @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3] - coordinates of a point
//              @lambda - lambda parameter calculated according 'giveLambda()' function
//              @direction - derivative direction: _x1_, _x2_, _x3_
//********************************************************************************************************
double eshelbySoluLambda::giveLambdaFirstDerivative( const double a[3], const double x[3], double lambda,
                             derivativeDirection direction )
{
  double a1 = a[0], a2 = a[1], a3 = a[2], x1 = x[0], x2 = x[1], x3 = x[2], derivative = 0.,
    var_a = 0, var_x = 0;

  //variable initiations
  switch( direction ){
  case _x1_ :
    var_x = x1;
    var_a = a1;
    break;
  case _x2_ :
    var_x = x2;
    var_a = a2;
    break;
  case _x3_ :
    var_x = x3;
    var_a = a3;
    break;
  default:
    _errorr( "giveLambdaFirstDerivative: invalid derivative direction ");
    break;
  }//end of switch: direction

  derivative = 2. * var_x / ( ( SQR( var_a ) + lambda ) * ( SQR( x1 )/SQR( SQR( a1 ) + lambda ) +
                                SQR( x2 )/SQR( SQR( a2 ) + lambda ) +
                                SQR( x3 )/SQR( SQR( a3 ) + lambda ) ) );
  return( derivative );
}//end of function: giveLambdaFirstDerivative()

//********************************************************************************************************
// description: Function gives the seccond derivative of Lambda
// last edit:   08. 09. 2009
// note:        Solution adopted from Toshio Mura's book (1982), page(s) 73
//              @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
//              @x[3] - coordinates of a point
//              @lambda - lambda parameter calculated according 'giveLambda()' function
//              @dLa - array with first derivatives of lambda
//              @direction - derivative direction: _xixj_; i,j = 1,2,3
//********************************************************************************************************
double eshelbySoluLambda::giveLambdaSeccondDerivative( const double a[3], const double x[3], double lambda,
                               derivativeDirection direction )
{
  double ai = a[0], aj = a[1], ak = a[2], xi = x[0], xj = x[1], xk = x[2], derivative = 0., dummy = 0.,
    dLa[3] = { this->giveLambdaFirstDerivative( a, x, lambda, _x1_ ),
           this->giveLambdaFirstDerivative( a, x, lambda, _x2_ ),
           this->giveLambdaFirstDerivative( a, x, lambda, _x3_ ) };
  double dLai = dLa[0], dLaj = dLa[1], dLak = dLa[2];
  derivativeDirection dirType;

  //variable initiations
  switch( direction ){
  case _x1x1_ :
    dirType = _HOMOGENEOUS_;
    //default set-up
    break;
  case _x1x2_ :
    dirType = _MIXED_;
    //default set-up
    break;
  case _x1x3_ :
    dirType = _MIXED_;
    SWAP( dLaj, dLak, dummy );
    SWAP( xj, xk, dummy );
    SWAP( aj, ak, dummy );
    break;
  case _x2x1_ :
    dirType = _MIXED_;
    //default set-up
    break;
  case _x2x2_ :
    dirType = _HOMOGENEOUS_;
    SWAP( dLai, dLaj, dummy );
    SWAP( xi, xj, dummy );
    SWAP( ai, aj, dummy );
    break;
  case _x2x3_ :
    dirType = _MIXED_;
    SWAP( dLai, dLak, dummy );
    SWAP( xi, xk, dummy );
    SWAP( ai, ak, dummy );
    break;
  case _x3x1_ :
    dirType = _MIXED_;
    SWAP( dLaj, dLak, dummy );
    SWAP( xj, xk, dummy );
    SWAP( aj, ak, dummy );
    break;
  case _x3x2_ :
    dirType = _MIXED_;
    SWAP( dLai, dLak, dummy );
    SWAP( xi, xk, dummy );
    SWAP( ai, ak, dummy );
    break;
  case _x3x3_ :
    dirType = _HOMOGENEOUS_;
    SWAP( dLai, dLak, dummy );
    SWAP( xi, xk, dummy );
    SWAP( ai, ak, dummy );
    break;
  default:
    _errorr( "giveLambdaSeccondDerivative: invalid derivative direction ");
    break;
  }//end of switch: direction

  //homogeneous or mixed derivative formula
  if( dirType == _MIXED_ ){
    derivative =
      ( -2. * xi * dLaj / SQR( SQR( ai ) + lambda )
    - 2. * xj * dLai / SQR( SQR( aj ) + lambda )
    + 2. * SQR( xi ) * dLai * dLaj / CUB( SQR( ai ) + lambda )
    + 2. * SQR( xj ) * dLai * dLaj / CUB( SQR( aj ) + lambda )
    + 2. * SQR( xk ) * dLai * dLaj / CUB( SQR( ak ) + lambda ) )
      /
      ( SQR( x[0] )/SQR( SQR( a[0] ) + lambda ) +
    SQR( x[1] )/SQR( SQR( a[1] ) + lambda ) +
    SQR( x[2] )/SQR( SQR( a[2] ) + lambda ) );
  }
  else if( dirType ==  _HOMOGENEOUS_ ){
    derivative =
      ( 2. / ( SQR( ai ) + lambda )
    - 4. * xi * dLai / SQR( SQR( ai ) + lambda )
    + 2. * SQR( xi ) * SQR( dLai ) / CUB( SQR( ai ) + lambda )
    + 2. * SQR( xj ) * SQR( dLai ) / CUB( SQR( aj ) + lambda )
    + 2. * SQR( xk ) * SQR( dLai ) / CUB( SQR( ak ) + lambda ) )
      /
      ( SQR( x[0] ) / SQR( SQR( a[0] ) + lambda ) +
    SQR( x[1] ) / SQR( SQR( a[1] ) + lambda ) +
    SQR( x[2] ) / SQR( SQR( a[2] ) + lambda ) );
  }else ;

  return( derivative );
}//end of function: giveLambdaSeccondDerivative()


} // end of namespace mumech

/*end of file*/