homogenization.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           homogenization.cpp
// author(s):      Jan Vorel, 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.
//********************************************************************************************************

// toto je kvuli WIN32_LEAN_AND_MEAN
#include "types.h"

#include "homogenization.h"
#include "problem.h"
#include "mesh.h"
#include "mnode.h"
#include "melement.h"
#include "elementaryFunctions.h"
#include "matrixOperations.h"


namespace mumech {


Homogenization :: Homogenization (long i, Problem *p)
{
  id = i;
  P = p;     if (p == NULL)  _errorr("Problem not defined");
  volume = -1;
  inside = new bool[p->noIncl];
  for (long i=0; i<p->noIncl; i++)  inside[i] = true;
}


Homogenization :: ~Homogenization ()
{
  P = NULL;    // do not deleta, this is only pointer to problem allocated upper
  delete [] inside;
}


void
Homogenization :: set_boundingBox (double x1, double y1, double x2, double y2)
{
  P->check_dim(2);
  bb1[0] = x1;  bb1[1] = y1;  bb1[2] = 0.0;
  bb2[0] = x2;  bb2[1] = y2;  bb2[2] = 0.0;
}
void
Homogenization :: set_boundingBox (double x1, double y1, double z1, double x2, double y2, double z2)
{
  P->check_dim(3);
  bb1[0] = x1;  bb1[1] = y1;  bb1[2] = z1;
  bb2[0] = x2;  bb2[1] = y2;  bb2[2] = z2;
}
void
Homogenization :: set_boundingBox (const double *p1, const double *p2)
{
  bb1[0] = p1[0];  bb1[1] = p1[1];  bb1[2] = p1[2];
  bb2[0] = p2[0];  bb2[1] = p2[1];  bb2[2] = p2[2];
}

void
Homogenization :: find_inclusions_in_BB (void)
{
  for (int i=0; i<P->noIncl; i++)
    inside[i] = P->give_inclusion(i)->is_inside_of_BB(bb1, bb2);
}


// Functions to approach Problem Properties
int
Homogenization :: giveNumberOfInclusions() const
{
  return P->noIncl;
}

ProblemType Homogenization :: giveProblemType() const
{
  if( P->give_twodim() )  return PT_2D;
  else                    return PT_3D;
}

int
Homogenization :: giveSizeOfSymVector()
{
  return P->give_VM_TENS_RANGE();
}

int
Homogenization :: giveSizeOfReducedMatrix()
{
  return P->give_ISO_C_RANGE();
}

int
Homogenization :: giveSizeOfFullMatrix()
{
  return P->give_VM_TENS_RANGE();
}

double
Homogenization :: giveTotalVolume (void) const
{
  if (volume < 0) {
    if (P->give_twodim())  ((Homogenization*)this)->volume = (bb2[0]-bb1[0])*(bb2[1]-bb1[1]);
    else                   ((Homogenization*)this)->volume = (bb2[0]-bb1[0])*(bb2[1]-bb1[1])*(bb2[2]-bb1[2]);
  }

  return volume;
}




// Functions to approach Inclusion Properties
void
Homogenization :: giveFullEshelbyMatrixOfInclusion(double** answer, const long inclusionNumber )
{
  P->give_inclusion(inclusionNumber)->give_EshelbyMatrixFull(answer);
}

void
Homogenization :: giveReducedEshelbyMatrixOfInclusion(double* answer, const long inclusionNumber )
{
  P->give_inclusion(inclusionNumber)->give_EshelbyMatrixReduced(answer);
}

void
Homogenization :: giveFullStiffnessMatrixOfInclusion( double* answer, const long inclusionNumber )
{
  P->give_inclusion(inclusionNumber)->give_StiffnessMatrixFull(answer);
}

void
Homogenization :: giveReducedStiffnessMatrixOfInclusion( double* answer, const long inclusionNumber )
{
  P->give_inclusion(inclusionNumber)->give_StiffnessMatrixReduced(answer);
}

double
Homogenization :: giveVolumeFractionOfInclusion( const long inclusionNumber )
{
  Inclusion *incl = this->P->give_inclusion(inclusionNumber);
  return incl->give_volume() / this->giveTotalVolume();
}

double
Homogenization :: giveTotalVolumeFractionOfInclusions (void) const
{
  double tf = 0.0;
  for (long i=0; i<P->noIncl; i++)
    if (inside[i])
      tf += P->inclusions[i]->give_volume();

  return tf  / this->giveTotalVolume();
}

void
Homogenization :: giveTransformationMatrixStressStrainG2L (double* answer, const long inclusionNumber)
{
  P->give_inclusion(inclusionNumber)->give_TeMatrix_G2L(answer);
}

void
Homogenization :: giveTransformationMatrixStressStrainL2G( double* answer, const long inclusionNumber )
{
  P->give_inclusion(inclusionNumber)->give_TeMatrix_L2G(answer);
}

void
Homogenization :: giveFullMatrixInGCSFromReducedMatrixInLCS(double* answer, double* ALoc, const long inclusionNumber )
{
  // answer = full matrix stored row by row
  // ALoc = full matrix in local coordinate system

  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  }

