inclusion.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           inclusion.cpp
// author(s):      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 "inclusion.h"

#include "macros.h"
#include "matrixOperations.h"
#include "problem.h"

#include "esuf.h"
#include "esuf_Sphere.h"
#include "esuf_OblateSpheroid.h"
#include "esuf_ProlateSpheroid.h"
#include "esuf_Penny.h"
#include "esuf_FlatEllipsoid.h"
#include "esuf_EllipticCylinder.h"
#include "esuf_Cylinder.h"

#include "esei.h"
#include "esei_Ellipsoid.h"
#include "esei_Sphere.h"
#include "esei_OblateSpheroid.h"
#include "esei_ProlateSpheroid.h"
#include "esei_Penny.h"
#include "esei_FlatEllipsoid.h"
#include "esei_EllipticCylinder.h"
#include "esei_Cylinder.h"

#include "transformTensors.h"
#include "stiffness.h"
#include "mesoface.h"

#include "matrix_vector.h"
#include "gelib.h"
#include "tixy2.h"


namespace mumech {

//using namespace transformTensor;

Inclusion :: Inclusion (long i, const Problem *p)
{
  //
  // INPUT DATA
  //
  id = i;
  this->P = p;
  
  shape = IS_VOID;
  E = nu = 0.0;
  
  origin = new double[3];  // 3 coordinates also for 2d in order to print 3 coordinates to vtk file
  
  int vr = p->give_VECT_RANGE();     // vector range, 2 or 3
  int tr = p->give_VM_TENS_RANGE();  // reduced vector range, 3 or 6
  int ir = p->give_ISO_C_RANGE();    // reduced eshelby tensor, 5 or 12

  a       = new double[vr];
  eAngles = new double[P->give_EA_RANGE()];
  //Remote_strains = NULL;
  
  //
  // PRE-COMPUTED DATA, when AF->char = true
  //
  actionRadius = SQRactionRadius = 0.0;
  
  noActingIncls = 0;
  actingIncls = NULL;
  
  Epsilon_Tau = NULL;
  
  n_approx_points=0;
  approx_points = NULL;
  approx_points_eps_tau = NULL;
  approx_coef = NULL;
  
  globPert_displc = globPert_strain = globPert_stress = NULL;
  
  rotated = true;
  volume = -1.0;
  
  C    = new double[ir];
  S    = new double[ir];
  SInv = new double[ir];

  //
  T    = new double[vr*vr];
  TInv = new double[vr*vr];
  Te    = new double[tr*tr];
  TeInv = new double[tr*tr];
}

Inclusion :: ~Inclusion ()
{
  delete [] origin;
  delete [] a;
  delete [] eAngles;
  
  // if (Remote_strains != NULL) {
  //   for (int i=0; i<P->give_nlc(); i++)
  //     delete [] Remote_strains[i];
  //   
  //   delete [] Remote_strains;
  // }
  
  delete [] actingIncls;
  
  if (C)    delete [] C;
  if (S)    delete [] S;
  if (SInv) delete [] SInv;

  if (T)    delete [] T;
  if (TInv) delete [] TInv;
  if (Te)    delete [] Te;
  if (TeInv) delete [] TeInv;
  
  DeleteArray2D (Epsilon_Tau, P->give_nLC());
  
  DeleteArray2D (globPert_displc, P->give_nLC());
  DeleteArray2D (globPert_strain, P->give_nLC());
  DeleteArray2D (globPert_stress, P->give_nLC());
}

bool Inclusion :: scan_eAngles_RAD (Stream *stream)
{
  return ST_scan_array (stream, P->give_EA_RANGE(), this->eAngles);
}
bool Inclusion :: scan_eAngles_DEG (Stream *stream)
{
  if (this->scan_eAngles_RAD(stream) == false)  return false;
  
  this->EAdeg2rad();
  return true;
}

void Inclusion :: allocate_nlc_fields (void)
{
  int nlc = P->give_nLC();
  if (nlc==0) _errorr("zero nlc");
  
  // if (Remote_strains != NULL)  _errorr("Remote_strains allocated already");
  // Remote_strains = new double*[nlc];
  // for (int i=0; i<nlc; i++)
  //   Remote_strains[i] = new double[P->give_TENS_RANGE()];

  if (globPert_displc == NULL)
    AllocateArray2D(globPert_displc, P->give_nLC(), P->give_VECT_RANGE());
  
  if (Epsilon_Tau == NULL) {
    Epsilon_Tau = AllocateArray2D (P->give_nLC(), P->give_VM_TENS_RANGE());
    // termitovo pokud je to mozne, tak to cleanovat az na miste prirazeni
    CleanArray2d (Epsilon_Tau, P->give_nLC(), P->give_VM_TENS_RANGE());
  }
}

void Inclusion :: ActingIncls_allocate (void)
{
  if (actingIncls) _errorr("nonzero actingIncls");
  
  if (noActingIncls == 0) return;
  
  actingIncls = new int[noActingIncls];
}

int Inclusion :: find_overlap (const Inclusion *incl) const
{
  double epsilon = 1.0e-10;

  // intersection of circumsribed spheres
  double axes = this->a[0] + incl->a[0];
  double dist;
  if (P->twodim) dist = DIST_POINTS_SQR_2D( this->origin, incl->origin );
  else           dist = DIST_POINTS_SQR_3D( this->origin, incl->origin );

  if (dist > axes * axes + epsilon)  return false;

  //
  if (P->twodim) { if (this->shape == IS_CIRCLE  &&  incl->shape == IS_CIRCLE)  return true; }
  else           { if (this->shape == IS_SPHERE  &&  incl->shape == IS_SPHERE)  return true; }

  if (P->twodim) return MesoFace::ellipses_overlap   ( this, incl );
  else           return MesoFace::ellipsoids_overlap ( this, incl );

  return 0;
}

bool Inclusion :: is_inside_of_BB (const double *bb1, const double *bb2) const
{
  if (this->shape != IS_CIRCLE  &&  this->shape != IS_SPHERE)  errol;

  double aa = (1 - 1e-4) * a[0];

  ;                        if ( (origin[0] - aa) < bb1[0]  ||  bb2[0] < (origin[0] + aa) ) return false;
  ;                        if ( (origin[1] - aa) < bb1[1]  ||  bb2[1] < (origin[1] + aa) ) return false;
  if (P->twodim == false)  if ( (origin[2] - aa) < bb1[2]  ||  bb2[2] < (origin[2] + aa) ) return false;

  return true;
}


//
// DIRECT API
//
void Inclusion :: set_centroids (double x, double y)           { P->check_dim(2); origin[0] = x; origin[1] = y; origin[2] = 0.0; } // 3. coordinate in order to print 3 coordinates to vtk file
void Inclusion :: set_centroids (double x, double y, double z) { P->check_dim(3); origin[0] = x; origin[1] = y; origin[2] = z; }
void Inclusion :: set_Semiaxes_dimensions (double x, double y)           { P->check_dim(2); a[0] = x; a[1] = y; }
void Inclusion :: set_Semiaxes_dimensions (double x, double y, double z) { P->check_dim(3); a[0] = x; a[1] = y; a[2] = z; }
void Inclusion :: set_Euller_angles_deg (double x)                     { P->check_dim(2); eAngles[0] = x;                                 this->EAdeg2rad(); }
void Inclusion :: set_Euller_angles_deg (double x, double y, double z) { P->check_dim(3); eAngles[0] = x; eAngles[1] = y; eAngles[2] = z; this->EAdeg2rad(); }
//
void Inclusion :: set_all_2d (double x, double y, double e, double n, double a1, double a2, double e1)
{
  this->set_centroids(x, y);
  this->set_Youngs_modulus(e);
  this->set_Poissons_ratio(n);
  this->set_Semiaxes_dimensions(a1, a2);
  this->set_Euller_angles_deg  (e1);
}
void Inclusion :: set_all_3d (double x, double y, double z, double e, double n,
                              double a1, double a2, double a3, double e1, double e2, double e3 )
{
  this->set_centroids(x, y, z);
  this->set_Youngs_modulus(e);
  this->set_Poissons_ratio(n);
  this->set_Semiaxes_dimensions(a1, a2, a3);
  this->set_Euller_angles_deg  (e1, e2, e3);
}



//void InclusionRecord :: set_Remote_strains_for_lc (int lc, const double *r1)
//{
//  CopyVector(r1, Remote_strains[lc], P->give_TENS_RANGE());
//}

void Inclusion :: input_data_initialize_and_check_consistency (void)
{
  if (P->is_converted_to_equivalent() == false) {
    // conversion infinity axes given in input file by -1.0 to INFTY
    for (int i=0; i<P->give_VECT_RANGE(); i++)
      if (a[i] < 0.0) {
    if (a[i] == -1.0)  a[i] = INFTY;
    else _errorr("negative axis is not equal to -1.0");
      }
    
    // detection of inclusion shape
    shape = this->giveInclusionGeometry(shape);
  }
  else {
    // conversion infinity axes given in input file by -1.0 to INFTY
    for (int i=0; i<P->give_VECT_RANGE(); i++)
      if (a[i] < 0.0)
    _errorr("negative axis");
    
    // detection of inclusion shape
    if (shape == IS_VOID)  _errorr("error");
    
    //initialize SQR action radii
    this->SQRactionRadius = SQR( this->actionRadius );

    TransformTensors::give_TInv (this->TInv,  this->T,  this->shape);
  }
  
  //
  rotated = true;
  if (this->shape == IS_SPHERE  ||  this->shape == IS_CIRCLE)
    rotated = false;
  else {
    if (eAngles[0] == 0.0)  rotated = false;
    if (P->twodim == false)
      if (eAngles[1] != 0.0  ||  eAngles[2] != 0.0)
        rotated = true;
  }

  // set the size of the derivative steps

  // standovo nastaveni
  // this->ndiff_1 = (this->a[0] + this->a[1] + this->a[2])/3.*0.005;
  // this->ndiff_2 = this->ndiff_1;
  // nebo
  // _NDIFF_H_MULTIPLIER_ = 0.008
  // if (this->a[2]>0) this->ndiff_1 = this->a[2] * _NDIFF_H_MULTIPLIER_;
  // else              this->ndiff_1 = this->a[1] * _NDIFF_H_MULTIPLIER_;

  // puvodni hodnoty
  this->ndiff_1 = 1.0e-4;
  this->ndiff_2 = this->ndiff_1;

  // allocate
  this->allocate_nlc_fields();
}

bool Inclusion :: point_is_inside (const double *coords, double epsilon) const
{
  double l;
  double loc[3];
  
  // GLOBAL->LOCAL transform of the point coordinates
  this->transformCoords_G2L(coords, loc);
  
  switch (this->shape) {
  case IS_ELLIPSOID :         
  case IS_SPHERE :
  //case IS_PENNY :
  //case IS_FLAT_ELLIPSOID :
  case IS_OBLATE_SPHEROID :   
  case IS_PROLATE_SPHEROID :  l = SQR( loc[0]/a[0] ) + SQR( loc[1]/a[1] ) + SQR( loc[2]/a[2] );  break;
  //case IS_ELLIPTIC_CYLINDER :
  //case IS_CYLINDER :          l =                      SQR( loc[1]/a[1] ) + SQR( loc[2]/a[2] );  break;
    
  case IS_ELLIPSE : 
  case IS_CIRCLE :            l = SQR( loc[0]/a[0] ) + SQR( loc[1]/a[1] )                     ;  break;
  
  default:
    _errorr( "InclusionRecord :: point_is_inside: unknown inclusion shape\n");
  }
  
  return  l <=  SQR(1 + epsilon);
}//end of function: point_is_inside ()

void Inclusion :: transformCoords_G2L (const double *glob, double *loc) const
{
  //swap of coordinate systems
  SubtractTwoVectors (glob, this->origin, loc, P->give_VECT_RANGE());

  //rotating of the coordinate system
  if (rotated)
    giveMatrixVectorProduct (T, loc, loc, P->give_VECT_RANGE());
}

void Inclusion :: transformCoords_L2G (const double *loc, double *glob) const
{
  // termitovo tato fce je docasne zbastlena pro ucely approx funkci
  
  if (rotated) errol;
  //MatrixOperations::giveMatrixVectorProduct(TInv, loc, glob, P->give_VECT_RANGE());
  
  if (loc == glob) {
    // termitovo switchnout na :
    // AddVector(this->origin, glob, P->give_VECT_RANGE());
    for (int i = 0; i < P->give_VECT_RANGE(); i++)
      glob[i] += this->origin[i];
  }
  else
    errol;
}

void Inclusion :: rotateDisplc_L2G (const double *loc, double *glob) const
{
  if (rotated)
    giveMatrixVectorProduct (this->TInv, loc, glob, P->give_VECT_RANGE());
  else
    if (glob != loc) CopyVector (loc, glob, P->give_VECT_RANGE());
}


void Inclusion :: rotateStrain_G2L (const double *glob, double *loc) const
{
  if (rotated)
    giveMatrixVectorProduct (this->Te, glob, loc, P->give_VM_TENS_RANGE());
  else
    if (glob != loc)  CopyVector (glob, loc, P->give_VM_TENS_RANGE());
}
void Inclusion :: rotateStrain_L2G (const double *loc, double *glob) const
{
  if (rotated)
    giveMatrixVectorProduct (this->TeInv, loc, glob, P->give_VM_TENS_RANGE() );
  else
    if (glob != loc) CopyVector (loc, glob, P->give_VM_TENS_RANGE());
}




double Inclusion :: give_volume( void )
{
  if( volume < 0.0 )
    this->computeVolume();

  return this->volume;
}

void Inclusion::update_approximations(int strainID)
{
  if (P->approximation == 1 && this->n_approx_points) {
    for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)
      this->Epsilon_Tau[strainID][j] = 0.;
    for (int i = 0; i < this->n_approx_points; i++)
      for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)
    this->Epsilon_Tau[strainID][j] += this->approx_points_eps_tau[i][j];
    for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)
      this->Epsilon_Tau[strainID][j] /= this->n_approx_points;
    
    
    for (int i = 0; i < P->give_dimension(); i++) {
      for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)
    this->approx_coef[strainID][i][j] = 1. / (this->a[i] * P->approx_point_pos1*2.)*(this->approx_points_eps_tau[2 * i + 2][j] - this->approx_points_eps_tau[2 * i + 1][j]);
    }
    /*for (int i = 0; i < P->give_dimension(); i++) {
      for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)
      printf("%lf ",this->approx_coef[strainID][i][j]);
      printf("\n");
      }
      enter();*/
  }
  else if (P->approximation == 2 && this->n_approx_points) {
    if (P->give_dimension() == 2) {
      vektor R(6);
      matice M(6, 6);
      double x, y, x2, y2, x3, y3, x4, y4;
      for (int i = 0; i < P->give_VM_TENS_RANGE(); i++) {
    R.vynulovat();
    M.vynulovat();
    /*
      R=[sum(z.*1);sum(z.*x);sum(z.*y);sum(z.*x.*x);sum(z.*y.*y);sum(z.*x.*y)];

      M=[     9           sum(x)          sum(y)          sum(x.*x)       sum(y.*y)       sum(x.*y);
      sum(x)      sum(x.*x)       sum(x.*y)       sum(x.*x.*x)    sum(x.*y.*y)    sum(x.*x.*y);
      sum(y)      sum(x.*y)       sum(y.*y)       sum(x.*x.*y)    sum(y.*y.*y)    sum(x.*y.*y);
      sum(x.*x)   sum(x.*x.*x)    sum(x.*x.*y)    sum(x.*x.*x.*x) sum(x.*x.*y.*y) sum(x.*x.*x.*y);
      sum(y.*y)   sum(x.*y.*y)    sum(y.*y.*y)    sum(x.*x.*y.*y) sum(y.*y.*y.*y) sum(x.*y.*y.*y);
      sum(x.*y)   sum(x.*x.*y)    sum(x.*y.*y)    sum(x.*x.*x.*y) sum(x.*y.*y.*y) sum(x.*x.*y.*y)
      ];    */
    for (int j = 0; j < this->n_approx_points; j++) {
      x = this->approx_points[j][0] - this->origin[0];x2 = x*x;x3 = x*x*x;x4 = x*x*x*x;
      y = this->approx_points[j][1] - this->origin[1];y2 = y*y;y3 = y*y*y;y4 = y*y*y*y;
      R.v[0] += this->approx_points_eps_tau[j][i];
      R.v[1] += this->approx_points_eps_tau[j][i] * x;
      R.v[2] += this->approx_points_eps_tau[j][i] * y;
      R.v[3] += this->approx_points_eps_tau[j][i] * x * x;
      R.v[4] += this->approx_points_eps_tau[j][i] * y * y;
      R.v[5] += this->approx_points_eps_tau[j][i] * x * y;
                                                                                
      M.g(0, 0) += 1.;
      M.g(0, 1) += x;
      M.g(0, 2) += y;
      M.g(0, 3) += x2;
      M.g(0, 4) += y2;
      M.g(0, 5) += x*y;

      M.g(1, 0) += x;
      M.g(1, 1) += x2;
      M.g(1, 2) += x*y;
      M.g(1, 3) += x3;
      M.g(1, 4) += x*y2;
      M.g(1, 5) += x2*y;

      M.g(2, 0) += y;
      M.g(2, 1) += x*y;
      M.g(2, 2) += y2;
      M.g(2, 3) += x2*y;
      M.g(2, 4) += y3;
      M.g(2, 5) += x*y2;

      M.g(3, 0) += x2;
      M.g(3, 1) += x3;
      M.g(3, 2) += x2*y;
      M.g(3, 3) += x4;
      M.g(3, 4) += x2*y2;
      M.g(3, 5) += x3*y;

      M.g(4, 0) += y2;
      M.g(4, 1) += x*y2;
      M.g(4, 2) += y3;
      M.g(4, 3) += x2*y2;
      M.g(4, 4) += y4;
      M.g(4, 5) += x*y3;

      M.g(5, 0) += x*y;
      M.g(5, 1) += x2*y;
      M.g(5, 2) += x*y2;
      M.g(5, 3) += x3*y;
      M.g(5, 4) += x*y3;
      M.g(5, 5) += x2*y2;
    }
                                                                
    //M.print(stdout);
    //R.print(stdout);
    gaussovaEliminace(M, R);
    //R.print(stdout);
    //enter();
    for (int j = 0; j < 5; j++)
      this->approx_coef[strainID][j][i] = R.v[j+1];
    this->Epsilon_Tau[strainID][i] = R.v[0];
      }
    }
  }
  else {
    _errorr("Not implemented yet.");
  }

  /*printf("\nIncl EpsTau:\n");
    for (int i = 0; i < this->n_approx_points; i++) {
    //for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)
    printf("%lf %lf %lf\n", this->approx_points_eps_tau[i][0], this->approx_points_eps_tau[i][1], this->approx_points_eps_tau[i][2]);
    }*/
  //enter();
}

