tests_tools.cpp

Go to the documentation of this file.

#include "tests.h"

#include "problem.h"
#include "mesh.h"
#include "homogenizationMethods.h"
#include "comparison.h"
#include "types.h"
#include "mnode.h"
#include "melement.h"
#include "vtk.h"
#include "mesoface.h"

#include "gelib.h"
#include "librw.h"

#include <algorithm>

using namespace mumech;

struct Bod {
  int ID;
  double coord[3];
  int mat[2];
  double rslt[2][6];
  double disp[4];
  bool last=false;
};

// *******************************************   ULOHY   ********************************************

// ***  2d ulohy pocitane 2d a 3d mumechem a ansysem, tisk rezu napeti

// Jedna inkluze.
void x2d_1I (void);
// Dve 2d elipticke inkluze na ose x, muze se menit jejich vzdalenost 'dx'.
void x2d_2I (void);
// Pravidelny grid kruhovych 2d inkluzi poctu "2n+1" x "2n+1".
void x2d_grid (void);


// ***  ulohy do clanku o mumechu

// vzorovy 3d priklad o dvou inkluzich
void clanek_3d_vzor (void);
// homogenizace - 2d uloha, pravidelny grid kruhovych inkluzi "2n+1" x "2n+1".
void clanek_homog_grid (void);


// ***  ruzne

// testovani visualizace
void visualization_1I_2d (void);
void visualization_1I_3d (void);
void visualization_197I_2d3d (void);


void intersection_2d (void);
void intersection_3d (void);

// 2d uloha pocitana 2d a 3d mumechem, postupna rotace inkluze, tisk hodnot v jednom bode; NUMERIKA je zubata!!
void x2Dx3D_rotation (void);

void penergy_2d (void);

// ---

void hrana_inkluze (void);

void Energy_x2D_1I_MMxFEM (void);

void ovlivneni (void);

// nove testy od standy, speedtests, zatim nezpracovane
void speedtest2Dx3D(void);
void createNodes(double a, double b, double c, double da, double db, double dc, double **&points, int &n);
void test3D();


// **************************************   DROBNE FUNKCE   **************************************
void modeset (Problem *p, int mode);
// **************************************   KOMPLEXNI FUNKCE   **************************************
void print_graf (const char *vtkGeomFile, const char *path, const char *suff, int lc, char axe, int comp, double coordMin, double coordMax, double coordDelta, char pt_flag, double dshift, DiffTypes dt, Bod * b = NULL, bool odd = false);
void compute_and_print_reg_mesh (const char *vtkGeomFile, const char *vtkMeshFileOutPert, double x, int n, char pt, DiffTypes dt);

// ************************************   ANSYS POSTPROCESSING   ************************************
Bod* read_ansys_strain_stress_plot (const char *path, int lc, char axe, int comp, double coordMin, double coordMax, bool strn, bool half=false, bool bod=false);
void read_ansys_displ_plot (const char *path, int lc, char axe, int comp, double coordMin, double coordMax);

// **********************************   GENEROVANI EQIV SOUBORU   ***********************************
void eqiv_gener_1x2dI_3d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double y, double E, double nu, double a1, double a2, double ea, DiffTypes dt);
void eqiv_gener_1x2dI_2d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double y, double E, double nu, double a1, double a2, double ea, DiffTypes dt);
void eqiv_gener_2x2dI_2d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double E, double nu, double a1, double a2, DiffTypes dt, int mode=0);
void eqiv_gener_2x2dI_3d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double E, double nu, double a1, double a2, DiffTypes dt, int mode=0);
void eqiv_gener_gx2dI_2d (const char *vtkAmndFile, int gx, int gy, double dx, double dy, double Em, double num, const double *rs, double E, double nu, double a1, double a2, DiffTypes dt, int mode);
void eqiv_gener_gx2dI_3d (const char *vtkAmndFile, int gx, int gy, double dx, double dy, double Em, double num, const double *rs, double E, double nu, double a1, double a2, DiffTypes dt, int mode);

void eqiv_gener_1x3dI (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double y, double z, double E, double nu,
               double a1, double a2, double a3, double ea1, double ea2, double ea3, DiffTypes dt);


void tests_tools (void)
{
  //x2d_1I ();
  x2d_2I ();
  //x2d_grid ();
  
  // visualization ();

  //intersection_2d();
  //intersection_3d ();

  //clanek_3d_vzor ();
  //clanek_homog_grid ();

  // ---

  //x2Dx3D_rotation();

  //penergy_2d ();

  

  // *** hrana inkluze ext x int
  //hrana_inkluze ();
  
  //print_incls ();
  //print_incls_2d ();
}

// **************************************************************************************************
// *******************************************   ULOHY   ********************************************
// **************************************************************************************************
void compare_grafs (const char *name1, const char *name2)
{
  FILE *file1 = fopen(name1, "r");
  FILE *file2 = fopen(name2, "r");
  char c;
  int n=0;
  double v1, v2;
  bool bre=false;
  
  while (true) {
    c = fgetc(file1); ungetc(c,file1);
    if (c == '\n' || c == EOF) break;
    
    fscanf (file1, "%lf", &v1);
    fscanf (file2, "%lf", &v2);
    if (v1 != v2) errol;
    
    fscanf (file1, "%lf", &v1);
    fscanf (file2, "%lf", &v2);
    if (! areSame(1e-9, 0.01, v1, v2)) {
      printf ("%lf   %lf \n", v1, v2);
      printf ("cd .. ; xmgrace -param /home/termit/.config/grace/link.par %s %s ; cd -\n", name1, name2);
      bre = true;
      break;//warningerrol;
    }
    
    FP_skip_line(file1);
    FP_skip_line(file2);
    n++;
  }
  if (bre==false && n<20)
    errol;
  
  fclose(file1);
  fclose(file2);
}



void penergy_2d (void)
{
  Comparison :: give_energy_MM ();

  return;
}




