homogenizationMethods.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 "homogenizationMethods.h"
#include "problem.h"
#include "mesh.h"
#include "mnode.h"
#include "melement.h"
#include "elementaryFunctions.h"
#include "matrixOperations.h"
#include "esuf.h"

#include <algorithm>


namespace mumech {

// ********************************************************************
// ************ DILUTE METHOD *****************************************
// ********************************************************************
void
Dilute :: giveHomogenizedStiffnessMatrix( double* answer )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for stiffness matrix must be allocated" );
#endif
  
  int mSize = this->giveSizeOfReducedMatrix();
  int mSizeFull = this->giveSizeOfFullMatrix()*this->giveSizeOfFullMatrix();
  double *Ai    = new double[mSize];
  double *Cm    = new double[mSize];
  double *CmInv = new double[mSize];
  double *Ci_Cm = new double[mSize];
  double *Aig   = new double[mSizeFull];
  double *dum1  = new double[mSizeFull];

  double ci_sum = 0.0;
  
  ProblemType pT = this->giveProblemType();
  
  P->matrix_giveReducedStiffMatrix(Cm);
  giveInverseOfReducedMatrix(CmInv, Cm, pT);
  
  
  // Dilute homogenization
  giveFullMatrixFromReducedMatrix(answer, Cm, pT); // answer = Cm;
  
  for( int i = 0; i < this->giveNumberOfInclusions(); i++) {
    if (! inside[i])
      continue;

    // concentration factor in local coordinate system
    giveDiluteReducedConcentrationFactorOfInclusionInLC(Ai, Ci_Cm, Cm, CmInv, i);

    // transformation to global coordinate system
    // Ci and Cm are isotropic => no rotation needed
    giveFullMatrixInGCSFromReducedMatrixInLCS(Aig, Ai, i);
    
    giveProductOfReducedMatrixRegularMatrix(dum1, Ci_Cm, Aig, pT, this->giveSizeOfFullMatrix()); // dum1 = Ci_Cm*Ai;
    
    AddVectorMultipledBy(giveVolumeFractionOfInclusion(i), dum1, answer, mSizeFull); // answer += ci*dum1 ( =(ci*(Ci-Cm)*Ai) )

    ci_sum += giveVolumeFractionOfInclusion(i);
  }
  
  if( ci_sum > 1. ) {
    _errorr1( "Volume of inclusions is greater then 1" );
  }

  delete [] Ai;
  delete [] Cm;
  delete [] CmInv;
  delete [] Ci_Cm;
  delete [] Aig;
  delete [] dum1;
}

void
Dilute :: giveDiluteReducedConcentrationFactorOfInclusionInLC( double* answer, double* matrix_inclStiffmat, double* matrixStiffMat, double* matrixComplMat, const long inclusionNumber )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for Eshelby tensor must be allocated" );
  if( matrix_inclStiffmat == NULL ) _errorr( "Field for matrix_inclStiffmat tensor must be allocated" );
#endif
  
  int mSize = this->giveSizeOfReducedMatrix();
  
  double *AiInv = new double[mSize];
  double *Ci    = new double[mSize];
  double *S     = new double[mSize];
  double *dum1  = new double[mSize];
  double *dum2  = new double[mSize];
  
  ProblemType pT = this->giveProblemType();
  
  giveReducedUnitMatrix(AiInv, pT); // AiInv = I
  
  giveReducedEshelbyMatrixOfInclusion(S, inclusionNumber);
  giveReducedStiffnessMatrixOfInclusion(Ci, inclusionNumber);
  
  CopyVector(Ci, matrix_inclStiffmat, mSize); // Ci_Cm = Ci;
  AddVectorMultipledBy(-1.0, matrixStiffMat, matrix_inclStiffmat, mSize); // Ci_Cm -= Cm;
  giveProductOfReducedMatrices(dum1, matrixComplMat, matrix_inclStiffmat, pT); // dum1 = CmInv*Ci_Cm;
  giveProductOfReducedMatrices(dum2, S, dum1, pT); // dum2 = S*dum1;
  
  AddVector(dum2, AiInv, mSize); // AiInv += dum2;
  giveInverseOfReducedMatrix(answer, AiInv, pT); // Ai = AiInv^{-1};
  
  delete [] AiInv;
  delete [] Ci;
  delete [] S;
  delete [] dum1;
  delete [] dum2;
}

void
Dilute :: giveFullDiluteConcentrationFactorOfInclusionInLC( double* answer, double* matrix_inclStiffmat, double* matrixStiffMat, double* matrixComplMat, const long inclusionNumber )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for Eshelby tensor must be allocated" );
  if( matrix_inclStiffmat == NULL ) _errorr( "Field for matrix_inclStiffmat tensor must be allocated" );
