problem.h

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           problem.h
// author(s):      Ladislav Svoboda, Jan Novak
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************
// *******************************************************************************************************
//  TADY BUDU ZATIM DRZET OBECNE INFORMACE O KODU
//  a) FORMATY
//     TENSOR - symetricke tensory druheho radu, 2x2 pro 2d a 3x3 pro 3d, hlavne strain a stress
//              - interne jsou ukladany (vetsinou) v 1d poli v TNTV_THEORETICAL_FEEP notaci
//              - I/O, tj. pri komunikaci s uzivatelem, jsou defaultne v 1d poli v TNTV_THEORETICAL_ROW notaci
//
//
// *******************************************************************************************************




#ifndef MUMECH_PROBLEM_H
#define MUMECH_PROBLEM_H

#include "types.h"
#include "macros.h"
#include "matrix.h"
#include "inclusion.h"



namespace mumech {

class Inclusion;
class Mesh;
class MatrixRecord;
class Homogenization;


class Problems
{
 protected:
  // INPUT VARIABLES
  bool twodim;                 
  
 public:
  Problems () {}
  virtual ~Problems () {}
  
  virtual classID give_classid() const = 0; //  { errol; return classVoid; } ///= 0; = 0;
  
  
  virtual bool is_converted_to_equivalent (void) const = 0;
  virtual int give_nLC (void) const = 0;
  
  int give_twodim (void) const { return twodim; }
  int give_dimension (void) const { return twodim ? 2 : 3; }
  
  int ndisplc_one (void) const { return this->give_VECT_RANGE()    ; }
  int nstrain_one (void) const { return this->give_VM_TENS_RANGE() ; }
  
  int give_VECT_RANGE    (void) const { return (twodim ? 2 : 3); }
  int give_EA_RANGE      (void) const { return (twodim ? 1 : 3); }
  int give_TENS_RANGE    (void) const { return (twodim ? 4 : 9); }
  int give_VM_TENS_RANGE (void) const { return (twodim ? 3 : 6); }
  int give_ISO_C_RANGE   (void) const { return (twodim ? 5 : 12); }
  int give_EL_POT_RANGE  (void) const { return (twodim ? 0 : 13); }
  int give_TRNSFM_MTRX_VEC_RANGE  (void) const { return (twodim ? 4 :  9); }
  int give_TRNSFM_MTRX_TENS_RANGE (void) const { return (twodim ? 9 : 36); }
  
  void set_dimension (int d)
  {
    if      (d == 2) twodim = true;
    else if (d == 3) twodim = false;
    else _errorr ("wrong dimension number, choose 2 or 3");
  }
  
};

class ProblemFEM : public Problems
{
private:
  int nLC;
  
public:
  ProblemFEM (int d, int n) : Problems () {
    this->set_dimension(d);
    this->nLC = n;
    
  }
  virtual ~ProblemFEM () {}
  
  virtual classID give_classid() const { return classProblemFEM; }
  
  virtual bool is_converted_to_equivalent (void) const { return true; }
  virtual int give_nLC (void) const { return nLC; }
 
  
};

class Problem : public Problems
{
  friend class Inclusion;
  friend class InclusionRecord3D;
  friend class InclusionRecord2D;
  friend class printBenchmarkResultsVtk;
  friend class selfBalanceAlgorithm;
  friend class eshelbySoluUniformField;
  friend class eshelbySoluEllipticIntegrals;
  friend class MatrixRecord;
  friend class Homogenization;
  friend class Dilute;
  friend class DifferentialScheme;
  friend class RegGrid;
  friend class Mesh;
  friend class Point;
  friend class mNode;
  friend class mElement;
  
  
  // /** All variables are private. */
  // ***************************************************************************************************
  // /**   ***   @name   ***                   PRIVATE VARIABLES                   ***   */*   ***   ***
  // *************************************************-*************************************************
  // @{
 private:
  // STATE VARIABLES
  bool data_set;               
  bool data_equivalent;        
  
  // INPUT VARIABLES
  //bool twodim;                 ///<  2 dimension problem; 3d is default; twodim == true - 2d, twodim == false - 3d
  
  MatrixRecord* matrix;        
  int noIncl;                  
  Inclusion **inclusions;      
  
  SBAtype SBA_type;            
  bool SBA_optimized;          
  int SBA_maxiters;            
  int SBA_reqiters;            
  
  int intFieldsShape;          
  int approximation;           
  double approx_point_pos1;
  
  DiffTypes diffType;          
  