void x2d_1I (void)
{
  Bod* b = NULL;
  //char suff[100];
  char vtkAmndFile2d [] = "x2d_1I/amnd.2d.vtk";
  char vtkAmndFile3d [] = "x2d_1I/amnd.3d.vtk";
  
  // **************************************
  // ZDE URCIT, CO SE BUDE POCITAT/TISKNOUT
  // **************************************
  double angle = 30.0; // 0.0 or 30.0  !! Kdyz zde zmenis uhel a chces porovnat s ANSYSem, tak zmen link viz vyse.
  bool pert = false;
  bool totl = true;
  char xy[] = "XY"; // nebo treba "X-" "-Y"
  int nlc = 3;      if (nlc !=1 && nlc != 3) errol;  // nlc musi byt jen 1 nebo 3 aby se to mohlo automaticky generovat, tj. rucne zadany rs23 nebo vsechno
  int ncomp = 3;    if (ncomp<1 || ncomp>3)  errol;
  
  double rs2 [] = { 1.0, 0.0, 0.0, 0.0 };
  double rs3 [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  
  bool m2d = true;
  bool m3d = true;
  
  // *******
  // COMPUTE
  // *******
  
  // vytvoreni 2d a 3d prikladu
  if(m2d) eqiv_gener_1x2dI_2d (vtkAmndFile2d, 1.0, 0.2, nlc == 1 ? rs2 : NULL, 0.0, 0.0, 10.0, 0.4, 1.0, 0.5, angle, DT_NUMERICAL);
  if(m3d) eqiv_gener_1x2dI_3d (vtkAmndFile3d, 1.0, 0.2, nlc == 1 ? rs3 : NULL, 0.0, 0.0, 10.0, 0.4, 1.0, 0.5, angle, DT_NUMERICAL);
  
  //char name1[100], name2[100], name3[100];
  // loop over nlc (number of load cases)
  for (int i=0; i<nlc; i++) {
    // loop over ncomp (number of components)
    for (int j=0; j<ncomp; j++) {
      // loop over axes
      for (int k=0; k<2; k++) {
    if (xy[k] == '-') continue;
    
    // GRAFY TISK - anasys
        ;    b = read_ansys_strain_stress_plot("x2d_1I/", i, xy[k], j==2 ? 3 : j, -20.0, 20.0, true, false, true);
        ;        read_ansys_strain_stress_plot("x2d_1I/", i, xy[k], j==2 ? 3 : j, -20.0, 20.0, false);
        if (j<2) read_ansys_displ_plot        ("x2d_1I/", i, xy[k], j==2 ? 3 : j, -20.0, 20.0);
    
    // GRAFY TISK - mumech
    if (pert) {  print_graf(vtkAmndFile2d, "x2d_1I/", ".2d.pert",            i, xy[k],            j, -20.0, 20.0, 0.01, 'p', 0, DT_NUMERICAL, b); }
        if (totl) {  print_graf(vtkAmndFile2d, "x2d_1I/", ".2d"     ,            i, xy[k],            j, -20.0, 20.0, 0.01, 't', 0, DT_NUMERICAL, b); }
    
    if (pert) {  print_graf(vtkAmndFile3d, "x2d_1I/", ".3d.pert", i==2 ? 3 : i, xy[k], j==2 ? 3 : j, -20.0, 20.0, 0.01, 'p', 0, DT_ANALITICAL, b); }
    if (totl) {  print_graf(vtkAmndFile3d, "x2d_1I/", ".3d"     , i==2 ? 3 : i, xy[k], j==2 ? 3 : j, -20.0, 20.0, 0.01, 't', 0, DT_ANALITICAL, b); }
    
    // // COMPARE grafs
    // sprintf (name1, "%sgraf_lc%ld_%c_displc_%d%s.dat", "x2d_1I/", i, xy[k], j, ".2d");
    // sprintf (name2, "%sgraf_lc%ld_%c_displc_%d%s.dat", "x2d_1I/", i, xy[k], j, ".3d");
    // sprintf (name3, "%sgraf_lc%ld_%c_displc_%d%s.dat", "x2d_1I/", i, xy[k], j, ".ansys");
    // if (j < 2) { compare_grafs(name1, name2); compare_grafs(name2, name3); }
    // 
    // sprintf (name1, "%sgraf_lc%ld_%c_strain_%d%s.dat", "x2d_1I/", i, xy[k], j, ".2d");
    // sprintf (name2, "%sgraf_lc%ld_%c_strain_%d%s.dat", "x2d_1I/", i, xy[k], j, ".3d");
    // sprintf (name3, "%sgraf_lc%ld_%c_strain_%d%s.dat", "x2d_1I/", i, xy[k], j, ".ansys");
    // compare_grafs(name1, name2); compare_grafs(name2, name3);
    // 
    // sprintf (name1, "%sgraf_lc%ld_%c_stress_%d%s.dat", "x2d_1I/", i, xy[k], j, ".2d");
    // sprintf (name2, "%sgraf_lc%ld_%c_stress_%d%s.dat", "x2d_1I/", i, xy[k], j, ".3d");
    // sprintf (name3, "%sgraf_lc%ld_%c_stress_%d%s.dat", "x2d_1I/", i, xy[k], j, ".ansys");
    // compare_grafs(name1, name2); compare_grafs(name2, name3);
      }
    }
  }
}

void read_number_from_file (double &dx, const char *name)
{
  FILE *file = fopen(name, "r");
  fscanf (file,"%lf", &dx);
  fclose(file);
}

// Dve 2d elipticke inkluze na ose x, muze se menit jejich vzdalenost 'dx'.
// Inkluze jsou porad stejne, protoze v ANSYS reseni se da momentalne menit jen parametr 'dx' a hustota site.
void x2d_2I (void)
{
  Bod* b = NULL;
  char suff[100];
  char vtkAmndFile2d [] = "x2d_2I/amnd.2d.vtk";
  char vtkAmndFile3d [] = "x2d_2I/amnd.3d.vtk";
  
  // **************************************
  // ZDE URCIT, CO SE BUDE POCITAT/TISKNOUT
  // **************************************
  //double angle = 0.0;//                               -> u 2 inkluzi je uhel vzdy 0.0
  char xy[] = "X-";    // nebo treba "X-" "-Y"          -> zatim jen X
  int nlc = 1;    if (nlc !=1 && nlc != 1) errol;  //   -> zatim jen 1 lc, 3 nlc viz x2d_1I
  int ncomp = 3;  if (ncomp<1 || ncomp>3)  errol;
  
  double dx = 1.25;
  
  double Em =  1.0, num = 0.2;   // 1.0 0.2
  double Ei = 10.0, nui = 0.3;   // 1.0 0.4
  double rs2 [] = { 1.0, 0.0, 0.0, 0.0 };
  double rs3 [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  
  read_number_from_file (dx, "x2d_2I/dx.txt");
  
  bool m2d = true;
  bool m3d = false;
  bool ANSYS = true; //        void      standa add
  bool mode[] = { false, false, 1, false, 0,     0,   false, false };  int cm = 8;
  
  // OK SET - jeden vybrany set, ktery je ulozen, pouze pro kontrolu s novym resenim, ze se nic nezmenilo/neposralo
  m2d = true;
  m3d = false;
  ANSYS = false;
  mode[0] = mode[2] = mode[4] = mode[5] = true;
  mode[1] = mode[3] = mode[6] = mode[7] = false;
  
  
  // *******
  // COMPUTE
  // *******
  
  for (int m=0; m<cm; m++) {
    if (! mode[m])  continue;
    
    // vytvoreni 2d a 3d prikladu
    if(m2d) eqiv_gener_2x2dI_2d (vtkAmndFile2d, Em, num, rs2, dx, Ei, nui, 1.0, 0.5, DT_NUMERICAL , m);
    if(m3d) eqiv_gener_2x2dI_3d (vtkAmndFile3d, Em, num, rs3, dx, Ei, nui, 1.0, 0.5, DT_ANALITICAL, m);
    
    // GRAFY TISK - anasys
    if (ANSYS)  b = read_ansys_strain_stress_plot("x2d_2I/", 0, 'X', 0, -20.0, 20.0, true, false, true);
    if (ANSYS)      read_ansys_strain_stress_plot("x2d_2I/", 0, 'X', 0, -20.0, 20.0, false);
    if (ANSYS)      read_ansys_displ_plot        ("x2d_2I/", 0, 'X', 0, -20.0, 20.0);
    
    // GRAFY TISK - mumech
    if(m2d) { sprintf (suff, ".%02d.2d", m);  print_graf(vtkAmndFile2d, "x2d_2I/", suff, 0, 'X', 0, -20.0, 20.0, 0.01, 't', 0, DT_NUMERICAL , b); }
    if(m3d) { sprintf (suff, ".%02d.3d", m);  print_graf(vtkAmndFile3d, "x2d_2I/", suff, 0, 'X', 0, -20.0, 20.0, 0.01, 't', 0, DT_ANALITICAL, b); }
    
    delete [] b; b = NULL;
    
    //Comparison::cmp_mm_ansys_total_energy("x2d_2I/", vtkAmndFile2d, 0);
    
    
    // // debug
    // Problem *p = new Problem;
    // p->read_inclusions_plus_initialize_and_print (vtkAmndFile2d, NULL, DT_NUMERICAL);
    //
    // double coords[] = { 4.5, 0.0, 0.0};
    // double **strain  = AllocateArray2D(1,3);
    // double **stress  = AllocateArray2D(1,3);
    // double **stress2 = AllocateArray2D(1,3);
    // p->giveFieldsOfPoint(NULL,strain,stress,coords,'t',0,1,PFCM_FULL,-3, TVRN_THEORETICAL_FEEP);
    // double C[5];
    // p->give_inclusion(1)->give_StiffnessMatrixReduced(C);
    // MatrixOperations::giveTVproduct_3is3x3to5and3(stress2[0], C, strain[0]);
    //
    // coords[1] = 0.0; // kvuli break pointu
    //
    // DeleteArray2D(strain,1);
    // DeleteArray2D(stress,1);
    // DeleteArray2D(stress2,1);
    // // debug

  }
}


// 2d uloha, pravidelny grid kruhovych inkluzi "2n+1" x "2n+1".
void x2d_grid (void)
{
  Bod* b = NULL;
  char suff[100];
  char vtkAmndFile2d [] = "x2d_grid/amnd.2d.vtk";
  char vtkAmndFile3d [] = "x2d_grid/amnd.3d.vtk";
  
  // **************************************
  // ZDE URCIT, CO SE BUDE POCITAT/TISKNOUT
  // **************************************
  //double angle = 0.0;//                               -> u gridu inkluzi je uhel vzdy 0.0
  char xy[] = "XY";    // nebo treba "X-" "-Y"          -> zatim jen X
  int nlc = 1;    if (nlc !=1 && nlc != 1) errol;  //   -> zatim jen 1 lc, 3 nlc viz x2d_1I
  int ncomp = 3;  if (ncomp<0 || ncomp>3)  errol;  //   -> obvykle 3, at se tisknou vsechny slozky
  
  int nx = 2;
  int ny = 2;
  double dx = 1.5;
  
  double Em =  1.0, num = 0.2;
  double Ei = 10.0, nui = 0.3; // 10.0 0.3
  double rs2 [] = { 1.0, 0.0, 0.0, 0.0 };
  double rs3 [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

  //read_number_from_file (dx, "x2d_grid/dx.txt");
  //read_number_from_file (Ei, "x2d_grid/Ei.txt");
  
  bool m2d = true;
  bool m3d = false;
  bool ANSYS = true; //         void         standa addext
  bool mode[] = { false, false, true, false,  false, false, false, false  };  int cm = 8;
  
  // *******
  // COMPUTE
  // *******
  for (int m=0; m<cm; m++) {
    if (! mode[m])  continue;
    
    // vytvoreni 2d a 3d prikladu
    if(m2d) eqiv_gener_gx2dI_2d (vtkAmndFile2d, 1+2*nx, 1+2*ny, 2*dx, 2*dx, Em, num, nlc == 1 ? rs2 : NULL, Ei, nui, 1.0, 1.0, DT_NUMERICAL , m);
    if(m3d) eqiv_gener_gx2dI_3d (vtkAmndFile3d, 1+2*nx, 1+2*ny, 2*dx, 2*dx, Em, num, nlc == 1 ? rs3 : NULL, Ei, nui, 1.0, 1.0, DT_ANALITICAL, m);
    
    for (int i=0; i<nlc; i++) {
      if (0) {
        //double e2, e3;
        Comparison::cmp_mm_ansys_total_energy("x2d_grid/", vtkAmndFile2d, 0);
        //if (d3) e3 = Comparison::cmp_mm_ansys_total_energy("x2d_grid/ansys.IB", vtkAmndFile3d, "x2d_grid/ansys.IB.mesh.vtu", "x2d_grid/ansys.IB.rslt.vtu");

        //double vf = (PI*1.0*1.0) / (2*dx * 2*dx) * 100.0;
        //vf = (vf < 10) ? ((int)vf*10)/10 : vf;

        //if (d2) { FILE *file = fopen("x2d_grid/energy_error_2d.dat", "a");  fprintf(file, "%lf %lf %g\n", 2*dx, e2, vf);  fclose(file); }
        //if (d3) { FILE *file = fopen("x2d_grid/energy_error_3d.dat", "a");  fprintf(file, "%lf %lf %g\n", 2*dx, e3, vf);  fclose(file); }
      }

      if (1) {
        // loop over ncomp (number of components)
        for (int j=0; j<ncomp; j++) {
          // loop over axes
          for (int k=0; k<2; k++) {
            if (xy[k] == '-') continue;

            // GRAFY TISK - anasys
            if (ANSYS)    b = read_ansys_strain_stress_plot("x2d_grid/", i, xy[k], j==2 ? 3 : j, -20.0, 20.0, true, false, true);
            if (ANSYS)        read_ansys_strain_stress_plot("x2d_grid/", i, xy[k], j==2 ? 3 : j, -20.0, 20.0, false);
            if (ANSYS && j<2) read_ansys_displ_plot        ("x2d_grid/", i, xy[k], j==2 ? 3 : j, -20.0, 20.0);

            // GRAFY TISK - mumech
            if(m2d) { sprintf (suff, ".%02d.2d", m);  print_graf(vtkAmndFile2d, "x2d_grid/", suff,            i, xy[k],            j, -20.0, 20.0, 0.01, 't', 0, DT_NUMERICAL , b, true); }
            if(m3d) { sprintf (suff, ".%02d.3d", m);  print_graf(vtkAmndFile3d, "x2d_grid/", suff, i==2 ? 3 : i, xy[k], j==2 ? 3 : j, -20.0, 20.0, 0.01, 't', 0, DT_ANALITICAL, b, true); }

            // one
            //eqiv_gener_1x2dI_3d (vtkAmndFile3d, 1.0, 0.2, rs3, 0.0, 0.0, 10.0, 0.4, 1.0, 1.0, 0.0, DT_ANALITICAL);
            //print_graf(vtkAmndFile3d, "x2d_grid/", ".3dONE"     , 'x', 0, -dx, dx, 0.01, 't', 20, DT_ANALITICAL);
          }
        }
      }
    }
  }
}



// --------------------------------------------------------







// vzorovy 3d priklad o dvou inkluzich
void clanek_3d_vzor (void)
{
  Problem *p = new Problem;
  p->read_input_file("clanek/inhomogeneity.vtk");
  p->input_data_initialize_and_check_consistency();
  p->convert_to_equivalent_problem();
  p->print_equivalent_problem("clanek/equivalent.vtk");
  delete p;
  p = new Problem;
  p->read_input_file("clanek/equivalent.vtk");
  p->input_data_initialize_and_check_consistency();
  double coords[] = { 5.0,  2.0, 0.0};
  double **stress = AllocateArray2D(2,9);
  p->giveFieldsOfPoint(NULL,NULL,stress,coords,'p',0,2);
  p->printFieldsOnMeshVTK("clanek/results.vtk","clanek/mesh.vtk",'t',2,1);


  const double p1[] = {-4,-3,-3};
  const double p2[] = { 5, 4, 3};
  const long   nn[] = {100,100,100};
  p->printFieldsOnMeshGrid("clanek/results_grid.vtk", p1, p2, nn, 't', 0, 1, PFCM_FULL);
  //p->compute_and_print_fields_on_regular_mesh((const double[]){-4,-3,-3}, (const double[]){5,4,3}, (const long[]){20,20,20}, "clanek/fields_reg.vtk", 0, 1, PFCM_OPTIMIZED, 't');

  // int hID = p->give_new_homogenization_id();
  // p->homogenization(hID)->compute_fields_??? (TypeOfHomogenization);
  // const double **field = p->homogenization(hID)->give_stress();
  //p->print_equivalent_problem(vtkAmndFile);


  p->print_visualization("clanek/visualization.vtk",5);
  delete p;
  DeleteArray2D(stress,2);
}

void visualization_1I_2d (void)
{
  char vtkAmndFile [] = "print/amnd2d.vtk";
  char vtkMeshFileF [] = "print/meshF.2d.vtk";
  //char vtkMeshFileX [] = "print/meshX.2d.vtk";
  
  double rs2 [] = { 1.0, 0.0, 0.0, 0.0 };
  eqiv_gener_1x2dI_2d (vtkAmndFile, 1.0, 0.2, rs2, 0.0, 0.0, 10.0, 0.4, 1.0, 0.5, 0.0, DT_NUMERICAL);
  
  Problem *p = new Problem;
  p->read_inclusions_plus_initialize_and_print (vtkAmndFile, NULL, DT_NUMERICAL);
  
  p->print_visualization(vtkMeshFileF, 10, 2, true);
  //p->print_geom(vtkMeshFileX, 10, 2, false);
  delete p;
}

void visualization_1I_3d (void)
{
  char vtkAmndFile [] = "print/amnd.vtk";
  char vtkMeshFileF [] = "print/meshF.vtk";
  char vtkMeshFileX [] = "print/meshX.vtk";
  
  double rs [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  
  // ellipsoid                   Em   num        x    y    z   E    nu   sa1  sa2  sa3  ea1  ea2  ea3
  eqiv_gener_1x3dI (vtkAmndFile, 1.0, 0.4, rs, -2.0, 0.0, 0.0, 5.0, 0.4, 4.0, 0.9, 0.8, 0.0, 0.0, 0.0, DT_NUMERICAL);
  
  Problem *p = new Problem;
  p->read_inclusions_plus_initialize_and_print (vtkAmndFile, NULL, DT_ANALITICAL);
  
  p->print_visualization(vtkMeshFileF, 20, 3, true);
  p->print_visualization(vtkMeshFileX, 20, 3, false);
  delete p;
}

// The most simple way to visualize input geometry.
void visualization_197I_2d3d (void)
{
  Problem *p =new Problem;
  p->set_matrix_E_nu(1.0, 0.25);
  p->read_inclusions_plus_initialize_and_print("visualization/stackSmashing.geom.vtk", NULL, DT_ANALITICAL);
  p->print_visualization("visualization/stackSmashing.topo.vtk", 4, 3);
  delete p;
}



// homogenizace - 2d uloha, pravidelny grid kruhovych inkluzi "2n+1" x "2n+1".
void clanek_homog_grid (void)
{
  char vtkAmndFile2d [] = "clanek/grid.amnd.2d.vtk";

  int n = 3;
  double dx = 1.6; // 1.2 1.6

  // vytvoreni 2d prikladu
  eqiv_gener_gx2dI_2d (vtkAmndFile2d, 1+2*n, 1+2*n, 2*dx, 2*dx, 1.0, 0.2, NULL, 10.0, 0.3, 1.0, 1.0, DT_NUMERICAL , 5);

  Problem *p = new Problem;
  Homogenization *h;

  double sm[3*3];
  double p1[3];  p1[0] = -dx;  p1[1] = -dx;  p1[2] = 0.0;
  double p2[3];  p2[0] =  dx;  p2[1] =  dx;  p2[2] = 0.0;


  // ******
  //   2D
  // ******
  // PROBLEM
  p->read_input_file (vtkAmndFile2d);
  p->input_data_initialize_and_check_consistency();
  //p->print_visualization("clanek/grid.2d.vtk", 5);
  file_delete_if_exist(vtkAmndFile2d);

  // MESHES
  // muMECH solution
  long   nn[3];  nn[0] = nn[1] = 80;   nn[2] = 0;

  Mesh *meshMM = p->give_new_mesh();
  meshMM->generate_regular_mesh (p1, p2, nn);
  meshMM->compute_element_fields ('t', false, true, true, 0, 0, 0, PFCM_FULL); // false
  //meshMM->print_geometry_file_vtk("clanek/bounding_boxde.vtk",0,1); // debug


  // ANSYS solution
  Mesh *meshAS = p->give_new_mesh();
  meshAS->read_geometry_file_vtk ("clanek/homog_ansys/list.IB.rslt.lc0.vtu", 0, 0, false);
  meshAS->read_geometry_file_vtk ("clanek/homog_ansys/list.IB.rslt.lc1.vtu", 0, 1, true);
  meshAS->read_geometry_file_vtk ("clanek/homog_ansys/list.IB.rslt.lc2.vtu", 0, 2, true);

  FILE *file = fopen("clanek/grid.2d.txt", "w");

  // *** STANDARDNI HOMOGENIZACE ***
  p->give_inclusion(0)->give_StiffnessMatrixFull(sm);                                                                                                  fprintf (file, "%s:\t",   "Inclusion"       ); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  h = p->give_new_homogenization (HT_VoightBound); h->set_boundingBox (p1, p2);  h->find_inclusions_in_BB();  h->giveHomogenizedStiffnessMatrix (sm);  fprintf (file, "%s:\t",   h->giveClassName()); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  h = p->give_new_homogenization (HT_RegGrid);     ((RegGrid*)h)->set_mesh(meshMM);                           h->giveHomogenizedStiffnessMatrix (sm);  fprintf (file, "%s:\t",   h->giveClassName()); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  h = p->give_new_homogenization (HT_RegGrid);     ((RegGrid*)h)->set_mesh(meshAS);                           h->giveHomogenizedStiffnessMatrix (sm);  fprintf (file, "%s:\t",   h->giveClassName()); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  h = p->give_new_homogenization (HT_MoriTanaka);  h->set_boundingBox (p1, p2);  h->find_inclusions_in_BB();  h->giveHomogenizedStiffnessMatrix (sm);  fprintf (file, "%s:\t",   h->giveClassName()); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  h = p->give_new_homogenization (HT_Dilute);      h->set_boundingBox (p1, p2);  h->find_inclusions_in_BB();  h->giveHomogenizedStiffnessMatrix (sm);  fprintf (file, "%s:\t\t", h->giveClassName()); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  h = p->give_new_homogenization (HT_ReussBound);  h->set_boundingBox (p1, p2);  h->find_inclusions_in_BB();  h->giveHomogenizedStiffnessMatrix (sm);  fprintf (file, "%s:\t",   h->giveClassName()); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");
  p->matrix_giveFullStiffMatrix(sm);                                                                                                                   fprintf (file, "%s:\t\t", "Matrix"          ); for (int i=0; i<9; i++)  fprintf (file, " %15.8lf", sm[i]);   fprintf (file, "\n");

  double fv = h->giveTotalVolumeFractionOfInclusions();


  fprintf (file, "Fraction volume: %6.2lf\n", 100*fv);
  fclose(file);

  //
  delete p;
}






void hrana_inkluze (void)
{
  char vtkAmndFile [] = "hrana_inkluze/amnd.vtk";
  
  double rs [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  
  // ellipsoid                   Em   num        x    y    z   E    nu   sa1  sa2  sa3  ea1  ea2  ea3
  eqiv_gener_1x3dI (vtkAmndFile, 1.0, 0.4, rs, -2.0, 0.0, 0.0, 5.0, 0.4, 1.0, 0.9, 0.8, 0.0, 0.0, 0.0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.ellipsoid.numer.3d", 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.ellipsoid.anal.3d" , 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_ANALITICAL);
  
  // sphere                      Em   num        x    y    z   E    nu   sa1  sa2  sa3  ea1  ea2  ea3
  eqiv_gener_1x3dI (vtkAmndFile, 1.0, 0.4, rs, -2.0, 0.0, 0.0, 5.0, 0.4, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.sphere.numer.3d", 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.sphere.anal.3d" , 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_ANALITICAL);
  
  // oblate spheroid             Em   num        x    y    z   E    nu   sa1  sa2  sa3  ea1  ea2  ea3
  eqiv_gener_1x3dI (vtkAmndFile, 1.0, 0.4, rs, -2.0, 0.0, 0.0, 5.0, 0.4, 1.0, 1.0, 0.9, 0.0, 0.0, 0.0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.oblatesph.numer.3d", 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.oblatesph.anal.3d" , 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_ANALITICAL);
  
  // prolate spheroid            Em   num        x    y    z   E    nu   sa1  sa2  sa3  ea1  ea2  ea3
  eqiv_gener_1x3dI (vtkAmndFile, 1.0, 0.4, rs, -2.0, 0.0, 0.0, 5.0, 0.4, 1.0, 0.9, 0.9, 0.0, 0.0, 0.0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.prolatesph.numer.3d", 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_NUMERICAL);
  print_graf(vtkAmndFile, "hrana_inkluze/", ".1I.prolatesph.anal.3d" , 0, 'x', 0, -3.12, -2.88, 0.01, 'p', 0, DT_ANALITICAL);
}




// ****************************************   ENERGY   **********************************************



// **************************************************************************************************
// **************************************   KOMPLEXNI FUNKCE   **************************************
// **************************************************************************************************
void compute_and_print_reg_mesh (const char *vtkGeomFile, const char *vtkMeshFileOutPert, double x, int n, char pt, DiffTypes dt)
{
  Problem *p = new Problem;
  p->read_inclusions_plus_initialize_and_print (vtkGeomFile, NULL, dt);
  
  // regular mesh
  double p1[3] = {-x, -x, 0 };
  double p2[3] = { x,  x, 0 };
  long   nn[3] = { n,  n, 0 };
  
  p->printFieldsOnMeshGrid(vtkMeshFileOutPert, p1, p2, nn, pt, 0, 0, PFCM_FULL);
  
  delete p;
}



// **************************************************************************************************
// **********************************   GENEROVANI EQIV SOUBORU   ***********************************
// **************************************************************************************************

// generovani eqiv souboru s jednou 3d inkluzi
void eqiv_gener_1x3dI (const char *vtkAmndFile, double Em, double num, const double *rs,
               double x, double y, double z, double E, double nu,
               double a1, double a2, double a3, double ea1, double ea2, double ea3, DiffTypes dt)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set();
  p->set_dimension(3);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  //
  p->set_numberOfRemoteStrains(1);
  p->set_RemoteStrain(0, rs);
  
  // inclusions
  p->set_number_of_inclusions(1);
  // first inclusion
  p->set_inclusion_centroids(0, x, y, z);
  p->set_inclusion_E_nu(0, E, nu);
  p->set_inclusion_semiaxesDimensions(0, a1, a2, a3);
  p->set_inclusion_EulerAnglesDEG(0, ea1, ea2, ea3);
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  
  p->convert_to_equivalent_problem();
  
  p->print_equivalent_problem(vtkAmndFile);
  
  delete p;
}


// generovani eqiv souboru s jednou 2d inkluzi pomoci 2d mumechu
void eqiv_gener_1x2dI_2d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double y, double E, double nu, double a1, double a2, double ea, DiffTypes dt)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set();
  p->set_dimension(2);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  if (rs) { p->set_numberOfRemoteStrains(1);
            p->set_RemoteStrain(0, rs); }
  else      p->set_UnitRemoteStrains();
  
  // inclusions
  p->set_number_of_inclusions(1);
  // first inclusion
  p->set_inclusion_centroids(0, x, y);
  p->set_inclusion_E_nu(0, E, nu);
  p->set_inclusion_semiaxesDimensions(0, a1, a2);
  p->set_inclusion_EulerAnglesDEG(0, ea);
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  
  p->convert_to_equivalent_problem();
  p->print_visualization("vis.vtk",4);
  p->print_equivalent_problem(vtkAmndFile);
  
  delete p;
}

// generovani eqiv souboru s jednou 2d inkluzi pomoci 3d mumechu
void eqiv_gener_1x2dI_3d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double y, double E, double nu, double a1, double a2, double ea, DiffTypes dt)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set();
  p->set_dimension(3);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  if (rs) { p->set_numberOfRemoteStrains(1);
            p->set_RemoteStrain(0, rs); }
  else      p->set_UnitRemoteStrains();
  
  // inclusions
  p->set_number_of_inclusions(1);
  // first inclusion
  p->set_inclusion_centroids(0, x, y, 0.0);
  p->set_inclusion_E_nu(0, E, nu);
  p->set_inclusion_semiaxesDimensions(0, 100.0*a1, a1, a2);
  p->set_inclusion_EulerAnglesDEG(0, ea, 90.0, 90.0);
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  
  p->convert_to_equivalent_problem();
  
  p->print_equivalent_problem(vtkAmndFile);
  
  delete p;
}


void eqiv_gener_gx2dI_2d (const char *vtkAmndFile, int gx, int gy, double dx, double dy, double Em, double num, const double *rs, double E, double nu, double a1, double a2, DiffTypes dt, int mode)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set();
  p->set_dimension(2);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  //
  if (rs) { p->set_numberOfRemoteStrains(1);
            p->set_RemoteStrain(0, rs); }
  else      p->set_UnitRemoteStrains();
  
  // inclusions
  p->set_number_of_inclusions(gx*gy);
  long i, j, id;
  id = 0;
  for (i=0; i<gx; i++) {
    for (j=0; j<gy; j++) {
      p->set_inclusion_centroids(id, dx*(0.5*(gx-1)-i), dy*(0.5*(gy-1)-j));
      p->set_inclusion_E_nu(id, E, nu);
      p->set_inclusion_semiaxesDimensions(id, a1, a2);
      p->set_inclusion_EulerAnglesDEG(id, 0.0);
      id++;
    }
  }
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  
  p->set_SBA_optimized(false);
  
  modeset (p, mode);

  int na = 6;
  p->set_problem_size( -dx*(0.5*gx - 0.5) - na*a1, -dy*(0.5*gy - 0.5) - na*a2, 0.0,     dx*(0.5*gx - 0.5) + na*a1, dy*(0.5*gy - 0.5) + na*a2, 0.0,     0.05);

  p->convert_to_equivalent_problem();
  p->print_equivalent_problem(vtkAmndFile);
  p->print_visualization("vis2d.vtk", 4, 2);
  
  delete p;
}
//
void eqiv_gener_gx2dI_3d (const char *vtkAmndFile, int gx, int gy, double dx, double dy, double Em, double num, const double *rs, double E, double nu, double a1, double a2, DiffTypes dt, int mode)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set();
  p->set_dimension(3);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  //
  if (rs) { p->set_numberOfRemoteStrains(1);
            p->set_RemoteStrain(0, rs); }
  else      p->set_UnitRemoteStrains();
  
  // inclusions
  p->set_number_of_inclusions(gx*gy);
  long i, j, id;
  id = 0;
  for (i=0; i<gx; i++) {
    for (j=0; j<gy; j++) {
      p->set_inclusion_centroids(id, dx*(0.5*(gx-1)-i), dy*(0.5*(gy-1)-j), 0.0);
      p->set_inclusion_E_nu(id, E, nu);
      p->set_inclusion_semiaxesDimensions(id, 100*a1, a1, a2);
      p->set_inclusion_EulerAnglesDEG(id, 0.0, 90.0, 90.0);
      id++;
    }
  }
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  
  p->set_SBA_optimized(false);
  
  modeset (p, mode);
  
  p->convert_to_equivalent_problem();
  p->print_equivalent_problem(vtkAmndFile);
  //p->print_visualization("vis3d.vtk", 4, 3);
  
  delete p;
}

