aster/src/mgopt_posttraitement.cpp

#include "gestionversion.h"
#include "mgopt_posttraitement.h"
#include "fem_solution.h"
#include "fem_maillage.h"
#include "mg_maillage.h"
#include "mg_gestionnaire.h"
#include "mg_triangle_peau.h"
#include "mg_geometrie_outils.h"
#include "fem_maillage_outils.h"
#include "tpl_fonctions_generiques.h"
#include "ot_algorithme_geometrique.h"
#include "mg_maillage_outils.h"
#include "mailleur_analyse.h"
#include "mg_face_element.h"
#include "mg_export.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>
#include <complex>
//ajout pour pouvoir utiliser les fonction venat de mailleur_analyse.cpp
#include "mailleur_analyse.h"
#include "m3d_triangle.h"
#include "fct_taille.h"
#include "mg_segment.h"
//---------------------------------------------------------------------------


MGOPT_POSTTRAITEMENT::MGOPT_POSTTRAITEMENT():affichageactif(0)
{
for (int i=0;i<9;i++) etat[i]=0;
params.ajouter("decoupage",0.,OT_PARAMETRES::DOUBLE,"0. découpage en conservant les mailles entieres 1. découpage par isodensite");
params.ajouter("seuil",0.5,OT_PARAMETRES::DOUBLE,"Valeur du seuil de densité pour les 2 découpages");
params.ajouter("nomdensextra","densiteextra.sol",OT_PARAMETRES::STRING,"Nom du fichier comportant la solution de la densité extrapolé au noeud");
params.ajouter("consimpose",1.,OT_PARAMETRES::DOUBLE,"0. Maillage non design non inclus dans le decoupage par mailles entieres 1. Maillage design inclus");
params.ajouter("consoptimise",1.,OT_PARAMETRES::DOUBLE,"0. Maillage optimise non inclus dans le decoupage par mailles entieres 1. Maillage optimise inclus");
params.ajouter("consnooptimise",0.,OT_PARAMETRES::DOUBLE,"0. Maillage non optimise non inclus dans le decoupage par mailles entieres 1. Maillage non optimise inclus");
params.ajouter("reactivation",1.,OT_PARAMETRES::DOUBLE,"1. Reactivation de mailles isolées sur le long des aretes 0. non reactivation");
params.ajouter("chen2005_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application du lissage de chen2005 0. méthode non utilisee");
params.ajouter("chen2005_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage de chen2005");
params.ajouter("chen2005_epsilon",0.01,OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2005");
params.ajouter("chen2005_sigma",0.1,OT_PARAMETRES::DOUBLE,"Écart type du filtre de lissage de chen2005");
params.ajouter("chen2005_filtre",1.,OT_PARAMETRES::DOUBLE,"Choix du filtre de lissage de chen2005 1. Gaussian 2.Laplacien 3.elfallahford");
params.ajouter("chen2008_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application deu lissage de chen2008 0. méthode non utilisee");
params.ajouter("chen2008_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage de chen2008");
params.ajouter("chen2008_epsilon",0.01,OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2008");
params.ajouter("chen2008_sigma",0.1,OT_PARAMETRES::DOUBLE,"Écart type du filtre de lissage de chen2008");
params.ajouter("chen2008_gamma",1.,OT_PARAMETRES::DOUBLE,"Vitesse d'avancement de déplacement du lissage de chen2008");
params.ajouter("chen2008_filtre",1.,OT_PARAMETRES::DOUBLE,"Choix du filtre de lissage de chen2008 1. Gaussian 2.Laplacien 3.elfallahford");
params.ajouter("jiao2012_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application deu lissage de jiao2012 0. méthode non utilisee");
params.ajouter("jiao2012_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage de jioa2012");
params.ajouter("McKenzie2016_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application du lissage de McKenzie2016 0. méthode non utilisee");
params.ajouter("McKenzie2016_lambda",0.33,OT_PARAMETRES::DOUBLE,"Paramète lambda, la valeur doit être positive");
params.ajouter("McKenzie2016_nu",-0.34,OT_PARAMETRES::DOUBLE,"Paramète nu, la valeur doit être negative et respecter 0<lambda<(-nu)");
params.ajouter("McKenzie2016_epsilon", 0.01000, OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2008");
params.ajouter("McKenzie2016_sigma", 0.10000, OT_PARAMETRES::DOUBLE,"Écart type du filtre de lissage de chen2008");
params.ajouter("McKenzie2016_gamma_", 1.0, OT_PARAMETRES::DOUBLE,"Vitesse d'avancement de déplacement du lissage de chen2008");
params.ajouter("McKenzie2016_filtre", 1.0, OT_PARAMETRES::DOUBLE,"Choix du filtre de lissage de chen2008 1. Gaussian 2.Laplacien 3.elfallahford");
params.ajouter("McKenzie2016_itertaubmax",100.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales de la partie taubin du lissage de McKenzie2016");
params.ajouter("McKenzie2016_itermax",100.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales de la partie chen 2008 du lissage de McKenzie2016");
params.ajouter("Taubin1995_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application du lissage de Taubin1995 0. méthode non utilisee");
params.ajouter("Taubin1995_lambda",0.33,OT_PARAMETRES::DOUBLE,"Paramète lambda, la valeur doit être positive et moins de 1");
params.ajouter("Taubin1995_nu",-0.34,OT_PARAMETRES::DOUBLE,"Paramète nu, la valeur doit être negative et respecter 0<lambda<(-nu)");
params.ajouter("Taubin1995_itermax",100.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage Taubin1995");
}


MGOPT_POSTTRAITEMENT::~MGOPT_POSTTRAITEMENT()
{
for (int i=0;i<lst_peau.size();i++)
        delete lst_peau[i];
}

void MGOPT_POSTTRAITEMENT::active_affichage(void (*fonc)(char*))
{
affiche=fonc;
affichageactif=1;
}


void MGOPT_POSTTRAITEMENT::posttraite(char *fichier)
{
affiche((char*)"");
affiche((char*)"***********************************************************");
affiche((char*)"Post-traitement d'un calcul d'optimisation de topologie");
affiche((char*)"***********************************************************");
affiche((char*)"");
affiche((char*)"");
affiche((char*)"");
affiche((char*)"Écriture d'un fichier de parametres par défaut");
params.enregistrer(fichier);
}


void MGOPT_POSTTRAITEMENT::conserve(int origine)
{
etat[(origine-1000)/10]=1;
}


void MGOPT_POSTTRAITEMENT::copieorigine(FEM_MAILLAGE* mai)
{
LISTE_FEM_ELEMENT3::iterator   it_tetra;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it_tetra);tet!=NULL;tet=mai->get_suivant_element3(it_tetra))
        {
        MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
        mgtet->change_origine(mgtet->get_origine());
        }
}


void MGOPT_POSTTRAITEMENT::gain_poids(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
{
LISTE_MG_TETRA::iterator   it_tetra;
double volume_initial = 0.;
double volume_final = 0.;
for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
        {
        int originetet=mgtet->get_origine();
        MG_NOEUD* noeud1 = mgtet->get_noeud1();
        double x1 = noeud1->get_x();
        double y1 = noeud1->get_y();
        double z1 = noeud1->get_z();
        MG_NOEUD* noeud2 = mgtet->get_noeud2();
        double x2 = noeud2->get_x();
        double y2 = noeud2->get_y();
        double z2 = noeud2->get_z();
        MG_NOEUD* noeud3 = mgtet->get_noeud3();
        double x3 = noeud3->get_x();
        double y3 = noeud3->get_y();
        double z3 = noeud3->get_z();
        MG_NOEUD* noeud4 = mgtet->get_noeud4();
        double x4 = noeud4->get_x();
        double y4 = noeud4->get_y();
        double z4 = noeud4->get_z();
        double volume_tet = (x2-x1)*((y3-y1)*(z4-z1)-(z3-z1)*(y4-y1))-(y2-y1)*((x3-x1)*(z4-z1)-(z3-z1)*(x4-x1)) + (z2-z1)*((x3-x1)*(y4-y1)-(y3-y1)*(x4-x1));
        volume_tet = volume_tet/6.;
        volume_initial = volume_initial + volume_tet;
        if (originetet != MAGIC::ORIGINE::MAILLEUR_AUTO)
                {
                volume_final = volume_final + volume_tet;
                }
        }
}

void MGOPT_POSTTRAITEMENT::reactivation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
{
cout << "   Reactivation d'elements manquants sur aretes vives et faces planes" << endl;
LISTE_MG_TETRA::iterator   it_tetra;
for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
        {
        int originetet=mgtet->get_origine();
        int compteur = 0;
        if (originetet == MAGIC::ORIGINE::MAILLEUR_AUTO)
                {
                MG_TRIANGLE* tri1 = mgtet->get_triangle1();
                MG_TRIANGLE* tri2 = mgtet->get_triangle2();
                MG_TRIANGLE* tri3 = mgtet->get_triangle3();
                MG_TRIANGLE* tri4 = mgtet->get_triangle4();
                int nb_tet1 = tri1->get_lien_tetra()->get_nb();
                for (int k = 0;k<nb_tet1;k++)
                        {
                        MG_TETRA* tet1 = tri1->get_lien_tetra()->get(k);
                        int ori1 = tet1->get_origine();
                        if ((tet1 != mgtet) && (ori1 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
                        }
                int nb_tet2 = tri2->get_lien_tetra()->get_nb();
                for (int k = 0;k<nb_tet2;k++)
                        {
                        MG_TETRA* tet2 = tri2->get_lien_tetra()->get(k);
                        int ori2 = tet2->get_origine();
                        if ((tet2 != mgtet) && (ori2 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
                        }
                int nb_tet3 = tri3->get_lien_tetra()->get_nb();
                for (int k = 0;k<nb_tet3;k++)
                        {
                        MG_TETRA* tet3 = tri3->get_lien_tetra()->get(k);
                        int ori3 = tet3->get_origine();
                        if ((tet3 != mgtet) && (ori3 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
                        }
                int nb_tet4 = tri4->get_lien_tetra()->get_nb();
                for (int k = 0;k<nb_tet4;k++)
                        {
                        MG_TETRA* tet4 = tri4->get_lien_tetra()->get(k);
                        int ori4 = tet4->get_origine();
                        if ((tet4 != mgtet) && (ori4 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
                        }
                if (compteur == 0) mgtet->change_origine(MAGIC::ORIGINE::OPTIMISE);
                }
        }
}


MG_MAILLAGE* MGOPT_POSTTRAITEMENT::extract_skin_par_decoupage(FEM_SOLUTION* sol,double limit,MG_GESTIONNAIRE& gest2,std::string nom)
{
affiche((char*)"Extraction de l'enveloppe par découpage");
sol->active_solution(0);
FEM_MAILLAGE* mai=sol->get_maillage();
affiche((char*)"    Extrapolation de la densité aux noeuds");
LISTE_FEM_NOEUD::iterator it;
for (FEM_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
    {
    double nume=0.;
    double deno=0.;
    int nb=no->get_lien_element3()->get_nb();
    int passe=0;
    for (int i=0;i<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];
          int col,ligne;
          double volume=tet->get_jacobien(jac,jac,ligne,col,1.);
          passe=1;
          nume=nume+tet->get_solution()*volume;
          deno=deno+volume;
          }
      }
     if (passe==1) no->change_solution(nume/deno); else no->change_solution(1.);               
    }
if (nom!="")
  {
  affiche((char*)"    Enregistrement de la densité aux noeuds");
  MG_GESTIONNAIRE *gest=mai->get_mg_maillage()->get_gestionnaire();
  std::string chemin=nom+".sol";
  FEM_SOLUTION* femdens=new FEM_SOLUTION(mai,1,(char*)chemin.c_str(),mai->get_nb_fem_element3(),"Extrapolation aux noeuds",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
  gest->ajouter_fem_solution(femdens);
  femdens->change_legende(0,"Densite");
  LISTE_FEM_ELEMENT3::iterator it3;
  int i=0;
  for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
      {
      MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
      if (tet2->get_origine()==MAGIC::ORIGINE::IMPOSE)
          {
          femdens->ecrire(1.,i,0,0);
          femdens->ecrire(1.,i,0,1);
          femdens->ecrire(1.,i,0,2);
          femdens->ecrire(1.,i,0,3);
          }
      else
          {
          femdens->ecrire(tet->get_fem_noeud(0)->get_solution(),i,0,0);
          femdens->ecrire(tet->get_fem_noeud(1)->get_solution(),i,0,1);
          femdens->ecrire(tet->get_fem_noeud(2)->get_solution(),i,0,2);
          femdens->ecrire(tet->get_fem_noeud(3)->get_solution(),i,0,3);
          }
      i++;
      }
  MG_EXPORT exp;
  exp.gmsh(mai,nom);
  gest->supprimer_fem_solution(gest->get_nb_fem_solution()-1);
  }
gest2.vide();
MG_GEOMETRIE* geo=new MG_GEOMETRIE((char*)"VIRTUEL",(char*)"VIRTUEL");
geo->change_valeur_unite(mai->get_mg_geometrie()->get_valeur_unite());
gest2.ajouter_mg_geometrie(geo);
MG_FACE_ELEMENT* face=new MG_FACE_ELEMENT();
geo->ajouter_mg_face(face);
MG_MAILLAGE* mgmai=new MG_MAILLAGE(geo);
gest2.ajouter_mg_maillage(mgmai);

MG_MAILLAGE* mgfem=mai->get_mg_maillage();
LISTE_MG_TRIANGLE::iterator itmg;
for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
  tri->change_nouveau_numero(0);
affiche((char*)"    Decoupage des tetra optimises");
LISTE_FEM_ELEMENT3::iterator it3;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
{
MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
if (tet2->get_origine()==MAGIC::ORIGINE::IMPOSE)
    {
    tet2->get_triangle1()->change_origine(MAGIC::ORIGINE::IMPOSE);
    tet2->get_triangle2()->change_origine(MAGIC::ORIGINE::IMPOSE);
    tet2->get_triangle3()->change_origine(MAGIC::ORIGINE::IMPOSE);
    tet2->get_triangle4()->change_origine(MAGIC::ORIGINE::IMPOSE);
    tet2->get_noeud1()->change_solution(tet->get_fem_noeud(0)->get_solution());
    tet2->get_noeud2()->change_solution(tet->get_fem_noeud(1)->get_solution());
    tet2->get_noeud3()->change_solution(tet->get_fem_noeud(2)->get_solution());
    tet2->get_noeud4()->change_solution(tet->get_fem_noeud(3)->get_solution());
    tet2->get_noeud1()->change_nouveau_numero(tet->get_fem_noeud(0)->get_id());
    tet2->get_noeud2()->change_nouveau_numero(tet->get_fem_noeud(1)->get_id());
    tet2->get_noeud3()->change_nouveau_numero(tet->get_fem_noeud(2)->get_id());
    tet2->get_noeud4()->change_nouveau_numero(tet->get_fem_noeud(3)->get_id());
    if (tet2->get_triangle1()->get_nouveau_numero()==0) 
      tet2->get_triangle1()->change_nouveau_numero(tet2->get_noeud4()->get_id());
    else 
      tet2->get_triangle1()->change_nouveau_numero(-1);
    if (tet2->get_triangle2()->get_nouveau_numero()==0) 
      tet2->get_triangle2()->change_nouveau_numero(tet2->get_noeud3()->get_id());
    else 
      tet2->get_triangle2()->change_nouveau_numero(-1);
    if (tet2->get_triangle3()->get_nouveau_numero()==0) 
      tet2->get_triangle3()->change_nouveau_numero(tet2->get_noeud1()->get_id());
    else 
      tet2->get_triangle3()->change_nouveau_numero(-1);
    if (tet2->get_triangle4()->get_nouveau_numero()==0) 
      tet2->get_triangle4()->change_nouveau_numero(tet2->get_noeud2()->get_id());
    else 
      tet2->get_triangle4()->change_nouveau_numero(-1);
    }
}

TPL_LISTE_ENTITE<MG_TRIANGLE*> tri_impose_interne;
for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
{
if (tri->get_origine()==MAGIC::ORIGINE::IMPOSE)
  if (tri->get_lien_topologie()->get_dimension()==3)
    if (tri->get_nouveau_numero()>0)
       tri_impose_interne.ajouter(tri);
}
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
    {
    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);
          std::vector<MG_NOEUD*> tab;
          interpole_segment(mai,no1,no2,&tab,limit,mgmai);
          interpole_segment(mai,no1,no3,&tab,limit,mgmai);
          interpole_segment(mai,no1,no4,&tab,limit,mgmai);
          interpole_segment(mai,no2,no3,&tab,limit,mgmai);
          interpole_segment(mai,no2,no4,&tab,limit,mgmai);
          interpole_segment(mai,no3,no4,&tab,limit,mgmai);
          int nb=tab.size();
          FEM_NOEUD* noext;
          if (nb>0)
            {
            if (no1->get_solution()<limit) noext=no1;
            if (no2->get_solution()<limit) noext=no2;
            if (no3->get_solution()<limit) noext=no3;
            if (no4->get_solution()<limit) noext=no4;
            }
          if (nb==3)
                  {
                  MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                  oriente_tri(tri,noext->get_coord());
                  }
          if (nb==4)
                  {
                  if (test_du_point_milieu(tab[0],tab[1],tet)==1)
                          {
                          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[1],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          oriente_tri(tri,noext->get_coord());
                          oriente_tri(tri2,noext->get_coord());
                          }
                  else if (test_du_point_milieu(tab[0],tab[2],tet)==1)
                          {
                          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[2],tab[1],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[2],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          oriente_tri(tri,noext->get_coord());
                          oriente_tri(tri2,noext->get_coord());}
                  else if (test_du_point_milieu(tab[0],tab[3],tet)==1)
                          {
                          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[3],tab[1],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[3],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          oriente_tri(tri,noext->get_coord());
                          oriente_tri(tri2,noext->get_coord());}
                  else if (test_du_point_milieu(tab[1],tab[2],tet)==1)
                          {
                          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[1],tab[2],tab[0],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          oriente_tri(tri,noext->get_coord());
                          oriente_tri(tri2,noext->get_coord());
                          }
                  else if (test_du_point_milieu(tab[1],tab[3],tet)==1)
                          {
                          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[1],tab[3],tab[0],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[3],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          oriente_tri(tri,noext->get_coord());
                          oriente_tri(tri2,noext->get_coord());}
                  else
                          {
                          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[2],tab[3],tab[0],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[2],tab[3],tab[1],mgmai,MAGIC::ORIGINE::TRIANGULATION);
                          oriente_tri(tri,noext->get_coord());
                          oriente_tri(tri2,noext->get_coord());
                          }
                  
                  }
          }
    }
affiche((char*)"    Decoupage des triangles non design interne au volume");
for (int i=0;i<tri_impose_interne.get_nb();i++)
    {
     MG_TRIANGLE* tri=tri_impose_interne.get(i);
    MG_NOEUD* no1=tri->get_noeud1();
    MG_NOEUD* no2=tri->get_noeud2();
    MG_NOEUD* no3=tri->get_noeud3();
     std::vector<MG_NOEUD*> tab;
     int num1=-1;
     int num2=-1;
     int num3=-1;
     if (no1->get_solution()<limit)
            {
            MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no1->get_nouveau_numero()),mgmai);
            tab.push_back(no);
            num1=1;
            }
    if (no2->get_solution()<limit)
            {
            MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no2->get_nouveau_numero()),mgmai);
            tab.push_back(no);
            num2=1;
            }
    if (no3->get_solution()<limit)
            {
            MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no3->get_nouveau_numero()),mgmai);
            tab.push_back(no);
            num3=1;
            }
    if (num1*num2<0) interpole_segment(mai,mai->get_fem_noeudid(no1->get_nouveau_numero()),mai->get_fem_noeudid(no2->get_nouveau_numero()),&tab,limit,mgmai);
    if (num1*num3<0) interpole_segment(mai,mai->get_fem_noeudid(no1->get_nouveau_numero()),mai->get_fem_noeudid(no3->get_nouveau_numero()),&tab,limit,mgmai);
    if (num2*num3<0) interpole_segment(mai,mai->get_fem_noeudid(no2->get_nouveau_numero()),mai->get_fem_noeudid(no3->get_nouveau_numero()),&tab,limit,mgmai);
    MG_NOEUD* noint=mgfem->get_mg_noeudid(tri->get_nouveau_numero());
    int nb=tab.size();
      if (nb==3)
          {
          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::IMPOSE);
          oriente_tri(tri,noint->get_coord());
          tri->inverse_sens();
          }
      if (nb==4)
          {
          MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::IMPOSE);
          MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,MAGIC::ORIGINE::IMPOSE);
          oriente_tri(tri,noint->get_coord());
          tri->inverse_sens();
          oriente_tri(tri2,noint->get_coord());
          tri2->inverse_sens();
          }
    }
affiche((char*)"    Decoupage des triangles sur la frontiere du volume");
LISTE_FEM_ELEMENT2::iterator it2;
for (FEM_ELEMENT2* tri=mai->get_premier_element2(it2);tri!=NULL;tri=mai->get_suivant_element2(it2))
    {
    std::vector<MG_NOEUD*> tab;
    FEM_NOEUD *no1=tri->get_fem_noeud(0);
    FEM_NOEUD *no2=tri->get_fem_noeud(1);
    FEM_NOEUD *no3=tri->get_fem_noeud(2);
    int ori=((MG_FACE*)tri->get_lien_topologie())->get_mg_coface(0)->get_orientation();
    OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
    OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
    OT_VECTEUR_3D nor=n1n3&n1n2;
    double xyzext[3];
    xyzext[0]=no1->get_x()+nor.get_x()*ori;
    xyzext[1]=no1->get_y()+nor.get_y()*ori;
    xyzext[2]=no1->get_z()+nor.get_z()*ori;
    if (tri->get_mg_element_maillage()->get_origine()==MAGIC::ORIGINE::IMPOSE)
      {
      MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai);
      MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai);
      MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai);
      int num1=interpole_segment(mai,no1,no2,&tab,limit,mgmai,0);
      int num2=interpole_segment(mai,no1,no3,&tab,limit,mgmai,0);
      int num3=interpole_segment(mai,no2,no3,&tab,limit,mgmai,0);
      int nb=tab.size();
      if (nb==0)
          {
          MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno2,mgno3,mgmai,MAGIC::ORIGINE::IMPOSE);
          oriente_tri(tri,xyzext);
          }
      if (nb==1)
        {
            if (num1==1)
                {
                MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
                MG_TRIANGLE_PEAU* tri2=insere_triangle(face,mgno2,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
                oriente_tri(tri,xyzext);
                oriente_tri(tri2,xyzext);
                }
            if (num2==1)
                {
                MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno2,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
                MG_TRIANGLE_PEAU* tri2=insere_triangle(face,mgno2,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
                oriente_tri(tri,xyzext);
                oriente_tri(tri2,xyzext);
                }
             if (num3==1)
                {
                MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno2,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
                MG_TRIANGLE_PEAU* tri2=insere_triangle(face,mgno1,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
                oriente_tri(tri,xyzext);
                oriente_tri(tri2,xyzext);
                }
        }
      if (nb==2)
        {
            if (num1==0)
              {
              MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno2,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
              MG_TRIANGLE_PEAU* tri2=insere_triangle(face,mgno2,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
              MG_TRIANGLE_PEAU* tri3=insere_triangle(face,mgno3,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
              oriente_tri(tri,xyzext);
              oriente_tri(tri2,xyzext);
              oriente_tri(tri3,xyzext);
              }
            if (num2==0)
              {
              MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
              MG_TRIANGLE_PEAU* tri2=insere_triangle(face,mgno3,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
              MG_TRIANGLE_PEAU* tri3=insere_triangle(face,mgno2,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
              oriente_tri(tri,xyzext);
              oriente_tri(tri2,xyzext);
              oriente_tri(tri3,xyzext);
              }
            if (num2==0)
              {
              MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno2,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
              MG_TRIANGLE_PEAU* tri2=insere_triangle(face,mgno3,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
              MG_TRIANGLE_PEAU* tri3=insere_triangle(face,mgno1,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
              oriente_tri(tri,xyzext);
              oriente_tri(tri2,xyzext);
              oriente_tri(tri3,xyzext);
              }
        }
      }
    else if ((no1->get_solution()>=limit) && (no2->get_solution()>=limit) && (no3->get_solution()>=limit)) 
          {
          MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai);
          MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai);
          MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai);
          MG_TRIANGLE_PEAU* tri=insere_triangle(face,mgno1,mgno2,mgno3,mgmai,MAGIC::ORIGINE::TRIANGULATION);
          oriente_tri(tri,xyzext);
          }
    else if (!((no1->get_solution()<limit) && (no2->get_solution()<limit) && (no3->get_solution()<limit)))
      {
      std::vector<MG_NOEUD*> tab;
      int num1=-1;
      int num2=-1;
      int num3=-1;
      if (no1->get_solution()>=limit)
            {
            MG_NOEUD* no=get_noeud_peau(no1,mgmai);
            tab.push_back(no);
            num1=1;
            }
      if (no2->get_solution()>=limit)
            {
            MG_NOEUD* no=get_noeud_peau(no2,mgmai);
            tab.push_back(no);
            num2=1;
            }
      if (no3->get_solution()>=limit)
            {
            MG_NOEUD* no=get_noeud_peau(no3,mgmai);
            tab.push_back(no);
            num3=1;
            }
      if (num1*num2<0) interpole_segment(mai,no1,no2,&tab,limit,mgmai);
      if (num1*num3<0) interpole_segment(mai,no1,no3,&tab,limit,mgmai);
      if (num2*num3<0) interpole_segment(mai,no2,no3,&tab,limit,mgmai);
      int nb=tab.size();
      if (nb==3)
          {
          MG_TRIANGLE_PEAU* tri=insere_triangle(face,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
          oriente_tri(tri,xyzext);
          }
      if (nb==4)
          {
          MG_TRIANGLE_PEAU* tri=insere_triangle(face,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
          MG_TRIANGLE_PEAU* tri2=insere_triangle(face,tab[1],tab[2],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
          oriente_tri(tri,xyzext);
          oriente_tri(tri2,xyzext);
          }
      }
    }
LISTE_MG_TRIANGLE::iterator it_tri;
for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
        {
        MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
        MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
        MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
        MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
        tripeau->change_voisin1(voisin1);
        tripeau->change_voisin2(voisin2);
        tripeau->change_voisin3(voisin3);
        tripeau->change_nouveau_numero(0);
        }
int fin;
do
        {
        fin=1;
        for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
                if (tripeau->get_nouveau_numero()==0)
                        {
                        fin=0;
                        std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
                        lst_peau.push_back(peau);
                        tripeau->change_nouveau_numero(1);
                        peau->push_back(tripeau);
                        determine_peau(peau);
                        }
                }
        }
while (fin==0);
return mgmai;
}


void MGOPT_POSTTRAITEMENT::suppression_peaux_isoles(MG_MAILLAGE* mgmai)
{
affiche((char*)"Suppression des peaux isolées"); 

char message[500];
for (int cas=0;cas<2;cas++)
{
if (cas==0) affiche((char*)"  Analyse initiale des peaux");
if (cas==1) affiche((char*)"  Analyse finale des peaux");
int nbisole=0;
for (int i=0;i<lst_peau.size();i++)
  {
  int isole=1;
  for (int j=0;j<lst_peau[i]->size();j++)
        if ((*lst_peau[i])[j]->get_origine()==MAGIC::ORIGINE::IMPOSE) {isole=0;break;}
  if (isole==1)
    {
    nbisole++;
    for (int j=0;j<lst_peau[i]->size();j++)
      {
      mgmai->supprimer_mg_triangleid((*lst_peau[i])[j]->get_id());
      }
    lst_peau[i]->clear();  
    }
  }
sprintf(message,"     %d peaux, %d non isoles, %d isoles",(int)lst_peau.size(),(int)lst_peau.size()-nbisole,nbisole);
affiche(message);
for (std::vector<std::vector<MG_TRIANGLE_PEAU*> *>::iterator j=lst_peau.end()-1;j>lst_peau.begin();j--)
  if ((*j)->size()==0)
      lst_peau.erase(j);
}
LISTE_MG_TRIANGLE::iterator itmg;
for (MG_TRIANGLE* tri=mgmai->get_premier_triangle(itmg);tri!=NULL;tri=mgmai->get_suivant_triangle(itmg))
  {
  if (tri->get_origine()==MAGIC::ORIGINE::IMPOSE)
        {
        tri->get_noeud1()->change_origine(MAGIC::ORIGINE::IMPOSE);
        tri->get_noeud2()->change_origine(MAGIC::ORIGINE::IMPOSE);
        tri->get_noeud3()->change_origine(MAGIC::ORIGINE::IMPOSE);
        }
  }

}


void MGOPT_POSTTRAITEMENT::oriente_tri(MG_TRIANGLE_PEAU* tri,double *xyz)
{
MG_NOEUD* no1=tri->get_noeud1();
MG_NOEUD* no2=tri->get_noeud2();
MG_NOEUD* no3=tri->get_noeud3();
OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
OT_VECTEUR_3D normal=n1n3&n1n2;
OT_VECTEUR_3D dir(no1->get_coord(),xyz);
dir.norme();
normal.norme();
double ps=normal*dir;
if (ps<0) tri->inverse_sens();
}


MG_NOEUD* MGOPT_POSTTRAITEMENT::get_noeud_peau(FEM_NOEUD* no,MG_MAILLAGE* mai)
{
static std::map<std::string,MG_NOEUD*> correspond;
unsigned long id=no->get_id();
char message[255];
sprintf(message,"_%lu_",id);
std::string key=message;
MG_NOEUD* mgno=correspond[key];
if (mgno==NULL)
    {
    int origine=MAGIC::ORIGINE::TRIANGULATION;
    if (((MG_NOEUD*)no->get_mg_element_maillage())->get_lien_topologie()->get_dimension()==1) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
    mgno=new MG_NOEUD(NULL,no->get_x(),no->get_y(),no->get_z(),origine);
    mai->ajouter_mg_noeud(mgno);
    correspond[key]=mgno;
    }
return mgno;
}


int MGOPT_POSTTRAITEMENT::test_du_point_milieu(MG_NOEUD* no1,MG_NOEUD* no2,FEM_ELEMENT3* tet)
{
double *xyz1=no1->get_coord();
double *xyz2=no2->get_coord();
double xyz[3];
xyz[0]=0.5*(xyz1[0]+xyz2[0]);
xyz[1]=0.5*(xyz1[1]+xyz2[1]);
xyz[2]=0.5*(xyz1[2]+xyz2[2]);
FEM_MAILLAGE_OUTILS ot;
if (((ot.estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
//if (((MG_MAILLAGE_OUTILS::estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
return 0;
}


int MGOPT_POSTTRAITEMENT::interpole_segment(FEM_MAILLAGE* fem,FEM_NOEUD* no1,FEM_NOEUD* no2,std::vector<MG_NOEUD*> *tab,double limit,MG_MAILLAGE* mai,int creation)
{
static std::map<std::string,MG_NOEUD*> correspond;
unsigned long id1=no1->get_id();
unsigned long id2=no2->get_id();
char message[255];
sprintf(message,"_%lu_%lu_",id1,id2);
std::string key1=message;
sprintf(message,"_%lu_%lu_",id2,id1);
std::string key2=message;
double sol1=no1->get_solution();
double sol2=no2->get_solution();
MG_NOEUD* n1=(MG_NOEUD*)no1->get_mg_element_maillage();
MG_NOEUD* n2=(MG_NOEUD*)no2->get_mg_element_maillage();
MG_SEGMENT *seg=fem->get_mg_maillage()->get_mg_segment(n1->get_id(),n2->get_id());
int origine=MAGIC::ORIGINE::TRIANGULATION;
if (seg->get_lien_topologie()->get_dimension()==2) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
if (seg->get_lien_topologie()->get_dimension()==1) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;

if (fabs(sol1-limit)<0.0000001)
      {
        sprintf(message,"%lu",id1);
        std::string key=message;
        MG_NOEUD* no=correspond[key];
        if (no==NULL)
          {
          double *xyz1=no1->get_coord();
          no=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],origine);
          mai->ajouter_mg_noeud(no);
          correspond[key]=no;
          }
        int present=0;
        for (int i=0;i<tab->size();i++)
              if ((*tab)[i]==no) present=1;
        if (present==0) 
          {
          tab->push_back(no);
          return 1;
          }
        else return 0;
      }
if (fabs(sol2-limit)<0.0000001)
   {
        sprintf(message,"%lu",id2);
        std::string key=message;
        MG_NOEUD* no=correspond[key];
        if (no==NULL)
          {
          double *xyz2=no2->get_coord();
          no=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],origine);
          mai->ajouter_mg_noeud(no);
          correspond[key]=no;
          }
        int present=0;
        for (int i=0;i<tab->size();i++)
              if ((*tab)[i]==no) present=1;
        if (present==0) 
          {
          tab->push_back(no);
          return 1;
          }
        else return 0;
          
      }
      
if (((sol1<limit) && (sol2>limit))||((sol1>limit) && (sol2<limit)))
    {
    MG_NOEUD* no=correspond[key1];
    if (no==NULL) no=correspond[key2];
    if ((no==NULL) && (creation==1))
      {
      double t=(limit-sol1)/(sol2-sol1);
      double xyz[3];
      double *xyz1=no1->get_coord();
      double *xyz2=no2->get_coord();
      xyz[0]=xyz1[0]+t*(xyz2[0]-xyz1[0]);
      xyz[1]=xyz1[1]+t*(xyz2[1]-xyz1[1]);
      xyz[2]=xyz1[2]+t*(xyz2[2]-xyz1[2]);
      no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],origine);
      mai->ajouter_mg_noeud(no);
      correspond[key1]=no;
      }
    if (no!=NULL) 
      {
      tab->push_back(no);
      return 1;
      }
    }
return 0;
}

void MGOPT_POSTTRAITEMENT::lire_params(char *fichier)
{
params.lire(fichier);
}

void MGOPT_POSTTRAITEMENT::posttraite(FEM_SOLUTION* sol, MG_GESTIONNAIRE& gest2,char *nomparam)
{
if (nomparam!=NULL) lire_params(nomparam);
affiche((char*)"Extraction du maillage de surface");
FEM_MAILLAGE* fem=sol->get_maillage();
int decoup=params.get_valeur("decoupage");
MG_MAILLAGE* mgmai;
if (decoup==0)
  {
  affiche((char*)"       Découpage par mailles 3D entieres");
  adapte_seuil(fem,sol); 
  mgmai=extract_skin_maille_entiere(fem,gest2);
  }
if (decoup==1)
  {
  affiche((char*)"       Découpage par isodensité");
  double seuil=params.get_valeur("seuil");
  std::string nomfich=params.get_nom("nomdensextra");
  mgmai=extract_skin_par_decoupage(sol,seuil,gest2,nomfich.c_str());
  }
suppression_peaux_isoles(mgmai);
affiche((char*)"Procedure de lissage");
int chen2005_debut=params.get_valeur("chen2005_debut");
int chen2008_debut=params.get_valeur("chen2008_debut");
int jiao2012_debut=params.get_valeur("jiao2012_debut");
int McKenzie2016_debut=params.get_valeur("McKenzie2016_debut");
int Taubin1995_debut=params.get_valeur("Taubin1995_debut");
int chen2005fait=0;
if (chen2005_debut==0) chen2005fait=1;
int chen2008fait=0;
if (chen2008_debut==0) chen2008fait=1;
int jiao2012fait=0;
if (jiao2012_debut==0) jiao2012fait=1;
int McKenzie2016fait=0;
if (McKenzie2016_debut==0) McKenzie2016fait=1;
int Taubin1995fait=0;
if (Taubin1995_debut==0) Taubin1995fait=1;

int fin=0;
int iteration=1;
do      
  {     
  if ((iteration>=chen2005_debut) && (chen2005fait==0))
      {
      char message[300];
      sprintf(message, "   Iteration %d : lissage de Chen 2005",iteration);
      affiche(message);
      double epsilon=params.get_valeur("chen2005_epsilon");
      double sigma=params.get_valeur("chen2005_sigma");
      int iter_max=params.get_valeur("chen2005_itermax");
      int nbiter=lissage_chen2005(mgmai,gest2,epsilon,sigma,iter_max);
      if (nbiter==iter_max) sprintf(message,"      Arrêt de la procédure de lissage de Chen2005 après %d itérations",nbiter);
                else sprintf(message,"      Convergence de la procédure de lissage de Chen2005 après %d itérations",nbiter);
      affiche(message);
      chen2005fait=1;
      iteration=iteration+nbiter;
      }
  else if ((iteration>=chen2008_debut) && (chen2008fait==0))
      {
      char message[300];
      sprintf(message, "   Iteration %d : lissage de Chen 2008",iteration);
      affiche(message);
      double epsilon=params.get_valeur("chen2008_epsilon");
      double sigma=params.get_valeur("chen2008_sigma");
      double gamma=params.get_valeur("chen2008_gamma");
      int iter_max=params.get_valeur("chen2008_itermax");
      int nbiter=lissage_chen2008(mgmai,gest2,sigma,gamma,epsilon,iter_max);
      if (nbiter==iter_max) sprintf(message,"      Arrêt de la procédure de lissage de Chen2008 après %d itérations",nbiter);
                else sprintf(message,"      Convergence de la procédure de lissage de Chen2008 après %d itérations",nbiter);
      affiche(message);
      chen2008fait=1;
      iteration=iteration+nbiter;
      }
  else if ((iteration>=jiao2012_debut) && (jiao2012fait==0))
      {
      char message[300];
      sprintf(message, "   Iteration %d : lissage de Jiao 2012",iteration);
      affiche(message);
      int iter_max=params.get_valeur("jiao2012_itermax");
      int nbiter=lissage_jiao2012(mgmai,gest2,iter_max);
      if (nbiter==iter_max) sprintf(message,"      Arrêt de la procédure de lissage de Jiao2012 après %d itérations",nbiter);
                else sprintf(message,"      Convergence de la procédure de lissage de Jiao2012 après %d itérations",nbiter);
      affiche(message);
      jiao2012fait=1;
      iteration=iteration+nbiter;
      }
  else if ((iteration>=McKenzie2016_debut) && (McKenzie2016fait==0))
      {
      char message[300];
      sprintf(message, "Lissage de McKenzie2016");
      affiche(message);
      double lambda=params.get_valeur("McKenzie2016_lambda");
      double nu=params.get_valeur("McKenzie2016_nu");
      double epsilon=params.get_valeur("McKenzie2016_epsilon");
      double sigma=params.get_valeur("McKenzie2016_sigma");
      double gamma_ =params.get_valeur("McKenzie2016_sigma");
      int filtre=params.get_valeur("McKenzie2016_filtre");
      int iter_max=params.get_valeur("McKenzie2016_itermax");
      int itr_taubin=params.get_valeur("McKenzie2016_itertaubmax");
      int nbiter=lissage_McKenzie2016(mgmai,gest2,lambda,nu,epsilon,sigma,gamma_,filtre, iter_max, itr_taubin);
       if (nbiter==iter_max) sprintf(message,"      Arrêt de la procédure de lissage de McKenzie2016 après %d itérations",nbiter);
                else sprintf(message,"      Convergence de la procédure de lissage de McKenzie2016 après %d itérations",nbiter);
      affiche(message);
      McKenzie2016fait=1;
      iteration=iteration+nbiter;
      }
      
       else if ((iteration>=Taubin1995_debut) && (Taubin1995fait==0))
      {
      char message[300];
      sprintf(message, "Lissage de Taubin1995");
      affiche(message);
      double lambda=params.get_valeur("Taubin1995_lambda");
      double nu=params.get_valeur("Taubin1995_nu");
      int iter_max=params.get_valeur("Taubin1995_itermax");
      int nbiter=lissage_Taubin1995(mgmai,gest2,lambda,nu, iter_max);
       if (nbiter==iter_max) sprintf(message,"      Arrêt de la procédure de lissage de Taubin1995 après %d itérations",nbiter);
                else sprintf(message,"      Convergence de la procédure de lissage de Taubin1995 après %d itérations",nbiter);
      affiche(message);
      Taubin1995fait=1;
      iteration=iteration+nbiter;
      }
      
   else iteration++;
   if (chen2005fait==1)
    if (chen2008fait==1)
      if (jiao2012fait==1)
          if (McKenzie2016fait==1)
            if (Taubin1995fait==1)
          fin=1;
  }
while (fin==0);


}  
  

MG_MAILLAGE* MGOPT_POSTTRAITEMENT::extract_skin_maille_entiere(FEM_MAILLAGE* mai,MG_GESTIONNAIRE &gest2)
{  
//Menu de la methode de lissage
affiche((char*)"Extraction du maillage de surface");
int coderesu = 0;
int mai2_id;
int imp=params.get_valeur("consimpose");
int opti=params.get_valeur("consoptimise");
int m_auto=params.get_valeur("consnooptimise");
if (opti==1) conserve(MAGIC::ORIGINE::OPTIMISE);
if (imp==1) conserve(MAGIC::ORIGINE::IMPOSE);
if (m_auto==1) conserve(MAGIC::ORIGINE::MAILLEUR_AUTO); 
copieorigine(mai);
MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
//gain_poids(mg_mai,gest2);
int reactiv=params.get_valeur("reactivation");
if (reactiv == 1)
        {
        reactivation(mg_mai,gest2);
        }
affiche((char*)"   Analyse des cas non manifold");
do
        {
        int nbmaniare,nbmanino,nbpeau;
        coderesu = extract_skin(mg_mai,gest2,nbpeau,nbmaniare,nbmanino,&mai2_id);
        char message[500];
        sprintf(message,"             %d peaux. %d manifold par arete. %d manifold par noeud",nbpeau,nbmaniare,nbmanino);
        affiche(message);
//      gain_poids(mg_mai,gest2);
        }
while (coderesu == 0);

MG_MAILLAGE* mg_mai2=gest2.get_mg_maillageid(mai2_id);

return mg_mai2;
  
}


void MGOPT_POSTTRAITEMENT::adapte_seuil(class FEM_MAILLAGE* fem,FEM_SOLUTION* solution)
{
double seuil=params.get_valeur("seuil");
solution->active_solution(0);
LISTE_FEM_ELEMENT3::iterator it;
for (FEM_ELEMENT3 *tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
        {
        if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=MAGIC::ORIGINE::IMPOSE) 
                if (tet->get_solution()>seuil) 
                        ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::OPTIMISE); 
                else
                        ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::MAILLEUR_AUTO); 
        
        }
}


