problem.cpp

Go to the documentation of this file.

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

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

//included libraries
#include "problem.h"

#include "matrixOperations.h"
//#include "elementaryFunctions.h"
#include "macros.h"
#include "inclusion.h"
#include "transformTensors.h"
#include "esuf.h"
#include "vtk.h"
#include "mesh.h"
#include "homogenization.h"
#include "mnode.h"
#include "melement.h"
#include "selfBalanceAlgorithm.h"
//#include "matrix.h"
#include "homogenizationMethods.h"

#include "tixy2.h"
#include "gelib.h"
#include "librw.h"
#include "mathlib.h"

#include <stdio.h>
#include <math.h>


using namespace tinyxml2;

namespace mumech {


//********************************************************************************************************
//                                            PUBLIC FUNCTIONS
//********************************************************************************************************
#define maxNumberOfMeshes 1000
Problem :: Problem (void) : Problems ()
{
  data_set = false;
  data_equivalent = false;
  
  twodim = false;
  
  matrix = new MatrixRecord (this);
  noIncl = 0;
  inclusions = NULL;
  
  SBA_optimized = true;
  SBA_maxiters =  999;
  SBA_reqiters = 1000;
  
  intFieldsShape = 2; // now 2 is default -> all the tests give the same results
  approximation = 0;
  approx_point_pos1 = 0.99;
  
  SBA_type = SBAT_VOID;

  diffType = DT_NUMERICAL;

  CleanVector(problem_size_A, 3);
  CleanVector(problem_size_B, 3);
  node_dist = 0.;
  
  verbose = 0;
  
  min_axes_diff = _MIN_AXES_DIFF;
  max_axes_diff = _MAX_AXES_DIFF;
  
  min_axes_diff_change = true;
  max_axes_diff_change = false;
  
  volumeOfIncls = 0.0;
  totalNoActingIncls = 0;
  
  nmeshes = 0;
  meshes = new Mesh*[maxNumberOfMeshes];
  memset (meshes, 0, maxNumberOfMeshes * sizeof(Mesh*));
  
  nhomogs = 0;
  homogs = new Homogenization*[maxNumberOfMeshes];
  memset (homogs, 0, maxNumberOfMeshes * sizeof(Homogenization*));
}

Problem :: ~Problem ()
{
  for (int i=0; i<noIncl; i++)  delete inclusions[i];
  delete [] inclusions;
  
  if (meshes) { for (int i=0; i<maxNumberOfMeshes; i++)  delete meshes[i];  delete [] meshes; }
  if (homogs) { for (int i=0; i<maxNumberOfMeshes; i++)  delete homogs[i];  delete [] homogs; }
  
  delete matrix;
} // end of function: destructor

void Problem :: input_data_initialize_and_check_consistency (void)
{
  for (int i=0; i<noIncl; i++)
    inclusions[i]->input_data_initialize_and_check_consistency();
  
}


void Problem :: set_number_of_inclusions (long n)
{
  noIncl = n; // number of inclusions
  if (twodim) { inclusions = new Inclusion*[noIncl];  for (int i=0; i<noIncl; i++) inclusions[i] = new InclusionRecord2D(i, this); }
  else        { inclusions = new Inclusion*[noIncl];  for (int i=0; i<noIncl; i++) inclusions[i] = new InclusionRecord3D(i, this); }
}




// allocate and generate a single 2d inclusion 2d problem
void Problem :: 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)
{
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  this->set_data_set();
  this->set_dimension(2);
  this->set_diffType(dt);
  this->set_matrix_E_nu(Em, num);
  if (rs) {
    this->set_numberOfRemoteStrains(1);
    this->set_RemoteStrain(0, rs);
  }
  else
    this->set_UnitRemoteStrains();
  
  // inclusions
  this->set_number_of_inclusions(1);
  // first inclusion
  this->set_inclusion_centroids(0, x, y);
  this->set_inclusion_E_nu(0, E, nu);
  this->set_inclusion_semiaxesDimensions(0, a1, a2);
  this->set_inclusion_EulerAnglesDEG(0, ea);
  
  this->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  
  this->convert_to_equivalent_problem();
  
}

void Problem :: set_approximation (int val)
{
  if (val < 0 || val>2)_errorr("Wrong approximation value.");
  if (this->approximation > 0)_errorr("Another approximation is already set. It can be set only once.");
  int n_points = 0;
  int n_approximations = 0;
  
  this->approximation = val;
  //2D
  // x y | xx yy xy
  //3D
  // x y z | xx yy zz yz xy xy
  
  if (val == 1) {
    n_points = this->give_dimension() * 2 + 1;
    n_approximations = this->give_dimension();
    for (int i = 0; i < noIncl; i++) {
      this->inclusions[i]->n_approx_points = n_points;
      this->inclusions[i]->approx_points = AllocateArray2D(n_points, 3);
      this->inclusions[i]->approx_points_eps_tau = AllocateArray2D(n_points, this->give_VM_TENS_RANGE());
      this->inclusions[i]->approx_coef = new double**[this->give_nLC()];
      for (int j = 0; j < this->give_nLC(); j++) {
    this->inclusions[i]->approx_coef[j] = new double*[n_approximations];
    for (int k = 0; k < n_approximations; k++) {
      this->inclusions[i]->approx_coef[j][k] = new double[this->give_VM_TENS_RANGE()];
      for (int ii = 0; ii < this->give_VM_TENS_RANGE(); ii++)this->inclusions[i]->approx_coef[j][k][ii] = 0.;
    }
      }
      CleanArray2d(this->inclusions[i]->approx_points, n_points, 3);
      CleanArray2d(this->inclusions[i]->approx_points_eps_tau, n_points, this->give_VM_TENS_RANGE());
      
      this->inclusions[i]->approx_points[1][0] -= approx_point_pos1*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[2][0] += approx_point_pos1*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[3][1] -= approx_point_pos1*this->inclusions[i]->a[1];
      this->inclusions[i]->approx_points[4][1] += approx_point_pos1*this->inclusions[i]->a[1];
      if (this->give_dimension() == 3) {
    this->inclusions[i]->approx_points[5][2] -= approx_point_pos1*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[6][2] += approx_point_pos1*this->inclusions[i]->a[2];
      }
      for (int j = 0; j < n_points; j++)this->inclusions[i]->transformCoords_L2G(this->inclusions[i]->approx_points[j], this->inclusions[i]->approx_points[j]);
    }
  }
  if (val == 2) {
    if (this->give_dimension() == 2) {
      n_points = 9;
      n_approximations = 5;
    }
    else {
      n_points = 19;
      n_approximations = 9;
    }
    for (int i = 0; i < noIncl; i++) {
      this->inclusions[i]->n_approx_points = n_points;
      this->inclusions[i]->approx_points = AllocateArray2D(n_points, 3);
      this->inclusions[i]->approx_points_eps_tau = AllocateArray2D(n_points, this->give_VM_TENS_RANGE());
      this->inclusions[i]->approx_coef = new double**[this->give_nLC()];
      for (int j = 0; j < this->give_nLC(); j++) {
    this->inclusions[i]->approx_coef[j] = new double*[n_approximations];
    for (int k = 0; k < n_approximations; k++) {
      this->inclusions[i]->approx_coef[j][k] = new double[this->give_VM_TENS_RANGE()];
      for (int ii = 0; ii < this->give_VM_TENS_RANGE(); ii++)this->inclusions[i]->approx_coef[j][k][ii] = 0.;
    }
      }
      CleanArray2d(this->inclusions[i]->approx_points, n_points, 3);
      CleanArray2d(this->inclusions[i]->approx_points_eps_tau, n_points, this->give_VM_TENS_RANGE());
      
      this->inclusions[i]->approx_points[1][0] -= approx_point_pos1*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[2][0] += approx_point_pos1*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[3][1] -= approx_point_pos1*this->inclusions[i]->a[1];
      this->inclusions[i]->approx_points[4][1] += approx_point_pos1*this->inclusions[i]->a[1];
      
      this->inclusions[i]->approx_points[5][0] -= 0.5*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[5][1] -= 0.5*this->inclusions[i]->a[1];
      this->inclusions[i]->approx_points[6][0] += 0.5*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[6][1] -= 0.5*this->inclusions[i]->a[1];
      this->inclusions[i]->approx_points[7][0] -= 0.5*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[7][1] += 0.5*this->inclusions[i]->a[1];
      this->inclusions[i]->approx_points[8][0] += 0.5*this->inclusions[i]->a[0];
      this->inclusions[i]->approx_points[8][1] += 0.5*this->inclusions[i]->a[1];
      
      if (this->give_dimension() == 3) {
    this->inclusions[i]->approx_points[9][2] -= approx_point_pos1*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[10][2] += approx_point_pos1*this->inclusions[i]->a[2];
    
    this->inclusions[i]->approx_points[11][0] -= 0.5*this->inclusions[i]->a[0];
    this->inclusions[i]->approx_points[11][2] -= 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[12][0] += 0.5*this->inclusions[i]->a[0];
    this->inclusions[i]->approx_points[12][2] -= 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[13][0] -= 0.5*this->inclusions[i]->a[0];
    this->inclusions[i]->approx_points[13][2] += 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[14][0] += 0.5*this->inclusions[i]->a[0];
    this->inclusions[i]->approx_points[14][2] += 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[15][1] -= 0.5*this->inclusions[i]->a[1];
    this->inclusions[i]->approx_points[15][2] -= 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[16][1] += 0.5*this->inclusions[i]->a[1];
    this->inclusions[i]->approx_points[16][2] -= 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[17][1] -= 0.5*this->inclusions[i]->a[1];
    this->inclusions[i]->approx_points[17][2] += 0.5*this->inclusions[i]->a[2];
    this->inclusions[i]->approx_points[18][1] += 0.5*this->inclusions[i]->a[1];
    this->inclusions[i]->approx_points[18][2] += 0.5*this->inclusions[i]->a[2];
    
      }
      
      
      for (int j = 0; j < n_points; j++)this->inclusions[i]->transformCoords_L2G(this->inclusions[i]->approx_points[j], this->inclusions[i]->approx_points[j]);
    }
  }
  
  /*for (int i = 0; i < noIncl; i++) {
    printf("\nIncl n. %d\n(%d points)\n", i, this->inclusions[i]->n_approx_points);
    for (int j = 0; j < this->inclusions[i]->n_approx_points; j++) {
    printf("%lf %lf %lf\n",this->inclusions[i]->approx_points[j][0], this->inclusions[i]->approx_points[j][1], this->inclusions[i]->approx_points[j][2]);
    }
    }
    enter();*/
}


// /**
//  * @brief scan_DATA_field_head
//  *
//  * @param stream
//  * @param strm
//  * @param str
//  * @param type
//  * @param nexc
//  * @param format
//  * @param name
//  */
// void scan_DATA_field_head (Stream *stream, Stream *strm, const char *str, char type, int nexc, char format, const char *name)
// {
//   char flint[8];
//   if      (format == 'i') { if (stream->isFile()) sprintf(flint,"int");    else sprintf(flint,"Int32");   }
//   else if (format == 'f') { if (stream->isFile()) sprintf(flint,"float");  else sprintf(flint,"Float32"); }
//   else errol;
//   
//   char auxs[255];
//   sprintf (auxs, "Error in section %s", name);
//   
//   if (stream->isFile()) {
//     SP_scan_expected_word_exit   (str, flint, auxs, CASE);
//     if (type == 'S') {
//       SP_scan_expected_number_exit (str, nexc, auxs);
//       FP_scan_expected_word_exit   (stream->file(), "LOOKUP_TABLE default", auxs, CASE);
//       FP_skip_line (stream->file());
//     }
//     strm->redefine( stream->file() );
//   }
//   else {
//     strm->redefine( XP_giveDAtext  (stream->tixel(), 0, false, "ascii", flint, name, &nexc));
//   }
// }
void scan_FIELD_head (FILE *stream, long noincl, int nexc, char format, const char *name)
{
  char flint[8];
  if      (format == 'i')  sprintf(flint,"int");
  else if (format == 'f')  sprintf(flint,"float");
  else errol;
  
  char auxs[255];
  sprintf (auxs, "Error in section %s", name);
  
  FP_scan_expected_number_exit (stream, noincl, auxs);
  FP_scan_expected_number_exit (stream, nexc, auxs);
  FP_scan_expected_word_exit   (stream, flint, auxs, CASE);
}