//
void eqiv_gener_2x2dI_3d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double E, double nu, double a1, double a2, DiffTypes dt, int mode)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set();
  p->set_dimension(3);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  //
  p->set_numberOfRemoteStrains(1);
  p->set_RemoteStrain(0, rs);
  
  // loop over all inclusions
  p->set_number_of_inclusions(2);
  // first inclusion
  p->set_inclusion_centroids(0, -x, 0.0, 0.0);
  p->set_inclusion_E_nu(0, E, nu);
  p->set_inclusion_semiaxesDimensions(0, 100*a1, a1, a2);
  p->set_inclusion_EulerAnglesDEG(0, 0.0, 90.0, 90.0);
  // second inclusion
  p->set_inclusion_centroids(1,  x, 0.0, 0.0);
  p->set_inclusion_E_nu(1, E, nu);
  p->set_inclusion_semiaxesDimensions(1, 100*a1, a1, a2);
  p->set_inclusion_EulerAnglesDEG(1, 0.0, 90.0, 90.0);
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT   ******************
  
  p->set_SBA_optimized(false);

  modeset (p, mode);

  p->convert_to_equivalent_problem();
  
  p->print_equivalent_problem(vtkAmndFile);
  
  delete p;
}

//
void eqiv_gener_2x2dI_2d (const char *vtkAmndFile, double Em, double num, const double *rs, double x, double E, double nu, double a1, double a2, DiffTypes dt, int mode)
{
  Problem *p = new Problem;
  
  // **************   DATA INPUT   ******************
  // START - DATA INPUT VIA FUNCTIONS
  // do not change the order of command
  
  // set problem
  p->set_data_set();
  p->set_dimension(2);
  p->set_diffType(dt);
  p->set_matrix_E_nu(Em, num);
  //
  p->set_numberOfRemoteStrains(1);
  p->set_RemoteStrain(0, rs);
  
  // loop over all inclusions
  p->set_number_of_inclusions(2);
  // first inclusion
  p->set_inclusion_centroids(0, -x, 0.0);
  p->set_inclusion_E_nu(0, E, nu);
  p->set_inclusion_semiaxesDimensions(0, a1, a2);
  p->set_inclusion_EulerAnglesDEG(0, 0.0);
  // second inclusion
  p->set_inclusion_centroids(1,  x, 0.0);
  p->set_inclusion_E_nu(1, E, nu);
  p->set_inclusion_semiaxesDimensions(1, a1, a2);
  p->set_inclusion_EulerAnglesDEG(1, 0.0);
  
  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT   ******************
  
  modeset (p, mode);
  p->set_SBA_optimized(false);
  p->convert_to_equivalent_problem();
  
  p->print_equivalent_problem(vtkAmndFile);
  
  delete p;
}