void Inclusion :: Eps02EpsTau (double *e_tau, const double *e_0)
{
  if (P->give_dimension() == 2) {
    ((InclusionRecord2D*)this)->pointInside = 1;
    ((InclusionRecord2D*)this)->x[0] = 0.0;
    ((InclusionRecord2D*)this)->x[1] = 0.0;
    ((InclusionRecord2D*)this)->compute_matrix_D();
  }
  else
    _errorr("3D not implemented yet.");
  
  double m[5], n[5], Cs[5], B[5], e_z[3], e_z_t[3];
  
  SubtractTwoVectors(C, P->matrix->C, Cs, 5);
  
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(n, Cs, S);
  AddVector(P->matrix->C, n, 5);
  giveInverseMatrix3x3to5(m, n);
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(B, m, Cs);
  MultiplyVector(B, 5, -1.0);
  //for (int i = 0; i < 3; i++)e_z_t[i] = P->matrix->remote_strains()[strainID][i] + e_z_add[i];
  rotateStrain_G2L(e_0, e_z);
  MatrixOperations::giveTVproduct_3is3x3to5and3(e_tau, B, e_z);
}


void Inclusion :: EAdeg2rad (void)
{
  for (int i=0; i<P->give_EA_RANGE(); i++)
    eAngles[i] = DEG2RAD( eAngles[i] );
}


//local constant and macros definitions
//*************************************
#define _ACT_RAD_MULT_2d     6. //must be type double
#define _ACT_RAD_MULT_3d     6. //must be type double

//Macros of the function: giveTransformationStrainOperator:
//--------------------------------------------------------
#define Power( a, b )      ( ( b == 2 )? SQR( a ) : pow( a, b ) )//because of MATHEMATICA CForm[] syntax
//Components of stiffness tensor C
/* 1st row */
#define L1111              C[__1111_]
#define L1122              C[__1122_]
#define L1133              C[__1133_]
/* 2nd row */
#define L2211              C[__2211_]
#define L2222              C[__2222_]
#define L2233              C[__2233_]
/* 3rd row */
#define L3311              C[__3311_]
#define L3322              C[__3322_]
#define L3333              C[__3333_]
/* 4th row */
#define L1212              C[__1212_]
/* 5th row */
#define L2323              C[__2323_]
/* 6th row */
#define L1313              C[__1313_]
//Components of stiffness tensor C1
/* 1st row */
#define Ll1111             C1[__1111_]
#define Ll1122             C1[__1122_]
#define Ll1133             C1[__1133_]
/* 2nd row */
#define Ll2211             C1[__2211_]
#define Ll2222             C1[__2222_]
#define Ll2233             C1[__2233_]
/* 3rd row */
#define Ll3311             C1[__3311_]
#define Ll3322             C1[__3322_]
#define Ll3333             C1[__3333_]
/* 4th row */
#define Ll1212             C1[__1212_]
/* 5th row */
#define Ll2323             C1[__2323_]
/* 6th row */
#define Ll1313             C1[__1313_]
//Components of the Eshelby tensor
/* 1st row */
#define S1111              S[__1111_]
#define S1122              S[__1122_]
#define S1133              S[__1133_]
/* 2nd row */
#define S2211              S[__2211_]
#define S2222              S[__2222_]
#define S2233              S[__2233_]
/* 3rd row */
#define S3311              S[__3311_]
#define S3322              S[__3322_]
#define S3333              S[__3333_]
/* 4th row */
#define S1212              S[__1212_]
/* 5th row */
#define S2323              S[__2323_]
/* 6th row */
#define S1313              S[__1313_]