bool Problem :: file_type_s2e (const char *s)
{
  if      (_STRCMP(s, "2D"     ) == 0) { twodim = true;  return false; }
  else if (_STRCMP(s, "2Dequiv") == 0) { twodim = true;  return true;  }
  else if (_STRCMP(s, "3D"     ) == 0) { twodim = false; return false; }
  else if (_STRCMP(s, "3Dequiv") == 0) { twodim = false; return true;  }
  else _errorr2 ("Invalid structure of the given VTK file, the key word %s at the second line is not supported", s);
  return false;
}

void Problem :: read_input_file (const char *filename)
{
  if (data_set) _errorr("data_set already given");
  data_set = true;
  
  if (this->verbose > 0)
    printf( "Reading of VTK file of multiple inclusion, please wait ...\n" );
  
  // open stream with automatic detection VTK legaci/XML file format
  char *fn = new char[255];
  strcpy(fn, filename);
  Stream *stream = new Stream();
  stream->open(STRM_void, "r", (const char*&)fn); // do not delete fn, it is deleted inside of the function
  
  //
  char *WORD=NULL;
  const XMLElement *xnel = NULL;
  
  bool lgc = stream->isFile();
  bool amd = false;
  
  if (lgc) WORD = new char[255];
  
  
  // HEAD OD FILE
  if (lgc) {
    // 1. line
    FP_skip_line (stream->file());
    
    // 2. line
    FP_scan_word (stream->file(), WORD);
    FP_skip_line (stream->file());
    
    amd = this->file_type_s2e (WORD);
    
    // 3. line
    FP_scan_expected_word_exit (stream->file(), "ASCII", "Invalid structure of the given VTK file", true);
    FP_skip_line (stream->file());
    // 4. line
    FP_scan_expected_word_exit (stream->file(), "DATASET UNSTRUCTURED_GRID", "Invalid structure of the given VTK file", true);
    FP_skip_line (stream->file());
  }
  else {
    // check head of xml file and read possible ctrl section
    // return XMLElement "Piece"
    const XMLNode *parent = stream->tixnod();
    
    //* declaration - delete if present
    if (parent->FirstChild()->ToDeclaration() != NULL)  stream->tixnod()->DeleteChild( stream->tixnod()->FirstChild() );
    
    //* comment - delete if present
    while (true)
      if (parent->FirstChild()->ToComment() != NULL)  stream->tixnod()->DeleteChild( stream->tixnod()->FirstChild() );
      else                                            break;
    
    //
    parent = XP_give_unique_expected_elem (parent, "VTKFile");
    stream->relink_downF();
    
    XP_check_expected_attribute (parent, "type", "UnstructuredGrid");
    XP_check_expected_attribute (parent, "version", "0.1");
    XP_check_expected_attribute (parent, "byte_order", "LittleEndian");
    
    //
    parent = XP_give_unique_expected_elem (parent, "UnstructuredGrid", false);
    stream->relink_downF();
    
    // read AppendedData section
    parent = parent->NextSibling();
    if (! parent)
      _errorr("No AppendedData section");
    
    if (_STRCMP(parent->Value(), "AppendedData") != 0)     _errorr ("name of xelement is not \"AppendedData\"");
    
    if (! (parent = parent->FirstChild ()))  _errorr("No child in AppendedData");
    if (   parent->Value()[0] != '_'      )  _errorr("There is no character \'_\' at the start of section AppendedData");
    parent = parent->NextSibling();
    
    if (_STRCMP(parent->Value(), "Input_Type") != 0)  _errorr("name of xelement is not Input_Type");
    amd = this->file_type_s2e (parent->ToElement()->GetText());
    
    while (parent->NextSibling()) {
      parent = parent->NextSibling();
      const char *text = parent->Value();
      long auxl;
      
      if (_STRCMP(text, "SBA_mode") == 0) {
        sscanf (parent->ToElement()->GetText(), "%ld", &auxl);
        SBA_optimized = (bool)(auxl-1);
      }
      else if (_STRCMP(text, "conversion_mode") == 0) {
        sscanf (parent->ToElement()->GetText(), "%d", &intFieldsShape);
      }
      else if (_STRCMP(text, "Matrix_record") == 0) {
        double e, nu;
        sscanf (parent->ToElement()->GetText(), "%lf %lf", &e, &nu);
        matrix->set_E_nu(e, nu);
      }
      else if (_STRCMP(text, "Remote_strains") == 0) {
        long nlc;
        const char *str = parent->ToElement()->GetText();
        SP_scan_expected_number_exit (str, this->give_TENS_RANGE(), "wrong length of Remote_strains tensor, expected 4 or 9");
    SP_scan_number(str, nlc);
    this->matrix->set_nlc(nlc);
    // load cases
    if (this->matrix->scan_Remote_strains(str) != true)   _errorr2 ("There is few numbers in section %s", WORD);
      }
      else if (_STRCMP(text, "Min_axes_diff") == 0) { double m;  sscanf (parent->ToElement()->GetText(), "%lf", &m);  this->set_semiaxes_min_difference(m); }
      else if (_STRCMP(text, "Max_axes_diff") == 0) { double m;  sscanf (parent->ToElement()->GetText(), "%lf", &m);  this->set_semiaxes_max_difference(m); }
    }
    
    // go to Piece
    XP_give_unique_expected_elem (stream->tixnod(), "Piece");
    stream->relink_downF();
    xnel = stream->tixel();
  }
  
  
  // variable declaration
  Stream auxstrm;
  // 
  bool status = true, prnt;
  int auxi, les, ges;
  long i, j, nn, auxl, noel, cno, offset;
  char ed='-';
  const char *data=NULL, *offs=NULL, *typs, *str;
  char LINE[1023];
  int lLINE=1023;
  VTKrange range = VTKR_void;
  int nLC;
  
  cno = 0;
  prnt = true;
  les = ges = 0;
  while (true) {
    //* key WORD
    if (lgc) {
      if (fgets (LINE, lLINE, stream->file()) == NULL)  break;
      str = LINE;
      SP_scan_word (str, WORD);
    }
    else {
      if (prnt) { if (! stream->relink_downF())  errol; else prnt = false; }
      else      { if (! stream->relink_next())   break; }
      WORD = (char *)stream->tixel()->Value();
    }
    
    //*  ******************
    //*  ***   POINTS   ***
    //*  ******************
    if      (_STRCMP(WORD, "Points") == 0) {
      if (range == VTKR_void)  range = VTKR_points;  else errol;
      
      // number of nodes
      if (lgc) {
    SP_scan_number (str, noel);
    SP_scan_expected_word2_exit  (str, "double", "float", "Invalid structure of the given VTK file, unsupported datatype of POINTS", CASE);
      }
      else {
    xnel->ToElement()->QueryLongAttribute("NumberOfPoints", &noel);
    data = XP_giveDAtext (stream->tixel(), 1, true, "ascii", "Float32", NULL, &(auxi=3));
      }
      
      noIncl = noel; // number of inclusions
      if (twodim) { inclusions = new Inclusion*[noIncl];  for (int i=0; i<noIncl; i++) inclusions[i] = new InclusionRecord2D(i, this); }
      else        { inclusions = new Inclusion*[noIncl];  for (int i=0; i<noIncl; i++) inclusions[i] = new InclusionRecord3D(i, this); }
      
      // for (i=0; i<noIncl; i++)  inclRec[i].set_id_problem(i, this);
      
      // coords
      noIncl = 0;
      if (twodim) {
    while (noel) { noel--;
          if (lgc) { status += ( fscanf (stream->file(), "%lf %lf %lf", inclusions[noIncl]->origin, inclusions[noIncl]->origin+1, inclusions[noIncl]->origin+2) == 3 );  if (inclusions[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0"); }
          else     { status += ( sscanf (data,           "%lf %lf %lf", inclusions[noIncl]->origin, inclusions[noIncl]->origin+1, inclusions[noIncl]->origin+2) == 3 );  if (inclusions[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
      noIncl ++;
    }
      }
      else {
    while (noel) { noel--;
          if (lgc) { status += ( fscanf (stream->file(), "%lf %lf %lf", inclusions[noIncl]->origin, inclusions[noIncl]->origin+1, inclusions[noIncl]->origin+2) == 3 ); }
          else     { status += ( sscanf (data,           "%lf %lf %lf", inclusions[noIncl]->origin, inclusions[noIncl]->origin+1, inclusions[noIncl]->origin+2) == 3 );  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
      noIncl ++;
    }
      }
      
      if (lgc)  FP_skip_line (stream->file());
    }
    //*  *************************************
    //*  ***   UNSTRUCTURED GRID - CELLS   ***
    //*  *************************************
    else if (_STRCMP(WORD, "CELLS") == 0) {
      if (range == VTKR_points)  range = VTKR_cells;  else errol;
      
      // number of cells
      if (lgc) {
        SP_scan_number (str, noel);
        SP_scan_number (str, cno);
      }
      else {
        xnel->ToElement()->QueryLongAttribute("NumberOfCells" , &noel);
        data = XP_giveDAtext (stream->tixel(), 1, false, "ascii", "Int32", "connectivity", NULL);
        offs = XP_giveDAtext (stream->tixel(), 2, false, "ascii", "Int32", "offsets", NULL);
        typs = XP_giveDAtext (stream->tixel(), 3, true,  "ascii", "UInt8", "types", NULL);
      }
      
      if (!lgc) offset = 0;
      // elem nodes
      for (i=0; i<noel; i++) {
        // nn
        if (lgc) { fscanf (stream->file(),"%ld",&nn); cno = cno - nn - 1; }
        else     { nn = offset;  status += SP_scan_number (offs, offset);  nn = offset - nn; }
        //
    if (nn != 1) _errorr3 ("The number of CELL(%ld) nodes is %ld, it should be 1", i, nn);
        //
    if (lgc) fscanf (stream->file(),"%ld", &auxl);
    else     status += SP_scan_number (data, auxl);
    if (auxl != i) _errorr4 ("The CELL(%ld) node ID is %ld, it should be %ld", i, auxl, i);
      }
      
      if (lgc) FP_skip_line (stream->file(), 1);
      
      // types of elems
      if (lgc) {
    fgets (LINE, lLINE, stream->file());  str=LINE;
        SP_scan_expected_word_exit   (str, "CELL_TYPES", "Invalid structure of the given VTK file", CASE);
        SP_scan_expected_number_exit (str,         noel, "Invalid structure of the given VTK file, the number following CELL_TYPES");
      }
      
      while (noel) { noel--;
        if (lgc)  fscanf (stream->file(),"%ld",&auxl);
        else      status += SP_scan_number (typs, auxl);
        if (auxl != 1) _errorr2 ("Unsupported CELL_TYPE %ld, only \"1\" is supported", auxl);
      }
      
      if (lgc) FP_skip_line (stream->file(), 1);
    }
    //*  *******************************
    //*  ***   POLYDATA - ELEMENTS   ***
    //*  *******************************
    else if (_STRCMP(WORD, "LINES") == 0  ||  _STRCMP(WORD, "POLYGONS" ) == 0  ||  _STRCMP(WORD, "POLYS") == 0  ||  _STRCMP(WORD, "Vertices") == 0  ||  _STRCMP(WORD, "Verts") == 0  ||  _STRCMP(WORD, "Strips") == 0) {
      _errorr2("Invalid structure of the given VTK file, the key word %s is not supported", WORD);
    }
    //*  **********************
    //*  ***   POINT_DATA   ***
    //*  **********************
    else if (_STRCMP(WORD, "POINT_DATA") == 0  ||  _STRCMP(WORD, "PointData") == 0) {
      if (range == VTKR_cells  || range == VTKR_points  ||  range == VTKR_fields)  range = VTKR_pd;  else errol;
      ed = 'p';
      // head
      if (lgc) { SP_scan_expected_number_exit (str, noIncl, "The number following POINT_DATA is not equal to number of inclusions"); }
      else     { prnt = true; }
    }
    //*  *********************
    //*  ***   CELL_DATA   ***
    //*  *********************
    else if (_STRCMP(WORD, "CELL_DATA") == 0  ||  _STRCMP(WORD, "CellData") == 0) {
      if (range == VTKR_cells || range == VTKR_data)  range = VTKR_cd;  else errol;
      ed = 'c';
       // head
      if (lgc) { SP_scan_expected_number_exit (str, noIncl, "The number following CELL_DATA is not equal to number of inclusions"); }
      else     { prnt = true; }
    }
    //*  ******************
    //*  ***   FIELDS   ***
    //*  ******************
    else if (_STRCMP(WORD, "FIELD") == 0) {
      if (range == VTKR_cells || range == VTKR_data)  range = VTKR_fields;  else errol;
      
      if (!lgc) _errorr("FIELD in nonlegacy");
      
      // head
      if (lgc)   SP_scan_word (str, WORD);
      else     { WORD = (char *)stream->tixel()->ToElement()->Attribute("Name");  if (!WORD)  _errorr ("there is no attribute \"Name\""); }
      
      if (_STRCMP(WORD, "unstructured_data") != 0)  _errorr("...");
      
      SP_scan_number (str, noel);
      
      for (i=0; i<noel; i++) {
    FP_scan_word (stream->file(), WORD);
    
    if (_STRCMP(WORD, "Acting_inclusions_list") == 0) {
          auxl = this->totalNoActingIncls;
      scan_FIELD_head (stream->file(), 1, auxl, 'i', WORD);
      if (auxl)
        for (j=0; j<noIncl; j++)
              if (FP_scan_array  (stream->file(), inclusions[i]->noActingIncls, inclusions[i]->actingIncls) != true)  _errorr2 ("There is few numbers in section %s", WORD);
    }
    else if (_STRCMP(WORD, "Remote_strains") == 0) {
      // length of tensor
      FP_scan_expected_number_exit (stream->file(), this->give_TENS_RANGE(), "wrong length of Remote_strains tensor, expected 4 or 9");
      // number of load cases
      fscanf (stream->file(), "%d", &nLC);
          this->matrix->set_nlc(nLC);
      FP_scan_expected_word_exit (stream->file(), "float", "float expected in remote strain", CASE);
      // load cases
      if (this->matrix->scan_Remote_strains(stream->file()) != true)   _errorr2 ("There is few numbers in section %s", WORD);
    }
        else if (_STRCMP(WORD, "SBA_mode") == 0) {
          scan_FIELD_head (stream->file(), 1, 1, 'i', WORD);
          fscanf (stream->file(), "%ld", &auxl);
          SBA_optimized = (bool)(auxl-1);
        }
        else if (_STRCMP(WORD, "conversion_mode") == 0) {
          scan_FIELD_head (stream->file(), 1, 1, 'i', WORD);
          fscanf (stream->file(), "%d", &intFieldsShape);
        }
        else if (_STRCMP(WORD, "Matrix_record") == 0) {
      scan_FIELD_head (stream->file(), 1, 2, 'f', WORD);
      double e, nu;
      fscanf (stream->file(), "%lf %lf", &e, &nu);
      
      this->matrix->set_E_nu(e, nu);
    }
    else
      _errorr3 ("Invalid name of %d -th FIELD \"%s\"", i, WORD);
      }
    }
    //*  **************************
    //*  ***   DATA - SCALARS   ***
    //*  **************************
    else if (_STRCMP(WORD, "SCALARS") == 0  ||  _STRCMP(WORD, "VECTORS") == 0  ||  _STRCMP(WORD, "TENSORS") == 0  ||  _STRCMP(WORD, "DataArray") == 0) {
      if (range == VTKR_pd || range == VTKR_cd || range == VTKR_data)  range = VTKR_data;  else errol;
      
      char type = WORD[0]; // save type of data - Scalar, Vector, Tensor, DataArray
      // head
      if (lgc)   SP_scan_word (str, WORD);
      else     { WORD = (char *)stream->tixel()->ToElement()->Attribute("Name");  if (!WORD)  _errorr ("there is no attribute \"Name\""); }
      
      //*  *** DATA for POINTS ***
      if (ed == 'p') {
    // bacic
        if      (_STRCMP(WORD, "Inclusion_shape"                            ) == 0) { auxl = 1;                                   scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'i', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_number (&auxstrm,       auxl                          )     != true)  _errorr2 ("There is few numbers in section %s", WORD); else inclusions[j]->shape = (InclusionGeometry)auxl; }
        else if (_STRCMP(WORD, "Youngs_modulus"                             ) == 0) { auxl = 1;                                   scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_number (&auxstrm,       inclusions[j] ->E                )     != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Poissons_ratio"                             ) == 0) { auxl = 1;                                   scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_number (&auxstrm,       inclusions[j] ->nu               )     != true)  _errorr2 ("There is few numbers in section %s", WORD); }
    
        else if (_STRCMP(WORD, "Semiaxes_dimensions"                        ) == 0) { auxl = this->give_VECT_RANGE();         scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, inclusions[j] ->a                )     != true)  _errorr2 ("There is few numbers in section %s", WORD); }
    
        else if (_STRCMP(WORD, "Euller_angles_deg"                          ) == 0) { auxl = this->give_EA_RANGE();               scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (inclusions[j]->scan_eAngles_DEG(&auxstrm)                              != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Euller_angles_rad"                          ) == 0) { auxl = this->give_EA_RANGE();               scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (inclusions[j]->scan_eAngles_RAD(&auxstrm)                              != true)  _errorr2 ("There is few numbers in section %s", WORD); }
    
    // else if (_STRNCMP(WORD, "Remote_strains_", 15                       ) == 0) {
    //   auxl = atoi(WORD+18);
    //   if (this->give_nLC() == 0) { this->set_nlc_allocate_nlc_fields(auxl);  lc = les = ges = 0; }
    //   else                if (auxl != nLC)  _errorr("Number of load cases mismatch");
    //   WORD[17] = '\0';
    //   auxl = atoi(WORD+15) - 1;
    //   pokud tady bude cteni, tak se musi udelat pres fci 'scan_...' uvnitr ktere se udela zmena notace z tensorove na feep
    //   if (auxl != lc) _errorr("Load case ID mismatch");                           auxl = this->give_TENS_RANGE();             scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, inclRec[j]->Remote_strains[lc])  != true)  _errorr2 ("There is few numbers in section %s", WORD); lc++; }
    
    // computed values
        else if (_STRCMP(WORD, "Action_radius"                              ) == 0) { auxl = 1;                                   scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_number (&auxstrm,       inclusions[j]->actionRadius     )       != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Acting_inclusions_number"                   ) == 0) { auxl = 1;                                   scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'i', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_number (&auxstrm,       inclusions[j]->noActingIncls    )       != true)  _errorr2 ("There is few numbers in section %s", WORD); }
    
    else if (_STRNCMP(WORD, "Local_Eq_strain_", 16                      ) == 0) {
      auxl = atoi(WORD+19);      if (auxl != nLC) _errorr("Local_Eq_strain: Number of load cases mismatch");  WORD[18] = '\0';
      auxl = atoi(WORD+16) - 1;  if (auxl != les) _errorr("Local_Eq_strain ID mismatch");
                                                                                      auxl = this->give_VM_TENS_RANGE();          scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (                     inclusions[j]->scan_locEigStrain_LC(&auxstrm, les) != true)  _errorr2 ("There is few numbers in section %s", WORD); les++; }
    else if (_STRNCMP(WORD, "Global_Eq_strain_", 17                     ) == 0) {
      auxl = atoi(WORD+20);      if (auxl != nLC) _errorr("Global_Eq_strain: Number of load cases mismatch");  WORD[19] = '\0';
      auxl = atoi(WORD+17) - 1;  if (auxl != ges) _errorr("Global_Eq_strain ID mismatch");
                                                                                      auxl = this->give_VM_TENS_RANGE();          scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (((InclusionRecord3D*)inclusions[j])->scan_globEigStrain_LC(&auxstrm, ges) != true)  _errorr2 ("There is few numbers in section %s", WORD); ges++; }
    
        else if (_STRCMP(WORD, "Stiffnesses_C"                              ) == 0) { auxl = this->give_ISO_C_RANGE();            scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, ((InclusionRecord3D*)inclusions[j])->C    )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Elliptic_potentials"                        ) == 0) { auxl = this->give_EL_POT_RANGE();           scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, ((InclusionRecord3D*)inclusions[j])->eInt )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Eshelby_tensor"                             ) == 0) { auxl = this->give_ISO_C_RANGE();            scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, ((InclusionRecord3D*)inclusions[j])->S    )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Eshelby_tensor_inverse"                     ) == 0) { auxl = this->give_ISO_C_RANGE();            scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, ((InclusionRecord3D*)inclusions[j])->SInv )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Global_2_local_transform_matrix_for_vectors") == 0) { auxl = this->give_TRNSFM_MTRX_VEC_RANGE();  scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, (                    inclusions[j])->T    )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
      //else if (_STRCMP(WORD, "Local_2_global_transform_matrix_for_vectors") == 0) { auxl = this->give_TRNSFM_MTRX_VEC_RANGE();  scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, (                    inclRec[j])->TInv )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Global_2_local_transform_matrix_for_tensors") == 0) { auxl = this->give_TRNSFM_MTRX_TENS_RANGE(); scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, ((InclusionRecord3D*)inclusions[j])->Te   )                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
        else if (_STRCMP(WORD, "Local_2_global_transform_matrix_for_tensors") == 0) { auxl = this->give_TRNSFM_MTRX_TENS_RANGE(); scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);  for (j=0; j<noIncl; j++) if (ST_scan_array  (&auxstrm, auxl, ((InclusionRecord3D*)inclusions[j])->TeInv)                   != true)  _errorr2 ("There is few numbers in section %s", WORD); }
    
        else
      _errorr2 ("Invalid name of SCALARS/DataArray for POINT_DATA/PointData \"%s\"", WORD);
    
    //
    if (_STRCMP(WORD, "Acting_inclusions_number") == 0)
      for (j=0; j<noIncl; j++) {
            inclusions[j]->ActingIncls_allocate();
            this->totalNoActingIncls += inclusions[j]->noActingIncls;
      }
      }
      
      //*  *** DATA for CELLS ***
      if (lgc) FP_skip_line (stream->file(), 1);
      else { if (stream->tixel()->NextSibling() == NULL)  stream->relink_up(); }
    }
    else if (_STRCMP(WORD, "") == 0) {;}
    else _errorr2 ("Invalid structure of the given VTK file, unexpected key word \"%s\"", WORD);
    
    if (lgc) { if (cno) _errorr2 ("cno(%ld) is not zero", cno); }
    else {
      if (data && data[0] != '\0')  _errorr2("CellData of name \"%s\" has too many numbers", WORD);
      if (offs && offs[0] != '\0')  _errorr("error");
    }
  }
  
  if (lgc) delete [] WORD;
  if (status == false)  _errorr("error");
  
  //if (lc != nLC) _errorr("few Remote_strains in input file");
  if (this->data_equivalent) {
    if (les != nLC) _errorr("few Local_Eq_strain in input file");
    if (ges != nLC) _errorr("few Global_Eq_strain in input file");
  }
  
  //* close
  stream->close();
  delete stream;
  
  if (amd) data_equivalent = true;
  
  // initialization
  if (this->data_equivalent) {
    if (matrix == NULL)          _errorr("Matrix record is not given");
    //if (SelfBalAlgorithm == SBAT_VOID) _errorr("SelfBalAlgorithm is not given");
  }
  
} // end of function Problem :: read_vtk_file (const char *vtkTopologyFile)


void Problem :: convert_to_equivalent_problem (void)
{
  if (!this->data_set)     _errorr("set");
  if ( this->data_equivalent) _errorr("initialized");
  
  if (matrix == NULL)          _errorr("Matrix record is not given");
  //if (SelfBalAlgorithm == SBAT_VOID) _errorr("SelfBalAlgorithm is not given");
  
  
  if (this->verbose > 0)
    printf( "Creating inclusion records, please wait ...\n" );
  
  Inclusion *prevInclRec = NULL;
  
  for (int i=0; i<noIncl; i++) {
    inclusions[i]->initialize(prevInclRec);
    prevInclRec = inclusions[i];
  }
  
  // inclusion dependencies, e.g. searching for adjacent inclusion with respect to acting radius
  if (SBA_optimized)  this->findAffectedInclusions();
  
  // incorporation of mutual interactions
  this->updateEigenstrainsBySBalgorithm();
  
  //
  this->data_equivalent = true;
}

void Problem :: updateEigenstrainsBySBalgorithm (void)
{
  clock_t selfUpdateTimeStart = clock( ); // time measurement 'start'
  
  long i, j;
  for (i=0; i<this->give_nLC(); i++) {
    if (this->verbose > 0)
      printf( "Transformation strain update of LoadCase no. %ld in progress, please wait ... \n", i );
    
    //
    for (j=0; j<noIncl; j++)
      inclusions[j]->SBA_computeInitialStrains(i);
    
    if (this->twodim) {
      this->updateEigenstrainsBySBalgorithm_2D (i);
    }
    else {
      selfBalanceAlgorithm::totalEigStrainInInclCentroidsUpdate (this);
      //
      for (j=0; j<noIncl; j++)
        ((InclusionRecord3D*)inclusions[j])->init_EigStrain_LC(i);
    }
  }
  
  //
  clock_t selfUpdateTimeEnd = clock( ); // time measurement 'end'
  // printing of computional time status
  if (this->verbose > 0)
    printf( "updateEigenstrainsBySBalgorithm: multiple self-balance algorithm time: %.5e sec\n",
        ( selfUpdateTimeEnd - selfUpdateTimeStart )/ ( double ) CLOCKS_PER_SEC );
  
}

void Problem :: updateEigenstrainsBySBalgorithm_2D (int strainID)
{
  if (SBA_type == SBAT_STANDA) {
    //#####################################################################################
    //##################################################### STANDOVO ######################
    //#####################################################################################
    double **last_eps_tau = NULL;
    
    double norma;
    
    int step = 0;
    
    do {
      // maximum number of interations reached - non convergence
      if (step >= this->get_SBA_maxiters())  _errorr("convergence not reached");
      // required number of interations reached - break
      if (step >= this->get_SBA_reqiters())  break;
      
      step++;
      norma = updateEpsTauInSBal(strainID, last_eps_tau);
      
    } while (norma > _SELF_BALANCE_NORM_LIMIT_);
    
    if (last_eps_tau != NULL) {
      for (int i = 0; i < this->give_nLC(); i++)delete[]last_eps_tau[i];
      delete[]last_eps_tau;
    }
    
    if (this->verbose > 0)
      printf("\nBallancing : %d steps\n", step);
    
    FILE *f;
    f = fopen("./EpsilonyTau.txt", "wt");
    for (int i = 0; i < this->noIncl; i++) {
      fprintf(f,"Inkluze %d\n", i);
      for (int j = 0; j < this->inclusions[i]->n_approx_points; j++) {
    for (int k = 0; k < this->give_VM_TENS_RANGE(); k++)fprintf(f,"%lf ", this->inclusions[i]->approx_points_eps_tau[j][k]);
    fprintf(f,"\n");
      }
    }
    fprintf(f, "Lokalni souradnice\n");
    for (int j = 0; j < this->inclusions[0]->n_approx_points; j++) {
      fprintf(f, "%d] %lf %lf\n", j, this->inclusions[0]->approx_points[j][0] - this->inclusions[0]->origin[0], this->inclusions[0]->approx_points[j][1] - this->inclusions[0]->origin[1]);
    }
    
    fclose(f);
    
    //for (int i = 0; i < this->noIncl; i++)printf("%lf %lf %lf\n", this->inclusions[i]->Epsilon_Tau[strainID][0], this->inclusions[i]->Epsilon_Tau[strainID][1], this->inclusions[i]->Epsilon_Tau[strainID][2]);
    //e(111);
  }
  else if (SBA_type == SBAT_ORIGINAL) {
    //#####################################################################################
    //##################################################### ORIGINAL ######################
    //#####################################################################################
    // odkomentoval jsem to, protoze to mam v testovacich prikladkach
    
    int i, j, jj;
      vektor pricinek_strain(3);
      double norma_delta_e_s;
      
      int step = 0;
      
      vektor *suma_e_s = new vektor[noIncl];
      for (i = 0; i < noIncl; i++)
      suma_e_s[i].define(3);
      
      do {
      // maximum number of interations reached - non convergence
      if (step >= this->get_SBA_maxiters())  _errorr("convergence not reached");
      // required number of interations reached - break
      if (step >= this->get_SBA_reqiters())  break;
      
      step++;
      for (i = 0; i < noIncl; i++)
      suma_e_s[i].vynulovat();
      
      if (SBA_optimized) {
      for (i = 0; i < noIncl; i++) {
      for (j = 0; j < inclusions[i]->noActingIncls; j++) {
      jj = inclusions[i]->actingIncls[j];
      ((InclusionRecord2D*)inclusions[jj])->compute_from_eps_tau(inclusions[i]->origin, pricinek_strain, ((InclusionRecord2D*)inclusions[jj])->delta_e_t);
      suma_e_s[i].pricist(pricinek_strain);
      }
      }
      }
      else {
      for (i = 0; i < noIncl; i++) {
      for (j = 0; j < noIncl; j++) {
      if (i != j) {
      ((InclusionRecord2D*)inclusions[j])->compute_from_eps_tau(inclusions[i]->origin, pricinek_strain, ((InclusionRecord2D*)inclusions[j])->delta_e_t);
      suma_e_s[i].pricist(pricinek_strain);
      }
      }
      }
      }
      
      for (i = 0; i < noIncl; i++)
      ((InclusionRecord2D*)inclusions[i])->delta_e_s.rovna_se(suma_e_s[i]);
      
      for (i = 0; i < noIncl; i++) {
      ((InclusionRecord2D*)inclusions[i])->compute_eps_tau_back(((InclusionRecord2D*)inclusions[i])->delta_e_s, ((InclusionRecord2D*)inclusions[i])->delta_e_t);

        AddVector(((InclusionRecord2D*)inclusions[i])->delta_e_t.v , inclusions[i]->Epsilon_Tau[strainID] , this->give_VM_TENS_RANGE());
      }
      
      norma_delta_e_s = 0.0;
      for (i = 0; i < noIncl; i++)
      norma_delta_e_s += ((InclusionRecord2D*)inclusions[i])->delta_e_s.norma();
      
      } while (norma_delta_e_s > _SELF_BALANCE_NORM_LIMIT_);
      
      delete[] suma_e_s;
      
      if (this->verbose > 0)
      printf("\nBallancing : %d steps\n", step);
  }
  else if (SBA_type == SBAT_MESHLESS) {
    //#####################################################################################
    //##################################################### MESHLESS ######################
    //#####################################################################################
    printf("Ballancing loadcase %d\n", strainID);
    double **tuhost = new double*[3];
    for (int i = 0; i < 3; i++)tuhost[i] = new double[3];
    for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)tuhost[i][j] = 0.;
    tuhost[0][0] = this->matrix->C[0];
    tuhost[0][1] = this->matrix->C[1];
    tuhost[1][0] = this->matrix->C[2];
    tuhost[1][1] = this->matrix->C[3];
    tuhost[2][2] = this->matrix->C[4];
    //for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)tuhost[i][j] *= 1.2;
    double AvSA[4], AvSB[4];
    
    // SET OF THE WAY OF COMPUTING
    int way = 0;
    bool onlyMatrixC = false;
    bool onlyVxV = false;
    int plus3 = 0;
    bool onlyExternal = true;
                                
    int celk_cyklu = this->noIncl * 3, cyklu = 0;
    int celk_cyklu_integrace = 0;
    for (int iA = 0; iA < noIncl; iA++)for (int iAl = 0; iAl < 3; iAl++)for (int iB = iA; iB < noIncl; iB++)for (int iBl = 0; iBl < 3; iBl++)celk_cyklu_integrace++;
    
    
    double **ballance_matrix;
    double *right_side;
    double **load;
    load = new double*[3];
    for (int i = 0; i < 3; i++) load[i] = new double[3];
    load[0][0] = 1.; load[0][1] = 0.; load[0][2] = 0.;
    load[1][0] = 0.; load[1][1] = 1.; load[1][2] = 0.;
    load[2][0] = 0.; load[2][1] = 0.; load[2][2] = 1.;
    int q = this->noIncl * 3 + plus3;
    int nnodes[2];
    FILE *ft;
    
    int nnn, pointInInclusion;
    double res, resv[3];
    double matC[3][3];
    for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)matC[i][j] = 0.;
    
    double ***node_values;
    double **nodes;
    double *strain;
    strain = new double[3];
    right_side = new double[q];
    ballance_matrix = new double *[q];
    for (int i = 0; i < q; i++)ballance_matrix[i] = new double[q];
    
    double problem_size[3];
    //update the problem size to be suitable for the mesh
    for (int i = 0; i < 2; i++) {
      if (floor(this->problem_size_A[i] / this->node_dist)*this->node_dist != this->problem_size_A[i])this->problem_size_A[i] = floor(this->problem_size_A[i] / this->node_dist)*this->node_dist;
      if (floor(this->problem_size_B[i] / this->node_dist)*this->node_dist != this->problem_size_B[i])this->problem_size_B[i] = floor(this->problem_size_B[i] / this->node_dist)*this->node_dist;
    }
    for (int i = 0; i < 3; i++)problem_size[i] = problem_size_B[i] - problem_size_A[i];
    if (problem_size[0] * problem_size[1] <= 0 || problem_size[0] < 0 || problem_size[1] < 0)_errorr("Wrong problem size.");
    for (int i = 0; i < 2; i++) {
      nnodes[i] = floor(problem_size[i] / this->node_dist) + 1;
    }
    int nnnodes = nnodes[0] * nnodes[1];
    node_values = new double**[this->noIncl * 3];
    for (int i = 0; i < this->noIncl * 3; i++) {
      node_values[i] = new double*[nnnodes];
      for (int j = 0; j < nnnodes; j++)node_values[i][j] = new double[3];
    }
    nodes = new double*[nnnodes];
    for (int i = 0; i < nnnodes; i++)nodes[i] = new double[3];
    nnn = 0;
    for (int i = 0; i < nnodes[1]; i++)for (int j = 0; j < nnodes[0]; j++) {
    nodes[nnn][0] = problem_size_A[0] + this->node_dist*j;
    nodes[nnn][1] = problem_size_A[1] + this->node_dist*i;
    nodes[nnn][2] = 0.;
    nnn++;
      }
    
    
    double dsds = this->node_dist*this->node_dist;
    double volume = problem_size[0] * problem_size[1];
    double dsds_div_volume = dsds / volume;
    
    int iA, iB;
    printf("Pocitam hodnoty strainu...\n");
    // compute all the strain values
    for (int inclA = 0; inclA < this->noIncl; inclA++)
      for (int inclAload = 0; inclAload < 3; inclAload++) {
    iA = inclA * 3 + inclAload;
    for (int nodeID = 0; nodeID < nnnodes; nodeID++) {
      if (((InclusionRecord2D*)inclusions[inclA])->point_is_inside(nodes[nodeID], 0.) == true)
        ((InclusionRecord2D*)inclusions[inclA])->compute_strain_pert_from_eps_zero_int(nodes[nodeID], strain, load[inclAload]);
      else
        ((InclusionRecord2D*)inclusions[inclA])->compute_strain_pert_from_eps_zero_ext(nodes[nodeID], strain, load[inclAload]);
      
      /*if (((InclusionRecord2D*)inclusions[inclA])->point_is_inside(nodes[nodeID], 0.) == true)
        for (int i = 0; i < 3; i++)node_values[iA][nodeID][i] = 0.;
        else*/
      for (int i = 0; i < 3; i++)node_values[iA][nodeID][i] = strain[i];
    }
    cyklu++;
                                                                printf("%3.1lf %%\n", 100.*cyklu / celk_cyklu);
      }
    
    //e(951);
    
    // export do souboru..
    /*printf("Zapisuju do souboru...\n");
      FILE *f;
      FILE *tf;
      tf = fopen("./bal/inside.txt", "wt");
      char nazev_souboru[100];
      f = fopen("./bal/coords.txt", "wt");
      for (int nodeID = 0; nodeID < nnnodes; nodeID++)
      fprintf(f, "%lf %lf\n", nodes[nodeID][0], nodes[nodeID][1]);
      fclose(f);
      
      for (int inclA = 0; inclA < this->noIncl; inclA++)
      for (int inclAload = 0; inclAload < 3; inclAload++) {
      sprintf(nazev_souboru, "./bal/incl_%03d_load_%d.txt",inclA,inclAload);
      f = fopen(nazev_souboru, "wt");
      iA = inclA * 3 + inclAload;
      for (int nodeID = 0; nodeID < nnnodes; nodeID++) {
      fprintf(f, "%lf %lf %lf\n", node_values[iA][nodeID][0], node_values[iA][nodeID][1], node_values[iA][nodeID][2]);
      if(iA==0)fprintf(tf,"%d\n", this->give_inclusion_with_point_inside(nodes[nodeID]));
      }
      fclose(f);
      }
      fclose(tf);
      e(69);*/
                                
                                


    printf("Integruju...\n");
    cyklu = 0;
    // and now integrate
    if (onlyVxV) {
      for (int inclA = 0; inclA < this->noIncl; inclA++) {
    for (int inclAload = 0; inclAload < 3; inclAload++) {
      iA = inclA * 3 + inclAload;

      // compute the integrals for vector
      res = 0;
      for (int nodeID = 0; nodeID < nnnodes; nodeID++) {
        for (int i = 0; i < 3; i++)res -= node_values[iA][nodeID][i] * this->matrix->Remote_strain[strainID][i];
      }
      res *= dsds;
      right_side[inclA * 3 + inclAload] = res;

      for (int inclB = inclA; inclB < this->noIncl; inclB++) {
        for (int inclBload = 0; inclBload < 3; inclBload++) {
          iB = inclB * 3 + inclBload;

          // compute the integrals for matrix
          res = 0;
          for (int nodeID = 0; nodeID < nnnodes; nodeID++) {
        for (int i = 0; i < 3; i++)res += node_values[iA][nodeID][i] * node_values[iB][nodeID][i];
          }
          res *= dsds;
          ballance_matrix[inclA * 3 + inclAload][inclB * 3 + inclBload] = res;
          ballance_matrix[inclB * 3 + inclBload][inclA * 3 + inclAload] = res;

          cyklu++;

        }
        printf("%3.1lf %%\n", 100.*cyklu / celk_cyklu_integrace);
      }
    }
      }
    }
    else {
      for (int inclA = 0; inclA < this->noIncl; inclA++) {
    for (int inclAload = 0; inclAload < 3; inclAload++) {
      iA = inclA * 3 + inclAload;

      //MD
      for (int i = 0; i < 3; i++)AvSA[i] = 0.;

      // compute the integrals for vector
      res = 0;
      for (int nodeID = 0; nodeID < nnnodes; nodeID++) {
        //MD
        for (int i = 0; i < 3; i++)AvSA[i] += node_values[iA][nodeID][i];

        pointInInclusion = this->give_inclusion_with_point_inside(nodes[nodeID]);
        //if (pointInInclusion != inclA && pointInInclusion != inclA)pointInInclusion = -1;
        if (!onlyExternal || (onlyExternal && pointInInclusion == -1)) {
          if (pointInInclusion == -1 || onlyMatrixC) {
        matC[0][0] = this->matrix->C[0];
        matC[0][1] = this->matrix->C[1];
        matC[1][0] = this->matrix->C[2];
        matC[1][1] = this->matrix->C[3];
        matC[2][2] = this->matrix->C[4];
          }
          else {
        // toto jsme sem davali spolu, tak jsem to zakomentoval a pro jistotu zde nechal
        // matC[0][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[0] - this->matrix->C[0];
        // matC[0][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[1] - this->matrix->C[1];
        // matC[1][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[2] - this->matrix->C[2];
        // matC[1][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[3] - this->matrix->C[3];
        // matC[2][2] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[4] - this->matrix->C[4];
        
        matC[0][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[0];
        matC[0][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[1];
        matC[1][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[2];
        matC[1][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[3];
        matC[2][2] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[4];
          }
          for (int i = 0; i < 3; i++)resv[i] = 0.;
          for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)
                       resv[i] += matC[j][i] * node_values[iA][nodeID][j];
          for (int i = 0; i < 3; i++)res -= resv[i] * this->matrix->Remote_strain[strainID][i];

          //MD
          if (way == 2) {
        for (int i = 0; i < 3; i++)resv[i] = 0.;
        for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)
                         resv[i] += node_values[iA][nodeID][j] * tuhost[j][i];
        for (int i = 0; i < 3; i++)res += resv[i] * this->matrix->Remote_strain[strainID][i];//+
          }
        }
      }
      //e(111);
      //MD
      for (int i = 0; i < 3; i++)AvSA[i] *= dsds_div_volume;

      res *= dsds;

      //MD
      if (way == 1) {
        for (int i = 0; i < 3; i++)resv[i] = 0.;
        for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)
                     resv[i] += tuhost[j][i] * AvSA[j];
        for (int i = 0; i < 3; i++)res += resv[i] * this->matrix->Remote_strain[strainID][i] * volume;
      }

      right_side[inclA * 3 + inclAload] = res;

      for (int inclB = inclA; inclB < this->noIncl; inclB++) {
        for (int inclBload = 0; inclBload < 3; inclBload++) {
          iB = inclB * 3 + inclBload;

          //MD
          for (int i = 0; i < 3; i++)AvSB[i] = 0.;

          // compute the integrals for matrix
          res = 0;
          for (int nodeID = 0; nodeID < nnnodes; nodeID++) {
        //MD
        for (int i = 0; i < 3; i++)AvSB[i] += node_values[iB][nodeID][i];

        pointInInclusion = this->give_inclusion_with_point_inside(nodes[nodeID]);
        //if (pointInInclusion != inclA && pointInInclusion != inclA)pointInInclusion = -1;
        if (!onlyExternal || (onlyExternal && pointInInclusion == -1)) {
          if (pointInInclusion == -1 || onlyMatrixC) {
            matC[0][0] = this->matrix->C[0];
            matC[0][1] = this->matrix->C[1];
            matC[1][0] = this->matrix->C[2];
            matC[1][1] = this->matrix->C[3];
            matC[2][2] = this->matrix->C[4];
          }
          else {
            // toto jsme sem davali spolu, tak jsem to zakomentoval a pro jistotu zde nechal
            // matC[0][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[0] - this->matrix->C[0];
            // matC[0][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[1] - this->matrix->C[1];
            // matC[1][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[2] - this->matrix->C[2];
            // matC[1][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[3] - this->matrix->C[3];
            // matC[2][2] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[4] - this->matrix->C[4];
            
            matC[0][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[0];
            matC[0][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[1];
            matC[1][0] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[2];
            matC[1][1] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[3];
            matC[2][2] = ((InclusionRecord2D*)inclusions[pointInInclusion])->C[4];
          }
          for (int i = 0; i < 3; i++)resv[i] = 0.;
          for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)
                           resv[i] += node_values[iA][nodeID][j] * matC[j][i];
          for (int i = 0; i < 3; i++)res += resv[i] * node_values[iB][nodeID][i];

          //MD
          if (way == 2) {
            for (int i = 0; i < 3; i++)resv[i] = 0.;
            for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)
                         resv[i] += node_values[iA][nodeID][j] * tuhost[j][i];
            for (int i = 0; i < 3; i++)res -= resv[i] * node_values[iB][nodeID][i];//-
          }
        }

          }
          //MD
          for (int i = 0; i < 3; i++)AvSB[i] *= dsds_div_volume;

          res *= dsds;

          if (way == 1) {
        //MD
        for (int i = 0; i < 3; i++)resv[i] = 0.;
        for (int i = 0; i < 3; i++)for (int j = 0; j < 3; j++)
                         resv[i] += tuhost[j][i] * AvSA[j];
        for (int i = 0; i < 3; i++)res -= resv[i] * AvSB[i] * problem_size[0] * problem_size[1];
          }

          ballance_matrix[inclA * 3 + inclAload][inclB * 3 + inclBload] = res;
          ballance_matrix[inclB * 3 + inclBload][inclA * 3 + inclAload] = res;

          cyklu++;

        }
        printf("%3.1lf %%\n", 100.*cyklu / celk_cyklu_integrace);
      }
    }
      }
    }
    if (plus3) {
      for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) {
      ;
    }
      }
    }


    //printf("\nmatice sestavena\n");
    matice M(noIncl * 3, noIncl * 3);
    vektor V(noIncl * 3);

    for (int i = 0; i < noIncl * 3; i++)for (int j = 0; j < noIncl * 3; j++)M.g(i, j) = ballance_matrix[i][j];
    for (int i = 0; i < noIncl * 3; i++)V.v[i] = right_side[i];

    /*FILE *fo;
      fo = fopen("./matlab_mat.txt", "wt");
      M.print(fo);
      fclose(fo);
      fo = fopen("./matlab_vec.txt", "wt");
      V.print(fo);
      fclose(fo);*/

    printf("Resim soustavu rovnic...\n");
    gaussovaEliminace(M, V);
    V.print(stdout);
    e(555);
    double computed_load[3];
    for (int inclID = 0; inclID < noIncl; inclID++) {
      for (int j = 0; j < 3; j++)
    computed_load[j] = V.v[inclID * 3 + j];
      ((InclusionRecord2D*)inclusions[inclID])->compute_eps_tau(strainID, this->inclusions[inclID]->Epsilon_Tau[strainID], computed_load, false);
    }

    for (int i = 0; i < this->noIncl; i++)printf("%lf %lf %lf\n", this->inclusions[i]->Epsilon_Tau[strainID][0], this->inclusions[i]->Epsilon_Tau[strainID][1], this->inclusions[i]->Epsilon_Tau[strainID][2]);
    e(111);
  }
  else if (SBA_type == SBAT_VOID) {
    ;
  }
  else _errorr("Unknown SBAmode");
}