//---------------------------- Lissage Chen 2005 Gilles Philippe -------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::lissage_chen2005(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
{
//Methode de lissage inspiree de Chen(2005) par Gilles Philipe Picher Martel
double un_sur_pi = 1./M_PI;
int compteur = 0;
int fin = 0;


do
        {
        TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
        TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
        TPL_LISTE_ENTITE<double> liste_wij;
        LISTE_MG_TRIANGLE::iterator   it_tri;
        int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
        for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
                OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
                normal_f_i.norme();
                liste_normales2.ajouter(normal_f_i);
                //Remplissage de la liste des voisins du triangle i
                TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
                MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
                double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins1;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_1);
                        }
                MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
                double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins2;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_2);
                        }
                MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
                double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins3;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_3);
                        }
                liste_voisins.supprimer(mgtri_i);
                int nb_voisins = liste_voisins.get_nb();
                double w_ij = 1./nb_voisins;
                double phi_i_min = 10.;
                double s_i = 0.0;
                double phi_im = 0.0;
                double *phi_ij = new double[nb_voisins];
                OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
                OT_VECTEUR_3D eta_i(0.,0.,0.);
                TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
                int j = 0;
                for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
                        {
                        OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
                        //1-Calculer la normale moyenne pour chaque triangle 
                        normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j; 
                        //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
                        double prod_scalaire = normal_f_i*normal_f_j;
                        if (prod_scalaire > 1.)
                                { 
                                prod_scalaire = 1.;
                                }
                        if (prod_scalaire < -1.)
                                {
                                prod_scalaire = -1.;
                                }
                        phi_ij[j] = acos(prod_scalaire)*un_sur_pi; 
                        //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
                        if (phi_ij[j] < phi_i_min)
                                {
                                phi_i_min = phi_ij[j]; 
                                eta_i = normal_f_j;
                                }
                        //3.1-On calcul l'angle moyen phi_im
                        phi_im = phi_im + w_ij*phi_ij[j];       
                        j++;
                        }
                normal_f_i_mean.norme();
                j = 0;
                for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
                        {
                        //3.2-Calcul de s_i selon la variance
                        s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
                        j++;
                        }
                delete[] phi_ij;
                //4-On calcule une nouvelle normale pour chaque triangle
                double pond;
                int num=params.get_valeur("chen2005_filtre");
                if (num==1) pond=ponderation_gaussian(s_i,sigma);
                else if (num==2) pond=ponderation_laplacian(s_i,sigma);
                else if (num==3) pond=ponderation_elfallahford(s_i,sigma);
                OT_VECTEUR_3D normal_f_i_new = pond*normal_f_i_mean + (1. - pond)*eta_i;
                normal_f_i_new.norme();
                liste_normales.ajouter(normal_f_i_new);
                liste_wij.ajouter(w_ij);
                mgtri->change_nouveau_numero(k);
                k++;
                }
        
        LISTE_MG_NOEUD::iterator   it_no;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
                double w_ij_prime = 0.0;
                OT_VECTEUR_3D v_temp(0.,0.,0.);
                OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
                for(int j=0;j<nb_voisins_j;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
                        //On calcule le centroide cj du triangle mgtri_j
                        MG_NOEUD* n1 = mgtri_j->get_noeud1();
                        MG_NOEUD* n2 = mgtri_j->get_noeud2();
                        MG_NOEUD* n3 = mgtri_j->get_noeud3();
                        double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
                        double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
                        double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
                        //On forme le vecteur vi_cj
                        OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
                        OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
                        // w_ij_prime correspond a la somme des aires des triangles voisins du noeuds
                        OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
                        OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
                        OT_VECTEUR_3D prodvect = AB&AC;
                        double w_ij = 0.5*prodvect.get_longueur();
                        w_ij_prime = w_ij_prime + w_ij;
                        v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new; 
                        }
                //5-On met a jour la position des noeuds
                v_i = v_i + v_temp/w_ij_prime;
                int origine = noeud_i->get_origine();
                if (origine != MAGIC::ORIGINE::IMPOSE)
                        {
                        noeud_i->change_x(v_i.get_x());
                        noeud_i->change_y(v_i.get_y());
                        noeud_i->change_z(v_i.get_z());
                        }
                }
        //Critere d'arret de l'algorithme
        int l=0;
        int nb_tri = mg_mai->get_nb_mg_triangle();
        for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
                OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
                OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
                double critere = 1. - normal_f_i*normal_f_i_new;
                if (critere <= epsilon) l++;
                }
        double tolerance = 0.01*nb_tri;
        if (nb_tri - l <= 0) fin = 1;
        compteur++;
        }