void modeset (Problem *p, int mode)
{
  if      (mode == 0) {  p->set_intFieldsShape(4);  p->set_SBAM(SBAT_ORIGINAL); }  // puvodni balancing
  else if (mode == 1) {  p->set_intFieldsShape(3);  p->set_SBAM(SBAT_ORIGINAL); }  // puvodni balancing s addExt
  else if (mode == 2) {  p->set_intFieldsShape(4);  p->set_SBAM(SBAT_VOID); }      // zadny balancing
  else if (mode == 3) {  p->set_intFieldsShape(3);  p->set_SBAM(SBAT_VOID); }      // zadny balancing s addExt
  else if (mode == 4) {  p->set_intFieldsShape(4);  p->set_SBAM(SBAT_STANDA); }
  else if (mode == 5) {  p->set_intFieldsShape(0);  p->set_SBAM(SBAT_STANDA); }
  else if (mode == 6) {  p->set_intFieldsShape(2);  p->set_SBAM(SBAT_MESHLESS); }
  else if (mode == 7) {  p->set_intFieldsShape(0);  p->set_SBAM(SBAT_MESHLESS); }
  else errol;
}


void print_graf (const char *vtkGeomFile, const char *path, const char *suff, int lc, char axe, int comp, double coordMin, double coordMax, double coordDelta, char pt_flag, double dshift, DiffTypes dt, Bod *b, bool odd)
{
  bool adapt = true;
  
  if (coordMin >= coordMax) errol;
  
  double *d = new double[9];  CleanVector(d, 9);
  double *e = new double[9];
  double *s = new double[9];
  
  Problem *p = new Problem;
  p->read_inclusions_plus_initialize_and_print (vtkGeomFile, NULL, dt);

  //
  int cind;
  if      (axe == 'X') cind = 0;
  else if (axe == 'Y') cind = 1;
  else if (axe == 'Z') cind = 2;
  else errol;
  
  double coords[] = { 0.0, 0.0, 0.0 };
  coords[cind] = coordMin;
  long n=0;
  
  char name[100];
  
  // 2D: predpoklad natvrdo, ze graf se tiskne na 2d uloze (pocitana muze byt i 3d mumechem)
  bool D =  comp >1 ? false : true;
  
  // 2D: predpoklad natvrdo, ze graf se tiskne na 2d uloze (pocitana muze byt i 3d mumechem)
  FILE * fileD;  sprintf(name, "%sgraf_lc%d_%c_displc_%d%s.dat", path, lc == 3 ? 2 : lc, axe,                 comp, suff);  if (D) OpenFile(fileD, name, _WRITE_);
  FILE * fileE;  sprintf(name, "%sgraf_lc%d_%c_strain_%d%s.dat", path, lc == 3 ? 2 : lc, axe, comp == 3 ? 2 : comp, suff);         OpenFile(fileE, name, _WRITE_);
  FILE * fileS;  sprintf(name, "%sgraf_lc%d_%c_stress_%d%s.dat", path, lc == 3 ? 2 : lc, axe, comp == 3 ? 2 : comp, suff);         OpenFile(fileS, name, _WRITE_);
  
  if (b) {
    long i= 0;
    bool out = ! odd;
    do {
      if (b[i].coord[cind] < coordMin  ||  coordMax < b[i].coord[cind])  continue;
      for (int aa=0; aa<2; aa++) {
        p->giveFieldsOfPointOneRS(d, e, s, b[i].coord, pt_flag, lc, PFCM_FULL, out ? -1 : -2, STRN_VOIGT_FEEP);
    
    if (aa==0 && D) fprintf (fileD, "%lf\t%lf\n", b[i].coord[cind], d[comp] + dshift);
    ;               fprintf (fileE, "%lf\t%lf\n", b[i].coord[cind], e[comp]);
    ;               fprintf (fileS, "%lf\t%lf\n", b[i].coord[cind], s[comp]);
    
        if (b[i].mat[0] != 0 && b[i].mat[1] != 0  &&  b[i].mat[0] != b[i].mat[1]) { // pokud bod obsahuje dva materialy
          if (b[i].mat[0] > 9 && b[i].mat[1] > 9) break;  // oba materialy jsou matrice, tj. maji cisla vetsi nez 9
      if(aa==0) out = out ? false : true;
    }
    else {
      break;
    }
      }
    }
    while (b[i++].last != true);
  }
  else {
    double toler = 1.e-3;
    if (3*toler > coordDelta)  toler = 0.2 * coordDelta;  
    long pic_old, pic, picreq; // parent incl id
    pic_old = -1;
    picreq = -2;
    bool doublep = false;
    //double c1;
    
    while (coords[cind] - coordMax < 0.1*coordDelta ) {
      if (!doublep)
    pic = p->give_inclusion_with_point_inside (coords, ((pic_old == -1) ? 1.0 : -1.0) * toler);
      else
    doublep = false;
      
      picreq = pic;
      // prekrocil jsem hranici inkluze ...
      if (pic != pic_old) {
    // ... smerem dovnitr inkluze
    if (pic_old == -1) {
          if ( p->give_inclusion(pic    )->point_is_inside (coords, -1 * toler) == false) {
        picreq = pic_old;
        doublep = true;
      }
    }
    // ... smerem ven z inkluze
    else {
          if ( p->give_inclusion(pic_old)->point_is_inside (coords,      toler) == true ) {
        picreq = pic_old;
        doublep = true;
      }
    }
      }
      pic_old = pic;
      
      p->giveFieldsOfPointOneRS(d, e, s, coords, pt_flag, lc, PFCM_FULL, picreq, STRN_VOIGT_FEEP);

      if (D) fprintf (fileD, "%lf\t%lf\n", coords[cind], d[comp] + dshift);
      ;      fprintf (fileE, "%lf\t%lf\n", coords[cind], e[comp]);
      ;      fprintf (fileS, "%lf\t%lf\n", coords[cind], s[comp]);
      
      if (!doublep) {
    if (!adapt) coords[cind] = coordMin + (++n)*coordDelta;
    else        coords[cind] += coordDelta;
      }
      //c1 = coords[0];
      
      picreq = -2;
    }
  }
  
  if (D) CloseFile(fileD);
  CloseFile(fileE);
  CloseFile(fileS);
  
  delete p;
  
  delete [] d;
  delete [] e;
  delete [] s;
}

// **************************************************************************************************
// ************************************   ANSYS POSTPROCESSING   ************************************
// **************************************************************************************************

bool compare_x(Bod const& first, Bod const& second) { return first.coord[0] < second.coord[0]; }
bool compare_y(Bod const& first, Bod const& second) { return first.coord[1] < second.coord[1]; }

