selfBalanceAlgorithm.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                         ###  ###  #####  ####  ##  ##
//                       ##  ##  ########  ##    ##  ## ##  ##
//                       ##  ##  ## ## ##  ###   ##     ######
//                       ##  ##  ##    ##  ##    ##  ## ##  ##
//                       #####   ##    ##  #####  ####  ##  ##
//                       ##
//
// name:           selfBalanceAlgorithm.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 <string.h>
#include "macros.h"
#include "esuf.h"
#include "esei.h"
#include "eshelbySoluLambda.h"
#include "transformTensors.h"
#include "selfBalanceAlgorithm.h"
#include "matrixOperations.h"
#include "elementaryFunctions.h"
#include "problem.h"


namespace mumech {


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

void selfBalanceAlgorithm :: copyGlobTotalEigStrainsToAuxVector (InclusionRecord3D **inclRec, double *auxVect, int n)
{
  for( int i = 0; i < n; i++ ){
    CopyVector( inclRec[i]->globEigStrainTotal, ( auxVect + i * 6 ), 6 );
  }
  
  return ;
}//end of function: copyGlobTotalEigStrainsToAuxVector( )


void selfBalanceAlgorithm :: totalEigStrainInInclCentroidsUpdate (Problem *problem)
{
  if (problem->SBA_type == SBAT_STANDA) {
    ;
  }
  else if (problem->SBA_type == SBAT_ORIGINAL) {
    // test whether all inclusions are 3D
    InclusionRecord3D **inclRec = new InclusionRecord3D*[problem->noIncl];
    for (int i=0; i<problem->noIncl; i++) {
      inclRec[i] = dynamic_cast<InclusionRecord3D*>(problem->inclusions[i]);
      if (inclRec[i] == NULL) _errorr("");
    }
    
    // info print
    if (problem->verbose > 0)  printf( "Transformation strain update in progress, please wait ... \n" );
    
    Point *point = NULL;
    double * auxStrains = NULL;
    double norm;
    int step, i, j;
    
    // auxiliary 'Point' memory allocation
    point = new Point[problem->noIncl];
    for( i = 0; i < problem->noIncl; i++ ){
      point[i].set_nlc(problem, 1);
      CopyVector( inclRec[i]->origin, point[i].x, 3 );
    }
    
    // auxiliary memory allocation to save the total global eigenstrains of previous 'balance' step
    auxStrains = new double[problem->noIncl * 6];
    long noActIncls, indx;
    
    // Self-balance iterative algorithm.
    step = 0;
    do {
      // maximum number of interations reached - non convergence
      if (step >= problem->get_SBA_maxiters())  _errorr("convergence not reached");
      // required number of interations reached - break
      if (step >= problem->get_SBA_reqiters())  break;
      step++;
      
      //update of auxGlobEigenstrain
      selfBalanceAlgorithm::copyGlobTotalEigStrainsToAuxVector( inclRec, auxStrains, problem->noIncl );
      
      for (i = 0; i < problem->noIncl; i++)
    CleanVector(inclRec[i]->locEigStrainTot, 6);
      
      for (i = 0; i < problem->noIncl; i++) {
    // number of active inclusions
        if (problem->SBA_optimized == false) noActIncls = problem->noIncl;
        else                                 noActIncls = inclRec[i]->noActingIncls;
    
    for (j=0; j<noActIncls; j++) {
          if (problem->SBA_optimized == false) indx = j;
          else                                 indx = inclRec[i]->actingIncls[j];
      
      if (indx != i) {
        inclRec[i]->esuf->giveEshelbyStrainOfOnePoint( &point[indx]);
        inclRec[indx]->SBA_updateGlobAndLocStrains(point+indx);
      }
      
    }
      }
      
      for (i=0; i<problem->noIncl; i++)
    inclRec[i]->addtot();
      
      //quadratic norm calculation
      norm = selfBalanceAlgorithm::giveQuadNormMultEigstrain( inclRec, auxStrains, problem->noIncl );
      
      //
      if (problem->verbose > 0)
    printf( "self-balance norm[%d] = %.6e\n", ( step ), norm );
      
    } while( norm > _SELF_BALANCE_NORM_LIMIT_ );
    
    if (problem->verbose > 0)
      printf( "steps: %d\n", step );
    
    // Final local eigenstrains update
    selfBalanceAlgorithm::localEigStrainUpdateTotal( inclRec, problem->noIncl );
    
    // DEBUG PRINT
    if (problem && problem->give_verbose() > 2) {
      //Total global transformation strain printing
      //selfBalanceAlgorithm::printTotEigenStrains( inclRec, problem->noIncl, "Updated transformation strains:\n" );
    }
    
    //auxiliary memory de-allocation
    delete [] point;
    delete [] auxStrains;
    
    delete [] inclRec;
  }
  else if (problem->SBA_type== SBAT_VOID) {
    ;
  }
  else _errorr("Unknown SBAmode");
  
  return ;
}//end of function: totalEigStrainInInclCentroidsUpdate( )


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


// //********************************************************************************************************
// // description: Function copies strain perturbatios in inclusion origins into an auxiliary vector
// // last edit:   02. 03. 2010
// // -------------------------------------------------------------------------------------------------------
// //              @point - point data record, defined in: types.h
// //              @auxVect - auxiliary vector
// //              @n - total number of inclusions
// //********************************************************************************************************
// void selfBalanceAlgorithm::copyPertStrainsToAuxVect( Point * point, double * auxVect, int n )
// {
//   for( int i = 0; i < n; i++ ){
//     CopyVector( point[i].strain[0], ( auxVect + i * 6 ), 6 );
//   }
// 
//   return ;
// }//end of function: copyPertStrainsToAuxVect( )

double selfBalanceAlgorithm :: giveQuadNormMultEigstrain (InclusionRecord3D **inclRec, double * oldStrain, int n)
{
  double norm = 0.;

  // calculation of quadratic norm of eigenstrains in global coordinates
  while (n--)
    norm += give_quadNormTwoVectors ( oldStrain + n * 6, inclRec[n]->globEigStrainTotal, 6);
  
  return norm;
}

//********************************************************************************************************
// description: Function updates total local eigenstrains from their total global counterparts
// last edit:   02. 03. 2010
// -------------------------------------------------------------------------------------------------------
//              @inclRec - inclusion record (data structure of given inclusions: types.h)
//              @n - total number of inclusions
//********************************************************************************************************
void selfBalanceAlgorithm::localEigStrainUpdateTotal( InclusionRecord3D **inclRec, int n )
{
  for (int i=0; i<n; i++)
    inclRec[i]->rotateStrain_G2L(inclRec[i]->globEigStrainTotal, inclRec[i]->locEigStrain);
}

//********************************************************************************************************
// description: Function prints resulting total transformation strains in inclusion centroids
// last edit:   02. 03. 2010
// -------------------------------------------------------------------------------------------------------
//              @inclRec - inclusion record (data structure of given inclusions: types.h)
//              @n - total number of inclusions
//********************************************************************************************************
void selfBalanceAlgorithm::printTotEigenStrains( InclusionRecord3D **inclRec, int n, const char * notice )
{
  char auxName[200];
  
  printf( "%s", notice );
  
  for( int i = 0; i < n; i++ ){
    sprintf( auxName, "Inclusion %d:\n", ( i + 1 ) );
    printVectorRowForm( inclRec[i]->globEigStrainTotal, 6, ( const char * ) auxName );
    //CleanString( auxName, 200 );
  }
}//end of function: printTotEigenStrains( )


} // end of namespace mumech

/*end of file*/