inclusion.h

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           inclusion.h
// 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.
//********************************************************************************************************

#ifndef MUMECH_INCLUSION_H
#define MUMECH_INCLUSION_H

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


using namespace gelibspace;

namespace mumech {

class Problem;
class Point;

class Inclusion
{
  friend class Homogenization;

 public:
  //
  // INPUT DATA
  // 
  long id;                    
  const Problem *P;           
  
  InclusionGeometry shape;    
  double E;                   
  double nu;                  
  
  double *origin;             
  double *a;                  
  double *eAngles;            
  
  
  //
  // *** CHARACTERISTIC VALUES ***
  //
  bool rotated;
  double volume;

  double *C;                  
  double *S;                  
  double *SInv;               

  // termitovo: TInv = T^T - pak je urcite zbytecny TInv ukladat do rec file, ztrata mista
  // TInv se pouziva na jednom miste jako TInv*vector, muze se udelat fce T *^T vector, tim padem se pouzije T a TInv je uplne zbytecne
  double *T;                  
  double *TInv;               

  // Strain and stress reduced vectors are in theoretical notation => trasformation matrices are same.
  double *Te;                 
  double *TeInv;              
  
 public:
  
  //
  // PRE-COMPUTED DATA, when AF->char = true
  //
  // inclusion data structure - inclusion record previously calculated
  
  double actionRadius;           
  double SQRactionRadius;        
  double ndiff_1;                
  double ndiff_2;                
  
  int noActingIncls;             
  int * actingIncls;             
  
  // termitovo Epsilon_Tau by se casem melo ztotoznit s asi locEigStrain_LC, ale opatrne
  double **Epsilon_Tau;          
  // approximation
  // toto je tu pro pocitani nekonstantnich poli v inkluzi, jsou to integracni/approximacni body
  // nejsou deallokovany
  int n_approx_points;
  double **approx_points;
  double **approx_points_eps_tau;
  double ***approx_coef;
  
  double **globPert_displc;      
  double **globPert_strain;      
  double **globPert_stress;      
  
 public:
  Inclusion (long i, const Problem *p);
  virtual ~Inclusion();
  
  bool scan_eAngles_RAD (Stream *stream);
  bool scan_eAngles_DEG (Stream *stream);
  virtual bool scan_locEigStrain_LC (Stream *stream, int lc) = 0;

  virtual void allocate_nlc_fields (void);
  
  virtual InclusionGeometry giveInclusionGeometry (InclusionGeometry shp_o) const = 0;
  virtual void input_data_initialize_and_check_consistency (void);
  virtual void initialize (const Inclusion *prevInclRec) = 0;
  virtual void SBA_computeInitialStrains (int lc) = 0;
  virtual bool point_is_affected (const double *point) const = 0;
  
  bool point_is_inside (const double *coords, double epsilon = 0.0) const;
  
  void transformCoords_G2L (const double *glob, double *loc) const;
  void transformCoords_L2G (const double *loc, double *glob) const;
  
  void rotateStrain_G2L (const double *glob, double *loc) const;
  void rotateStrain_L2G (const double *loc, double *glob) const;
  void rotateDisplc_L2G (const double *loc, double *glob) const;

  void ActingIncls_allocate (void);

  int find_overlap (const Inclusion *incl) const;
  bool is_inside_of_BB (const double *bb1, const double *bb2) const;
  
  // DIRECT API
  void set_centroids (double x, double y);
  void set_centroids (double x, double y, double z);
  void set_Inclusion_shape (InclusionGeometry val) { shape = val; }
  void set_Youngs_modulus (double val) { E = val; }
  void set_Poissons_ratio (double val) { nu = val; }
  void set_Semiaxes_dimensions (double x, double y);
  void set_Semiaxes_dimensions (double x, double y, double z);
  void set_Euller_angles_deg (double x);
  void set_Euller_angles_deg (double x, double y, double z);
  //void set_Remote_strains_for_lc (int lc, const double *r1);

  // agregates
  void set_all_2d (double x, double y, double e, double n,
                   double a1, double a2, double e1);
  void 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 );



  double give_volume( void );

  //double** pointer_globPertStrain(void) { return this->globPert_strain; };
  //double** pointer_globPertStress(void) { return this->globPert_stress; };
  //double** pointer_globPertDisplc(void) { return this->globPert_displc; };
  
  void update_approximations(int lc);
  
  void Eps02EpsTau(double *e_tau, const double *e_0);
  
 protected:
  void EAdeg2rad (void);
  
  void giveTransformationStrainOperator (double oper[12], const double C[12], const double C1[12], const double S[12]);
  void give_EshelbyMatrixFull (double** answer) const;
  void give_EshelbyMatrixReduced (double* answer) const;
 public:
  void give_StiffnessMatrixFull (double* answer) const;
  void give_StiffnessMatrixReduced (double* answer) const;
 protected:
  void give_TeMatrix_G2L (double* answer) const;
  void give_TeMatrix_L2G (double* answer) const;

 public:
  const double** give_EshelbyPertDisplc_internal (                 int lc, int nlc, const double *coords);
  //void           give_EshelbyPertDisplc_internal (double **displc, int lc, int nlc, const double *coords);
  const double** give_EshelbyPertStrain_internal (                 int lc, int nlc);
  //void           give_EshelbyPertStrain_internal (double **strain, int lc, int nlc);
  const double** give_EshelbyPertStress_internal (                 int lc, int nlc);
  //void           give_EshelbyPertStress_internal (double **stress, int lc, int nlc);
  
