matrixOperations.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           matrixOperations.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 "matrixOperations.h"
#include "elementaryFunctions.h"
#include "macros.h"
#include <math.h>


namespace mumech {


double MatrixOperations :: giveTTproduct_1is2x2to3and2x2to3_SS  (const double *T1, const double *T2)
{
  return T1[0]*T2[0] + T1[1]*T2[1] + 2.0*T1[2]*T2[2];
}


//********************************************************************************************************
// description: function prints a field which is saved as 2D array
// last edit:   25. 01. 2010
// -------------------------------------------------------------------------------------------------------
// note:        @field - '@n x @m' matrix
//              @n - number of rows of printed field
//              @m - number of columns of printed field
//********************************************************************************************************
void MatrixOperations :: printField( double ** field, int n, int m, const char * notice )
{
  printf( "%s", notice );
  for( int i = 0; i < n; i++ ){
    for( int j = 0; j < m; j++ ){
      if( ABS( field[i][j] ) > MACHINE_EPS_PRN ){
    printf( "%.6e ", field[i][j] );
      }
      else{
    printf( "%.6e ", 0. );
      }
    }printf( "\n" );
  }

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

void MatrixOperations :: giveIsoTensAndVectorProduct (const double *tens, const double *vect, double *result, int dim)
{
  if( vect == result ){//in-place product
    double *auxVec = new double[dim];
    //product calculation
    if (dim == 3) {
      auxVec[0] = tens[0] * vect[0] + tens[3] * vect[1] + tens[5] * vect[2];
      auxVec[1] = tens[3] * vect[0] + tens[1] * vect[1] + tens[4] * vect[2];
      auxVec[2] = tens[5] * vect[0] + tens[4] * vect[1] + tens[2] * vect[2];
    }
    else if (dim == 2) {
      auxVec[0] = tens[0] * vect[0] + tens[2] * vect[1];
      auxVec[1] = tens[2] * vect[0] + tens[1] * vect[1];
    }
    else _errorr("unsupported dimension");
    
    //result initialization
    CopyVector(auxVec, result, dim);
    delete [] auxVec;
  }
  else{//ou-of-place product
    //product calculation and result initialization
    if (dim == 3) {
      result[0] = tens[0] * vect[0] + tens[3] * vect[1] + tens[5] * vect[2];
      result[1] = tens[3] * vect[0] + tens[1] * vect[1] + tens[4] * vect[2];
      result[2] = tens[5] * vect[0] + tens[4] * vect[1] + tens[2] * vect[2];
    }
    else if (dim == 2) {
      result[0] = tens[0] * vect[0] + tens[2] * vect[1];
      result[1] = tens[2] * vect[0] + tens[1] * vect[1];
    }
    else _errorr("unsupported dimension");
  }
  
  return ;
}//end of function: giveIsoTensAndVectorProduct( )

void MatrixOperations :: giveTVproduct_6is6x6to12and6 (double *result, const double *T, const double *vect)
{
  if( vect == result ){//in-place product
    double auxVect[6];
    //product calculation
    auxVect[0] = T[0] * vect[0] + T[1] * vect[1] + T[2] * vect[2];
    auxVect[1] = T[3] * vect[0] + T[4] * vect[1] + T[5] * vect[2];
    auxVect[2] = T[6] * vect[0] + T[7] * vect[1] + T[8] * vect[2];
    auxVect[3] = T[9] * vect[3];
    auxVect[4] = T[10] * vect[4];
    auxVect[5] = T[11] * vect[5];
    //result initialization
    CopyVector(auxVect, result, 6);
  }
  else{//ou-of-place product
    //product calculation and result initialization
    result[0] = T[0] * vect[0] + T[1] * vect[1] + T[2] * vect[2];
    result[1] = T[3] * vect[0] + T[4] * vect[1] + T[5] * vect[2];
    result[2] = T[6] * vect[0] + T[7] * vect[1] + T[8] * vect[2];
    result[3] = T[9] * vect[3];
    result[4] = T[10] * vect[4];
    result[5] = T[11] * vect[5];
  }

  return ;
}

void MatrixOperations :: giveTVproduct_3is3x3to5and3 (double *result, const double *T, const double *vect)
{
  if( vect == result ){//in-place product
    double auxVect[3];
    //product calculation
    auxVect[0] = T[0] * vect[0] + T[1] * vect[1];
    auxVect[1] = T[2] * vect[0] + T[3] * vect[1];
    auxVect[2] = T[4] * vect[2];
    //result initialization
    CopyVector(auxVect, result, 3);
  }
  else{//ou-of-place product
    //product calculation and result initialization
    result[0] = T[0] * vect[0] + T[1] * vect[1];
    result[1] = T[2] * vect[0] + T[3] * vect[1];
    result[2] = T[4] * vect[2];
  }
}

double MatrixOperations :: giveVTVproduct_1is3and3x3to5and3  (const double *T, const double *vect)
{
  return
    vect[0] * ( T[0] * vect[0] + T[1] * vect[1]) +
    vect[1] * ( T[2] * vect[0] + T[3] * vect[1]) +
    vect[2] * ( T[4] * vect[2]);
}

void MatrixOperations :: giveTTproduct_3x3to5is3x3to5and3x3to5 (double *result, const double *T1, const double *T2)
{
  if (result == T1  ||  result == T2) { //in-place product
    _errorr("in-place product not implemented yet");
  }
  else { //out-of-place product
    result[0] = T1[0] * T2[0]  +  T1[1] * T2[2];    result[1] = T1[0] * T2[1]  +  T1[1] * T2[3];
    result[2] = T1[2] * T2[0]  +  T1[3] * T2[2];    result[3] = T1[2] * T2[1]  +  T1[3] * T2[3];

    result[4] = T1[4] * T2[4];
  }
}







//********************************************************************************************************
// description: function prints a field which is saved as 2D array
// last edit:   24. 02. 2010
// -------------------------------------------------------------------------------------------------------
// note:        @field - '@n x @m' matrix
//              @n - number of rows of printed field
//              @m - number of columns of printed field
//              @lineBreak - a number after which the line break is expected (formating purposes)
//********************************************************************************************************
void MatrixOperations :: printField( double ** field, int n, int m, int lineBreak, const char * notice )
{
  int i = 0, j = 0, k = 0;

  printf( "%s", notice );
  for( i = 0; i < n; i++ ){
    for( j = 0; j < m; j++, k++ ){
      if( k == lineBreak ){
    printf( "\n" );
    k = 0;
      }else ;
      if( ABS( field[i][j] ) > MACHINE_EPS_PRN ){
    ( k == 0 )? printf( "%.6e", field[i][j] ) : printf( " %.6e", field[i][j] );
      }
      else{
    ( k == 0 )? printf( "%.6e", 0. ) : printf( " %.6e", 0. );
      }
    }printf( "\n\n" ); k = 0;
  }

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

//********************************************************************************************************
// description: function prints an anisotropic field of type 'double' which is saved as 2D anisotropic
//              array
// last edit:   25. 02. 2010
// -------------------------------------------------------------------------------------------------------
// note:        @field - '@n x @m' matrix
//              @n - number of rows of printed field
//              @m - pointer to a vecotor containig the length (number of entries) of each line
//              @lineBreak - a number after which the line break is expected (formating purposes)
//********************************************************************************************************
void MatrixOperations :: printAnisoField( double ** field, int n, int * m, int lineBreak,
                    const char * notice )
{
  int i = 0, j = 0, k = 0;

  printf( "%s", notice );
  for( i = 0; i < n; i++ ){
    for( j = 0; j < m[i]; j++, k++ ){
      if( k == lineBreak ){
    printf( "\n" );
    k = 0;
      }else ;
      if( ABS( field[i][j] ) > MACHINE_EPS_PRN ){
    ( k == 0 )? printf( "%.6e", field[i][j] ) : printf( " %.6e", field[i][j] );
      }
      else{
    ( k == 0 )? printf( "%.6e", 0. ) : printf( " %.6e", 0. );
      }
    }printf( "\n\n" ); k = 0;
  }

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

//********************************************************************************************************
// description: function prints an anisotropic field ot type 'int' which is saved as 2D anisotropic array
// last edit:   25. 02. 2010
// -------------------------------------------------------------------------------------------------------
// note:        @field - '@n x @m' matrix
//              @n - number of rows of printed field
//              @m - pointer to a vecotor containig the length (number of entries) of each line
//              @lineBreak - a number after which the line break is expected (formating purposes)
//********************************************************************************************************
void MatrixOperations :: printAnisoField( int ** field, int n, int * m, int lineBreak,
                    const char * notice )
{
  int i = 0, j = 0, k = 0;

  printf( "%s", notice );
  for( i = 0; i < n; i++ ){
    for( j = 0; j < m[i]; j++, k++ ){
      if( k == lineBreak ){
    printf( "\n" );
    k = 0;
      }else ;
      if( ABS( field[i][j] ) > MACHINE_EPS_PRN ){
    ( k == 0 )? printf( "%d", field[i][j] ) : printf( " %d", field[i][j] );
      }
      else{
    ( k == 0 )? printf( "%d", 0 ) : printf( " %d", 0 );
      }
    }printf( "\n\n" ); k = 0;
  }

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



void ludcmp(double **a, int n, int *indx, double *d)
{
  int i,imax,j,k;
  double big,dum,temp=0;
  double *vv;

  vv = new double [n];
  *d=1.0;
  for (i=0;i<n;i++) {
    big=0.0;
    for (j=0;j<n;j++) {
      temp=fabs(a[i][j]);
      if (temp > big) big=temp;
    }
    if (big == 0.0) _errorr1( "Singular matrix in routine ludcmp\n");
    vv[i]=1.0/big;
  }
  
  for (j=0;j<n;j++) {
    for (i=0;i<j;i++)
      {
    for (k=0;k<i;k++)
      a[i][j] -= a[i][k]*a[k][j];
      }
    big=0.0;
    for (i=j;i<n;i++)
      {
    for (k=0;k<j;k++)
      a[i][j] -= a[i][k]*a[k][j];
    if((dum=vv[i]*fabs(a[i][j])) >= big)
      {
        big=dum;
        imax=i;
      }
      }
    if (j != imax)
      {
    for (k=0;k<n;k++)
      {
        dum=a[imax][k];
        a[imax][k]=a[j][k];
        a[j][k]=dum;
      }
    *d = -(*d);
    vv[imax]=vv[j];
      }
    indx[j]=imax;
    if (a[j][j] == 0.0)
      a[j][j] = 1.0e-20;
    if (j != n-1)
      {
    dum=1.0/(a[j][j]);
    for (i=j+1;i<n;i++)
      a[i][j] *= dum;
      }
  }
  delete [] vv;
}

void lubksb(double **a, int n, int *indx, double b[])
{
  int i,ii=-1,ip,j;
  double sum;
  
  for (i=0;i<n;i++) {
    ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];
    if (ii>=0.0)
      for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
    else if (sum) ii=i;
    b[i]=sum;
  }
  for (i=n-1;i>-1;i--) {
    sum=b[i];
    for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}


void
MatrixOperations :: giveInverseMatrix (double *mtrx, int n)
{
#ifdef DEBUG
  if (mtrx == NULL)  _errorr( "No matrix for inversion defined" );
  if (n    == 0   )  _errorr( "Zero size of the matrix for inversion" );
#endif

  double **A = NULL;
  AllocateArray2D (A,n,n);
  copyMatrix_1Dto2D (mtrx, A, n);

  giveInverseMatrix_core (&mtrx, A, n, true);

  DeleteArray2D(A, n);
}

void
MatrixOperations :: giveInverseMatrix (double **mtrx, int n)
{
#ifdef DEBUG
  if (mtrx == NULL)  _errorr( "No matrix for inversion defined" );
  if (n    == 0   )  _errorr( "Zero size of the matrix for inversion" );
#endif

  double **A = NULL;
  AllocateArray2D (A,n,n);
  CopyArray2D(mtrx, A, n, n);

  giveInverseMatrix_core (mtrx, A, n, true);

  DeleteArray2D(A, n);
}

void
MatrixOperations :: giveInverseMatrix_core (double **answer, double **mtrx, int n, bool vect)
{
#ifdef DEBUG
  if (answer == NULL)  _errorr ("Field for answer must be allocated");
  if (mtrx   == NULL)  _errorr( "No matrix for inversion defined" );
  if (n      == 0   )  _errorr( "Zero size of the matrix for inversion" );
#endif
  
  int i, j;
  double d;
  
  double *col  = new double [n];
  int    *indx = new int [n];
  
  ludcmp( mtrx, n, indx, &d );
  
  for ( j=0; j<n; j++ ) {
    for ( i=0; i<n; i++ ) col[i] = 0.0;
    col[j] = 1.0;
    lubksb( mtrx, n, indx, col );
    if (vect) for ( i=0; i<n; i++ ) answer[0][i*n + j] = col[i];
    else      for ( i=0; i<n; i++ ) answer   [i   ][j] = col[i];
  }
  
  delete [] col;
  delete [] indx;
}







void MatrixOperations :: CopyTensorAndReduceToFeep (const double *a, double *b)
{
  b[0] = a[0];
  b[1] = a[4];
  b[2] = a[8];
  b[3] = a[1];
  b[4] = a[5];
  b[5] = a[2];
  
  return ;
}//end of function: CopyTensorAndReduceToFeep( )

void MatrixOperations :: convert2Dtensor_TRtoROW_notation (double *t)
{
  t[3] = t[1];
  t[1] = t[2];
}
void MatrixOperations :: convert3Dtensor_TFEEPtoROW_notation (double *t)
{
  t[8] = t[2];
  t[7] = t[4];
  t[2] = t[5];
  t[6] = t[5];
  t[5] = t[7];
  t[4] = t[1];
  t[1] = t[3];
}
void MatrixOperations :: convert2Dtensor_TRtoVR_notation (double *t)
{
  t[2] *= 2.0;
}
void MatrixOperations :: convert3Dtensor_TRtoVR_notation (double *t)
{
  t[3] *= 2.0;
  t[4] *= 2.0;
  t[5] *= 2.0;
}

void MatrixOperations :: convert3Dtensor_9ROWto6FEEP_notation (double *t)
{
  if (t[1] != t[3] || t[2] != t[6] || t[5] != t[7])  _errorr("tensor is not symmetrical");
  t[1] = t[4];
  t[2] = t[8];
  t[4] = t[7];
  t[5] = t[6];
  t[6] = t[7] = t[8] = 0.0;
}
void MatrixOperations :: convert2Dtensor_4ROWto3FEEP_notation (double *t)
{
  if (t[1] != t[2])  _errorr("tensor is not symmetrical");
  t[1] = t[3];
  t[3] = 0.0;
}

void MatrixOperations :: copy3Dto2Dtensors_FEEPreduced (const double *a, double *b)
{
  b[0] = a[0];
  b[1] = a[1];
  b[2] = a[3];
}
void MatrixOperations :: copy2Dto3Dtensors_FEEPreduced (const double *a, double *b)
{
  b[0] = a[0];
  b[1] = a[1];
  b[2] = 0.0;
  b[3] = a[2];
  b[4] = 0.0;
  b[5] = 0.0;
}

void MatrixOperations :: copy3DeshelbyTensor_reduced2full (const double *a, double **b)
{
  CleanArray2d (b, 6, 6);

  b[0][0] = a[0];
  b[0][1] = a[1];
  b[0][2] = a[2];

  b[1][0] = a[3];
  b[1][1] = a[4];
  b[1][2] = a[5];

  b[2][0] = a[6];
  b[2][1] = a[7];
  b[2][2] = a[8];

  b[3][3] = a[9];
  b[4][4] = a[10];
  b[5][5] = a[11];
}
void MatrixOperations :: copy2DeshelbyTensor_reduced2full (const double *a, double **b)
{
  CleanArray2d (b, 3, 3);

  b[0][0] = a[0];
  b[0][1] = a[1];

  b[1][0] = a[2];
  b[1][1] = a[3];

  b[2][2] = a[4];
}



void MatrixOperations :: copy3DeshelbyTensor_full2reduced (const double *a, double *b)
{
  b[0] = a[0];
  b[1] = a[1];
  b[2] = a[2];
  
  b[3] = a[6];
  b[4] = a[7];
  b[5] = a[8];
  
  b[6] = a[12];
  b[7] = a[13];
  b[8] = a[14];
  
  b[ 9] = a[21];
  b[10] = a[28];
  b[11] = a[35];
}

void MatrixOperations :: copy3DeshelbyTensor_reduced2full (const double *a, double *b)
{
  CleanVector(b, 6*6);

  b[ 0] = a[0];
  b[ 1] = a[1];
  b[ 2] = a[2];

  b[ 6] = a[3];
  b[ 7] = a[4];
  b[ 8] = a[5];

  b[12] = a[6];
  b[13] = a[7];
  b[14] = a[8];

  b[21] = a[9];
  b[28] = a[10];
  b[35] = a[11];
}

void MatrixOperations :: copy2DeshelbyTensor_full2reduced (const double *a, double *b)
{
  b[0] = a[0];
  b[1] = a[1];
  b[2] = a[3];
  b[3] = a[4];
  b[4] = a[8];
}
void MatrixOperations :: copy2DeshelbyTensor_reduced2full (const double *a, double *b)
{
  CleanVector(b, 9);
  
  b[0] = a[0];
  b[1] = a[1];
  b[3] = a[2];
  b[4] = a[3];
  b[8] = a[4];
}



// /**
//  * Function prints an Isotropic tensor having only 12 nonzero entries in matrix form.
//  * @param s[12]  Eshelby tensor
//  * @param flag   _MATEMATHICA_SYNTAX_ or _STANDARD_SYNTAX_
//  */
// void MatrixOperations :: printSMatrixReduced (double S[12], const char * notice, unsigned flag)
// {
//   if( flag == _STANDARD_SYNTAX_ ){
//     printf( "%s", notice );
//     printf( "%e %e %e %e %e %e\n", S[0], S[1], S[2], 0., 0., 0. );
//     printf( "%e %e %e %e %e %e\n", S[3], S[4], S[5], 0., 0., 0. );
//     printf( "%e %e %e %e %e %e\n", S[6], S[7], S[8], 0., 0., 0. );
//     printf( "%e %e %e %e %e %e\n", 0., 0., 0., S[9], 0., 0. );
//     printf( "%e %e %e %e %e %e\n", 0., 0., 0., 0., S[10], 0. );
//     printf( "%e %e %e %e %e %e\n", 0., 0., 0., 0., 0., S[11] );
//   }
//   else if( flag == _MATHEMATICA_SYNTAX_ ){//suitable for debuging
//     printf( "%s", notice );
//     printf( "{{%.6lf,%.6lf,%.6lf,%.6lf,%.6lf,%.6lf},", S[0], S[1], S[2], 0., 0., 0. );
//     printf( "{%.6lf,%.6lf,%.6lf,%.6lf,%.6lf,%.6lf},", S[3], S[4], S[5], 0., 0., 0. );
//     printf( "{%.6lf,%.6lf,%.6lf,%.6lf,%.6lf,%.6lf},", S[6], S[7], S[8], 0., 0., 0. );
//     printf( "{%.6lf,%.6lf,%.6lf,%.6lf,%.6lf,%.6lf},", 0., 0., 0., S[9], 0., 0. );
//     printf( "{%.6lf,%.6lf,%.6lf,%.6lf,%.6lf,%.6lf},", 0., 0., 0., 0., S[10], 0. );
//     printf( "{%.6lf,%.6lf,%.6lf,%.6lf,%.6lf,%.6lf}}\n", 0., 0., 0., 0., 0., S[11] );
//   }
//   else{
//     _errorr( "printIsoTensor: unknown syntax flag");
//   }
// 
//   return ;
// }//end of function: printIsoTensor( )




} // end of namespace mumech

/*end of file*/