while ((fin == 0) && (compteur < iter_max));

return compteur;
} 

//---------------------------- Lissage Jiao 2012 -------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::lissage_jiao2012(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, int iter_max)
{
// Lissage global avec methode de lissage inspiree de JIAO IMR2012 par Jean Christophe sept 2013
int compteur = 0;

do
        {
        vector<double> nouv_position_x;
        vector<double> nouv_position_y;
        vector<double> nouv_position_z;
        LISTE_MG_NOEUD::iterator   it_no;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                  if (origine != MAGIC::ORIGINE::IMPOSE)
                  {  
                  int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
                  double wij = 0.0;
                  double wij_sommej = 0.0;
                  OT_VECTEUR_3D v_temp(0.,0.,0.);
                  OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
                  for(int j=0;j<nb_voisins_j;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
                        //On calcule le centroide cj du triangle mgtri_j
                        MG_NOEUD* n1 = mgtri_j->get_noeud1();
                        MG_NOEUD* n2 = mgtri_j->get_noeud2();
                        MG_NOEUD* n3 = mgtri_j->get_noeud3();
                        double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
                        double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
                        double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
                        //On forme le vecteur vi_cj et le vecteur cj
                        OT_VECTEUR_3D cj(cj_x ,cj_y,cj_z);
                        OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
                        wij=vi_cj.get_longueur();
                        wij_sommej=wij_sommej+wij;
                        v_temp = v_temp + (wij*cj); 
                        }
                  //5-On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
                  v_i =v_temp/wij_sommej;
                  nouv_position_x.push_back(v_i.get_x());
                  nouv_position_y.push_back(v_i.get_y());
                  nouv_position_z.push_back(v_i.get_z());
                  }
                }
        //On actualise la position des noeuds
        int ind_noeud=0;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                  if (origine != MAGIC::ORIGINE::IMPOSE)
                  {  
                  noeud_i->change_x(nouv_position_x[ind_noeud]);
                  noeud_i->change_y(nouv_position_y[ind_noeud]);
                  noeud_i->change_z(nouv_position_z[ind_noeud]);  
                  ind_noeud++;
                  }
                }
                //Critere d'arret de l'algorithme
        compteur++;
        }