void Problem :: findAffectedInclusions (void)
{
  int *auxActingIncls = NULL;
  int i, j;
  
  // auxiliary memory allocation
  auxActingIncls = new int[noIncl];
  
  // finding the dependencies
  for (i=0; i<noIncl; i++) {
    // initial initialization of the number of acting inclusions
    inclusions[i]->noActingIncls = 0;
    for (j=0; j<noIncl; j++) {
      if (i != j  &&  inclusions[j]->point_is_affected(inclusions[i]->origin)) {
        auxActingIncls[inclusions[i]->noActingIncls] = j;
        inclusions[i]->noActingIncls++;
      }
    }
    
    if (this->inclusions[i]->noActingIncls != 0) {
      inclusions[i]->ActingIncls_allocate();
      CopyVector( auxActingIncls, inclusions[i]->actingIncls, inclusions[i]->noActingIncls );
    }
  }
  
  // deallocation of auxiliary memory
  delete [] auxActingIncls;
}

void Problem :: give_ovlivneni (const double *coords, double *result)
{
  // auxiliary 'Point' memory allocation
  Point point;
  point.set_nlc (this, 1);
  point.x[0] = coords[0];
  point.x[1] = coords[1];
  point.x[2] = coords[2];
  
  if (inclusions[0]->point_is_inside(point.x, 1.e-7)) {
    CopyVector(((InclusionRecord3D*) inclusions[0])->locEigStrain, result, 6);
  }
  else {
    CleanVector(((InclusionRecord3D*) inclusions[0])->locEigStrainTot, 6);
    ((InclusionRecord3D*) inclusions[0])->esuf->giveEshelbyStrainOfOnePoint(&point);
    ((InclusionRecord3D*) inclusions[0])->SBA_updateGlobAndLocStrains(&point);
    
    CopyVector(((InclusionRecord3D*) inclusions[0])->locEigStrainTot, result, 6);
  }
}