#endif

  int mSize = this->giveSizeOfFullMatrix();
  int nElements = mSize * mSize;

  double *AiInv = new double[nElements];
  double *Ci    = new double[nElements];
  double *S     = new double[nElements];
  double *dum1  = new double[nElements];
  double *dum2  = new double[nElements];

  ProblemType pT = this->giveProblemType();

  if( pT != PT_3D ) {
    _errorr( "Unsupported problem type, must be added" );
  }

  giveFullUnitMatrix(AiInv, pT); // AiInv = I

  // Eshelby ttensor of inclusion in anisotropic matrix
  InclusionRecord3D *incl = dynamic_cast <InclusionRecord3D*>(P->give_inclusion(inclusionNumber));
  incl->esuf->giveSijkl(S, incl->a, matrixStiffMat, 15, 40);

  giveFullStiffnessMatrixOfInclusion(Ci, inclusionNumber);

  CopyVector(Ci, matrix_inclStiffmat, nElements); // Ci_Cm = Ci;
  AddVectorMultipledBy(-1.0, matrixStiffMat, matrix_inclStiffmat, nElements); // Ci_Cm -= Cm;
  giveMatrixMatrixProduct(matrixComplMat, matrix_inclStiffmat, dum1, mSize); // dum1 = CmInv*Ci_Cm;
  giveMatrixMatrixProduct(S, dum1, dum2, mSize); // dum2 = S*dum1;

  AddVector(dum2, AiInv, nElements); // AiInv += dum2;

  CopyVector(AiInv, answer, nElements);
  MatrixOperations::giveInverseMatrix(answer, mSize);

  delete [] AiInv;
  delete [] Ci;
  delete [] S;
  delete [] dum1;
  delete [] dum2;
}


// ********************************************************************
// ************ MORI-TANAKA METHOD ************************************
// ********************************************************************
void
MoriTanaka :: giveHomogenizedStiffnessMatrix( double* answer )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for stiffness matrix must be allocated" );
#endif
  
  int mSize = this->giveSizeOfReducedMatrix();
  int mSizeFull = this->giveSizeOfFullMatrix()*this->giveSizeOfFullMatrix();
  double *Ai    = new double[mSize];
  double *Cm    = new double[mSize];
  double *CmInv = new double[mSize];
  double *Ci_Cm = new double[mSize];
  double *Aig   = new double[mSizeFull];
  double *dum1  = new double[mSizeFull];
  double *O     = new double[mSizeFull];
  double *I     = new double[mSizeFull];
  double *U     = new double[mSizeFull];

  double ci_sum = 0.0;

  ProblemType pT = this->giveProblemType();

  P->matrix_giveReducedStiffMatrix(Cm);
  giveInverseOfReducedMatrix(CmInv, Cm, pT);


  // Mori-Tanaka homogenization
  giveFullMatrixFromReducedMatrix(answer, Cm, pT); // answer = Cm;

  CleanVector(O, mSizeFull);
  CleanVector(U, mSizeFull);
  for( int i = 0; i < this->giveNumberOfInclusions(); i++) {
    if (! inside[i])
      continue;

    // concentration factor in local coordinate system
    Dilute::giveDiluteReducedConcentrationFactorOfInclusionInLC(Ai, Ci_Cm, Cm, CmInv, i);

    // transformation to global coordinate system
    // Ci and Cm are isotropic => no rotation needed
    giveFullMatrixInGCSFromReducedMatrixInLCS(Aig, Ai, i);

    AddVectorMultipledBy(giveVolumeFractionOfInclusion(i), Aig, O, mSizeFull);

    giveProductOfReducedMatrixRegularMatrix(dum1, Ci_Cm, Aig, pT, this->giveSizeOfFullMatrix()); // dum1 = Ci_Cm*Ai;
    AddVectorMultipledBy(giveVolumeFractionOfInclusion(i), dum1, U, mSizeFull); // U += ci*dum1 ( =(ci*(Ci-Cm)*Ai) )

    ci_sum += giveVolumeFractionOfInclusion(i);
  }

  if( ci_sum > 1. ) {
    _errorr1( "Volume of inclusions is greater then 1" );
  }

  giveFullUnitMatrix(I, this->giveProblemType());
  AddVectorMultipledBy((1.-ci_sum), I, O, mSizeFull); // O += cm*I

  MatrixOperations :: giveInverseMatrix(O, this->giveSizeOfFullMatrix());
  
  giveMatrixMatrixProduct(U, O, dum1, this->giveSizeOfFullMatrix());
  
  AddVector(dum1, answer, mSizeFull);
  
  delete [] Ai;
  delete [] Cm;
  delete [] CmInv;
  delete [] Ci_Cm;
  delete [] Aig;
  delete [] dum1;
  delete [] O;
  delete [] I;
  delete [] U;
}