while (compteur < iter_max);

return compteur;
}

//---------------------------- Lissage Chen 2008 -------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::lissage_chen2008(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigma, double gamma, double epsilon, int iter_max)
{
//Lissage global inspire de Chen(2008) corrige par Jean Christophe sept 2013
double un_sur_pi = 1./M_PI;
int compteur = 0;
int fin = 0;
compteur = 0;

do
        {
        vector<double> nouv_position_x;
        vector<double> nouv_position_y;
        vector<double> nouv_position_z;  
        TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
        TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
        TPL_LISTE_ENTITE<double> liste_wij;
        LISTE_MG_TRIANGLE::iterator   it_tri;
        int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
        for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
                OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
                normal_f_i.norme();
                liste_normales2.ajouter(normal_f_i);
                //Remplissage de la liste des voisins du triangle i
                TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
                MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
                double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins1;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_1);
                        }
                MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
                double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins2;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_2);
                        }
                MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
                double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins3;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_3);
                        }
                liste_voisins.supprimer(mgtri_i);
                int nb_voisins = liste_voisins.get_nb();
                double w_ij = 1./nb_voisins;
                double phi_i_min = 10.;
                double s_i = 0.0;
                double phi_im = 0.0;
                double *phi_ij = new double[nb_voisins];
                OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
                normal_f_i_mean = normal_f_i;
                OT_VECTEUR_3D eta_i(0.,0.,0.);
                TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
                int j = 0;
                for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
                        {
                        OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
                        //1-Calculer la normale moyenne pour chaque triangle 
                        normal_f_i_mean = normal_f_i_mean + normal_f_j; 
                        //2.1-On calcule l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
                        double prod_scalaire = normal_f_i*normal_f_j;
                        if (prod_scalaire > 1.)
                                { 
                                prod_scalaire = 1.;
                                }
                        if (prod_scalaire < -1.)
                                {
                                prod_scalaire = -1.;
                                }
                        phi_ij[j] = acos(prod_scalaire)*un_sur_pi; 
                        //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
                        if (phi_ij[j] < phi_i_min)
                                {
                                phi_i_min = phi_ij[j]; 
                                eta_i = normal_f_j;
                                }
                        //3.1-On calcule l'angle moyen phi_im
                        phi_im = phi_im + w_ij*phi_ij[j];       
                        j++;
                        }
                normal_f_i_mean.norme();
                j = 0;
                for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
                        {
                        //3.2-Calcul de s_i selon la variance
                        s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
                        j++;
                        }
                delete[] phi_ij;
                //4-On calcule une nouvelle normale pour chaque triangle
                double pond;
                int num=params.get_valeur("chen2008_filtre");
                if (num==1) pond=ponderation_gaussian(s_i,sigma);
                else if (num==2) pond=ponderation_laplacian(s_i,sigma);
                else if (num==3) pond=ponderation_elfallahford(s_i,sigma);
                OT_VECTEUR_3D normal_f_i_new = pond*normal_f_i_mean + (1. - pond)*eta_i;
                normal_f_i_new.norme();
                liste_normales.ajouter(normal_f_i_new);
                liste_wij.ajouter(w_ij);
                mgtri->change_nouveau_numero(k);
                k++;
                }
        
        LISTE_MG_NOEUD::iterator   it_no;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                  if (origine != MAGIC::ORIGINE::IMPOSE)
                  {
                  int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
                  double w_ij_prime = 0.0;
                  OT_VECTEUR_3D v_temp(0.,0.,0.);
                  OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
                        for(int j=0;j<nb_voisins_j;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
                        //On calcule le centroide cj du triangle mgtri_j
                        MG_NOEUD* n1 = mgtri_j->get_noeud1();
                        MG_NOEUD* n2 = mgtri_j->get_noeud2();
                        MG_NOEUD* n3 = mgtri_j->get_noeud3();
                        double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
                        double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
                        double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
                        //On forme le vecteur vi_cj
                        OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
                        OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
                        v_temp = v_temp + (vi_cj*normal_f_i_new)*normal_f_i_new; 
                        }
                //5-On met a jour la position des noeuds
                  v_i = v_i + (gamma/(2*nb_voisins_j))*v_temp;
                  nouv_position_x.push_back(v_i.get_x());
                  nouv_position_y.push_back(v_i.get_y());
                  nouv_position_z.push_back(v_i.get_z());
                  }
                }

        //On actualise la position des noeuds
        int ind_noeud=0;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                  if (origine != MAGIC::ORIGINE::IMPOSE)
                  {  
                  noeud_i->change_x(nouv_position_x[ind_noeud]);
                  noeud_i->change_y(nouv_position_y[ind_noeud]);
                  noeud_i->change_z(nouv_position_z[ind_noeud]);  
                  ind_noeud++;
                  }
                }
        
        //Critere d'arret de l'algorithme
        int l=0;
        int nb_tri = mg_mai->get_nb_mg_triangle();
        for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
                OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
                OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
                double critere = 1. - normal_f_i*normal_f_i_new;
                if (critere <= epsilon) l++;
                }
        double tolerance = 0.01*nb_tri;
        if (nb_tri - l <= 0) fin = 1;
        compteur++;
        }