  int mSize = this->giveSizeOfFullMatrix();
  double *dumM = new double[mSize*mSize];

  giveFullMatrixFromReducedMatrix(dumM, ALoc, this->giveProblemType());
  giveFullMatrixInGCSFromFullMatrixInLCS(answer, dumM, inclusionNumber);

  delete [] dumM;
}

void
Homogenization :: giveFullMatrixInGCSFromFullMatrixInLCS(double* answer, double* ALoc, const long inclusionNumber )
{
  // answer = full matrix stored row by row
  // ALoc = full matrix in local coordinate system

  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  }

  int mSize = this->giveSizeOfFullMatrix();
  double *dumM = new double[mSize*mSize];
  double *dumT = new double[mSize*mSize];
  
  giveTransformationMatrixStressStrainL2G(dumT, inclusionNumber);

  giveMatrixMatrixProduct(dumT, ALoc, dumM, mSize);

  giveTransformationMatrixStressStrainG2L(dumT, inclusionNumber);

  giveMatrixMatrixProduct(dumM, dumT, answer, mSize);

  delete [] dumM;
  delete [] dumT;
}

void
Homogenization :: giveFullMatrixInLCSFromFullMatrixInGCS(double* answer, double* AGlob, const long inclusionNumber )
{
  // answer = full matrix stored row by row
  // AGlob = full matrix in global coordinate system

  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  }

  int mSize = this->giveSizeOfFullMatrix();
  double *dumM = new double[mSize*mSize];
  double *dumT = new double[mSize*mSize];

  giveTransformationMatrixStressStrainG2L(dumT, inclusionNumber);

  giveMatrixMatrixProduct(dumT, AGlob, dumM, mSize);

  giveTransformationMatrixStressStrainL2G(dumT, inclusionNumber);

  giveMatrixMatrixProduct(dumM, dumT, answer, mSize);

  delete [] dumM;
  delete [] dumT;
}


// General Functions
void
Homogenization :: giveInverseOfReducedMatrix( double* answer, const double* rM, ProblemType pT )
{
  // rM = reduced matrix
  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  }

  switch (pT) {
  case PT_2D: giveInverseMatrix3x3to5 (answer, rM);   break;
  case PT_3D: giveInverseMatrix6x6to12(answer, rM);   break;
  default:
    _errorr( "Unsupported problem type" );
  }
}

void
Homogenization :: giveReducedUnitMatrix( double* answer, ProblemType pT )
{
  if( answer == NULL )  _errorr( "Field for answer must be allocated" );

  switch( pT ) {
  case PT_2D:  MatrixOperations::giveUnitSMatrixReduced_2d (answer);  break;
  case PT_3D:  MatrixOperations::giveUnitSMatrixReduced_3d (answer);  break;
  default: _errorr( "Unsupported problem type" );
  }
}

