transformTensors.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           transformTensors.cpp
// author(s):      Jan Novak, Ladislav Svoboda
// 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 <math.h>
#include "macros.h"
#include "transformTensors.h"
#include "matrixOperations.h"
#include "elementaryFunctions.h"
#include "problem.h"


namespace mumech {

void TransformTensors :: give_T (double *T, const double *eullerAngles, EullerRotations flag, InclusionGeometry inclGeometry)
{
  switch ( inclGeometry ) {
  case IS_ELLIPSOID :
  case IS_OBLATE_SPHEROID :
  case IS_PROLATE_SPHEROID :
  //case IS_ELLIPTIC_CYLINDER :
  //case IS_CYLINDER :
  //case IS_PENNY :
  //case IS_FLAT_ELLIPSOID :
    TransformTensors::give_T_3d (T, eullerAngles, flag);
    break;
  case IS_SPHERE :
    giveUnitMatrix3x3(T);   //unit matrix, not used furthermore anyway (hopefully)
    break;
  case IS_CIRCLE :
    giveUnitMatrix2x2(T);   //unit matrix, not used furthermore anyway (hopefully)
    break;
  case IS_ELLIPSE :
    T[0] =  cos(eullerAngles[0]);  T[1] = sin(eullerAngles[0]);
    T[2] = -T[1];                  T[3] = T[0];
    break;
  default:
    _errorr2( "Unsupported shape of inclusion \"%s\"", IST_e2s (inclGeometry));
  }
  
  return ;
}

void TransformTensors :: give_Te (double *Te, const double *T, InclusionGeometry inclGeometry)
{
  switch ( inclGeometry ) {
    case IS_ELLIPSOID :
    case IS_OBLATE_SPHEROID :
    case IS_PROLATE_SPHEROID :
    //case IS_ELLIPTIC_CYLINDER :
    //case IS_CYLINDER :
    //case IS_PENNY :
    //case IS_FLAT_ELLIPSOID :
      TransformTensors::give_Te_3d (Te, T);
      break;
    case IS_SPHERE :
      MatrixOperations::giveUnitSMatrixFull_3d (Te);  // unit matrix, not used furthermore anyway (hopefully)
      break;
    case IS_CIRCLE :
      MatrixOperations::giveUnitSMatrixFull_2d (Te);  // unit matrix, not used furthermore anyway (hopefully)
      break;
    case IS_ELLIPSE :
      // first row
      Te[0] = SQR( T[0] );
      Te[1] = SQR( T[1] );
      Te[2] = 2.0 * T[0] * T[1];

      // second row
      Te[3] = SQR( T[2] );
      Te[4] = SQR( T[3] );
      Te[5] = 2.0 * T[2] * T[3];

      // third row
      Te[6] = T[0] * T[2];
      Te[7] = T[1] * T[3];
      Te[8] = T[0] * T[3] + T[1] * T[2];

      break;
    default:
      _errorr2( "Unsupported shape of inclusion \"%s\"", IST_e2s (inclGeometry));
  }
}//end of function: give_Te( )

void TransformTensors :: give_TInv (double *TInv, const double *T, InclusionGeometry inclGeometry)
{
  switch ( inclGeometry ) {

  case IS_ELLIPSOID :
  //case IS_ELLIPTIC_CYLINDER :
  //case IS_CYLINDER :
  //case IS_PENNY :
  //case IS_FLAT_ELLIPSOID :
  case IS_OBLATE_SPHEROID :
  case IS_PROLATE_SPHEROID :     giveTransposedMatrix (T, TInv, 3);  break;
  case IS_ELLIPSE :              giveTransposedMatrix (T, TInv, 2);  break;

  // unit matrix, not used furthermore anyway (hopefully)
  case IS_SPHERE :    giveUnitMatrix3x3(TInv);  break;
  case IS_CIRCLE :    giveUnitMatrix2x2(TInv);  break;

  default: _errorr2( "Unsupported shape of inclusion \"%s\"", IST_e2s (inclGeometry));
  }

  return ;
}

void TransformTensors :: give_TeInv (double *TeInv, const double *Te, InclusionGeometry inclGeometry)
{
  switch ( inclGeometry ) {
    case IS_ELLIPSOID :
    case IS_OBLATE_SPHEROID :
    case IS_PROLATE_SPHEROID :
    //case IS_ELLIPTIC_CYLINDER :
    //case IS_CYLINDER :
    //case IS_PENNY :
    //case IS_FLAT_ELLIPSOID :
      // ALTERNATIVE - more consumption then "transposition plus some 2.0 multiplication"
      // this->Te( TeInv, TInv );

      giveTransposedMatrix (Te, TeInv, 6);

      int i, j;

      for (i=0; i<3; i++)
        for (j=3; j<6; j++)
          TeInv[i*6+j] *= 2.0;

      for (i=3; i<6; i++)
        for (j=0; j<3; j++)
          TeInv[i*6+j] /= 2.0;

      break;

    case IS_ELLIPSE :
      giveTransposedMatrix (Te, TeInv, 3);

      TeInv[2] *= 2.0;  TeInv[5] *= 2.0;
      TeInv[6] /= 2.0;  TeInv[7] /= 2.0;

      break;

    case IS_SPHERE :
      MatrixOperations::giveUnitSMatrixFull_3d (TeInv);  // unit matrix, not used furthermore anyway (hopefully)
      break;
    case IS_CIRCLE :
      MatrixOperations::giveUnitSMatrixFull_2d (TeInv);  // unit matrix, not used furthermore anyway (hopefully)
      break;
    default:
      _errorr2( "Unsupported shape of inclusion \"%s\"", IST_e2s (inclGeometry));
  }

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


void TransformTensors :: give_T_3d (double T[9], const double eullerAngles[3], EullerRotations flag)
{
  double phi = eullerAngles[0], nu = eullerAngles[1], psi = eullerAngles[2];
  double c1 = cos( phi ), c2 = cos( nu ), c3 = cos( psi );
  double s1 = sin( phi ), s2 = sin( nu ), s3 = sin( psi );

  switch( flag ){
  case _XZX_:
    //1st row
    T[0] = c2;
    T[1] = c1 * s2;
    T[2] = s1 * s2;
    //2nd row
    T[3] = -c3 * s2;
    T[4] = c1 * c2 * c3 - s1 * s3;
    T[5] = c2 * c3 * s1 + c1 * s3;
    //3rd row
    T[6] = s2 * s3;
    T[7] = -c3 * s1 - c1 * c2 * s3;
    T[8] = c1 * c3 - c2 * s1 * s3;
    break;
  case _XYX_:
    //1st row
    T[0] = c2;
    T[1] = s1 * s2;
    T[2] = -c1 * s2;
    //2nd row
    T[3] = s2 * s3;
    T[4] = c1 * c3 - c2 * s1 * s3;
    T[5] = c3 * s1 + c1 * c2 * s3;
    //3rd row
    T[6] = c3 * s2;
    T[7] = -c2 * c3 * s1 - c1 * s3;
    T[8] = c1 * c2 * c3 - s1 * s3;
    break;
  case _YXY_:
    //1st row
    T[0] = c1 * c3 - c2 * s1 * s3;
    T[1] = s2 * s3;
    T[2] = -c3 * s1 - c1 * c2 * s3;
    //2nd row
    T[3] = s1 * s2;
    T[4] = c2;
    T[5] = c1 * s2;
    //3rd row
    T[6] = c2 * c3 * s1 + c1 * s3;
    T[7] = -c3 * s2;
    T[8] = c1 * c2 * c3 - s1 * s3;
    break;
  case _YZY_:
    //1st row
    T[0] = c1 * c2 * c3 - s1 * s3;
    T[1] = c3 * s2;
    T[2] = -c2 * c3 * s1 - c1 * s3;
    //2nd row
    T[3] = -c1 * s2;
    T[4] = c2;
    T[5] = s1 * s2;
    //3rd row
    T[6] = c3 * s1 + c1 * c2 * s3;
    T[7] = s2 * s3;
    T[8] = c1 * c3 - c2 * s1 * s3;
    break;
  case _ZYZ_:
    //1st row
    T[0] = c1 * c2 * c3 - s1 * s3;
    T[1] = c2 * c3 * s1 + c1 * s3;
    T[2] = -c3 * s2;
    //2nd row
    T[3] = -c3 * s1 - c1 * c2 * s3;
    T[4] = c1 * c3 - c2 * s1 * s3;
    T[5] = s2 * s3;
    //3rd row
    T[6] = c1 * s2;
    T[7] = s1 * s2;
    T[8] = c2;
    break;
  case _ZXZ_:
    //1st row
    T[0] = c1 * c3 - c2 * s1 * s3;
    T[1] = c3 * s1 + c1 * c2 * s3;
    T[2] = s2 * s3;
    //2nd row
    T[3] = -c2 * c3 * s1 - c1 * s3;
    T[4] = c1 * c2 * c3 - s1 * s3;
    T[5] = c3 * s2;
    //3rd row
    T[6] = s1 * s2;
    T[7] = -c1 * s2;
    T[8] = c2;
    break;
  case _XZY_:
    //1st row
    T[0] = c2 * c3;
    T[1] = c1 * c3 * s2 + s1 * s3;
    T[2] = c3 * s1 * s2 - c1 * s3;
    //2nd row
    T[3] = -s2;
    T[4] = c1 * c2;
    T[5] = c2 * s1;
    //3rd row
    T[6] = c2 * s3;
    T[7] = -c3 * s1 + c1 * s2 * s3;
    T[8] = c1 * c3 + s1 * s2 * s3;
    break;
  case _XYZ_:
    //1st row
    T[0] = c2 * c3;
    T[1] = c3 * s1 * s2 + c1 * s3;
    T[2] = -c1 * c3 * s2 + s1 * s3;
    //2nd row
    T[3] = -c2 * s3;
    T[4] = c1 * c3 - s1 * s2 * s3;
    T[5] = c3 * s1 + c1 * s2 * s3;
    //3rd row
    T[6] = s2;
    T[7] = -c2 * s1;
    T[8] = c1 * c2;
    break;
  case _YXZ_:
    //1st row
    T[0] = c1 * c3 + s1 * s2 * s3;
    T[1] = c2 * s3;
    T[2] = -c3 * s1 + c1 * s2 * s3;
    //2nd row
    T[3] = c3 * s1 * s2 - c1 * s3;
    T[4] = c2 * c3;
    T[5] = c1 * c3 * s2 + s1 * s3;
    //3rd row
    T[6] = c2 * s1;
    T[7] = -s2;
    T[8] = c1 * c2;
    break;
  case _YZX_:
    //1st row
    T[0] = c1 * c2;
    T[1] = s2;
    T[2] = -c2 * s1;
    //2nd row
    T[3] = -c1 * c3 * s2 + s1 * s3;
    T[4] = c2 * c3;
    T[5] = c3 * s1 * s2 + c1 * s3;
    //3rd row
    T[6] = c3 * s1 + c1 * s2 * s3;
    T[7] = -c2 * s3;
    T[8] = c1 * c3 - s1 * s2 * s3;
    break;
  case _ZYX_:
    //1st row
    T[0] = c1 * c2;
    T[1] = c2 * s1;
    T[2] = -s2;
    //2nd row
    T[3] = -c3 * s1 + c1 * s2 * s3;
    T[4] = c1 * c3 + s1 * s2 * s3;
    T[5] = c2 * s3;
    //3rd row
    T[6] = c1 * c3 * s2 + s1 * s3;
    T[7] = c3 * s1 * s2 - c1 * s3;
    T[8] = c2 * c3;
    break;
  case _ZXY_:
    //1st row
    T[0] = c1 * c3 - s1 * s2 * s3;
    T[1] = c3 * s1 + c1 * s2 * s3;
    T[2] = -c2 * s3;
    //2nd row
    T[3] = -c2 * s1;
    T[4] = c1 * c2;
    T[5] = s2;
    //3rd row
    T[6] = c3 * s1 * s2 + c1 * s3;
    T[7] = -c1 * c3 * s2 + s1 * s3;
    T[8] = c2 * c3;
    break;
  default:
    _errorr( "T: Invalid successive rotations of coordinate system");
    break;
  }//end of switch: flag

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


#define _T( i, j )  T[( i - 1 ) * 3 + ( j - 1 )]

void TransformTensors :: give_Te_3d (double Te[36], const double T[9])
{
  //1st row
  Te[0] = SQR( _T( 1, 1 ) );
  Te[1] = SQR( _T( 1, 2 ) );
  Te[2] = SQR( _T( 1, 3 ) );
  //T[3] <-> T[5] in TNTV_THEORETICAL_FEEP notation
  Te[5] = 2. * _T( 1, 1 ) * _T( 1, 3 );//s_13 transform
  Te[4] = 2. * _T( 1, 2 ) * _T( 1, 3 );//s_23 transform
  Te[3] = 2. * _T( 1, 1 ) * _T( 1, 2 );//s_12 transform

  //2nd row
  Te[6] = SQR( _T( 2, 1 ) );
  Te[7] = SQR( _T( 2, 2 ) );
  Te[8] = SQR( _T( 2, 3 ) );
  //T[9] <-> T[11] in TNTV_THEORETICAL_FEEP notation
  Te[11] = 2. * _T( 2, 1 ) * _T( 2, 3 );//s_13 transform
  Te[10] = 2. * _T( 2, 2 ) * _T( 2, 3 );//s_23 transform
  Te[9] =  2. * _T( 2, 1 ) * _T( 2, 2 );//s_12 transform

  //3rd row
  Te[12] = SQR( _T( 3, 1 ) );
  Te[13] = SQR( _T( 3, 2 ) );
  Te[14] = SQR( _T( 3, 3 ) );
  //T[15] <-> T[17] in TNTV_THEORETICAL_FEEP notation
  Te[17] = 2. * _T( 3, 1 ) * _T( 3, 3 );//s_13 transform
  Te[16] = 2. * _T( 3, 2 ) * _T( 3, 3 );//s_23 transform
  Te[15] = 2. * _T( 3, 1 ) * _T( 3, 2 );//s_12 transform

  //4th row (in TNTV_THEORETICAL_FEEP notation 4th row <-> 6th row)
  Te[18] = _T( 2, 1 ) * _T( 1, 1 );
  Te[19] = _T( 2, 2 ) * _T( 1, 2 );
  Te[20] = _T( 2, 3 ) * _T( 1, 3 );
  //T[21] <-> T[23] in TNTV_THEORETICAL_FEEP notation
  Te[23] = _T( 2, 3 ) * _T( 1, 1 ) + _T( 2, 1 ) * _T( 1, 3 );//s_13 transform
  Te[22] = _T( 2, 3 ) * _T( 1, 2 ) + _T( 2, 2 ) * _T( 1, 3 );//s_23 transform
  Te[21] = _T( 1, 1 ) * _T( 2, 2 ) + _T( 2, 1 ) * _T( 1, 2 );//s_12 transform

  //5th row
  Te[24] = _T( 3, 1 ) * _T( 2, 1 );
  Te[25] = _T( 3, 2 ) * _T( 2, 2 );
  Te[26] = _T( 3, 3 ) * _T( 2, 3 );
  //T[33] <-> T[35] in TNTV_THEORETICAL_FEEP notation
  Te[29] = _T( 3, 3 ) * _T( 2, 1 ) + _T( 3, 1 ) * _T( 2, 3 );//s_13 transform
  Te[28] = _T( 3, 2 ) * _T( 2, 3 ) + _T( 3, 3 ) * _T( 2, 2 );//s_23 transform
  Te[27] = _T( 3, 1 ) * _T( 2, 2 ) + _T( 3, 2 ) * _T( 2, 1 );//s_12 transform

  //6th row (in TNTV_THEORETICAL_FEEP notation 6th row <-> 4th row)
  Te[30] = _T( 3, 1 ) * _T( 1, 1 );
  Te[31] = _T( 3, 2 ) * _T( 1, 2 );
  Te[32] = _T( 3, 3 ) * _T( 1, 3 );
  //T[27] <-> T[29] in TNTV_THEORETICAL_FEEP notation
  Te[35] = _T( 3, 3 ) * _T( 1, 1 ) + _T( 3, 1 ) * _T( 1, 3 );//s_13 transform
  Te[34] = _T( 3, 2 ) * _T( 1, 3 ) + _T( 3, 3 ) * _T( 1, 2 );//s_23 transform
  Te[33] = _T( 3, 2 ) * _T( 1, 1 ) + _T( 3, 1 ) * _T( 1, 2 );//s_12 transform

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



} // end of namespace mumech

/*end of file*/