while ((fin == 0) && (compteur < iter_max));

return compteur;
  
}


//---------------------------- Cycle Taubin 1995 par Bryan McKnezie 2016-------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::cycle_taubin1995(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double lambda, double nu)
{
  
      /////partie lambda du lissage
        vector<double> nouv_position_x;
        vector<double> nouv_position_y;
        vector<double> nouv_position_z;
        LISTE_MG_NOEUD::iterator   it_no;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                if ((origine != MAGIC::ORIGINE::IMPOSE)&& (noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&(noeud_i->get_lien_topologie()==NULL))///////////////////////////////////////////////////////////////////////////////
                  {  
                  TPL_MAP_ENTITE<MG_NOEUD*> liste_voisin;
                  MG_MAILLAGE_OUTILS ot;
                  ot.get_noeud_voisin(noeud_i,liste_voisin,2);
                    
                  int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
                  double delta_x = 0.0;
                  double delta_y = 0.0;
                  double delta_z = 0.0;
                  double v_x = noeud_i->get_x();//v_1->get_x();
                  double v_y = noeud_i->get_y();//v_1->get_y();
                  double v_z = noeud_i->get_z();//v_1->get_z();
                  double multiplicateur = (1.0/nb_voisins_j);
                        
                  for(int j=0;j<nb_voisins_j;j++)
                        {
                        double v_i_x=liste_voisin.get(j)->get_x();
                        double v_i_y=liste_voisin.get(j)->get_y();
                        double v_i_z=liste_voisin.get(j)->get_z();
                        delta_x = delta_x+(multiplicateur*(v_i_x-v_x));
                        delta_y = delta_y+(multiplicateur*(v_i_y-v_y));
                        delta_z = delta_z+(multiplicateur*(v_i_z-v_z));

                        }
                  double x_prime = (v_x)+(lambda*delta_x);
                  double y_prime = (v_y)+(lambda*delta_y);
                  double z_prime = (v_z)+(lambda*delta_z);
                        
                  //On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
                  nouv_position_x.push_back(x_prime);
                  nouv_position_y.push_back(y_prime);
                  nouv_position_z.push_back(z_prime);
                  }
                }
        //On actualise la position des noeuds
        int ind_noeud=0;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                  if ((origine !=  MAGIC::ORIGINE::IMPOSE)&& (noeud_i->get_lien_topologie()==NULL))//((origine !=  MAGIC::ORIGINE::IMPOSE)&& (origine !=  MAGIC::ORIGINE::TRIANGULATION_ARETE)&&(noeud_i->get_lien_topologie()==NULL))///////////////////////////////////////////////////////////////////////////////
                  {  
                  noeud_i->change_x(nouv_position_x[ind_noeud]);
                  noeud_i->change_y(nouv_position_y[ind_noeud]);
                  noeud_i->change_z(nouv_position_z[ind_noeud]); 
                  ind_noeud++;
                  }
                }
                
              
                //partie nu du lissage
              nouv_position_x.clear();
              nouv_position_y.clear();
              nouv_position_z.clear();
              for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                if ((origine != MAGIC::ORIGINE::IMPOSE)&& (noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&(noeud_i2->get_lien_topologie()==NULL))///////////////////////////////////////////////////////////////////////////////
                  {  
                  TPL_MAP_ENTITE<MG_NOEUD*> liste_voisin;
                  MG_MAILLAGE_OUTILS ot;
                  ot.get_noeud_voisin(noeud_i,liste_voisin,2);
                    
                  int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
                  double delta_x = 0.0;
                  double delta_y = 0.0;
                  double delta_z = 0.0;
                  double v_x = noeud_i->get_x();//v_1->get_x();
                  double v_y = noeud_i->get_y();//v_1->get_y();
                  double v_z = noeud_i->get_z();//v_1->get_z();
                  double multiplicateur = (1.0/nb_voisins_j);
                        
                  for(int j=0;j<nb_voisins_j;j++)
                        {
                        double v_i_x=liste_voisin.get(j)->get_x();
                        double v_i_y=liste_voisin.get(j)->get_y();
                        double v_i_z=liste_voisin.get(j)->get_z();
                        delta_x = delta_x+(multiplicateur*(v_i_x-v_x));
                        delta_y = delta_y+(multiplicateur*(v_i_y-v_y));
                        delta_z = delta_z+(multiplicateur*(v_i_z-v_z));

                        }
                  double x_prime = (v_x)+(nu*delta_x);
                  double y_prime = (v_y)+(nu*delta_y);
                  double z_prime = (v_z)+(nu*delta_z);
                        
                  //On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
                  nouv_position_x.push_back(x_prime);
                  nouv_position_y.push_back(y_prime);
                  nouv_position_z.push_back(z_prime);
                }
        }
        //On actualise la position des noeuds
        ind_noeud=0;
        for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
                {
                int origine = noeud_i->get_origine();
                  if ((origine != MAGIC::ORIGINE::IMPOSE)&& ( noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&( noeud_i2->get_lien_topologie()==NULL))/////////////////////////////////////////////////////////////////////////////
                  {  
                  noeud_i->change_x(nouv_position_x[ind_noeud]);
                  noeud_i->change_y(nouv_position_y[ind_noeud]);
                  noeud_i->change_z(nouv_position_z[ind_noeud]);  
                  ind_noeud++;
                  }
                }
return 1;
  
}

//---------------------------- Lissage Taubin 1995 par Bryan McKenzie-------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::lissage_Taubin1995(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double lambda, double nu, int iter_max)
{
int compt = 0;
int ok = 0;
char mess[200];
  do
  {
   int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
   sprintf(mess, "cycle taubin faite"); affiche(mess);
    compt = compt+1;
    sprintf(mess, "incrémentation compt %d",compt); affiche(mess);
  }while (compt < iter_max);
  return compt;
  sprintf(mess, "fin Taubin1995"); affiche(mess);
}
  

//---------------------------- Varience/Ecart type/ pour McKenzie2016-------------------------------------------------------------------------
void MGOPT_POSTTRAITEMENT::varience_McKenzie2016(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double *vari0,double *vari1,double *vari2)
{
char mess[300]; 
int ecart_type = 0;
double s_i = 0.0;
double phi_im = 0.0;
double ecart_type_i = 0.;
ecart_type = 0;
do
{
  
        TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
        TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
        TPL_LISTE_ENTITE<double> liste_wij;
        LISTE_MG_TRIANGLE::iterator   it_tri;
        int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
        for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
                OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
                normal_f_i.norme();
                liste_normales2.ajouter(normal_f_i);
                //Remplissage de la liste des voisins du triangle i
                TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
                MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
                double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins1;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_1);
                        }
                MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
                double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins2;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_2);
                        }
                MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
                double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
                for (int j = 0;j<nb_voisins3;j++)
                        {
                        MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
                        liste_voisins.ajouter(mgtri_3);
                        }
                liste_voisins.supprimer(mgtri_i);
                int nb_voisins = liste_voisins.get_nb();
                double w_ij = 1./nb_voisins;
                double phi_i_min = 10.;
                s_i = 0.0;
                phi_im = 0.0;
                double *phi_ij = new double[nb_voisins];
                OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
                normal_f_i_mean = normal_f_i;
                OT_VECTEUR_3D eta_i(0.,0.,0.);
                TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
                int j = 0;
                 double un_sur_pi = 1./M_PI;  
                for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
                        {
                        OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
                        //1-Calculer la normale moyenne pour chaque triangle 
                        normal_f_i_mean = normal_f_i_mean + normal_f_j; 
                        //2.1-On calcule l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
                        double prod_scalaire = normal_f_i*normal_f_j;
                        if (prod_scalaire > 1.)
                                { 
                                prod_scalaire = 1.;
                                }
                        if (prod_scalaire < -1.)
                                {
                                prod_scalaire = -1.;
                                }
                        phi_ij[j] = acos(prod_scalaire)*un_sur_pi; 
                        //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
                        if (phi_ij[j] < phi_i_min)
                                {
                                phi_i_min = phi_ij[j]; 
                                eta_i = normal_f_j;
                                }
                        //3.1-On calcule l'angle moyen phi_im
                        phi_im = phi_im + w_ij*phi_ij[j];       
                        j++;
                        }
                normal_f_i_mean.norme();
                j = 0;
                for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
                        {
                        //3.2-Calcul de s_i selon la variance
                        s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
                        j++;
                        }
                delete[] phi_ij;
                }
                ecart_type_i = sqrt(s_i);
                ecart_type = 1;
                }       
while (ecart_type == 0);        
ecart_type = 0;//reinitialisation d'ecart_type pour la prochainepasse
*vari0 = phi_im;
*vari1 = s_i;
*vari2 = ecart_type_i;
//return  &varience[0];  
  
}