double Problem :: compute_supplement_energy (int lc)
{
  double E = 0.0;
  for (int i=0; i<noIncl; i++)
    E += this->inclusions[i]->compute_supplement_energy(lc);

  return E;
}



// ********************************************************************************************************
// *** *** *** ***                             PRINT FUNCTIONS                              *** *** *** ***
// ********************************************************************************************************

//*  ***   ELEMENTS   ***
void print_VTK_elems_head (Stream *stream, long noIncl)
{
  long i;
  
  if (stream->isFile()) {
    fprintf (stream->file(), "CELLS %ld %ld\n", noIncl, 2*noIncl);
    for (i=0; i<noIncl; i++)
      fprintf (stream->file(), "1 %ld\n", i);
    
    fprintf (stream->file(), "CELL_TYPES %ld\n", noIncl);
    for (i=0; i<noIncl; i++)
      fprintf (stream->file(), "1\n");
  }
  else {
    stream->tixel()->SetAttribute("NumberOfCells", (int)noIncl);
    
    XMLElement *cells = stream->tix_doc()->NewElement("Cells");   stream->tixel()->InsertEndChild(cells);
    
    XMLElement *dac = stream->tix_doc()->NewElement("DataArray");  cells->InsertEndChild(dac);
    dac->SetAttribute("format", "ascii");
    dac->SetAttribute("type", "Int32");
    dac->SetAttribute("Name", "connectivity");
    
    XMLElement *dao = (XMLElement*)dac->ShallowClone(NULL);   cells->InsertEndChild(dao);
    dao->SetAttribute("Name", "offsets");
    
    XMLElement *dat = (XMLElement*)dac->ShallowClone(NULL);   cells->InsertEndChild(dat);
    dat->SetAttribute("type", "UInt8");
    dat->SetAttribute("Name", "types");
    
    char auxs[7];
    sprintf (auxs,"0");
    for (long i=1; i<=noIncl; i++) {
      ;                         dac->InsertEndChild(stream->tix_doc()->NewText(auxs));
      sprintf (auxs,"%ld", i);  dao->InsertEndChild(stream->tix_doc()->NewText(auxs));
      ;                         dat->InsertEndChild(stream->tix_doc()->NewText("1" ));
    }
  }
}



