mailleur_fem.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/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mailleur_fem.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  Derniere modification par francois
//####//  jeu 25 jui 2024 18:16:33 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "gestionversion.h"
#include "mailleur_fem.h"
#include "mg_gestionnaire.h"
#include "fem_segment2.h"
#include "fem_mini_segment2.h"
#include "fem_segment3.h"
#include "fem_triangle3.h"
#include "fem_triangle6.h"
#include "fem_quadrangle4.h"
#include "fem_quadrangle8.h"
#include "fem_tetra4.h"
#include "fem_tetra10.h"
#include "fem_hexa8.h"
#include "fem_hexa20.h"
#include "fem_penta6.h"
#include "fem_penta15.h"
#include "fem_graph_noeud.h"
#include "fem_maillage_quadratique_outils.h"
#include "ot_decalage_parametre.h"
#include "mg_arete_element.h"
#include "mg_face_element.h"
#include <math.h>
 
 
 
MAILLEUR_FEM::MAILLEUR_FEM(OT_CPU* comp):MAILLEUR(false,comp)
{
 
}
 
MAILLEUR_FEM::MAILLEUR_FEM(MAILLEUR_FEM &mdd):MAILLEUR(mdd)
{
ini_param();
}
 
 
MAILLEUR_FEM::~MAILLEUR_FEM()
{
}
 
 
int MAILLEUR_FEM::maille(FEM_MAILLAGE* fem,bool courbure_discrete,int num)
{
    int degre=fem->get_degre();
    MG_MAILLAGE* maillage=fem->get_mg_maillage();
    MG_GEOMETRIE* mg_geometrie=maillage->get_mg_geometrie();
    TPL_MAP_ENTITE<MG_SEGMENT*> mini;
    if (mg_geometrie!=NULL)
      if (mg_geometrie->get_nb_mg_poutre()!=0)
            if (mg_geometrie->get_nb_mg_volume()!=0)
                recherche_connexion_multidimension(mini,maillage,mg_geometrie);
    if (degre==1) construire_lineaire(fem,mini,maillage,mg_geometrie);
    if (degre==2) construire_quadratique(fem,mini,maillage,mg_geometrie,courbure_discrete);
 
    if (num>0) optimise_numerotation(fem);
    if (num>1) copie_numerotation_opt(fem);
 
    BOITE_3D b;
    int nx,ny,nz;
    bool structure=maillage->get_param_structure(b,nx,ny,nz);
    if (structure==true)
      fem->change_param_structure(b,nx,ny,nz);
    return OK;
}
 
 
void MAILLEUR_FEM::recherche_connexion_multidimension(TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie)
{
LISTE_MG_NOEUD::iterator it;
for (MG_NOEUD* no=maillage->get_premier_noeud(it);no!=NULL;no=maillage->get_suivant_noeud(it))
    {
    if ((no->get_lien_hexa()->get_nb()>0)  || (no->get_lien_tetra()->get_nb()))
      if (no->get_lien_segment()->get_nb()>0)  
        {
        bool estjoint=false;
        for (int i=0;i<no->get_lien_segment()->get_nb();i++)
            {
            LISTE_MG_POUTRE::iterator itp;
            for (MG_POUTRE* pou=mg_geometrie->get_premier_poutre(itp);pou!=NULL;pou=mg_geometrie->get_suivant_poutre(itp))
              if (pou->contient_element(no->get_lien_segment()->get(i)))  {estjoint=true;}
            }
        if (estjoint) 
            {
              for (int i=0;i<no->get_lien_segment()->get_nb();i++)
              {
                bool bord=false;
        LISTE_MG_POUTRE::iterator itp;
                int nbele2=no->get_lien_segment()->get(i)->get_lien_triangle()->get_nb();
                for (int j=0;j<nbele2;j++)
                  if (no->get_lien_segment()->get(i)->get_lien_triangle()->get(j)->get_lien_tetra()->get_nb()>0) 
                  if (no->get_lien_segment()->get(i)->get_lien_triangle()->get(j)->get_lien_tetra()->get_nb()<2) bord=true;
                nbele2=no->get_lien_segment()->get(i)->get_lien_quadrangle()->get_nb();
                for (int j=0;j<nbele2;j++)
                  if (no->get_lien_segment()->get(i)->get_lien_quadrangle()->get(j)->get_lien_hexa()->get_nb()>0)
                  if (no->get_lien_segment()->get(i)->get_lien_quadrangle()->get(j)->get_lien_hexa()->get_nb()<2) bord=true;
                if (bord==true) mini.ajouter(no->get_lien_segment()->get(i));
              }
            //std::cout << no->get_id() << " est un joint 1D-3D : " << mini.get_nb() << " mini "  << std::endl;
            
            }
        }
        
    }
LISTE_MG_TRIANGLE::iterator itt;
for (MG_TRIANGLE* tri=maillage->get_premier_triangle(itt);tri!=NULL;tri=maillage->get_suivant_triangle(itt))
    if (tri->get_origine()==MAGIC::ORIGINE::SECTION)
     {
     mini.ajouter(tri->get_segment1());
     mini.ajouter(tri->get_segment2());
     mini.ajouter(tri->get_segment3());
    }
     
LISTE_MG_QUADRANGLE::iterator itq;
for (MG_QUADRANGLE* qua=maillage->get_premier_quadrangle(itq);qua!=NULL;qua=maillage->get_suivant_quadrangle(itq))
    if (qua->get_origine()==MAGIC::ORIGINE::SECTION)
     {
     mini.ajouter(qua->get_segment1());
     mini.ajouter(qua->get_segment2());
     mini.ajouter(qua->get_segment3());
     mini.ajouter(qua->get_segment4());
    }
 
  
}
 
 
 
 
int MAILLEUR_FEM::construire_lineaire(FEM_MAILLAGE* fem,TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie)
{
    LISTE_MG_NOEUD::iterator it_noeud;
    MG_NOEUD* mgnoeud=maillage->get_premier_noeud(it_noeud);
    std::vector<FEM_NOEUD*> lstnoeuddirect;
    unsigned int i=0;
    while (mgnoeud)
    {
        mgnoeud->change_nouveau_numero(i);
        FEM_NOEUD* femnoeud=new FEM_NOEUD(mgnoeud);
        femnoeud->change_numero(i+1);
        i++;
        fem->ajouter_fem_noeud(femnoeud);
        lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
        if (mgnoeud->get_lien_topologie()!=NULL)
          if (mgnoeud->get_lien_topologie()->get_dimension()==0)
              {
              FEM_ELEMENT0* femele=new FEM_ELEMENT0(mgnoeud->get_lien_topologie(),mgnoeud,&femnoeud);
             fem-> ajouter_fem_element0(femele);
              }
         mgnoeud=maillage->get_suivant_noeud(it_noeud);
         
    }
    int dimsansgeo=0;
    if (mg_geometrie==NULL)
    {
        if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
        else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
        else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
    }
    if (mg_geometrie!=NULL)
    if (strcmp(mg_geometrie->get_type_geometrie(),"VIRTUEL")==0)
    {
        if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
        else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
        else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
    }
    
    if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==1)))
    {
        LISTE_MG_SEGMENT::iterator it_seg;
        MG_SEGMENT * mgseg = maillage->get_premier_segment(it_seg);
        while (mgseg)
        {
            if (mgseg->get_lien_topologie()!=NULL) 
              if (mgseg->get_lien_topologie()->get_dimension()!=1) 
                    {
                    mgseg=maillage->get_suivant_segment(it_seg);
                    continue;
                    }
            if (mg_geometrie!=NULL)
              if (mgseg->get_lien_topologie()==NULL) 
                if (dimsansgeo!=1)
                    {
                    mgseg=maillage->get_suivant_segment(it_seg);
                    continue;
                    }
                FEM_NOEUD *tabnoeud[2];
                tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
                tabnoeud[1]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
                FEM_SEGMENT2* seg=new FEM_SEGMENT2(mgseg,tabnoeud);
                fem->ajouter_fem_element1(seg);
            mgseg=maillage->get_suivant_segment(it_seg);
        }
        TPL_MAP_ENTITE<MG_SEGMENT*>::ITERATEUR itm;
        for (MG_SEGMENT* mgseg=mini.get_premier(itm);mgseg!=NULL;mgseg=mini.get_suivant(itm))
          {
          FEM_NOEUD *tabnoeud[2];
                tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
                tabnoeud[1]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
                FEM_SEGMENT2* seg=new FEM_MINI_SEGMENT2(mgseg,tabnoeud);
                fem->ajouter_fem_element1(seg);  
          }
    }
    if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==2)))
    {
        LISTE_MG_TRIANGLE::iterator it_tri;
        MG_TRIANGLE* mgtri=maillage->get_premier_triangle(it_tri);
        while (mgtri)
        {
            if (mgtri->get_lien_topologie()!=NULL)      if (mgtri->get_lien_topologie()->get_dimension()!=2) 
                    {
                    mgtri=maillage->get_suivant_triangle(it_tri);
                    continue;
                    }
            if (mg_geometrie!=NULL)
              if (mgtri->get_lien_topologie()==NULL) 
                if (dimsansgeo!=2)
                    {
                    mgtri=maillage->get_suivant_triangle(it_tri);
                    continue;
                    }
            FEM_NOEUD *tabnoeud[3];
            tabnoeud[0]=lstnoeuddirect[mgtri->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgtri->get_noeud2()->get_nouveau_numero()];
            tabnoeud[2]=lstnoeuddirect[mgtri->get_noeud3()->get_nouveau_numero()];
            FEM_TRIANGLE3* tri=new FEM_TRIANGLE3(mgtri,tabnoeud);
            fem->ajouter_fem_element2(tri);
            mgtri=maillage->get_suivant_triangle(it_tri);
        }
        LISTE_MG_QUADRANGLE::iterator it_quad;
        MG_QUADRANGLE* mgquad=maillage->get_premier_quadrangle(it_quad);
        while (mgquad)
        {
            if (mgquad->get_lien_topologie()!=NULL)     if (mgquad->get_lien_topologie()->get_dimension()!=2) {
                    mgquad=maillage->get_suivant_quadrangle(it_quad);
                    continue;
                }
            if (mg_geometrie!=NULL)
              if (mgquad->get_lien_topologie()==NULL) 
                if (dimsansgeo!=2)
                    {
                    mgquad=maillage->get_suivant_quadrangle(it_quad);
                    continue;
                    }
            FEM_NOEUD *tabnoeud[4];
            tabnoeud[0]=lstnoeuddirect[mgquad->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgquad->get_noeud2()->get_nouveau_numero()];
            tabnoeud[2]=lstnoeuddirect[mgquad->get_noeud3()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[mgquad->get_noeud4()->get_nouveau_numero()];
            FEM_QUADRANGLE4* quad=new FEM_QUADRANGLE4(mgquad,tabnoeud);
            fem->ajouter_fem_element2(quad);
            mgquad=maillage->get_suivant_quadrangle(it_quad);
        }
    }
    if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==3)))
    {
        LISTE_MG_TETRA::iterator it_tetra;
        MG_TETRA* mgtetra=maillage->get_premier_tetra(it_tetra);
        while (mgtetra)
        {
            FEM_NOEUD *tabnoeud[4];
            tabnoeud[0]=lstnoeuddirect[mgtetra->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgtetra->get_noeud2()->get_nouveau_numero()];
            tabnoeud[2]=lstnoeuddirect[mgtetra->get_noeud3()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[mgtetra->get_noeud4()->get_nouveau_numero()];
            FEM_TETRA4* tet=new FEM_TETRA4(mgtetra,tabnoeud);
            fem->ajouter_fem_element3(tet);
            mgtetra=maillage->get_suivant_tetra(it_tetra);
        }
        LISTE_MG_HEXA::iterator it_hex;
        MG_HEXA* mghex=maillage->get_premier_hexa(it_hex);
        while (mghex)
        {
            FEM_NOEUD *tabnoeud[8];
            tabnoeud[0]=lstnoeuddirect[mghex->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mghex->get_noeud2()->get_nouveau_numero()];
            tabnoeud[2]=lstnoeuddirect[mghex->get_noeud3()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[mghex->get_noeud4()->get_nouveau_numero()];
            tabnoeud[4]=lstnoeuddirect[mghex->get_noeud5()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[mghex->get_noeud6()->get_nouveau_numero()];
            tabnoeud[6]=lstnoeuddirect[mghex->get_noeud7()->get_nouveau_numero()];
            tabnoeud[7]=lstnoeuddirect[mghex->get_noeud8()->get_nouveau_numero()];
            FEM_HEXA8* hex=new FEM_HEXA8(mghex,tabnoeud);
            fem->ajouter_fem_element3(hex);
            mghex=maillage->get_suivant_hexa(it_hex);
        }
        LISTE_MG_PENTA::iterator it_pen;
        MG_PENTA* mgpen=maillage->get_premier_penta(it_pen);
        while (mgpen)
        {
            FEM_NOEUD *tabnoeud[8];
            tabnoeud[0]=lstnoeuddirect[mgpen->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgpen->get_noeud2()->get_nouveau_numero()];
            tabnoeud[2]=lstnoeuddirect[mgpen->get_noeud3()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[mgpen->get_noeud4()->get_nouveau_numero()];
            tabnoeud[4]=lstnoeuddirect[mgpen->get_noeud5()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[mgpen->get_noeud6()->get_nouveau_numero()];
            FEM_PENTA6* pen=new FEM_PENTA6(mgpen,tabnoeud);
            fem->ajouter_fem_element3(pen);
            mgpen=maillage->get_suivant_penta(it_pen);
        }
    }
    return OK;
}
 
 
int MAILLEUR_FEM::recal_element_quadratique(FEM_MAILLAGE* fem)
{
  std::map<long,FEM_ELEMENT3*> map_tet_replace;
  int size_map_tet_replace=map_tet_replace.size();
  int nbpas=(int)param.get_valeur("Quadratisation_pas");
  double dismin=param.get_valeur("Quadratisation_dis");
  int nbpassemax=(int)param.get_valeur("Quadratisation_nbpassemax");
  int nbpasse=1;
  for (int passe=0;passe<nbpasse;passe++)
  {
    char mess[255];
    sprintf(mess,"        Passe %d",passe+1);
    affiche(mess);
    LISTE_FEM_ELEMENT3::iterator itele3;
    //bool recale=false;
    for (FEM_ELEMENT3* tet=fem->get_premier_element3(itele3);tet;tet=fem->get_suivant_element3(itele3))
    {
      FEM_MAILLAGE_QUADRATIQUE_OUTILS ou;
      double dis=ou.get_distorsion2(tet);
      double jmin=ou.get_jmin(tet);;
      if ((dis<dismin-1e-12) || (jmin<0.))
        {
        OT_VECTEUR_3D vec[10];
        OT_VECTEUR_3D depart[10];
        for (int i=1;i<9;i++)
            {
            if ((i==2) || (i==4)) continue;
            double* xyzori=tet->get_fem_noeud(i)->get_coord_ori();  
            double* xyz=tet->get_fem_noeud(i)->get_coord();
            vec[i].change_x(xyz[0]-xyzori[0]);
            vec[i].change_y(xyz[1]-xyzori[1]);
            vec[i].change_z(xyz[2]-xyzori[2]);
            depart[i].change_x(xyz[0]);
            depart[i].change_y(xyz[1]);
            depart[i].change_z(xyz[2]);
            }
        int pas=0;
        while ((dis<dismin-1e-12)||(jmin<0.))
            {
            for (int i=1;i<9;i++)
                {
                if ((i==2) || (i==4)) continue;
                double* xyz=tet->get_fem_noeud(i)->get_coord();
                tet->get_fem_noeud(i)->change_x(depart[i].get_x()-pas*1./1./nbpas*vec[i].get_x());
                tet->get_fem_noeud(i)->change_y(depart[i].get_y()-pas*1./1./nbpas*vec[i].get_y());
                tet->get_fem_noeud(i)->change_z(depart[i].get_z()-pas*1./1./nbpas*vec[i].get_z());
                }
            dis=ou.get_distorsion2(tet);
            jmin=ou.get_jmin(tet);
            pas++;
            }
        map_tet_replace.insert(std::pair<long,FEM_ELEMENT3*>(tet->get_id(),tet));
        if (passe==0) sprintf(mess,"           Tetra %lu replacé à %.2lf%% de sa position optimale\t jmin=%le distorsion=%le",tet->get_id(),100.-(pas-1.)*1./nbpas*100.,jmin,dis);
        else sprintf(mess,"           Tetra %lu replacé à %.2lf%% de sa position precedente\t jmin=%le  distorsion=%le",tet->get_id(),100.-(pas-1.)*1./nbpas*100.,jmin,dis);
        affiche(mess);
        }
        
    }
    if (size_map_tet_replace!=map_tet_replace.size()) nbpasse++;
    size_map_tet_replace=map_tet_replace.size();
    if (nbpasse>nbpassemax) nbpasse=nbpassemax;
  }
  if (fem->get_nb_fem_element3()!=0)
  {
  long nb_tetra_replace = map_tet_replace.size();
  char mess[255];
  sprintf(mess,"        Nb tetra replace : %li (%lf%%)",nb_tetra_replace,((double)nb_tetra_replace/fem->get_nb_fem_element3())*(double)100.0);
  affiche(mess);
  }
  return OK;
}
 
 
int MAILLEUR_FEM::construire_quadratique(FEM_MAILLAGE* fem,TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie,bool courbure_discrete)
{
    int recal=(int)param.get_valeur("Quadratisation_jmin");
    int dimsansgeo=0;
    bool virtuel =false;
    if (mg_geometrie!=NULL)
        if (strcmp(mg_geometrie->get_type_geometrie(),"VIRTUEL")==0) virtuel=true;
        if (virtuel)
        {
            if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
            else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
            else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
        }
    LISTE_MG_NOEUD::iterator it_nd;
    MG_NOEUD* mgnoeud=maillage->get_premier_noeud(it_nd);
    if (courbure_discrete) prepare_quad_discrete(fem->get_mg_maillage(),virtuel);
    std::vector<FEM_NOEUD*> lstnoeuddirect;
    unsigned int i=0;
    while (mgnoeud)
    {
        mgnoeud->change_nouveau_numero(i);
        FEM_NOEUD* femnoeud=new FEM_NOEUD(mgnoeud);
        femnoeud->change_numero(i+1);
        i++  ;
        fem->ajouter_fem_noeud(femnoeud);
        //mgnoeud=maillage->get_suivant_noeud(it_nd);
        lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
        if (mgnoeud->get_lien_topologie()!=NULL)
            if (mgnoeud->get_lien_topologie()->get_dimension()==0)
            {
                FEM_ELEMENT0* femele=new FEM_ELEMENT0(mgnoeud->get_lien_topologie(),mgnoeud,&femnoeud);
                fem->ajouter_fem_element0(femele);
            }
            mgnoeud=maillage->get_suivant_noeud(it_nd);
        
    }
    if (mg_geometrie==NULL)
    {
        if (maillage->get_nb_mg_tetra()!=0) dimsansgeo=3;
        else if (maillage->get_nb_mg_triangle()!=0) dimsansgeo=2;
        else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
    }
        int nbmgnoeud = maillage->get_nb_mg_noeud();
    
 
    LISTE_MG_SEGMENT::iterator it_seg;
    MG_SEGMENT* mgseg=maillage->get_premier_segment(it_seg);
 
    FEM_MAILLAGE_QUADRATIQUE_OUTILS ou ;
    i=0;
    while (mgseg)
    {
        mgseg->change_nouveau_numero(i++);
        MG_NOEUD* no1t=mgseg->get_noeud1();
        MG_NOEUD* no2t=mgseg->get_noeud2();
        MG_NOEUD *noeud1=mgseg->get_noeud1();
        MG_NOEUD *noeud2=mgseg->get_noeud2();
        double *xyz1=noeud1->get_coord();
        double *xyz2=noeud2->get_coord();
        double xori=0.5*(xyz1[0]+xyz2[0]);
        double yori=0.5*(xyz1[1]+xyz2[1]);
        double zori=0.5*(xyz1[2]+xyz2[2]);
        FEM_NOEUD* femnoeud;
        if (mgseg->get_lien_topologie()!=NULL)
            if (mgseg->get_lien_topologie()->get_dimension()==1)
            {
                double L3=0,L4=0,j=0.0; 
                double xyz[3]={xori,yori,zori};
                if (!virtuel)
                {
                    MG_ARETE* arete=(MG_ARETE*)(mgseg->get_lien_topologie());
                    double t1,t2,tm;
                    arete->inverser(t1,xyz1);
                    arete->inverser(t2,xyz2);
                    if (arete->get_courbe()->est_periodique())
                        if (t1>t2) t2=t2+arete->get_courbe()->get_periode();
                    double t=0.5*(t1+t2);
                    double xyz3[3];
                    arete->evaluer(t,xyz);      
                    double L1= ou.get_distance_curviligne(0,0.5,xyz1,xyz,xyz2);
                    double L2= ou.get_distance_curviligne(0.5,1,xyz1,xyz,xyz2); 
                    if (L1/L2 >1.01)                      
                    {            
                        do                           
                        {
                            j+=0.0001;               
                            tm=(-t+t1)*j+t;      
                            arete->evaluer(tm,xyz3);             
                            L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);      
                            L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
                            
                        }while (0.5-(L4/(L3+L4))>=0.0001);
                        
                        femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
                    }
                    
                    if (L2/L1 >1.01)  
                    {
                        do                   
                        {
                            j+=0.0001;
                            tm=(-t+t2)*j+t;      
                            arete->evaluer(tm,xyz3);             
                            L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);      
                            L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
                            
                        }while (0.5-(L3/(L3+L4))>=0.0001);
                        
                        femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
                    }
                }
                if(j == 0) // geo virtuel
                {
                    if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
                    femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
                }
                
            }
        if (mgseg->get_lien_topologie()!=NULL)
            if (mgseg->get_lien_topologie()->get_dimension()==2)
                {
                    double L3=0,L4=0,j=0;  
                    double xyz[3]={xori,yori,zori};
                    if (!virtuel)
                    {
                        MG_FACE* face=(MG_FACE*)(mgseg->get_lien_topologie());
                        double uv1[2],uv2[2];
                        OT_DECALAGE_PARAMETRE decalage(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
                        face->inverser(uv1,xyz1);
                        face->inverser(uv2,xyz2);
                        if (face->get_surface()->est_periodique_u()==1)
                        {
                            double eps=1e-10*face->get_surface()->get_periode_u();
                            if (uv1[0]<-eps) uv1[0]=uv1[0]+face->get_surface()->get_periode_u();
                            if (uv1[0]>face->get_surface()->get_periode_u()-eps) uv1[0]=uv1[0]-face->get_surface()->get_periode_u();
                            if (uv2[0]<-eps) uv2[0]=uv2[0]+face->get_surface()->get_periode_u();
                            if (uv2[0]>face->get_surface()->get_periode_u()-eps) uv2[0]=uv2[0]-face->get_surface()->get_periode_u();
                        }         
                        if (face->get_surface()->est_periodique_v()==1)
                        {
                            double eps=1e-10*face->get_surface()->get_periode_v();
                            if (uv1[1]<-eps) uv1[1]=uv1[1]+face->get_surface()->get_periode_v();
                            if (uv1[1]>face->get_surface()->get_periode_v()-eps) uv1[1]=uv1[1]-face->get_surface()->get_periode_v();
                            if (uv2[1]<-eps) uv2[1]=uv2[1]+face->get_surface()->get_periode_v();
                            if (uv2[1]>face->get_surface()->get_periode_v()-eps) uv2[1]=uv2[1]-face->get_surface()->get_periode_v();
                        }        
                        double du=decalage.calcul_decalage_parametre_u(uv1[0]);
                        double dv=decalage.calcul_decalage_parametre_v(uv1[1]);
                        double u1=decalage.decalage_parametre_u(uv1[0],du);
                        double v1=decalage.decalage_parametre_v(uv1[1],dv);
                        double u2=decalage.decalage_parametre_u(uv2[0],du);
                        double v2=decalage.decalage_parametre_v(uv2[1],dv);
                        double uv[2],uv3[2],xyz3[3];
                        double milieu[3];
                        milieu[0]=0.5*(xyz1[0]+xyz2[0]);
                        milieu[1]=0.5*(xyz1[1]+xyz2[1]);
                        milieu[2]=0.5*(xyz1[2]+xyz2[2]);
                        OT_VECTEUR_3D vec(xyz1,xyz2);
                        face->inverser(uv,milieu,0.25*vec.get_longueur());
                        face->evaluer(uv,xyz);              
                        double u=decalage.decalage_parametre_u(uv[0],du);
                        double v=decalage.decalage_parametre_v(uv[1],dv);
                        
                        double L1= ou.get_distance_curviligne(0,0.5,xyz1,xyz,xyz2);
                        double L2= ou.get_distance_curviligne(0.5,1,xyz1,xyz,xyz2);
                        if (L1/L2 >1.01)                  
                        {                  
                            do                       
                            {
                                j+=0.0001;
                                if ((u1<u) && (u2<u) && (u1<u2)) u1=u1+face->get_surface()->get_periode_u();
                                if ((v1<v) && (v2<v) && (v1<v2)) v1=v1+face->get_surface()->get_periode_v();
                                
                                double u3=(-u+u1)*j+u;
                                double v3=(-v+v1)*j+v;   
                                uv3[0]=decalage.decalage_parametre_u(u3,-du);
                                uv3[1]=decalage.decalage_parametre_v(v3,-dv);
                                face->evaluer(uv3,xyz3);             
                                L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);            
                                L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);           
                                
                            }while ( 0.5-(L4/(L3+L4)) >=0.0001);                  
                            
                            femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
                        }
                        
                        if (L2/L1 >1.01) 
                        {                       
                            do                       
                            {
                                j+=0.0001;
                                if ((u1<u) && (u2<u) && (u2<u1)) u2=u2+face->get_surface()->get_periode_u();
                                if ((v1<v) && (v2<v) && (v2<v1)) v2=v2+face->get_surface()->get_periode_v();    
                                
                                double u3=(-u+u2)*j+u;
                                double v3=(-v+v2)*j+v;
                                uv3[0]=decalage.decalage_parametre_u(u3,-du);
                                uv3[1]=decalage.decalage_parametre_v(v3,-dv);
                                
                                face->evaluer(uv3,xyz3);                
                                L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
                                L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
                                
                            }while (0.5-(L3/(L3+L4)) >=0.0001);            
                            
                            femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
                        }   
                    }
                    if (j == 0.0)    // geo virtuel
                    {
                        if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
                        femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
                    }
                    
                }
        if (mgseg->get_lien_topologie()!=NULL)
                {
                if (mgseg->get_lien_topologie()->get_dimension()==3)
                    {
                        double x=0.5*(xyz1[0]+xyz2[0]);
                        double y=0.5*(xyz1[1]+xyz2[1]);
                        double z=0.5*(xyz1[2]+xyz2[2]);
                        femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),x,y,z,x,y,z);
                    }
                }
                else   // pas de geo
               {
                    double xyz[3];
                    xyz[0]=xori;
                    xyz[1]=yori;
                    xyz[2]=zori;
                    if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
 
                    femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
                }
                femnoeud->change_numero(i+nbmgnoeud+1);
                fem->ajouter_fem_noeud(femnoeud);
                mgseg=maillage->get_suivant_segment(it_seg);
                lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
    }
 
    if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==1)))
    {
        
        mgseg=maillage->get_premier_segment(it_seg);
        while (mgseg)
        {
            if (mgseg->get_lien_topologie()!=NULL) 
                if (mgseg->get_lien_topologie()->get_dimension()!=1) 
                {
                    mgseg=maillage->get_suivant_segment(it_seg);
                    continue;
                }
            if (mg_geometrie!=NULL)
                if (mgseg->get_lien_topologie()==NULL) 
                     if (dimsansgeo!=1)
                        {
                            mgseg=maillage->get_suivant_segment(it_seg);
                            continue;
                        }
            FEM_NOEUD *tabnoeud[3];
            tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgseg->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[2]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
            FEM_SEGMENT3* seg=new FEM_SEGMENT3(mgseg,tabnoeud);
            fem->ajouter_fem_element1(seg);
            mgseg=maillage->get_suivant_segment(it_seg);
        }
    }
    if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==2)))
    {
        LISTE_MG_TRIANGLE::iterator it_tri;
        MG_TRIANGLE* mgtri=maillage->get_premier_triangle(it_tri);
        while (mgtri)
        {
            if (mgtri->get_lien_topologie()!=NULL)      if (mgtri->get_lien_topologie()->get_dimension()!=2) 
            {
                mgtri=maillage->get_suivant_triangle(it_tri);
                continue;
            }
            if (mg_geometrie!=NULL)
                if (mgtri->get_lien_topologie()==NULL) 
                    if (dimsansgeo!=2)
                    {
                        mgtri=maillage->get_suivant_triangle(it_tri);
                        continue;
                    }
            FEM_NOEUD *tabnoeud[6];
            tabnoeud[0]=lstnoeuddirect[mgtri->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgtri->get_segment1()->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[2]=lstnoeuddirect[mgtri->get_noeud2()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[mgtri->get_segment2()->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[4]=lstnoeuddirect[mgtri->get_noeud3()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[mgtri->get_segment3()->get_nouveau_numero()+nbmgnoeud];
            FEM_TRIANGLE6* tri=new FEM_TRIANGLE6(mgtri,tabnoeud);
            fem->ajouter_fem_element2(tri);
            mgtri=maillage->get_suivant_triangle(it_tri);
        }
        LISTE_MG_QUADRANGLE::iterator it_quad;
        MG_QUADRANGLE* mgquad=maillage->get_premier_quadrangle(it_quad);
        while (mgquad)
        {
            if (mgquad->get_lien_topologie()!=NULL)     if (mgquad->get_lien_topologie()->get_dimension()!=2) {
                mgquad=maillage->get_suivant_quadrangle(it_quad);
                continue;
            }
            if (mg_geometrie!=NULL)
                if (mgquad->get_lien_topologie()==NULL) 
                    if (dimsansgeo!=2)
                    {
                        mgquad=maillage->get_suivant_quadrangle(it_quad);
                        continue;
                    }
                    FEM_NOEUD *tabnoeud[8];
                tabnoeud[0]=lstnoeuddirect[mgquad->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[mgquad->get_segment1()->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[2]=lstnoeuddirect[mgquad->get_noeud2()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[mgquad->get_segment2()->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[4]=lstnoeuddirect[mgquad->get_noeud3()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[mgquad->get_segment3()->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[6]=lstnoeuddirect[mgquad->get_noeud4()->get_nouveau_numero()];
            tabnoeud[7]=lstnoeuddirect[mgquad->get_segment4()->get_nouveau_numero()+nbmgnoeud];
            FEM_QUADRANGLE8* quad=new FEM_QUADRANGLE8(mgquad,tabnoeud);
            fem->ajouter_fem_element2(quad);
            mgquad=maillage->get_suivant_quadrangle(it_quad);
        }
    }
    if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==3)))
    {
        
        LISTE_MG_TETRA::iterator it_tetra;
        MG_TETRA* mgtetra=maillage->get_premier_tetra(it_tetra);
        while (mgtetra)
        {
            FEM_NOEUD *tabnoeud[10];
            tabnoeud[0]=lstnoeuddirect[mgtetra->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud1()->get_id(),mgtetra->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[2]=lstnoeuddirect[mgtetra->get_noeud2()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud2()->get_id(),mgtetra->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[4]=lstnoeuddirect[mgtetra->get_noeud3()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud3()->get_id(),mgtetra->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[6]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud1()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud2()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud3()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[9]=lstnoeuddirect[mgtetra->get_noeud4()->get_nouveau_numero()];
            FEM_TETRA10* tet=new FEM_TETRA10(mgtetra,tabnoeud);
            fem->ajouter_fem_element3(tet);
            mgtetra=maillage->get_suivant_tetra(it_tetra);
        }
        if (recal==1) if(recal_element_quadratique(fem)==FAIL) return FAIL;
        
        
        
        LISTE_MG_HEXA::iterator it_hexa;
        MG_HEXA* mghex=maillage->get_premier_hexa(it_hexa);
        while (mghex)
        {
            FEM_NOEUD *tabnoeud[20];
            tabnoeud[0]=lstnoeuddirect[mghex->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud1()->get_id(),mghex->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[2]=lstnoeuddirect[mghex->get_noeud2()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud2()->get_id(),mghex->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[4]=lstnoeuddirect[mghex->get_noeud3()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud3()->get_id(),mghex->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[6]=lstnoeuddirect[mghex->get_noeud4()->get_nouveau_numero()];
            tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud4()->get_id(),mghex->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud1()->get_id(),mghex->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[9]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud2()->get_id(),mghex->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[10]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud3()->get_id(),mghex->get_noeud7()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[11]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud4()->get_id(),mghex->get_noeud8()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[12]=lstnoeuddirect[mghex->get_noeud5()->get_nouveau_numero()];
            tabnoeud[13]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud5()->get_id(),mghex->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[14]=lstnoeuddirect[mghex->get_noeud6()->get_nouveau_numero()];
            tabnoeud[15]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud6()->get_id(),mghex->get_noeud7()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[16]=lstnoeuddirect[mghex->get_noeud7()->get_nouveau_numero()];
            tabnoeud[17]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud7()->get_id(),mghex->get_noeud8()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[18]=lstnoeuddirect[mghex->get_noeud8()->get_nouveau_numero()];
            tabnoeud[19]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud8()->get_id(),mghex->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
            FEM_HEXA20* hex=new FEM_HEXA20(mghex,tabnoeud);
            fem->ajouter_fem_element3(hex);
            mghex=maillage->get_suivant_hexa(it_hexa);
        }
        LISTE_MG_PENTA::iterator it_pen;
        MG_PENTA* mgpen=maillage->get_premier_penta(it_pen);
        while (mgpen)
        {
            FEM_NOEUD *tabnoeud[15];
            tabnoeud[0]=lstnoeuddirect[mgpen->get_noeud1()->get_nouveau_numero()];
            tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud1()->get_id(),mgpen->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[2]=lstnoeuddirect[mgpen->get_noeud2()->get_nouveau_numero()];
            tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud2()->get_id(),mgpen->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[4]=lstnoeuddirect[mgpen->get_noeud3()->get_nouveau_numero()];
            tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud3()->get_id(),mgpen->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[6]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud1()->get_id(),mgpen->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud2()->get_id(),mgpen->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud3()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[9]=lstnoeuddirect[mgpen->get_noeud4()->get_nouveau_numero()];
            tabnoeud[10]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud4()->get_id(),mgpen->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[11]=lstnoeuddirect[mgpen->get_noeud5()->get_nouveau_numero()];
            tabnoeud[12]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud5()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
            tabnoeud[13]=lstnoeuddirect[mgpen->get_noeud6()->get_nouveau_numero()];
            tabnoeud[14]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud4()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
            FEM_PENTA15* pen=new FEM_PENTA15(mgpen,tabnoeud);
            fem->ajouter_fem_element3(pen);
            mgpen=maillage->get_suivant_penta(it_pen);
        }
    }
    return OK;
}
 
 
 
 
void MAILLEUR_FEM::optimise_numerotation(FEM_MAILLAGE* fem)
{
    FEM_NOEUD* noeud=fem->get_fem_noeud(0);
    FEM_GRAPH_NOEUD *graph;
    graph=new FEM_GRAPH_NOEUD(noeud,fem);
    int e=graph->get_excentricite();
    TPL_MAP_ENTITE<FEM_NOEUD*> dernierniveau=graph->get_dernier_niveau();
    int nbdernierniveau=dernierniveau.get_nb();
    for (int i=0;i<nbdernierniveau;i++)
    {
        FEM_GRAPH_NOEUD* graphtmp=new FEM_GRAPH_NOEUD(dernierniveau.get(i),fem);
        int etmp=graphtmp->get_excentricite();
        if (etmp>e)
        {
            delete graph;
            graph=graphtmp;
            dernierniveau=graph->get_dernier_niveau();
            nbdernierniveau=dernierniveau.get_nb();
            i=-1;
            e=etmp;
        }
        else delete graphtmp;
    }
    int numero=fem->get_nb_fem_noeud();
    for (int i=0;i<e+1;i++)
    {
        TPL_MAP_ENTITE<FEM_NOEUD*> niveau=graph->get_niveau(i);
        int nbnoeud=niveau.get_nb();
        for (int j=0;j<nbnoeud;j++)
        {
            FEM_NOEUD* noeud=niveau.get(j);
            noeud->change_numero_opt(numero);
            numero--;
        }
    }
 
 
    delete graph;
}
 
 
void MAILLEUR_FEM::copie_numerotation_opt(FEM_MAILLAGE* fem)
{
    int nbnoeud=fem->get_nb_fem_noeud();
    for (int i=0;i<nbnoeud;i++)
    {
        FEM_NOEUD* noeud=fem->get_fem_noeud(i);
        noeud->change_numero(noeud->get_numero_opt());
    }
 
}
 
void MAILLEUR_FEM::prepare_quad_discrete(MG_MAILLAGE* mai,bool virtuel)
{
MG_GEOMETRIE *geo=mai->get_mg_geometrie();
liste_normale.clear();
segaprojete.clear();
segnonprojete.clear();
LISTE_MG_TETRA::iterator it;
LISTE_MG_NOEUD::iterator itn;
LISTE_MG_TRIANGLE::iterator itt;
LISTE_MG_SEGMENT::iterator its;
for (MG_TETRA* tet=mai->get_premier_tetra(it);tet!=NULL;tet=mai->get_suivant_tetra(it))
    {
    tet->change_nouveau_numero(0);
    tet->get_triangle1()->change_nouveau_numero(0);
    tet->get_triangle2()->change_nouveau_numero(0);
    tet->get_triangle3()->change_nouveau_numero(0);
    tet->get_triangle4()->change_nouveau_numero(0);
    }
for (MG_TETRA* tet=mai->get_premier_tetra(it);tet!=NULL;tet=mai->get_suivant_tetra(it))
    {
    if (tet->get_lien_topologie()!=NULL)
        {
        if (tet->get_triangle1()->get_nouveau_numero()==0)
            tet->get_triangle1()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
        else if (tet->get_triangle1()->get_nouveau_numero()>2)
            {
            if (tet->get_triangle1()->get_nouveau_numero()==tet->get_lien_topologie()->get_id()) 
                tet->get_triangle1()->change_nouveau_numero(2);
            else
                tet->get_triangle1()->change_nouveau_numero(1);
            }
        
        }
    else tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
    
    if (tet->get_lien_topologie()!=NULL)
        {
        if (tet->get_triangle2()->get_nouveau_numero()==0)
            tet->get_triangle2()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
        else if (tet->get_triangle2()->get_nouveau_numero()>2)
            {
            if (tet->get_triangle2()->get_nouveau_numero()==tet->get_lien_topologie()->get_id()) 
                tet->get_triangle2()->change_nouveau_numero(2);
            else
                tet->get_triangle2()->change_nouveau_numero(1);
            }
        
        }
    else tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
    
    if (tet->get_lien_topologie()!=NULL)
        {
        if (tet->get_triangle3()->get_nouveau_numero()==0)
            tet->get_triangle3()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
        else if (tet->get_triangle3()->get_nouveau_numero()>2)
            {
            if (tet->get_triangle3()->get_nouveau_numero()==tet->get_lien_topologie()->get_id()) 
                tet->get_triangle3()->change_nouveau_numero(2);
            else
                tet->get_triangle3()->change_nouveau_numero(1);
            }
        
        }
    else tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
    
    if (tet->get_lien_topologie()!=NULL)
        {
        if (tet->get_triangle4()->get_nouveau_numero()==0)
            tet->get_triangle4()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
        else if (tet->get_triangle4()->get_nouveau_numero()>2)
            {
            if (tet->get_triangle4()->get_nouveau_numero()==tet->get_lien_topologie()->get_id()) 
                tet->get_triangle4()->change_nouveau_numero(2);
            else
                tet->get_triangle4()->change_nouveau_numero(1);
            }
        
        }
    else tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
            
    }
if (geo!=NULL)
    {
    LISTE_MG_FACE::iterator itf;
    for (MG_FACE* face=geo->get_premier_face(itf);face!=NULL;face=geo->get_suivant_face(itf))
        {
        int nbele=face->get_lien_maillage()->get_nb();
        for (int i=0;i<nbele;i++)
            {
            MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get(i);
            ele->change_nouveau_numero(1);
            }
            
        }
    }
for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
    seg->change_nouveau_numero(0);
for (MG_TRIANGLE* tri=mai->get_premier_triangle(itt);tri!=NULL;tri=mai->get_suivant_triangle(itt))
    {
    if (tri->get_nouveau_numero()>2) tri->change_nouveau_numero(1);
    if (tri->get_lien_topologie()!=NULL)
            if (virtuel)
                if ((((MG_FACE_ELEMENT*)tri->get_lien_topologie())->get_nb_contrainte()>0))
                    {
                    segnonprojete[tri->get_segment1()]=tri->get_segment1();    
                    segnonprojete[tri->get_segment2()]=tri->get_segment2();    
                    segnonprojete[tri->get_segment3()]=tri->get_segment3();    
                    noeudcourburenulle[tri->get_noeud1()]=tri->get_noeud1();
                    noeudcourburenulle[tri->get_noeud2()]=tri->get_noeud2();
                    noeudcourburenulle[tri->get_noeud3()]=tri->get_noeud3();
                    continue;
                    }
    if (tri->get_nouveau_numero()==1)
                    {
                    segaprojete[tri->get_segment1()]=tri->get_segment1();    
                    segaprojete[tri->get_segment2()]=tri->get_segment2();    
                    segaprojete[tri->get_segment3()]=tri->get_segment3();    
                    }
    }
}
 
bool MAILLEUR_FEM::get_courbure_noeud(MG_NOEUD* no,MG_NOEUD* dirno,OT_VECTEUR_3D &n,double &r)
{
if (noeudcourburenulle[no]!=NULL)
        {
        r=0.;
        n={0.,0.,0.};
        return true;
        }
OT_VECTEUR_3D PP1(no->get_coord(),dirno->get_coord());
if (liste_normale.count(no)==0)
    {
    double deno=0.;
    OT_VECTEUR_3D n0(0.,0.,0.);
    n=n0;
    int passe=0;
    std::map<MG_NOEUD*,std::pair<int,int> > etoilesegment;
    std::map<MG_NOEUD*,std::pair<MG_TRIANGLE*,MG_TRIANGLE*> > etoilesegmenttriangle;
 
    for (int i=0;i<no->get_lien_triangle()->get_nb();i++)
        {
        MG_TRIANGLE *tri=no->get_lien_triangle()->get(i);
        if (tri->get_nouveau_numero()!=1) continue;
        MG_NOEUD* not1=tri->get_noeud1();    
        MG_NOEUD* not2=tri->get_noeud2();    
        MG_NOEUD* not3=tri->get_noeud3();    
        passe++;
        MG_NOEUD *nobase1,*nobase2;
        if (not1==no) {nobase1=not2;nobase2=not3;}
        if (not2==no) {nobase1=not3;nobase2=not1;}
        if (not3==no) {nobase1=not1;nobase2=not2;}
        if (etoilesegment.count(nobase1)==0) etoilesegment[nobase1].first=1; else etoilesegment[nobase1].second=1; 
        if (etoilesegment.count(nobase2)==0) etoilesegment[nobase2].first=-1; else etoilesegment[nobase2].second=-1; 
        if (etoilesegmenttriangle.count(nobase1)==0) etoilesegmenttriangle[nobase1].first=tri; else etoilesegmenttriangle[nobase1].second=tri; 
        if (etoilesegmenttriangle.count(nobase2)==0) etoilesegmenttriangle[nobase2].first=tri; else etoilesegmenttriangle[nobase2].second=tri; 
        }
    if (passe==0) return false;
 
    std::map<MG_NOEUD*,std::pair<int,int> > :: iterator itseg=etoilesegment.begin();
    MG_NOEUD* segdepart=itseg->first;
    MG_NOEUD* segcour=itseg->first;
    MG_TRIANGLE* tricour=NULL;
    int sens=1;
    do
        {
        if (etoilesegmenttriangle[segcour].first!=tricour) tricour=etoilesegmenttriangle[segcour].first; else tricour=etoilesegmenttriangle[segcour].second;
        if (etoilesegment[segcour].first*etoilesegment[segcour].second>0) sens=-sens;
        MG_NOEUD* not1=tricour->get_noeud1();    
        MG_NOEUD* not2=tricour->get_noeud2();    
        MG_NOEUD* not3=tricour->get_noeud3();    
        OT_VECTEUR_3D vec1(not1->get_coord(),not2->get_coord());
        OT_VECTEUR_3D vec2(not1->get_coord(),not3->get_coord());
        OT_VECTEUR_3D pvec=vec1&vec2;     
        double aire=pvec.get_longueur(); 
        pvec.norme();
        pvec=sens*pvec;
        n=n+aire*pvec;
        deno=deno+aire;
        if ((tricour->get_noeud1()!=no) && (tricour->get_noeud1()!=segcour)) segcour=tricour->get_noeud1();
        else if ((tricour->get_noeud2()!=no) && (tricour->get_noeud2()!=segcour)) segcour=tricour->get_noeud2();
        else if ((tricour->get_noeud3()!=no) && (tricour->get_noeud3()!=segcour)) segcour=tricour->get_noeud3();
        }
    while (segcour!=segdepart);
    n=(1./deno)*n;
    n.norme();
    liste_normale[no]=n;
    }
else n=liste_normale[no];
r=PP1*PP1/2./(n*PP1);
return true;               
}
 
void MAILLEUR_FEM::get_quad_seg(MG_SEGMENT *seg,bool virtuel,double *xyz)
{
    if (virtuel)
        if (seg->get_lien_topologie()!=NULL)
            if (((MG_ARETE_ELEMENT*)seg->get_lien_topologie())->get_nb_contrainte()>0) return;
    if (segaprojete[seg]==NULL) 
        return;
    if (segnonprojete[seg]!=NULL) 
        return;
    MG_NOEUD* no1=seg->get_noeud1();
    MG_NOEUD* no2=seg->get_noeud2();
    OT_VECTEUR_3D n1,n2;
    double r1,r2;
    bool def1=get_courbure_noeud(no1,no2,n1,r1);
    bool def2=get_courbure_noeud(no2,no1,n2,r2);
    if (def1==false) return;
    if (def2==false) return;
    double r=0.5*(r1+r2);
    OT_VECTEUR_3D n=0.5*(n1+n2);
    n.norme();
    double dist=0.;
    if (fabs(r)<1000.*seg->get_longueur())
        {
        dist=r*r-0.25*seg->get_longueur()*seg->get_longueur();
        dist=fabs(r)-sqrt(dist);
        if (r>0) dist=-dist;
        }
    xyz[0]=xyz[0]+dist*n.get_x();
    xyz[1]=xyz[1]+dist*n.get_y();
    xyz[2]=xyz[2]+dist*n.get_z();
}