//---------------------------- Lissage McKenzie2016 par Bryan Mckenzie 2016-------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::lissage_McKenzie2016(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double lambda, double nu, double epsilon,double sigma,double gamma_,int filtre, int iter_max, int itr_taubin)
{
// Lissage global avec methode de lissage inspiree de Taubin 1995 et Chen 2008 par Bryan McKenzie 2016
  
  
//sauvegarde des differents position de noeuds pour conservation des position
vector<double> coord_x;
vector<double> coord_y;
vector<double> coord_z;
vector<double> vecteur_original;
   LISTE_MG_NOEUD::iterator it_no; 
   int ii=0;
  for (MG_NOEUD* noeud=mg_mai->get_premier_noeud(it_no);noeud!=NULL;noeud=mg_mai->get_suivant_noeud(it_no))
        {
        MG_NOEUD* tri=(MG_NOEUD*)noeud;
        if(tri->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO) 
                {
                coord_x.push_back(noeud->get_x());
                coord_y.push_back(noeud->get_y());
                coord_z.push_back(noeud->get_z());
                vecteur_original.push_back(sqrt(coord_x[ii]*coord_x[ii] + coord_y[ii]*coord_y[ii] + coord_z[ii]*coord_z[ii]));
                ii++;
                }

        }
        

char mess[300]; 
sprintf(mess, "Depart Taubin1995"); affiche(mess);
double v_initial = 0.0;
double v_final = 0.0;
double v_taubin = 0.0;
double v_taubin_prime = 0.0;
double d_initial = 0.0;
double d_final = 0.0;
double d_taubin = 0.0;
double nu_ = nu;
double sigma_ = sigma;//variable utiliser pour sassuré que la pièce ne grossie pas dan taubin_chen 2
int ecart_type = 0;
double s_i = 0.0;
double phi_im = 0.0;
double ecart_type_i = 0.;
int compteur_triangle_retourne =0;
int compteur_total = 0;
int compteur_taubin_initial = 0;
int compteur_chen = 0;
int compteur_taubin_chen = 0;
int compteur_taubin_final = 0;
int compteur_temp = 0; 
int compteur_taubin2D = 0; 
int triangle_final_retourne =0;
int iter_taubin = itr_taubin;


int activationPI = 1; // variable qui active l'adjustement du volume dans taubin1995
int activationPI2 = 1; // variable qui active l'adjustement du volume dans taubin_chen
double condition_de_varience = 1000.0; //variable de sortie pour les prenières iteration taubin pour s'assuré que la surface a une certaine surface lisse
int taubin_premiere_passe = 0;

  /// Calcul du volume initial
  MG_MAILLAGE_OUTILS ot;
  v_initial=ot.get_volume(mg_mai);
  
        
//******************parti taubin pur***********************************//
                
        
//calcule de la varience et l'écart type pour la surface du maillaige au début
double varience_ori0 = 10.0;
double varience_ori1 = 10.0;
double varience_ori2 = 10.0;
varience_McKenzie2016(mg_mai,gest2, &varience_ori0, &varience_ori1, &varience_ori2);
      sprintf(mess,"      Moyenne des normes des triangles de la surface : %f",varience_ori0); affiche(mess);
      sprintf(mess,"      Varience des normes des triangles de la surface : %f",varience_ori1); affiche(mess);
      sprintf(mess,"      Ecart type des normes des triangles de la surface : %f",varience_ori2); affiche(mess);
v_taubin=v_initial; //ajout vincent
do
        {
          v_taubin_prime = v_taubin;//assignation du volume de l'itération précédente pour fair la loi de controle
          int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
                
/// Calcul du volume 
        MG_MAILLAGE_OUTILS ot;
        v_taubin=ot.get_volume(mg_mai);
                
                //implemantion d'un PI por controler le volume
        if (activationPI==1)
          {
                if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))<-0.1)//fine adjustment makes part bigger if new iteration of part is smaller than the preceeding pass and volume is less than 10% bigger than original
                {
                  nu_ =nu_-0.00005;
                }
                if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))>0.1)//fine adjustment makes part smaller if part is much bigger than preceeding pass and the volume is larger than original
                {
                  nu_ = nu_+0.00005;
                }
                if (fabs(v_taubin)>(fabs(v_initial)*1.015))//gross adjustment makes part smaller if part is 3% bigger than original make smaller
                {
                  nu_ =nu_+0.0005;
                }
               if (fabs(v_taubin)<(fabs(v_initial)*0.985))//gross adjustment makes part bigger if it is smaller than the original
                {
                  nu_ =nu_-0.0005;
                } 
          }
         //double un_sur_pi = 1./M_PI;      
        
 
        //code qui vérifie la qualité des different trianges
MAILLEUR_ANALYSE m3d(mg_mai);
m3d.change_eps_angle_retourne(0.03000);
m3d.change_borne(0.1,0.2,0.5);
double qualmin,qualmax,qualmoy;
int tab[7];
m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
int nbretourne=tab[0];
int itri=0;


  if (tab[0]>0)
  {
    compteur_triangle_retourne++;
    compteur_temp--;
    sprintf(mess, "Triangle retourné");
    affiche(mess);
  }
        compteur_taubin_initial++;
        compteur_total++;
        compteur_temp++;
        
        
//calcule de la varience et l'écart type pour la surface du maillaige au début
double varience0 = 0.0;
double varience1 =0.0;
double varience2 =0.0;
varience_McKenzie2016(mg_mai,gest2, &varience0, &varience1, &varience2);

        if (((compteur_temp >= iter_taubin)&&(varience1<0.0275))||(compteur_taubin_initial >(2*iter_taubin)))
        {
        taubin_premiere_passe = 1;
        }
        
        
        }       
while (taubin_premiere_passe == 0);
        

//code d'utilisation de chen pour garder la forme et de taubin pour retirer les triangles retournés
        sprintf(mess, "Début de Chen pour %d iterations avec Taubin qui résoud les triangles retournés",iter_max); affiche(mess);

  do
  {
    
int nbiter=lissage_chen2008(mg_mai,gest2,sigma,gamma_,epsilon,1);

// code de vérifiaction de triangle retourner
MAILLEUR_ANALYSE m3d(mg_mai);
m3d.change_eps_angle_retourne(0.03000);
m3d.change_borne(0.1,0.2,0.5);
double qualmin,qualmax,qualmoy;
int tab[7];
m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
int triangleretourner =tab[0];
      
//code de retourenage de triangle inverser avec taubin 1995
if (triangleretourner >0)       
{
do
        {
          v_taubin_prime = v_taubin;//assignation du volume de l'itération précédente pour fair la loi de controle
          int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
                
/// Calcul du volume 
        MG_MAILLAGE_OUTILS ot;
        v_taubin=ot.get_volume(mg_mai);
                
                //implemantion d'un PI por controler le volume
                if (activationPI2 = 1)
                {
                
                if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))<-0.1)//fine adjustment makes part bigger if new iteration of part is smaller than the preceeding pass and volume is less than 10% bigger than original
                {
                  nu_ =nu_-0.00005;
                }
                
                if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))>0.1)//fine adjustment makes part smaller if part is much bigger than preceeding pass and the volume is larger than original
                {
                  nu_ = nu_+0.00005;
                }

                if (fabs(v_taubin)>(fabs(v_initial)*1.015))//gross adjustment makes part smaller if part is 3% bigger than original make smaller
                {
                  nu_ =nu_+0.0005;
                }
                if (fabs(v_taubin)<(fabs(v_initial)*0.985))//gross adjustment makes part bigger if it is smaller than the original
                {
                  nu_ =nu_-0.0005;
                } 
                }
        
//code qui vérifie la qualité des different trianges
m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
if (tab[0]>0)
  {
    triangleretourner = 1;
    sprintf(mess, "Trianlge retourner pas régler"); affiche(mess);

  } 
  else
  {
    triangleretourner = 0;
    sprintf(mess, "Trianlge retourner régler"); affiche(mess);
  }
compteur_triangle_retourne++;
compteur_taubin_chen++; 
compteur_total++;
}
        
while ((triangleretourner > 0));
        
}       
compteur_total++;
compteur_chen++;


}
while (compteur_chen < iter_max);

sprintf(mess, "Début du retour de volume avec Taubin1995 qui résoud les triangles retournés"); affiche(mess);
// passe final de taubin pour revenir au volume initial
v_taubin=ot.get_volume(mg_mai);

if ((fabs(fabs(v_initial)-fabs(v_taubin)))>(0.02*(fabs(v_initial))))
{
do
        {
          v_taubin_prime = v_taubin;
          int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
                
/// Calcul du volume 
        v_taubin=ot.get_volume(mg_mai);
                
                //implemantion d'un PDF pour controler le volume

                if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))<-0.1)//fine adjustment makes part bigger if new iteration of part is smaller than the preceeding pass and volume is less than 10% bigger than original
                {
                  nu_ =nu_-0.00002;
                  sprintf(mess, "ajustement mineur 1"); affiche(mess);
                                }
                
                if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))>0.1)//fine adjustment makes part smaller if part is much bigger than preceeding pass and the volume is larger than original
                {
                  nu_ = nu_+0.00002;
                  sprintf(mess, "ajustement mineur 2"); affiche(mess);
                }

                if (fabs(v_taubin)>(fabs(v_initial)*1.015))//gross adjustment makes part smaller if part is 1.5% bigger than original make smaller
                {
                  nu_ =nu_+0.0005;
                  sprintf(mess, "ajustement majeur 1"); affiche(mess);
                                }
               if (fabs(v_taubin)<(fabs(v_initial)*0.985))//gross adjustment makes part bigger if it is more than 1.5% smaller than the original
                {
                  nu_ =nu_-0.0005;
                  sprintf(mess, "ajustement majeur 2"); affiche(mess);
                } 
                sprintf(mess, "Volume: %f\n old vol: %f\n nu_: %f\n",v_taubin,v_taubin_prime, nu_); affiche(mess);    
               
               
        //code qui vérifie la qualité des different trianges
MAILLEUR_ANALYSE m3d(mg_mai);
m3d.change_eps_angle_retourne(0.03000);
m3d.change_borne(0.1,0.2,0.5);
double qualmin,qualmax,qualmoy;
int tab[7];
m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
int triangleretourner =tab[0];
if (tab[0]>0)
  {
    compteur_triangle_retourne++;
    sprintf(mess, "Triangle retourné");
    affiche(mess);
    triangle_final_retourne =1;
  }
  else{triangle_final_retourne =0;}
        
        compteur_taubin_final++;
        compteur_total++;
        
        }
        
while ((triangle_final_retourne ==0)&&(fabs(v_initial-v_taubin)>(fabs(0.02*v_initial))));//((compteur_ < 1000));

}

////////////////////////////calcul du mouvement des noeuds par rapport a l'original
/******************code qui calcule la moyenne de pourcentage de déplacement des differents noeuds*/   
   //LISTE_MG_NOEUD::iterator it_no; 
  ii=0;
  int numdepla=0;
  double distot=0;
  double distotx=0;
  double distoty=0;
  double distotz=0;
  double dismin=1e300,dismax=-1e300;
  double disminx=1e300,dismaxx=-1e300;
  double disminy=1e300,dismaxy=-1e300;
  double disminz=1e300,dismaxz=-1e300;
  double xmin=1e300,xmax=-1e300;
  double ymin=1e300,ymax=-1e300;
  double zmin=1e300,zmax=-1e300;
  for (MG_NOEUD* noeud=mg_mai->get_premier_noeud(it_no);noeud!=NULL;noeud=mg_mai->get_suivant_noeud(it_no))
        {
        if(noeud->get_origine()!=MAGIC::ORIGINE::IMPOSE)////////
                {
                double x=noeud->get_x();
                double y=noeud->get_y();
                double z=noeud->get_z();
                if (x<xmin) xmin=x;  
                if (y<ymin) ymin=y;  
                if (z<zmin) zmin=z;  
                if (x>xmax) xmax=x;  
                if (y>ymax) ymax=y;  
                if (z>zmax) zmax=z;  
                double dx=fabs(noeud->get_x()-coord_x[ii]); 
                double dy=fabs(noeud->get_y()-coord_y[ii]); 
                double dz=fabs(noeud->get_z()-coord_z[ii]); 
                double dis=sqrt(dx*dx+dy*dy+dz*dz);  
                if (dx>dismaxx) dismaxx=dx;
                if (dx<disminx) disminx=dx;
                if (dy>dismaxy) dismaxy=dy;
                if (dy<disminy) disminy=dy;
                if (dz>dismaxz) dismaxz=dz;
                if (dz<disminz) disminz=dz;
                if (dis>dismax) dismax=dis;
                if (dis<dismin) dismin=dis;
                distot=distot+dis;  
                distotx=distoty+dx;  
                distoty=distoty+dy;  
                distotz=distotz+dz;  
                numdepla++;
                }
        ii++;
        }
