mailleur3d_optimisation.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mailleur3d_optimisation.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:55 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "gestionversion.h"
#include "mailleur3d_optimisation.h"
#include "ot_mathematique.h"
#include "m3d_tetra.h"
#include "m3d_triangle.h"
#include "m3d_noeud.h"
#include "mg_maillage.h"
#include <math.h>
 
 
 
 
MAILLEUR3D_OPTIMISATION::MAILLEUR3D_OPTIMISATION(MG_MAILLAGE* mgmai,int niv):MAILLEUR(false),mg_maillage(mgmai),niveau_optimisation(niv)
{
o3d_data();
o3d_data2();
LISTE_MG_TETRA::iterator it;
TPL_MAP_ENTITE<M3D_TETRA*> listetet3d;
for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet!=NULL;tet=mg_maillage->get_suivant_tetra(it))
    if (tet->get_type_entite()!=MAGIC::TYPE_ENTITE::IDM3D_TETRA)    
        {
        M3D_TETRA* tet3d=new M3D_TETRA(tet->get_id(),tet->get_lien_topologie(),tet->get_noeud1(),tet->get_noeud2(),tet->get_noeud3(),tet->get_noeud4(),tet->get_triangle1(),tet->get_triangle2(),tet->get_triangle3(),tet->get_triangle4(),tet->get_origine());
        listetet3d.ajouter(tet3d);    
        }
 
TPL_MAP_ENTITE<M3D_TETRA*>::ITERATEUR it2;
 
for (M3D_TETRA* tet=listetet3d.get_premier(it2);tet!=NULL;tet=listetet3d.get_suivant(it2))
  {
  mg_maillage->supprimer_mg_tetraid(tet->get_id());
  mg_maillage->ajouter_mg_tetra(tet);  
  }
}
 
MAILLEUR3D_OPTIMISATION::~MAILLEUR3D_OPTIMISATION()
{
}
 
 
 
double MAILLEUR3D_OPTIMISATION::get_volume(MG_TETRA* tet)
{
    double *xyz1=tet->get_noeud1()->get_coord();
    double *xyz2=tet->get_noeud2()->get_coord();
    double *xyz3=tet->get_noeud3()->get_coord();
    double *xyz4=tet->get_noeud4()->get_coord();
    OT_VECTEUR_3D ab(xyz1,xyz2);
    OT_VECTEUR_3D ac(xyz1,xyz3);
    OT_VECTEUR_3D ad(xyz1,xyz4);
    double vol=(ab&ac)*ad;
    vol=vol/6.;
    return vol;
}
 
double MAILLEUR3D_OPTIMISATION::get_volume(double *xyz1,double *xyz2,double *xyz3,double *xyz4)
{
    OT_VECTEUR_3D ab(xyz1,xyz2);
    OT_VECTEUR_3D ac(xyz1,xyz3);
    OT_VECTEUR_3D ad(xyz1,xyz4);
    double vol=(ab&ac)*ad;
    vol=vol/6.;
    return vol;
}
 
void MAILLEUR3D_OPTIMISATION::optimise(MG_VOLUME* mgvol)
{
    char mess[255];
    sprintf(mess,"     Niveau d'optimisation : %d",niveau_optimisation);
    if (affichageactif==1) affiche(mess);
    int nbaoptimiser;
    int nbaoptimiserapres=mg_maillage->get_nb_mg_tetra();
    do
    {
        nbaoptimiser=nbaoptimiserapres;
        optimise(mgvol,nbaoptimiserapres);
        if (nbaoptimiserapres==0) 
        {
          nbaoptimiser=0;
          sprintf(mess,"     Tout est optimise");
          if (affichageactif==1) affiche(mess);
        }
          
    }
    while (nbaoptimiserapres!=nbaoptimiser);
}
 
 
 