// ********************************************************************
// ************ REUSS BOUND *******************************************
// ********************************************************************
void
ReussBound :: giveHomogenizedStiffnessMatrix( double* answer )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for stiffness matrix must be allocated" );
#endif
  
  int mSize = this->giveSizeOfReducedMatrix();
  double *C    = new double[mSize];
  double *CInv = new double[mSize];
  double *answerRedInv = new double[mSize];
  ProblemType pT = this->giveProblemType();

  double ci_sum = 0.0;


  CleanVector(answerRedInv, mSize);
  for( int i = 0; i < this->giveNumberOfInclusions(); i++) {
    if (! inside[i])
      continue;

    // inclusion stiffness matrix
    giveReducedStiffnessMatrixOfInclusion(C, i);

    // transformation to global coordinate system
    // inclusion stiffness matrix is isotropic => no rotation needed
    giveInverseOfReducedMatrix(CInv, C, pT);

    AddVectorMultipledBy(giveVolumeFractionOfInclusion(i), CInv, answerRedInv, mSize);
    ci_sum += giveVolumeFractionOfInclusion(i);
  }

  if( ci_sum > 1. ) {
    _errorr1( "Volume of inclusions is greater then 1" );
  }

  P->matrix_giveReducedStiffMatrix(C);
  giveInverseOfReducedMatrix(CInv, C, pT);

  AddVectorMultipledBy((1.-ci_sum), CInv, answerRedInv, mSize);

  // stiffness matrix from compliance matrix
  giveInverseOfReducedMatrix(C, answerRedInv, pT);

  giveFullMatrixFromReducedMatrix(answer, C, pT);

  delete [] C;
  delete [] CInv;
  delete [] answerRedInv;
}


// ********************************************************************
// ************ VOIGHT BOUND ******************************************
// ********************************************************************
void
VoightBound :: giveHomogenizedStiffnessMatrix( double* answer )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for stiffness matrix must be allocated" );
#endif
  
  int mSize = this->giveSizeOfReducedMatrix();
  double *C    = new double[mSize];
  double *answerRed = new double[mSize];
  ProblemType pT = this->giveProblemType();

  double ci_sum = 0.0;


  CleanVector(answerRed, mSize);
  for( int i = 0; i < this->giveNumberOfInclusions(); i++) {
    if (! inside[i])
      continue;

    // inclusion stiffness matrix
    giveReducedStiffnessMatrixOfInclusion(C, i);

    // transformation to global coordinate system
    // inclusion stiffness matrix is isotropic => no rotation needed

    AddVectorMultipledBy(giveVolumeFractionOfInclusion(i), C, answerRed, mSize);
    ci_sum += giveVolumeFractionOfInclusion(i);
  }

  if( ci_sum > 1. ) {
    _errorr1( "Volume of inclusions is greater then 1" );
  }
  
  P->matrix_giveReducedStiffMatrix(C);
  
  AddVectorMultipledBy((1.-ci_sum), C, answerRed, mSize);

  giveFullMatrixFromReducedMatrix(answer, answerRed, pT);

  delete [] C;
  delete [] answerRed;
}


// ********************************************************************
// ************ REGULAR GRID ******************************************
// ********************************************************************

void
RegGrid :: giveHomogenizedStiffnessMatrix (double *answer)
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for stiffness matrix must be allocated" );
#endif
  
  if (this->mesh) {
    // CHECK
    if (P->give_matrix()->remote_strain_is_unit() == false)  _errorr("Remote strain is not equal to \"unit set\".");
    if (P->is_converted_to_equivalent()           == false)  _errorr("Problem is not converted to equivalent.");
    // tady be se melo zkontrolovat, ze v meshi jsou spocitany vsechny vysledky, tj. vsechny jednotkove stavy
  }
  else {
    // tady se predpoklada, ze nejou zatezovaci stavy, neni prevedeno na equivalent problem
    // to vse je treba dat do cajku, zkontrolovat, ze to je spocitano, popripade spocitat...
    
    // NO MESH, mesh has to be genereated inside of bounding box
    //P->give_matrix()->set_unit_remote_strain();
    //P->input_data_initialize_and_check_consistency();
    //P->convert_to_equivalent_problem ();
    
    long n[3];
    n[0] = n[1] = n[2] = 40;
    if (P->give_twodim())  n[2] = 0;
    
    Mesh *bbmesh = P->give_new_mesh();
    bbmesh->generate_regular_mesh (bb1, bb2, n);
    bbmesh->compute_element_fields ('t', false, true, true, 0, 0, 0, PFCM_FULL);
    mesh = bbmesh;
    bbmesh = NULL;
  }
  
  // debug
  //mesh->print_geometry_file_vtk("homog/bounding_box.vtk",5,1);
  
  int tr = P->give_VM_TENS_RANGE();
    
  //
  double *strain = Allocate1Ddz(tr*tr);
  double *stress = Allocate1Ddz(tr*tr);

  bool rmsh = mesh->mtA == MT_REGULAR;
  double voli, vol = 0.0;

  long i, j, k;
  for (i=0; i<mesh->nElems; i++) {
    if (rmsh) voli = 1.0;
    else      voli = mesh->Elems[i]->give_volume();

    vol += voli;

    // loop over load cases
    for (j=0; j<tr; j++)
      // loop over stress vector
      for (k=0; k<tr; k++) {
        strain[tr*k+j] += voli * mesh->Elems[i]->strain[0][j][k];
        stress[tr*k+j] += voli * mesh->Elems[i]->stress[0][j][k];
      }
  }

  DivideVector (strain, tr*tr, vol);
  DivideVector (stress, tr*tr, vol);
  
  MatrixOperations :: giveInverseMatrix (strain, tr);
  // tady to musim prenasobit inversnim strainem, pac prumerny strain nebude jednotkovy, pac homogenizuju na prilis male oblasti, ne na nekonecnu
  if (1)  giveMatrixMatrixProduct (stress, strain, answer, tr);
  else    CopyVector(stress, answer, tr*tr);
  
  //
  delete [] strain;
  delete [] stress;
  
}