void Inclusion :: giveTransformationStrainOperator (double oper[12], const double C[12], const double C1[12], const double S[12])
{
  //S[__1133_] = S[__2233_] = S[__3311_] = S[__3322_] = S[__2323_] = S[__1313_] = S[__3333_]  = 0.0;
  //C[__1133_] = C[__2233_] = C[__3311_] = C[__3322_] = C[__2323_] = C[__1313_]  = C[__3333_]= 0.0;
  
  double
    B = Power(L1111,3) + 2*Power(L1122,3) - Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*
    (S1133*S2222*S3311 - S1122*S2233*S3311 - S1133*S2211*S3322 + S1111*S2233*S3322 +
     S1122*S2211*S3333 - S1111*S2222*S3333) +
    Power(L1111,2)*(Ll1122*(S1122 + S1133 + S2211 + S2233 + S3311 + S3322) +
            Ll1111*(S1111 + S2222 + S3333)) +
    Power(L1122,2)*(Ll1111*(-S1111 + S1122 + S1133 + S2211 - S2222 + S2233 + S3311 + S3322 - S3333) +
            2*Ll1122*(S1111 + S2222 + S3333)) +
    L1122*(Ll1111 - Ll1122)*(Ll1111*(S1122*S2233 + S1122*S3311 - S2222*S3311 + S2233*S3311 +
                     S2211*S3322 + S1133*(S2211 - S2222 + S3322) -
                     S1111*(S2233 + S3322) - (S1122 + S2211)*S3333) +
                 2*Ll1122*(S1122*S2211 + S1133*S3311 + S2233*S3322 - S2222*S3333 -
                       S1111*(S2222 + S3333))) -
    L1111*(3*Power(L1122,2) + L1122*(Ll1111*(S1122 + S1133 + S2211 + S2233 + S3311 + S3322) +
                     Ll1122*(2*S1111 + S1122 + S1133 + S2211 + 2*S2222 + S2233 +
                         S3311 + S3322 + 2*S3333)) +
       (Ll1111 - Ll1122)*(Ll1111*(S1122*S2211 + S1133*S3311 + S2233*S3322 -
                      S2222*S3333 - S1111*(S2222 + S3333)) +
                  Ll1122*(-(S2222*S3311) + S2233*S3311 + S2211*S3322 +
                      S2233*S3322 + S1133*(S2211 - S2222 + S3311 + S3322) +
                      S1122*(S2211 + S2233 + S3311 - S3333) - (S2211 + S2222)*S3333 -
                      S1111*(S2222 + S2233 + S3322 + S3333))));
  //end of euxiliary variable definitions
  
  //----------------------------------------------------------------------------------------------------
  //Components of the @oper tensor:
  /* 1st row */
  oper[__1111_] = (-(Power(L1111,2)*Ll1111) + Power(L1122,2)*(Ll1111 - 2*Ll1122) +
           Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S2233*S3322 - S2222*S3333) +
           L1122*(Ll1111 - Ll1122)*(Ll1111*(S2233 + S3322) + 2*Ll1122*(S2222 + S3333)) +
           L1111*(2*L1122*Ll1122 - (Ll1111 - Ll1122)*(Ll1111*(S2222 + S3333) +
                                  Ll1122*(S2222 + S2233 + S3322 + S3333))))/B;
  
  oper[__1122_] = -((Power(L1122,2)*Ll1111 + Power(L1111,2)*Ll1122 - L1111*
             (L1122*(Ll1111 + Ll1122) + (Ll1111 - Ll1122)*
              (Ll1111*S1122 + Ll1122*(S1122 + S1133 + S3322 - S3333))) +
             L1122*(Ll1111 - Ll1122)*(2*Ll1122*S1122 + Ll1111*(S1133 + S3322 - S3333)) +
             Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S1133*S3322 - S1122*S3333))/B);
  
  oper[__1133_] = (-(Power(L1122,2)*Ll1111) - Power(L1111,2)*Ll1122 +
           Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S1133*S2222 - S1122*S2233) -
           L1122*(Ll1111 - Ll1122)*(2*Ll1122*S1133 + Ll1111*(S1122 - S2222 + S2233)) +
           L1111*(L1122*(Ll1111 + Ll1122) + (Ll1111 - Ll1122)*
              (Ll1111*S1133 + Ll1122*(S1122 + S1133 - S2222 + S2233))))/B;
  
  //----------------------------------------------------------------------------------------------------
  /* 2nd row */
  oper[__2211_] = -((Power(L1122,2)*Ll1111 + Power(L1111,2)*Ll1122 -
             L1111*(L1122*(Ll1111 + Ll1122) + (Ll1111 - Ll1122)*
                (Ll1111*S2211 + Ll1122*(S2211 + S2233 + S3311 - S3333))) +
             L1122*(Ll1111 - Ll1122)*(2*Ll1122*S2211 + Ll1111*(S2233 + S3311 - S3333)) +
             Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S2233*S3311 - S2211*S3333))/B);

  oper[__2222_] = -((Power(L1111,2)*Ll1111 - Power(L1122,2)*(Ll1111 - 2*Ll1122) -
             Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S1133*S3311 - S1111*S3333) -
             L1122*(Ll1111 - Ll1122)*(Ll1111*(S1133 + S3311) + 2*Ll1122*(S1111 + S3333)) +
             L1111*(-2*L1122*Ll1122 + (Ll1111 - Ll1122)*
                (Ll1111*(S1111 + S3333) + Ll1122*(S1111 + S1133 + S3311 + S3333))))/B);

  oper[__2233_] = -((Power(L1122,2)*Ll1111 + Power(L1111,2)*Ll1122 -
             L1122*(Ll1111 - Ll1122)*(Ll1111*(S1111 - S1133 - S2211) - 2*Ll1122*S2233) +
             Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S1133*S2211 - S1111*S2233) -
             L1111*(L1122*(Ll1111 + Ll1122) + (Ll1111 - Ll1122)*
                (Ll1111*S2233 + Ll1122*(-S1111 + S1133 + S2211 + S2233))))/B);
  
  //----------------------------------------------------------------------------------------------------
  /* 3rd row */
  oper[__3311_] = (-(Power(L1122,2)*Ll1111) - Power(L1111,2)*Ll1122 +
           Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S2222*S3311 - S2211*S3322) -
           L1122*(Ll1111 - Ll1122)*(2*Ll1122*S3311 + Ll1111*(S2211 - S2222 + S3322)) +
           L1111*(L1122*(Ll1111 + Ll1122) + (Ll1111 - Ll1122)*
              (Ll1111*S3311 + Ll1122*(S2211 - S2222 + S3311 + S3322))))/B;

  oper[__3322_] = (-(Power(L1122,2)*Ll1111) - Power(L1111,2)*Ll1122 +
           L1122*(Ll1111 - Ll1122)*(Ll1111*(S1111 - S1122 - S3311) - 2*Ll1122*S3322) -
           Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S1122*S3311 - S1111*S3322) +
           L1111*(L1122*(Ll1111 + Ll1122) + (Ll1111 - Ll1122)*
              (Ll1111*S3322 + Ll1122*(-S1111 + S1122 + S3311 + S3322))))/B;

  oper[__3333_] = (-(Power(L1111,2)*Ll1111) + Power(L1122,2)*(Ll1111 - 2*Ll1122) +
           Power(Ll1111 - Ll1122,2)*(Ll1111 + 2*Ll1122)*(S1122*S2211 - S1111*S2222) +
           L1122*(Ll1111 - Ll1122)*(Ll1111*(S1122 + S2211) + 2*Ll1122*(S1111 + S2222)) +
           L1111*(2*L1122*Ll1122 - (Ll1111 - Ll1122)*
              (Ll1111*(S1111 + S2222) + Ll1122*(S1111 + S1122 + S2211 + S2222))))/B;
  
  //----------------------------------------------------------------------------------------------------
  /* 4th row */
  oper[__1212_] = -(Ll1212/(L1212 + Ll1212*S1212));
  
  //----------------------------------------------------------------------------------------------------
  /* 5th row */
  oper[__2323_] = -(Ll1212/(L1212 + Ll1212*S2323));
  
  //----------------------------------------------------------------------------------------------------
  /* 6th row */
  oper[__1313_] = -(Ll1212/(L1212 + Ll1212*S1313));
  
  //----------------------------------------------------------------------------------------------------
  
#ifdef TOM_DEBUG
  //oper[__1133_] = oper[__2233_] = oper[__3311_] = oper[__3322_] = oper[__2323_] = oper[__1313_] = oper[__3333_] = 0.0;
  
  std::cout << "Matice C" << std::endl;
  for (int i = 0; i < 12; i++)
    std::cout << "C[" << i << "] = " << C[i] << std::endl;
  std::cout << std::endl;
  
  std::cout << "Matice C* = C1 - C" << std::endl;
  for (int i = 0; i < 12; i++)
    std::cout << "C*[" << i << "] = " << C1[i] << std::endl;
  std::cout << std::endl;

  std::cout << "Matice S" << std::endl;
  for (int i = 0; i < 12; i++)
    std::cout << "S[" << i << "] = " << S[i] << std::endl;
  std::cout << std::endl;

  std::cout << "Matice B" << std::endl;
  for (int i = 0; i < 12; i++)
    std::cout << "B[" << i << "] = " << oper[i] << std::endl;
  std::cout << std::endl;
#endif
  
  return ;
} // end of function: giveTransformationStrainOperator( )

void Inclusion :: give_EshelbyMatrixFull (double** answer) const
{
#ifdef DEBUG
  if (answer == NULL)  _errorr( "Field for Eshelby tensor must be allocated" );
#endif

  if (P->twodim) MatrixOperations::copy2DeshelbyTensor_reduced2full (this->S, answer);
  else           MatrixOperations::copy3DeshelbyTensor_reduced2full (this->S, answer);
}

void Inclusion :: give_EshelbyMatrixReduced (double* answer) const
{
#ifdef DEBUG
  if (answer == NULL)  _errorr( "Field for Eshelby tensor must be allocated" );
#endif

  CopyVector(this->S, answer, P->give_ISO_C_RANGE());
}

void Inclusion :: give_StiffnessMatrixFull (double* answer) const
{
#ifdef DEBUG
  if (answer == NULL)  _errorr( "Field for Eshelby tensor must be allocated" );
#endif

  if (P->twodim) MatrixOperations::copy2DeshelbyTensor_reduced2full( this->C, answer);
  else           MatrixOperations::copy3DeshelbyTensor_reduced2full( this->C, answer);
}

void Inclusion :: give_StiffnessMatrixReduced (double* answer) const
{
#ifdef DEBUG
  if (answer == NULL)  _errorr( "Field for Eshelby tensor must be allocated" );
#endif

  CopyVector (this->C, answer, P->give_ISO_C_RANGE());
}

void Inclusion :: give_TeMatrix_G2L (double* answer) const
{
#ifdef DEBUG
  if (answer == NULL)  _errorr( "Field for Eshelby tensor must be allocated" );
#endif

  CopyVector(this->Te, answer, P->give_VM_TENS_RANGE() * P->give_VM_TENS_RANGE());
}
void Inclusion :: give_TeMatrix_L2G (double* answer) const
{
#ifdef DEBUG
  if (answer == NULL)  _errorr( "Field for Eshelby tensor must be allocated" );
#endif

  CopyVector(this->TeInv, answer, P->give_VM_TENS_RANGE() * P->give_VM_TENS_RANGE());
}







const double** Inclusion :: give_EshelbyPertDisplc_internal (int lc, int nlc, const double *coords)
{
  if (P->twodim)  ((InclusionRecord2D*)this)->getDisplacement(coords, globPert_displc, lc, nlc, 1);
  else            ((InclusionRecord3D*)this)->esuf->giveEshelbyDisplacementOfOnePoint (globPert_displc, coords, lc, nlc);
  
  return (const double **)globPert_displc;
}
//void Inclusion :: give_EshelbyPertDisplc_internal (double **displc, int lc, int nlc, const double *coords)
//{
//  CopyArray2D(this->give_EshelbyPertDisplc_internal(lc, nlc, coords), displc, nlc, P->give_VECT_RANGE());
//}



const double **Inclusion :: give_EshelbyPertStrain_internal (int lc, int nlc)
{
  if (globPert_strain == NULL) {
    globPert_strain = new double*[P->give_nLC()];
    memset (globPert_strain, 0, P->give_nLC() * sizeof(double*));
  }
  
  for (int i = lc; i<lc + nlc; i++)
    if (globPert_strain[i] == NULL) {
      globPert_strain[i] = new double[P->give_VM_TENS_RANGE()];
      
      if (P->twodim) MatrixOperations::giveTVproduct_3is3x3to5and3(globPert_strain[i], this->S, ((InclusionRecord2D*)this)->Epsilon_Tau[i]);
      else           MatrixOperations::giveTVproduct_6is6x6to12and6(globPert_strain[i], this->S, ((InclusionRecord3D*)this)->locEigStrain_LC[i]);
      
      this->rotateStrain_L2G(globPert_strain[i], globPert_strain[i]);
    }
  
  return (const double **)globPert_strain + lc;
}
//void Inclusion :: give_EshelbyPertStrain_internal (double **strain, int lc, int nlc)
//{
//  CopyArray2D(this->give_EshelbyPertStrain_internal(lc, nlc), strain, nlc, P->give_VM_TENS_RANGE());
//}