void MAILLEUR3D_OPTIMISATION::optimise(MG_VOLUME* mgvol,int& nbmauvais)
{
    int passe=1;
    //passe++;
    int nbtet=0;
    if (mgvol!=NULL)
    {
      TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
      for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
      {
        if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
        M3D_TETRA* mtet=(M3D_TETRA*)tet;
        double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
        mtet->change_qualite(qual);
        if (qual<0.1*niveau_optimisation) nbtet++;
        ajouter_ordre_tetra(mtet,0);
      }
    }
   else
   {
    LISTE_MG_TETRA::iterator it;
    for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
    {
        if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
        M3D_TETRA* mtet=(M3D_TETRA*)tet;
        double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
        mtet->change_qualite(qual);
        if (qual<0.1*niveau_optimisation) nbtet++;
        ajouter_ordre_tetra(mtet,0);
    }
   }
    
    if (nbtet==0) return;
    for (int phase=0;phase<2;phase++)
    {
        char mess[255];
        if (phase==0) sprintf(mess,"     Phase %d : %d a optimiser",2*passe-1+phase,nbtet);
        if (phase==1) sprintf(mess,"     Phase %d",2*passe-1+phase);
        if (affichageactif==1) affiche(mess);
        while (lst_tetra[phase].size()>0)
        {
            ORDRE_TETRA::iterator i=lst_tetra[phase].begin();
            M3D_TETRA* tet=(*i).second;
            supprimer_ordre_tetra(tet);
            double crit_avant=tet->get_qualite();
            if (crit_avant<0.1*niveau_optimisation)
            {
                double crit[4]={-1.,-1.,-1.,-1.};
                MG_NOEUD* no[4];
                no[0]=(MG_NOEUD*)tet->get_noeud1();
                no[1]=(MG_NOEUD*)tet->get_noeud2();
                no[2]=(MG_NOEUD*)tet->get_noeud3();
                no[3]=(MG_NOEUD*)tet->get_noeud4();
                double x[4],y[4],z[4];
                bouge_point(mgvol,no[0],crit[0],x[0],y[0],z[0]);
                bouge_point(mgvol,no[1],crit[1],x[1],y[1],z[1]);
                bouge_point(mgvol,no[2],crit[2],x[2],y[2],z[2]);
                bouge_point(mgvol,no[3],crit[3],x[3],y[3],z[3]);
                std::vector<double> vecteur_crit(crit,crit+4);
                std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
                double crit_opt=*it1;
                if (crit_opt>crit_avant)
                {
                    int num=it1-vecteur_crit.begin();
                    no[num]->change_x(x[num]);
                    no[num]->change_y(y[num]);
                    no[num]->change_z(z[num]);
                    int nb=no[num]->get_lien_tetra()->get_nb();
                    for (int j=0;j<nb;j++)
                    {
                        M3D_TETRA* mtet=(M3D_TETRA*)no[num]->get_lien_tetra()->get(j);
                        double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
                        mtet->change_qualite(qual);
                    }
                    crit_avant=crit_opt;
                }
            }
            
            int coquille_non_optimise=1;
            if (crit_avant<0.1*niveau_optimisation)
            {
                COQUILLE coque[6];
                double crit[6]={-1,-1,-1,-1,-1,-1};
                remaille_coquille(tet->get_noeud1(),tet->get_noeud2(),crit[0],coque[0]);
                remaille_coquille(tet->get_noeud1(),tet->get_noeud3(),crit[1],coque[1]);
                remaille_coquille(tet->get_noeud1(),tet->get_noeud4(),crit[2],coque[2]);
                remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[3],coque[3]);
                remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[4],coque[4]);
                remaille_coquille(tet->get_noeud3(),tet->get_noeud4(),crit[5],coque[5]);
                std::vector<double> vecteur_crit(crit,crit+6);
                std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
                double crit_opt=*it1;
                if (crit_opt>crit_avant)
                {
                    coquille_non_optimise=0;
                    int num=it1-vecteur_crit.begin();
                    for (int i=0;i<coque[num].taille;i++)
                    {
                        M3D_TETRA *tettmp=(M3D_TETRA*)coque[num].tet[i];
                        if (tet!=tettmp) supprimer_ordre_tetra(tettmp);
                        mg_maillage->supprimer_mg_tetraid(coque[num].tet[i]->get_id());
                    }
                    for (int i=0;i<2*coque[num].taille-4;i++)
                    {
                        MG_NOEUD* noeud1=coque[num].new_tetra[4*i];
                        MG_NOEUD* noeud2=coque[num].new_tetra[4*i+1];
                        MG_NOEUD* noeud3=coque[num].new_tetra[4*i+2];
                        MG_NOEUD* noeud4=coque[num].new_tetra[4*i+3];
                        MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
                        MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
                        MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
                        MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
                        if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
                        if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
                        if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
                        if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
                        M3D_TETRA* mtet=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4,MAGIC::ORIGINE::MAILLEUR_AUTO);
                        mg_maillage->ajouter_mg_tetra(mtet);
                        double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
                        mtet->change_qualite(qual);
                        if (phase==0)
                            ajouter_ordre_tetra(mtet,1);
                    }
 
 
                }
 
            }
            if ((phase==0) && (coquille_non_optimise==1))
                ajouter_ordre_tetra(tet,1);
        }
 
    }
    int nb=0;
    //for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
    if (mgvol!=NULL)
    {
      TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
      for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
        {
        if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
        M3D_TETRA* mtet=(M3D_TETRA*)tet;
        if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
        }
    }
    else
        {
        LISTE_MG_TETRA::iterator it;
        for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
                {
                if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
                M3D_TETRA* mtet=(M3D_TETRA*)tet;
                if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
                }
        }
 
    nbmauvais=nb;
 
}
 
 
 