  double problem_size_A[3];    
  double problem_size_B[3];    
  double node_dist;            
  
  int verbose;                 
  
  double min_axes_diff;        
  double max_axes_diff;        
  
  bool min_axes_diff_change;   
  bool max_axes_diff_change;   
  
  
  // COMPUTED VARIABLES
  double volumeOfIncls;        
  int totalNoActingIncls;      
  
  //
  int nmeshes;                 
  Mesh **meshes;               
  
  int nhomogs;                 
  Homogenization **homogs;     
  // @}
  
  
  // ***************************************************************************************************
  // /**   ***   @name   ***                   PUBLIC FUNCTIONS                    ***   */*   ***   ***
  // *************************************************-*************************************************
  // @{
 public:
  Problem (void);
  virtual ~Problem ();
  
  virtual classID give_classid() const { return classProblem; }
  
  void read_input_file (const char *filename);

  void input_data_initialize_and_check_consistency (void);
  
  void convert_to_equivalent_problem (void);
  
  void print_equivalent_problem (const char *filename);
  
  long giveFieldsOfPoint (double **displc, double **strain, double **stress, const double *coords,
                          char ptFlag, int rs, int nrs,
                          PFCmode pfcMode=PFCM_OPTIMIZED, long reqIncl = -3, STRNotation tn = STRN_THEORETICAL_ROW) const;

  long giveFieldsOfPointOneRS (double *displc, double *strain, double *stress, const double *coords,
                               char ptFlag, int rs,
                               PFCmode pfcMode=PFCM_OPTIMIZED, long reqIncl = -3, STRNotation tn = STRN_THEORETICAL_ROW) const;
  
  void printFieldsOnMeshVTK (const char *mesh_file_out, const char *mesh_file,
                             char ptFlag, int rs, int nrs,
                             PFCmode pfcMode=PFCM_OPTIMIZED) const;

  void printFieldsOnMeshGrid (const char *mesh_file_out, const double *p1, const double *p2, const long *n,
                              char ptFlag, int rs, int nrs,
                              PFCmode pfcMode=PFCM_OPTIMIZED) const;

  void print_visualization (const char *filename, int n, int dim=0, bool refined=false);

  void check_overlap (void);

  void matrix_giveReducedStiffMatrix (double *m) const { this->give_matrix()->giveReducedStiffMatrix(m); }
  void matrix_giveFullStiffMatrix    (double *m) const { this->give_matrix()->giveFullStiffMatrix(m); }


  virtual int give_nLC (void) const { if (matrix->nRS==0) _errorr("Number of load cases have not been set");  return matrix->nRS; }
  DiffTypes give_diffType (void) const { return diffType; }
  //SBAmode give_SBA_mode (void) { return SBA_mode; }

  // int set_verbose (void) const { return verbose; }

  Mesh* give_new_mesh(void);
  Homogenization* give_new_homogenization (HomogenizationType ht);
  
  // AGREGATED FUNCTIONS
  
  Problem* read_inclusions_plus_initialize_and_print (const char *inclusion_file, const char *equiv_file, DiffTypes dt);
  
  // @}
  
  
  
  // /** Public functions, obsolete or debugging functions. */
  // ***************************************************************************************************
  // /**   ***   @name   ***         PUBLIC FUNCTIONS - obsolete debug             ***   */*   ***   ***
  // *************************************************-*************************************************
  // @{
 public:
  //
  void set_SBA_optimized                 (bool val) { SBA_optimized  = val; }
  void set_SBA_maximumNumberOfIterations  (int val) { SBA_maxiters   = val; }
  void set_SBA_requiredNumberOfIterations (int val) { SBA_reqiters   = val; }
  void set_intFieldsShape                 (int val) { intFieldsShape = val; }
  void set_SBAM (SBAtype val) { SBA_type = val; }
  
  void set_approximation(int val);
  
  void set_problem_size(double a1, double a2, double a3, double b1, double b2, double b3, double node_distance) {
    problem_size_A[0] = a1;
    problem_size_A[1] = a2;
    problem_size_A[2] = a3;
    problem_size_B[0] = b1;
    problem_size_B[1] = b2;
    problem_size_B[2] = b3;
    node_dist = node_distance;
  }
  