const double** Inclusion :: give_EshelbyPertStress_internal (int lc, int nlc)
{
  if (globPert_stress == NULL) {
    globPert_stress = new double*[P->give_nLC()];
    memset (globPert_stress, 0, P->give_nLC() * sizeof(double*));
  }
  
  // termitovo tady udela vypocet pres this->C, je to rychlejsi  mozna
  for (int i=lc; i<lc+nlc; i++)
    if (globPert_stress[i] == NULL) {
      globPert_stress[i] = new double[P->give_VM_TENS_RANGE()];
      
      if (P->twodim) {
    vektor eps(3);
    
    eps.vynulovat();
    //eps.odecist(((InclusionRecord2D*)this)->Epsilon_Tau[i]);
    for (int j = 0; j < P->give_VM_TENS_RANGE(); j++)eps.v[j] -= ((InclusionRecord2D*)this)->Epsilon_Tau[i][j];
    rotate_vector(eps, -((InclusionRecord2D*)this)->eAngles[0]);

        AddVector (*this->give_EshelbyPertStrain_internal(i, 1), eps.v, 3);
    
    MatrixOperations::giveTVproduct_3is3x3to5and3(this->globPert_stress[i], P->matrix->C, eps.v);
      }
      else {
    double effStrain[6];
    SubtractTwoVectors (*this->give_EshelbyPertStrain_internal(i, 1), ((InclusionRecord3D*)this)->globEigStrain_LC[i], effStrain, 6);
    
    MatrixOperations::giveTVproduct_6is6x6to12and6 (globPert_stress[i], P->matrix->C, effStrain);
      }
    }
  
  return (const double **)globPert_stress+lc;
}
// ///
// void Inclusion :: give_EshelbyPertStress_internal (double **stress, int lc, int nlc)
// {
//   CopyArray2D(this->give_EshelbyPertStress_internal(lc, nlc), stress, nlc, P->give_VM_TENS_RANGE());
// }



//  ///
//  void Inclusion::add_EshelbyPertDisplc_internal_SIFCM(double **displc, int lc, int nlc, const double *coords, double **e_t_SIFCM) {
//    errol;
//  }


void Inclusion::add_EshelbyPertStrain_internal_SIFCM(double **strain, int lc, int nlc, double **e_t_SIFCM)
{
  double *globPert_strain_SIFCM = new double[P->give_VM_TENS_RANGE()];
  
  for (int i = 0; i < nlc; i++) {
    if (P->twodim) MatrixOperations::giveTVproduct_3is3x3to5and3 (globPert_strain_SIFCM, this->S, e_t_SIFCM[i]);
    else           MatrixOperations::giveTVproduct_6is6x6to12and6(globPert_strain_SIFCM, this->S, e_t_SIFCM[i]);
    
    this->rotateStrain_L2G (globPert_strain_SIFCM, globPert_strain_SIFCM);
    AddVector(globPert_strain_SIFCM, strain[i], P->give_VM_TENS_RANGE());
  }
  
  delete [] globPert_strain_SIFCM;
}
void Inclusion::add_EshelbyPertStress_internal_SIFCM(double **stress, int lc, int nlc, double **e_t_SIFCM, double **e_s_ext_add)
{
  double *globPert_stress_SIFCM = new double[P->give_VM_TENS_RANGE()];
  
  for (int i = 0; i < nlc; i++) {
    
    if (P->twodim) {
      vektor eps(3);
      for (int j = 0; j < 3; j++) eps.v[j] = -e_t_SIFCM[i][j];

      rotate_vector (eps, -((InclusionRecord2D*)this)->eAngles[0]);
      
      for (int j = 0; j < 3; j++) eps.v[j] += e_s_ext_add[i][j];
      
      MatrixOperations::giveTVproduct_3is3x3to5and3 (globPert_stress_SIFCM, P->matrix->C, eps.v);

      AddVector(globPert_stress_SIFCM, stress[i], P->give_VM_TENS_RANGE());
      
    }
    else {
      errol;  // tady to neni spravne, ne?
      //double effStrain[6];
      //SubtractTwoVectors(*this->give_EshelbyPertStrain_internal(i, 1), ((InclusionRecord3D*)this)->globEigStrain_LC[i], effStrain, 6);
      
      //MatrixOperations::giveTVproduct_6is6x6to12and6(globPert_stress_SIFCM, P->matrix->C, effStrain);
    }
    
  }
}



void Inclusion :: give_EshelbyPertFields_external (Point *point, int lc, int nlc, bool disp, bool strn)
{
  if (P->twodim) {
    if (strn)  ((InclusionRecord2D*)this)->giveExtStrainPert (point->x, point->strain, lc, nlc);
    if (disp)  ((InclusionRecord2D*)this)->getDisplacement   (point->x, point->displacement, lc, nlc, 0);
  }
  else
    ((InclusionRecord3D*)this)->esuf->giveEshelbyFieldsOfOnePoint(point, lc, nlc, disp, strn);

}

double Inclusion :: compute_supplement_energy (int lc)
{
  // ted si tady pocitam nejaky picovinky, Mura ma na tuto energii vzorecek
  
  //  Suplement of the isotropic elastic stiffness tensor of an inclusion stored in reduced form.
  double *Cs = new double[P->give_ISO_C_RANGE()];
  SubtractTwoVectors (this->C, P->matrix->C, Cs, P->give_ISO_C_RANGE());

  double *stress = new double[P->nstrain_one()];

  if (P->twodim) MatrixOperations::giveTVproduct_3is3x3to5and3  (stress, this->C, P->matrix->remote_strains()[lc]);
  else           MatrixOperations::giveTVproduct_6is6x6to12and6 (stress, this->C, P->matrix->remote_strains()[lc]);

  
  double E = give_VectorVectorProduct (P->matrix->remote_strains()[lc], stress, P->nstrain_one());
  E *= this->give_volume();

  delete [] Cs;

  return E;
}


//  //********************************************************************************************************
//  // description: function prints the record of an inclusion
//  // last edit:   02. 03. 2010
//  // -------------------------------------------------------------------------------------------------------
//  // note:        @inclRec - inclusion record of one inclusion
//  //********************************************************************************************************
//  void inclRecordInit::printInclRecord( InclusionRecord inclRec, const char * notice )
//  {
//    matrixOperations mtrxOper;
//    eshelbySoluUniformField eshSolu;
//  
//    /*
//   if verbose
//    double auxMtrx[36] = { 0., 0., 0., 0., 0., 0.,
//               0., 0., 0., 0., 0., 0.,
//               0., 0., 0., 0., 0., 0.,
//               0., 0., 0., 0., 0., 0.,
//               0., 0., 0., 0., 0., 0.,
//               0., 0., 0., 0., 0., 0.,
//    };
//    double auxVector[6] = { 0., 0., 0., 0., 0., 0. };
//  
//    mtrxOper.giveMatrixMatrixProduct( inclRec.TeInv, inclRec.Te, auxMtrx, 6 );
//    mtrxOper.giveMatrixVectorProduct( inclRec.TeInv, inclRec.locRemoteStrain, auxVector, 6 );
//  #endif
//    */
//  
//    printf( "%s", notice );
//    printf( "Poisson's ratio: %e\n", inclRec.nu );
//    printf( "Young's modulus: %e\n", inclRec.E );
//    printf( "Shape:           %d\n", ( int ) inclRec.shape );
//    mtrxOper.printVectorRowForm( inclRec.origin, 3, "Origin:\n" );
//    mtrxOper.printVectorRowForm( inclRec.a, 3, "Semiaxes dimensions:\n" );
//    mtrxOper.printVectorRowForm( inclRec.sort_a, 3, "Sorted semiaxes dimensions:\n" );
//    mtrxOper.printVectorRowForm( inclRec.eAngles, 3, "Euller angles:\n" );
//    mtrxOper.printVectorRowForm( inclRec.globRemoteStrain[0], 6, "Global remote strain field:\n" );
//    mtrxOper.printVectorRowForm( inclRec.locRemoteStrain[0], 6, "Local remote strain field:\n" );
//    eshSolu.printIsoTensor( inclRec.SInv, "Inverse Eshelby tensor:\n", _STANDARD_SYNTAX_ );
//    mtrxOper.printMatrix( inclRec.T, 3, "Global->local vector transformation matrix:\n" );
//    mtrxOper.printMatrix( inclRec.TInv, 3, "Local->global vector transformation matrix:\n" );
//    mtrxOper.printMatrix( inclRec.Te, 6, "Global->local strain tensor transformation matrix:\n" );
//    mtrxOper.printMatrix( inclRec.TeInv, 6, "Local->global strain tensor transformation matrix:\n" );
//    /*
//      if  verbose
//      mtrxOper.printMatrix( auxMtrx, 6, "debug: Te * TeInv (suppose to be unit matrix):\n" );
//      mtrxOper.printVectorRowForm( auxVector, 6, "debug: Global remote strain field:\n" );
//      #endif
//    */
//    printf( "Action radius: %e\n", inclRec.actionRadius );
//    printf( "Number of acting inclusions: %d\n", inclRec.noActingIncls );
//    mtrxOper.printVectorRowForm( inclRec.actingIncls, inclRec.noActingIncls,
//                     "Numbers of acting inclusions:\n" );
//    eshSolu.printIsoTensor( inclRec.C, "Total stifness matrix:\n", _STANDARD_SYNTAX_ );
//    mtrxOper.printVectorRowForm( inclRec.locEigStrain[0], 6, "Initial local equivalent eigenstrain:\n" );
//    mtrxOper.printVectorRowForm( inclRec.globEigStrain[0], 6, "Initial global equivalent eigenstrain:\n" );
//    mtrxOper.printVectorRowForm( inclRec.globEigStrainTotal[0], 6, "Initial total global equivalent eigenstrain:\n" );
//  
//    return ;
//  }//end of function: printInclRecord( )









// ***********************************************************************************************************
// ***                                         3D INCLUSION RECORD                                         ***
// ***********************************************************************************************************

InclusionRecord3D :: InclusionRecord3D (long i, const Problem *p) : Inclusion(i, p)
{
  locEigStrain_LC = globEigStrain_LC = NULL;
  ellInt = NULL;
  esuf = NULL;
}

InclusionRecord3D :: ~InclusionRecord3D ()
{
  DeleteArray2D ( locEigStrain_LC, P->give_nLC());
  DeleteArray2D (globEigStrain_LC, P->give_nLC());

  delete ellInt;
  delete esuf;
}

bool InclusionRecord3D :: scan_locEigStrain_LC (Stream *stream, int lc)
{
  if (locEigStrain_LC == NULL)  this->allocate_nlc_fields ();
  return ST_scan_array (stream, P->give_VM_TENS_RANGE(), this->locEigStrain_LC[lc]);
}
bool InclusionRecord3D :: scan_globEigStrain_LC (Stream *stream, int lc)
{
  if (globEigStrain_LC == NULL)  this->allocate_nlc_fields ();
  return ST_scan_array (stream, P->give_VM_TENS_RANGE(), this->globEigStrain_LC[lc]);
}

void InclusionRecord3D :: computeVolume (void)
{
  switch(shape) {
  case IS_ELLIPSOID :
  case IS_SPHERE :
  //case IS_PENNY :
  //case IS_FLAT_ELLIPSOID :
  case IS_OBLATE_SPHEROID :
  case IS_PROLATE_SPHEROID :
    this->volume = 4.0/3.0 * PI * a[0] * a[1] * a[2];
    break;
  //case IS_ELLIPTIC_CYLINDER :
  //case IS_CYLINDER :
  //  _errorr( "Not yet defined" );
  //  break;
  default:
    _errorr( "Undefined shape" );
    break;
  }
}