int MAILLEUR3D_OPTIMISATION::bouge_point(MG_VOLUME* mgvol,MG_NOEUD* noeud,double& crit,double& x,double& y, double& z)
{
    if (noeud->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO) return 0;
 //   if (noeud->get_type_entite()!=IDM3D_NOEUD) return 0;
   
     if (mgvol!=NULL) 
      if (noeud->get_lien_topologie()!=mgvol) return 0;
    if (mgvol==NULL)
        if (noeud->get_type_entite()==MAGIC::TYPE_ENTITE::IDM3D_NOEUD)
            if (((M3D_NOEUD*)noeud)->get_etat()==MAGIC::MAILLEURFRONTALETAT::INACTIF) return 0;
    double xopt=0.;
    double yopt=0.;
    double zopt=0.;
    double qual_dep=1.;
    int nb_tet=noeud->get_lien_tetra()->get_nb();
    double tab_coord[3000];
    double gamma=0;
    for (int i=0;i<nb_tet;i++)
    {
        M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
        qual_dep=std::min(qual_dep,tet->get_qualite());
        MG_NOEUD* no1=tet->get_noeud1();
        MG_NOEUD* no2=tet->get_noeud2();
        MG_NOEUD* no3=tet->get_noeud3();
        MG_NOEUD* no4=tet->get_noeud4();
        MG_NOEUD *noeud1,*noeud2,*noeud3;
        if (no1==noeud)
        {
            noeud1=no2;
            noeud2=no4;
            noeud3=no3;
        }
        if (no2==noeud)
        {
            noeud1=no1;
            noeud2=no3;
            noeud3=no4;
        }
        if (no3==noeud)
        {
            noeud1=no1;
            noeud2=no4;
            noeud3=no2;
        }
        if (no4==noeud)
        {
            
            noeud1=no1;
            noeud2=no2;
            noeud3=no3;
        }
        double xyzi[3];
        double *xyz1=noeud1->get_coord();
        double *xyz2=noeud2->get_coord();
        double *xyz3=noeud3->get_coord();
        xyzi[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
        xyzi[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
        xyzi[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
        OT_VECTEUR_3D v12(xyz1,xyz2);
        OT_VECTEUR_3D v13(xyz1,xyz3);
        OT_VECTEUR_3D v23(xyz2,xyz3);
        OT_VECTEUR_3D normal=v12&v13;
        double perimetre=v12.get_longueur()+v13.get_longueur()+v23.get_longueur();
        double hauteur = (perimetre/3.) * 0.8 ;
        normal.norme();
        tab_coord[3*i]=xyzi[0]+normal.get_x()*hauteur;
        tab_coord[3*i+1]=xyzi[1]+normal.get_y()*hauteur;
        tab_coord[3*i+2]=xyzi[2]+normal.get_z()*hauteur;
        if (tet->get_qualite()> 1e-10) gamma=gamma+1./tet->get_qualite()/tet->get_qualite();
    }
    gamma=1./gamma;
    for (int i=0;i<nb_tet;i++)
    {
        M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
        if (tet->get_qualite()> 1e-10)
        {
            double alpha=gamma/tet->get_qualite()/tet->get_qualite();
            xopt=xopt+alpha*tab_coord[3*i];
            yopt=yopt+alpha*tab_coord[3*i+1];
            zopt=zopt+alpha*tab_coord[3*i+2];
        }
 
    }
    double delta[3]={xopt-noeud->get_x(),yopt-noeud->get_y(),zopt-noeud->get_z()};
    double bmin=0.,bmax=1.;
    double vieuxx=noeud->get_x();
    double vieuxy=noeud->get_y();
    double vieuxz=noeud->get_z();
    double qualini=0.;
    for (int iteration=0;iteration<5;iteration++)
    {
        double alpha=(bmin+bmax)*0.5;
        noeud->change_x(vieuxx+alpha*delta[0]);
        noeud->change_y(vieuxy+alpha*delta[1]);
        noeud->change_z(vieuxz+alpha*delta[2]);
        double qualcour=1.;
        for (int i=0;i<nb_tet;i++)
        {
            M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
            double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
            qualcour=std::min(qualcour,qual);
        }
        double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
        noeud->change_x(vieuxx+alpha_eps*delta[0]);
        noeud->change_y(vieuxy+alpha_eps*delta[1]);
        noeud->change_z(vieuxz+alpha_eps*delta[2]);
        double qualcour_eps=1.;
        for (int i=0;i<nb_tet;i++)
        {
            M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
            double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
            qualcour_eps=std::min(qualcour_eps,qual);
        }
        if (qualcour>qualcour_eps) bmin =alpha;
        else bmax=alpha;
        qualini=std::max(qualini,qualcour);
    }
    noeud->change_x(vieuxx);
    noeud->change_y(vieuxy);
    noeud->change_z(vieuxz);
    if (qualini>qual_dep)
    {
        x=vieuxx+(bmin+bmax)*0.5*delta[0];
        y=vieuxy+(bmin+bmax)*0.5*delta[1];
        z=vieuxz+(bmin+bmax)*0.5*delta[2];
        crit=qualini;
        return 1;
    }
        
    return 0;
}
 
 
 
void MAILLEUR3D_OPTIMISATION::remaille_coquille(MG_NOEUD* noeud1,MG_NOEUD* noeud2, double&  crit, COQUILLE& coque)
{
    MG_NOEUD* noeud;
    if (noeud1->get_id()>noeud2->get_id()) noeud=noeud2;
    else noeud=noeud1;
    int nb=noeud->get_lien_petit_segment()->get_nb();
    for (int i=0;i<nb;i++)
    {
        MG_SEGMENT* seg=noeud->get_lien_petit_segment()->get(i);
        if  ( ((seg->get_noeud1()==noeud1) && (seg->get_noeud2()==noeud2)) || ((seg->get_noeud1()==noeud2) && (seg->get_noeud2()==noeud1)))
        {
          //if (mg_maillage->get_mg_geometrie()==NULL)
          //{
            if (seg->get_dimension_topo_null()==2)
              {
                coque.taille=0;
                crit=0.;
                return;
              }
          //}
          //else 
          if (seg->get_lien_topologie()!=NULL) if ((seg->get_lien_topologie()->get_dimension()!=3) || (seg->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO))
            {
                coque.taille=0;
                crit=0.;
                return;
            }
        }
    }
    int nb1=noeud1->get_lien_tetra()->get_nb();
    int nb2=noeud2->get_lien_tetra()->get_nb();
    MG_TETRA* coq[100];
    int nb_coq=0;
    for (int i=0;i<nb1;i++)
        for (int j=0;j<nb2;j++)
            if (noeud1->get_lien_tetra()->get(i)==noeud2->get_lien_tetra()->get(j))
            {
                coq[nb_coq]=noeud1->get_lien_tetra()->get(i);
                nb_coq++;
            }
    if ((nb_coq<4) || (nb_coq>10))
    {
        coque.taille=0;
        crit=0.;
        return;
    }
    MG_NOEUD* tet_noeud[50];
    int nb_tet_noeud=0;
    coque.taille=nb_coq;
    for (int i=0;i<nb_coq;i++)
    {
        coque.tet[i]=coq[i];
        
        if ((coq[i]->get_noeud1()!=noeud1) && (coq[i]->get_noeud1()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud1();
        if ((coq[i]->get_noeud2()!=noeud1) && (coq[i]->get_noeud2()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud2();
        if ((coq[i]->get_noeud3()!=noeud1) && (coq[i]->get_noeud3()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud3();
        if ((coq[i]->get_noeud4()!=noeud1) && (coq[i]->get_noeud4()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud4();
        coque.volume=coque.volume+fabs(get_volume(coq[i]));
    }
    MG_NOEUD* polygone[20];
    int nb_poly=2;
    polygone[0]=tet_noeud[0];
    tet_noeud[0]=NULL;
    polygone[1]=tet_noeud[1];
    tet_noeud[1]=NULL;
    while (nb_poly!=nb_coq+1)
    {
        for (int j=0;j<nb_coq;j++)
        {
            if (tet_noeud[2*j]==polygone[nb_poly-1])
            {
                polygone[nb_poly++]=tet_noeud[2*j+1];
                tet_noeud[2*j]=NULL;
                tet_noeud[2*j+1]=NULL;
            }
            if (tet_noeud[2*j+1]==polygone[nb_poly-1])
            {
                polygone[nb_poly++]=tet_noeud[2*j];
                tet_noeud[2*j]=NULL;
                tet_noeud[2*j+1]=NULL;
            }
 
        }
    }
    double vol1=get_volume(polygone[0]->get_coord(),polygone[1]->get_coord(),polygone[nb_poly-1]->get_coord(),noeud1->get_coord());
    double vol2=get_volume(polygone[0]->get_coord(),polygone[nb_poly-1]->get_coord(),polygone[1]->get_coord(),noeud2->get_coord());
    if (vol1*vol2<0.)
    {
        crit=0.;
        coque.taille=0;
        return;
    }
    MG_NOEUD *noeuda=noeud1;
    MG_NOEUD *noeudb=noeud2;
    if (vol1<0.)
    {
        noeuda=noeud2;
        noeudb=noeud1;
    }
    int nb_solution,nb_triangle;
    
    if (nb_coq==4)      {
        nb_solution=2;
        nb_triangle=4;
    }
    if (nb_coq==5)      {
        nb_solution=5;
        nb_triangle=10;
    }
    if (nb_coq==6)      {
        nb_solution=14;
        nb_triangle=20;
    }
    if (nb_coq==7)      {
        nb_solution=42;
        nb_triangle=35;
    }
    if (nb_coq==8)      {
        nb_solution=132;
        nb_triangle=56;
    }
    if (nb_coq==9)      {
        nb_solution=429;
        nb_triangle=84;
    }
    if (nb_coq==10)     {
        nb_solution=1430;
        nb_triangle=120;
    }
    double crit_triangle[120];
    for (int i=0;i<nb_triangle;i++)
    {
        double crit1=OPERATEUR::qualite_tetra(polygone[tab_face[nb_coq-4][i][0]]->get_coord(),polygone[tab_face[nb_coq-4][i][1]]->get_coord(),polygone[tab_face[nb_coq-4][i][2]]->get_coord(),noeuda->get_coord());
        double crit2=OPERATEUR::qualite_tetra(polygone[tab_face[nb_coq-4][i][0]]->get_coord(),polygone[tab_face[nb_coq-4][i][2]]->get_coord(),polygone[tab_face[nb_coq-4][i][1]]->get_coord(),noeudb->get_coord());
        crit_triangle[i]=std::min(crit1,crit2);
    }
    /* examen de chaque solution  */
    double crit_opt=0.;
    int numero_solution= -1;
    double crit_solution[1430];
    for (int i=0;i<nb_solution;i++)
    {
        double volume=0.;
        for (int j=0;j<nb_coq-2;j++)
        {
            MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][0]];
            MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][1]];
            MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][2]];
            MG_NOEUD* nod=noeuda;
            volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
            nod=noeudb;
            volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
        }
        double eps=0.0018*pow(volume,0.666666666);
        if (OPERATEUR::egal(volume,coque.volume,eps))
        {
            crit_solution[i]=1.;
            for (int j=0;j<nb_coq-2;j++)
                crit_solution[i]=std::min(crit_solution[i],crit_triangle[tab_solution[nb_coq-4][i][j]]);
        }
        else crit_solution[i]=0.;
        if (crit_opt<crit_solution[i])
        {
            crit_opt=crit_solution[i];
            numero_solution=i;
        }
 
    }
    if (numero_solution==(-1))
    {
        crit=0.;
        coque.taille=0;
        return;
    }
    crit=crit_opt;
    for (int j=0;j<nb_coq-2;j++)
    {
        MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][0]];
        MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][1]];
        MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][2]];
        coque.new_tetra[8*j]=noa;
        coque.new_tetra[8*j+1]=nob;
        coque.new_tetra[8*j+2]=noc;
        coque.new_tetra[8*j+3]=noeuda;
        coque.new_tetra[8*j+4]=noa;
        coque.new_tetra[8*j+5]=noc;
        coque.new_tetra[8*j+6]=nob;
        coque.new_tetra[8*j+7]=noeudb;
    }
}
 
 
void MAILLEUR3D_OPTIMISATION::ajouter_ordre_tetra(M3D_TETRA* tet,int num)
{
    double val;
    if (num==0) val=tet->get_qualite();
    if (num==1) val=1./tet->get_qualite();
    std::pair<double,M3D_TETRA*> tmp(val,tet);
    ORDRE_TETRA::iterator p=lst_tetra[num].insert(tmp);
    lst_tetraid[num][tet->get_id()]=p;
}
 