void
Homogenization :: giveFullUnitMatrix( double* answer, ProblemType pT )
{
  if( answer == NULL )  _errorr( "Field for answer must be allocated" );

  switch( pT ) {
  case PT_2D:  MatrixOperations::giveUnitSMatrixFull_2d (answer);  break;
  case PT_3D:  MatrixOperations::giveUnitSMatrixFull_3d (answer);  break;
  default: _errorr( "Unsupported problem type" );
  }
}

void
Homogenization :: giveProductOfReducedMatrices( double* answer, double* A, double* B, ProblemType pT )
{
  // answer = A*B;
  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  }

  CleanVector(answer, giveSizeOfReducedMatrix());

  switch( pT )
  {
  case PT_2D: {
    answer[0] = A[0]*B[0] + A[1]*B[2];
    answer[1] = A[0]*B[1] + A[1]*B[3];
    answer[2] = A[2]*B[0] + A[3]*B[2];
    answer[3] = A[2]*B[1] + A[3]*B[3];
    answer[4] = A[4]*B[4];

    break;
  }
  case PT_3D: {
    answer[0] = A[0]*B[0] + A[1]*B[3] + A[2]*B[6];
    answer[1] = A[0]*B[1] + A[1]*B[4] + A[2]*B[7];
    answer[2] = A[0]*B[2] + A[1]*B[5] + A[2]*B[8];
    answer[3] = A[3]*B[0] + A[4]*B[3] + A[5]*B[6];
    answer[4] = A[3]*B[1] + A[4]*B[4] + A[5]*B[7];
    answer[5] = A[3]*B[2] + A[4]*B[5] + A[5]*B[8];
    answer[6] = A[6]*B[0] + A[7]*B[3] + A[8]*B[6];
    answer[7] = A[6]*B[1] + A[7]*B[4] + A[8]*B[7];
    answer[8] = A[6]*B[2] + A[7]*B[5] + A[8]*B[8];
    answer[9] = A[9]*B[9];
    answer[10] = A[10]*B[10];
    answer[11] = A[11]*B[11];

    break;
  }
  default:
    _errorr( "Unsupported problem type" );
  }
}

void
Homogenization :: giveFullMatrixFromReducedMatrix( double* answer, const double* A, ProblemType pT )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for answer must be allocated" );
#endif
  
  CleanVector(answer, giveSizeOfFullMatrix()*giveSizeOfFullMatrix());
  
  switch( pT )
  {
  case PT_2D: {
    MatrixOperations::copy2DeshelbyTensor_reduced2full (A, answer);
    break;
  }
  case PT_3D: {
    MatrixOperations::copy3DeshelbyTensor_reduced2full (A, answer);
    break;
  }
  default:
    _errorr( "Unsupported problem type" );
  }
}

void
Homogenization :: giveProductOfRegularMatrixReducedMatrix( double* answer, double* A, double* redB , int n, ProblemType pT )
{
  // answer = A[n,k]*B[k,l];
  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  } else if( giveSizeOfFullMatrix() != n ) {
    _errorr( "Matrix mismatch" );
  }
  
  double *dumB = new double[n*n];
  
  giveFullMatrixFromReducedMatrix(dumB, redB, pT);
  
  giveMatrixMatrixProduct(A, dumB, answer, n);
  
  delete [] dumB;
}

void
Homogenization :: giveProductOfReducedMatrixRegularMatrix( double* answer, double* redA, double* B, ProblemType pT, int n )
{
  // answer = A[n,k]*B[k,l];
  if( answer == NULL ) {
    _errorr( "Field for answer must be allocated" );
  } else if( giveSizeOfFullMatrix() != n ) {
    _errorr( "Matrix mismatch" );
  }
  
  double *dumA = new double[n*n];
  
  giveFullMatrixFromReducedMatrix(dumA, redA, pT);
  
  giveMatrixMatrixProduct(dumA, B, answer, n);
  
  delete [] dumA;
}

} // end of namespace mumech

/*end of file*/