mstruct_generateur_polycristaux.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mstruct_generateur_polycristaux.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2025
//####//  Derniere modification par francois
//####//  mer 14 mai 2025 17:54:57 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
  #include "gestionversion.h"
  #include "mg_file.h"
  #include "mstruct_generateur_polycristaux.h"
  #include "tpl_grille.h"
  #ifdef PROJECT_POLYCRISTAUX
  #include "polycristal.h"
  #endif
  #include "occ_import.h"
  #include "fct_taille_fem_solution_generateur_gradient.h"
  #include "mailleur0d.h"
  #include "mailleur1d.h"
  #include "mailleur2d.h"
  #include "mailleur3d.h"
  #include "mailleur_fem.h"
  #include "mailleur_analyse.h"
  #include "mg_export.h"
  #include "mstruct_outils.h"
  #include "mstruct_definition.h"
  #include "mgaster.h"
  #include "lc_point.h"
  #include "mailleur_analyse.h"
  #include <fstream>
  #include <math.h>
  #include <random>
  #include "mstruct_borne.h"
#include <sys/stat.h>
#include <mutex>
 
 
MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT::MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT()
{
nbcristaux=-1;
degre=-1;
dg=-1;
KDH=-1.;
GDH=-1.;
EDH=-1.;
KCH=-1.;
GCH=-1.;
ECH=-1.;
nbphase=0;
}
 
MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT::MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT(MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT& mdd):nbcristaux(mdd.nbcristaux),degre(mdd.degre),dg(mdd.dg),KDH(mdd.KDH),GDH(mdd.GDH),EDH(mdd.EDH),KCH(mdd.KCH),GCH(mdd.GCH),ECH(mdd.ECH),nbphase(mdd.nbphase),nbcristauxphase(mdd.nbcristauxphase),masse(mdd.masse),volume(mdd.volume),cdm(mdd.cdm)
{
memcpy(tens1,mdd.tens1,sizeof(double)*9);
memcpy(tens2,mdd.tens2,sizeof(double)*9);
memcpy(tens3,mdd.tens3,sizeof(double)*9);
        
}
 
MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT::~MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT()
{
}
 
 
MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE::MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE() 
{       
Evoigt=-1;
Ereuss=-1;
Kvoigt=-1;
Kreuss=-1;
Gvoigt=-1;
Greuss=-1;
nuvoigt=-1;
nureuss=-1;
}
 
MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE::MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE(MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE& mdd):Evoigt(mdd.Evoigt),Ereuss(mdd.Ereuss),Kvoigt(mdd.Kvoigt),Kreuss(mdd.Kreuss),Gvoigt(mdd.Gvoigt),Greuss(mdd.Greuss),nuvoigt(mdd.nuvoigt),nureuss(mdd.nureuss)
{
}
 
MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE::~MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE() 
{
}
 
 
 
  MSTRUCT_GENERATEUR_POLYCRISTAUX::MSTRUCT_GENERATEUR_POLYCRISTAUX()
  {
  ini_param(NULL);
  }
 
 
  MSTRUCT_GENERATEUR_POLYCRISTAUX::MSTRUCT_GENERATEUR_POLYCRISTAUX(char* nomparam)
  {
  ini_param(nomparam);
  }
  
  MSTRUCT_GENERATEUR_POLYCRISTAUX::MSTRUCT_GENERATEUR_POLYCRISTAUX(char* fichierin, char* fichierout)
  {
          MG_FILE gest(fichierin);
          MG_GEOMETRIE *mggeo=gest.get_mg_geometrie(0);
          FILE *out=fopen(fichierout,"wt");
          LISTE_MG_POINT::iterator it;
          for (MG_POINT* pt=mggeo->get_premier_point(it);pt!=NULL;pt=mggeo->get_suivant_point(it))
                 if (pt->get_type()==MG_ELEMENT_GEOMETRIQUE::TYPE_ELEMENT_GEOMETRIQUE::LC_POINT)
                 {
                   double xyz[3];
                   pt->evaluer(xyz);
                   fprintf(out,"%e %e %e\n",xyz[0],xyz[1],xyz[2]);  
                 }
           fclose(out);
  }
 
 
  MSTRUCT_GENERATEUR_POLYCRISTAUX::MSTRUCT_GENERATEUR_POLYCRISTAUX(MSTRUCT_GENERATEUR_POLYCRISTAUX& mdd):param(mdd.param)
  {
 
  }
 
  MSTRUCT_GENERATEUR_POLYCRISTAUX::~MSTRUCT_GENERATEUR_POLYCRISTAUX()
  {
 
  }
 
  void MSTRUCT_GENERATEUR_POLYCRISTAUX::active_affichage(fonction_affiche *fonc)
  {
      affiche2=fonc;
      affichageactif=1;
  }
 
  void MSTRUCT_GENERATEUR_POLYCRISTAUX::affiche(char *mess)
  {
    if (affichageactif==1) affiche2(mess);
  }
 
  void MSTRUCT_GENERATEUR_POLYCRISTAUX::ini_param(char *nom)
  {
  param.ajouter("nommateriau","Cuivre",OT_PARAMETRES::STRING,"nom du matériaux");
  //param.ajouter("bornereuss",109e9,OT_PARAMETRES::DOUBLE,"Valeur de la borne de Reuss en Pa");
  //param.ajouter("bornevoigt",143e9,OT_PARAMETRES::DOUBLE,"Valeur de la borne de Voigt en Pa");
  param.ajouter("boitexmin",0,OT_PARAMETRES::DOUBLE,"valeur x min de l'échantillon");
  param.ajouter("boiteymin",0,OT_PARAMETRES::DOUBLE,"valeur y min de l'échantillon");
  param.ajouter("boitezmin",0,OT_PARAMETRES::DOUBLE,"valeur z min de l'échantillon");
  param.ajouter("boitexmax",1,OT_PARAMETRES::DOUBLE,"valeur x max de l'échantillon");
  param.ajouter("boiteymax",1,OT_PARAMETRES::DOUBLE,"valeur y max de l'échantillon");
  param.ajouter("boitezmax",1,OT_PARAMETRES::DOUBLE,"valeur z max de l'échantillon");
  param.ajouter("nb_cristaux",100,OT_PARAMETRES::DOUBLE,"Nombre de cristaux à tirer aléatoirement");
  param.ajouter("rayon",0.05,OT_PARAMETRES::DOUBLE,"Distance séparant deux points tirés aléatoirement en % de la longueur totale");
  param.ajouter("nomfichier","nomcao",OT_PARAMETRES::STRING,"Nom du fichier de l'étude sans extension");
  param.ajouter("enregistrementstep",0.,OT_PARAMETRES::DOUBLE,"0. pas de fichier step 1. avec fichier step");
  param.ajouter("epsfusion",1e-6,OT_PARAMETRES::DOUBLE,"epsilon de la fusion des noeuds des Voronoi");
  param.ajouter("dg",0.1,OT_PARAMETRES::DOUBLE,"ecart nodale pour la carte de taille");
  param.ajouter("gradient",0.25,OT_PARAMETRES::DOUBLE,"gradient max visé pour la carte de taille");
  param.ajouter("bornemin_carte",0.1,OT_PARAMETRES::DOUBLE,"Borne minimal de l'écart nodal");
  param.ajouter("type_bornemin",1,OT_PARAMETRES::DOUBLE,"0. borne min absolue 1. borne min relative");
  param.ajouter("Nbphase",2.0,OT_PARAMETRES::DOUBLE,"Nombre de phase du polycristal");
  std::vector<double> lstphase;
  lstphase.push_back(0.68);
  lstphase.push_back(0.32);
  param.ajouter("Proportionphase",lstphase,OT_PARAMETRES::DOUBLE,"Proportion de chacune des phases en masse");
  std::vector<double> lstdens;
  lstdens.push_back(5.);
  lstdens.push_back(12.);
  param.ajouter("Densite",lstdens,OT_PARAMETRES::DOUBLE,"Densité de chacune des phases");
  std::vector<double> lst1,lst2,lst3;
  lst1.push_back(66.69e9);
  lst1.push_back(120.1e9);
  lst2.push_back(0.4189);
  lst2.push_back(0.3);
  lst3.push_back(75.4e9);
  lst3.push_back(85.2e9);
  
  param.ajouter("E1",lst1,OT_PARAMETRES::DOUBLE,"Module d'Young dans la direction 1 du polycristal");
  param.ajouter("E2",lst1,OT_PARAMETRES::DOUBLE,"Module d'Young dans la direction 2 du polycristal");
  param.ajouter("E3",lst1,OT_PARAMETRES::DOUBLE,"Module d'Young dans la direction 3 du polycristal");
  param.ajouter("n1",lst2,OT_PARAMETRES::DOUBLE,"Coefficent de Poisson dans le plan 12 du polycristal");
  param.ajouter("n2",lst2,OT_PARAMETRES::DOUBLE,"Coefficent de Poisson dans le plan 23 du polycristal");
  param.ajouter("n3",lst2,OT_PARAMETRES::DOUBLE,"Coefficent de Poisson dans le plan 13 du polycristal");
  param.ajouter("G1",lst3,OT_PARAMETRES::DOUBLE,"Module de Cisaillement dans le plan 12 du polycristal");
  param.ajouter("G2",lst3,OT_PARAMETRES::DOUBLE,"Module de Cisaillement dans le plan 23 du polycristal");
  param.ajouter("G3",lst3,OT_PARAMETRES::DOUBLE,"Module de Cisaillement dans le plan 13 du polycristal");
  param.ajouter("Degre",1,OT_PARAMETRES::DOUBLE,"1 maillage lineaire 2 maillage quadratique");
  param.ajouter("fichiercumul","cumul",OT_PARAMETRES::STRING,"Nom du fichier sans extension du cumul des resultats");
  param.ajouter("enregistrerres","1.0",OT_PARAMETRES::DOUBLE,"0. pas d'enregistrement des résultats et statistiques 1. enregistrement");
  param.ajouter("paramaster","param_aster.txt",OT_PARAMETRES::STRING,"Nom du fichier des parametres aster");
  param.ajouter("reprise_geo",0.,OT_PARAMETRES::DOUBLE,"0. nouvelle géométrie 1. reprise de la dernière géométrie");
  param.ajouter("reprise_ori",0.,OT_PARAMETRES::DOUBLE,"0. nouvelle orientation 1. reprise des dernières orientations");
  param.ajouter("reprise_mai",0.,OT_PARAMETRES::DOUBLE,"0. nouveau maillage 1. reprise du dernier maillage");
  param.ajouter("operationaeffectuer",3.,OT_PARAMETRES::DOUBLE,"0. rien 1. generation 2. maillage 3. calcul");
  if (nom!=NULL) param.lire(nom);
  }
 
 
  void MSTRUCT_GENERATEUR_POLYCRISTAUX::cree_param(char* nom,char *nomexe)
  {
  affiche((char*)"\nCreation du fichier de parametre");
  param.change_valeur("operationaeffectuer", 2.);
  param.change_valeur("enregistrerres", 0.);
  param.enregistrer(nom);
  
  char nom2[500];
  strcpy(nom2, nom);
  strcat(nom2, "2");
  FILE *paraout2=fopen(nom2,"wt");
  fprintf(paraout2,"-> %s\n",nom);
  fprintf(paraout2,"operationaeffectuer = 3.000000\n");
  fprintf(paraout2,"reprise_geo = 1.000000\n");
  fprintf(paraout2,"reprise_ori = 1.000000\n");
  fprintf(paraout2,"reprise_mai = 1.000000\n");
  fprintf(paraout2,"enregistrerres = 1.000000\n");
  
  
  fclose(paraout2);
  
  
  
  
  FILE *out = fopen("vasy", "wt");
  fprintf(out, "#!/bin/bash\n");
fprintf(out, "for ((c=1;c<=$1;c++))\n");
fprintf(out, "do\n");
fprintf(out, "timeout 900s %s  -polycristaux -out cristal.magic -param %s\n", nomexe,nom);
fprintf(out, "VAR=$?\n");
fprintf(out, "if (( $VAR )) ; then\n"); 
fprintf(out, "          echo \"Timeout\"\n");
fprintf(out, "else\n"); 
fprintf(out, "                  %s  -polycristaux -out cristal.magic -param %s\n", nomexe,nom2);
fprintf(out, "fi\n");
fprintf(out, "done\n");
 
  
  fclose(out);
  out = fopen("vasy_par", "wt");
  fprintf(out, "#!/bin/bash\n");
fprintf(out, "VAR=1\n");
fprintf(out, "while (($VAR))\n");
fprintf(out, "do\n");
bool chemindanspath=false;
if ((nomexe[0]>='a') && (nomexe[0]<='z')) chemindanspath=true;
if ((nomexe[0]>='A') && (nomexe[0]<='Z')) chemindanspath=true;
if (chemindanspath)
        fprintf(out, "timeout 1800s %s  -polycristaux -out cristal.magic -param %s\n", nomexe,nom);
else 
        fprintf(out, "timeout 1800s ../%s  -polycristaux -out cristal.magic -param %s\n", nomexe,nom);
fprintf(out, "VAR=$?\n");
fprintf(out, "if (( $VAR )) ; then\n"); 
fprintf(out, "          echo \"Timeout\"\n");
fprintf(out, "else\n"); 
if (chemindanspath)
        fprintf(out, "                  %s  -polycristaux -out cristal.magic -param %s\n", nomexe,nom2);
else
        fprintf(out, "                  ../%s  -polycristaux -out cristal.magic -param %s\n", nomexe,nom2);
fprintf(out, "fi\n");
fprintf(out, "done\n");
 
  
  fclose(out);
  chmod ("vasy",strtol("0777", NULL,8));
  chmod ("vasy_par",strtol("0777", NULL,8));
  MG_EXPORT exp;
  exp.change_param_aster((char*)"Noeud_ele",2.);
  exp.ecrire_params_aster((char*)"param_aster.txt");
  
  }
 
 
        void MSTRUCT_GENERATEUR_POLYCRISTAUX::genere(char *nommagic)
        {
        int reprendgeo=(int)param.get_valeur("reprise_geo");
        int reprendmai=(int)param.get_valeur("reprise_mai");
        int reprendori=(int)param.get_valeur("reprise_ori");
 
          
 
 
        #ifdef PROJECT_POLYCRISTAUX
    MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT statres;
        MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE statborne;
        MG_GESTIONNAIRE *gest;
        int afaire=(int)param.get_valeur("operationaeffectuer");
        MG_GEOMETRIE* mggeo;
        double dg=param.get_valeur("dg");
    statres.dg=dg;
        double epsfusion=param.get_valeur("epsfusion");
    std::string nomfichier=param.get_nom("nomfichier");
        std::string nomfichierbrep=nomfichier+".brep";
        std::string nomfichierocaf=nomfichier+".ocaf";
        std::string nomfichiercarte=nomfichier+".sol";
        std::string nomfichierborne=nomfichier+".borne";
        double xmin=param.get_valeur("boitexmin");
        double ymin=param.get_valeur("boiteymin");
        double zmin=param.get_valeur("boitezmin");
        double xmax=param.get_valeur("boitexmax");
        double ymax=param.get_valeur("boiteymax");
        double zmax=param.get_valeur("boitezmax");
        BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
        BOITE_3D boite2=boite;
        boite2.change_grosseur(1.01);
    if (afaire>0) 
        {
        std::random_device seed;
        std::mt19937_64 generateur(seed());  
          affiche((char*)"\nInitialisation");
            
          if (reprendgeo==0)
          {
          gest=new MG_GESTIONNAIRE;
          TPL_GRILLE<POLYCRISTAUX::POINT_TIRE*> grille;
          grille.initialiser(boite2.get_xmin(),boite2.get_ymin(),boite2.get_zmin(),boite2.get_xmax(),boite2.get_ymax(),boite2.get_zmax(),100,100,100);
          /*POLYCRISTAUX::POINT_TIRE p1(0,0,0,0);grille.inserer(&p1);
          POLYCRISTAUX::POINT_TIRE p2(1,1000,0,0);grille.inserer(&p2);
          POLYCRISTAUX::POINT_TIRE p3(2,1000,1000,0);grille.inserer(&p3);
          POLYCRISTAUX::POINT_TIRE p4(3,0,1000,0);grille.inserer(&p4);
          POLYCRISTAUX::POINT_TIRE p5(4,0,0,0);grille.inserer(&p5);
          POLYCRISTAUX::POINT_TIRE p6(5,1000,0,0);grille.inserer(&p6);
          POLYCRISTAUX::POINT_TIRE p7(6,1000,1000,0);grille.inserer(&p7);
          POLYCRISTAUX::POINT_TIRE p8(7,0,1000,0);grille.inserer(&p8);*/
          double dx=boite.get_xmax()-boite.get_xmin();  
          double dy=boite.get_ymax()-boite.get_ymin();  
          double dz=boite.get_zmax()-boite.get_zmin();  
          affiche((char*)"Generation de cristaux aléatoires"); 
          double rayon=param.get_valeur("rayon")*0.3333333333333333333*(dx+dy+dz);
          std::uniform_real_distribution<double> distribution_position_x(boite.get_xmin(),boite.get_xmax());
          std::uniform_real_distribution<double> distribution_position_y(boite.get_ymin(),boite.get_ymax());
          std::uniform_real_distribution<double> distribution_position_z(boite.get_zmin(),boite.get_zmax());
          int nbinsere=0;
          int nbainsere=param.get_valeur("nb_cristaux");
          std::vector<POLYCRISTAUX::POINT_TIRE*> listepointainserer;
          std::vector<double> listedepointvoro;
          listedepointvoro.push_back(boite.get_xmin());listedepointvoro.push_back(boite.get_ymin());listedepointvoro.push_back(boite.get_zmin());
          listedepointvoro.push_back(boite.get_xmin());listedepointvoro.push_back(boite.get_ymax());listedepointvoro.push_back(boite.get_zmin());
          listedepointvoro.push_back(boite.get_xmax());listedepointvoro.push_back(boite.get_ymax());listedepointvoro.push_back(boite.get_zmin());
          listedepointvoro.push_back(boite.get_xmax());listedepointvoro.push_back(boite.get_ymin());listedepointvoro.push_back(boite.get_zmin());
          listedepointvoro.push_back(boite.get_xmin());listedepointvoro.push_back(boite.get_ymin());listedepointvoro.push_back(boite.get_zmax());
          listedepointvoro.push_back(boite.get_xmin());listedepointvoro.push_back(boite.get_ymax());listedepointvoro.push_back(boite.get_zmax());
          listedepointvoro.push_back(boite.get_xmax());listedepointvoro.push_back(boite.get_ymax());listedepointvoro.push_back(boite.get_zmax());
          listedepointvoro.push_back(boite.get_xmax());listedepointvoro.push_back(boite.get_ymin());listedepointvoro.push_back(boite.get_zmax());
        
          while (nbinsere<nbainsere)
            {
            double x=distribution_position_x(generateur);
            double y=distribution_position_y(generateur);
            double z=distribution_position_z(generateur);
            TPL_MAP_ENTITE<POLYCRISTAUX::POINT_TIRE*> lsttrouve;
            grille.rechercher(x,y,z,rayon,lsttrouve);
            TPL_MAP_ENTITE<POLYCRISTAUX::POINT_TIRE*>::ITERATEUR it;
            bool bonpoint=true;
            for (POLYCRISTAUX::POINT_TIRE* pt=lsttrouve.get_premier(it);pt!=NULL;pt=lsttrouve.get_suivant(it))
                {
                OT_VECTEUR_3D vec(x-pt->x,y-pt->y,z-pt->z);
                if (vec.get_longueur()<rayon) bonpoint=false;
                }
            if (bonpoint)
                {
                listedepointvoro.push_back(x);
                listedepointvoro.push_back(y);
                listedepointvoro.push_back(z);
                POLYCRISTAUX::POINT_TIRE *pt=new POLYCRISTAUX::POINT_TIRE(nbinsere+8,x,y,z);
                grille.inserer(pt);
                listepointainserer.push_back(pt);
                nbinsere++;
                }
            }
            
            for (int i=0;i<listepointainserer.size();i++)
              delete listepointainserer[i];
            
        char message[1000];
            sprintf(message,"   Construction des Voronoi avec borne de fusion de  %lf",epsfusion);
        affiche(message);
            Polycristal poly(listedepointvoro,(char*)nomfichierbrep.c_str(),epsfusion);
            if (affichageactif) poly.active_affichage(affiche2);
            bool avecstep=false; 
            int valeur=(int)param.get_valeur("enregistrementstep");
            if (valeur==0) avecstep=false;
            if (valeur==1) avecstep=true;
            poly.construit(avecstep);
            OCC_IMPORT occ_import;
            mggeo = occ_import.importer_fichier_brep_V2017(*gest,(char*)nomfichierbrep.c_str(),(char*)nomfichierocaf.c_str(),FICHIEROCC,1e-6,false);
        occ_import.importer_triangulation_V2017(*gest,mggeo,0.1,1);
            mggeo->change_valeur_unite(1.);
            int nb=listedepointvoro.size()/3;
            for (int i=8;i<nb;i++)
                  {
                  double x=listedepointvoro[3*i];  
                  double y=listedepointvoro[3*i+1];  
                  double z=listedepointvoro[3*i+2];
                  LC_POINT* pt=new LC_POINT(x,y,z);
                  mggeo->ajouter_mg_point(pt);
                  }
            
          }
          else 
            {
            affiche((char*)"Reprise de calcul d'un ancien polycristal"); 
            gest=new MG_FILE(nommagic);
            mggeo=gest->get_mg_geometrie(0);
            
            }
          statres.nbcristaux=mggeo->get_nb_mg_volume();
 
      char mess[500];
          sprintf(mess,"   Nombre de cristaux MAGiC assemblés = %d ", mggeo->get_nb_mg_volume());;
          affiche(mess);
        double tens1[9],tens2[9],tens3[9];
        for (int i=0;i<9;i++)
                {
                tens1[i]=0.;  
                tens2[i]=0.;  
                tens3[i]=0.;  
                }
        if (reprendori==0)
        {
          std::uniform_real_distribution<double> distribution_angle1(0.,360.);
          std::uniform_real_distribution<double> distribution_angle2(-1.,1.);
          std::uniform_real_distribution<double> distribution_angle3(0.,360.);
          affiche((char*)"Conditions aux limites");
          LISTE_MG_VOLUME::iterator it;
          int nbphase=(int)param.get_valeur("Nbphase");
      statres.nbphase=nbphase;
      statres.volume.resize(statres.nbphase+1);
      statres.masse.resize(statres.nbphase+1);
      statres.nbcristauxphase.resize(statres.nbphase+1);
      statres.cdm.resize(statres.nbphase+1);
      //std::discrete_distribution<int> d({20, 80});
      //std::map<int, int> m;
      std::vector<double> volume(nbphase);
      std::vector<double> tirage(nbphase);
      std::vector<OT_VECTEUR_3D> cdm(nbphase);
      std::vector<double> cible(nbphase);
      std::vector<int> ordrecible(nbphase);
      std::vector<double> densite(nbphase);
      std::vector<int> nbcristaux(nbphase);
      //OT_TENSEUR rigtot(6,6),soutot(6,6);
          //MSTRUCT_BORNE  borne;
          for (int i=0;i<nbphase;i++)
        {
        volume[i]=0.;
        tirage[i]=0.;
        cdm[i].change_x(0.);
        cdm[i].change_y(0.);
        cdm[i].change_z(0.);
        cible[i]=param.get_valeur("Proportionphase",i);
        ordrecible[i]=i;
        densite[i]=param.get_valeur("Densite",i);
        nbcristaux[i]=0;
                /*double E1=param.get_valeur("E1",i);
                double E2=param.get_valeur("E2",i);
                double E3=param.get_valeur("E3",i);
                double n1=param.get_valeur("n1",i);
                double n2=param.get_valeur("n2",i);
                double n3=param.get_valeur("n3",i);
                double G1=param.get_valeur("G1",i);
                double G2=param.get_valeur("G2",i);
                double G3=param.get_valeur("G3",i);
                OT_TENSEUR rig(6,6),rigphase(6,6),souphase(6,6);
                borne.get_tenseur_rigidite(E1/1e9,E2/1e9,E3/1e9,2.*G1/1e9,2.*G2/1e9,2.*G3/1e9,n1,n2,n3,rig);
                borne.get_VoigtReuss_poly(rig,rigphase,souphase);
                rigphase=cible[i]*rigphase;
        souphase=cible[i]*souphase;
                rigtot=rigtot+rigphase;
        soutot=soutot+souphase;*/
        }
      /*double Ev,Gv,Kv,nuv;  
          double Er,Gr,Kr,nur;  
          borne.get_VoigtReuss_poly(rigtot,soutot,Kv,Gv,Ev,nuv,Kr,Gr,Er,nur);
          statborne.Evoigt=Ev*1e9;
      statborne.Kvoigt=Kv/3.*1e9;
          statborne.Gvoigt=Gv/2.*1e9;
          statborne.nuvoigt=nuv;
          statborne.Ereuss=Er*1e9;
          statborne.Kreuss=Kr/3.*1e9;
          statborne.Greuss=Gr/2.*1e9;
          statborne.nureuss=nur;
          
          FILE *inborne=fopen(nomfichierborne.c_str(),"wb");
          fwrite(&statborne,sizeof(statborne),1,inborne);
          fclose(inborne);*/
          
      std::shuffle(ordrecible.begin(),ordrecible.end(),generateur);
      std::vector<MG_VOLUME*> listevolume;
      for (MG_VOLUME* vol=mggeo->get_premier_volume(it);vol!=NULL;vol=mggeo->get_suivant_volume(it))
          listevolume.push_back(vol);
          std::shuffle(listevolume.begin(),listevolume.end(),generateur);
      double voltot=(boite.get_xmax()-boite.get_xmin())*(boite.get_ymax()-boite.get_ymin())*(boite.get_zmax()-boite.get_zmin());
      double nume=voltot;
      std::vector<double> deno(nbphase);
      for (int i=0;i<nbphase;i++)
          deno[i]=1.;
      for (int i=0;i<nbphase;i++)
        {
        nume=nume*densite[i];
        for (int j=0;j<nbphase;j++)
            if (i!=j) deno[j]=deno[j]*densite[i];
                else deno[j]=deno[j]*cible[i];
        }
      double massetot=0.;
      for (int i=0;i<nbphase;i++)
          massetot=massetot+deno[i];
      massetot=nume/massetot;
          
    
      for (int i=0;i<listevolume.size();i++)
          {
      MG_VOLUME *vol=listevolume[i];    
      int tir=-1,j=0;
      while (tir==-1)
        {
        if (tirage[ordrecible[j]]/massetot<cible[ordrecible[j]]) tir=ordrecible[j];
        j++;
        //if (j>=nbphase) tir=j-1;
        }
      double volu,masse;
      OT_VECTEUR_3D cdm2;
      OT_MATRICE_3D m1,m2;
      double dens=densite[tir];
      vol->get_propriete_massique(gest->get_mg_maillage(0),volu,masse,cdm2,m1,m2,dens,-1);    
      volume[tir]=volume[tir]+volu;  
      tirage[tir]=tirage[tir]+masse;  
      cdm[tir]=cdm[tir]+masse*cdm2;
      nbcristaux[tir]=nbcristaux[tir]+1;
      
      
          for (int i=0;i<vol->get_nb_ccf();i++) 
            vol->supprimer_ccf(0);
          double E1=param.get_valeur("E1",tir);
          double E2=param.get_valeur("E2",tir);
          double E3=param.get_valeur("E3",tir);
          double n1=param.get_valeur("n1",tir);
          double n2=param.get_valeur("n2",tir);
          double n3=param.get_valeur("n3",tir);
          double G1=param.get_valeur("G1",tir);
          double G2=param.get_valeur("G2",tir);
          double G3=param.get_valeur("G3",tir);
          vol->ajouter_ccf((char*)"Ph",(double)tir);
          vol->ajouter_ccf((char*)"E1",E1);
          vol->ajouter_ccf((char*)"E2",E2);
          vol->ajouter_ccf((char*)"E3",E3);
          vol->ajouter_ccf((char*)"n1",n1);
          vol->ajouter_ccf((char*)"n2",n2);
          vol->ajouter_ccf((char*)"n3",n3);
          vol->ajouter_ccf((char*)"G1",G1);
          vol->ajouter_ccf((char*)"G2",G2);
          vol->ajouter_ccf((char*)"G3",G3);
          vol->ajouter_ccf((char*)"Ro",dens);
          double  a1=distribution_angle1(generateur);
          double  a3=distribution_angle3(generateur);
          double  a2=distribution_angle2(generateur);
          a2=acos(a2)/M_PI*180.;
          vol->ajouter_ccf((char*)"a1",a1);
          vol->ajouter_ccf((char*)"a2",a2);
          vol->ajouter_ccf((char*)"a3",a3);
          double x1[3],x2[3],x3[3];
          double psi=a1/180.*M_PI;
          double teta=a2/180.*M_PI;
          double phi=a3/180.*M_PI;
          x1[0]=cos(phi)*cos(psi)-sin(phi)*cos(teta)*sin(psi);
          x1[1]=cos(phi)*sin(psi)+sin(phi)*cos(teta)*cos(psi);
          x1[2]=sin(phi)*sin(teta);
          x2[0]=-sin(phi)*cos(psi)-cos(phi)*cos(teta)*sin(psi);
          x2[1]=-sin(phi)*sin(psi)+cos(phi)*cos(teta)*cos(psi);
          x2[2]=cos(phi)*sin(teta);
          x3[0]=sin(teta)*sin(psi);
          x3[1]=-sin(teta)*cos(psi);
          x3[2]=cos(teta);
          for (int i=0;i<3;i++)
            for (int j=0;j<3;j++)
              {
              tens1[i*3+j]+=x1[i]*x1[j];
              tens2[i*3+j]+=x2[i]*x2[j];
              tens3[i*3+j]+=x3[i]*x3[j];
              }
              
            }
      double massetotale=0.;
      double volumetotal=0.;
      for (int i=0;i<nbphase;i++)
          {
          volumetotal=volumetotal+volume[i];
          massetotale=massetotale+tirage[i];
          }
      char mess[500];
      sprintf(mess,"   Masse estimée %lf  calculée %lf",massetot,massetotale);
      affiche(mess);
      statres.masse[nbphase]=massetotale;  
      statres.volume[nbphase]=volumetotal;  
      for (int i=0;i<nbphase;i++)
        {
        char mess[500];
        sprintf(mess,"   Phase %d (%d cristaux) =  tirage %d \n                            masse %lf ou %lf%% \n                            volume %lf ou %lf%%", i+1,nbcristaux[i],ordrecible[i]+1,tirage[i],tirage[i]/massetotale*100.,volume[i],volume[i]/volumetotal*100.);
        affiche(mess);
        sprintf(mess,"                               Centre de masse  x = %lf ",cdm[i].get_x()/tirage[i]);
        affiche(mess);
        sprintf(mess,"                                                y = %lf ",cdm[i].get_y()/tirage[i]);
        affiche(mess);
        sprintf(mess,"                                                z = %lf ",cdm[i].get_z()/tirage[i]);
        affiche(mess);
        statres.masse[i]=tirage[i];
        statres.volume[i]=volume[i];
        statres.cdm[i]=cdm[i]/tirage[i];
        statres.cdm[nbphase]=statres.cdm[nbphase]+cdm[i];
        statres.nbcristauxphase[i]=nbcristaux[i];
        }
     statres.cdm[statres.nbphase]=(1./statres.masse[statres.nbphase])*statres.cdm[statres.nbphase];
            sprintf(mess,"                               Centre de masse  global x = %lf ",statres.cdm[statres.nbphase].get_x());
        affiche(mess);
        sprintf(mess,"                                                       y = %lf ",statres.cdm[statres.nbphase].get_y());
        affiche(mess);
        sprintf(mess,"                                                       z = %lf ",statres.cdm[statres.nbphase].get_z());
        affiche(mess);
       }
          else 
            {
            affiche((char*)"Reprise conditions aux limites");
            LISTE_MG_VOLUME::iterator it;
        statres.nbphase=(int)param.get_valeur("Nbphase");
        statres.volume.resize(statres.nbphase+1);
        statres.masse.resize(statres.nbphase+1);
        statres.nbcristauxphase.resize(statres.nbphase+1);
        statres.cdm.resize(statres.nbphase+1);
        for (int i=0;i<statres.nbphase;i++)
            {
            statres.volume[i]=0.;
            statres.masse[i]=0.;
            statres.nbcristauxphase[i]=0;
            statres.cdm[i].change_x(0.);
            statres.cdm[i].change_y(0.);
            statres.cdm[i].change_z(0.);
            }
            for (MG_VOLUME* vol=mggeo->get_premier_volume(it);vol!=NULL;vol=mggeo->get_suivant_volume(it))
              {
              double a1,a2,a3;
              vol->get_valeur_ccf((char*)"a1",a1);
              vol->get_valeur_ccf((char*)"a2",a2);
              vol->get_valeur_ccf((char*)"a3",a3);
              double x1[3],x2[3],x3[3];
              double psi=a1/180.*M_PI;
              double teta=a2/180.*M_PI;
              double phi=a3/180.*M_PI;
              x1[0]=cos(phi)*cos(psi)-sin(phi)*cos(teta)*sin(psi);
              x1[1]=cos(phi)*sin(psi)+sin(phi)*cos(teta)*cos(psi);
              x1[2]=sin(phi)*sin(teta);
              x2[0]=-sin(phi)*cos(psi)-cos(phi)*cos(teta)*sin(psi);
              x2[1]=-sin(phi)*sin(psi)+cos(phi)*cos(teta)*cos(psi);
              x2[2]=cos(phi)*sin(teta);
              x3[0]=sin(teta)*sin(psi);
              x3[1]=-sin(teta)*cos(psi);
              x3[2]=cos(teta);
              for (int i=0;i<3;i++)
                for (int j=0;j<3;j++)
                  {
                  tens1[i*3+j]+=x1[i]*x1[j];
                  tens2[i*3+j]+=x2[i]*x2[j];
                  tens3[i*3+j]+=x3[i]*x3[j];
                  }
              double val;
              vol->get_valeur_ccf((char*)"Ph",val);
              int numphase=(int)val;
          vol->get_valeur_ccf((char*)"Ro",val);
          double dens=val;
              double volu,mass;
          OT_VECTEUR_3D cdm2;
          OT_MATRICE_3D m1,m2;
          vol->get_propriete_massique(gest->get_mg_maillage(0),volu,mass,cdm2,m1,m2,dens,-1);    
          statres.volume[numphase]=statres.volume[numphase]+volu;
          statres.masse[numphase]=statres.masse[numphase]+mass;
          statres.nbcristauxphase[numphase]++;
          statres.cdm[numphase]=statres.cdm[numphase]+mass*cdm2;
              }
              statres.masse[statres.nbphase]=0.;
          statres.volume[statres.nbphase]=0.;
          statres.cdm[statres.nbphase].change_x(0.);
          statres.cdm[statres.nbphase].change_y(0.);
          statres.cdm[statres.nbphase].change_z(0.);
          for (int i=0;i<statres.nbphase;i++)
            {
            statres.volume[statres.nbphase]=statres.volume[statres.nbphase]+statres.volume[i];
            statres.masse[statres.nbphase]=statres.masse[statres.nbphase]+statres.masse[i];
            statres.cdm[i]=(1./statres.masse[i])*statres.cdm[i];
            statres.cdm[statres.nbphase]=statres.cdm[statres.nbphase]+statres.masse[i]*statres.cdm[i];
            }
         statres.cdm[statres.nbphase]=(1./statres.masse[statres.nbphase])*statres.cdm[statres.nbphase];
         for (int i=0;i<statres.nbphase;i++)
                {
                char mess[500];
                sprintf(mess,"   Phase %d  (%d cristaux) :    masse %lf ou %lf%% \n                               volume %lf ou %lf%%", i+1,statres.nbcristauxphase[i],statres.masse[i],statres.masse[i]/statres.masse[statres.nbphase]*100.,statres.volume[i],statres.volume[i]/statres.volume[statres.nbphase]*100.);
            
                affiche(mess);
         
         sprintf(mess,"                               Centre de masse  x = %lf ",statres.cdm[i].get_x());
        affiche(mess);
        sprintf(mess,"                                                y = %lf ",statres.cdm[i].get_y());
        affiche(mess);
        sprintf(mess,"                                                z = %lf ",statres.cdm[i].get_z());
        affiche(mess);
        
                }
        sprintf(mess,"                               Centre de masse  global x = %lf ",statres.cdm[statres.nbphase].get_x());
        affiche(mess);
        sprintf(mess,"                                                       y = %lf ",statres.cdm[statres.nbphase].get_y());
        affiche(mess);
        sprintf(mess,"                                                       z = %lf ",statres.cdm[statres.nbphase].get_z());
        affiche(mess);
        
            }
          
          for (int i=0;i<9;i++)
                {
                statres.tens1[i]=tens1[i]/mggeo->get_nb_mg_volume();  
                statres.tens2[i]=tens2[i]/mggeo->get_nb_mg_volume();   
                statres.tens3[i]=tens3[i]/mggeo->get_nb_mg_volume();   
                }
    }
        int degre=(int)param.get_valeur("Degre");
    statres.degre=degre;
        FEM_MAILLAGE* fem;
        if (afaire>1) 
        {
 
        if (reprendmai==0)
        {
          gest->supprimer_tout_fem_solution();
          gest->supprimer_tout_fem_maillage();
          for (int i=gest->get_nb_mg_maillage()-1;i>0;i--)
                gest->supprimer_mg_maillage(i);
          //FCT_TAILLE_FEM_SOLUTION_GENERATEUR_CONSTANT carte(gest,mggeo,dg,1,1,1,1,(char*)nomfichiercarte.c_str());
          affiche((char*)"Carte de taille");
          
          std::vector<double> lstpoint;
          int nbarete=mggeo->get_nb_mg_arete();
          double bornerafmin=param.get_valeur("bornemin_carte");
      int typeborne=(int)param.get_valeur("type_bornemin");
      if (typeborne==1) bornerafmin=bornerafmin*dg;
      for (int i=0;i<nbarete;i++)
              {
              MG_ARETE* are=mggeo->get_mg_arete(i);
              double val=are->get_longueur(are->get_tmin(),are->get_tmax());
                if (val<epsfusion) std::cout<< " pb" << std::endl;
              if (val<dg) 
                  {
                  double xyz[3];
                  are->get_cosommet1()->get_sommet()->get_point()->evaluer(xyz);
                  lstpoint.push_back(xyz[0]);
                  lstpoint.push_back(xyz[1]);
                  lstpoint.push_back(xyz[2]);
                  lstpoint.push_back(std::max(val,bornerafmin));
                  }
            
                }
      char message[500];
          sprintf(message,"   Écart nodal demande : %lf Écart nodal minimal : %lf",dg,bornerafmin);
          affiche(message);
          sprintf(message,"   %d points de raffinement",(int)(lstpoint.size())/4);
          affiche(message);
          //affiche((char*)"Enregistrement");
          //gest->enregistrer(nommagic);
      double gradient=param.get_valeur("gradient");
      FCT_TAILLE_FEM_SOLUTION_GENERATEUR_GRADIENT carte(gest,mggeo,dg,2,50,50,50,(char*)nomfichiercarte.c_str(),lstpoint,gradient);
          //carte.active_affichage(affiche2);
      carte.construit();
          carte.enregistrer(nommagic);
          
          MG_MAILLAGE *mai=new MG_MAILLAGE(mggeo);
          gest->ajouter_mg_maillage(mai);
          affiche((char*)"Maillage 0D");
          MAILLEUR0D m0d(mai,mggeo);
          m0d.active_affichage(affiche2);
          m0d.maille();
          affiche((char*)"Maillage 1D");
          MAILLEUR1D m1d(mai,mggeo,&carte);
          //m1d.active_affichage(affiche2);
          m1d.maille();
          affiche((char*)"Maillage 2D");
          MAILLEUR2D m2d(mai,mggeo,&carte);
          //m2d.active_affichage(affiche2);
          m2d.maille();
          affiche((char*)"Maillage 3D");
          MAILLEUR3D m3d(mai,mggeo,&carte,false);
          //m3d.active_affichage(affiche2);
          LISTE_MG_VOLUME::iterator it;int ii=0;
          for (MG_VOLUME* vol=mggeo->get_premier_volume(it);vol!=NULL;vol=mggeo->get_suivant_volume(it))
                        {
                        ii++;
                        char mess[5000];
                        if (ii!=1) sprintf(mess,"\033[A\033[K    Polycristal %d / %d soit %.2lf%%",ii,mggeo->get_nb_mg_volume(),ii*100.0/mggeo->get_nb_mg_volume());
            else sprintf(mess,"    Polycristal %d / %d soit %.2lf%%",ii,mggeo->get_nb_mg_volume(),ii*100.0/mggeo->get_nb_mg_volume());
            affiche(mess);
                        m3d.maille(vol);
                        }
      //MAILLEUR_ANALYSE mana(mai);
      //mana.active_affichage(affiche2);
      //mana.analyse();
      //optimise_maillage(mai);
          //mana.analyse();
      fem=new FEM_MAILLAGE(mggeo,mai,degre);
          gest->ajouter_fem_maillage(fem);
      char mess[400];
      sprintf(mess,"Maillage FEM degre %d",degre);
      affiche(mess);
          MAILLEUR_FEM mfem;
          mfem.maille(fem,0);
        }
        else
          {
          affiche((char*)"Reprise de calcul d'un ancien maillage"); 
          fem=gest->get_fem_maillage(1);
          }
          
          //verification volume
          double volume=0.;
          LISTE_FEM_ELEMENT3::iterator it;
          for (FEM_ELEMENT3* ele=fem->get_premier_element3(it);ele!=NULL;ele=fem->get_suivant_element3(it))
            {
            MG_TETRA* tet=(MG_TETRA*)ele->get_mg_element_maillage();
            MG_NOEUD* n1=tet->get_noeud1();
            MG_NOEUD* n2=tet->get_noeud2();
            MG_NOEUD* n3=tet->get_noeud3();
            MG_NOEUD* n4=tet->get_noeud4();
            OT_VECTEUR_3D vec1(n1->get_coord(),n2->get_coord());
            OT_VECTEUR_3D vec2(n1->get_coord(),n3->get_coord());
            OT_VECTEUR_3D vec3(n1->get_coord(),n4->get_coord());
            OT_VECTEUR_3D pvec=vec1&vec2;
            double vol=fabs(pvec*vec3);
            volume=volume+vol;
            }
          volume=volume/6.;
          char message[500];
          sprintf(message,"   Volume des mailles = %lf",volume);
          affiche(message);
          
          /*
          MAILLEUR_ANALYSE mana(fem);
          mana.active_affichage(affiche2);
          mana.analyse(NULL);
          */}
          if (afaire>2)
          {
          double epsx,epsy,epsz,epsxy,epsxz,epsyz;
          double sigx,sigy,sigz,sigxy,sigxz,sigyz;
          affiche((char*)"\n---> 1. CLDHS\n");
          MSTRUCT_OUTILS::appliquer_conditions_limites(mggeo,boite,MSTRUCT::TYPE_ETUDE::MECANIQUE,MSTRUCT::TYPE_CONDITIONS_LIMITES::DEFORMATION_HOMOGENE,MSTRUCT::TYPE_CHARGEMENT::SPHERIQUE,0.0001,0,0.0001);
          int erreur=calcule_cacteristique_mecanique(gest,fem,epsx,epsy,epsz,epsxy,epsxz,epsyz,sigx,sigy,sigz,sigxy,sigxz,sigyz);
          if (erreur==0)
            {
            statres.KDH=1./3.*(sigx+sigy+sigz)/(epsx+epsy+epsz);
            mggeo->efface_ccf(false,true,true,true);
            char message[500];
        sprintf(message,"          KDH=%lf",statres.KDH);
            affiche(message);
            }
        else exit(EXIT_FAILURE);
          
          if (erreur==0)
          {
          
          affiche((char*)"\n---> 2. CLDHD\n");
          MSTRUCT_OUTILS::appliquer_conditions_limites(mggeo,boite,MSTRUCT::TYPE_ETUDE::MECANIQUE,MSTRUCT::TYPE_CONDITIONS_LIMITES::DEFORMATION_HOMOGENE,MSTRUCT::TYPE_CHARGEMENT::DEVIATORIQUE,0.0001,0,0.0001);
          erreur=calcule_cacteristique_mecanique(gest,fem,epsx,epsy,epsz,epsxy,epsxz,epsyz,sigx,sigy,sigz,sigxy,sigxz,sigyz);
          if (erreur==0)
            {
            statres.GDH=1./6.*(sigxy/epsxy+sigyz/epsyz+sigxz/epsxz);
            statres.EDH=9*statres.KDH*statres.GDH/(3*statres.KDH+statres.GDH);
            mggeo->efface_ccf(false,true,true,true);
            char message[500];
            sprintf(message,"          GDH=%lf\n          EDH=%lf\n",statres.GDH,statres.EDH);
            affiche(message);
            }
        else exit(EXIT_FAILURE);
          }
          
          if (erreur==0)
          {
          affiche((char*)"\n---> 3. CLCHS\n");
          MSTRUCT_OUTILS::appliquer_conditions_limites(mggeo,boite,MSTRUCT::TYPE_ETUDE::MECANIQUE,MSTRUCT::TYPE_CONDITIONS_LIMITES::CONTRAINTE_HOMOGENE,MSTRUCT::TYPE_CHARGEMENT::SPHERIQUE,0.0001,0,0.0001);
          erreur=calcule_cacteristique_mecanique(gest,fem,epsx,epsy,epsz,epsxy,epsxz,epsyz,sigx,sigy,sigz,sigxy,sigxz,sigyz);
          if (erreur==0)
            {
            statres.KCH=1./3.*(sigx+sigy+sigz)/(epsx+epsy+epsz);
            mggeo->efface_ccf(false,true,true,true);
            char message[500];
            sprintf(message,"          KCH=%lf\n\n",statres.KCH);
            affiche(message);
            }
       else exit(EXIT_FAILURE);
 
          }
          
          if (erreur==0)
          {
          affiche((char*)"\n---> 4. CLCHD\n");
          MSTRUCT_OUTILS::appliquer_conditions_limites(mggeo,boite,MSTRUCT::TYPE_ETUDE::MECANIQUE,MSTRUCT::TYPE_CONDITIONS_LIMITES::CONTRAINTE_HOMOGENE,MSTRUCT::TYPE_CHARGEMENT::DEVIATORIQUE,0.0001,0,0.0001);
          calcule_cacteristique_mecanique(gest,fem,epsx,epsy,epsz,epsxy,epsxz,epsyz,sigx,sigy,sigz,sigxy,sigxz,sigyz);
          if (erreur==0)
              {
              statres.GCH=1./6.*(sigxy/epsxy+sigyz/epsyz+sigxz/epsxz);
              statres.ECH=9*statres.KCH*statres.GCH/(3*statres.KCH+statres.GCH);
              mggeo->efface_ccf(false,true,true,true);
          char message[500];
              sprintf(message,"          GCH=%lf\n          ECH=%lf\n",statres.GCH,statres.ECH);
              affiche(message);
              }
      else exit(EXIT_FAILURE);
 
          //double GDH,EDH,KCH,GCH,ECH;
          affiche((char*)"\n\n Résultats \n\n");
      char message[500];
          sprintf(message,"     KDH=%lf;GDH=%lf;EDH=%lf;",statres.KDH,statres.GDH,statres.EDH);
          affiche(message);
          sprintf(message,"     KCH=%lf;GCH=%lf;ECH=%lf;\n",statres.KCH,statres.GCH,statres.ECH);
          affiche(message);
          gest->enregistrer(nommagic);
          
      
      
      }  
     
        }
        int enregistreresultat=(int)param.get_valeur("enregistrerres");
    if (enregistreresultat==1)
        {
        /*FILE *outborne=fopen(nomfichierborne.c_str(),"rb");
                fread(&statborne,sizeof(statborne),1,outborne);
                fclose(outborne);*/
          calcule_borne(statborne,statres);
          std::mutex mtx; 
          std::string nomcumul=param.get_nom("fichiercumul");
          std::string nomcumul1 = nomcumul+".txt";
          std::string nomcumul2 = nomcumul+".gnu";
          std::string nomcumul3 = nomcumul+".svg";
          FILE *outtest=fopen((char*)nomcumul1.c_str(),"rt");
          bool existe =false;
          if (outtest!=NULL) {existe=true;fclose(outtest);}
          mtx.lock();
          FILE *out=fopen((char*)nomcumul1.c_str(),"at");
          if (!existe) 
      { fprintf(out,"Nbcristaux;degre;dg;KDH;GDH;EDH;KCH;GCH;ECH;ori1_11;ori1_12;ori1_13;ori1_21;ori1_22;ori1_23;ori1_31;ori1_32;ori1_33;ori2_11;ori2_12;ori2_13;ori2_21;ori2_22;ori2_23;ori2_31;ori2_32;ori2_33;ori3_11;ori3_12;ori3_13;ori3_21;ori3_22;ori3_23;ori3_31;ori3_32;ori3_33;Nbphase;");
      for (int i=0;i<statres.nbphase;i++)
          fprintf(out,"N%d;VPh%d;MPh%d;CDMx%d;CDMy%d;CDMz%d;",i+1,i+1,i+1,i+1,i+1,i+1);
      fprintf(out,"VPhtot;MPhtot;CDMxtot;CDMytot;CDMztot,Kv,Gv,Ev,nuv,Kr,Gr,Er,nur\n");
      }
          fprintf(out,"%d;%d;%lf;",statres.nbcristaux,statres.degre,statres.dg);
          fprintf(out,"%lf;%lf;%lf;",statres.KDH,statres.GDH,statres.EDH);
          fprintf(out,"%lf;%lf;%lf;",statres.KCH,statres.GCH,statres.ECH);
          for (int i=0;i<9;i++)
            fprintf(out,"%lf;",statres.tens1[i]);
          for (int i=0;i<9;i++)
            fprintf(out,"%lf;",statres.tens2[i]);
          for (int i=0;i<9;i++)
            fprintf(out,"%lf;",statres.tens3[i]);
          fprintf(out,"%d;",statres.nbphase);
      for (int i=0;i<statres.nbphase;i++)
          fprintf(out,"%d;%lf;%lf;%lf;%lf;%lf;",statres.nbcristauxphase[i],statres.volume[i],statres.masse[i],statres.cdm[i].get_x(),statres.cdm[i].get_y(),statres.cdm[i].get_z());
      fprintf(out,"%lf;%lf;%lf;%lf;%lf;",statres.volume[statres.nbphase],statres.masse[statres.nbphase],statres.cdm[statres.nbphase].get_x(),statres.cdm[statres.nbphase].get_y(),statres.cdm[statres.nbphase].get_z());
      fprintf(out,"%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf;",statborne.Kvoigt,statborne.Gvoigt,statborne.Evoigt,statborne.nuvoigt,statborne.Kreuss,statborne.Greuss,statborne.Ereuss,statborne.nureuss);
      fprintf(out,"\n");
          
        fclose(out);
        mtx.unlock();
            /*outtest=fopen((char*)nomcumul2.c_str(),"rt");
            existe =false;
              if (outtest!=NULL) {existe=true;fclose(outtest);}
              if (!existe) 
                  {
                  mtx.lock();
              FILE *out=fopen((char*)nomcumul2.c_str(),"wt");  
                  fprintf(out, "set terminal svg\n");
                  fprintf(out, "set output \"%s\"\n", nomcumul3.c_str());
                  fprintf(out, "set datafile separator \";\"\n");
                  fprintf(out, "set key on  bmargin horizontal\n");
                  fprintf(out, "set xlabel 'Numéro d''essai'\n");
                  fprintf(out, "set ylabel 'Module d''Young (GPa)'\n");
                  fprintf(out, "set y2label 'Écart à l''isotropie'\n");
                  fprintf(out, "set pointsize 1\n");
                  fprintf(out, "set autoscale x\n");
                  fprintf(out, "set autoscale y\n");
                  fprintf(out, "set y2tics 0.5,0.1,1\n");
                  fprintf(out, "set y2r [0.5:1]\n");
                  fprintf(out, "set ytics nomirror\n");
                  fprintf(out, "set xtics nomirror\n");
                  double facteur = (param.get_valeur("bornevoigt")-param.get_valeur("bornereuss"))/3./1e9;
                  fprintf(out, "set offsets graph 0, 0, %lf,0\n", facteur);
                  fprintf(out, "#\n");
                  fprintf(out, "d(x)=dd\n");
                  fprintf(out, "fit d(x) '%s' u :($6/1e9) via dd\n",nomcumul1.c_str());
                  fprintf(out, "c(x)=cc\n");
                  fprintf(out, "fit c(x) '%s' u :($9/1e9) via cc\n",nomcumul1.c_str());
                  fprintf(out, "#\n");
                  fprintf(out, "set title \"%s\\n %d cristaux\"\n",param.get_nom("nommateriau").c_str(), (int)param.get_valeur("nb_cristaux"));
                  fprintf(out, "plot '%s' u :($6/1e9) title 'EDH' pointtype 5 lc rgb \"blue\", '%s' u :($9/1e9) title 'ECH' pointtype 5 lc rgb \"red\",(%lf/1e9) title 'Borne Voigt' w l lc rgb \"violet\" lw 3,(%lf/1e9) title 'Borne Reuss' w l lc rgb \"dark-orange\" lw 3,  '%s'  u :(1.-sqrt(1./9.*(($10-1./3)*($10-1./3)+($14-1./3.)*($14-1./3)+($18-1./3.)*($18-1./3.)+($19-1./3)*($19-1./3)+($23-1./3.)*($23-1./3)+($27-1./3.)*($27-1./3.)+($28-1./3)*($28-1./3)+($32-1./3.)*($32-1./3)+($36-1./3.)*($36-1./3.)))) axis x1y2 title 'Isotropie' lc rgb \"green\" pointtype 1,c(x) title sprintf(\"moyenne %%.1f\",cc) lc rgb \"red\" dt 2 ,d(x) title sprintf(\"moyenne %%.1f\",dd) lc rgb \"blue\" dt 2", nomcumul1.c_str(),nomcumul1.c_str(),statres.Evoigt,statres.Ereuss,nomcumul1.c_str() );
                  
                  
                  fclose(out);
          mtx.unlock();
                  }*/
          }
          
    
        
        affiche((char*)"Enregistrement");
    gest->enregistrer(nommagic);
    delete gest;
    
  #endif 
  }
 
 
 
  int MSTRUCT_GENERATEUR_POLYCRISTAUX::calcule_cacteristique_mecanique(MG_GESTIONNAIRE *gest,FEM_MAILLAGE* fem, double& epsx, double& epsy, double& epsz, double& epsxy, double& epsxz, double& epsyz, double& sigx, double& sigy, double& sigz, double& sigxy, double& sigxz, double& sigyz)
  {
 
    //affiche((char*)"Calcul aster");
    MGASTER ef;
    ef.active_affichage(affiche2);
    std::string nomparamaster=param.get_nom("paramaster");
    std::string nomfichier=param.get_nom("nomfichier");
    int erreur=ef.calcule((char*)nomparamaster.c_str(),fem,(char*)nomfichier.c_str(),MAGIC::CALCUL_ASTER::ELASTIQUE,(char*)"00011111",false);
    if (erreur!=0) return erreur;
    int nb=gest->get_nb_fem_solution();
    for (int i=0;i<nb;i++)
    {
    
    FEM_SOLUTION* sol=gest->get_fem_solution(i);
    std::string nom=sol->get_nom();
    nom.erase(2,1000);
    if (nom=="GE")
          {
          epsx=sol->get_moyenne_volumique_champs(NULL,0); 
          epsy=sol->get_moyenne_volumique_champs(NULL,1); 
          epsz=sol->get_moyenne_volumique_champs(NULL,2); 
          epsxy=sol->get_moyenne_volumique_champs(NULL,3); 
          epsxz=sol->get_moyenne_volumique_champs(NULL,4); 
          epsyz=sol->get_moyenne_volumique_champs(NULL,5); 
          }
    if ((nom=="GS") && (sol->get_nb_champ()==6))
          {
          sigx=sol->get_moyenne_volumique_champs(NULL,0); 
          sigy=sol->get_moyenne_volumique_champs(NULL,1); 
          sigz=sol->get_moyenne_volumique_champs(NULL,2); 
          sigxy=sol->get_moyenne_volumique_champs(NULL,3); 
          sigxz=sol->get_moyenne_volumique_champs(NULL,4); 
          sigyz=sol->get_moyenne_volumique_champs(NULL,5); 
          }
    }
    return erreur;
  }
  
  
  
  
  
  
  void MSTRUCT_GENERATEUR_POLYCRISTAUX::optimise_maillage(MG_MAILLAGE* mai)
{  
double bornerafmin=param.get_valeur("bornemin_carte");
int typeborne=(int)param.get_valeur("type_bornemin"); 
double dg=param.get_valeur("dg");
if (typeborne==1) bornerafmin=bornerafmin*dg;
TPL_MAP_ENTITE<MG_TETRA*> coq;
std::multimap<MG_NOEUD*,MG_NOEUD*> lst1;
std::map<MG_NOEUD*,MG_NOEUD*> lst2;
LISTE_MG_SEGMENT::iterator it;
for (MG_SEGMENT* seg=mai->get_premier_segment(it);seg!=NULL;seg=mai->get_suivant_segment(it))
    {
    if (seg->get_longueur()<bornerafmin)     
        {
        int nb1=seg->get_noeud1()->get_lien_tetra()->get_nb();
        int nb2=seg->get_noeud2()->get_lien_tetra()->get_nb();
        MG_NOEUD* no1=seg->get_noeud1();
        MG_NOEUD* no2=seg->get_noeud2();
        if (no1>no2)
            {
            MG_NOEUD* tmp=no1;
            no1=no2;
            no2=tmp;
            }
        std::map<MG_NOEUD*,MG_NOEUD*>::iterator itn=lst2.find(no1);
        if (itn!=lst2.end())
            no1=itn->second;    
        std::pair<MG_NOEUD*,MG_NOEUD*> tmp1(no1,no2);
        std::pair<MG_NOEUD*,MG_NOEUD*> tmp2(no2,no1);
        std::pair<std::map<MG_NOEUD*,MG_NOEUD*>::iterator,bool> ret=lst2.insert(tmp2);
        if (ret.second==false) continue;
        lst1.insert(tmp1);
        for (int i=0;i<nb1;i++)
            for (int j=0;j<nb2;j++)
                if (seg->get_noeud1()->get_lien_tetra()->get(i)==seg->get_noeud2()->get_lien_tetra()->get(j))
                    {
                    MG_TETRA* tet=seg->get_noeud1()->get_lien_tetra()->get(i);
                    coq.ajouter(tet);
                    }
        }
    }
TPL_MAP_ENTITE<MG_TETRA*>::ITERATEUR it2;
for (MG_TETRA* tet=coq.get_premier(it2);tet!=NULL;tet=coq.get_suivant(it2))
    mai->supprimer_mg_tetraid(tet->get_id());
for (std::multimap<MG_NOEUD*,MG_NOEUD*>::iterator itn=lst1.begin();itn!=lst1.end();itn++)
    {
    TPL_LISTE_ENTITE<MG_NOEUD*> lsttetaconstruire;
    TPL_LISTE_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lsttetaconstruire2;
    coq.vide();
    MG_NOEUD* no1=itn->first;    
    MG_NOEUD* no2=itn->second;    
    int nb=no2->get_lien_tetra()->get_nb();
    for (int i=0;i<nb;i++)
        {
        MG_TETRA* tet=no2->get_lien_tetra()->get(i);
        if (tet->get_noeud1()==no2) lsttetaconstruire.ajouter(no1); else lsttetaconstruire.ajouter(tet->get_noeud1());
        if (tet->get_noeud2()==no2) lsttetaconstruire.ajouter(no1); else lsttetaconstruire.ajouter(tet->get_noeud2());
        if (tet->get_noeud3()==no2) lsttetaconstruire.ajouter(no1); else lsttetaconstruire.ajouter(tet->get_noeud3());
        if (tet->get_noeud4()==no2) lsttetaconstruire.ajouter(no1); else lsttetaconstruire.ajouter(tet->get_noeud4());
        lsttetaconstruire2.ajouter(tet->get_lien_topologie());
        coq.ajouter(tet);
        }
       
for (int i=0;i<lsttetaconstruire.get_nb()/4;i++)
    {
    double qual=OPERATEUR::qualite_tetra(lsttetaconstruire.get(4*i)->get_coord(),lsttetaconstruire.get(4*i+1)->get_coord(),lsttetaconstruire.get(4*i+2)->get_coord(),lsttetaconstruire.get(4*i+3)->get_coord());
    if (qual<1e-14) std::cout << qual << std::endl;
        mai->ajouter_mg_tetra(lsttetaconstruire2.get(i),lsttetaconstruire.get(4*i),lsttetaconstruire.get(4*i+1),lsttetaconstruire.get(4*i+2),lsttetaconstruire.get(4*i+3),MAGIC::ORIGINE::MAILLEUR_AUTO);
    }
        
    }
for (MG_TETRA* tet=coq.get_premier(it2);tet!=NULL;tet=coq.get_suivant(it2))
    mai->supprimer_mg_tetraid(tet->get_id());
                        
    
    
    
    
}
 
 void MSTRUCT_GENERATEUR_POLYCRISTAUX::calcule_borne(MSTRUCT_GENERATEUR_POLYCRISTAUX_BORNE& statborne,MSTRUCT_GENERATEUR_POLYCRISTAUX_RESULTAT &echantillon)
{
  OT_TENSEUR rigtot(6,6),soutot(6,6);
  MSTRUCT_BORNE  borne;
          for (int i=0;i<echantillon.nbphase;i++)
        {
        double E1=param.get_valeur("E1",i);
                double E2=param.get_valeur("E2",i);
                double E3=param.get_valeur("E3",i);
                double n1=param.get_valeur("n1",i);
                double n2=param.get_valeur("n2",i);
                double n3=param.get_valeur("n3",i);
                double G1=param.get_valeur("G1",i);
                double G2=param.get_valeur("G2",i);
                double G3=param.get_valeur("G3",i);
                OT_TENSEUR rig(6,6),rigphase(6,6),souphase(6,6);
                borne.get_tenseur_rigidite(E1/1e9,E2/1e9,E3/1e9,2.*G1/1e9,2.*G2/1e9,2.*G3/1e9,n1,n2,n3,rig);
                borne.get_VoigtReuss_poly(rig,rigphase,souphase);
                double proportion=echantillon.masse[i]/echantillon.masse[echantillon.nbphase];
                rigphase=proportion*rigphase;
        souphase=proportion*souphase;
                rigtot=rigtot+rigphase;
        soutot=soutot+souphase;
        }
      double Ev,Gv,Kv,nuv;  
          double Er,Gr,Kr,nur;  
          borne.get_VoigtReuss_poly(rigtot,soutot,Kv,Gv,Ev,nuv,Kr,Gr,Er,nur);
          statborne.Evoigt=Ev*1e9;
      statborne.Kvoigt=Kv/3.*1e9;
          statborne.Gvoigt=Gv/2.*1e9;
          statborne.nuvoigt=nuv;
          statborne.Ereuss=Er*1e9;
          statborne.Kreuss=Kr/3.*1e9;
          statborne.Greuss=Gr/2.*1e9;
          statborne.nureuss=nur;
        
}