void InclusionRecord3D :: allocate_nlc_fields (void)
{
  Inclusion :: allocate_nlc_fields();
  
  int nlc = P->give_nLC();
  if (nlc==0) _errorr("zero nlc");
  
  if (locEigStrain_LC == NULL) {
    // termitovo predelat na fci AllocateArray..
    locEigStrain_LC  = new double*[nlc];
    for (int i=0; i<nlc; i++)
      locEigStrain_LC [i] = new double[P->give_VM_TENS_RANGE()];
  }
  
  if (globEigStrain_LC == NULL) {
    globEigStrain_LC = new double*[nlc];
    for (int i=0; i<nlc; i++)
      globEigStrain_LC [i] = new double[P->give_VM_TENS_RANGE()];
  }
}

InclusionGeometry InclusionRecord3D :: giveInclusionGeometry (InclusionGeometry shp_o) const
{
  InclusionGeometry shp_n;
  
  double min = P->give_semiaxes_min_difference();
  double max = P->give_semiaxes_max_difference();
  
  bool chmin = P->give_semiaxis_min_difference_change();
  //termitovo ne co chmax??    bool chmax = P->give_semiaxis_max_difference_change();
  
  //semiaxes definition and initialization: a1 > a2 > a3
  double a1 = a[0], a2 = a[1], a3 = a[2];
  
  if (a1 < a2  ||  a2 < a3)  _errorr("semiaxes not sorted");
  if (a1 < 0.0)              _errorr("negative semiaxis");
  
  if (shp_o) {
    switch (shp_o) {
    case IS_ELLIPSOID:         if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
                               if (!( a2 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 2. axis is infinite",               this->id+1, IST_e2s (shp_o));
                   if (!( a3 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 3. axis is infinite",               this->id+1, IST_e2s (shp_o));
                   if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                   if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                   if (!( a3 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 3. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                               if (!( a1 > a2))     _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not bigger then 2. one", this->id+1, IST_e2s (shp_o));
                               if (!( a2 > a3))     _errorr3("inclusion %ld with prescribed %s shape: 2. axis is not bigger then 3. one", this->id+1, IST_e2s (shp_o));
      break;
    case IS_SPHERE:            if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
                   if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                   if (!( a1 == a2))    _errorr3("inclusion %ld with prescribed %s shape: 1. and 2. axes are not equal",      this->id+1, IST_e2s (shp_o));
                               if (!( a2 == a3))    _errorr3("inclusion %ld with prescribed %s shape: 2. and 3. axes are not equal",      this->id+1, IST_e2s (shp_o));
      break;
    case IS_OBLATE_SPHEROID:   if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
                   if (!( a3 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 3. axis is infinite",               this->id+1, IST_e2s (shp_o));
                   if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                   if (!( a3 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 3. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                   if (!( a1 == a2))    _errorr3("inclusion %ld with prescribed %s shape: 1. and 2. axes are not equal",      this->id+1, IST_e2s (shp_o));
                               if (!( a1 > a3))     _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not bigger then 3. one", this->id+1, IST_e2s (shp_o));
      break;
    case IS_PROLATE_SPHEROID:  if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
                               if (!( a2 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 2. axis is infinite",               this->id+1, IST_e2s (shp_o));
                   if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                   if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
                               if (!( a1 > a2))     _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not bigger then 2. one", this->id+1, IST_e2s (shp_o));
                   if (!( a2 == a3))    _errorr3("inclusion %ld with prescribed %s shape: 2. and 3. axes are not equal",      this->id+1, IST_e2s (shp_o));
     break;
    //case IS_ELLIPTIC_CYLINDER: if (!( a1 == INFTY)) _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not infinite",           this->id+1, IST_e2s (shp_o));
    //                 if (!( a2 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 2. axis is infinite",               this->id+1, IST_e2s (shp_o));
    //                 if (!( a3 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 3. axis is infinite",               this->id+1, IST_e2s (shp_o));
    //                 if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
    //                 if (!( a3 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 3. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
    //                           if (!( a2 > a3))     _errorr3("inclusion %ld with prescribed %s shape: 2. axis is not bigger then 3. one", this->id+1, IST_e2s (shp_o));
    //  break;
    //case IS_CYLINDER:          if (!( a1 == INFTY)) _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not infinite",           this->id+1, IST_e2s (shp_o));
    //                           if (!( a2 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 2. axis is infinite",               this->id+1, IST_e2s (shp_o));
    //                 if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
    //                           if (!( a2 == a3))    _errorr3("inclusion %ld with prescribed %s shape: 2. and 3. axes are not equal",      this->id+1, IST_e2s (shp_o));
    //  break;
    //case IS_PENNY:             if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
    //                           if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
    //                 if (!( a1 == a2))    _errorr3("inclusion %ld with prescribed %s shape: 1. and 2. axes are not equal",      this->id+1, IST_e2s (shp_o));
    //                 if (!( a3 == 0.0))   _errorr3("inclusion %ld with prescribed %s shape: 3. axis is not 0.0",                this->id+1, IST_e2s (shp_o));
    //  break;
    //case IS_FLAT_ELLIPSOID:    if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
    //                 if (!( a2 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 2. axis is infinite",               this->id+1, IST_e2s (shp_o));
    //                 if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
    //                 if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
    //                           if (!( a1 > a2))     _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not bigger then 2. one", this->id+1, IST_e2s (shp_o));
    //                 if (!( a3 == 0.0))   _errorr3("inclusion %ld with prescribed %s shape: 3. axis is not 0.0",                this->id+1, IST_e2s (shp_o));
    //  break;
    default: _errorr2("Not supported inclusion shape %s", IST_e2s (shp_o));
    }
  }
  
  // detection of the shape with stable analysis
  // a1 = 0
  if (a1 == 0.0)  _errorr2("inclusion %ld: zero 1. axis - unsupported", this->id+1);
  // a1 = infty
  else if (a1 == INFTY) {
    // a2 = 0
    if (a2 == 0.0)    _errorr2("inclusion %ld: infinity 1. axis and zero 2. axis - unsupported ", this->id+1);
    // a2 = a1
    else if (a2 == INFTY)  _errorr2("inclusion %ld: infinity 1. axis and infinity 2. axis - unsupported", this->id+1);
    // a2 = any
    else {
      // a3 = 0
      if (a3 < max*a2)  _errorr2("inclusion %ld: infinity 1. axis and zero 3. axis - unsupported ", this->id+1);
      // a3 = a2
      else if (a2 - a3 < min*a2)  errol; //shp_n = IS_CYLINDER;
      // a3 = any
      else                        errol; //shp_n = IS_ELLIPTIC_CYLINDER;
    }
  }
  // a1 = any
  else {
    // a2 = 0.0
    if (a2 < max*a1)  _errorr2("inclusion %ld: any 1. axis and zero 2. axis - unsupported ", this->id+1);
    // a2 = a1
    else if (a1 - a2 < min*a1) {
      // a3 = 0
      if (a3 < max*a2)            errol; //shp_n = IS_PENNY;
      // a3 = a2
      else if (a2 - a3 < min*a2)  shp_n = IS_SPHERE;
      // a3 = any
      else                        shp_n = IS_OBLATE_SPHEROID;
    }
    // a2 = any
    else {
      // a3 = 0
      if (a3 < max*a2)            shp_n = IS_ELLIPSOID;
      // a3 = a2
      else if (a2 - a3 < min*a2)  shp_n = IS_PROLATE_SPHEROID;
      // a3 = any
      else                        shp_n = IS_ELLIPSOID;
    }
  }
  
  // Toto se musi promyslet, ktere konkretni kombinace tvaru mohou tvarove degenerovat.
  //
  if (shp_o) {
    if (shp_o != shp_n) {
      switch (shp_n) {
      case IS_SPHERE:
    if (chmin) {
      if (shp_o == IS_ELLIPSOID) {
        _warningg3("inclusion %ld with prescribed %s shape: shape is different - it was changed", this->id+1, IST_e2s (shp_o));
        return shp_n;
      }
      else _errorr3("inclusion %ld with prescribed %s shape: unsupported shape conversion", this->id+1, IST_e2s (shp_o));
    }
    else { _warningg("shape is different"); return shp_o; }
        break;
      default: _errorr2("Not supported inclusion shape %s", IST_e2s (shp_o));
      }
    }
    else
      return shp_o;
  }
  else
    return shp_n;
  
  _errorr("return have to be used above");
  return IS_VOID;
}

void InclusionRecord3D :: input_data_initialize_and_check_consistency (void)
{
  Inclusion :: input_data_initialize_and_check_consistency ();
  
  if (P->is_converted_to_equivalent() == false) {
    ;
  }
  else {
    ;
  }
  
  if (this->ellInt) errol;
  switch (this->shape) {
  case IS_ELLIPSOID :         ellInt = new eshelbySoluEllipticIntegralsEllipsoid(this);        break;
  case IS_SPHERE :            ellInt = new eshelbySoluEllipticIntegralsSphere(this);           break;
  case IS_OBLATE_SPHEROID :   ellInt = new eshelbySoluEllipticIntegralsOblateSpheroid(this);   break;
  case IS_PROLATE_SPHEROID :  ellInt = new eshelbySoluEllipticIntegralsProlateSpheroid(this);  break;
  //case IS_ELLIPTIC_CYLINDER : ellInt = new eshelbySoluEllipticIntegralsEllipticCylinder(this); break;
  //case IS_CYLINDER :          ellInt = new eshelbySoluEllipticIntegralsCylinder(this);         break;
  //case IS_PENNY :             ellInt = new eshelbySoluEllipticIntegralsPenny(this);            break;
  //case IS_FLAT_ELLIPSOID :    ellInt = new eshelbySoluEllipticIntegralsFlatEllipsoid(this);    break;
  default:
    _errorr( "Unknown shape of inclusion " );
  }
  
  if (this->esuf) errol;
  switch (this->shape) {
  case IS_ELLIPSOID :         esuf = new eshelbySoluUniformField(this);        break;
  case IS_SPHERE :            esuf = new eshelbySoluUniformFieldSphere(this);           break;
  case IS_OBLATE_SPHEROID :   esuf = new eshelbySoluUniformFieldOblateSpheroid(this);   break;
  case IS_PROLATE_SPHEROID :  esuf = new eshelbySoluUniformFieldProlateSpheroid(this);  break;
  //case IS_ELLIPTIC_CYLINDER : esuf = new eshelbySoluUniformFieldEllipticCylinder(this); break;
  //case IS_CYLINDER :          esuf = new eshelbySoluUniformFieldCylinder(this);         break;
  //case IS_PENNY :             esuf = new eshelbySoluUniformFieldPenny(this);            break;
  //case IS_FLAT_ELLIPSOID :    esuf = new eshelbySoluUniformFieldFlatEllipsoid(this);    break;
  default:
    _errorr( "Unknown shape of inclusion " );
  }
}

void InclusionRecord3D :: initialize (const Inclusion *prevInclRec)
{
  const InclusionRecord3D *prevInclRec3D;
  
  if (prevInclRec) {
    prevInclRec3D = dynamic_cast<const InclusionRecord3D*>(prevInclRec);
    if (!prevInclRec3D) _errorr("");
  }
  else
    prevInclRec3D = NULL;
  
  // elliptical Ferers-Dyson potentials of internal points
  this->ellInt->giveEllipticIntegrals (this->eInt, 0., true);
  
  // Eshelby tensor
  //eshSolu.giveEshelbyTensor( this->inclRec[i].S, this->inclRec[i].sort_a, this->inclRec[i].nu,
  //this->inclRec[i].eInt, this->inclRec[i].shape ); // ORIGINAL CODE. Bug?
  this->esuf->giveEshelbyTensor (this->S, this->eInt);
  
  // inverse of the Eshelby tensor
  this->esuf->giveEshelbyTensorInverse (this->SInv, this->S);
  // global->local transformation matrices
  TransformTensors::give_T (this->T,  this->eAngles, _ZXZ_, this->shape);
  TransformTensors::give_Te(this->Te, this->T,              this->shape);
  
  //local->global transformation matrices
  TransformTensors::give_TInv (this->TInv,  this->T,  this->shape);
  TransformTensors::give_TeInv(this->TeInv, this->Te, this->shape);
  
  // initialize action radii
  this->actionRadius = _ACT_RAD_MULT_3d * this->a[0];
  this->SQRactionRadius = SQR( this->actionRadius );
  
  // total and complementary stifness of inclusion
  if (prevInclRec3D && prevInclRec3D->E == this->E && prevInclRec3D->nu == this->nu)
    // convenient if identical stiffnesses in few sets
    CopyVector (prevInclRec3D->C,  this->C,  12);
  else
    Stiffness::giveReducedIsoStiffMatrix (this->C, this->E, this->nu, false);

#ifdef TOM_DEBUG
  std::cout << "Matice B" << std::endl;
  std::cout << inclRec[0].transfTensOper[4] << std::endl;
  std::cout << inclRec[0].transfTensOper[5] << std::endl;
  std::cout << inclRec[0].transfTensOper[7] << std::endl;
  std::cout << inclRec[0].transfTensOper[8] << std::endl;
  std::cout << inclRec[0].transfTensOper[10] << std::endl;
  std::cout << std::endl;
  
  std::cout << "LocRemoteStrain" << std::endl;
  for (int i = 0; i < 6; i++)
    std::cout << "[" << i << "]" << inclRec[0].locRemoteStrain[0][i] << std::endl;
  std::cout << std::endl;
  
  std::cout << "LocEigStrain" << std::endl;
  for (int i = 0; i < 6; i++)
    std::cout << "[" << i << "]" << inclRec[0].locEigStrain[0][i] << std::endl;
  std::cout << std::endl;
  std::cout << std::endl;
#endif

}


void InclusionRecord3D :: SBA_computeInitialStrains (int lc)
{
  double locRemoteStrain[6];    

  // Local remote strain initialization.
  this->rotateStrain_G2L(P->matrix->Remote_strain[lc], locRemoteStrain);
  

  double Cs[12];                 
  SubtractTwoVectors (this->C, P->matrix->C, Cs, 12);

  double B[12];   


  this->giveTransformationStrainOperator (B, P->matrix->C, Cs, this->S);



  // Local transformation eigenstrain this->locEigStrain initialization.
  // A transformation eigenstrain is given by the remote strain applied to an infinite matrix surrounding each inclusion.
  // Matrix-vector multiplication B(6x6) * eps(6x1). Matrix B is in the same form as material stiffness matrix.
  MatrixOperations::giveTVproduct_6is6x6to12and6 (this->locEigStrain, B, locRemoteStrain);
  
  // Global transformation eigenstrain initialization.
  this->rotateStrain_L2G(this->locEigStrain, this->globEigStrainTotal);
  
  return ;
}//end of function: updateStrainsInInclRecord( )


void InclusionRecord3D :: SBA_updateGlobAndLocStrains (Point *point)
{
  double auxStrain[6] = { 0., 0., 0., 0., 0., 0. };

  //local strain preturbation
  this->rotateStrain_G2L(point->strain[0], auxStrain);
  
  //calculation of local transformation strain
  MatrixOperations::giveTVproduct_6is6x6to12and6(auxStrain, SInv, auxStrain);
  
  AddVector(auxStrain, locEigStrainTot, 6);
  
  return;
}//end of function: updateGlobAndLocStrainsInInclRecord( )

void InclusionRecord3D :: addtot ()
{
  CopyVector(locEigStrainTot, locEigStrain, 6);

  //calculation of global transformation strain
  this->rotateStrain_L2G(locEigStrainTot, locEigStrainTot);

  //final update of global total transformation strain
  SumTwoVectors (globEigStrainTotal, locEigStrainTot, globEigStrainTotal, 6);
}


bool InclusionRecord3D :: point_is_affected (const double *point) const
{
  double rx = this->origin[0] - point[0];
  double ry = this->origin[1] - point[1];
  double rz = this->origin[2] - point[2];
  
  if ( ABS( rx ) < this->actionRadius &&
       ABS( ry ) < this->actionRadius &&
       ABS( rz ) < this->actionRadius ) {
    double r = VEC_SQR_NORM_3D( rx, ry, rz );
    if (r <= this->SQRactionRadius)
      return true;
  }
  
  return false;
}


void InclusionRecord3D :: init_EigStrain_LC (int lc)
{
  CopyVector (locEigStrain,       locEigStrain_LC [lc], P->give_VM_TENS_RANGE());
  CopyVector (globEigStrainTotal, globEigStrain_LC[lc], P->give_VM_TENS_RANGE());
}


// ***********************************************************************************************************
// ***                2D ONLY - START          2D INCLUSION RECORD                                         ***
// ***********************************************************************************************************

InclusionRecord2D :: InclusionRecord2D (long i, const Problem *p) : Inclusion(i, p)
{
  // for standa: dusledne tady vse vynulovat, a porovnat jako v .h
  x[0]=0.0;
  x[1]=0.0;
  x_last_derivation[0] = 0.0;
  x_last_derivation[1] = 0.0;
  x_last_IS[0] = 0.0;
  x_last_IS[1] = 0.0;
  lambda=0.0;
  h=0.0;
  H=0.0;
  D.define(3,3);
  L.define(2, 3);
  Sext.define(3,3);
  e_t.define(3);
  
  delta_e_s.define(3);
  delta_e_t.define(3);
  any_derivative_preparation_computed = false;
  any_IS_computed = false;
}

InclusionRecord2D :: ~InclusionRecord2D ()
{
  
}

bool InclusionRecord2D :: scan_locEigStrain_LC (Stream *stream, int lc)
{
  if (Epsilon_Tau == NULL)  this->allocate_nlc_fields ();
  return ST_scan_array (stream, P->give_VM_TENS_RANGE(), this->Epsilon_Tau[lc]);
}


void InclusionRecord2D :: allocate_nlc_fields (void)
{
  Inclusion :: allocate_nlc_fields();
}


InclusionGeometry InclusionRecord2D :: giveInclusionGeometry (InclusionGeometry shp_o) const
{
  InclusionGeometry shp_n;
  
  double min = P->give_semiaxes_min_difference();
  double max = P->give_semiaxes_max_difference();
  
  bool chmin = P->give_semiaxis_min_difference_change();
  // termitovo: na co toto? bool chmax = P->give_semiaxis_max_difference_change();
  
  //semiaxes definition and initialization: a1 > a2 > a3
  double a1 = a[0], a2 = a[1];
  
  if (a1 < a2)   _errorr("semiaxes not sorted");
  if (a1 < 0.0)  _errorr("negative semiaxis");
  
  if (shp_o) {
    switch (shp_o) {
    case IS_ELLIPSE: if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite",               this->id+1, IST_e2s (shp_o));
                     if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0",                    this->id+1, IST_e2s (shp_o));
             if (!( a1 > a2))     _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not bigger then 2. one", this->id+1, IST_e2s (shp_o));
      break;
    case IS_CIRCLE:  if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite", this->id+1, IST_e2s (shp_o));
             if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0",      this->id+1, IST_e2s (shp_o));
             if (!( a1 == a2))    _errorr3("inclusion %ld with prescribed %s shape: 1. and 2. axes are not equal", this->id+1, IST_e2s (shp_o));
      break;
    //case IS_STRIP:   if (!( a1 == INFTY)) _errorr3("inclusion %ld with prescribed %s shape: 1. axis is not infinite", this->id+1, IST_e2s (shp_o));
    //           if (!( a2 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 2. axis is infinite", this->id+1, IST_e2s (shp_o));
    //           if (!( a2 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 2. axis is 0.0", this->id+1, IST_e2s (shp_o));
    //  break;
    //case IS_CRACK:   if (!( a1 < INFTY))  _errorr3("inclusion %ld with prescribed %s shape: 1. axis is infinite", this->id+1, IST_e2s (shp_o));
    //                 if (!( a1 > 0.0))    _errorr3("inclusion %ld with prescribed %s shape: 1. axis is 0.0", this->id+1, IST_e2s (shp_o));
    //           if (!( a2 == 0.0))   _errorr3("inclusion %ld with prescribed %s shape: 3. axis is not 0.0", this->id+1, IST_e2s (shp_o));
    //  break;
    default: _errorr2("Not supported inclusion shape %s", IST_e2s (shp_o));
    }
  }
  
  // detection of the shape with stable analysis
  // a1 = 0
  if (a1 == 0.0)  _errorr2("inclusion %ld: zero 1. axis - unsupported", this->id+1);
  // a1 = infty
  else if (a1 == INFTY) {
    // a2 = 0
    if (a2 == 0.0)         _errorr2("inclusion %ld: infinity 1. axis and zero 2. axis - unsupported ", this->id+1);
    // a2 = a1
    else if (a2 == INFTY)  _errorr2("inclusion %ld: infinity 1. axis and infinity 2. axis - unsupported", this->id+1);
    // a2 = any
    else errol; //shp_n = IS_STRIP;
  }
  // a1 = any
  else {
    // a2 = 0
    if (a2 < max*a1)            errol; //shp_n = IS_CRACK;
    // a2 = a1
    else if (a1 - a2 < min*a1)  shp_n = IS_CIRCLE;
    // a2 = any
    else                        shp_n = IS_ELLIPSE;
  }
  
  // Toto se musi promyslet, ktere konkretni kombinace tvaru mohou tvarove degenerovat.
  //
  if (shp_o) {
    if (shp_o != shp_n) {
      switch (shp_n) {
      case IS_CIRCLE:
    if (chmin) {
      if (shp_o == IS_ELLIPSE) {
        _warningg3("inclusion %ld with prescribed %s shape: shape is different - it was changed", this->id+1, IST_e2s (shp_o));
        return shp_n;
      }
      else _errorr3("inclusion %ld with prescribed %s shape: unsupported shape conversion", this->id+1, IST_e2s (shp_o));
    }
    else { _warningg("shape is different"); return shp_o; }
        break;
      default: _errorr2("Not supported inclusion shape %s", IST_e2s (shp_o));
      }
    }
    else
      return shp_o;
  }
  else
    return shp_n;
  
  _errorr("");
  return IS_VOID;
}

void InclusionRecord2D :: input_data_initialize_and_check_consistency (void)
{
  Inclusion :: input_data_initialize_and_check_consistency ();
  
  if (P->is_converted_to_equivalent() == false) {
    ;
  }
  else {
    ;
  }
  
  this->h = this->derivative_step(100000000);
  this->H = this->derivative_step(10000);
}

void InclusionRecord2D :: initialize (const Inclusion *prevInclRec)
{
  const InclusionRecord2D *prevInclRec2D;
  
  if (prevInclRec) {
    prevInclRec2D = dynamic_cast<const InclusionRecord2D*>(prevInclRec);
    if (!prevInclRec2D) _errorr("");
  }
  else
    prevInclRec2D = NULL;
  
  
  actionRadius = _ACT_RAD_MULT_2d * a[0];
  
  //
  if (prevInclRec2D != NULL && prevInclRec2D->E == this->E && prevInclRec2D->nu == this->nu){
    CopyVector (prevInclRec2D->C,  this->C,  5);
  }
  else
    Stiffness::giveReducedIsoStiffMatrix (this->C, this->E, this->nu, true);

  this->compute_S_int();

  // global->local transformation matrices
  TransformTensors::give_T (this->T,  this->eAngles, _NONE_, this->shape);
  TransformTensors::give_Te(this->Te, this->T,               this->shape);

  //local->global transformation matrices
  TransformTensors::give_TInv (this->TInv,  this->T,    this->shape);
  TransformTensors::give_Te    (this->TeInv, this->TInv, this->shape);
}

void InclusionRecord2D :: SBA_computeInitialStrains (int lc)
{
  this->compute_initial_eps_tau(lc);
        CopyVector(this->e_t.v, this->Epsilon_Tau[lc], P->give_VM_TENS_RANGE());
  this->delta_e_t.rovna_se(this->e_t);
}


bool InclusionRecord2D :: point_is_affected (const double *point) const
{
  return ( VEC_SQR_NORM_2D ( origin[0] - point[0], origin[1] - point[1] )  <=   this->actionRadius * this->actionRadius );
}

double InclusionRecord2D:: S_give(int index)
{
  if (x_last_IS[0] != x[0] || x_last_IS[1] != x[1] || any_IS_computed == false){
    compute_I_S();
    x_last_IS[0] = x[0];
    x_last_IS[1] = x[1];
    any_IS_computed = true;
  }
  if(index==1111)return Sext.g(0,0);
  if(index==2222)return Sext.g(1,1);
  if(index==1122)return Sext.g(0,1);
  if(index==2211)return Sext.g(1,0);
  if(index==1212)return Sext.g(2,2);
  return 0.0;
}

double InclusionRecord2D :: I(int index)
{
  if (x_last_IS[0] != x[0] || x_last_IS[1] != x[1] || any_IS_computed == false){
    compute_I_S();
    x_last_IS[0] = x[0];
    x_last_IS[1] = x[1];
    any_IS_computed = true;
  }
  if(index==1)return I1;
  if(index==2)return I2;
  if(index==11)return I11;
  if(index==22)return I22;
  if(index==12 || index==21)return I12;
  return 0.0;
}

double InclusionRecord2D :: derivative_step(double delitel)
{
  if(delitel==0)delitel=1000000;
  if(a[0]<a[1]){return a[0]/delitel;}else{return a[1]/delitel;}
}

void InclusionRecord2D :: derivative_preparation()
{
  int i;
  P_a[0].x[0]=x[0];     P_a[0].x[1]=x[1]+h; 
  P_a[1].x[0]=x[0]+h;   P_a[1].x[1]=x[1];
  P_b[0].x[0]=x[0]-H;   P_b[0].x[1]=x[1]+H; 
  P_b[1].x[0]=x[0];     P_b[1].x[1]=x[1]+H;
  P_b[2].x[0]=x[0]+H;   P_b[2].x[1]=x[1]+H;
  P_b[3].x[0]=x[0]-H;   P_b[3].x[1]=x[1];
  P_b[4].x[0]=x[0]+H;   P_b[4].x[1]=x[1];
  P_b[5].x[0]=x[0]-H;   P_b[5].x[1]=x[1]-H;
  P_b[6].x[0]=x[0];     P_b[6].x[1]=x[1]-H;
  P_b[7].x[0]=x[0]+H;   P_b[7].x[1]=x[1]-H;
  for(i=0;i<8;i++){
    P_b[i].a[0]=a[0];
    P_b[i].a[1]=a[1];
  }
  for(i=0;i<2;i++){
    P_a[i].a[0]=a[0];
    P_a[i].a[1]=a[1];
  }
  for(i=0;i<8;i++)P_b[i].spocitat();
  for(i=0;i<2;i++)P_a[i].spocitat();
}

//  Pa:      0   
//           X  1
//     
//  Pb:   0  1  2
//        3  X  4
//        5  6  7

double InclusionRecord2D :: dI(int index,int smer1,int smer2)
{
  if (P->diffType != DT_NUMERICAL)_errorr("Only numerical derivation is available.");
  if (pointInside)return 0.0;
  if (x_last_derivation[0] != x[0] || x_last_derivation[1] != x[1] || any_derivative_preparation_computed == false){
    derivative_preparation();
    x_last_derivation[0] = x[0];
    x_last_derivation[1] = x[1];
    any_derivative_preparation_computed = true;
  }
  if(smer2==0){
    if(smer1==1)return (P_a[1].I(index)-I(index))/(h);
    if(smer1==2)return (P_a[0].I(index)-I(index))/(h);
  }
  else{
    if(smer1==smer2){
      if(smer1==1)return ( P_b[4].I(index) - 2.*I(index) + P_b[3].I(index) )/(H*H);
      if(smer1==2)return ( P_b[1].I(index) - 2.*I(index) + P_b[6].I(index) )/(H*H);
    }
    else{
      return (P_b[2].I(index)+P_b[5].I(index)-P_b[0].I(index)-P_b[7].I(index))/(4*H*H);
    }
  }
  return 0; 
}

void InclusionRecord2D :: compute_I_S()
{
  if (a[0] == a[1]){
    if (pointInside){
      lambda = 0.;
      I1 = I2 = 2.*PI;
      I11 = I22 = I12 = PI / pow(a[0], 2);
      Sext.g(1, 1) = Sext.g(0, 0) = 3. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * I11 + (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*I1;
      Sext.g(1, 0) = Sext.g(0, 1) = 1. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * I12 - (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*I1;
      Sext.g(2, 2) = 1. / (1. - P->matrix->nu())*((a[0] * a[0] + a[1] * a[1]) / (2.*(a[0] + a[1])*(a[0] + a[1])) + (0.5 - P->matrix->nu()));
    }
    else{
      lambda = x[0] * x[0] + x[1] * x[1] - a[0] * a[0];
      I1 = I2 = 2.*PI*a[0] * a[0] / (a[0] * a[0] + lambda);
      I11 = I22 = I12 = PI*a[0] * a[0] / pow(a[0] * a[0] + lambda, 2);
      Sext.g(1, 1) = Sext.g(0, 0) = 1. / 8. / PI / (1. - P->matrix->nu())*((2. * P->matrix->nu()*I1 - I1 + a[0] * a[0] * I11) + 2. * (a[0] * a[0] * I11 - I1 + (1. - P->matrix->nu())*I1*2.));
      Sext.g(1, 0) = Sext.g(0, 1) = 1. / 8. / PI / (1. - P->matrix->nu())*((2. * P->matrix->nu()*I1 - I2 + a[0] * a[0] * I12));
      Sext.g(2, 2) = 2. / 8. / PI / (1. - P->matrix->nu())*(a[0] * a[0] * I12 - I2 + (1. - P->matrix->nu())*(I1 + I2));
    }
  }
  else{
    if (pointInside){
      lambda = 0.;
      I1 = 4.*PI*a[1] / (a[0] + a[1]);
      I2 = 4.*PI*a[0] / (a[0] + a[1]);
      I12 = 4.*PI / (a[0] * a[0] + 2.*a[0] * a[1] + a[1] * a[1]);
      I11 = 1. / 3. * (4. * PI / (a[0] * a[0]) - I12);
      I22 = 1. / 3. * (4. * PI / (a[1] * a[1]) - I12);
      Sext.g(0, 0) = 3. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * I11 + (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*I1;
      Sext.g(1, 1) = 3. / (8. * PI*(1. - P->matrix->nu()))*a[1] * a[1] * I22 + (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*I2;
      Sext.g(0, 1) = 1. / (8. * PI*(1. - P->matrix->nu()))*a[1] * a[1] * I12 - (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*I1;
      Sext.g(1, 0) = 1. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * I12 - (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*I2;
      Sext.g(2, 2) = 2. / (2. - 2.*P->matrix->nu())*((a[0] * a[0] + a[1] * a[1]) / (2.*(a[0] + a[1])*(a[0] + a[1])) + (1. - 2.*P->matrix->nu()) / 2.);
    }
    else{
      double b = a[0] * a[0] + a[1] * a[1] - x[0] * x[0] - x[1] * x[1];
      double c = a[0] * a[0] * a[1] * a[1] - x[0] * x[0] * a[1] * a[1] - x[1] * x[1] * a[0] * a[0];
      double Di = b*b - 4. * c;
      lambda = 0.5*(-b + sqrt(Di));
      I1 = 4.*PI*a[0] * a[1] / (a[1] * a[1] - a[0] * a[0])*(sqrt((a[1] * a[1] + lambda) / (a[0] * a[0] + lambda)) - 1.);
      I2 = 4.*PI*a[0] * a[1] / (a[0] * a[0] - a[1] * a[1])*(sqrt((a[0] * a[0] + lambda) / (a[1] * a[1] + lambda)) - 1.);
      I12 = 1.*(I2 - I1) / (a[0] * a[0] - a[1] * a[1]);
      I11 = 1. / 3. * ((I1 + I2) / (a[0] * a[0] + lambda) - I12);
      I22 = 1. / 3. * ((I1 + I2) / (a[1] * a[1] + lambda) - I12);
      Sext.g(0, 0) = 1. / 8. / PI / (1. - P->matrix->nu())*((2 * P->matrix->nu()*I1 - I1 + a[0] * a[0] * I11) + 2. * (a[0] * a[0] * I11 - I1 + (1. - P->matrix->nu())*(I1 + I1)));
      Sext.g(1, 1) = 1. / 8. / PI / (1. - P->matrix->nu())*((2 * P->matrix->nu()*I2 - I2 + a[1] * a[1] * I22) + 2. * (a[1] * a[1] * I22 - I2 + (1. - P->matrix->nu())*(I2 + I2)));
      Sext.g(0, 1) = 1. / 8. / PI / (1. - P->matrix->nu())*((2 * P->matrix->nu()*I1 - I2 + a[0] * a[0] * I12));
      Sext.g(1, 0) = 1. / 8. / PI / (1. - P->matrix->nu())*((2 * P->matrix->nu()*I2 - I1 + a[1] * a[1] * I12));
      Sext.g(2, 2) = 2. / 8. / PI / (1. - P->matrix->nu())*(a[0] * a[0] * I12 - I2 + (1. - P->matrix->nu())*(I1 + I2));
    }
  }
}

void InclusionRecord2D :: compute_S_int()
{
  double _I1,_I2,_I11,_I22,_I12;
  if (a[0] == a[1]){
    _I1 = _I2 = 2.*PI;
    _I11 = _I22 = _I12 = PI / (a[0] * a[0]);
    S[0] = S[3] = 3. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * _I11 + (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*_I1;
    S[1] = S[2] = 1. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * _I12 - (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*_I1;
    S[4]        = 1. / (1. - P->matrix->nu())*((a[0] * a[0] + a[1] * a[1]) / (2.*(a[0] + a[1])*(a[0] + a[1])) + (0.5 - P->matrix->nu()));
  }
  else{
    _I1 = 4.*PI*a[1] / (a[0] + a[1]);
    _I2 = 4.*PI*a[0] / (a[0] + a[1]);
    _I12 = 4.*PI / ((a[0] + a[1])*(a[0] + a[1]));
    _I11 = 1. / 3. * (4. * PI / (a[0] * a[0]) - _I12);
    _I22 = 1. / 3. * (4. * PI / (a[1] * a[1]) - _I12);
    S[0] = 3. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * _I11 + (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*_I1;
    S[3] = 3. / (8. * PI*(1. - P->matrix->nu()))*a[1] * a[1] * _I22 + (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*_I2;
    S[1] = 1. / (8. * PI*(1. - P->matrix->nu()))*a[1] * a[1] * _I12 - (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*_I1;
    S[2] = 1. / (8. * PI*(1. - P->matrix->nu()))*a[0] * a[0] * _I12 - (1. - 2. * P->matrix->nu()) / (8. * PI*(1. - P->matrix->nu()))*_I2;
    S[4] = 1. / (1. - P->matrix->nu())*((a[0] * a[0] + a[1] * a[1]) / (2.*(a[0] + a[1])*(a[0] + a[1])) + (0.5 - P->matrix->nu()));
  }
}
double InclusionRecord2D :: r(int i,int j)
{
  if(i==j)return 1.0;
  return 0.0;
}

double InclusionRecord2D :: compute_D(int i,int j,int k,int l)
{
  return 1.0/(8.0*PI*(1.0-P->matrix->nu()))*(
                  8.0*PI*(1.0-P->matrix->nu())*S_give(1000*i+100*j+10*k+l)+2.0*P->matrix->nu()*r(k,l)*x[i-1]*dI(i,j,0)
                  +(1.0-P->matrix->nu())*(  r(i,l)*x[k-1]*dI(k,j,0)  +  r(j,l)*x[k-1]*dI(k,i,0)  +  r(i,k)*x[l-1]*dI(l,j,0)   +   r(j,k)*x[l-1]*dI(l,i,0)   )
                  -r(i,j)*x[k-1]*(dI(k,l,0)-pow(a[i-1],2)*dI(10*k+i,l,0)) - (r(i,k)*x[j-1]+r(j,k)*x[i-1])*(dI(j,l,0)-pow(a[i-1],2)*dI(10*i+j,l,0))
                  -(r(i,l)*x[j-1]+r(j,l)*x[i-1])*(dI(j,k,0)-pow(a[i-1],2)*dI(10*i+j,k,0))-x[i-1]*x[j-1]*(dI(j,l,k)-pow(a[i-1],2)*dI(10*i+j,l,k))
                  );
}

void InclusionRecord2D :: compute_matrix_D()
{
  D.g(0, 0) = compute_D(1, 1, 1, 1);
  D.g(1, 1) = compute_D(2, 2, 2, 2);
  D.g(0, 1) = compute_D(1, 1, 2, 2);
  D.g(1, 0) = compute_D(2, 2, 1, 1);
  D.g(2, 2) = 2.*compute_D(1, 2, 1, 2) - Sext.g(2, 2);
  D.g(0, 2) = 2.*compute_D(1, 1, 1, 2);
  D.g(1, 2) = 2.*compute_D(2, 2, 1, 2);
  D.g(2, 0) = compute_D(1, 2, 1, 1);
  D.g(2, 1) = compute_D(1, 2, 2, 2);
}
void InclusionRecord2D::compute_eps_tau_back(vektor &e_s_, vektor &e_t_)
{
  pointInside = 1;
  x[0] = 0.0;
  x[1] = 0.0;
  compute_matrix_D();
  matice D_inv(3, 3);
  D_inv.inverze_3x3(D);
  vektor es(3);
  rotateStrain_G2L(e_s_.v, es.v);
  MaticexVektor(e_t_, D_inv, es);
}
void InclusionRecord2D :: compute_strain_pert_from_eps_zero_ext(const double *point, double *es, double *ez)
{
  pointInside = 1;
  x[0] = 0.0;
  x[1] = 0.0;
  compute_matrix_D();
  //int i;
  
  double m[5], n[5], Cs[5], B[5], e_z[3]; //, e_z_t[3];
  vektor vet(3),ves(3);
  SubtractTwoVectors(C, P->matrix->C, Cs, 5);
  
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(n, Cs, S);
  AddVector(P->matrix->C, n, 5);
  giveInverseMatrix3x3to5(m, n);
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(B, m, Cs);
  MultiplyVector(B, 5, -1.0);
  rotateStrain_G2L(ez, e_z);
  MatrixOperations::giveTVproduct_3is3x3to5and3(vet.v, B, e_z);
  
  pointInside = 0;
  this->transformCoords_G2L(point, x);
  compute_matrix_D();
  
  MaticexVektor(ves, D, vet);
  rotateStrain_L2G(ves.v, es);
}
void InclusionRecord2D::compute_strain_pert_from_eps_zero_int(const double *point, double *es, double *ez)
{
  pointInside = 1;
  x[0] = 0.0;
  x[1] = 0.0;
  compute_matrix_D();
  //int i;
  
  double m[5], n[5], Cs[5], B[5], e_z[3]; //, e_z_t[3];
  vektor vet(3), ves(3);
  SubtractTwoVectors(C, P->matrix->C, Cs, 5);
  
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(n, Cs, S);
  AddVector(P->matrix->C, n, 5);
  giveInverseMatrix3x3to5(m, n);
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(B, m, Cs);
  MultiplyVector(B, 5, -1.0);
  rotateStrain_G2L(ez, e_z);
  MatrixOperations::giveTVproduct_3is3x3to5and3(vet.v, B, e_z);
  
  MaticexVektor(ves, D, vet);
  rotateStrain_L2G(ves.v, es);
}

void InclusionRecord2D :: compute_from_eps_tau(const double *point, vektor &es, vektor &et)
{
  pointInside = 0;
  this->transformCoords_G2L(point, x);
  compute_matrix_D();

  MaticexVektor(es, D, et);
  rotateStrain_L2G(es.v, es.v);
}

void InclusionRecord2D :: compute_initial_eps_tau(int strainID)
{
  // termitovo:
  // zeptat se standy jestli je e_t v globalu nebo lokalu
  
  pointInside = 1;
  x[0]=0.0;
  x[1]=0.0;
  compute_matrix_D();

  double m[5], n[5], Cs[5], B[5], e_z[3];

  SubtractTwoVectors (C, P->matrix->C, Cs, 5);

  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(n, Cs, S);
  AddVector(P->matrix->C, n, 5);
  giveInverseMatrix3x3to5(m, n);
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(B, m, Cs);
  MultiplyVector(B, 5, -1.0);
  rotateStrain_G2L(P->matrix->remote_strains()[strainID], e_z);
  MatrixOperations::giveTVproduct_3is3x3to5and3(e_t.v, B, e_z);
}

void InclusionRecord2D :: update_eps_tau (int strainID, const double *e_z_add)
{
  pointInside = 1;
  x[0] = 0.0;
  x[1] = 0.0;
  compute_matrix_D();
  int i;
  double m[5], n[5], Cs[5], B[5], e_z[3],e_z_t[3];
  
  SubtractTwoVectors(C, P->matrix->C, Cs, 5);
  
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(n, Cs, S);
  AddVector(P->matrix->C, n, 5);
  giveInverseMatrix3x3to5(m, n);
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(B, m, Cs);
  MultiplyVector(B, 5, -1.0);
  for (i = 0; i < 3; i++)e_z_t[i] = P->matrix->remote_strains()[strainID][i]+ e_z_add[i];
  rotateStrain_G2L(e_z_t, e_z);
  MatrixOperations::giveTVproduct_3is3x3to5and3(Epsilon_Tau[strainID], B, e_z);
}

void InclusionRecord2D :: compute_eps_tau(int strainID,double *e_tau, double *e_z_add,bool add)
{
  // termitovo tady se pocita porad stejne D znovu pro vsechny body na inkluzi, nebo se mi tozda
  pointInside = 1;
  x[0] = 0.0;
  x[1] = 0.0;
  compute_matrix_D();
  int i;
  double m[5], n[5], Cs[5], B[5], e_z[3], e_z_t[3];
  
  SubtractTwoVectors(C, P->matrix->C, Cs, 5);
  
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(n, Cs, S);
  AddVector(P->matrix->C, n, 5);
  giveInverseMatrix3x3to5(m, n);
  MatrixOperations::giveTTproduct_3x3to5is3x3to5and3x3to5(B, m, Cs);
  MultiplyVector(B, 5, -1.0);
  if(add) for (i = 0; i < 3; i++)e_z_t[i] = P->matrix->remote_strains()[strainID][i] + e_z_add[i];
  else    for (i = 0; i < 3; i++)e_z_t[i] = e_z_add[i];
  rotateStrain_G2L(e_z_t, e_z);
  MatrixOperations::giveTVproduct_3is3x3to5and3(e_tau, B, e_z);
}

void InclusionRecord2D :: giveExtStrainPert (const double *point, double **strain, int lc, int nlc)
{
  pointInside = 0;
  this->transformCoords_G2L(point, x);
  compute_matrix_D();
  vektor epsTau(3);
  
  
  vektor v(3);
  for (int i = lc; i < lc + nlc; i++) {
    
    for (int j = 0; j < 3; j++)
      epsTau.v[j] = this->Epsilon_Tau[i][j];
    //new begin
    // now it is just a FAKE polynomial eigenstrain solution...
    if (P->approximation >= 1)
      for (int j = 0; j < P->give_VM_TENS_RANGE(); j++) {
        epsTau.v[j] += 1.0*(this->approx_coef[i][0][j] * x[0]) / ((1. + this->lambda));
        epsTau.v[j] += 1.0*(this->approx_coef[i][1][j] * x[1]) / ((1. + this->lambda));
        if(P->give_dimension()==3)epsTau.v[j] += 1.0*(this->approx_coef[i][2][j] * x[2]) / ((1. + this->lambda));
      }
    if (P->approximation >= 2) {
      if (P->give_dimension() == 2) {
        for (int j = 0; j < P->give_VM_TENS_RANGE(); j++) {
          epsTau.v[j] += 1.0*(this->approx_coef[i][2][j] * x[0] * x[0]) / (SQR(1. + this->lambda));
          epsTau.v[j] += 1.0*(this->approx_coef[i][3][j] * x[1] * x[1]) / (SQR(1. + this->lambda));
          epsTau.v[j] += 1.0*(this->approx_coef[i][4][j] * x[0] * x[1]) / (SQR(1. + this->lambda));
        }
      }
      else {
        
      }
      
    }
    
    
    
    
    //new end
    MaticexVektor(v, D, epsTau);
    rotateStrain_L2G(v.v, strain[i - lc]);
  }
}

void InclusionRecord2D::getDisplacement(const double *point, double **displacement, int lc, int nlc, int _pointInside)
{
  pointInside = _pointInside;
  this->transformCoords_G2L(point, x);

  double _2nu = 2.0 * P->matrix->nu();
  double _1min2nu = 1. - _2nu;
  double multTRN = 1. / (8. * PI * (1. - P->matrix->nu()));
  
  L.g(0, 0) = multTRN * (x[0] * (_1min2nu * I(1) + 3. * a[0] * a[0] * I(11)));
  if (!pointInside)L.g(0, 0) += multTRN * x[0] * x[0] * (-dI(1, 1, 0) + a[0] * a[0] * dI(11, 1, 0));
  L.g(1, 1) = multTRN * (x[1] * (_1min2nu * I(2) + 3. * a[1] * a[1] * I(22)));
  if (!pointInside)L.g(1, 1) += multTRN * x[1] * x[1] * (-dI(2, 2, 0) + a[1] * a[1] * dI(22, 2, 0));
  L.g(0, 1) = multTRN * (x[0] * (_2nu     * I(1) - I(2) + a[0] * a[0] * I(12)));
  if (!pointInside)L.g(0, 1) += multTRN * x[0] * x[1] * (-dI(2, 2, 0) + a[0] * a[0] * dI(12, 2, 0));
  L.g(1, 0) = multTRN * (x[1] * (_2nu     * I(2) - I(1) + a[1] * a[1] * I(21)));
  if (!pointInside)L.g(1, 0) += multTRN * x[0] * x[1] * (-dI(1, 1, 0) + a[1] * a[1] * dI(21, 1, 0));
  L.g(0, 2) = 2.*multTRN * (x[1] * (_1min2nu * I(2) + a[0] * a[0] * I(12)));
  if (!pointInside)L.g(0, 2) += 2.*multTRN * x[0] * x[0] * (-dI(1, 2, 0) + a[0] * a[0] * dI(11, 2, 0));
  L.g(1, 2) = 2.*multTRN * (x[0] * (_1min2nu * I(1) + a[1] * a[1] * I(21)));
  if (!pointInside)L.g(1, 2) += 2.*multTRN * x[0] * x[1] * (-dI(1, 2, 0) + a[1] * a[1] * dI(21, 2, 0));
  
  vektor vect(2), etset(3);
  for (int i = lc; i < lc + nlc; i++) {
    for (int j = 0; j < 3; j++)
      etset.v[j] = this->Epsilon_Tau[i][j];
    
    MaticexVektor(vect, L, etset);
    this->rotateDisplc_L2G (vect.v, displacement[i - lc]);
  }
}
void InclusionRecord2D :: computeVolume (void)
{
  switch( shape ) {
  case IS_ELLIPSE:
  case IS_CIRCLE:
    this->volume = PI * a[0] * a[1];
    break;
  //case IS_STRIP:
  //case IS_CRACK:
  //  _errorr( "Not yet defined" );
  //  break;
  default:
    _errorr( "Undefined shape" );
  }
}

// ***********************************************************************************************************
// ****************   2D ONLY -  END    **********************************************************************
// ***********************************************************************************************************


} // end of namespace mumech

/*end of file*/