//  ***   SCALAR / VECTOR / TENSOR  body   ***
XMLElement* print_VTK_field_head (Stream *stream, const char *name, char format, int slen, int noincl)
{
  XMLElement *da = NULL;
  
  char flint[8];
  if      (format == 'i') { if (stream->isFile()) sprintf(flint,"int");    else sprintf(flint,"Int32");   }
  else if (format == 'f') { if (stream->isFile()) sprintf(flint,"float");  else sprintf(flint,"Float32"); }
  else errol;
  
  if (stream->isFile()) {
    fprintf (stream->file(), "%s %d %d %s\n", name, noincl, slen, flint);
  }
  else {
    errol;
    //da = stream->tix_doc()->NewElement("DataArray");  stream->tixel()->InsertEndChild(da);
    //da->SetAttribute("format", "ascii");
    //da->SetAttribute("type", flint);
    //da->SetAttribute("Name", name);
    //if      (type == 's')  da->SetAttribute("NumberOfComponents", "1");
    //else if (type == 'v')  da->SetAttribute("NumberOfComponents", "3");
    //else if (type == 't')  da->SetAttribute("NumberOfComponents", "9");
    //else    _errorr ("error");
  }
  
  return da;
}



// ///
// void print_values (char *auxs, int n, const double *values, int rv)
// {
//   if (n==0) {
//     _errorr("");
//     auxs[0] = '\0';
//     return;
//   }
//   
//   //int n1 = 15;
//   int n1 = 14; // debug
//   int pos = 0;
//   
//   for (int i=0; i<n; i++) {
//     sprintf (auxs+pos, "%*.*e", n1, 8, values[i]);  pos += n1;
//     
//     if (auxs[pos] != '\0') pos++;
//     if (auxs[pos] != '\0') errol;
//     
//     if (i+1<n) {
//       if ((i+1)%rv == 0)  sprintf (auxs+pos, "\n");
//       else                sprintf (auxs+pos, " ");
//       pos++;
//     }
//   }
// }