Bod* read_ansys_strain_stress_plot (const char *path, int lc, char axe, int comp, double coordMin, double coordMax, bool strn, bool half, bool bod)
{
  char name[1000];
  long i, j, nNodes;
  Bod * Res;
  FILE *file;
  
  // nacti souradnice
  sprintf (name, "%sANSYS/rslts_lc%d/list.%c.nodes.txt", path, lc, axe);
  file = _openFileN("r", "ansys coords file", name);
  
  FP_skip_line(file);
  FP_skip_expected_string(file, " *GET  NNODE     FROM  NODE  ITEM=COUN       VALUE=", true);
  fscanf(file, "%ld", &nNodes);
  FP_skip_line(file, 3);
  
  Res = new Bod[nNodes];
  memset(Res, 0, nNodes*sizeof(Bod));

  for (i = 0; i < nNodes; i++) {
    if (i % 50 == 0)
      FP_skip_line(file, 11);

    fscanf(file, "%d", &Res[i].ID);
    fscanf(file, "%lf %lf %lf", &Res[i].coord[0], &Res[i].coord[1], &Res[i].coord[2]);
  }
  
  fclose(file);
  
  
  // nacti vysledky
  sprintf (name, "%sANSYS/rslts_lc%d/list.%c.%s.txt", path, lc, axe, strn ? "strain" : "stress");
  file = _openFileN("r", "ansys result file", name);
  
  int current_mat;
  
  FP_skip_to_line_starting_with(file, "NODAL RESULTS ARE FOR MATERIAL", true);
  fscanf(file, "%d", &current_mat);
  FP_skip_line(file, 5);
  
  int current_node;
  char c;
  
  while (true) { // podle poctu uzlu
    c = fgetc(file);
    if (c == ' ') {
      fscanf(file, "%d", &current_node);  // kontroluju jestli radek zacina cislem
      for (i = 0; i < nNodes; i++) { // cyklus pres vsechny body
        if (current_node == Res[i].ID) { // kdyz se shoduje ID
      int pos = 0;
          if (Res[i].mat[0] != 0) {
            if (Res[i].mat[1] != 0) errol;
        pos = 1;
      }
      
          Res[i].mat[pos] = current_mat;
          for (j = 0; j < 6; j++)
            fscanf(file, "%lf", &Res[i].rslt[pos][j]);
      break;
    }
      }
      if (i == nNodes)
    errol;
      
      FP_skip_line(file);
    }
    else if (c == '1') {
      FP_skip_to_line_starting_with(file, "NODAL RESULTS ARE FOR MATERIAL", true);
      fscanf(file, "%d", &current_mat);
      FP_skip_line(file, 5);
    }
    else
      break;
  }
  
  fclose(file);
  
  int coord_index;
  if      (axe == 'X') { coord_index = 0;  std::sort(Res, Res + nNodes, &compare_x); }
  else if (axe == 'Y') { coord_index = 1;  std::sort(Res, Res + nNodes, &compare_y); }
  else errol;
  
  // tiskni graf
  sprintf (name, "%sgraf_lc%d_%c_%s_%d.ansys.dat", path, lc, axe, strn ? "strain" : "stress", comp == 3 ? 2 : comp);
  file = _openFileN("w", "plot file", name);
  
  int cind;
  if      (axe == 'X') cind = 0;
  else if (axe == 'Y') cind = 1;
  else if (axe == 'Z') cind = 2;
  else errol;
  
  for (i = 0; i < nNodes; i++) {
    if (Res[i].coord[cind] < coordMin  ||  coordMax < Res[i].coord[cind])  continue;
    
    if (Res[i].mat[0] != 0 && Res[i].mat[1] != 0) { // pokud bod obsahuje dva materialy
      
      // kdyz predchozi bod.material na prvni pozici neni nula vytisknu nejdriv pro prvni material a pak pro druhy
      if (Res[i - 1].mat[0] == Res[i].mat[0]) {
        fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].rslt[0][comp]);
        fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].rslt[1][comp]);
      }
      // kdyz predchozi bod.material na druhe pozici neni nula vytisknu nejdriv pro druhy material a pak pro prvni
      else if (Res[i - 1].mat[0] == Res[i].mat[1]) {
        fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].rslt[1][comp]);
        fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].rslt[0][comp]);
      }
      else errol;
    }
    else { // kdyz bod neobsahuje dva materialy
      if (Res[i].mat[0]) // kdyz je to material 1
        fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].rslt[0][comp]);
      else if (Res[i].mat[1]) // kdyz je to material 2
        fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].rslt[1][comp]);
      else
    errol;
    }
  }
  
  fclose(file);
  
  Bod *b = NULL;
  int nn;
  if (bod) {
    if (half) nn = nNodes-1;
    else      nn = 0;

    b = new Bod[nNodes + nn];
    memcpy (b+nn, Res, sizeof(Bod));
    
    for (i = 1; i < nNodes; i++) {
      ;           memcpy (b+nn+i, Res+i, sizeof(Bod));
      if (half) { memcpy (b+nn-i, Res+i, sizeof(Bod));
                  b[nn-i].coord[coord_index] *= -1; }
    }
    
    b[nn + nNodes-1].last = true;
  }
  
  delete [] Res;
  
  return b;
}

void read_ansys_displ_plot (const char *path, int lc, char axe, int comp, double coordMin, double coordMax)
{
  char name[1000];
  long i, j, nNodes;
  Bod * Res;
  FILE *file;
  
  // nacti souradnice
  sprintf (name, "%sANSYS/rslts_lc%d/list.%c.nodes.txt", path, lc, axe);
  file = _openFileN("r", "ansys coords file", name);
  
  FP_skip_line(file);
  FP_skip_expected_string(file, " *GET  NNODE     FROM  NODE  ITEM=COUN       VALUE=", true);
  fscanf(file, "%ld", &nNodes);
  FP_skip_line(file, 3);
  
  Res = new Bod[nNodes];
  memset(Res, 0, nNodes*sizeof(Bod));

  for (i = 0; i < nNodes; i++) {
    if (i % 50 == 0)
      FP_skip_line(file, 11);

    fscanf(file, "%d", &Res[i].ID);
    fscanf(file, "%lf %lf %lf", &Res[i].coord[0], &Res[i].coord[1], &Res[i].coord[2]);
  }
  
  fclose(file);
  
  
  // nacti vysledky
  sprintf (name, "%sANSYS/rslts_lc%d/list.%c.%s.txt", path, lc, axe, "displ");
  file = _openFileN("r", "ansys result file", name);
  
  FP_skip_to_line_starting_with(file, "NODE      UX          UY          UZ          USUM  ", true);
  FP_skip_line(file);
  
  int current_node;
  char c;
  
  while (true) { // podle poctu uzlu
    c = fgetc(file);
    if (c == ' ') {
      fscanf(file, "%d", &current_node);  // kontroluju jestli radek zacina cislem
      for (i = 0; i < nNodes; i++) { // cyklus pres vsechny body
        if (current_node == Res[i].ID) { // kdyz se shoduje ID
      
      for (j = 0; j < 4; j++)
            fscanf(file, "%lf", &Res[i].disp[j]);
      break;
    }
      }
      if (i == nNodes)
    errol;
      
      FP_skip_line(file);
    }
    else if (c == '1') {
      FP_skip_to_line_starting_with(file, "NODE      UX          UY          UZ          USUM  ", true);
      FP_skip_line(file);
    }
    else
      break;
  }
  
  fclose(file);
  
  int coord_index;
  if      (axe == 'X') { coord_index = 0;  std::sort(Res, Res + nNodes, &compare_x); }
  else if (axe == 'Y') { coord_index = 1;  std::sort(Res, Res + nNodes, &compare_y); }
  else errol;
  
  // tiskni graf
  sprintf (name, "%sgraf_lc%d_%c_%s_%d.ansys.dat", path, lc, axe, "displc", comp);
  file = _openFileN("w", "plot file", name);
  
  int cind;
  if      (axe == 'X') cind = 0;
  else if (axe == 'Y') cind = 1;
  else if (axe == 'Z') cind = 2;
  else errol;
  
  for (i = 0; i < nNodes; i++) {
    if (Res[i].coord[cind] < coordMin  ||  coordMax < Res[i].coord[cind])  continue;
    
    fprintf(file, "%lf\t%lf\n", Res[i].coord[coord_index], Res[i].disp[comp]);
  }
  
  fclose(file);
  
  delete [] Res;
  
  return;
}