  /*void give_inside_strain_of_incl0_on_other_affect(double *strain_others,double *strain) {
    ((InclusionRecord2D*)inclusions[0])->update_eps_tau(0, strain_others);
    
    double **str = AllocateArray2D(1, give_VM_TENS_RANGE());
    for (int i = 0; i < give_VM_TENS_RANGE(); i++)str[0][i] = 0.;
    inclusions[0]->give_EshelbyPertStrain_internal(str, 0, 1);
    for (int i = 0; i < give_VM_TENS_RANGE(); i++)strain[i] = str[0][i]+ strain_others[i];
    }*/
  
  
  void set_semiaxes_min_difference (double val) { min_axes_diff = val; }
  void set_semiaxes_max_difference (double val) { max_axes_diff = val; if (val > INFTY)  _errorr2("Set max_axes_diff greather then mumech infinity value %lf", INFTY); }

  Inclusion* give_inclusion (long i) const { return inclusions[i]; }
  long give_inclusion_with_point_inside (const double *coords, double epsilon = 0.0) const;
  void give_ovlivneni (const double *coords, double *result);
  void compute_fields_and_cmp_with_FEM_on_mesh (const char *mesh_file, const char *fem_results, int lc, int nlc, PFCmode pfcMode);
  double compute_supplement_energy (int lc);

  // @}
  
  
  
  // /** Public functions for direct data input, obligatory for correct input. */
  // ***************************************************************************************************
  // /**   ***   @name   ***            DIRECT INPUT API - obligatory              ***   */*   ***   ***
  // *************************************************-*************************************************
  // @{
 public:
  void set_data_set (void) { data_set = true; }
  //void set_dimension (int d);
  void set_matrix_E_nu (double E, double nu) { this->matrix->set_E_nu(E, nu); }
  void set_numberOfRemoteStrains (int n) { this->matrix->set_nlc(n); }
  void set_RemoteStrain (int id, const double *rs) { this->matrix->set_Remote_strain_for_lc(id, rs); }
  void set_number_of_inclusions (long n);
  
  void set_inclusion_centroids          (long id, double x, double y, double z=0.0)    { if (twodim)  this->inclusions[id]->set_centroids (x,y);  else  this->inclusions[id]->set_centroids (x,y,z); }
  void set_inclusion_E_nu               (long id, double E, double nu)                 { this->inclusions[id]->set_Youngs_modulus(E); this->inclusions[id]->set_Poissons_ratio(nu); }
  void set_inclusion_semiaxesDimensions (long id, double a1, double a2, double a3=0.0) { if (twodim)  this->inclusions[id]->set_Semiaxes_dimensions (a1,a2);  else  this->inclusions[id]->set_Semiaxes_dimensions (a1,a2,a3); }
  void set_inclusion_EulerAnglesDEG (long id, double e1, double e2=0.0, double e3=0.0) { if (twodim)  this->inclusions[id]->set_Euller_angles_deg (e1);  else  this->inclusions[id]->set_Euller_angles_deg (e1,e2,e3); }
  // @}
  
  
  // /** Public functions for direct input of the voluntary data, default values are used if not set. */
  // ***************************************************************************************************
  // /**   ***   @name   ***            DIRECT INPUT API - voluntary              ***   */*   ***   ***
  // *************************************************-*************************************************
  // @{
 public:
  void set_inclusion_shape (long id, InclusionGeometry shape) { this->inclusions[id]->set_Inclusion_shape(shape); }
  void set_diffType (DiffTypes val) { diffType = val; }
  void set_DataDimDiff (int dim, DiffTypes dt) { this->set_data_set(); this->set_dimension(dim); this->set_diffType(dt); }
  void set_UnitRemoteStrains (void) { this->matrix->set_unit_remote_strain(); }
  void set_inclusion_all (long id, double x, double y, double e, double n, double a1, double a2, double e1)
  { this->inclusions[id]->set_all_2d(x, y, e, n, a1, a2, e1); }
  void set_inclusion_all (long id, double x, double y, double z, double e, double n, double a1, double a2, double a3, double e1, double e2, double e3)
  { this->inclusions[id]->set_all_3d(x, y, z, e, n, a1, a2, a3, e1, e2, e3); }

  // allocate and generate a single 2d inclusion 2d problem
  void generate_equiv_1x2dI_2d (double Em, double num, const double *rs, double x, double y, double E, double nu, double a1, double a2, double ea, DiffTypes dt);
  // @}
  
  
  
  // /**  */
  // ***************************************************************************************************
  // /**   ***   @name   ***                  PRIVATE FUNCTIONS                    ***   */*   ***   ***
  // *************************************************-*************************************************
  // @{
 private:
  