  //void           add_EshelbyPertDisplc_internal_SIFCM(double **displc, int lc, int nlc, const double *coords, double **e_t_SIFCM);
  void           add_EshelbyPertStrain_internal_SIFCM(double **strain, int lc, int nlc, double **e_t_SIFCM);
  void           add_EshelbyPertStress_internal_SIFCM(double **stress, int lc, int nlc, double **e_t_SIFCM, double **e_s_ext_add);
  
  void give_EshelbyPertFields_external (Point *point, int lc, int nlc, bool disp, bool strn);

  double compute_supplement_energy (int lc);

private:
  virtual void computeVolume() = 0;


};//end of class declaration

// ***********************************************************************************************************
// ***                                         3D INCLUSION RECORD                                         ***
// ***********************************************************************************************************
class eshelbySoluEllipticIntegrals;
class eshelbySoluUniformField;

class InclusionRecord3D : public Inclusion
{
 public:
  // 
  // *** CLASSES ***
  //
  eshelbySoluEllipticIntegrals *ellInt;
  eshelbySoluUniformField      *esuf;
  
  //
  // *** CHARACTERISTIC MATRICES ***
  //
  double eInt[13];             

  //
  // PRE-COMPUTED DATA - previously calculated inclusion record, when P->data_equivalent == true
  //
  double locEigStrain[6];        
  double locEigStrainTot[6];     
  double globEigStrainTotal[6];  
  
  double **locEigStrain_LC;      
  double **globEigStrain_LC;     
  
 public:
  InclusionRecord3D (long i, const Problem *p);
  virtual ~InclusionRecord3D();
  
  bool scan_locEigStrain_LC (Stream *stream, int lc);
  bool scan_globEigStrain_LC (Stream *stream, int lc);
  
  void allocate_nlc_fields (void);
  
  InclusionGeometry giveInclusionGeometry (InclusionGeometry shp_o) const;
  void input_data_initialize_and_check_consistency (void);
  void initialize (const Inclusion *prevInclRec);
  void SBA_computeInitialStrains (int lc);
  bool point_is_affected (const double *point) const;

  void SBA_updateGlobAndLocStrains (Point *point);

  void addtot ();
  
  // *** PRINT ***
  //void init_locEigStrain_Singl (NotationType notationFlag);
  //void init_locEigStrain_Multi (NotationType notationFlag);
  
  void init_EigStrain_LC (int lc);


 private:
  virtual void computeVolume();

};



// ***********************************************************************************************************
// ****************   2D ONLY - START   **********************************************************************
// ***********************************************************************************************************
class bod_pomocny{
public:
    double x[2];
    double a[2];
    double lambda;
    double I1;
    double I2;
    double I11;
    double I12;
    double I22;
    
    double I(int index){
    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;
  }
    void spocitat(){
    if (a[0] == a[1]){
      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] / (a[0] * a[0] * a[0] * a[0] + 2.*a[0] * a[0] * lambda + lambda*lambda);
    }
    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] - a[0] * a[0] * x[1] * x[1];
      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])*(pow(a[1] * a[1] + lambda, 0.5) / pow(a[0] * a[0] + lambda, 0.5) - 1.);
      I2 = 4.*PI*a[0] * a[1] / (a[0] * a[0] - a[1] * a[1])*(pow(a[0] * a[0] + lambda, 0.5) / pow(a[1] * a[1] + lambda, 0.5) - 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);
    }
    }
};




class InclusionRecord2D : public Inclusion
{
 public:
  // spocitano - vybalancovano
  vektor delta_e_t;
  vektor delta_e_s;
  
  // body
  double x[2];
  double x_last_derivation[2];
  double x_last_IS[2];
  double lambda;
  double I1;
  double I2;
  double I11;
  double I12;
  double I22;
  
  matice Sext;
  matice D;
  matice L;
  
  vektor e_t;                     // epsilon tau
  
  
  bod_pomocny P_a[2];             // auxiliary points for the first derivatives
  bod_pomocny P_b[8];             // auxiliary points for the second derivatives
  
  int pomocne_body_spocitane;
  double h,H;                     // dx pri derivaci
  bool any_derivative_preparation_computed;
  bool any_IS_computed;
  int pointInside;
  
 public:
  InclusionRecord2D (long i, const Problem *p);
  virtual ~InclusionRecord2D();
  
  bool scan_locEigStrain_LC (Stream *stream, int lc);
  
  void allocate_nlc_fields (void);
  
  InclusionGeometry giveInclusionGeometry (InclusionGeometry shp_o) const;
  void input_data_initialize_and_check_consistency (void);
  void initialize (const Inclusion *prevInclRec);
  void SBA_computeInitialStrains (int lc);
  bool point_is_affected (const double *point) const;


  double S_give(int index);
  double I(int index);
  double derivative_step(double delitel);
  void derivative_preparation();
  double dI(int index,int smer1,int smer2);
  void compute_I_S();
  void compute_S_int();
  double r(int i,int j);
  double compute_D(int i,int j,int k,int l);
  void compute_matrix_D();
  void compute_eps_tau_back(vektor &e_s_, vektor &e_t_);
  void compute_strain_pert_from_eps_zero_int(const double *point, double *es, double *ez);
  void compute_strain_pert_from_eps_zero_ext(const double *point, double *es, double *ez);
  void compute_from_eps_tau(const double *point, vektor &es, vektor &et);
  void giveExtStrainPert (const double *point, double **strain, int lc, int nlc);
  void compute_initial_eps_tau(int strainID);
  
  void compute_eps_tau(int strainID, double *e_tau, double *e_z_add,bool add);
  void update_eps_tau(int strainID, const double *e_z_add);
  
  void getDisplacement(const double *point, double **displacement, int lc, int nlc, int _pointInside);

 private:
  virtual void computeVolume();

};

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

} // end of namespace mumech

#endif // end of MUMECH_INCLUSION_H

/*end of file*/