affiche((char*)"");
sprintf(mess, "Déplacement absolu minimum                            en X : %f   en Y : %f  en Z : %f  en amplitude : %f",disminx,disminy,disminz,dismin); affiche(mess);
sprintf(mess, "Déplacement absolu maximum                            en X : %f   en Y : %f  en Z : %f  en amplitude : %f",dismaxx,dismaxy,dismaxz,dismax); affiche(mess);
sprintf(mess, "Déplacement absolu moyen (tous les noeuds)            en X : %f   en Y : %f  en Z : %f  en amplitude : %f",distotx/ii,distoty/ii,distotz/ii,distot/ii); affiche(mess);
sprintf(mess, "Déplacement absolu moyen (les noeuds du non design)   en X : %f   en Y : %f  en Z : %f  en amplitude : %f",distotx/numdepla,distoty/numdepla,distotz/numdepla,distot/numdepla); affiche(mess);
double dx=xmax-xmin;
double dy=ymax-ymin;
double dz=zmax-zmin;
double dmax=dx;
if (dy>dmax) dmax=dy;
if (dz>dmax) dmax=dz;
sprintf(mess, "Déplacement relatif moyen (tous les noeuds)           en X : %f%%  en Y : %f%% en Z : %f%% en amplitude : %f%%",distotx/ii/dx*100.,distoty/ii/dy*100.,distotz/ii/dz*100.,distot/ii/dmax*100.); affiche(mess);
sprintf(mess, "Déplacement relatif moyen (les noeuds du non design)  en X : %f%%  en Y : %f%% en Z : %f%% en amplitude : %f%%",distotx/numdepla/dx*100.,distoty/numdepla/dy*100.,distotz/numdepla/dz*100.,distot/numdepla/dmax*100.); affiche(mess);


affiche((char*)"");

affiche((char*)"");


/// Calcul du volume final
        v_final=ot.get_volume(mg_mai);
        