  virtual bool is_converted_to_equivalent (void) const { return data_equivalent; }
  //int give_twodim (void) const { return twodim; }
  int is_twodim (void) const { return twodim; }
  //virtual int give_dimension (void) const { return twodim ? 2 : 3; }
  int get_SBA_maxiters (void) { return SBA_maxiters; }
  int get_SBA_reqiters (void) { return SBA_reqiters; }

  void check_dim (int i) const { if ( (i != 2 && i != 3) || (i == 2 && twodim==false) || (i == 3 && twodim==true)) _errorr("problem dimension error"); }
  
  int give_verbose (void) const { return verbose; }
  double give_semiaxes_min_difference (void) const { return min_axes_diff; }
  double give_semiaxes_max_difference (void) const { return max_axes_diff; }
  bool give_semiaxis_min_difference_change (void) const { return min_axes_diff_change; }
  bool give_semiaxis_max_difference_change (void) const { return max_axes_diff_change; }
  MatrixRecord* give_matrix (void) const { return matrix; }
  
  // /// Gives vector range, same as dimension.
  // int give_VECT_RANGE    (void) const { return (twodim ? 2 : 3); }
  // /// Gives vector range, same as dimension.
  // int give_EA_RANGE      (void) const { return (twodim ? 1 : 3); }
  // /// Gives tensor range, same as dimension^2. Theoretical range of a second order tensor.
  // int give_TENS_RANGE    (void) const { return (twodim ? 4 : 9); }
  // /// Gives range of a second order tensor in Voigt-Mandel notation.
  // int give_VM_TENS_RANGE (void) const { return (twodim ? 3 : 6); }
  // /// Gives range of ...
  // int give_ISO_C_RANGE   (void) const { return (twodim ? 5 : 12); }
  // /// Gives range of ...
  // int give_EL_POT_RANGE  (void) const { return (twodim ? 0 : 13); }
  // /// Gives range of ...
  // int give_TRNSFM_MTRX_VEC_RANGE  (void) const { return (twodim ? 4 :  9); }
  // /// Gives range of ...
  // int give_TRNSFM_MTRX_TENS_RANGE (void) const { return (twodim ? 9 : 36); }
  
  
  //int ndisplc_one (void) const { return this->give_VECT_RANGE()    ; }
  //int nstrain_one (void) const { return this->give_VM_TENS_RANGE() ; }
  int ndisplc_all (void) const { return this->give_VECT_RANGE()    * matrix->nRS; }
  int nstrain_all (void) const { return this->give_VM_TENS_RANGE() * matrix->nRS; }

  bool file_type_s2e (const char *s);
  
  int check_lc_nlc (int lc, int nlc) const { return this->check_lc_nlc (lc, nlc, matrix->nRS); }
  static int check_lc_nlc (int lc, int nlc, int nRS)
  {
    if      (nlc <  0)  _errorr("small nlc");
    else if (nlc == 0)  nlc = nRS;
    else if (nlc > nRS) _errorr("big nlc");
    if (lc < 0)  _errorr("small lc");
    else if (lc+nlc > nRS)
      _errorr("big lc+nlc");
    return nlc;
  }
  
  PointPositionFlag givePointPosition (const double loc_x[3], const double sort_a[3], InclusionGeometry shape) const;
  
  void printVtkFileCompleteInclRec (const char *rsltFileName);
  
  void findAffectedInclusions (void);
  void updateEigenstrainsBySBalgorithm (void);
  void updateEigenstrainsBySBalgorithm_2D (int strainID);
    
  // 
  void updateStrainsInInclRecord (int noUpdate);
  //void giveTransfEigenstrains();
  
  long give_EshelbyPertFieldsOnePoint (const double *coords, double **disp, double **strain, double **stress,
                                       int rs, int nrs, PFCmode pfcMode, long reqIncl = -3) const;
  
  //long giveEshelbyTotalFields_OnePoint (const double *coords, double **disp, double **strain, double **stress,
  //                                      int rs, int nrs, PFCmode pfcMode, long reqIncl = -3, TensorNotationTypeInVector tn = TNTV_THEORETICAL_ROW) const;
  
  void give_EshelbyPertFieldsOnePoint_external (const double *coords, double **disp, double **strain, double **stress,
                           int lc, int nlc, PFCmode pfcMode, long parent_incl) const;
  
  long giveActingInclusions (long * actIncls, const double *coords) const;
  
  
  double updateEpsTauInSBal(int strainID, double **&last_eps_tau);
  
  // @}
  
}; // end of class declaration

} // end of namespace mumech

#endif

/*end of file*/