void speedtest2Dx3D(void){
  int i, j;
  double ax[3];
  ax[0] = 50.;
  ax[1] = 1.;
  ax[2] = 0.5;
  double r1 = ax[1] * 1.01;
  double r2 = ax[1] * 4.5;
  double dr = ax[1] * 0.005;
  int boduu = round((r2 - r1) / dr + 1);
  int bodu = 360 * boduu;
  double **points2D;
  double **points3D;
  points2D = new double*[bodu];
  points3D = new double*[bodu];
  for (i = 0; i < bodu; i++)points2D[i] = new double[3];
  for (i = 0; i < bodu; i++)points3D[i] = new double[3];
  for (i = 0; i < 360; i++)for (j = 0; j < boduu; j++){
    points2D[i * boduu + j][0] = (r1 + dr*j)*cos(PI*i / 180);
    points2D[i * boduu + j][1] = (r1 + dr*j)*sin(PI*i / 180);
    points2D[i * boduu + j][2] = 0.;
    points3D[i * boduu + j][0] = 0.;
    points3D[i * boduu + j][1] = (r1 + dr*j)*cos(PI*i / 180);
    points3D[i * boduu + j][2] = (r1 + dr*j)*sin(PI*i / 180);
  }
  /*double a = ax[1]*3.;
  double da = ax[1]*0.01;
  int bodu = round((2 * a / da + 1)*(2 * a / da + 1));
  int bodusqrt = round(2 * a / da + 1);
  double **points2D;
  double **points3D;
  points2D = new double*[bodu];
  points3D = new double*[bodu];
  for (i = 0; i < bodu; i++)points2D[i] = new double[3];
  for (i = 0; i < bodu; i++)points3D[i] = new double[3];
  for (i = 0; i < bodusqrt; i++)for (j = 0; j < bodusqrt; j++){
  points2D[i*bodusqrt + j][0] = da*(i - a / da);
  points2D[i*bodusqrt + j][1] = da*(j - a / da);
  points2D[i*bodusqrt + j][2] = 0.;
  points3D[i*bodusqrt + j][0] = 0.;
  points3D[i*bodusqrt + j][1] = da*(i - a / da);
  points3D[i*bodusqrt + j][2] = da*(j - a / da);
  }*/
  Problem *analFce2D, *analFce3D;

  int nlc = 1;

  double **strain_2D;
  strain_2D = new double*[nlc];
  for (i = 0; i < nlc; i++)strain_2D[i] = new double[4];
  double **stress_2D;
  stress_2D = new double*[nlc];
  for (i = 0; i < nlc; i++)stress_2D[i] = new double[4];
  double **disp_2D;
  disp_2D = new double*[nlc];
  for (i = 0; i < nlc; i++)disp_2D[i] = new double[2];

  double **strain_3D;
  strain_3D = new double*[nlc];
  for (i = 0; i < nlc; i++)strain_3D[i] = new double[9];
  double **stress_3D;
  stress_3D = new double*[nlc];
  for (i = 0; i < nlc; i++)stress_3D[i] = new double[9];
  double **disp_3D;
  disp_3D = new double*[nlc];
  for (i = 0; i < nlc; i++)disp_3D[i] = new double[3];
  //int _2D[] = { 0, 3, 1 };
  //int _3D[] = { 4, 8, 5 };
  clock_t ta, tb;
  double times[4];
  double rs2D[] = { 1.0, 0.0, 0.0, 0.0 };
  double rs3D[] = { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
  //rozjezd

  //2D
  printf("Starting the measurement...\n");
  analFce2D = new Problem();

  analFce2D->set_dimension(2);
  analFce2D->set_number_of_inclusions(1);
  analFce2D->set_inclusion_shape(0, IS_ELLIPSE);
  analFce2D->set_inclusion_centroids(0, 0, 0);
  analFce2D->set_inclusion_E_nu(0, 2.4, 0.3);
  analFce2D->set_inclusion_EulerAnglesDEG(0, 30);
  analFce2D->set_inclusion_semiaxesDimensions(0, ax[1], ax[2]);
  //analFce2D->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 1., 0.5, 0, 30, 0, 0);
  analFce2D->set_numberOfRemoteStrains(1);
  analFce2D->set_RemoteStrain(0, rs2D);
  analFce2D->set_data_set();

  analFce2D->set_SBA_optimized(true);
  analFce2D->set_diffType(DT_NUMERICAL);
  analFce2D->input_data_initialize_and_check_consistency();
  analFce2D->convert_to_equivalent_problem();
  for (i = 0; i < bodu; i++)analFce2D->giveFieldsOfPoint(disp_2D, strain_2D, stress_2D, points2D[i], 'p', 0, nlc, PFCM_OPTIMIZED);
  delete analFce2D;

  //zacatek mereni
  
  //2D
  printf("Measuring 2D numerical diff...\n");
  analFce2D = new Problem();
  
  analFce2D->set_dimension(2);
  analFce2D->set_number_of_inclusions(1);
  analFce2D->set_inclusion_shape(0, IS_ELLIPSE);
  analFce2D->set_inclusion_centroids(0, 0, 0);
  analFce2D->set_inclusion_E_nu(0, 2.4, 0.3);
  analFce2D->set_inclusion_EulerAnglesDEG(0, 30);
  analFce2D->set_inclusion_semiaxesDimensions(0, ax[1], ax[2]);
  //analFce2D->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 1., 0.5, 0, 30, 0, 0);
  analFce2D->set_numberOfRemoteStrains(1);
  analFce2D->set_RemoteStrain(0, rs2D);
  analFce2D->set_data_set();

  analFce2D->set_SBA_optimized(true);
  analFce2D->set_diffType(DT_NUMERICAL);
  analFce2D->input_data_initialize_and_check_consistency();
  analFce2D->convert_to_equivalent_problem();
  ta = clock();
  for (i = 0; i < bodu; i++)analFce2D->giveFieldsOfPoint(disp_2D, strain_2D, stress_2D, points2D[i], 'p', 0, nlc, PFCM_OPTIMIZED);
  tb = clock();
  delete analFce2D;
  times[0] = (tb - ta) / (double)CLOCKS_PER_SEC;
  printf("It took %.2lf seconds.\n", times[0]);

  //3D
  printf("Measuring 3D numerical diff...\n");
  analFce3D = new Problem();

  analFce3D->set_dimension(3);
  analFce3D->set_number_of_inclusions(1);
  analFce3D->set_inclusion_shape(0, IS_ELLIPSOID);
  analFce3D->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, ax[0], ax[1], ax[2], 0, 30, 0);
  analFce3D->set_numberOfRemoteStrains(1);
  analFce3D->set_RemoteStrain(0, rs3D);
  analFce3D->set_data_set();

  analFce3D->set_SBA_optimized(true);
  analFce3D->set_diffType(DT_NUMERICAL);
  analFce3D->input_data_initialize_and_check_consistency();
  analFce3D->convert_to_equivalent_problem();
  ta = clock();
  for (i = 0; i < bodu; i++)analFce3D->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points3D[i], 'p', 0, nlc, PFCM_OPTIMIZED);
  tb = clock();
  delete analFce3D;
  times[1] = (tb - ta) / (double)CLOCKS_PER_SEC;
  printf("It took %.2lf seconds.\n", times[1]);

  //3D
  printf("Measuring 3D analytical diff...\n");
  analFce3D = new Problem();

  analFce3D->set_dimension(3);
  analFce3D->set_number_of_inclusions(1);
  analFce3D->set_inclusion_shape(0, IS_ELLIPSOID);
  analFce3D->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, ax[0], ax[1], ax[2], 0, 30, 0);
  analFce3D->set_numberOfRemoteStrains(1);
  analFce3D->set_RemoteStrain(0, rs3D);
  analFce3D->set_data_set();

  analFce3D->set_SBA_optimized(true);
  analFce3D->set_diffType(DT_ANALITICAL);
  analFce3D->input_data_initialize_and_check_consistency();
  analFce3D->convert_to_equivalent_problem();
  ta = clock();
  for (i = 0; i < bodu; i++)analFce3D->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points3D[i], 'p', 0, nlc, PFCM_OPTIMIZED);
  tb = clock();
  delete analFce3D;
  times[2] = (tb - ta) / (double)CLOCKS_PER_SEC;
  printf("It took %.2lf seconds.\n", times[2]);

  //3D
  printf("Measuring 3D complex diff...\n");
  analFce3D = new Problem();

  analFce3D->set_dimension(3);
  analFce3D->set_number_of_inclusions(1);
  analFce3D->set_inclusion_shape(0, IS_ELLIPSOID);
  analFce3D->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, ax[0], ax[1], ax[2], 0, 30, 0);
  analFce3D->set_numberOfRemoteStrains(1);
  analFce3D->set_RemoteStrain(0, rs3D);
  analFce3D->set_data_set();

  analFce3D->set_SBA_optimized(true);
  analFce3D->set_diffType(DT_COMPLEX);
  analFce3D->input_data_initialize_and_check_consistency();
  analFce3D->convert_to_equivalent_problem();
  ta = clock();
  for (i = 0; i < bodu; i++)analFce3D->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points3D[i], 'p', 0, nlc, PFCM_OPTIMIZED);
  tb = clock();
  delete analFce3D;
  times[3] = (tb - ta) / (double)CLOCKS_PER_SEC;
  printf("It took %.2lf seconds.\n", times[3]);

  //konec mereni -> prehled vysledku
  double ref = times[2];
  printf("\n2D ND: %.1lf\n3D ND: %.1lf\n3D AD: %.1lf\n3D CD: %.1lf\n", times[0] / ref * 100, times[1] / ref * 100, times[2] / ref * 100, times[3] / ref * 100);
}
void createNodes(double a, double b, double c, double da, double db, double dc,double **&points,int &n) {
    int i,ia,ib,ic;
    int na, nb, nc;
    na = round((2.*a / da + 1));
    nb = round((2.*b / db + 1));
    nc = round((2.*c / dc + 1));
    n = na*nb*nc;
    points = new double*[n];
    for (i = 0; i < n; i++)points[i] = new double[3];
    int id = 0;
    for (ia = 0; ia < na; ia++)for (ib = 0; ib < nb; ib++)for (ic = 0; ic < nc; ic++) {
        points[id][0] = -a + ia*da;
        points[id][1] = -b + ib*db;
        points[id][2] = -c + ic*dc;
        id++;
    }
}
void test3D(void) {

    

    int i, j, bodu;
    double **points;
    Problem *analFce;
    int nlc = 1;
    double a = 10., b = 10., c = 10., da = .5, db = .5, dc = .5;
    // enum DiffTypes { DT_Void, DT_ANALITICAL, DT_NUMERICAL, DT_COMPLEX };
    DiffTypes difftype = DT_NUMERICAL;
    
    

    double **strain_3D;
    strain_3D = new double*[nlc];
    for (i = 0; i < nlc; i++)strain_3D[i] = new double[9];
    double **stress_3D;
    stress_3D = new double*[nlc];
    for (i = 0; i < nlc; i++)stress_3D[i] = new double[9];
    double **disp_3D;
    disp_3D = new double*[nlc];
    for (i = 0; i < nlc; i++)disp_3D[i] = new double[3];
    clock_t ta, tb;
        //double times[4];
    double rs3D[] = { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
    
    if (difftype == DT_NUMERICAL)printf("NUMERICAL DIFF\n");
    if (difftype == DT_ANALITICAL)printf("ANALYTICAL DIFF\n");
    if (difftype == DT_COMPLEX)printf("COMPLEX DIFF\n");

    //3D Ellipsoid
    createNodes(a,b,c,da,db,dc, points, bodu);
    printf("Computing 3D Ellipsoid...\n");
    analFce = new Problem();
    analFce->set_dimension(3);
    analFce->set_number_of_inclusions(1);
    analFce->set_inclusion_shape(0, IS_ELLIPSOID);
    analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 8., 5., 4., 20., 30., 50.);
    analFce->set_numberOfRemoteStrains(1);
    analFce->set_RemoteStrain(0, rs3D);
    analFce->set_data_set();
        analFce->set_SBA_optimized(true);
    analFce->set_diffType(difftype);
    analFce->input_data_initialize_and_check_consistency();
    analFce->convert_to_equivalent_problem();
    for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
    delete analFce;
    for (i = 0; i < bodu; i++)delete []points[i];
    delete[]points;

    //3D Sphere
    createNodes(a, b, c, da, db, dc, points, bodu);
    printf("Computing 3D Sphere...\n");
    analFce = new Problem();
    analFce->set_dimension(3);
    analFce->set_number_of_inclusions(1);
    analFce->set_inclusion_shape(0, IS_SPHERE);
    analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 5., 5., 5., 20., 30., 50.);
    analFce->set_numberOfRemoteStrains(1);
    analFce->set_RemoteStrain(0, rs3D);
    analFce->set_data_set();
        analFce->set_SBA_optimized(true);
    analFce->set_diffType(difftype);
    analFce->input_data_initialize_and_check_consistency();
    analFce->convert_to_equivalent_problem();
    for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
    delete analFce;
    for (i = 0; i < bodu; i++)delete[]points[i];
    delete[]points;

        //createNodes(a, b, c, da, db, dc, points, bodu);
        //printf("Computing 3D Ellyptic cylinder...\n");
        //analFce = new Problem();
        //analFce->set_dimension(3);
        //analFce->set_number_of_inclusions(1);
        //analFce->set_inclusion_shape(0, IS_ELLIPTIC_CYLINDER);
        //analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, INFTY, 5., 2., 20., 30., 50.);
        //analFce->set_numberOfRemoteStrains(1);
        //analFce->set_RemoteStrain(0, rs3D);
        //analFce->set_data_set();
        //analFce->set_SBA_optimized(true);
        //analFce->set_diffType(difftype);
        //analFce->input_data_initialize_and_check_consistency();
        //analFce->convert_to_equivalent_problem();
        //for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
        //delete analFce;
        //for (i = 0; i < bodu; i++)delete[]points[i];
        //delete[]points;

        //createNodes(a, b, c, da, db, dc, points, bodu);
        //printf("Computing 3D Cylinder...\n");
        //analFce = new Problem();
        //analFce->set_dimension(3);
        //analFce->set_number_of_inclusions(1);
        //analFce->set_inclusion_shape(0, IS_CYLINDER);
        //analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, INFTY, 5., 5., 20., 30., 50.);
        //analFce->set_numberOfRemoteStrains(1);
        //analFce->set_RemoteStrain(0, rs3D);
        //analFce->set_data_set();
        //analFce->set_SBA_optimized(true);
        //analFce->set_diffType(difftype);
        //analFce->input_data_initialize_and_check_consistency();
        //analFce->convert_to_equivalent_problem();
        //for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
        //delete analFce;
        //for (i = 0; i < bodu; i++)delete[]points[i];
        //delete[]points;

    //3D Oblate spheroid
    createNodes(a, b, c, da, db, dc, points, bodu);
    printf("Computing 3D Oblate spheroid...\n");
    analFce = new Problem();
    analFce->set_dimension(3);
    analFce->set_number_of_inclusions(1);
    analFce->set_inclusion_shape(0, IS_OBLATE_SPHEROID);
    analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 5., 5., 2., 20., 30., 50.);
    analFce->set_numberOfRemoteStrains(1);
    analFce->set_RemoteStrain(0, rs3D);
    analFce->set_data_set();
        analFce->set_SBA_optimized(true);
    analFce->set_diffType(difftype);
    analFce->input_data_initialize_and_check_consistency();
    analFce->convert_to_equivalent_problem();
    for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
    delete analFce;
    for (i = 0; i < bodu; i++)delete[]points[i];
    delete[]points;

    //3D Prolate spheroid
    createNodes(a, b, c, da, db, dc, points, bodu);
    printf("Computing 3D Prolate spheroid...\n");
    analFce = new Problem();
    analFce->set_dimension(3);
    analFce->set_number_of_inclusions(1);
    analFce->set_inclusion_shape(0, IS_PROLATE_SPHEROID);
    analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 8., 5., 5., 20., 30., 50.);
    analFce->set_numberOfRemoteStrains(1);
    analFce->set_RemoteStrain(0, rs3D);
    analFce->set_data_set();
        analFce->set_SBA_optimized(true);
    analFce->set_diffType(difftype);
    analFce->input_data_initialize_and_check_consistency();
    analFce->convert_to_equivalent_problem();
    for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
    delete analFce;
    for (i = 0; i < bodu; i++)delete[]points[i];
    delete[]points;
    
        //createNodes(a, b, c, da, db, dc, points, bodu);
        //printf("Computing 3D Flat ellipsoid...\n");
        //analFce = new Problem();
        //analFce->set_dimension(3);
        //analFce->set_number_of_inclusions(1);
        //analFce->set_inclusion_shape(0, IS_FLAT_ELLIPSOID);
        //analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 5., 3., 0., 20., 30., 50.);
        //analFce->set_numberOfRemoteStrains(1);
        //analFce->set_RemoteStrain(0, rs3D);
        //analFce->set_data_set();
        //analFce->set_SBA_optimized(true);
        //analFce->set_diffType(difftype);
        //analFce->input_data_initialize_and_check_consistency();
        //analFce->convert_to_equivalent_problem();
        //for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
        //delete analFce;
        //for (i = 0; i < bodu; i++)delete[]points[i];
        //delete[]points;

        //createNodes(a, b, c, da, db, dc, points, bodu);
        //printf("Computing 3D Penny...\n");
        //analFce = new Problem();
        //analFce->set_dimension(3);
        //analFce->set_number_of_inclusions(1);
        //analFce->set_inclusion_shape(0, IS_PENNY);
        //analFce->set_inclusion_all(0, 0, 0, 0, 2.4, 0.3, 5., 5., 0., 20., 30., 50.);
        //analFce->set_numberOfRemoteStrains(1);
        //analFce->set_RemoteStrain(0, rs3D);
        //analFce->set_data_set();
        //analFce->set_SBA_optimized(true);
        //analFce->set_diffType(difftype);
        //analFce->input_data_initialize_and_check_consistency();
        //analFce->convert_to_equivalent_problem();
        //for (i = 0; i < bodu; i++)analFce->giveFieldsOfPoint(disp_3D, strain_3D, stress_3D, points[i], 'p', 0, nlc, PFCM_OPTIMIZED);
        //delete analFce;
        //for (i = 0; i < bodu; i++)delete[]points[i];
        //delete[]points;

    /*
    IS_ELLIPSOID = 1, IS_SPHERE = 2,
    IS_OBLATE_SPHEROID = 8, IS_PROLATE_SPHEROID = 9,
    IS_ELLIPTIC_CYLINDER = 3, IS_CYLINDER = 4,
    IS_PENNY = 5,
    IS_FLAT_ELLIPSOID = 7,
    */
}












