mgopt_simp_adaptative.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/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mgopt_simp_adaptative.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:57 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
#include "gestionversion.h"
#include "ot_mathematique.h"
#include "mg_file.h"
#include "mg_export.h"
#include "fem_solution.h"
#include "fem_maillage.h"
#include "mg_gestionnaire.h"
#include "mgopt_simp_adaptative.h"
#include "mgopt_simp.h"
#include "ot_cpu.h"
#include <fct_generateur_fichier.h>
#include <fct_taille_fem_solution_fonction_isotrope.h>
#include <tpl_octree.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <fstream>
 
void affiche(char* message)
{
    std::cout << message << std::endl;
}
 
 
 
MGOPT_SIMP_ADAPTATIVE::MGOPT_SIMP_ADAPTATIVE(bool save):MGOPT_SIMP(save)
{
params.ajouter("carte_taille_structure_adapt",20.,OT_PARAMETRES::DOUBLE,"Utilisation d'un support de carte de taille structure 0-non nombre=nombre de subdivision en une direction");
params.ajouter("iter_adaptative_max",10.,OT_PARAMETRES::DOUBLE,"Nombre maximal d'itérations de la SIMP adaptative");
params.ajouter("seuil_convergence_adaptative",0.1,OT_PARAMETRES::DOUBLE,"Critere de convergence minimal dans les iterations de SIMP lors d'adaptation");
params.ajouter("enm",7.,OT_PARAMETRES::DOUBLE,"Ecart nodal maximal de la carte de taille de la SIMP adaptative");
params.ajouter("khi",0.,OT_PARAMETRES::DOUBLE,"Parametre de controle de la vitesse de raffinement toujours > 0 (si 0. maillage constant de taille enm)");
params.ajouter("coef_critere_convergence",1.,OT_PARAMETRES::DOUBLE,"Coefficient de multiplication du critere_convergence apres chaque cycle");
params.ajouter("seuil_lissage_gradient",0.3,OT_PARAMETRES::DOUBLE,"Seuil de lissage du gradient");
}
 
MGOPT_SIMP_ADAPTATIVE::MGOPT_SIMP_ADAPTATIVE(MGOPT_SIMP_ADAPTATIVE &mdd):MGOPT_SIMP(mdd)
{
}
 