// ********************************************************************
// ************ DIFFERENTIAL SCHEME ***********************************
// ********************************************************************
void
DifferentialScheme :: giveHomogenizedStiffnessMatrix( double* answer )
{
#ifdef DEBUG
  if( answer == NULL )  _errorr( "Field for stiffness matrix must be allocated" );
#endif

  int mSize = this->giveSizeOfFullMatrix();
  int nElements = mSize*mSize;

  double *Ceff       = new double[nElements];
  double *CeffInv    = new double[nElements];
  double *CeffLoc    = new double[nElements];
  double *CeffInvLoc = new double[nElements];
  double *Ai         = new double[nElements];
  double *Ci_Cm      = new double[nElements];
  double *AiLoc      = new double[nElements];
  double *Ci_CmLoc   = new double[nElements];
  double *dum1       = new double[nElements];

  double ci_sum = 0.0;

  // sorting of inclusion based on their volume
  Inclusion **incl = new Inclusion*[this->giveNumberOfInclusions()];
  for( int i = 0; i < this->giveNumberOfInclusions(); i++) {
    incl[i] = this->P->give_inclusion(i);
  }
  std::sort(incl, incl+this->giveNumberOfInclusions(), &ascendingSort);


  P->matrix_giveFullStiffMatrix(Ceff);  // stifness of the matrix (isotropic)
  for( int i = 0; i < this->giveNumberOfInclusions(); i++) {

    // each inclusion is threated separately and effective stiffness is determined when new inclusion is added
    // Dilute scheme is utilized in this task (// Ceff_{i} = Ceff_{i-1} + (ci*(Ci-Cm)*Ai))

    // id of the inclusion - based on the volume
    int ii = incl[i]->id;
    
    CopyVector(Ceff, CeffInv, nElements);
    MatrixOperations::giveInverseMatrix(CeffInv, mSize);
    if( i > 0 ) {
      // rotation to LC of the inclusion
      giveFullMatrixInLCSFromFullMatrixInGCS(CeffLoc, Ceff, ii);
      giveFullMatrixInLCSFromFullMatrixInGCS(CeffInvLoc, CeffInv, ii);
    } else {
      // Initial stiffness matrix is isotropic
      CopyVector(Ceff, CeffLoc, nElements);
      CopyVector(CeffInv, CeffInvLoc, nElements);
    }

    // concentration factor in local coordinate system
    giveFullDiluteConcentrationFactorOfInclusionInLC(AiLoc, Ci_CmLoc, CeffLoc, CeffInvLoc, ii);

    // transformation to global coordinate system
    giveFullMatrixInGCSFromFullMatrixInLCS(Ci_Cm, Ci_CmLoc, ii);
    giveFullMatrixInGCSFromFullMatrixInLCS(Ai, AiLoc, ii);


    giveMatrixMatrixProduct(Ci_Cm, Ai, dum1, mSize); // dum1 = Ci_Cm*Ai;

    AddVectorMultipledBy(giveVolumeFractionOfInclusion(ii), dum1, Ceff, nElements); // answer += ci*dum1 ( =(ci*(Ci-Cm)*Ai) )

    // calculation of total volume
    ci_sum += giveVolumeFractionOfInclusion(ii);
  }

  CopyVector(Ceff, answer, nElements);

  delete [] Ceff;
  delete [] CeffInv;
  delete [] CeffLoc;
  delete [] CeffInvLoc;
  delete [] Ai;
  delete [] Ci_Cm;
  delete [] AiLoc;
  delete [] Ci_CmLoc;
  delete [] dum1;
  delete [] incl;
}

} // end of namespace mumech

/*end of file*/