// ***********************************************************************************************************************
// ***********************************   INTERSECTION   ******************************************************************
// ***********************************************************************************************************************
void intersection_3d_core (double c1x, double c1y, double c1z,  double a1x, double a1y, double a1z,  double e1x, double e1y, double e1z,
                           double c2x, double c2y, double c2z,  double a2x, double a2y, double a2z,  double e2x, double e2y, double e2z,
                           bool expected, bool ee, const char *s)
{
  Problem *p = new Problem;

  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set ();
  p->set_dimension (3);
  p->set_diffType (DT_NUMERICAL);
  p->set_matrix_E_nu (1.0, 0.4);
  double rs [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  p->set_numberOfRemoteStrains (1);
  p->set_RemoteStrain (0, rs);

  // inclusions
  p->set_number_of_inclusions (2); // id  coordinates     E    nu    semiaxes        eAngles
  p->set_inclusion_all               (0,  c1x, c1y, c1z,  5.0, 0.4,  a1x, a1y, a1z,  e1x, e1y, e1z);
  p->set_inclusion_all               (1,  c2x, c2y, c2z,  5.0, 0.4,  a2x, a2y, a2z,  e2x, e2y, e2z);

  p->input_data_initialize_and_check_consistency ();
  // **************   DATA INPUT VIA FUNCTIONS   ******************

  p->convert_to_equivalent_problem ();
  p->print_visualization("intersection_3d_visualization_DELETE.vtk", 12);

  int result = MesoFace::ellipsoids_overlap ( p->give_inclusion(0), p->give_inclusion(1) );

  if (result == -1) {
    if (ee == false  ||  (ee == true  &&  expected == false)) {
      printf (" NOT CONVERGED %s", s);
      return;
    }
  }
  else {
    if ((bool)result != expected) {
      printf ("  \nERRORR - expectation of overlap is %s\n", expected ? "YES" : "NO");
      return;
    }
  }

  printf (" OK%s", s);

  //get_confirm ("run visualization", true);

  delete p;
}

// 2d uloha pocitana 2d a 3d mumechem, postupna rotace inkluze, tisk hodnot v jednom bode
// Pokud je 3d uloha derivovana numericky,  vysledna krivka je zubata.
// Pokud je 3d uloha derivovana analyticky, vysledna krivka je hladka.
void x2Dx3D_rotation(void)
{
  int i, k;
  Problem *analFce2D, *analFce3D;

  const double rs2d[] = { 1.0, 0.0, 0.0, 0.0 };
  const double rs3d[] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

  double *point = new double[3];
  point[0] = 3.5;
  point[1] = 1.5;
  point[2] = 0.;

  double *displc_2D = new double[2];   double *displc_3D = new double[3];
  double *strain_2D = new double[4];   double *strain_3D = new double[9];
  double *stress_2D = new double[4];   double *stress_3D = new double[9];

  int _2D[] = { 0, 3, 1 };
  int _3D[] = { 0, 4, 1 };

  double **results = NULL;
  AllocateArray2D(results, 181, 16);

  DiffTypes dt2d = DT_NUMERICAL;
  DiffTypes dt3d = DT_NUMERICAL;

  char vtkAmndFile2d [] = "rotation/amnd.2d.vtk";
  char vtkAmndFile3d [] = "rotation/amnd.3d.vtk";
  for (int ii = 0; ii <= 180; ii++){
    // vytvoreni prikladu
    eqiv_gener_1x2dI_2d (vtkAmndFile2d, 1.0, 0.4, rs2d, 0.0, 0.0, 2.4, 0.3, 1.0, 0.5, ii*1.0, dt2d);
    eqiv_gener_1x2dI_3d (vtkAmndFile3d, 1.0, 0.4, rs3d, 0.0, 0.0, 2.4, 0.3, 1.0, 0.5, ii*1.0, dt3d);

    analFce2D = (new Problem())->read_inclusions_plus_initialize_and_print (vtkAmndFile2d, NULL, dt2d);
    analFce3D = (new Problem())->read_inclusions_plus_initialize_and_print (vtkAmndFile3d, NULL, dt3d);

    analFce2D->giveFieldsOfPointOneRS(displc_2D, strain_2D, stress_2D, point, 'p', 0, PFCM_OPTIMIZED);
    analFce3D->giveFieldsOfPointOneRS(displc_3D, strain_3D, stress_3D, point, 'p', 0, PFCM_OPTIMIZED);

    for (k = 0; k < 3; k++){
      results[ii][2 * k     ] = strain_2D[_2D[k]];
      results[ii][2 * k +  1] = strain_3D[_3D[k]];
    }
    for (k = 0; k < 3; k++){
      results[ii][2 * k +  6] = stress_2D[_2D[k]];
      results[ii][2 * k +  7] = stress_3D[_3D[k]];
    }
    for (k = 0; k < 2; k++){
      results[ii][2 * k + 12] = displc_2D[k];
      results[ii][2 * k + 13] = displc_3D[k];
    }
    delete analFce2D;
    delete analFce3D;
  }
  int index;
  char name[10], filename2d[25], filename3d[25];
  // 0-2 strain, 3-5 stress, 6-7 displacement
  for (index = 0; index < 8; index++) {
    if (index == 0) sprintf(name, "strain1");
    if (index == 1) sprintf(name, "strain2");
    if (index == 2) sprintf(name, "strain3");
    if (index == 3) sprintf(name, "stress1");
    if (index == 4) sprintf(name, "stress2");
    if (index == 5) sprintf(name, "stress3");
    if (index == 6) sprintf(name, "disp1");
    if (index == 7) sprintf(name, "disp2");
    FILE *o2, *o3;
    sprintf (filename2d, "./rotation/%s.2d.dat", name);
    sprintf (filename3d, "./rotation/%s.3d.dat", name);
    o2 = fopen(filename2d, "wt");
    o3 = fopen(filename3d, "wt");
    for (i = 0; i <= 180; i++) {
      fprintf(o2, "%d %.14lf\n", i, results[i][2 * index    ]);
      fprintf(o3, "%d %.14lf\n", i, results[i][2 * index + 1]);
    }
    fclose(o2);
    fclose(o3);
  }
}


void intersection_3d_core_6 (double cx, double cy, double cz,
                             double a1x, double a1y, double a1z,  double g1x, double g1y, double g1z,  double e1x, double e1y, double e1z,
                             double a2x, double a2y, double a2z,  double g2x, double g2y, double g2z,  double e2x, double e2y, double e2z)
{
  double e = 0.01 * (a1y+a2y)/2.0;

  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx + (g1x+g2x  ), cy              , cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  true , true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx + (g1x+g2x+e), cy              , cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  false, true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx - (g1x+g2x  ), cy              , cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  true , true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx - (g1x+g2x+e), cy              , cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  false, true, " ");  // endik

  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy + (g1y+g2y  ), cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  true , true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy + (g1y+g2y+e), cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  false, true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy - (g1y+g2y  ), cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  true , true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy - (g1y+g2y+e), cz              ,  a2x,a2y,a2z,   e2x,e2y,e2z,  false, true, " ");  // endik

  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy              , cz + (g1z+g2z  ),  a2x,a2y,a2z,   e2x,e2y,e2z,  true , true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy              , cz + (g1z+g2z+e),  a2x,a2y,a2z,   e2x,e2y,e2z,  false, true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy              , cz - (g1z+g2z  ),  a2x,a2y,a2z,   e2x,e2y,e2z,  true , true, "");  // endik
  intersection_3d_core (cx              , cy              , cz              ,  a1x,a1y,a1z,   e1x,e1y,e1z,
                        cx              , cy              , cz - (g1z+g2z+e),  a2x,a2y,a2z,   e2x,e2y,e2z,  false, true, " \n");  // endik
}
// Prvni inkluze ma dane pootoceni, pootoceni druhe je zde generovano, pokud predpokladam uhel 90, tak je pouze 6 pozic.
void intersection_3d_core_6x6 (double cx, double cy, double cz,
                               double a1x, double a1y, double a1z,  double g1x, double g1y, double g1z,  double e1x, double e1y, double e1z,
                               double a2x, double a2y, double a2z)
{
  intersection_3d_core_6 (cx,cy,cz, a1x,a1y,a1z, g1x,g1y,g1z, e1x,e1y,e1z, a2x,a2y,a2z,  a2x, a2y, a2z,    0.0,  0.0,  0.0);
  intersection_3d_core_6 (cx,cy,cz, a1x,a1y,a1z, g1x,g1y,g1z, e1x,e1y,e1z, a2x,a2y,a2z,  a2y, a2x, a2z,    0.0,  0.0, 90.0);
  intersection_3d_core_6 (cx,cy,cz, a1x,a1y,a1z, g1x,g1y,g1z, e1x,e1y,e1z, a2x,a2y,a2z,  a2y, a2z, a2x,    0.0, 90.0, 90.0);
  intersection_3d_core_6 (cx,cy,cz, a1x,a1y,a1z, g1x,g1y,g1z, e1x,e1y,e1z, a2x,a2y,a2z,  a2x, a2z, a2y,    0.0, 90.0,  0.0);
  intersection_3d_core_6 (cx,cy,cz, a1x,a1y,a1z, g1x,g1y,g1z, e1x,e1y,e1z, a2x,a2y,a2z,  a2z, a2x, a2y,   90.0, 90.0,  0.0);
  intersection_3d_core_6 (cx,cy,cz, a1x,a1y,a1z, g1x,g1y,g1z, e1x,e1y,e1z, a2x,a2y,a2z,  a2z, a2y, a2x,   90.0, 90.0, 90.0);
}
// Prvni inkluze pootoceni je zde generovano, pokud predpokladam uhel 90, tak je pouze 6 pozic.
void intersection_3d_core_6x6x6 (double cx, double cy, double cz,
                                 double a1x, double a1y, double a1z,
                                 double a2x, double a2y, double a2z)
{
  intersection_3d_core_6x6 (cx,cy,cz, a1x,a1y,a1z,  a1x, a1y, a1z,    0.0,  0.0,  0.0,  a2x,a2y,a2z);
  intersection_3d_core_6x6 (cx,cy,cz, a1x,a1y,a1z,  a1y, a1x, a1z,    0.0,  0.0, 90.0,  a2x,a2y,a2z);
  intersection_3d_core_6x6 (cx,cy,cz, a1x,a1y,a1z,  a1y, a1z, a1x,    0.0, 90.0, 90.0,  a2x,a2y,a2z);
  intersection_3d_core_6x6 (cx,cy,cz, a1x,a1y,a1z,  a1x, a1z, a1y,    0.0, 90.0,  0.0,  a2x,a2y,a2z);
  intersection_3d_core_6x6 (cx,cy,cz, a1x,a1y,a1z,  a1z, a1x, a1y,   90.0, 90.0,  0.0,  a2x,a2y,a2z);
  intersection_3d_core_6x6 (cx,cy,cz, a1x,a1y,a1z,  a1z, a1y, a1x,   90.0, 90.0, 90.0,  a2x,a2y,a2z);
}