double delta_vol = v_final-v_initial;
double vol_pourcent = (v_final/v_initial)*100;
double delta_vol_pourcent = vol_pourcent-100;
sprintf(mess, "Volume initial: %f",v_initial);affiche(mess);
sprintf(mess, "Volume final: %f",v_final);affiche(mess);
sprintf(mess, "Volume final pourcentage: %f%%",vol_pourcent);affiche(mess);
sprintf(mess, "Delta Volume: %f",delta_vol);affiche(mess);
sprintf(mess, "Delta Volume pourcentage: %f%%",delta_vol_pourcent);affiche(mess);
sprintf(mess, " Iteration inital Taubin1995: %i",compteur_taubin_initial); affiche(mess);
sprintf(mess, " Iteration Chen2008: %i",compteur_chen); affiche(mess);
sprintf(mess, "         Iteration Taubin1995 dans Chen2008: %i",compteur_taubin_chen); affiche(mess);
sprintf(mess, "         epsilon: %f",epsilon); affiche(mess);
sprintf(mess, "         sigma: %f",sigma); affiche(mess);
sprintf(mess, "         gamma_: %f",gamma_); affiche(mess);
sprintf(mess, "         filtre: %i",filtre); affiche(mess);
sprintf(mess, " Iteration final Taubin1995: %i",compteur_taubin_final); affiche(mess);
sprintf(mess, "         lambda(depart): %f",lambda); affiche(mess);
sprintf(mess, "         nu(depart):%f",nu); affiche(mess);
sprintf(mess, "         nu(fin):%f",nu_); affiche(mess);
sprintf(mess, " Iteration triangle retourné: %i",compteur_triangle_retourne); affiche(mess);
sprintf(mess, " Iteration Total: %i",compteur_total); affiche(mess);
                        
        
//verification de nombre de triange retournée
MAILLEUR_ANALYSE m3d(mg_mai);
double eps_angle_retourne=0.03;
m3d.change_eps_angle_retourne(eps_angle_retourne);
double borne1=0.1;
double borne2=0.2;
double borne3=0.5;
m3d.change_borne(borne1,borne2,borne3);
double qualmin,qualmax,qualmoy;
int tab[7];
m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);              
sprintf(mess,"\n\nANALYSEUR DE MG MAILLAGE 2D"); affiche(mess);
  sprintf(mess,"      Nombre de triangles avec des noeuds fusionnés : %d",tab[6]); affiche(mess); 
  sprintf(mess,"      Nombre de triangles avec 1 voisin manquant : %d",tab[5]); affiche(mess); 
  sprintf(mess,"      critere d'angle limite pour l'inversion : %lf",eps_angle_retourne); affiche(mess); 
  sprintf(mess,"      qualite moyenne des triangles de frontiere : %lf",qualmoy); affiche(mess); 
  sprintf(mess,"      qualite min des triangles de frontiere : %lf",qualmin); affiche(mess); 
  sprintf(mess,"      qualite max des triangles de frontiere : %lf",qualmax); affiche(mess); 
  sprintf(mess,"      nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./mg_mai->get_nb_mg_triangle()); affiche(mess); 
  sprintf(mess,"      nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./mg_mai->get_nb_mg_triangle()); affiche(mess); 
  sprintf(mess,"      nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./mg_mai->get_nb_mg_triangle()); affiche(mess); 
  sprintf(mess,"      nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./mg_mai->get_nb_mg_triangle()); affiche(mess); 
  sprintf(mess,"      nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab[0],tab[0]*100./mg_mai->get_nb_mg_triangle()); affiche(mess);
                
//calcule de la varience et l'écart type pour la surface du maillaige
double varience00 =0.0;
double varience11 =0.0;
double varience22 =0.0;
varience_McKenzie2016(mg_mai,gest2, &varience00, &varience11, &varience22);
      sprintf(mess,"\n      Moyenne des normes des triangles de la surface : %f",varience00); affiche(mess);
      sprintf(mess,"      Varience des normes des triangles de la surface : %f",varience11); affiche(mess);
      sprintf(mess,"      Ecart type des normes des triangles de la surface : %f",varience22); affiche(mess);
                
return compteur_total;
}


//---------------------------- Extract skin -------------------------------------------------------------------------
int MGOPT_POSTTRAITEMENT::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,int &nbpeau,int &nbmaniare,int &nbmanino, int *mai2_id)
{

//etape 0 - On commence par mettre a zero tous les nouveau_numero des triangles et des noeuds du maillage
LISTE_MG_TRIANGLE::iterator   it_tri;
for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
        {
        mgtri->change_nouveau_numero(0);
        mgtri->change_origine(MAGIC::ORIGINE::MAILLEUR_AUTO);
        }
LISTE_MG_NOEUD::iterator   it_noeud;
for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
        {
        mgnoeud->change_nouveau_numero(0);
        }

//etape 1 - On boucle ensuite tous les tetraedres pour ajouter un compteur du nombre de fois qu'un triangle appartient a 1 tetra
LISTE_MG_TETRA::iterator   it_tetra;
for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
        {
        int origine = mgtet->get_origine();
        if (origine==MAGIC::ORIGINE::IMPOSE)
                {
                mgtet->get_triangle1()->change_origine(MAGIC::ORIGINE::IMPOSE);
                mgtet->get_triangle2()->change_origine(MAGIC::ORIGINE::IMPOSE);
                mgtet->get_triangle3()->change_origine(MAGIC::ORIGINE::IMPOSE);
                mgtet->get_triangle4()->change_origine(MAGIC::ORIGINE::IMPOSE);
                }
        if (origine==MAGIC::ORIGINE::OPTIMISE)
                {
                if  (mgtet->get_triangle1()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle1()->change_origine(MAGIC::ORIGINE::OPTIMISE);
                if  (mgtet->get_triangle2()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle2()->change_origine(MAGIC::ORIGINE::OPTIMISE);
                if  (mgtet->get_triangle3()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle3()->change_origine(MAGIC::ORIGINE::OPTIMISE);
                if  (mgtet->get_triangle4()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle4()->change_origine(MAGIC::ORIGINE::OPTIMISE);
                }
        
        if (((origine == MAGIC::ORIGINE::OPTIMISE) && (etat[(MAGIC::ORIGINE::OPTIMISE-1000)/10]==1) ) ||
           ((origine == MAGIC::ORIGINE::IMPOSE) && (etat[(MAGIC::ORIGINE::IMPOSE-1000)/10]==1) ) ||
        ((origine == MAGIC::ORIGINE::MAILLEUR_AUTO) && (etat[(MAGIC::ORIGINE::MAILLEUR_AUTO-1000)/10]==1) ) ||
        ((origine == MAGIC::ORIGINE::TRIANGULATION) && (etat[(MAGIC::ORIGINE::TRIANGULATION-1000)/10]==1) ) ||
        ((origine == MAGIC::ORIGINE::MODIFICATION) && (etat[(MAGIC::ORIGINE::MODIFICATION-1000)/10]==1) ) ||
        ((origine == MAGIC::ORIGINE::DUPLIQUER) && (etat[(MAGIC::ORIGINE::DUPLIQUER-1000)/10]==1) ) )

                {
                int num1 = mgtet->get_triangle1()->get_nouveau_numero();
                int num2 = mgtet->get_triangle2()->get_nouveau_numero();
                int num3 = mgtet->get_triangle3()->get_nouveau_numero();
                int num4 = mgtet->get_triangle4()->get_nouveau_numero();
                num1++;
                num2++;
                num3++;
                num4++;
                mgtet->get_triangle1()->change_nouveau_numero(num1);
                mgtet->get_triangle2()->change_nouveau_numero(num2);
                mgtet->get_triangle3()->change_nouveau_numero(num3);
                mgtet->get_triangle4()->change_nouveau_numero(num4);
                }
        }

//etape 2 - On boucle l'ensemble des triangles identifies a l'etape 1 pour identifier les noeuds leur appartenant
for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
        {
        int num = mgtri->get_nouveau_numero();
        if (num == 1) 
                {
                MG_NOEUD* noeud1 = mgtri->get_noeud1();
                MG_NOEUD* noeud2 = mgtri->get_noeud2();
                MG_NOEUD* noeud3 = mgtri->get_noeud3();
                noeud1->change_nouveau_numero(1);
                noeud2->change_nouveau_numero(1);
                noeud3->change_nouveau_numero(1);
                if (mgtri->get_origine()==MAGIC::ORIGINE::IMPOSE) 
                        {
                        noeud1->change_origine(MAGIC::ORIGINE::IMPOSE);
                        noeud2->change_origine(MAGIC::ORIGINE::IMPOSE);
                        noeud3->change_origine(MAGIC::ORIGINE::IMPOSE);
                        }
                }
        }

//etape 3 - On cree un nouveau maillage pour la peau
gest2.vide();
MG_GEOMETRIE* geo=new MG_GEOMETRIE((char*)"VIRTUEL",(char*)"VIRTUEL");
geo->change_valeur_unite(mg_mai->get_mg_geometrie()->get_valeur_unite());
gest2.ajouter_mg_geometrie(geo);
MG_FACE_ELEMENT* face=new MG_FACE_ELEMENT();
geo->ajouter_mg_face(face);
MG_MAILLAGE* mai2 = new MG_MAILLAGE(geo); 
gest2.ajouter_mg_maillage(mai2);

//etape 4 - On boucle l'ensemble des noeuds identifies a l'etape 2 pour les recreer dans le second maillage
for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
        {
        int num = mgnoeud->get_nouveau_numero();
        if (num == 1) 
                {
                double x = mgnoeud->get_x();
                double y = mgnoeud->get_y();
                double z = mgnoeud->get_z();
                int origine=MAGIC::ORIGINE::TRIANGULATION;
                int nb_tri=mgnoeud->get_lien_triangle()->get_nb();
                bool tri2d=false,tri3d=false;
                for (int i=0;i<nb_tri;i++)
                    {
                    if (mgnoeud->get_lien_triangle()->get(i)->get_nouveau_numero()==1)
                      {
                      if (mgnoeud->get_lien_triangle()->get(i)->get_lien_topologie()->get_dimension()==2) tri2d=true;
                      if (mgnoeud->get_lien_triangle()->get(i)->get_lien_topologie()->get_dimension()==3) tri3d=true;
                      }
                    }
                if ((tri2d)&&(tri3d)) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
                if (mgnoeud->get_lien_topologie()->get_dimension()==1) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
                if (mgnoeud->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
                MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
                mai2->ajouter_mg_noeud(noeud1);
                mgnoeud->change_nouveau_numero(noeud1->get_id());
                noeud1->change_nouveau_numero(mgnoeud->get_id());
                }
        }

// etape 5 - On boucle l'ensemble des triangles identifies a l'etape 1 pour les recreer dans le maillage 2
for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
        {
        int originetet=mgtet->get_origine();
        if (((originetet == MAGIC::ORIGINE::OPTIMISE) && (etat[(MAGIC::ORIGINE::OPTIMISE-1000)/10]==1) ) ||
           ((originetet == MAGIC::ORIGINE::IMPOSE) && (etat[(MAGIC::ORIGINE::IMPOSE-1000)/10]==1) ) ||
        ((originetet == MAGIC::ORIGINE::MAILLEUR_AUTO) && (etat[(MAGIC::ORIGINE::MAILLEUR_AUTO-1000)/10]==1) ) ||
        ((originetet == MAGIC::ORIGINE::TRIANGULATION) && (etat[(MAGIC::ORIGINE::TRIANGULATION-1000)/10]==1) ) ||
        ((originetet == MAGIC::ORIGINE::MODIFICATION) && (etat[(MAGIC::ORIGINE::MODIFICATION-1000)/10]==1) ) ||
        ((originetet == MAGIC::ORIGINE::DUPLIQUER) && (etat[(MAGIC::ORIGINE::DUPLIQUER-1000)/10]==1) ) )
                {
                MG_NOEUD* noeud1 = mgtet->get_noeud1();
                MG_NOEUD* noeud2 = mgtet->get_noeud2();
                MG_NOEUD* noeud3 = mgtet->get_noeud3();
                MG_NOEUD* noeud4 = mgtet->get_noeud4();
                MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
                MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
                MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
                MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
                MG_TRIANGLE* tri1=mgtet->get_triangle1();
                MG_TRIANGLE* tri2=mgtet->get_triangle2();
                MG_TRIANGLE* tri3=mgtet->get_triangle3();
                MG_TRIANGLE* tri4=mgtet->get_triangle4();
                if (tri1->get_nouveau_numero()==1)
                        {
                        int origine=MAGIC::ORIGINE::TRIANGULATION;
                        if (tri1->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
                        MG_FACE* topo=NULL;
                        if (tri1->get_lien_topologie()->get_dimension()==2) topo=face;
                        MG_TRIANGLE_PEAU* tripeau = insere_triangle(topo,node1,node2,node3,mai2,origine);
                        tripeau->change_nouveau_numero(0);
                        tri1->change_nouveau_numero(0);
                        }
                if (tri2->get_nouveau_numero()==1)
                        {
                        int origine=MAGIC::ORIGINE::TRIANGULATION;
                        if (tri2->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
                        MG_FACE* topo=NULL;
                        if (tri2->get_lien_topologie()->get_dimension()==2) topo=face;
                        MG_TRIANGLE_PEAU* tripeau = insere_triangle(topo,node1,node4,node2,mai2,origine);
                        tripeau->change_nouveau_numero(0);
                        tri2->change_nouveau_numero(0);
                        }
                if (tri3->get_nouveau_numero()==1)
                        {
                        int origine=MAGIC::ORIGINE::TRIANGULATION;
                        MG_FACE* topo=NULL;
                        if (tri3->get_lien_topologie()->get_dimension()==2) topo=face;
                        if (tri3->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
                        MG_TRIANGLE_PEAU* tripeau = insere_triangle(topo,node2,node4,node3,mai2,origine);
                        tripeau->change_nouveau_numero(0);
                        tri3->change_nouveau_numero(0);
                        }
                if (tri4->get_nouveau_numero()==1)
                        {
                        int origine=MAGIC::ORIGINE::TRIANGULATION;
                        MG_FACE* topo=NULL;
                        if (tri4->get_lien_topologie()->get_dimension()==2) topo=face;
                        if (tri4->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
                        MG_TRIANGLE_PEAU* tripeau = insere_triangle(topo,node1,node3,node4,mai2,origine);
                        tripeau->change_nouveau_numero(0);
                        tri4->change_nouveau_numero(0);
                        }
                }
}
// etape 6 recherche des vosins
for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
        {
        MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
        MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
        MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
        MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
        tripeau->change_voisin1(voisin1);
        tripeau->change_voisin2(voisin2);
        tripeau->change_voisin3(voisin3);
        }
// etape 7 - On recherche les peaux
int fin;
do
        {
        fin=1;
        for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
                {
                MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
                if (tripeau->get_nouveau_numero()==0)
                        {
                        fin=0;
                        std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
                        lst_peau.push_back(peau);
                        tripeau->change_nouveau_numero(1);
                        peau->push_back(tripeau);
                        determine_peau(peau);
                        }
                }
        }
while (fin==0);
nbpeau=lst_peau.size();
//test de manifold 
LISTE_MG_SEGMENT::iterator itseg;
for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
        seg->change_nouveau_numero(0);
TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeuddepuisarete;
nbpeau=lst_peau.size();
for (int i=0;i<nbpeau;i++)
        {
        int nbtri=lst_peau[i]->size();
        for (int j=0;j<nbtri;j++)
                {
                MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
                tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
                tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
                tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
                if (tripeau->get_segment1()->get_nouveau_numero()>2) 
                        nonmanifoldarete.ajouter(tripeau->get_segment1());
                if (tripeau->get_segment2()->get_nouveau_numero()>2) 
                        nonmanifoldarete.ajouter(tripeau->get_segment2());
                if (tripeau->get_segment3()->get_nouveau_numero()>2) 
                        nonmanifoldarete.ajouter(tripeau->get_segment3());
                }
        }

int nbnonmanifoldarete=nonmanifoldarete.get_nb();
for (int i=0;i<nbnonmanifoldarete;i++)
        {
        MG_SEGMENT* seg=nonmanifoldarete.get(i);
        MG_NOEUD* n1=seg->get_noeud1();
        MG_NOEUD* n2=seg->get_noeud2();
        nonmanifoldnoeuddepuisarete.ajouter(n1);
        nonmanifoldnoeuddepuisarete.ajouter(n2);
        }
LISTE_MG_NOEUD::iterator itnoeud;
for (MG_NOEUD* no=mai2->get_premier_noeud(itnoeud);no!=NULL;no=mai2->get_suivant_noeud(itnoeud))
        {
        if (nonmanifoldnoeuddepuisarete.existe(no)) continue;
        if (est_non_manifold(no)) nonmanifoldnoeud.ajouter(no); 
        }
        
nbmaniare = nonmanifoldarete.get_nb(); 
nbmanino  = nonmanifoldnoeud.get_nb(); 

//etape 8 - Traitement des cas de non manifold
//non manifold par arete
for (int i=0;i<nbnonmanifoldarete;i++)
        {
        MG_SEGMENT* segment=nonmanifoldarete.get(i);
        MG_NOEUD* noeud1 = mg_mai->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
        MG_NOEUD* noeud2 = mg_mai->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
        int nb_tetra = noeud1->get_lien_tetra()->get_nb();
        for (int j=0;j<nb_tetra;j++)
                {
                MG_TETRA* mgtet =noeud1->get_lien_tetra()->get(j);
                int originetet=mgtet->get_origine();
                if (originetet == MAGIC::ORIGINE::MAILLEUR_AUTO)
                        {
                        //On reactive le tetra si l'autre noeud du segment lui appartient aussi
                        MG_NOEUD* no1 = mgtet->get_noeud1();
                        MG_NOEUD* no2 = mgtet->get_noeud2();
                        MG_NOEUD* no3 = mgtet->get_noeud3();
                        MG_NOEUD* no4 = mgtet->get_noeud4();
                        if((no1 == noeud2) ||(no2 == noeud2) || (no3 == noeud2) || (no4 == noeud2))
                                mgtet->change_origine(MAGIC::ORIGINE::OPTIMISE);
                        }
                }       
        }

//non manifold par noeud
int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
for (int i=0;i<nbnonmanifoldnoeud;i++)
        {
        MG_NOEUD* noeud=mg_mai->get_mg_noeudid(nonmanifoldnoeud.get(i)->get_nouveau_numero());
        int nb_tetra = noeud->get_lien_tetra()->get_nb();
        for (int j = 0;j<nb_tetra;j++)
                { 
                MG_TETRA* mgtet =noeud->get_lien_tetra()->get(j);
                int originetet=mgtet->get_origine();
                if (originetet == MAGIC::ORIGINE::MAILLEUR_AUTO)
                        {
                        mgtet->change_origine(MAGIC::ORIGINE::OPTIMISE);        
                        }
                }
        }
*mai2_id = mai2->get_id();
if ((nbnonmanifoldarete != 0) || (nbnonmanifoldnoeud != 0))
        {
        for (int i=0;i<lst_peau.size();i++)
                {
                delete lst_peau[i];
                }
        lst_peau.clear();
        gest2.supprimer_mg_maillageid(*mai2_id);
        return 0;
        }

return 1;
}

//---------------------------- Determine peau -------------------------------------------------------------------------
void MGOPT_POSTTRAITEMENT::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
{
int num=0;
while (num!=peau->size())
        {
        MG_TRIANGLE_PEAU* p=(*peau)[num];
        if (p->get_voisin1()->get_nouveau_numero()==0)
                {
                peau->push_back(p->get_voisin1());
                p->get_voisin1()->change_nouveau_numero(1);
                }               
        if (p->get_voisin2()->get_nouveau_numero()==0)
                {
                peau->push_back(p->get_voisin2());
                p->get_voisin2()->change_nouveau_numero(1);
                }
        if (p->get_voisin3()->get_nouveau_numero()==0)
                {
                peau->push_back(p->get_voisin3());
                p->get_voisin3()->change_nouveau_numero(1);
                }
        num++;
        }
}

//---------------------------- Insere triangle -------------------------------------------------------------------------
MG_TRIANGLE_PEAU* MGOPT_POSTTRAITEMENT::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3,MG_MAILLAGE* mg_maillage,int origine)
{
MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
if (mgsegment1==NULL) 
        mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
if (mgsegment2==NULL) 
        mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
if (mgsegment3==NULL) 
        mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
mg_maillage->ajouter_mg_triangle(mtriangle);
return mtriangle;
}

//---------------------------- Recherche voisin -------------------------------------------------------------------------
MG_TRIANGLE_PEAU* MGOPT_POSTTRAITEMENT::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
{
MG_TRIANGLE_PEAU* trisol=NULL;
double angleref=4.*M_PI;
int nb1=mg_noeud1->get_lien_triangle()->get_nb();
int nb2=mg_noeud2->get_lien_triangle()->get_nb();
for (int i=0;i<nb1;i++)
for (int j=0;j<nb2;j++)
if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j)) 
        if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
                        {
                        double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
                        if (angle<angleref)
                                {
                                angleref=angle;
                                trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
                                }
                        } 
return trisol;
}

int MGOPT_POSTTRAITEMENT::est_non_manifold(MG_NOEUD* no)
{
static int compteur=0;
compteur++;
int nb_tri=no->get_lien_triangle()->get_nb();
TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
for (int i=0;i<nb_tri;i++)
        lst.ajouter(no->get_lien_triangle()->get(i));
MG_TRIANGLE_PEAU* premier_triangle=(MG_TRIANGLE_PEAU*)lst.get(0);
lst.supprimer(premier_triangle);
MG_TRIANGLE_PEAU* tricour=(MG_TRIANGLE_PEAU*)premier_triangle;
int i=0;
do
        {
        if (lst.existe(tricour->get_voisin1()) || ((tricour->get_voisin1()==premier_triangle)&& (i>1))) 
                {
                tricour=tricour->get_voisin1();
                lst.supprimer(tricour);         
                }
        else if (lst.existe(tricour->get_voisin2()) || ((tricour->get_voisin2()==premier_triangle)&& (i>1)))    
                {
                tricour=tricour->get_voisin2();
                lst.supprimer(tricour);         
                }
        else if (lst.existe(tricour->get_voisin3()) || ((tricour->get_voisin3()==premier_triangle)&& (i>1)))
                {
                tricour=tricour->get_voisin3();
                lst.supprimer(tricour);         
                }
        i++;
        }
while (tricour!=premier_triangle);
if (lst.get_nb()>0) 
        return 1;
return 0;
}

//---------------------------- Calcule angle -------------------------------------------------------------------------
double MGOPT_POSTTRAITEMENT::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
{
    MG_NOEUD* noeud1=ft1->get_noeud1();
    MG_NOEUD* noeud2=ft1->get_noeud2();
    MG_NOEUD* noeud3=ft1->get_noeud3();
    MG_NOEUD* noeud4=ft2->get_noeud1();
    MG_NOEUD* noeud5=ft2->get_noeud2();
    MG_NOEUD* noeud6=ft2->get_noeud3();
    double angle=get_angle_par_noeud<MG_NOEUD*>(noeud1,noeud2,noeud3,noeud4,noeud5,noeud6);
    return angle;
}

//---------------------------- Ponderation gaussian -------------------------------------------------------------------------
double MGOPT_POSTTRAITEMENT::ponderation_gaussian(double s,double sigma)
{
      double w_s;
      w_s = exp(-(s*s)/(2.*sigma*sigma));
      return w_s;
}

//---------------------------- Ponderation Laplacian -------------------------------------------------------------------------
double MGOPT_POSTTRAITEMENT::ponderation_laplacian(double s,double sigma)
{
      double w_s;
      w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
      return w_s;
}

//---------------------------- Ponderation elfallahford -------------------------------------------------------------------------
double MGOPT_POSTTRAITEMENT::ponderation_elfallahford(double s,double sigma)
{
      double w_s;
      w_s = 1./sqrt(1+pow((s/sigma),2));
      return w_s;
}
Revision:	805
Committed:	Mon Jun 20 16:37:35 2016 UTC (8 years, 10 months ago) by francois
File size:	102368 byte(s)
Log Message:	Optimisation de la méthode de lissage McKenzie2016 suppression de la méthode inconnue nan2015 dans le posttraitement des optimisation topologique