void print_valuesi (char *auxs, int n, const int *values, int rv)
{
  if (n==0) {
    auxs[0] = '\0';
    return;
  }
  
  int pos = 0;
  
  for (int i=0; i<n; i++) {
    sprintf (auxs+pos, "%15d", values[i]);  pos += 15;
    
    if (auxs[pos] != '\0') pos++;
    if (auxs[pos] != '\0') errol;
    
    if (i+1<n) {
      if ((i+1)%rv == 0)  sprintf (auxs+pos, "\n");
      else                sprintf (auxs+pos, " ");
      pos++;
    }
  }
}


//
void Problem :: printVtkFileCompleteInclRec (const char *rsltFileName)
{
  long i, j;
  XMLElement *da;
  long auxl;
  Stream *stream;
  bool lgc = true;
  char auxs[20001];   auxs[20000] =  '\0';
  
  // allocate stream
  if (lgc)  stream = new Stream(STRM_file);
  else      stream = new Stream(STRM_tixel);
  
  //* START
  char *filename = new char[255];
  strcpy(filename, rsltFileName);// do not delete filename, it is deleted inside of the function
  //
  print_VTK_START (stream, filename);
  
  //* HEAD
  char type[8];
  if (this->twodim) sprintf (type,"2Dequiv");
  else              sprintf (type,"3Dequiv");
  
  if (lgc)
    print_VTK_head (stream->file(), type);
  
  //* NODES
  da = print_VTK_nodes_head (stream, this->noIncl);
  
  for (i=0; i<this->noIncl; i++) {
    sprintf (auxs, "%15.8e %15.8e %15.8e", this->inclusions[i]->origin[0], this->inclusions[i]->origin[1], this->inclusions[i]->origin[2]);
    
    if (lgc) fprintf (stream->file(), "%s\n", auxs);
    else     da->InsertEndChild( stream->tix_doc()->NewText(auxs) );
  }
  
  // ELEMENTS
  print_VTK_elems_head (stream, this->noIncl);
  
  
  // FIELDS
  fprintf (stream->file(), "FIELD unstructured_data %d\n", 4);
  // SBA
  da = print_VTK_field_head (stream, "SBA_mode", 'i', 1, 1);
  sprintf (auxs, "%d", 1 + (int)this->SBA_optimized);
  print_auxs (lgc, stream, da, auxs);
  // intFieldsShape
  da = print_VTK_field_head (stream, "conversion_mode", 'i', 1, 1);
  sprintf (auxs, "%d", this->intFieldsShape);
  print_auxs (lgc, stream, da, auxs);
  // Matrix record
  da = print_VTK_field_head (stream, "Matrix_record", 'f', 2, 1);
  sprintf (auxs   ,  "%12.6e", this->matrix->E());
  sprintf (auxs+12, " %12.6e", this->matrix->nu());
  print_auxs (lgc, stream, da, auxs);
  // Remote_strains
  da = print_VTK_field_head (stream, "Remote_strains", 'f', this->give_nLC(), this->give_TENS_RANGE());
  double rs[9];
  for (i=0; i<this->give_nLC(); i++) {
    CopyVector (this->matrix->Remote_strain[i], rs, this->give_TENS_RANGE());
    if (twodim)  MatrixOperations::convert2Dtensor_TRtoROW_notation (rs);
    else         MatrixOperations::convert3Dtensor_TFEEPtoROW_notation (rs);
    
    for (j=0; j<this->give_TENS_RANGE(); j++)
      sprintf (auxs+j*13, " %12.6e", rs[j]);
    
    print_auxs (lgc, stream, da, auxs);
  }
  
  // POINT DATA
  print_VTK_init_point_data (stream, this->noIncl);
  
  auxl = 1;                                    da = print_VTK_data_head (stream, "Inclusion_shape",                             's', 'i', auxl);  for (i=0; i<this->noIncl; i++) { sprintf (auxs, "%d",                         this->inclusions[i]->shape);          print_auxs (lgc, stream, da, auxs); }
  auxl = 1;                                    da = print_VTK_data_head (stream, "Youngs_modulus",                              's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { sprintf (auxs, "%.8e",                       this->inclusions[i] ->E      );       print_auxs (lgc, stream, da, auxs); }
  auxl = 1;                                    da = print_VTK_data_head (stream, "Poissons_ratio",                              's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { sprintf (auxs, "%.8e",                       this->inclusions[i] ->nu     );       print_auxs (lgc, stream, da, auxs); }
  auxl = this->give_VECT_RANGE();              da = print_VTK_data_head (stream, "Semiaxes_dimensions",                         's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs,  auxl,                   this->inclusions[i] ->a      );       print_auxs (lgc, stream, da, auxs); }
  auxl = this->give_EA_RANGE();                da = print_VTK_data_head (stream, "Euller_angles_rad",                           's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs,  auxl,                   this->inclusions[i]->eAngles );       print_auxs (lgc, stream, da, auxs); }
  auxl = 1 ;                                   da = print_VTK_data_head (stream, "Action_radius",                               's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { sprintf (auxs, "%.8e",                       this->inclusions[i]->actionRadius);   print_auxs (lgc, stream, da, auxs); }
  
  if (! this->twodim) {                        da = print_VTK_data_head (stream, "Acting_inclusions_number",                    's', 'i',    1);  for (i=0; i<this->noIncl; i++) { sprintf (auxs, "%d",                         this->inclusions[i]->noActingIncls);  print_auxs (lgc, stream, da, auxs); } }
  
  auxl = this->give_ISO_C_RANGE();             da = print_VTK_data_head (stream, "Stiffnesses_C",                               's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl,                    this->inclusions[i]->C);              print_auxs (lgc, stream, da, auxs); }
  
  if (! this->twodim) {                        da = print_VTK_data_head (stream, "Elliptic_potentials",                         's', 'f',   13);  for (i=0; i<this->noIncl; i++) { print_values (auxs, 13, ((InclusionRecord3D*)this->inclusions[i])->eInt);          print_auxs (lgc, stream, da, auxs); } }
  
  auxl = this->give_ISO_C_RANGE();             da = print_VTK_data_head (stream, "Eshelby_tensor",                              's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl,                    this->inclusions[i]->S);              print_auxs (lgc, stream, da, auxs); }
  
  if (! this->twodim) {                        da = print_VTK_data_head (stream, "Eshelby_tensor_inverse",                      's', 'f',   12);  for (i=0; i<this->noIncl; i++) { print_values (auxs, 12, ((InclusionRecord3D*)this->inclusions[i])->SInv);          print_auxs (lgc, stream, da, auxs); } }
  
  auxl = this->give_TRNSFM_MTRX_VEC_RANGE();   da = print_VTK_data_head (stream, "Global_2_local_transform_matrix_for_vectors", 's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl,                    this->inclusions[i]->T);              print_auxs (lgc, stream, da, auxs); }
  //   = this->give_TRNSFM_MTRX_VEC_RANGE();   da = print_VTK_data_head (stream, "Local_2_global_transform_matrix_for_vectors", 's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl,                    this->inclusions[i]->TInv);           print_auxs (lgc, stream, da, auxs); }
  auxl = this->give_TRNSFM_MTRX_TENS_RANGE();  da = print_VTK_data_head (stream, "Global_2_local_transform_matrix_for_tensors", 's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl,                    this->inclusions[i]->Te);             print_auxs (lgc, stream, da, auxs); }
  auxl = this->give_TRNSFM_MTRX_TENS_RANGE();  da = print_VTK_data_head (stream, "Local_2_global_transform_matrix_for_tensors", 's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl,                    this->inclusions[i]->TeInv);          print_auxs (lgc, stream, da, auxs); }
  
  
  int nlc = this->give_nLC();

  if (this->twodim) {
    auxl = this->give_VM_TENS_RANGE();
    for (j=0; j<nlc; j++) { sprintf (auxs, "Local_Eq_strain_%02ld_%02d", j+1, nlc);   da = print_VTK_data_head (stream, auxs,   's', 'f', auxl);  for (i=0; i<this->noIncl; i++) { print_values (auxs, auxl, this->inclusions[i]->Epsilon_Tau[j]);         print_auxs (lgc, stream, da, auxs); }  }
  }
  else {
    for (j=0; j<nlc; j++) { sprintf (auxs, "Local_Eq_strain_%02ld_%02d", j+1, nlc);   da = print_VTK_data_head (stream, auxs,   's', 'f',    6);  for (i=0; i<this->noIncl; i++) { print_values (auxs,    6, ((InclusionRecord3D*)this->inclusions[i])->locEigStrain_LC[j]);   print_auxs (lgc, stream, da, auxs); }  }
    
    for (j=0; j<nlc; j++) { sprintf (auxs, "Global_Eq_strain_%02ld_%02d", j+1, nlc);  da = print_VTK_data_head (stream, auxs,   's', 'f',    6);  for (i=0; i<this->noIncl; i++) { print_values (auxs,    6, ((InclusionRecord3D*)this->inclusions[i])->globEigStrain_LC[j]);  print_auxs (lgc, stream, da, auxs); }  }
    
    // FIELDS
    fprintf (stream->file(), "FIELD unstructured_data %d\n", 1);
    // Acting inclusions
    auxl = this->totalNoActingIncls;
    da = print_VTK_field_head (stream, "Acting_inclusions_list",                      'i', auxl, 1);
    for (i=0; i<this->noIncl; i++) {
      print_valuesi (auxs, this->inclusions[i]->noActingIncls, this->inclusions[i]->actingIncls, 10);
      print_auxs (lgc, stream, da, auxs);
    }
  }
  
  // END
  print_VTK_FINISH (stream);
  delete stream;
  
}//end of function: printVtkFileCompleteInclRec()

long Problem :: give_inclusion_with_point_inside (const double *coords, double epsilon) const
{
  long i;
  for (i=0; i<noIncl; i++)
    if (inclusions[i]->point_is_inside(coords, epsilon) == true) {
      //printf("inside: %d\n",i);
      return i;
    }
  
  return -1;
}


void convertTensorStrain (bool twodim, double *tens, STRNotation oldNot, STRNotation newNot)
{
  //
  if (oldNot == STRN_THEORETICAL_FEEP) {
    if      (newNot == STRN_THEORETICAL_FEEP) { ; }
    else if (newNot == STRN_THEORETICAL_ROW) {
      if (twodim)  MatrixOperations::convert2Dtensor_TRtoROW_notation (tens);
      else         MatrixOperations::convert3Dtensor_TFEEPtoROW_notation (tens);
    }
    else if (newNot == STRN_VOIGT_FEEP) {
      if (twodim)  MatrixOperations::convert2Dtensor_TRtoVR_notation(tens);
      else         MatrixOperations::convert3Dtensor_TRtoVR_notation(tens);
    }
    else _errorr ("not supported notation");
  }
  else
    _errorr ("not supported notation");
}
void convertTensorStress (bool twodim, double *tens, STRNotation oldNot, STRNotation newNot)
{
  //
  if (oldNot == STRN_THEORETICAL_FEEP) {
    if      (newNot == STRN_THEORETICAL_FEEP) { ; }
    else if (newNot == STRN_THEORETICAL_ROW) {
      if (twodim)  MatrixOperations::convert2Dtensor_TRtoROW_notation (tens);
      else         MatrixOperations::convert3Dtensor_TFEEPtoROW_notation (tens);
    }
    else if (newNot == STRN_VOIGT_FEEP) { ; }
    else _errorr ("not supported notation");
  }
  else
    _errorr ("not supported notation");
}



long Problem :: giveFieldsOfPointOneRS (double *displc, double *strain, double *stress, const double *coords,
                                        char ptFlag, int rs, PFCmode pfcMode, long reqIncl, STRNotation tn) const
{
  return this->giveFieldsOfPoint ( displc ? &displc : NULL, strain ? &strain : NULL, stress ? &stress : NULL, coords, ptFlag, rs, 1, pfcMode, reqIncl, tn);
}

long Problem :: giveFieldsOfPoint (double **displc, double **strain, double **stress, const double *coords,
                                   char ptFlag, int rs, int nrs, PFCmode pfcMode, long reqIncl, STRNotation tn) const
{
  // perturbation fields
  long parent_incl = this->give_EshelbyPertFieldsOnePoint (coords, displc, strain, stress, rs, nrs, pfcMode, reqIncl);
  
  // add homogeneous fields in case of total fields
  if (ptFlag == 't') {
    for (int i=0; i<nrs; i++) {
      // pro kontrolu dilute
      //if (parent_incl == -1) {
      //  CleanVector(strain[i], this->give_VM_TENS_RANGE());
      //  CleanVector(stress[i], this->give_VM_TENS_RANGE());
      //}
      if (displc) AddVector (this->matrix->give_globHomog_Displc(i+rs, coords), displc[i], this->give_VECT_RANGE());
      if (strain) AddVector (this->matrix->give_globHomog_Strain(i+rs),         strain[i], this->give_VM_TENS_RANGE());
      if (stress) AddVector (this->matrix->give_globHomog_Stress(i+rs),         stress[i], this->give_VM_TENS_RANGE());
    }
  }
  else if (ptFlag != 'p') _errorr("wrong pt flag");
  
  // convert to required notation
  if (strain)  for (int i=0; i<nrs; i++)  convertTensorStrain(twodim, strain[i], STRN_THEORETICAL_FEEP, tn);
  if (stress)  for (int i=0; i<nrs; i++)  convertTensorStress(twodim, stress[i], STRN_THEORETICAL_FEEP, tn);
  
  return parent_incl;
}


long Problem :: give_EshelbyPertFieldsOnePoint (const double *coo, double **displc, double **strain, double **stress,
                                                int lc, int nlc, PFCmode pfcMode, long reqIincl) const
{
  double *coords;
  // shift of coordinates od inclusion set
  if (this->matrix->origin) {
    coords = new double[3];
    SubtractTwoVectors (coo, this->matrix->origin, coords, this->give_VECT_RANGE());
  }
  else
    coords = (double*)coo;
  
  nlc = this->check_lc_nlc(lc, nlc);
  
  // point position
  long parent_incl;
  
  // no supposal of point position
  if (reqIincl == -3) {
    parent_incl = this->give_inclusion_with_point_inside (coords, 0.0);
  }
  // point is supposed to lay on any inclusion
  else if (reqIincl == -2) {
    parent_incl = this->give_inclusion_with_point_inside (coords, 1.e-6);
    if (parent_incl == -1)
      _errorr("point doesnot lay on inclusion, nekde asi zustala -2, ted tam ma byt -3 pro neurcitou polohu");
  }
  // point is supposed to lay on matrix
  else if (reqIincl == -1) {
    parent_incl = this->give_inclusion_with_point_inside (coords, -1.e-6);
    if (parent_incl != -1)
      _errorr("point doesnot lay on matrix");
  }
  // point is supposed to lay on an inclucion
  else if (reqIincl >= 0) {
    parent_incl = reqIincl;
    if (inclusions[parent_incl]->point_is_inside (coords, 1.e-3) == false)
      _errorr("point doesnot lay inside of the supposed inclusion");
  }
  else errol;
  
  // *** EVALUATION OF PERTURBATION FIELDS ***
  if (displc != NULL  ||  strain != NULL  ||  stress != NULL) {
    //parent_incl = -1; // pro debug hrany NEMAZAT
    
    
    double **strain_r;
    double **e_t_updated;
    
    
    // termitovo, na toto si udelat prikladek o dvou inkluzich, a porovnat grafy
    // zaroven poradne pochopit standovo balancing
    // pak vykostit starej balancing, asi
    // intFieldsShape:
    // 0 = nekonstanta,   s pricinku vnejsich inkluzi   redukovana
    // 1 = polynomial - if available
    // 2 =   konstanta,   s pricinku vnejsich inkluzi   redukovana
    // 3 = nekonstanta,   s pricinku vnejsich inkluzi neredukovana
    // 4 =   konstanta, bez pricinku vnejsich inkluzi
    
    
    if (parent_incl == -1) {
      // EXTERIOR POINT
      this->give_EshelbyPertFieldsOnePoint_external (coords, displc, strain, stress, lc, nlc, pfcMode, parent_incl);
    }
    else {
      // INTERIOR POINT
      
      // SUM EXT AND RECALCULATE, THEN ADD INT
      if (intFieldsShape == 0) {
    // I need to compute strain even if strain==NULL
        if (strain) strain_r = strain;
    else        strain_r = AllocateArray2D(nlc, this->give_VM_TENS_RANGE());
        
    // add ext
    this->give_EshelbyPertFieldsOnePoint_external(coords, displc, strain_r, stress, lc, nlc, pfcMode, parent_incl);
        
        // field for recalculated eps_tau
    e_t_updated = AllocateArray2D(nlc, this->give_VM_TENS_RANGE());
    
    if (strain || stress) {
      // computing the eps_tau
          
          // termitovo a bude to tady vubec fungovat kdyz nebude spustenej balancing??
          // termitovo takze tady se v danem bode spocita strain_r jako vnejsi zatizeni a pak z predpokladu konstantniho prubehu se spocte epstau a z nej strain
          // ale kdyz neni vnejsi zatizeni konstanta, tak by tady spis mel byt polynom, ne?
          
          // termitovo proc je tady true, tim se pocita na inkluzi strain_r plus remote_strain, neni lepsi dat false a nakonci secist s strainGlob...
      for (int i = 0; i < nlc; i++) ((InclusionRecord2D*)inclusions[parent_incl])->compute_eps_tau(i+lc, e_t_updated[i], strain_r[i], true);
          // termitovo jaktoze je tady add, nema to neco spolecneho s CHYBOU kde nesouhlasi strain se steessem?
          // musim si prostudovat balancing
      this->inclusions[parent_incl]->add_EshelbyPertStrain_internal_SIFCM (strain_r, lc, nlc, e_t_updated);
          
      if (stress) {
        CleanArray2d(stress, nlc, this->give_VM_TENS_RANGE());
        this->inclusions[parent_incl]->add_EshelbyPertStress_internal_SIFCM(stress, lc, nlc, e_t_updated, strain_r);
      }
    }
        
        if (!strain) DeleteArray2D (strain_r, nlc);
        
        DeleteArray2D (e_t_updated, nlc);
        
        if (displc) {
      //CleanArray2d(displc, nlc, this->give_VECT_RANGE());
          // add bylo uvnitr fce, vytahl jsem to nahoru at je to videt
          // termitovo proc je tady add??
          AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertDisplc_internal(lc, nlc, coords), displc, nlc, this->give_VECT_RANGE());
    }
      }
      
      // POLYNOMIAL INT FIELDS
      if (intFieldsShape == 1) {
    _errorr("Not implemented yet.");
      }
      
      // CONST INT STRAIN AND STRESS
      if (intFieldsShape == 2) {
        this->give_EshelbyPertFieldsOnePoint_external(this->inclusions[parent_incl]->origin, displc, strain, stress, lc, nlc, pfcMode, parent_incl);
        
        if (strain)  AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertStrain_internal(lc, nlc)        , strain, nlc, this->give_VM_TENS_RANGE());
    if (stress)  AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertStress_internal(lc, nlc)        , stress, nlc, this->give_VM_TENS_RANGE());
    if (displc)  AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertDisplc_internal(lc, nlc, coords), displc, nlc, this->give_VECT_RANGE()   );
      }
      
      // SUM ALL
      if (intFieldsShape == 3) {
        this->give_EshelbyPertFieldsOnePoint_external(coords, displc, strain, stress, lc, nlc, pfcMode, parent_incl);
        
        if (strain)  AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertStrain_internal(lc, nlc)        , strain, nlc, this->give_VM_TENS_RANGE());
    if (stress)  AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertStress_internal(lc, nlc)        , stress, nlc, this->give_VM_TENS_RANGE());
    if (displc)  AddArray2D (this->inclusions[parent_incl]->give_EshelbyPertDisplc_internal(lc, nlc, coords), displc, nlc, this->give_VECT_RANGE()   );
      }
      
      // CONST INT STRAIN AND STRESS - only the current inclusion
      if (intFieldsShape == 4) {
        if (strain)  CopyArray2D (this->inclusions[parent_incl]->give_EshelbyPertStrain_internal(lc, nlc)        , strain, nlc, this->give_VM_TENS_RANGE());
    if (stress)  CopyArray2D (this->inclusions[parent_incl]->give_EshelbyPertStress_internal(lc, nlc)        , stress, nlc, this->give_VM_TENS_RANGE());
    if (displc)  CopyArray2D (this->inclusions[parent_incl]->give_EshelbyPertDisplc_internal(lc, nlc, coords), displc, nlc, this->give_VECT_RANGE()   );
      }
      
    }
    
  }
  
  
  // dealloc
  if (this->matrix->origin)  delete [] coords;
  else coords = NULL;
  
  return parent_incl;
}


void Problem :: give_EshelbyPertFieldsOnePoint_external (const double *coords, double **displc, double **strain, double **stress,
                                                         int lc, int nlc, PFCmode pfcMode, long parent_incl) const
{
  Inclusion *incl;
  
  // affected inclusions
  long noActIncls;
  long *actIncls = NULL;
  if (pfcMode == PFCM_FULL) {
    noActIncls = this->noIncl;
  }
  else if (pfcMode == PFCM_OPTIMIZED) {
    // creating of auxiliary fields
    actIncls = new long[this->noIncl];
    //searching for acting inclusions
    noActIncls = this->giveActingInclusions(actIncls, coords);
  }
  else  _errorr("giveEshelbyPertFieldsOfOnePoint: unknofn PFCmode");
  
  // tenhle point mozna alokovat na Problemu at se furt nealokuje znovu
  Point *point = new Point;
  point->set_nlc(this, nlc);
  CopyVector (coords, point->x, this->give_VECT_RANGE());
  
  bool delstrain = false;
  // initialization fields
  if (stress && strain == NULL) {
    delstrain = true;
    AllocateArray2D (strain, nlc, this->give_VM_TENS_RANGE());
  }
  
  if (displc)  CleanArray2d (displc, nlc, this->give_VECT_RANGE());
  if (strain)  CleanArray2d (strain, nlc, this->give_VM_TENS_RANGE());
  
  for (int i=0; i<noActIncls; i++) {
    if (parent_incl == i) continue;
    
    // affected inclusion
    if (pfcMode == PFCM_FULL) incl = this->inclusions[i];
    else                      incl = this->inclusions[actIncls[i]];
    
    // calculation of perturbation fields
    incl->give_EshelbyPertFields_external (point, lc, nlc, displc, strain);
    
    // superposition of perturbation fields from all inclusions
    if (displc) AddArray2D (point->displacement, displc, nlc, this->give_VECT_RANGE()   );
    if (strain) AddArray2D (point->strain,       strain, nlc, this->give_VM_TENS_RANGE());
  }
  
  // evaluation of stress perturbation
  if (stress) {
    for (int i=0; i<nlc; i++)
      if (twodim)  MatrixOperations::giveTVproduct_3is3x3to5and3  (stress[i], this->matrix->C, strain[i]);
      else         MatrixOperations::giveTVproduct_6is6x6to12and6 (stress[i], this->matrix->C, strain[i]);
  }
  
  // deallocation
  if (pfcMode == PFCM_OPTIMIZED)  delete [] actIncls;
  if (delstrain) DeleteArray2D (strain, nlc);
  delete point;
  
}//end of function: addEshelbyPertFieldsOfOnePoint( )

void Problem :: printFieldsOnMeshVTK (const char *mesh_file_out, const char *mesh_file,
                                   char ptFlag, int rs, int nrs,
                                   PFCmode pfcMode) const
{
  nrs = this->check_lc_nlc(rs, nrs);

  Mesh mesh(1, this);
  mesh.read_geometry_file_vtk (mesh_file, -1, -1, false);
  mesh.compute_node_fields(ptFlag, true, true, true, 0, rs, nrs, pfcMode);
  mesh.print_geometry_file_vtk (mesh_file_out, 0, rs, nrs);
}

void Problem :: printFieldsOnMeshGrid (const char *mesh_file_out, const double *p1, const double *p2, const long *n,
                                       char ptFlag, int rs, int nrs,
                                       PFCmode pfcMode) const
{
  nrs = this->check_lc_nlc(rs, nrs);

  Mesh mesh(1, this);
  mesh.generate_regular_mesh (p1, p2, n);
  mesh.compute_node_fields(ptFlag, true, true, true, 0, rs, nrs, pfcMode);
  mesh.print_geometry_file_vtk (mesh_file_out, 0, rs, nrs);
}






Problem *Problem::read_inclusions_plus_initialize_and_print (const char *inclusion_file, const char *equiv_file, DiffTypes dt)
{
  this->set_diffType (dt);
  this->read_input_file (inclusion_file);
  
  this->input_data_initialize_and_check_consistency();
  
  if (this->data_equivalent == false)
    this->convert_to_equivalent_problem ();
  
  if (equiv_file)
    this->print_equivalent_problem (equiv_file);

  return this;
}

Mesh* Problem :: give_new_mesh (void)
{
  meshes[nmeshes] = new Mesh(nmeshes, this);
  nmeshes++;
  return meshes[nmeshes-1];
}

Homogenization* Problem :: give_new_homogenization (HomogenizationType ht)
{
  switch( ht )
  {
  case HT_Dilute: {
    homogs[nhomogs] = new Dilute(nhomogs, this);
    break;
  }
  case HT_MoriTanaka: {
    homogs[nhomogs] = new MoriTanaka(nhomogs, this);
    break;
  }
  case HT_ReussBound: {
    homogs[nhomogs] = new ReussBound(nhomogs, this);
    break;
  }
  case HT_VoightBound: {
    homogs[nhomogs] = new VoightBound(nhomogs, this);
    break;
  }
  case HT_RegGrid: {
    homogs[nhomogs] = new RegGrid(nhomogs, this);
    break;
  }
  case HT_DifferentialScheme: {
    homogs[nhomogs] = new DifferentialScheme(nhomogs, this);
    break;
  }
  default:
    _errorr( "Unknown homogenization type" );

  }

  nhomogs++;
  return homogs[nhomogs-1];
}

void Problem :: print_equivalent_problem (const char *vtkTopologyFileRec)
{
  // gives the total value of Kn number. Kn is the numbre of inclusion within the action region of each inclusion
  this->totalNoActingIncls = 0;
  for (long i=0; i<noIncl; i++)
    this->totalNoActingIncls += inclusions[i]->noActingIncls;
  
  // printing the complete inclusion topology and geometry file containing equivalent strains
  this->printVtkFileCompleteInclRec (vtkTopologyFileRec);
}



//void Problem :: compute_fields_and_cmp_with_FEM_on_mesh (const char *mesh_file, const char *fem_results, int lc, int nlc, PFCmode pfcMode)
//{
 //   nlc = P->check_lc_nlc(lc, nlc);
  // 
  // Mesh mesh(1, this);
  // mesh.read_geometry_file_vtk (mesh_file);
  // mesh.compute_fields(lc, nlc, pfcMode);
  // mesh.print_geometry_file_vtk (mesh_file_out, lc, nlc);
//}





//********************************************************************************************************
//                                            PROTECTED FUNCTIONS
//********************************************************************************************************


//********************************************************************************************************
// description: function gives a set of acting inclusions to a given point. Function returns the number
//              of acting inclusions. Function is a clon of 'inclRecordInit::giveAffectedInclusions'
//              method.
// last edit:   17. 03. 2010
// -------------------------------------------------------------------------------------------------------
// note:        @actIncls - pointer to an empty array on noIncls size, i.e.: actIncls[noIncls]
//              @coords - coordinates of a point to be checked out
//********************************************************************************************************
long Problem :: giveActingInclusions (long *actIncls, const double *coords) const
{
  long noActIncls = 0;
  
  //finding the dependencies
  for (long i=0; i<this->noIncl; i++)
    if (inclusions[i]->point_is_affected(coords)) {
      actIncls[noActIncls] = i;
      noActIncls++;
    }
  
  return noActIncls;
}//end of function: giveActingInclusions( )


double Problem :: updateEpsTauInSBal (int strainID, double **&last_eps_tau)
{
  vektor suma_e_s(this->give_VM_TENS_RANGE());
  double point[3];
  double **strain = AllocateArray2D(1, this->give_VM_TENS_RANGE());
  for (int i = 0; i < noIncl; i++) {
    for (int ap = -1; ap < this->inclusions[i]->n_approx_points; ap++) {
      if (ap == -1)for (int k = 0; k < 3; k++) point[k] = inclusions[i]->origin[k];
      else for (int k = 0; k < 3; k++) point[k] = inclusions[i]->approx_points[ap][k];
      suma_e_s.vynulovat();
      if (SBA_optimized) {
    for (int j = 0; j < inclusions[i]->noActingIncls; j++) {                                                                            
      ((InclusionRecord2D*)inclusions[inclusions[i]->actingIncls[j]])->giveExtStrainPert(point, strain, strainID, 1);                                   
      for (int ii = 0; ii < this->give_VM_TENS_RANGE(); ii++)suma_e_s.v[ii] += strain[0][ii];
    }
      }
      else {
    for (int j = 0; j < noIncl; j++) {
      if (i != j) {
        ((InclusionRecord2D*)inclusions[j])->giveExtStrainPert(point, strain, strainID, 1);
        for (int ii = 0; ii < this->give_VM_TENS_RANGE(); ii++)suma_e_s.v[ii] += strain[0][ii];
      }
    }
      }
      for (int j = 0; j < this->give_VM_TENS_RANGE(); j++)suma_e_s.v[j] += this->matrix->remote_strains()[strainID][j];

      if (ap == -1)this->inclusions[i]->Eps02EpsTau(this->inclusions[i]->Epsilon_Tau[strainID], suma_e_s.v);
      else this->inclusions[i]->Eps02EpsTau(this->inclusions[i]->approx_points_eps_tau[ap], suma_e_s.v);
    }
    if (this->inclusions[i]->n_approx_points)this->inclusions[i]->update_approximations(strainID);
  }



  double norma = 0.;
  if (last_eps_tau == NULL) {
    last_eps_tau = AllocateArray2D(this->noIncl, this->give_VM_TENS_RANGE());
    for (int i = 0; i < this->noIncl; i++)for (int j = 0; j < this->give_VM_TENS_RANGE(); j++) {
    last_eps_tau[i][j] = this->inclusions[i]->Epsilon_Tau[strainID][j];
      }
    norma = 100.;
  }
  else {
    for (int i = 0; i < this->noIncl; i++) {
      for (int j = 0; j < this->give_VM_TENS_RANGE(); j++)
    norma += SQR(this->inclusions[i]->Epsilon_Tau[strainID][j] - last_eps_tau[i][j]);

      for (int j = 0; j < this->give_VM_TENS_RANGE(); j++)
    last_eps_tau[i][j] = this->inclusions[i]->Epsilon_Tau[strainID][j];
    }
  }
  return norma;
}


void Problem :: print_visualization (const char *name, int n, int dim, bool refined)
{
  if (this->data_equivalent == false)  _errorr("Inverse matrix TInv at inclusion have to be computed before calling print_geom function, for example call initialize_inclusion_records");
  
  if (dim == 0)  dim = this->give_dimension();

  if (refined)  _errorr("refined je nejaky pokus, je rozbitej momentalne, ale mozna jen ve 2d");

  if (dim != this->give_dimension())
    _errorr("error in visualization of 3d problem in xy plane, still not repaired");

  if (dim == 2) {
    // unit mesh
    Mesh umesh(1, NULL);
    umesh.generate_regularSphereMesh_2d(n, NULL);
    
    long nn = umesh.nNodes;
    long ne = umesh.nElems;
    
    // Total mesh
    Mesh mesh(1, this);
    mesh.set_coll(3);
    mesh.nNodes = noIncl * nn;  mesh.Nodes = new mNode*   [mesh.nNodes];
    mesh.nElems = noIncl * ne;  mesh.Elems = new mElement*[mesh.nElems];
    
    // auxiliary mesh of one inclusion
    Mesh *mesh1;
    const InclusionRecord2D *incl;
    
    // loop over all inclusions
    for (int i=0; i<noIncl; i++) {
      incl = (InclusionRecord2D*)inclusions[i];
      if (!refined) mesh1 = new Mesh(1, this, &umesh);
      else        { mesh1 = new Mesh(1, this);
        errol; // tady by melo byt asi 2d a ne 3d
        mesh1->generate_regularSphereMesh_3d(n, incl->a);
      }
      
      mesh1->scale(incl->a);
      mesh1->local2global(incl->origin, incl->TInv);
      mesh1->shift_id(i*nn, i*ne);
      
      for (int j=0; j<nn; j++) {
    mesh.Nodes[i*nn+j] = mesh1->Nodes[j];
    mesh.Nodes[i*nn+j]->M = &mesh;
      }
      delete [] mesh1->Nodes; mesh1->Nodes = NULL;
      
      for (int j=0; j<ne; j++) {
    mesh.Elems[i*ne+j] = mesh1->Elems[j];
    mesh.Elems[i*ne+j]->M = &mesh;
      }
      delete [] mesh1->Elems; mesh1->Elems = NULL;
      
      delete mesh1;
    }
    
    mesh.print_geometry_file_vtk (name, 0, 0, 1);
  }
  else {
    if (this->twodim) _errorr("print_geom problem is 3d, please use parameter dim 2");
    
    // unit mesh
    Mesh umesh(1, NULL);
    umesh.generate_regularSphereMesh_3d (n, NULL);
    
    long nn = umesh.nNodes;
    long ne = umesh.nElems;
    
    // Total mesh
    Mesh mesh(1, this);
    mesh.set_coll(0);
    mesh.nNodes = noIncl * nn;  mesh.Nodes = new mNode*   [mesh.nNodes];
    mesh.nElems = noIncl * ne;  mesh.Elems = new mElement*[mesh.nElems];
    
    // auxiliary mesh of one inclusion
    Mesh *mesh1;
    const InclusionRecord3D *incl;
    
    // loop over all inclusions
    for (int i=0; i<noIncl; i++) {
      incl = (InclusionRecord3D*)inclusions[i];
      if (!refined) mesh1 = new Mesh(1, this, &umesh);
      else        { mesh1 = new Mesh(1, this);
        mesh1->generate_regularSphereMesh_3d (n, incl->a);
      }
      
      mesh1->scale(incl->a);
      mesh1->local2global(incl->origin, incl->TInv);
      mesh1->shift_id(i*nn, i*ne);
      
      for (int j=0; j<nn; j++) {
    mesh.Nodes[i*nn+j] = mesh1->Nodes[j];
    mesh.Nodes[i*nn+j]->M = &mesh;
      }
      delete [] mesh1->Nodes; mesh1->Nodes = NULL;
      
      for (int j=0; j<ne; j++) {
    mesh.Elems[i*ne+j] = mesh1->Elems[j];
    mesh.Elems[i*ne+j]->M = &mesh;
      }
      delete [] mesh1->Elems; mesh1->Elems = NULL;
      
      delete mesh1;
    }
    
    mesh.print_geometry_file_vtk (name, 0, 0, 1);
  }
}

void Problem :: check_overlap (void)
{
  long i, j, status;

  for (i=0; i<this->noIncl; i++) {
    for (j=i+1; j<this->noIncl; j++) {
      status = this->inclusions[i]->find_overlap(this->inclusions[j]);

      if (status == -1) {
        //
        _errorr("check_overlap: NOT converged");
      }

      if (status == 1) {
        //
        _errorr("check_overlap: OVERLAP");
      }
    }
  }
}



} // end of namespace mumech

/*end of file*/