MGOPT_SIMP_ADAPTATIVE::~MGOPT_SIMP_ADAPTATIVE()
{
}
 
 
double MGOPT_SIMP_ADAPTATIVE::get_gradient_element(FEM_ELEMENT3* tet)
{ 
  double jac[9],j[9];
  if ((tet->get_nb_fem_noeud()==4))
        {
        FEM_NOEUD* n1=tet->get_fem_noeud(0);
        FEM_NOEUD* n2=tet->get_fem_noeud(1);
        FEM_NOEUD* n3=tet->get_fem_noeud(2);
        FEM_NOEUD* n4=tet->get_fem_noeud(3);
        double *xyz1=n1->get_coord();
        double *xyz2=n2->get_coord();
        double *xyz3=n3->get_coord();
        double *xyz4=n4->get_coord();
        jac[0]=(xyz2[0]-xyz1[0]);
        jac[1]=(xyz2[1]-xyz1[1]);
        jac[2]=(xyz2[2]-xyz1[2]);
        jac[3]=(xyz3[0]-xyz1[0]);
        jac[4]=(xyz3[1]-xyz1[1]);
        jac[5]=(xyz3[2]-xyz1[2]);
        jac[6]=(xyz4[0]-xyz1[0]);
        jac[7]=(xyz4[1]-xyz1[1]);
        jac[8]=(xyz4[2]-xyz1[2]);
        double detj=jac[0]*jac[4]*jac[8]+jac[3]*jac[7]*jac[2]+jac[6]*jac[1]*jac[5]-jac[2]*jac[4]*jac[6]-jac[5]*jac[7]*jac[0]-jac[8]*jac[1]*jac[3];
 
        j[0*3+0] =  (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
        j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
        j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
        j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
        j[1*3+1] =  (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
        j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
        j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
        j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
        j[2*3+2] =  (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
 
        //Gradient par element Design
        double d0 = n1->get_solution();
        double d1 = n2->get_solution();
        double d2 = n3->get_solution();
        double d3 = n4->get_solution();
 
        double gd_x = j[0*3+0]*(d1-d0)+j[0*3+1]*(d2-d0)+j[0*3+2]*(d3-d0);
        double gd_y = j[1*3+0]*(d1-d0)+j[1*3+1]*(d2-d0)+j[1*3+2]*(d3-d0);
        double gd_z = j[2*3+0]*(d1-d0)+j[2*3+1]*(d2-d0)+j[2*3+2]*(d3-d0);
 
        OT_VECTEUR_3D gd(gd_x,gd_y,gd_z);
        return(gd.get_longueur());
        }
  return 0.;
}
 
 
 
void MGOPT_SIMP_ADAPTATIVE::gradient_noeud(FEM_MAILLAGE* femmai,MG_GESTIONNAIRE* gest,char* nom_solution_gradient)
{
  double grdmax = 0.001;
  double grdmax_grd = 0.001;
  TPL_MAP_ENTITE<FEM_NOEUD*> lstnode;
  TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR in;
  TPL_MAP_ENTITE<FEM_NOEUD*> liste_noeuds_design;
  TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR itd;
  TPL_MAP_ENTITE<FEM_NOEUD*> liste_noeuds_nondesign;
  TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR itnd;
  TPL_MAP_ENTITE<FEM_ELEMENT3*> liste_tetras_design;
  TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR ittd;
  TPL_MAP_ENTITE<FEM_ELEMENT3*> liste_tetras_nondesign;
  TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR ittnd;
  
  affiche((char*)"Extrapolation de la densité aux noeuds ");
  LISTE_FEM_NOEUD::iterator it;
  for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
        no->change_solution(0.);
  
  LISTE_FEM_ELEMENT3::iterator itt;
  for (FEM_ELEMENT3* tet=femmai->get_premier_element3(itt);tet!=NULL;tet=femmai->get_suivant_element3(itt))
        {
        if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
                {
                FEM_NOEUD* n1=tet->get_fem_noeud(0);
                FEM_NOEUD* n2=tet->get_fem_noeud(1);
                FEM_NOEUD* n3=tet->get_fem_noeud(2);
                FEM_NOEUD* n4=tet->get_fem_noeud(3);
                liste_noeuds_design.ajouter(n1);
                liste_noeuds_design.ajouter(n2);
                liste_noeuds_design.ajouter(n3);
                liste_noeuds_design.ajouter(n4);
                liste_tetras_design.ajouter(tet);
                }
        if (tet->get_mg_element_maillage()->get_origine()==MAGIC::ORIGINE::IMPOSE)
                {
                FEM_NOEUD* n1=tet->get_fem_noeud(0);
                FEM_NOEUD* n2=tet->get_fem_noeud(1);
                FEM_NOEUD* n3=tet->get_fem_noeud(2);
                FEM_NOEUD* n4=tet->get_fem_noeud(3);
                liste_noeuds_nondesign.ajouter(n1);
                liste_noeuds_nondesign.ajouter(n2);
                liste_noeuds_nondesign.ajouter(n3);
                liste_noeuds_nondesign.ajouter(n4);
                liste_tetras_nondesign.ajouter(tet);
                }
        }
  TPL_LISTE_ENTITE<FEM_NOEUD*> lstasup;
  for (FEM_NOEUD* no1=liste_noeuds_design.get_premier(itd);no1!=NULL;no1=liste_noeuds_design.get_suivant(itd))
        for (FEM_NOEUD* no2=liste_noeuds_nondesign.get_premier(itnd);no2!=NULL;no2=liste_noeuds_nondesign.get_suivant(itnd))
                if (no1==no2) lstasup.ajouter(no2);
    for (int i=0;i<lstasup.get_nb();i++)
        liste_noeuds_nondesign.supprimer(lstasup.get(i));
 
  for (FEM_NOEUD* no=liste_noeuds_design.get_premier(itd);no!=NULL;no=liste_noeuds_design.get_suivant(itd))
        {
        double nume=0.;
        double deno=0.;
        for (int i=0;i<no->get_lien_element3()->get_nb();i++)
                {
                FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
                if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
                        {
                        double jac[9];
                        double volume=tet->get_jacobien(jac,jac,1.);
                        nume=nume+tet->get_solution()*volume;
                        deno=deno+volume;
                        }
                }
        no->change_solution(nume/deno);
        }
 
  for (FEM_ELEMENT3* tet=liste_tetras_design.get_premier(ittd);tet!=NULL;tet=liste_tetras_design.get_suivant(ittd))
        {
        if ((tet->get_nb_fem_noeud()==4)&&(tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE))
                tet->change_solution(get_gradient_element(tet));
        }
 
  affiche((char*)"Extrapolation du gradient aux noeuds ");
  for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
        no->change_solution(0.);
 
  for (FEM_NOEUD* no=liste_noeuds_design.get_premier(itd);no!=NULL;no=liste_noeuds_design.get_suivant(itd))
        {
        double vect=0.;
        double vol=0.;
        for (int i=0;i<no->get_lien_element3()->get_nb();i++)
                {
                FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
                if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
                        {
                        FEM_NOEUD* no1=tet->get_fem_noeud(0);
                        FEM_NOEUD* no2=tet->get_fem_noeud(1);
                        FEM_NOEUD* no3=tet->get_fem_noeud(2);
                        FEM_NOEUD* no4=tet->get_fem_noeud(3);
                        double *xyz1=no1->get_coord();
                        double *xyz2=no2->get_coord();
                        double *xyz3=no3->get_coord();
                        double *xyz4=no4->get_coord();
                        OT_VECTEUR_3D vec1(xyz1,xyz2);
                        OT_VECTEUR_3D vec2(xyz1,xyz3);
                        OT_VECTEUR_3D vec3(xyz1,xyz4);
                        OT_VECTEUR_3D pvec=vec1&vec2 ;
                        double volume=(fabs(vec3*pvec))/6.;
                        vect=vect+volume*tet->get_solution();
                        vol=vol+volume;
                        }
                }
        if(vect/vol > 0.001) no->change_solution(vect/vol);
        if(vect/vol > grdmax) grdmax = vect/vol;
        if(vect/vol >= 1.e-308) lstnode.ajouter(no);
        }
 
        char message[500];
        sprintf(message,"gradient maximal = %lf",grdmax);
        affiche(message);
 
 
  // Calcul du gradient (du gradient) par element de design
  for (FEM_ELEMENT3* tet=liste_tetras_design.get_premier(ittd);tet!=NULL;tet=liste_tetras_design.get_suivant(ittd))
        tet->change_solution(0.);
  for (FEM_ELEMENT3* tet=liste_tetras_design.get_premier(ittd);tet!=NULL;tet=liste_tetras_design.get_suivant(ittd))
        {
        if ((tet->get_nb_fem_noeud()==4)&&(tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE))
                tet->change_solution(get_gradient_element(tet));
        }
  // Calcul de la variation maximale du gradient aux noeuds
  for (FEM_NOEUD* no=liste_noeuds_design.get_premier(itd);no!=NULL;no=liste_noeuds_design.get_suivant(itd))
        {
        double vect=0.;
        double vol=0.;
        for (int i=0;i<no->get_lien_element3()->get_nb();i++)
                {
                FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
                if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
                        {
                        FEM_NOEUD* no1=tet->get_fem_noeud(0);
                        FEM_NOEUD* no2=tet->get_fem_noeud(1);
                        FEM_NOEUD* no3=tet->get_fem_noeud(2);
                        FEM_NOEUD* no4=tet->get_fem_noeud(3);
                        double *xyz1=no1->get_coord();
                        double *xyz2=no2->get_coord();
                        double *xyz3=no3->get_coord();
                        double *xyz4=no4->get_coord();
                        OT_VECTEUR_3D vec1(xyz1,xyz2);
                        OT_VECTEUR_3D vec2(xyz1,xyz3);
                        OT_VECTEUR_3D vec3(xyz1,xyz4);
                        OT_VECTEUR_3D pvec=vec1&vec2 ;
                        double volume=(fabs(vec3*pvec))/6.;
                        vect=vect+volume*tet->get_solution();
                        vol=vol+volume;
                        }
                }
        if(vect/vol > grdmax_grd) grdmax_grd = vect/vol;
        }
  sprintf(message,"gradient maximal du gradient = %lf",grdmax_grd);
  affiche(message);
  double seuil = params.get_valeur("seuil_lissage_gradient");
  
  // Lissage du gradient a partir de sa variation
  TPL_MAP_ENTITE<FEM_NOEUD*> stnode2;
  for (FEM_NOEUD* no=liste_noeuds_design.get_premier(itd);no!=NULL;no=liste_noeuds_design.get_suivant(itd))
        {
        double vect=0.;
        double vol=0.;
        for (int i=0;i<no->get_lien_element3()->get_nb();i++)
                {
                FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
                if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
                        {
                        FEM_NOEUD* no1=tet->get_fem_noeud(0);
                        FEM_NOEUD* no2=tet->get_fem_noeud(1);
                        FEM_NOEUD* no3=tet->get_fem_noeud(2);
                        FEM_NOEUD* no4=tet->get_fem_noeud(3);
                        double *xyz1=no1->get_coord();
                        double *xyz2=no2->get_coord();
                        double *xyz3=no3->get_coord();
                        double *xyz4=no4->get_coord();
                        OT_VECTEUR_3D vec1(xyz1,xyz2);
                        OT_VECTEUR_3D vec2(xyz1,xyz3);
                        OT_VECTEUR_3D vec3(xyz1,xyz4);
                        OT_VECTEUR_3D pvec=vec1&vec2 ;
                        double volume=(fabs(vec3*pvec))/6.;
                        vect=vect+volume*tet->get_solution();
                        vol=vol+volume;
                        }
                }
        if(vect/vol > seuil*grdmax_grd)
                {
                TPL_MAP_ENTITE<FEM_NOEUD*> stnode;
                TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR stt;
                double som=0.;
                for (int i=0;i<no->get_lien_element3()->get_nb();i++)
                        {
                        FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
                        if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
                                {
                                FEM_NOEUD* no1=tet->get_fem_noeud(0);
                                FEM_NOEUD* no2=tet->get_fem_noeud(1);
                                FEM_NOEUD* no3=tet->get_fem_noeud(2);
                                FEM_NOEUD* no4=tet->get_fem_noeud(3);
                                stnode.ajouter(no1);
                                stnode.ajouter(no2);
                                stnode.ajouter(no3);
                                stnode.ajouter(no4);
                                }
                        }
                stnode.supprimer(no);
                int nb=stnode.get_nb();
                for (FEM_NOEUD* nod=stnode.get_premier(stt);nod!=NULL;nod=stnode.get_suivant(stt))
                        som=som+nod->get_solution();
                no->change_solution(som/nb);
                stnode2.ajouter(no);
/*              for (FEM_NOEUD* nod=stnode.get_premier(stt);nod!=NULL;nod=stnode.get_suivant(stt))
                        stnode.supprimer(nod);*/
                }
        }
  grdmax=1e-308;
  for (FEM_NOEUD* no=liste_noeuds_design.get_premier(itd);no!=NULL;no=liste_noeuds_design.get_suivant(itd))
        if(no->get_solution() > grdmax) grdmax = no->get_solution();
  sprintf(message,"Nombre de noeuds lisses # %d",stnode2.get_nb());
  affiche(message);
 
 
  FEM_SOLUTION* sol3=new FEM_SOLUTION(femmai,1,nom_solution_gradient,femmai->get_nb_fem_noeud(),"Gradient de densite",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
  gest->ajouter_fem_solution(sol3);
  sol3->change_legende(0,"noeuds");
  int i=0;
  for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
        {
        sol3->ecrire(no->get_solution(),i,0);
        i++;
        }
}
 
 
 
 
void MGOPT_SIMP_ADAPTATIVE::adapte(char *nomgestd,char *nomgestnd,char *nomsortie,char *nomparam,char *nomparamaster)
{
affiche((char*)"");
affiche((char*)"************************************");
affiche((char*)"Optimisation adaptative de topologie");
affiche((char*)"************************************");
affiche((char*)"");
 
affiche((char*)"");
affiche((char*)"");  
time_t heuredeb = time(NULL);
tdebut = *localtime(&heuredeb);
if (nomparam!=NULL) lire_params(nomparam);
char *p=strchr(nomgestd,'.');
strncpy(nometude,nomgestd,p-nomgestd);
nometude[p-nomgestd]=0;
p=strchr(nomsortie,'.');
strncpy(nometudesortie,nomsortie,p-nomsortie);
nometudesortie[p-nomsortie]=0;
 
OT_CPU cpu;
cpu.ajouter_etape("debut SIMP Adaptative");
 
affiche((char*)"Itération 1 de SIMP");
affiche((char*)"Preparation du modèle");
gestd=new MG_FILE(nomgestd);
gestnd=new MG_FILE(nomgestnd);
unsigned long maille=(unsigned long)params.get_valeur("maille");
FEM_MAILLAGE* fem;
 
cpu.ajouter_etape("duree preparation modele");
 
if (maille==0) fem=prepare_model(gestd,gestnd);
else fem=gestd->get_fem_maillageid(maille);
 
cpu.ajouter_etape("temps de maillage SIMP1");
 
optimisation(1,fem,nomparamaster,1);
 
cpu.ajouter_etape("temps calcul SIMP1");
 
  char message2[500];
  int nbiter=1;
  int nbitermax=params.get_valeur("iter_adaptative_max");
  int fin=0;
  double gdmax=1e-308;
 
  do 
        {
        char nomsortieiter[5000];
        sprintf(nomsortieiter,"%s_iter%d_gradient.sol",nometudesortie,nbiter);
        gradient_noeud(fem,gestd,nomsortieiter);
 
        LISTE_FEM_NOEUD::iterator it;
        double tenseur_i[9];
        double taille_min=1e308;
        for (FEM_NOEUD* no=gestd->get_fem_maillage(0)->get_premier_noeud(it);no!=NULL;no=gestd->get_fem_maillage(0)->get_suivant_noeud(it))
                {
                double *xyz=no->get_coord();
                carte->evaluer(xyz,tenseur_i);  // ancienne carte
                if (1./sqrt(tenseur_i[0]) <= taille_min) taille_min = 1./sqrt(tenseur_i[0]);
                }
 
        char message[500];
        sprintf(message,"Enregistrement resultat iteration %d",nbiter);
        affiche(message);
        sprintf(nomsortieiter,"%s_iter%d.magic",nometudesortie,nbiter);
        gestd->enregistrer(nomsortieiter);
 
        sprintf(message,"Calcul de la carte de taille");
        affiche(message);
        sprintf(nomsortieiter,"%s_iter%dctt.sol",nometudesortie,nbiter);
        FCT_TAILLE_FEM_SOLUTION_FONCTION_ISOTROPE nvcarte(gestd,gestd->get_fem_maillage(0),nomsortieiter);
        FCT_TAILLE_FEM_SOLUTION_DEFINITION_FONCTION_ISOTROPE::constante[0]=params.get_valeur("enm");
        FCT_TAILLE_FEM_SOLUTION_DEFINITION_FONCTION_ISOTROPE::constante[1]=-params.get_valeur("khi")*FCT_TAILLE_FEM_SOLUTION_DEFINITION_FONCTION_ISOTROPE::constante[0];
        nvcarte.construit_depuis_scalaire(gestd->get_fem_solution(1),0,&FCT_TAILLE_FEM_SOLUTION_DEFINITION_FONCTION_ISOTROPE::exponentielle);
        int nbdiv=params.get_valeur("carte_taille_structure_adapt");
        if (nbdiv>0)
            {
            sprintf(message,"Calcul de la carte de taille sur un support structure de %d,%d,%d divisions",nbdiv,nbdiv,nbdiv);
            affiche(message);
            nvcarte.structure_support(nbdiv,nbdiv,nbdiv);
            }
        
        
        sprintf(nomsortieiter,"%s_iter%d",nometudesortie,nbiter);
        nvcarte.enregistrer(nomsortieiter);
 
        delete gestd;
        delete gestnd;
        delete carte;
        carte==NULL;
        gestd=new MG_FILE(nomgestd);
        gestnd=new MG_FILE(nomgestnd);
 
        sprintf(nomsortieiter,"%s_iter%d.ctt",nometudesortie,nbiter);
        params.change_valeur("typecarte",3.);
        params.change_valeur("fichiercarte",nomsortieiter);
        params.change_valeur("festimation",20.);
        double seuilmin=params.get_valeur("seuil_convergence_adaptative");
        double val=std::max(params.get_valeur("coef_critere_convergence")*params.get_valeur("critere_convergence"),seuilmin);
        params.change_valeur("critere_convergence",val);
        nbiter++;
        sprintf(message,"Itération %d de la SIMP",nbiter);
        affiche(message);
 
        affiche((char*)"Preparation du modèle");
        fem=prepare_model(gestd,gestnd); 
 
        cpu.ajouter_etape("temps nouvelle carte et maillage");
        
        optimisation(nbiter,fem,nomparamaster,nbiter);
        
        cpu.ajouter_etape("temps optimisation");
        
 
        affiche((char*)"Enregistrement resultat");
        sprintf(nomsortieiter,"%s_iter%d.magic",nometudesortie,nbiter);
        gestd->enregistrer(nomsortieiter);
        
 
        double tenseur_j[9];
        double taille_min2 = 1.e308;
        double taille_max = 0.;
        for (FEM_NOEUD* no=gestd->get_fem_maillage(0)->get_premier_noeud(it);no!=NULL;no=gestd->get_fem_maillage(0)->get_suivant_noeud(it))
                {
                double *xyz=no->get_coord();
                carte->evaluer(xyz,tenseur_j);  // nouvelle carte
                if (1./sqrt(tenseur_j[0]) <= taille_min2) taille_min2 = 1./sqrt(tenseur_j[0]);
                if (1./sqrt(tenseur_j[0]) > taille_max) taille_max = 1./sqrt(tenseur_j[0]);
                }
        sprintf(message,"taille_min (ancienne carte) = %lf ; taille_min2 (nouvelle carte) = %lf",taille_min,taille_min2);
        affiche(message);
        //if ((taille_min2/taille_min) >= 0.8) fin=1;
        
        
        cpu.ajouter_etape("temps enregistrement iteration");
        std::vector<std::string> tab;
        cpu.get_tabfinal(tab);
        for (int i=0;i<tab.size();i++)
          printf("%s\n",(char*)tab[i].c_str());
        }
 
  while ((fin == 0) && (nbiter < nbitermax));
  if (nbiter==nbitermax) {sprintf(message2,"      Arrêt de la procedure d'optimisation adaptative apres %d cycles de raffinement",nbiter-1);affiche(message2);}
  else {sprintf(message2,"      Convergence de la procedure d'optimisation adaptative après %d cycles de raffinement",nbiter-1);affiche(message2);}
 
  int sauvesimple=params.get_valeur("enregistrementsimple");
  if (sauvesimple==1)
        {
        affiche((char*)"Enregistrement resultat simplifie");
        char nomfichier[5000];
        sprintf(nomfichier,"%s.opt",nometudesortie);
        enregistre_simple(nomfichier,fem);
        }
 
  affiche((char*)"Enregistrement finaux resultat sous GMSH");
  MG_EXPORT exp;
  char nomfichier[5000];
  sprintf(nomfichier,"%s_mg",nometudesortie);
  exp.gmsh(fem->get_mg_maillage(),nomfichier);
  sprintf(nomfichier,"%s_fem",nometudesortie);
  exp.gmsh(fem,nomfichier);
  
}