void MAILLEUR3D_OPTIMISATION::supprimer_ordre_tetra(M3D_TETRA* tet)
{
    int num=0;
    ORDRE_TETRA_PARID::iterator it=lst_tetraid[num].find(tet->get_id());
    if (it==lst_tetraid[num].end())
    {
        num=1;
        it=lst_tetraid[num].find(tet->get_id());
    }
    if (it==lst_tetraid[num].end()) return;
    lst_tetra[num].erase(it->second);
    lst_tetraid[num].erase(it);
}
 
 
void MAILLEUR3D_OPTIMISATION::change_niveau_optimisation(int num)
{
    niveau_optimisation=num;
}
 
int MAILLEUR3D_OPTIMISATION::get_niveau_optimisation(void)
{
    return niveau_optimisation;
}
 
 
 
class MG_TRIANGLE* MAILLEUR3D_OPTIMISATION::insere_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,int origine)
{
    MG_TRIANGLE* triangle=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
    if (triangle!=NULL) return NULL;
    MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
    MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeud1->get_id(),noeud3->get_id());
    MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeud2->get_id(),noeud3->get_id());
    if (segment1==NULL) segment1=cree_segment(mgvol,noeud1,noeud2,origine);
    if (segment2==NULL) segment2=cree_segment(mgvol,noeud1,noeud3,origine);
    if (segment3==NULL) segment3=cree_segment(mgvol,noeud2,noeud3,origine);
    MG_TRIANGLE* tri=cree_triangle(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
    return tri;
}
 
 
 
 
 
 
MG_TRIANGLE* MAILLEUR3D_OPTIMISATION::cree_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_SEGMENT* segment1,MG_SEGMENT* segment2,MG_SEGMENT* segment3,int origine)
{
    M3D_TRIANGLE* tri=new M3D_TRIANGLE(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
    mg_maillage->ajouter_mg_triangle(tri);
    return tri;
}
 
 
 
MG_SEGMENT* MAILLEUR3D_OPTIMISATION::cree_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,int origine)
{
    MG_SEGMENT* seg=mg_maillage->ajouter_mg_segment(mgvol,noeud1,noeud2,origine);
    seg->change_dimension_topo_null(3);
    return seg;
}