void intersection_3d (void)
{
  printf ("\n TOOLS - 3D intersection:\n");

  double  cx = 0.0, cy = 0.0, cz = 0.0;     // shift origin both inclusions
  double  a1x = 3.0, a1y = 2.0, a1z = 1.0;  // 1. inclusion semiaxes
  double  a2x = 3.0, a2y = 2.0, a2z = 1.0;  // 2. inclusion semiaxes

  // 6x6x6 prikladu - dve elipsy se dotykaji nebo tesne nedotykaji ve vrcholech
  // Vypada to, ze pokud lokalni osa x jedne elipsy miri do stredu druhe a elipsy se dotykaji v bode tak to obcas nedokonverguje.
  // Toto se bude dit tak vzacne, asi, ze tam proste hodim error natvrdo a uvidime.
  // V tomto prikladu 6x6x6 se s tim pocita (ee == true) a chyba je potlacena, v jinych prikladech nedokonvergovani hodi error.
  intersection_3d_core_6x6x6(cx,cy,cz, a1x,a1y,a1z, a2x,a2y,a2z);

  //
  intersection_3d_core (cx        , cy      , cz      ,  3.0, 2.0, 1.0,   0.0,  0.0,  0.0,
                        cx + 4.24 , cy + 2.0, cz + 1.0,  3.0, 2.0, 1.0,   0.0,  0.0,  0.0,  true , false, "\n");  // endik
  intersection_3d_core (cx        , cy      , cz      ,  3.0, 2.0, 1.0,   0.0,  0.0,  0.0,
                        cx + 4.26 , cy + 2.0, cz + 1.0,  3.0, 2.0, 1.0,   0.0,  0.0,  0.0,  false, false, "\n");  // endik

  //
  intersection_3d_core (cx        , cy      , cz      ,  3.0, 2.0, 2.0,  30.0, 45.0, 17.0,
                        cx + 4.35 , cy + 2.0, cz + 1.0,  2.0, 2.0, 2.0,   0.0,  0.0,  0.0,  true , false, "\n");  // endik
  intersection_3d_core (cx        , cy      , cz      ,  3.0, 2.0, 2.0,  30.0, 45.0, 17.0,
                        cx + 4.40 , cy + 2.0, cz + 1.0,  2.0, 2.0, 2.0,   0.0,  0.0,  0.0,  false, false, "\n");  // endik


  //
  intersection_3d_core (cx        , cy      , cz      ,  3.0, 2.0, 1.0,  90.0, 90.0, 90.0,
                        cx + 1.68 , cy + 2.0, cz + 1.0,  3.0, 2.0, 1.0,  90.0, 90.0, 90.0,  true , false, "\n");  // endik
  intersection_3d_core (cx        , cy      , cz      ,  3.0, 2.0, 1.0,  90.0, 90.0, 90.0,
                        cx + 1.70 , cy + 2.0, cz + 1.0,  3.0, 2.0, 1.0,  90.0, 90.0, 90.0,  false, false, "\n");  // endik



  printf("\n\n");
}


void intersection_2d_core (double c1x, double c1y,  double a1x, double a1y,  double e1,
                           double c2x, double c2y,  double a2x, double a2y,  double e2,
                           bool expected, const char *s)
{
  Problem *p = new Problem;

  // **************   DATA INPUT VIA FUNCTIONS   ******************
  // problem
  p->set_data_set ();
  p->set_dimension (2);

  p->set_matrix_E_nu (2.0, 0.2);
  double rs [] = { 1.0, 0.0, 0.0, 0.0 };
  p->set_numberOfRemoteStrains (1);
  p->set_RemoteStrain (0, rs);

  // inclusions
  p->set_number_of_inclusions (2); // id  coords     E    nu    semiaxes   eAngle
  p->set_inclusion_all               (0,  c1x, c1y,  1.0, 0.1,  a1x, a1y, e1);
  p->set_inclusion_all               (1,  c2x, c2y,  1.0, 0.1,  a2x, a2y, e2);

  p->input_data_initialize_and_check_consistency ();
  // **************   DATA INPUT VIA FUNCTIONS   ******************

  p->convert_to_equivalent_problem ();
  p->print_visualization("intersection_2d_visualization_DELETE.vtk", 20);

  bool result = MesoFace::ellipses_overlap (p->give_inclusion (0), p->give_inclusion (1));

  if (result != expected) {
    printf ("  \nERRORR - expectation of overlap is %s\n", expected ? "YES" : "NO");
    return;
  }

  printf (" OK%s", s);

  delete p;
}
void intersection_2d_core_4 (double cx, double cy,
                             double a1x, double a1y,  double g1x, double g1y,  double e1,
                             double a2x, double a2y,  double g2x, double g2y,  double e2)
{
  double e = 0.01 * (a1y+a2y)/2.0;

  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx + (g1x+g2x  ), cy              ,  a2x,a2y,  e2,  true , "");  // endik
  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx + (g1x+g2x+e), cy              ,  a2x,a2y,  e2,  false, "");  // endik
  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx - (g1x+g2x  ), cy              ,  a2x,a2y,  e2,  true , "");  // endik
  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx - (g1x+g2x+e), cy              ,  a2x,a2y,  e2,  false, " ");  // endik

  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx              , cy + (g1y+g2y  ),  a2x,a2y,  e2,  true , "");  // endik
  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx              , cy + (g1y+g2y+e),  a2x,a2y,  e2,  false, "");  // endik
  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx              , cy - (g1y+g2y  ),  a2x,a2y,  e2,  true , "");  // endik
  intersection_2d_core (cx              , cy              ,  a1x,a1y,  e1,
                        cx              , cy - (g1y+g2y+e),  a2x,a2y,  e2,  false, " ");  // endik

}

//
void intersection_2d (void)
{
  printf ("\n TOOLS - 2D intersection:\n");

  double  cx = 0.0, cy = 0.0;    // shift origin both inclusions
  double  a1x = 5.0, a1y = 2.0;  // 1. inclusion semiaxes
  double  a2x = 4.0, a2y = 1.5;  // 2. inclusion semiaxes

  // 2x2x4 prikladu - dve elipsy se dotykaji nebo tesne nedotykaji ve vrcholech
  intersection_2d_core_4 (cx, cy,  a1x, a1y,  a1x, a1y,  0.0,  a2x, a2y,  a2x, a2y,  0.0);
  intersection_2d_core_4 (cx, cy,  a1x, a1y,  a1x, a1y,  0.0,  a2x, a2y,  a2y, a2x, 90.0);
  intersection_2d_core_4 (cx, cy,  a1x, a1y,  a1y, a1x, 90.0,  a2x, a2y,  a2x, a2y,  0.0);
  intersection_2d_core_4 (cx, cy,  a1x, a1y,  a1y, a1x, 90.0,  a2x, a2y,  a2y, a2x, 90.0);

  //
  intersection_2d_core (cx       , cy      ,  a1x,a1y,   0.0,
                        cx - 4.62, cy + 3.0,  a2x,a2y,   0.0,  true , "\n");  // endik
  intersection_2d_core (cx       , cy      ,  a1x,a1y,   0.0,
                        cx - 4.66, cy + 3.0,  a2x,a2y,   0.0,  false, "\n");  // endik

  intersection_2d_core (cx       , cy      ,  a1x,a1y,   90.0,
                        cx - 3.952, cy + 2.0,  a2x,a2y, 110.0,  true , "\n");  // endik
  intersection_2d_core (cx       , cy      ,  a1x,a1y,   90.0,
                        cx - 3.952, cy + 2.0,  a2x,a2y, 109.0,  false, "\n");  // endik

  intersection_2d_core (cx       , cy      ,  a1x,a1y,  90.0,
                        cx - 4.50, cy + 2.0,  a2x,a2y,  40.0,  true , "\n");  // endik
  intersection_2d_core (cx       , cy      ,  a1x,a1y,  90.0,
                        cx - 4.52, cy + 2.0,  a2x,a2y,  40.0,  false, "\n");  // endik

  printf("\n\n");
}

// XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX
// Tisk grafu ovlivneni, viz p->give_ovlivneni (x, result)
void ovlivneni (void)
{
  //char vtkAmndFile [] = "ovlivneni_jednaInkluze.amnd.vtk";

  Problem *p;
  p = new Problem();

  // **************   DATA INPUT   ******************
  // set problem
  p->set_data_set();
  p->set_dimension(3);
  p->set_diffType(DT_ANALITICAL);
  p->set_matrix_E_nu(1.0, 0.4);
  //
  p->set_numberOfRemoteStrains(1);
  double rs1 [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  p->set_RemoteStrain(0, rs1);

  // loop over all inclusions
  p->set_number_of_inclusions(1);

  // first inclusion
  p->give_inclusion(0)->set_centroids(0.0, 0.0, 0.0);
  p->give_inclusion(0)->set_Inclusion_shape(IS_PROLATE_SPHEROID);
  p->give_inclusion(0)->set_Youngs_modulus(10.0);
  p->give_inclusion(0)->set_Poissons_ratio(0.4);
  p->give_inclusion(0)->set_Semiaxes_dimensions(1.2, 1.0, 1.0);
  p->give_inclusion(0)->set_Euller_angles_deg(0.0, 90.0, 90.0);

  p->input_data_initialize_and_check_consistency();
  // **************   DATA INPUT   ******************

  p->set_SBA_optimized(false);
  p->convert_to_equivalent_problem();

  long i;

  // menit dle libosti
  double minCoord = -4.0;
  double maxCoord =  4.0;
  double step = 0.1;
  double x[3];
  double result[6] = { 0., 0., 0., 0., 0., 0. };

  FILE * out1 = NULL;
  FILE * out2 = NULL;
  char name[100];


  for (i = 0; i < 2; i++) { // over X & Y
    sprintf(name, "pokus_graf_2.0_X_%02ld.dat", i);    OpenFile(out1, name, _WRITE_);
    sprintf(name, "pokus_graf_2.0_Y_%02ld.dat", i);    OpenFile(out2, name, _WRITE_);

    double coord = minCoord;

    while (coord <= maxCoord) { // over points
      x[0] = x[1] = x[2] = 0.0;
      x[i] = coord;

      p->give_ovlivneni (x, result);

      fprintf(out1, "%e\t%e\n", x[i], result[0]);
      fprintf(out2, "%e\t%e\n", x[i], result[1]);

      coord += step;

    } // end of while loop over points

    CloseFile(out1);
    CloseFile(out2);
  } // end of for loop over X & Y

  delete p;

  return;
}

/*end of file*/