toIbrep/src/toibrep.cpp

#include "gestionversion.h"
#include "toibrep.h"
#include "mg_file.h"
#include "ot_mathematique.h"
#include "tpl_map_entite.h"
#include "tpl_relation_entite.h"
#include "tpl_octree.h"
#include "toibrep_point.h"
#include "vct.h"
#include <math.h>
#include <fstream>
#include <string>
#include "IBrep.h"
#include "ot_cpu.h"
#include "fem_tetra4.h"
#include "xfem_tetra4.h"


TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm,int nbpas,OT_CPU* compt):geo(ge),gest(g),mai(femm),NPAS(nbpas),affichageactif(0),compteur(compt)
{
}

TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm,int nbpas):geo(ge),gest(g),mai(femm),NPAS(nbpas),affichageactif(0),compteur(NULL)
{
}
TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm):geo(ge),gest(g),mai(femm),NPAS(50),affichageactif(0),compteur(NULL)
{
}
TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm,class OT_CPU* compt):geo(ge),gest(g),mai(femm),NPAS(50),affichageactif(0),compteur(compt)
{
}

TOIBREP::~TOIBREP()
{
}


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

int TOIBREP::compare_etat(int etat,int valeur)
{
if ((etat==1) && (valeur==INTERIEUR)) return 1;
if ((etat==3) && (valeur==INTERIEUR)) return 1;
if ((etat==3) && (valeur==STRICTINTERIEUR)) return 1;
if ((etat==1) && (valeur==SUR_FACE)) return 1;
return 0;
}


int TOIBREP::estdansletetra(FEM_ELEMENT3 *tet,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3,*xyz4;
if (tet->get_nb_fem_noeud()==4)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(1)->get_coord();
        xyz3=tet->get_fem_noeud(2)->get_coord();
        xyz4=tet->get_fem_noeud(3)->get_coord();
        }
if (tet->get_nb_fem_noeud()==10)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(2)->get_coord();
        xyz3=tet->get_fem_noeud(4)->get_coord();
        xyz4=tet->get_fem_noeud(9)->get_coord();
        }
OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]);
OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]);
OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]);
OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]);
OT_MATRICE_3D mat(v1,v2,v3);
OT_MATRICE_3D mat1(v4,v2,v3);
OT_MATRICE_3D mat2(v1,v4,v3);
OT_MATRICE_3D mat3(v1,v2,v4);
double det=mat.get_determinant();
double xsi=mat1.get_determinant()/det;
double eta=mat2.get_determinant()/det;
double dseta=mat3.get_determinant()/det;
int reponse1=1;
int reponse2=1;
double eps=0.0000001;
if (xsi<-eps) reponse1=0;
if (eta<-eps) reponse1=0;
if (dseta<-eps) reponse1=0;
if (xsi+eta+dseta>1.+eps) reponse1=0;
if (xsi<eps) reponse2=0;
if (eta<eps) reponse2=0;
if (dseta<eps) reponse2=0;
if (xsi+eta+dseta>1.-eps) reponse2=0;
return reponse1+2*reponse2;
}


double TOIBREP::calculdist(double *n,double x,double y,double z,FEM_NOEUD* noeud)
{
double* xyz=noeud->get_coord();
double dist=sqrt((xyz[0]-x)*(xyz[0]-x)+(xyz[1]-y)*(xyz[1]-y)+(xyz[2]-z)*(xyz[2]-z));
OT_VECTEUR_3D vec1(n[0],n[1],n[2]);
OT_VECTEUR_3D vec2(xyz[0]-x,xyz[1]-y,xyz[2]-z);
double ps=vec1*vec2;
if (ps<0.) dist=-dist;
return dist;
}

void TOIBREP::recherche_arete_tangeante(TPL_MAP_ENTITE<MG_ARETE*> &lst,TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> &lsttopo)
{

  
LISTE_MG_ARETE::iterator it;
for (MG_ARETE* are=geo->get_premier_arete(it);are!=NULL;are=geo->get_suivant_arete(it))
    {
    if (lsttopo.existe(are))
          if (ot.arete_virtuelle(are)) lst.ajouter(are);
    }
      
}


//void TOIBREP::importer(std::string nomfichier,class MagXchange* data,std::string nomfichierout)
IBrep TOIBREP::importer(std::string nomfichier,std::string nomfichieribrep,MG_GROUPE_TOPOLOGIQUE* mggt)
{
TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lst;
if (mggt!=NULL)
        {
        int nb=mggt->get_nb();
        for (int i=0;i<nb;i++)
                {
                lst.ajouter(mggt->get(i));
                mggt->get(i)->get_topologie_sousjacente(&lst);
                }
        }
else 
        {
        int nbvol=geo->get_nb_mg_volume();
        for (int i=0;i<nbvol;i++)
                {
                lst.ajouter(geo->get_mg_volume(i));
                geo->get_mg_volume(i)->get_topologie_sousjacente(&lst);
                }
        }
affiche((char*)"Recherche aretes tangeantes");
TPL_MAP_ENTITE<MG_ARETE*> lstaretetan;
recherche_arete_tangeante(lstaretetan,lst);
char mess[300];
sprintf(mess,"      Nombre d'aretes tangeantes %d",lstaretetan.get_nb());
affiche(mess);
if (compteur!=NULL) compteur->ajouter_etape("Aretes tangeantes");
affiche((char*)"Traitement de base");
int nbface=geo->get_nb_mg_face();
int nbvraiface=0;
TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*>::ITERATEUR it;
for (MG_ELEMENT_TOPOLOGIQUE* ele=lst.get_premier(it);ele!=NULL;ele=lst.get_suivant(it))
        if (ele->get_dimension()==2) nbvraiface++;
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base selection");
int nb_noeud=mai->get_nb_fem_noeud();
int nb_element3=mai->get_nb_fem_element3();
std::string nomfichier2=nomfichier+"_dis.sol";
FEM_SOLUTION* solution=new FEM_SOLUTION(mai,nbvraiface+lstaretetan.get_nb(),(char*)nomfichier2.c_str(),nb_noeud,"LS",ENTITE_NOEUD);
int i=0;
gest->ajouter_fem_solution(solution);
nomfichier2=nomfichier+"_ele.sol";
FEM_SOLUTION* solution_ele=new FEM_SOLUTION(mai,nbvraiface+lstaretetan.get_nb(),(char*)nomfichier2.c_str(),nb_element3,"A",ENTITE_ELEMENT3);
gest->ajouter_fem_solution(solution_ele);
for (MG_ELEMENT_TOPOLOGIQUE* ele=lst.get_premier(it);ele!=NULL;ele=lst.get_suivant(it))
        {
        if (ele->get_dimension()!=2) continue;
        char mess[255];
        sprintf(mess,"F %lu",ele->get_id());
        solution->change_legende(i,mess);
        solution_ele->change_legende(i,mess);
        ++i;
        }
TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR itare;
for (MG_ARETE* are=lstaretetan.get_premier(itare);are!=NULL;are=lstaretetan.get_suivant(itare))
        {
        char mess[255];
        sprintf(mess,"A %lu",are->get_id());
        solution->change_legende(i,mess);
        solution_ele->change_legende(i,mess);
        ++i;
        }


//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base creation solution");
double xmin=1e300,ymin=1e300,zmin=1e308;
double xmax=-1e300,ymax=-1e300,zmax=-1e308;
TPL_MAP_ENTITE<FEM_NOEUD*> lst_noeud;
LISTE_FEM_NOEUD::iterator it2;
i=0;
for (FEM_NOEUD* noeud=mai->get_premier_noeud(it2);noeud!=NULL;noeud=mai->get_suivant_noeud(it2))
        {
        double* xyz=noeud->get_coord();
        if (xyz[0]<xmin) xmin=xyz[0];
        if (xyz[1]<ymin) ymin=xyz[1];
        if (xyz[2]<zmin) zmin=xyz[2];
        if (xyz[0]>xmax) xmax=xyz[0];
        if (xyz[1]>ymax) ymax=xyz[1];
        if (xyz[2]>zmax) zmax=xyz[2];
        for (int j=0;j<nbvraiface+lstaretetan.get_nb();j++)
                solution->ecrire(i,j,1e300);
        lst_noeud.ajouter(noeud);
        i++;
        }
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base recherche boite");
octree.initialiser(&lst_noeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base ini octree");
LISTE_FEM_ELEMENT3::iterator it3;
i=0;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
  {
  octree.inserer(tet);
  for (int j=0;j<nbvraiface+lstaretetan.get_nb();j++)
      solution_ele->ecrire(i,j,0);
  i++;
  }
if (compteur!=NULL) compteur->ajouter_etape("Traitement de base");
for (int i=0;i<nbface;i++)
//for (int i=0;i<1;i++)
        {
        MG_FACE* face=geo->get_mg_face(i);
        //MG_FACE* face=geo->get_mg_faceid(757);
        if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"      Face %d numero %lu",i,face->get_id());
          affiche(mess);
          }
        if (!lst.existe(face)) continue;
        vector<MG_FACE*> lstface;
        lstface.push_back(face);
         /// ici topologie virtuelle


        /*TPL_LISTE_ENTITE<double> lst;
        int type=face->get_surface()->get_type_geometrique(lst);
        int idem=0;
        for (int j=0;j<i;j++)
                {
                TPL_LISTE_ENTITE<double> lst2;
                MG_FACE* face2=geo->get_mg_face(j);
                int type2=face2->get_surface()->get_type_geometrique(lst2);
                if (type==type2)
                        if (lst.get_nb()==lst2.get_nb())
                                {
                                int diff=0;
                                for (int i=0;i<lst.get_nb();i++)
                                        if (fabs(lst.get(i)-lst2.get(i))>0.000001) diff=1;
                                if (diff==0) idem=1;
                                }
                }
        if (!idem) lstface.push_back(geo->get_mg_face(i));
        for (int j=i+1;j<nbface;j++)
                {
                TPL_LISTE_ENTITE<double> lst2;
                MG_FACE* face2=geo->get_mg_face(j);
                int type2=face2->get_surface()->get_type_geometrique(lst2);
                if (type==type2)
                        if (lst.get_nb()==lst2.get_nb())
                                {
                                int diff=0;
                                for (int i=0;i<lst.get_nb();i++)
                                        if (fabs(lst.get(i)-lst2.get(i))>0.000001) diff=1;
                                if (diff==0) lstface.push_back(face2);
                                }
                }*/
        levelsetn(&lst,&lstface,solution,solution_ele,i);
        }

TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR itare2;
i=0;
for (MG_ARETE* are=lstaretetan.get_premier(itare2);are!=NULL;are=lstaretetan.get_suivant(itare2))
        {
        if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"      Arete tangeante %d numero %lu",i,are->get_id());
          affiche(mess);
          }
        traite_arete_tangeante(are,nbface+i,solution,solution_ele);
        i++;
        }


affiche((char*)"Exportation IBREP");
IBrep brep=exporter_IBrep(nomfichieribrep,solution,solution_ele,mggt);
if (compteur!=NULL) compteur->ajouter_etape("Exportation IBREP");
affiche((char*)"Decoupage maillage");
for (int i=0;i<lstaretetan.get_nb();i++)
    decoupe_maillage(nbface+i,solution,solution_ele,0);
for (int i=0;i<nbface;i++)
  decoupe_maillage(i,solution,solution_ele,1);
if (compteur!=NULL) compteur->ajouter_etape("Decoupage maillage");
affiche((char*)"Tag maillage");
cree_tag_maillage();
if (compteur!=NULL) compteur->ajouter_etape("Tag maillage");
return brep;
}


void TOIBREP::cree_tag_maillage(void)
{
/*vector<FEM_ELEMENT3*> couche;
LISTE_FEM_ELEMENT3::iterator ittetra;
for (FEM_ELEMENT3 *tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
    {
    if (tet->est_decoupe_arete_tangeante())
        {
        for (int i=tet->get_nb_xfem()-1;i>-1;i--)
          tet->supprimer_xfem(mai,tet->get_xfem(i));
        }
    if (tet->get_nb_xfem()>0) 
        {
          tet->change_etat(2);
          couche.push_back(tet);
        }
    else tet->change_etat(3);
    }
for (int i=0;i<couche.size();i++)
    {
    FEM_ELEMENT3* tet=couche[i];
    for (int j=0;j<tet->get_nb_fem_noeud();j++)
      {
      int propage=0;
      if (tet->get_etat()==2)
          {
          for (int k=0;k<tet->get_nb_xfem();k++)
              {
              XFEM_TETRA4* xtet=(XFEM_TETRA4*)tet->get_xfem(k);
              if (xtet->est_noeud_fem(tet->get_fem_noeud(j))==1) propage=1;
              }
          }
      if (tet->get_etat()!=2) propage=1;
      if (propage==1)
        {
        for (int k=0;k<tet->get_fem_noeud(j)->get_lien_element3()->get_nb();k++)
            {
            if (tet->get_fem_noeud(j)->get_lien_element3()->get(k)->get_etat()==3)
                {
                tet->get_fem_noeud(j)->get_lien_element3()->get(k)->change_etat(1);
                couche.push_back(tet->get_fem_noeud(j)->get_lien_element3()->get(k));
                }
            }
        }
                
      }
    }
double volume=0.;
for (FEM_ELEMENT3 *tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
  {
  if (tet->get_etat()==1) volume=volume+get_volume(tet);
  if (tet->get_etat()==2) 
    {
    for (int i=0;i<tet->get_nb_xfem();i++)
      {
      XFEM_TETRA4* xtet=(XFEM_TETRA4*)tet->get_xfem(i);
      volume=volume+get_volume(xtet);
      }
    }
  
  }
char chaine[500];
sprintf(chaine,"              Volume=%le",volume);
affiche(chaine);*/
}


double TOIBREP::get_volume(FEM_ELEMENT3* tet)
{
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);
OT_VECTEUR_3D vec12(no1->get_coord(),no2->get_coord());
OT_VECTEUR_3D vec13(no1->get_coord(),no3->get_coord());
OT_VECTEUR_3D vec14(no1->get_coord(),no4->get_coord());
double vol= (vec12&vec13)*vec14*1./6.;
return vol;
}


void TOIBREP::decoupe_maillage(int num,FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele,int avecsuppression)
{
  
solution->active_solution(num);
solution_ele->active_solution(num);
LISTE_FEM_ELEMENT3::iterator ittetra;
for (FEM_ELEMENT3 *tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
      {
      if ((tet->get_solution()>0.5) || (tet->get_nb_xfem()!=0))
        {
        tet->change_etat(1);
        int nb=tet->get_nb_xfem();
        if (nb==0) 
          {
          int nb2=decoupe_tetra(tet,tet,mai);
          if (nb2!=0)
            if (avecsuppression==1) tet->decoupe_frontiere(); else tet->decoupe_arete_tangeante();
          }
        else
          for (int i=nb-1;i>-1;i--)
              {
              calcul_valeur_sous_element(tet,tet->get_xfem(i));
              int nb2=decoupe_tetra(tet,tet->get_xfem(i),mai);
              if (nb2!=0) tet->supprimer_xfem(mai,tet->get_xfem(i));
              if (nb2!=0)
                if (avecsuppression==1) tet->decoupe_frontiere(); else tet->decoupe_arete_tangeante();
              }
        if (avecsuppression==1)
          {
          nb=tet->get_nb_xfem();
          for (int i=nb-1;i>-1;i--)
              {
              if (signe_tetra(tet->get_xfem(i))==-1) 
                {
                FEM_NOEUD* no[4];
                no[0]=tet->get_xfem(i)->get_fem_noeud(0);
                no[1]=tet->get_xfem(i)->get_fem_noeud(1);
                no[2]=tet->get_xfem(i)->get_fem_noeud(2);
                no[3]=tet->get_xfem(i)->get_fem_noeud(3);
                tet->supprimer_xfem(mai,tet->get_xfem(i));
                for (int j=0;j<4;j++)
                  if (no[j]->get_numero_opt()==-1)
                    if (no[j]->get_numero()==0) 
                        mai->supprimer_fem_noeud(no[j]->get_id());
                }
              }
          }
        }
      }

}

void TOIBREP::calcul_valeur_sous_element(FEM_ELEMENT3* ele,FEM_ELEMENT3 *xele)
{
double x1=ele->get_fem_noeud(0)->get_x();
double y1=ele->get_fem_noeud(0)->get_y();
double z1=ele->get_fem_noeud(0)->get_z();
double x2=ele->get_fem_noeud(1)->get_x();
double y2=ele->get_fem_noeud(1)->get_y();
double z2=ele->get_fem_noeud(1)->get_z();
double x3=ele->get_fem_noeud(2)->get_x();
double y3=ele->get_fem_noeud(2)->get_y();
double z3=ele->get_fem_noeud(2)->get_z();
double x4=ele->get_fem_noeud(3)->get_x();
double y4=ele->get_fem_noeud(3)->get_y();
double z4=ele->get_fem_noeud(3)->get_z();
double sol1=ele->get_fem_noeud(0)->get_solution();
double sol2=ele->get_fem_noeud(1)->get_solution();
double sol3=ele->get_fem_noeud(2)->get_solution();
double sol4=ele->get_fem_noeud(3)->get_solution();
int nb=xele->get_nb_fem_noeud();
for (int i=0;i<nb;i++)
  if (xele->get_fem_noeud(i)->get_numero_opt()==1)
    {
    double x=xele->get_fem_noeud(i)->get_x();
    double y=xele->get_fem_noeud(i)->get_y();
    double z=xele->get_fem_noeud(i)->get_z();
    double u = -(x3 * y1 * z4 - x3 * y1 * z - y1 * z4 * x - z3 * y1 * x4 + z3 * y1 * x + y1 * z * x4 - y3 * x1 * z4 - z1 * y * x4 + z1 * y3 * x4 - z3 * x1 * y - z1 * y3 * x + y3 * x1 * z - z * x1 * y4 + z * x3 * y4 + z4 * y3 * x + z4 * x1 * y - z4 * x3 * y - z3 * y4 * x + z1 * y4 * x - z1 * x3 * y4 - y3 * z * x4 + z3 * x1 * y4 + z3 * y * x4 + z1 * x3 * y) / (-z1 * y2 * x3 - z1 * x2 * y4 - z1 * y3 * x4 + z1 * x2 * y3 + z1 * x3 * y4 + z3 * y1 * x4 + y2 * z3 * x1 + x2 * y1 * z4 - y1 * z2 * x4 + y1 * z2 * x3 + z1 * y2 * x4 + y3 * z2 * x4 + x1 * z2 * y4 - x1 * z2 * y3 - x3 * y1 * z4 - x3 * z2 * y4 - z3 * x1 * y4 + y2 * x3 * z4 - x2 * y3 * z4 + x2 * z3 * y4 - x2 * z3 * y1 - y2 * z3 * x4 - y2 * x1 * z4 + y3 * x1 * z4);
    double v = (z1 * y4 * x + z1 * y2 * x4 - z1 * y2 * x + z1 * x2 * y - z1 * x2 * y4 - z1 * y * x4 + z4 * x1 * y - z * x1 * y4 + y2 * x1 * z - y1 * z4 * x + y1 * z * x4 + z * x2 * y4 - z * y2 * x4 + z4 * y2 * x - y4 * z2 * x + y * z2 * x4 - z4 * x2 * y + y1 * z2 * x - x2 * y1 * z - x1 * z2 * y + x2 * y1 * z4 - y1 * z2 * x4 + x1 * z2 * y4 - y2 * x1 * z4) / (-z1 * y2 * x3 - z1 * x2 * y4 - z1 * y3 * x4 + z1 * x2 * y3 + z1 * x3 * y4 + z3 * y1 * x4 + y2 * z3 * x1 + x2 * y1 * z4 - y1 * z2 * x4 + y1 * z2 * x3 + z1 * y2 * x4 + y3 * z2 * x4 + x1 * z2 * y4 - x1 * z2 * y3 - x3 * y1 * z4 - x3 * z2 * y4 - z3 * x1 * y4 + y2 * x3 * z4 - x2 * y3 * z4 + x2 * z3 * y4 - x2 * z3 * y1 - y2 * z3 * x4 - y2 * x1 * z4 + y3 * x1 * z4);
    double w = -(-y2 * x3 * z - y2 * z3 * x1 + y2 * z3 * x + y2 * x1 * z - z3 * y1 * x + y1 * z2 * x - y1 * z2 * x3 - y3 * x1 * z - y3 * z2 * x + x2 * y3 * z + x2 * z3 * y1 - x2 * z3 * y - x2 * y1 * z - x1 * z2 * y + x1 * z2 * y3 + z1 * x2 * y + x3 * z2 * y + z3 * x1 * y + z1 * y2 * x3 + z1 * y3 * x - z1 * y2 * x - z1 * x2 * y3 + x3 * y1 * z - z1 * x3 * y) / (-z1 * y2 * x3 - z1 * x2 * y4 - z1 * y3 * x4 + z1 * x2 * y3 + z1 * x3 * y4 + z3 * y1 * x4 + y2 * z3 * x1 + x2 * y1 * z4 - y1 * z2 * x4 + y1 * z2 * x3 + z1 * y2 * x4 + y3 * z2 * x4 + x1 * z2 * y4 - x1 * z2 * y3 - x3 * y1 * z4 - x3 * z2 * y4 - z3 * x1 * y4 + y2 * x3 * z4 - x2 * y3 * z4 + x2 * z3 * y4 - x2 * z3 * y1 - y2 * z3 * x4 - y2 * x1 * z4 + y3 * x1 * z4);
    double sol=(1.-u-v-w)*sol1+u*sol2+v*sol3+w*sol4;
    xele->get_fem_noeud(i)->change_solution(sol);
    }
}


int TOIBREP::signe_tetra(FEM_ELEMENT3* tet)
{
int signep=0;
int signen=0;
int signem=0;
double val1=tet->get_fem_noeud(0)->get_solution();
double val2=tet->get_fem_noeud(1)->get_solution();
double val3=tet->get_fem_noeud(2)->get_solution();
double val4=tet->get_fem_noeud(3)->get_solution();
if (val1<1e200)
  {
  if (val1>0.00000001) signep++; 
  else if (val1<-0.00000001) signem++;
  else signen++;
  }
if (val2<1e200)
  {
  if (val2>0.00000001) signep++; 
  else if (val2<-0.00000001) signem++;
  else signen++;
  }
if (val3<1e200)
  {
  if (val3>0.00000001) signep++; 
  else if (val3<-0.00000001) signem++;
  else signen++;
  }
if (val4<1e200)
  {
  if (val4>0.00000001) signep++; 
  else if (val4<-0.00000001) signem++;
  else signen++;
  }v
if ((signep==0) && (signem>0)) return -1;
if ((signem==0) && (signep>0))  return 1;
return 0;
}


void TOIBREP::traite_arete_tangeante(MG_ARETE* are,int numsol,FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele)
{
//NPAS=150;
solution->active_solution(numsol);
solution_ele->active_solution(numsol);
vector<TOIBREP_POINT*> lst;
LISTE_FEM_ELEMENT3::iterator ittetra;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
        tet->change_etat(0);
LISTE_FEM_NOEUD::iterator itnoeud;
for (FEM_NOEUD* no=mai->get_premier_noeud(itnoeud);no!=NULL;no=mai->get_suivant_noeud(itnoeud))
        no->change_etat(0);
double tmin=are->get_tmin();
double tmax=are->get_tmax();
double eps=1e-3*(tmax-tmin);
for (int i=0;i<NPAS+1;i++)
  {
   double t=tmin+i*1.0/NPAS*(tmax-tmin);
   double xyz[3],xyzplus[3],xyzmoins[3];
   are->evaluer(t,xyz);
   TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,NULL);
   lst.push_back(pt);
   if (t+eps<tmax) 
          {
          are->evaluer(t+eps,xyzplus);
          pt->ajoute_point_eps(xyzplus);
          }
   if (t-eps>tmin) 
          {
          are->evaluer(t-eps,xyzmoins);
          pt->ajoute_point_eps(xyzmoins);
          }
   
  }
for (int i=0;i<lst.size();i++)
    {
    TOIBREP_POINT* pt=lst[i];
    double xyz[3];
    pt->get_coord3(xyz);
    TPL_MAP_ENTITE<FEM_ELEMENT3*> liste;
    octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste);
    TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR it;
    for (FEM_ELEMENT3* tet=liste.get_premier(it);tet!=NULL;tet=liste.get_suivant(it))
          {
          if (tet->get_etat()==1) continue;
          int etat=estdansletetra(tet,xyz[0],xyz[1],xyz[2]);
          if ((compare_etat(etat,SUR_FACE)) && (pt->get_interieur()))
                {
                int nbeps=pt->get_nb_point_eps();
                for (int k=0;k<nbeps;k++)
                    {
                    double xyztmp[3];
                    pt->get_point_eps(k,xyztmp);
                    int etattmp=estdansletetra(tet,xyztmp[0],xyztmp[1],xyztmp[2]);
                    if (compare_etat(etattmp,SUR_FACE)) {etat=STRICTINTERIEUR;break;}
                    if (compare_etat(etattmp,INTERIEUR)) {etat=STRICTINTERIEUR;break;}
                    if (compare_etat(etattmp,STRICTINTERIEUR)) {etat=STRICTINTERIEUR;break;}
                    }
                }
          if (compare_etat(etat,STRICTINTERIEUR))
                {
                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                        {
                        double dist=calcul_distance_level_ortho(tet->get_fem_noeud(k),are,pt);
                        tet->get_fem_noeud(k)->change_solution(dist);
                        tet->get_fem_noeud(k)->change_numero(pt->get_id());
                        tet->change_etat(1);  //temp
                        tet->get_fem_noeud(k)->change_etat(1); //temp
                        }
                }
           }
    }
remplir_trou_tangeant(&lst,are,solution_ele);    
    
    
int i=0;
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
        {
        if (noeud->get_etat()==1) solution->ecrire(i,numsol,noeud->get_solution());
        else solution->ecrire(i,numsol,1e300); 
        ++i;
        }
i=0;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
        {
        if (tet->get_etat()==1) solution_ele->ecrire(i,numsol,1.);
        ++i;
        }
int nbpt=lst.size();
for (int i=0;i<nbpt;i++) delete lst[i];
TOIBREP_POINT::remisecompteurid();
}


void TOIBREP::echantillonne_sommets(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face)
{
TPL_MAP_ENTITE<MG_SOMMET*> lstsom;
int nbbou=face->get_nb_mg_boucle();
for (int i=0;i<nbbou;i++)
    {
    MG_BOUCLE* boucle=face->get_mg_boucle(i);
    int nb_arete=boucle->get_nb_mg_coarete();
    for (int j=0;j<nb_arete;j++)
      {
      MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete();
      lstsom.ajouter(are->get_cosommet1()->get_sommet());
      lstsom.ajouter(are->get_cosommet2()->get_sommet());
      }
    }  
TPL_MAP_ENTITE<MG_SOMMET*>::ITERATEUR it;      
for (MG_SOMMET* som=lstsom.get_premier(it);som!=NULL;som=lstsom.get_suivant(it))
    {
    double xyz[3];
    som->get_point()->evaluer(xyz);
    double uv[2];
    face->inverser(uv,xyz);
    int N=NPAS;
    double rayon=1e-2;
    for (int k=0;k<N+1;k++)
          {
          double uv2[2];
          uv2[0]=uv[0]+rayon*cos(k*1.0/N*2.*M_PI);
          uv2[1]=uv[1]+rayon*sin(k*1.0/N*2.*M_PI);
          if ((face->valide_parametre_u(uv2[0])) && (face->valide_parametre_v(uv2[1]))) 
                        {
                        double distance=ot.calcule_distance_contour_face_uv(uv2,face);
                        if (distance>0.)
                          {
                          double xyz[3];
                          face->evaluer(uv2,xyz);
                          TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],1,face);
                          lst.push_back(pt);
                          }
                        }
          }
      
    }
}
void TOIBREP::echantillonne_aretes(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face)
{
echantillonne_sommets(lst,face);
int nbbou=face->get_nb_mg_boucle();
for (int i=0;i<nbbou;i++)
    {
    MG_BOUCLE* boucle=face->get_mg_boucle(i);
    int nb_arete=boucle->get_nb_mg_coarete();
    for (int j=0;j<nb_arete;j++)
      {
      MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete();
      double tmin=are->get_tmin();
      double tmax=are->get_tmax();
      double N=NPAS;
      for (int k=0;k<N+1;k++)
          {
          double t=tmin+k*1.0/N*(tmax-tmin);
          double xyz[3];
          are->evaluer(t,xyz);
          TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
          lst.push_back(pt);
//        octree.inserer(pt);
          }
     /* double t=tmin+0.2/NPAS*(tmax-tmin);
      double xyz[3];
      are->evaluer(t,xyz);
      TOIBREP_POINT *pt1=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
      lst.push_back(pt1);
      t=tmax-0.2/NPAS*(tmax-tmin);
      are->evaluer(t,xyz);
      TOIBREP_POINT *pt2=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
      lst.push_back(pt2);*/
      
      
      }
    }
}

void TOIBREP::levelsetn(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lsttopo,vector<MG_FACE*> *lstface,class FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele,int numsol)
{
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"          Echantillonnage");
          affiche(mess);
          }
        
solution->active_solution(numsol);
solution_ele->active_solution(numsol);
int nbface=lstface->size();
if (nbface==0) return;
vector<TOIBREP_POINT*> lst;
BOITE_2D boite=ot.get_boite_2D((*lstface)[0]);
echantillonne_aretes(lst,(*lstface)[0]);
for (int iface=1;iface<nbface;iface++)
        {
        MG_FACE* face=(*lstface)[iface];
        BOITE_2D boiteface=ot.get_boite_2D(face);
        boite=boite+boiteface;
        echantillonne_aretes(lst,(*lstface)[iface]);
        }
MG_FACE* face=(*lstface)[0];
int trouve=0;
int orientation;
int j=0;
do 
        {
        MG_COFACE* coface=face->get_mg_coface(j);
        MG_VOLUME* vol=coface->get_coquille()->get_mg_volume();
        if (lsttopo->existe(vol)) 
                {
                orientation=coface->get_orientation();
                trouve=1;
                }
        j++;
        }
while (trouve==0);
LISTE_FEM_ELEMENT3::iterator ittetra;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
        tet->change_etat(0);
LISTE_FEM_NOEUD::iterator itnoeud;
for (FEM_NOEUD* no=mai->get_premier_noeud(itnoeud);no!=NULL;no=mai->get_suivant_noeud(itnoeud))
        no->change_etat(0);
double umin=boite.get_xmin();
double umax=boite.get_xmax();
double vmin=boite.get_ymin();
double vmax=boite.get_ymax();
if (face->get_surface()->est_periodique_u())
        {
        umin=0;
        umax=face->get_surface()->get_periode_u();
        }
if (face->get_surface()->est_periodique_v())
        {
        vmin=0;
        vmax=face->get_surface()->get_periode_v();
        }
double epsu=1e-3*(umax-umin);
double epsv=1e-3*(vmax-vmin);
for (int i=0;i<NPAS;i++)
        for (int j=0;j<NPAS;j++)
                {
                double uv[2];
                uv[0]=umin+(i*1.0+0.5)/NPAS*(umax-umin);
                uv[1]=vmin+(j*1.0+0.5)/NPAS*(vmax-vmin);
                double xyz[3];
                if ((face->valide_parametre_u(uv[0])) && (face->valide_parametre_v(uv[1]))) 
                        {
                        double distance=ot.calcule_distance_contour_face_uv(uv,face);
                        int inte=0;
                        if (distance>0.) inte=1.; else continue;
                        face->evaluer(uv,xyz);
                        TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],inte,face);
                        lst.push_back(pt);
                        if (uv[0]+epsu<umax)
                            {
                            double uvtmp[2]={uv[0]+epsu,uv[1]};
                            double distance=ot.calcule_distance_contour_face_uv(uvtmp,face);
                            if (distance>1e-8)
                                {
                                face->evaluer(uvtmp,xyz);
                                pt->ajoute_point_eps(xyz);
                                }
                            }  
                        if (uv[0]-epsu>umin)
                            {
                            double uvtmp[2]={uv[0]-epsu,uv[1]};
                            double distance=ot.calcule_distance_contour_face_uv(uvtmp,face);
                            if (distance>1e-8)
                                {
                                face->evaluer(uvtmp,xyz);
                                pt->ajoute_point_eps(xyz);
                                }
                            }
                        if (uv[1]+epsv<vmax)
                            {
                            double uvtmp[2]={uv[0],uv[1]+epsv};
                            double distance=ot.calcule_distance_contour_face_uv(uvtmp,face);
                            if (distance>1e-8)
                                {
                                face->evaluer(uvtmp,xyz);
                                pt->ajoute_point_eps(xyz);
                                }
                            }  
                        if (uv[1]-epsv>vmin)
                            {
                            double uvtmp[2]={uv[0],uv[1]-epsv};
                            double distance=ot.calcule_distance_contour_face_uv(uvtmp,face);
                            if (distance>1e-8)
                                {
                                face->evaluer(uvtmp,xyz);
                                pt->ajoute_point_eps(xyz);
                                }
                            }
//                      octree.inserer(pt);
                        }
                }
/* debut temp
{
MG_GESTIONNAIRE gestemp;
MG_MAILLAGE *mai=new MG_MAILLAGE(NULL);
gestemp.ajouter_mg_maillage(mai);
for (int i=0;i<lst.size();i++)
{
double xyz[3];
lst[i]->get_coord3(xyz);
MG_NOEUD* no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],TRIANGULATION);
mai->ajouter_mg_noeud(no);
if (i!=0) {MG_SEGMENT* seg=new MG_SEGMENT(NULL,no,mai->get_mg_noeud(mai->get_nb_mg_noeud()-2),TRIANGULATION);mai->ajouter_mg_segment(seg);}
}

gestemp.enregistrer("debug.magic");
}*/
// fin temp
if (compteur!=NULL) compteur->ajouter_etape("Echantillonnage");
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"          Preparation newton");
          affiche(mess);
          }
calcullevelsetpremierepasse(face,orientation,&lst,0,lst.size());
if (compteur!=NULL)  compteur->ajouter_etape("Preparation newton");

if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"          Newton");
          affiche(mess);
          }
calcullevelsetdeuxiemepasse(&lst);
if (compteur!=NULL) compteur->ajouter_etape("Newton");
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"          Remplissage des trous");
          affiche(mess);
          }
remplir_trou(&lst,face,orientation);
if (compteur!=NULL) compteur->ajouter_etape("Remplissage des trous");
//etendrelevelset(solution,numsol);
int i=0;
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
        {
        if (noeud->get_etat()==1) solution->ecrire(i,numsol,noeud->get_solution());
        else solution->ecrire(i,numsol,1e300); 
        ++i;
        }
i=0;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
        {
        if (tet->get_etat()==1) solution_ele->ecrire(i,numsol,1.);
        ++i;
        }
int nbpt=lst.size();
for (int i=0;i<nbpt;i++) delete lst[i];
TOIBREP_POINT::remisecompteurid();


}
void TOIBREP::remplir_trou_tangeant(std::vector<TOIBREP_POINT*> *lst,MG_ARETE* are,FEM_SOLUTION *solution_ele)
{
MG_FACE* face1=are->get_mg_coarete(0)->get_boucle()->get_mg_face();
MG_FACE* face2=are->get_mg_coarete(1)->get_boucle()->get_mg_face();
TPL_MAP_ENTITE<FEM_ELEMENT3*> lst_element3_concerne;
int num1,num2;
int nbchamp=solution_ele->get_nb_champ();
for (int i=0;i<nbchamp;i++)
        {
        std::string nom=solution_ele->get_legende(i);
        char nom2[2];
        unsigned long id;
        sscanf(nom.c_str(),"%s %lu",nom2,&id);
        if (id==face1->get_id()) num1=i;
        if (id==face2->get_id()) num2=i;
        }  
LISTE_FEM_ELEMENT3::iterator ittet;
int i=0;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittet);tet!=NULL;tet=mai->get_suivant_element3(ittet))
{
if (solution_ele->lire(i,num1)>0.000001) 
  lst_element3_concerne.ajouter(tet);
if (solution_ele->lire(i,num2)>0.000001) 
  lst_element3_concerne.ajouter(tet);
i++;
}


int couche_active;
int nb_couche=0;
do
{
nb_couche++;
TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR  ittetra;
std::vector<FEM_ELEMENT3*> couche;
couche_active=0;
for (FEM_ELEMENT3* tet=lst_element3_concerne.get_premier(ittetra);tet!=NULL;tet=lst_element3_concerne.get_suivant(ittetra))
      {
      if (tet->get_etat()!=0) continue;
      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);
      int nbactif=0;
      int nbpositif=0;
      if ((no1->get_etat()==1) || (no1->get_etat()==-1))
            {
            nbactif++;
            if (no1->get_solution()>0.) nbpositif++;
            }
      if ((no2->get_etat()==1) || (no2->get_etat()==-1))
            {
            nbactif++;
            if (no2->get_solution()>0.) nbpositif++;
            }
      if ((no3->get_etat()==1) || (no3->get_etat()==-1))
            {
            nbactif++;
            if (no3->get_solution()>0.) nbpositif++;
            }
      if ((no4->get_etat()==1) || (no4->get_etat()==-1))
            {
            nbactif++;
            if (no4->get_solution()>0.) nbpositif++;
            }
      if ((nbactif==3) && (nbpositif<3) && (nbpositif>0))
            {
            couche.push_back(tet);
            }
      }
for (int i=0;i<couche.size();i++)
      {
      FEM_ELEMENT3* tet=couche[i];
      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);
      TOIBREP_POINT* pt;
      FEM_NOEUD* noeudref;
      if (no1->get_etat()==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];}
      if (no2->get_etat()==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];}
      if (no3->get_etat()==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];}
      if (no4->get_etat()==0) {noeudref=no4;pt=(*lst)[no1->get_numero()];}
      MG_FACE* face=pt->get_mg_face();
      double t,xyz[3];
      double dist=calcul_distance_level_ortho(noeudref,are,pt,t,xyz);
      noeudref->change_solution(dist);
      noeudref->change_etat(2);
      TOIBREP_POINT *nvpt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,0.,1,NULL);
      lst->push_back(nvpt);
      noeudref->change_numero(nvpt->get_id()); 
      }
for (FEM_ELEMENT3* tet=lst_element3_concerne.get_premier(ittetra);tet!=NULL;tet=lst_element3_concerne.get_suivant(ittetra))
      {
      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);
      if (no1->get_etat()!=0)
      if (no2->get_etat()!=0)
      if (no3->get_etat()!=0)
      if (no4->get_etat()!=0)
      if (tet->get_etat()==0)
          {
          int reactive=0;
          for (int i=0;i<tet->get_nb_fem_noeud();i++)
              {
              if (reactive) continue;
              FEM_NOEUD* no=tet->get_fem_noeud(i);
              TOIBREP_POINT *pt=(*lst)[no->get_numero()];
              double t;
              pt->get_coord1(t);
              double eps=1e-4*(are->get_tmax()-are->get_tmin());
              if ((t>=are->get_tmin()-eps) && (t<are->get_tmax()+eps)) 
                {
                   tet->change_etat(1);
                   no1->change_etat(1);
                   no2->change_etat(1);
                   no3->change_etat(1);
                   no4->change_etat(1);
                   reactive=1;
                   couche_active=1;
                   
                 }
                
              }
           if (!reactive)
              {
              if (no1->get_etat()==2) no1->change_etat(-1);
              if (no2->get_etat()==2) no2->change_etat(-1);
              if (no3->get_etat()==2) no3->change_etat(-1);
              if (no4->get_etat()==2) no4->change_etat(-1);
              tet->change_etat(-1);
              }
          }
     }
}
while (couche_active==1);       
TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR  ittetra;
for (FEM_ELEMENT3* tet=lst_element3_concerne.get_premier(ittetra);tet!=NULL;tet=lst_element3_concerne.get_suivant(ittetra))
      {
      if (tet->get_etat()!=1) continue;
      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);
      int nbpositif=0;
      int nbnegatif=0;
      double epsdist=1e-6*tet->get_boite_3D().get_rayon();
      if (no1->get_solution()>epsdist) nbpositif++;
      if (no1->get_solution()<-epsdist) nbnegatif++;
      if (no2->get_solution()>epsdist) nbpositif++;
      if (no2->get_solution()<-epsdist) nbnegatif++;
      if (no3->get_solution()>epsdist) nbpositif++;
      if (no3->get_solution()<-epsdist) nbnegatif++;
      if (no4->get_solution()>epsdist) nbpositif++;
      if (no4->get_solution()<-epsdist) nbnegatif++;
      if (!(((nbnegatif>0) && (nbpositif>0)) ||(nbnegatif+nbpositif==1)))
         tet->change_etat(0);
      }
}


void TOIBREP::remplir_trou(std::vector<TOIBREP_POINT*> *lst,MG_FACE* face,int sens)
{
int castot=0;
int cas=0;
int couche_active;
int nb_couche=0;
LISTE_FEM_ELEMENT3::iterator ittetra;
do
{
nb_couche++;

std::vector<FEM_ELEMENT3*> couche;
couche_active=0;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
      {
      if (tet->get_etat()!=0) continue;
      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);
      int nbactif=0;
      int nbpositif=0;
      if ((no1->get_etat()==1) || (no1->get_etat()==-1))
            {
            nbactif++;
            if (no1->get_solution()>0.) nbpositif++;
            }
      if ((no2->get_etat()==1) || (no2->get_etat()==-1))
            {
            nbactif++;
            if (no2->get_solution()>0.) nbpositif++;
            }
      if ((no3->get_etat()==1) || (no3->get_etat()==-1))
            {
            nbactif++;
            if (no3->get_solution()>0.) nbpositif++;
            }
      if ((no4->get_etat()==1) || (no4->get_etat()==-1))
            {
            nbactif++;
            if (no4->get_solution()>0.) nbpositif++;
            }
      if ((nbactif==3) && (nbpositif<3) && (nbpositif>0))
            {
            couche.push_back(tet);
            }
      }
for (int i=0;i<couche.size();i++)
      {
        FEM_ELEMENT3* tet=couche[i];
        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);
      TOIBREP_POINT* pt;
         FEM_NOEUD* noeudref;
            if (no1->get_etat()==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];}
            if (no2->get_etat()==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];}
            if (no3->get_etat()==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];}
            if (no4->get_etat()==0) {noeudref=no4;pt=(*lst)[no1->get_numero()];}
            MG_FACE* face=pt->get_mg_face();
            double uv[2];
            double xyz[3];
            double dist=calcul_distance(noeudref,face,pt,xyz,uv);
            double normal[3];
            face->calcul_normale_unitaire(uv,normal);
            normal[0]=normal[0]*sens*(-1.);
            normal[1]=normal[1]*sens*(-1.);
            normal[2]=normal[2]*sens*(-1.);
            double dist2=calculdist(normal,xyz[0],xyz[1],xyz[2],noeudref);
            if (dist2<0.) dist=dist*(-1);
            noeudref->change_solution(dist);
            noeudref->change_etat(2);
            TOIBREP_POINT *nvpt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],1,face);
            lst->push_back(nvpt);
            noeudref->change_numero(nvpt->get_id()); 
            
      }
 for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
      {
      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);
      if (no1->get_etat()!=0)
      if (no2->get_etat()!=0)
      if (no3->get_etat()!=0)
      if (no4->get_etat()!=0)
      if (tet->get_etat()==0)
          {
          castot++;
          double lstptdecoup[5*3];
          int nb=0;
          decoupe_tetra_noeud(tet,lstptdecoup,&nb);
          int reactive=0;
          for (int i=0;i<nb;i++)
              {
              if (reactive) continue;
              double uv[2];
              face->inverser(uv,lstptdecoup+3*i);
              if ((face->valide_parametre_u(uv[0])) && (face->valide_parametre_v(uv[1]))) 
                {
                double distance=ot.calcule_distance_contour_face_uv(uv,face);
                if (distance>1e-6)
                  {
                  cas++;
                  tet->change_etat(1);
                  no1->change_etat(1);
                  no2->change_etat(1);
                  no3->change_etat(1);
                  no4->change_etat(1);
                  reactive=1;
                  couche_active=1;
                  }
                }
              }
           if (!reactive)
              {
              if (no1->get_etat()==2) no1->change_etat(-1);
              if (no2->get_etat()==2) no2->change_etat(-1);
              if (no3->get_etat()==2) no3->change_etat(-1);
              if (no4->get_etat()==2) no4->change_etat(-1);
              tet->change_etat(-1);
              }
          }
      }
}
while (couche_active==1);       
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittetra);tet!=NULL;tet=mai->get_suivant_element3(ittetra))
      {
      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);
     /* if (no1->get_etat()==1)
        if (no2->get_etat()==1)
          if (no3->get_etat()==1)
            if (no4->get_etat()==1)
              if (tet->get_etat()!=1)
                tet->change_etat(1);*/
        
      if (tet->get_etat()!=1) continue;
      int nbpositif=0;
      int nbnegatif=0;
      double epsdist=1e-6*tet->get_boite_3D().get_rayon();
      if (no1->get_solution()>epsdist) nbpositif++;
      if (no1->get_solution()<-epsdist) nbnegatif++;
      if (no2->get_solution()>epsdist) nbpositif++;
      if (no2->get_solution()<-epsdist) nbnegatif++;
      if (no3->get_solution()>epsdist) nbpositif++;
      if (no3->get_solution()<-epsdist) nbnegatif++;
      if (no4->get_solution()>epsdist) nbpositif++;
      if (no4->get_solution()<-epsdist) nbnegatif++;
      if (!(((nbnegatif>0) && (nbpositif>0)) ||(nbnegatif+nbpositif==1)))
         tet->change_etat(0);
      }
 double res=cas*100.0/castot;
char mess[255];
sprintf(mess,"                  %lf%% en %d couches",res,nb_couche);
affiche(mess);
}

int TOIBREP::decoupe_tetra(FEM_ELEMENT3* tet,FEM_ELEMENT3* xtet,FEM_MAILLAGE* mai)
{
  FEM_NOEUD* no1=xtet->get_fem_noeud(0);
  FEM_NOEUD* no2=xtet->get_fem_noeud(1);
  FEM_NOEUD* no3=xtet->get_fem_noeud(2);
  FEM_NOEUD* no4=xtet->get_fem_noeud(3);
  if (no1->get_solution()*no2->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no2->get_solution());
        double x=no1->get_x()+t*(no2->get_x()-no1->get_x());
        double y=no1->get_y()+t*(no2->get_y()-no1->get_y());
        double z=no1->get_z()+t*(no2->get_z()-no1->get_z());
        FEM_NOEUD* no=new FEM_NOEUD(tet->get_mg_element_maillage(),x,y,z);
        no->change_solution(0.0);
        no->change_numero(0);
        no->change_numero_opt(1);
        mai->ajouter_fem_noeud(no);
        FEM_NOEUD *lstno[4];
        lstno[0]=no1;
        lstno[1]=no;
        lstno[2]=no3;
        lstno[3]=no4;
        XFEM_TETRA4* tet1=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        lstno[0]=no;
        lstno[1]=no2;
        XFEM_TETRA4* tet2=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        int nb=decoupe_tetra(tet,tet1,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet1);
        nb=decoupe_tetra(tet,tet2,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet2);
        return 2;
        }
  if (no1->get_solution()*no3->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no3->get_solution());
        double x=no1->get_x()+t*(no3->get_x()-no1->get_x());
        double y=no1->get_y()+t*(no3->get_y()-no1->get_y());
        double z=no1->get_z()+t*(no3->get_z()-no1->get_z());
        FEM_NOEUD* no=new FEM_NOEUD(tet->get_mg_element_maillage(),x,y,z);
        no->change_solution(0.0);
        no->change_numero(0);
        no->change_numero_opt(1);
        mai->ajouter_fem_noeud(no);
        FEM_NOEUD *lstno[4];
        lstno[0]=no1;
        lstno[1]=no2;
        lstno[2]=no;
        lstno[3]=no4;
        XFEM_TETRA4* tet1=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        lstno[0]=no;
        lstno[2]=no3;
        XFEM_TETRA4* tet2=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        int nb=decoupe_tetra(tet,tet1,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet1);
        nb=decoupe_tetra(tet,tet2,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet2);
        return 2;
        }
   if (no1->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no4->get_solution());
        double x=no1->get_x()+t*(no4->get_x()-no1->get_x());
        double y=no1->get_y()+t*(no4->get_y()-no1->get_y());
        double z=no1->get_z()+t*(no4->get_z()-no1->get_z());
        FEM_NOEUD* no=new FEM_NOEUD(tet->get_mg_element_maillage(),x,y,z);
        no->change_solution(0.0);
        no->change_numero(0);
        no->change_numero_opt(1);
        mai->ajouter_fem_noeud(no);
        FEM_NOEUD *lstno[4];
        lstno[0]=no1;
        lstno[1]=no2;
        lstno[2]=no3;
        lstno[3]=no;
        XFEM_TETRA4* tet1=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        lstno[0]=no;
        lstno[3]=no4;
        XFEM_TETRA4* tet2=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        int nb=decoupe_tetra(tet,tet1,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet1);
        nb=decoupe_tetra(tet,tet2,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet2);
        return 2;
        }
  if (no2->get_solution()*no3->get_solution()<-0.00000001)
        {
        double t=no2->get_solution()/(no2->get_solution()-no3->get_solution());
        double x=no2->get_x()+t*(no3->get_x()-no2->get_x());
        double y=no2->get_y()+t*(no3->get_y()-no2->get_y());
        double z=no2->get_z()+t*(no3->get_z()-no2->get_z());
        FEM_NOEUD* no=new FEM_NOEUD(tet->get_mg_element_maillage(),x,y,z);
        no->change_solution(0.0);
        no->change_numero(0);
        no->change_numero_opt(1);
        mai->ajouter_fem_noeud(no);
        FEM_NOEUD *lstno[4];
        lstno[0]=no1;
        lstno[1]=no2;
        lstno[2]=no;
        lstno[3]=no4;
        XFEM_TETRA4* tet1=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        lstno[1]=no;
        lstno[2]=no3;
        XFEM_TETRA4* tet2=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        int nb=decoupe_tetra(tet,tet1,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet1);
        nb=decoupe_tetra(tet,tet2,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet2);
        return 2;
        }
    if (no2->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no2->get_solution()/(no2->get_solution()-no4->get_solution());
        double x=no2->get_x()+t*(no4->get_x()-no2->get_x());
        double y=no2->get_y()+t*(no4->get_y()-no2->get_y());
        double z=no2->get_z()+t*(no4->get_z()-no2->get_z());
        FEM_NOEUD* no=new FEM_NOEUD(tet->get_mg_element_maillage(),x,y,z);
        no->change_solution(0.0);
        no->change_numero_opt(1);
        no->change_numero(0);
        mai->ajouter_fem_noeud(no);
        FEM_NOEUD *lstno[4];
        lstno[0]=no1;
        lstno[1]=no2;
        lstno[2]=no3;
        lstno[3]=no;
        XFEM_TETRA4* tet1=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        lstno[1]=no;
        lstno[3]=no4;
        XFEM_TETRA4* tet2=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        int nb=decoupe_tetra(tet,tet1,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet1);
        nb=decoupe_tetra(tet,tet2,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet2);
        return 2;
        }
    if (no3->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no3->get_solution()/(no3->get_solution()-no4->get_solution());
        double x=no3->get_x()+t*(no4->get_x()-no3->get_x());
        double y=no3->get_y()+t*(no4->get_y()-no3->get_y());
        double z=no3->get_z()+t*(no4->get_z()-no3->get_z());
        FEM_NOEUD* no=new FEM_NOEUD(tet->get_mg_element_maillage(),x,y,z);
        no->change_solution(0.0);
        no->change_numero_opt(1);
        no->change_numero(0);
        mai->ajouter_fem_noeud(no);
        FEM_NOEUD *lstno[4];
        lstno[0]=no1;
        lstno[1]=no2;
        lstno[2]=no3;
        lstno[3]=no;
        XFEM_TETRA4* tet1=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        lstno[2]=no;
        lstno[3]=no4;
        XFEM_TETRA4* tet2=new XFEM_TETRA4(mai,(FEM_TETRA4*)tet,lstno);
        int nb=decoupe_tetra(tet,tet1,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet1);
        nb=decoupe_tetra(tet,tet2,mai);
        if (nb!=0) tet->supprimer_xfem(mai,tet2);
        return 2;
        }
return 0;
}


void TOIBREP::decoupe_tetra_noeud(FEM_ELEMENT3* tet,double* lst,int *nb)
{
  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);
  if (no1->get_solution()*no2->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no2->get_solution());
        double x=no1->get_x()+t*(no2->get_x()-no1->get_x());
        double y=no1->get_y()+t*(no2->get_y()-no1->get_y());
        double z=no1->get_z()+t*(no2->get_z()-no1->get_z());
        lst[3*(*nb)]=x;
        lst[3*(*nb)+1]=y;
        lst[3*(*nb)+2]=z;
        (*nb)++;
        }
  if (no1->get_solution()*no3->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no3->get_solution());
        double x=no1->get_x()+t*(no3->get_x()-no1->get_x());
        double y=no1->get_y()+t*(no3->get_y()-no1->get_y());
        double z=no1->get_z()+t*(no3->get_z()-no1->get_z());
        lst[3*(*nb)]=x;
        lst[3*(*nb)+1]=y;
        lst[3*(*nb)+2]=z;
        (*nb)++;
        }
  if (no1->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no4->get_solution());
        double x=no1->get_x()+t*(no4->get_x()-no1->get_x());
        double y=no1->get_y()+t*(no4->get_y()-no1->get_y());
        double z=no1->get_z()+t*(no4->get_z()-no1->get_z());
        lst[3*(*nb)]=x;
        lst[3*(*nb)+1]=y;
        lst[3*(*nb)+2]=z;
        (*nb)++;
        }
  if (no2->get_solution()*no3->get_solution()<-0.00000001)
        {
        double t=no2->get_solution()/(no2->get_solution()-no3->get_solution());
        double x=no2->get_x()+t*(no3->get_x()-no2->get_x());
        double y=no2->get_y()+t*(no3->get_y()-no2->get_y());
        double z=no2->get_z()+t*(no3->get_z()-no2->get_z());
        lst[3*(*nb)]=x;
        lst[3*(*nb)+1]=y;
        lst[3*(*nb)+2]=z;
        (*nb)++;
        }
  if (no2->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no2->get_solution()/(no2->get_solution()-no4->get_solution());
        double x=no2->get_x()+t*(no4->get_x()-no2->get_x());
        double y=no2->get_y()+t*(no4->get_y()-no2->get_y());
        double z=no2->get_z()+t*(no4->get_z()-no2->get_z());
        lst[3*(*nb)]=x;
        lst[3*(*nb)+1]=y;
        lst[3*(*nb)+2]=z;
        (*nb)++;
        }
  if (no3->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no3->get_solution()/(no3->get_solution()-no4->get_solution());
        double x=no3->get_x()+t*(no4->get_x()-no3->get_x());
        double y=no3->get_y()+t*(no4->get_y()-no3->get_y());
        double z=no3->get_z()+t*(no4->get_z()-no3->get_z());
        lst[3*(*nb)]=x;
        lst[3*(*nb)+1]=y;
        lst[3*(*nb)+2]=z;
        (*nb)++;
        }
}

void TOIBREP::calcullevelsetpremierepasse(MG_FACE* face,int sens,vector<TOIBREP_POINT*> *lst,int n1,int n2)
{
/*LISTE_FEM_ELEMENT3::iterator it3;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
{
double *xyz1,*xyz2,*xyz3,*xyz4;
if (tet->get_nb_fem_noeud()==4)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(1)->get_coord();
        xyz3=tet->get_fem_noeud(2)->get_coord();
        xyz4=tet->get_fem_noeud(3)->get_coord();
        }
if (tet->get_nb_fem_noeud()==10)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(2)->get_coord();
        xyz3=tet->get_fem_noeud(4)->get_coord();
        xyz4=tet->get_fem_noeud(9)->get_coord();
        }
double xcentre=0.25*(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0]);  
double ycentre=0.25*(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1]);  
double zcentre=0.25*(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2]);
double rayon1=(xyz1[0]-xcentre)*(xyz1[0]-xcentre)+(xyz1[1]-ycentre)*(xyz1[1]-ycentre)+(xyz1[2]-zcentre)*(xyz1[2]-zcentre);
double rayon2=(xyz2[0]-xcentre)*(xyz2[0]-xcentre)+(xyz2[1]-ycentre)*(xyz2[1]-ycentre)+(xyz2[2]-zcentre)*(xyz2[2]-zcentre);
double rayon3=(xyz3[0]-xcentre)*(xyz3[0]-xcentre)+(xyz3[1]-ycentre)*(xyz3[1]-ycentre)+(xyz3[2]-zcentre)*(xyz3[2]-zcentre);
double rayon4=(xyz4[0]-xcentre)*(xyz4[0]-xcentre)+(xyz4[1]-ycentre)*(xyz4[1]-ycentre)+(xyz4[2]-zcentre)*(xyz4[2]-zcentre);
double rayon=std::max(rayonETRA_ACTIF1,std::max(rayon2,std::max(rayon3,rayon4)));
rayon=sqrt(rayon);
TPL_MAP_ENTITE<TOIBREP_POINT*> liste;
octree.rechercher(xcentre,ycentre,zcentre,rayon,liste);
TPL_MAP_ENTITE<TOIBREP_POINT*>::ITERATEUR it;
for (TOIBREP_POINT* pt=liste.get_premier(it);pt!=NULL;pt=liste.get_suivant(it))
      {
      double xyzpt[3];
      pt->get_coord3(xyzpt);
      int etat=estdansletetra(tet,xyzpt[0],xyzpt[1],xyzpt[2]);
      if (compare_etat(etat,INTERIEUR))
                {
                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                        {
                        double normal[3];
                        double uv[2];
                        face->inverser(uv,xyzpt);
                        face->calcul_normale_unitaire(uv,normal);
                        normal[0]=normal[0]*sens*(-1.);
                        normal[1]=normal[1]*sens*(-1.);
                        normal[2]=normal[2]*sens*(-1.);
                        double dist=calculdist(normal,xyzpt[0],xyzpt[1],xyzpt[2],tet->get_fem_noeud(k));
                        if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution()))
                            {
                            tet->get_fem_noeud(k)->change_solution(dist);
                            tet->get_fem_noeud(k)->change_numero(pt->get_id());
                            pt->change_coord2(uv);
                            }
                        }
                }
      if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur()))
                {
                tet->change_etat(1);
                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                    tet->get_fem_noeud(k)->change_etat(1);
                }
      if (tet->get_etat()==1) break;
      }

}
octree.vide();
it}*/
/*LISTE_FEM_ELEMENT3::iterator it3;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
        octree.inserer(tet);*/
for (int i=n1;i<n2;i++)
                {
                TOIBREP_POINT* pt=(*lst)[i];
                double uv[2];
                double xyz[3];
                pt->get_coord3(xyz);
                double normal[3];
                face->inverser(uv,xyz);
                face->calcul_normale_unitaire(uv,normal);
                normal[0]=normal[0]*sens*(-1.);
                normal[1]=normal[1]*sens*(-1.);
                normal[2]=normal[2]*sens*(-1.);
                TPL_MAP_ENTITE<FEM_ELEMENT3*> liste;
                octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste);
                TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR it;
                for (FEM_ELEMENT3* tet=liste.get_premier(it);tet!=NULL;tet=liste.get_suivant(it))
                        {
                        if (tet->get_etat()==1) continue;
                        int etat=estdansletetra(tet,xyz[0],xyz[1],xyz[2]);
                        if ((compare_etat(etat,SUR_FACE)) && (pt->get_interieur()))
                          {
                          int nbeps=pt->get_nb_point_eps();
                          for (int k=0;k<nbeps;k++)
                              {
                              double xyztmp[3];
                              pt->get_point_eps(k,xyztmp);
                              int etattmp=estdansletetra(tet,xyztmp[0],xyztmp[1],xyztmp[2]);
                              if (compare_etat(etattmp,SUR_FACE)) {etat=STRICTINTERIEUR;break;}
                              if (compare_etat(etattmp,INTERIEUR)) {etat=STRICTINTERIEUR;break;}
                              if (compare_etat(etattmp,STRICTINTERIEUR)) {etat=STRICTINTERIEUR;break;}
                              }
                          }
                        if (compare_etat(etat,INTERIEUR))
                                {
                                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                                        {
                                        double dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_fem_noeud(k));
                                        if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution()))
                                                {
                                                tet->get_fem_noeud(k)->change_solution(dist);
                                                tet->get_fem_noeud(k)->change_numero(pt->get_id());
                                                pt->change_coord2(uv);
                                                }
                                        }
                                }
                        if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur()))
                                {
                                tet->change_etat(1);
                                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                                  {
                                  tet->get_fem_noeud(k)->change_etat(1);
                                  //octree.supprimer(tet);
                                  }
                                }
                        }
                }
//octree.vide();


}       

void TOIBREP::calcullevelsetdeuxiemepasse(vector<TOIBREP_POINT*> *lst)
{
LISTE_FEM_NOEUD::iterator itno;
for (FEM_NOEUD* no=mai->get_premier_noeud(itno);no!=NULL;no=mai->get_suivant_noeud(itno))
        if (no->get_solution()<1e250)
                {
                TOIBREP_POINT* pt=(*lst)[no->get_numero()];
                MG_FACE* face=pt->get_mg_face();
                double xyz[3],uv[2];
                double dist=calcul_distance(no,face,pt,xyz,uv);
                int signe=1;
                if (no->get_solution()<0) signe=-1;
                no->change_solution(signe*dist);
                }


}


double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_ARETE* are,TOIBREP_POINT* pt,double &tii,double *xyz,double precision)
{
double eps;
are->inverser(tii,noeud->get_coord());
int compteur=0;
OT_VECTEUR_3D Pt(noeud->get_x(),noeud->get_y(),noeud->get_z());
do
        {
        compteur++;
        double ti=tii;
        double xyz[3],dxyz[3],ddxyz[3];
        are->deriver_seconde(ti,ddxyz,dxyz,xyz);
        OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]);
        OT_VECTEUR_3D Ct_deriver(dxyz[0],dxyz[1],dxyz[2]);
        OT_VECTEUR_3D Ct_deriver_seconde(ddxyz[0],ddxyz[1],ddxyz[2]);
        OT_VECTEUR_3D Distance = Ct-Pt;
        tii=ti-Ct_deriver*Distance/(Ct_deriver_seconde*Distance+Ct_deriver.get_longueur2());
        eps=fabs(tii-ti);
        if (compteur>500) return 1e300;
        if (tii<are->get_tmin()) 
                {
                tii=are->get_tmin();
                eps=0.;
                }
        if (tii>are->get_tmax()) 
                {
                tii=are->get_tmax();
                eps=0.;
                }
        }
while (eps>precision);
are->evaluer(tii,xyz);
OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]);
double distance=(Ct-Pt).get_longueur();
return distance;
}


double TOIBREP::calcul_distance_level_ortho(FEM_NOEUD* no,MG_ARETE* are,TOIBREP_POINT* pt,double precision)
{
double t;
double xyz[3];
return calcul_distance_level_ortho(no,are,pt,t,xyz,precision);
}

double TOIBREP::calcul_distance_level_ortho(FEM_NOEUD* no,MG_ARETE* are,TOIBREP_POINT* pt,double &t,double *xyz,double precision)
{
calcul_distance(no,are,pt,t,xyz,precision);
double dir[3];
are->deriver(t,dir);
double uv[2];
MG_FACE* face=are->get_mg_coarete(0)->get_boucle()->get_mg_face();
face->inverser(uv,xyz);
double nor[3];
face->calcul_normale_unitaire(uv,nor);
OT_VECTEUR_3D dirv(dir);
OT_VECTEUR_3D norv(nor);
OT_VECTEUR_3D sensv=norv&dirv;
double d=-(sensv.get_x()*xyz[0]+sensv.get_y()*xyz[1],sensv.get_z()*xyz[2]);
double dist=fabs(sensv.get_x()*no->get_x()+sensv.get_y()*no->get_y()+sensv.get_z()*no->get_z()+d);
OT_VECTEUR_3D distv(no->get_x()-xyz[0],no->get_y()-xyz[1],no->get_z()-xyz[2]);
if (sensv*distv<0.) dist=-dist;
return dist;
}


double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_FACE* face,TOIBREP_POINT* pt,double *xyz,double *uv,double precision)
{
double uvii[2],eps;
pt->get_coord2(uvii);
int compteur=0;
OT_VECTEUR_3D Pt(noeud->get_x(),noeud->get_y(),noeud->get_z());
double delta_u,delta_v;
do
        {
        compteur++;
        double uvi[2];
        uvi[0]=uvii[0];
        uvi[1]=uvii[1];
        double xyzduu[3],xyzdvv[3],xyzduv[3],xyzdu[3],xyzdv[3],xyz[3];
        face->deriver_seconde(uvi,xyzduu,xyzduv,xyzdvv,xyz,xyzdu,xyzdv);
        OT_VECTEUR_3D S(xyz[0],xyz[1],xyz[2]);
        OT_VECTEUR_3D Su(xyzdu[0],xyzdu[1],xyzdu[2]);
        OT_VECTEUR_3D Sv(xyzdv[0],xyzdv[1],xyzdv[2]);
        OT_VECTEUR_3D Suu(xyzduu[0],xyzduu[1],xyzduu[2]);
        OT_VECTEUR_3D Suv(xyzduv[0],xyzduv[1],xyzduv[2]);
        OT_VECTEUR_3D Svv(xyzdvv[0],xyzdvv[1],xyzdvv[2]);
        OT_VECTEUR_3D Distance = S-Pt;
        double a[4],b[2];
        a[0]=Su.get_longueur2()+Distance*Suu;
        a[1]=Su*Sv+Distance*Suv;
        a[2]=Su*Sv+Distance*Suv;
        a[3]=Sv.get_longueur2()+Distance*Svv;
        b[0]=Distance*Su;b[0]=-b[0];
        b[1]=Distance*Sv;b[1]=-b[1];
        double deltau,deltav;
        double denominateur_delta=(a[0]*a[3]-a[2]*a[1]);
        if (a[0]<1E-12)
                deltau=0;
        else delta_u=(b[0]*a[3]-b[1]*a[1])/denominateur_delta;
        if (a[3]<1E-12)
                deltav=0;
        else delta_v=(a[0]*b[1]-a[2]*b[0])/denominateur_delta;
        /*if (fabs(denominateur_delta) < ( (fabs(a[0])+fabs(a[1])+fabs(a[2])+fabs(a[3]))*1e-12 )    )
            return 1e300;*/
        uvii[0]=uvi[0]+delta_u;
        uvii[1]=uvi[1]+delta_v;
        if (face->get_surface()->est_periodique_u()==1)
                {
                if(uvii[0]<0.) uvii[0]=face->get_surface()->get_periode_u()-uvii[0];
                if(uvii[0]>face->get_surface()->get_periode_u()) uvii[0]=uvii[0]-face->get_surface()->get_periode_u();
                }
        if (face->get_surface()->est_periodique_v()==1)
                {
                if(uvii[1]<0.) uvii[0]=face->get_surface()->get_periode_v()-uvii[1];
                if(uvii[1]>face->get_surface()->get_periode_v()) uvii[1]=uvii[1]-face->get_surface()->get_periode_v();
                }
        delta_u=uvii[0]-uvi[0];
        delta_v=uvii[1]-uvi[1];
        if (compteur>500) return 1e300;
        }

while ((fabs(delta_u)>precision)||(fabs(delta_v)>precision));
face->evaluer(uvii,xyz);
OT_VECTEUR_3D S(xyz[0],xyz[1],xyz[2]);
double distance=(S-Pt).get_longueur();
uv[0]=uvii[0];
uv[1]=uvii[1];
return distance;
}


void TOIBREP::etendrelevelset(FEM_SOLUTION* sol,int numsol)
{
sol->active_solution(numsol);


LISTE_FM know;
LISTE_FM trial;
LISTE_FM far;
LISTE_FM exterieur;
LISTE_FM_TRI trialtri;
LISTE_FM_TRI_ID trialtriid;


LISTE_FEM_ELEMENT3::iterator ittet;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittet);tet!=NULL;tet=mai->get_suivant_element3(ittet))
        {
        tet->change_solution(0.);
        int numsol=0;
        if (tet->get_fem_noeud(0)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(0)->change_numero(1);} else tet->get_fem_noeud(0)->change_numero(0);
        if (tet->get_fem_noeud(1)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(1)->change_numero(1);} else tet->get_fem_noeud(1)->change_numero(0);
        if (tet->get_fem_noeud(2)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(2)->change_numero(1);} else tet->get_fem_noeud(2)->change_numero(0);
        if (tet->get_fem_noeud(3)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(3)->change_numero(1);} else tet->get_fem_noeud(3)->change_numero(0);
        if (numsol==4)
                ajouter_liste(know,tet);
                
        else if (numsol==3)
                {
                if (tet->get_fem_noeud(0)->get_numero()==0) tet->change_numero(0);
                else if (tet->get_fem_noeud(1)->get_numero()==0) tet->change_numero(1);
                else if (tet->get_fem_noeud(2)->get_numero()==0) tet->change_numero(2);
                else tet->change_numero(3);
                ajouter_liste(trial,tet);
                } 
        else
                ajouter_liste(far,tet);
        
        }
for (LISTE_FM::iterator i=trial.begin();i!=trial.end();i++)
        {
        FEM_ELEMENT3* tet=*i;
        int signe;
        double sol=resoudgradT(tet,&signe);
        if (fabs(sol)>0.00000001)
                {
                if (fabs(sol)<fabs(tet->get_fem_noeud(tet->get_numero())->get_solution()))
                        {
                        tet->get_fem_noeud(tet->get_numero())->change_solution(sol);
                        ajouter_liste(trialtri,trialtriid,tet,fabs(tet->get_fem_noeud(tet->get_numero())->get_solution()));
                        }
                }
        else 
                {
                ajouter_liste(exterieur,tet);
                tet->change_solution(1e300);
                }
        }
int fin=0;
LISTE_FM_TRI::iterator itfin=trialtri.end();
itfin--;
double longref=(*itfin).first;
do
        {
        LISTE_FM_TRI::iterator it=trialtri.begin();
        FEM_ELEMENT3* tet=(*it).second;
        double longcourant=(*it).first;
        supprimer_liste(trialtri,trialtriid,tet);
        ajouter_liste(know,tet);
        FEM_NOEUD* noeud=tet->get_fem_noeud(tet->get_numero());
        noeud->change_numero(1);
        if (noeud->get_solution()>20000) 
                cout << "BUGGGGGGG" <<endl;
        int nbtetra=noeud->get_lien_element3()->get_nb();
        for (int i=0;i<nbtetra;i++)
                {
                FEM_ELEMENT3* tet2=noeud->get_lien_element3()->get(i);
                if (tet2==tet) continue;
                LISTE_FM_TRI_ID::iterator it=trialtriid.find(tet2->get_id());
                if (it!=trialtriid.end())
                        {
                        int signe;
                        double sol=resoudgradT(tet2,&signe);
                        double solution=tet2->get_fem_noeud(tet2->get_numero())->get_solution();
                        if (fabs(sol)>0.00000001)
                                if (!((solution<1e200)&&(sol*solution<0)))
                                        if (fabs(sol)<fabs(solution))
                                        {
                                        supprimer_liste(trialtri,trialtriid,tet2);
                                        ajouter_liste(trialtri,trialtriid,tet2,fabs(sol));
                                        tet2->get_fem_noeud(tet2->get_numero())->change_solution(sol);
                                        }
                        }
                LISTE_FM::iterator it2=find(far.begin(),far.end(),tet2);
                if (it2!=far.end())
                        {
                        int numsol=0;
                        if (tet2->get_fem_noeud(0)->get_numero()==1) numsol++;
                        if (tet2->get_fem_noeud(1)->get_numero()==1) numsol++;
                        if (tet2->get_fem_noeud(2)->get_numero()==1) numsol++;
                        if (tet2->get_fem_noeud(3)->get_numero()==1) numsol++;
                        //if (numsol==4) 
                                //cout << " BUG " <<endl;       
                        if (numsol==3)
                                {
                                if (tet2->get_fem_noeud(0)->get_numero()==0) tet2->change_numero(0);
                                else if (tet2->get_fem_noeud(1)->get_numero()==0) tet2->change_numero(1);
                                else if (tet2->get_fem_noeud(2)->get_numero()==0) tet2->change_numero(2);
                                else tet2->change_numero(3);
                                int signe;
                                double sol=resoudgradT(tet2,&signe);
                                double ancsol=tet2->get_fem_noeud(tet2->get_numero())->get_solution();
                                if (fabs(sol)>0.00000001)
                                        {
                                        if (!((ancsol<1e200) && (ancsol*sol<0.) ))
                                                {
                                                tet2->get_fem_noeud(tet2->get_numero())->change_solution(sol);
                                                supprimer_liste(far,tet2);
                                                ajouter_liste(trialtri,trialtriid,tet2,fabs(sol));
                                                }
                                        }
                                else 
                                        {
                                        tet2->change_solution(1e300);
                                        ajouter_liste(exterieur,tet2);
                                        }
                                }
                        }
                }
        if (trialtri.size()==0) fin=1;
        if (exterieur.size()>0)
        if (fin==0)
        if (longcourant>longref)
                {
                int nombre=exterieur.size();
                for (int i=0;i<nombre;i++)
                        {
                        FEM_ELEMENT3* tet2=(*(exterieur.begin()));
                        supprimer_liste(exterieur,tet2);
                        ajouter_liste(know,tet2);
                        for (int nd=0;nd<4;nd++)
                                {
                                FEM_NOEUD* noeud=tet2->get_fem_noeud(nd);
                                int nbtetra=noeud->get_lien_element3()->get_nb();
                                for (int i=0;i<nbtetra;i++)
                                        {
                                        FEM_ELEMENT3* tet3=noeud->get_lien_element3()->get(i);
                                        if (tet2==tet3) continue;
                                        LISTE_FM::iterator it2=find(far.begin(),far.end(),tet3);
                                        if (it2!=far.end())
                                                {
                                                supprimer_liste(far,tet3);
                                                ajouter_liste(exterieur,tet3);
                                                tet3->change_solution(1e300);
                                                }
                                        }
                                }
                        }
                LISTE_FM_TRI::iterator itfin=trialtri.end();
                itfin--;
                longref=(*itfin).first;
                }
        }
while (fin==0);
LISTE_FEM_NOEUD::iterator itnoeud;
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
        noeud->change_numero(0);
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittet);tet!=NULL;tet=mai->get_suivant_element3(ittet))
        {
        if (tet->get_solution()>1e200) continue;
        tet->get_fem_noeud(0)->change_numero(1);
        tet->get_fem_noeud(1)->change_numero(1);
        tet->get_fem_noeud(2)->change_numero(1);
        tet->get_fem_noeud(3)->change_numero(1);
        }
int i=0;        
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
        {
        if (noeud->get_numero()==1) sol->ecrire(i,numsol,noeud->get_solution()); else sol->ecrire(i,numsol,1e300);
        ++i;
        }
}
IBrep TOIBREP::exporter_IBrep(string chemin,FEM_SOLUTION* solution,FEM_SOLUTION *solution_ele,MG_GROUPE_TOPOLOGIQUE* mggt)
{
TPL_MAP_ENTITE<MG_VOLUME*> lst;
if (mggt==NULL)
        {
        int nbvol=geo->get_nb_mg_volume();
        for (int i=0;i<nbvol;i++)
                lst.ajouter(geo->get_mg_volume(i));
        }
else 
        {
        int nbmggt=mggt->get_nb();
        for (int i=0;i<nbmggt;i++)
                if (mggt->get(i)->get_dimension()==3) 
                        lst.ajouter((MG_VOLUME*)mggt->get(i));
        }
IBrep brep(100);
int nbvolexp=lst.get_nb();
for (int i=0;i<nbvolexp;i++)
        {
        MG_VOLUME* vol=lst.get(i);
        int nbcoq=vol->get_nb_mg_coquille();
        IVolumeN ivoltmp(vol->get_id());
        pair<IVolume*,bool> pair_ivol=brep.AddVolume(ivoltmp);
        if (vol->get_idoriginal()!="") brep.SetName(pair_ivol.first,vol->get_idoriginal());
        IVolume* ivol=pair_ivol.first;
        for (int j=0;j<nbcoq;j++)
                {
                MG_COQUILLE* coq=vol->get_mg_coquille(j);
                int nbface=coq->get_nb_mg_coface();
                IShell ishell(nbface);
                for (int k=0;k<nbface;k++)
                        {
                        MG_COFACE* coface=coq->get_mg_coface(k);
                        ishell[k]=coface->get_id();
                        MG_FACE* face=coface->get_face();
                        int sens=coface->get_orientation();
                        ICoFaceN icoface(coface->get_id(),face->get_id(),sens);
                        brep.AddCoFace(icoface);
                        IFace* iface=brep.GetFace(face->get_id());
                        if (iface==NULL)
                                {
                                IFaceN ifacenew(face->get_id());
                                pair<IFace*,bool> pair_iface=brep.AddFace(ifacenew);
                                if (face->get_idoriginal()!="") brep.SetName(pair_iface.first,face->get_idoriginal());
                                iface=pair_iface.first;
                                int nbbou=face->get_nb_mg_boucle();
                                for (int l=0;l<nbbou;l++)
                                        {
                                        MG_BOUCLE* bou=face->get_mg_boucle(l);
                                        int nbare=bou->get_nb_mg_coarete();
                                        ILoop iloop(nbare);
                                        for (int m=0;m<nbare;m++)
                                                {
                                                MG_COARETE* coare=bou->get_mg_coarete(m);
                                                iloop[m]=coare->get_id();
                                                MG_ARETE* are=coare->get_arete();
                                                int sens=coare->get_orientation();
                                                ICoEdgeN icoedge(coare->get_id(),are->get_id(),sens,face->get_id());
                                                brep.AddCoEdge(icoedge);
                                                IEdge* iedge=brep.GetEdge(are->get_id());
                                                IVertex *ivertex1,*ivertex2;
                                                ICoVertex *icover1,*icover2;
                                                if (iedge==NULL)
                                                        {
                                                        MG_COSOMMET* cover1=are->get_cosommet1();
                                                        MG_COSOMMET* cover2=are->get_cosommet2();
                                                        MG_SOMMET* ver1=cover1->get_sommet();
                                                        MG_SOMMET* ver2=cover2->get_sommet();
                                                        ICoVertexN icovertmp1(cover1->get_id(),are->get_id(),ver1->get_id(),cover1->get_t());
                                                        ICoVertexN icovertmp2(cover2->get_id(),are->get_id(),ver2->get_id(),cover2->get_t());
                                                        pair<ICoVertex*,bool> pair_icover1=brep.AddCoVertex(icovertmp1);
                                                        pair<ICoVertex*,bool> pair_icover2=brep.AddCoVertex(icovertmp2);
                                                        icover1=pair_icover1.first;
                                                        icover2=pair_icover2.first;
                                                        IEdgeN iedgetmp(are->get_id(),cover1->get_id(),cover2->get_id());
                                                        ivertex1=brep.GetVertex(ver1->get_id());
                                                        ivertex2=brep.GetVertex(ver2->get_id());
                                                        if (ivertex1==NULL)
                                                                {
                                                                MG_POINT* pt=ver1->get_point();
                                                                double xyz[3];
                                                                pt->evaluer(xyz);
                                                                IVertexN ivertex(ver1->get_id(),xyz);
                                                                pair<IVertex*,bool> pair_ivertex1=brep.AddVertex(ivertex);
                                                                if (ver1->get_idoriginal()!="") brep.SetName(pair_ivertex1.first,ver1->get_idoriginal());
                                                                ivertex1=pair_ivertex1.first;
                                                                }
                                                        if (ivertex2==NULL)
                                                                {
                                                                MG_POINT* pt=ver2->get_point();
                                                                double xyz[3];
                                                                pt->evaluer(xyz);
                                                                IVertexN ivertex(ver2->get_id(),xyz);
                                                                pair<IVertex*,bool> pair_ivertex2=brep.AddVertex(ivertex);
                                                                if (ver2->get_idoriginal()!="") brep.SetName(pair_ivertex2.first,ver2->get_idoriginal());
                                                                ivertex2=pair_ivertex2.first;
                                                                }
                                                        ivertex1->AddCoVertex(cover1->get_id());
                                                        ivertex2->AddCoVertex(cover2->get_id());
                                                        pair<IEdge*,bool> pair_iedge=brep.AddEdge(iedgetmp);
                                                        if (are->get_idoriginal()!="") brep.SetName(pair_iedge.first,are->get_idoriginal());
                                                        iedge=pair_iedge.first;
                                                        }
                                                iedge->AddCoEdge(coare->get_id(),brep);
                                                if (sens==1)
                                                        {
                                                        icover1=brep.GetCoVertex(iedge->FromCoVertex());
                                                        ivertex1=brep.GetVertex(icover1->Vertex());
                                                        ivertex1->AddFace(face->get_id());
                                                        }
                                                else
                                                        {
                                                        icover2=brep.GetCoVertex(iedge->ToCoVertex());
                                                        ivertex2=brep.GetVertex(icover2->Vertex());
                                                        ivertex2->AddFace(face->get_id());
                                                        }
                                                }
                                        iface->AddLoop(iloop);
                                        }
                                }
                        iface->AddCoFace(coface->get_id());
                        }
                ivol->AddShell(ishell);
                }
        }
int nbsol=solution->get_nb_champ();
//unsigned long *correspondface=new unsigned long[nbsol];
for (int i=0;i<nbsol;i++)
        {
        std::string nom=solution->get_legende(i);
        char nom2[2];
        unsigned long id;
        sscanf(nom.c_str(),"%s %lu",nom2,&id);
        std::pair<IField *,bool> fld=brep.AddField(id);
        if (nom2[0]!='F')
          {
          IField *f=fld.first;
          f->info.tag=1;
          }
        }
LISTE_FEM_NOEUD::iterator itnoeud;
int i=0;
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
        {
        int nbsolnoeud=0;
        for (int j=0;j<nbsol;j++) 
                {
                if (solution->lire(i,j)<1e200) nbsolnoeud++;    
                }
        double *xyz=noeud->get_coord();
        INodeN inoeudtmp(noeud->get_id(),xyz[0],xyz[1],xyz[2],nbsolnoeud);
        pair<INode*,bool> pair_inoeud=brep.AddNode(inoeudtmp);
        INode* inoeud=pair_inoeud.first;
        int num=0;
        IBrep::iterator_fields itf=brep.begin_fields();
        for (int j=0;j<nbsol;j++) 
                {
                if (solution->lire(i,j)<1e200)
                        {
                        inoeud->Fields()[num]=(*itf).num;
                        inoeud->Values()[num]=solution->lire(i,j);
                        num++;
                        }
                itf++;
                }
        i++;
        }
LISTE_FEM_ELEMENT3::iterator ittet;
int num=0;
i=0;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(ittet);tet!=NULL;tet=mai->get_suivant_element3(ittet))
        {
        int nbsolele=0;
        for (int j=0;j<nbsol;j++) 
                {
                if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001) nbsolele++;  
                }
        IElementN xfemele(IElement::TETRAHEDRON,nbsolele);
        num++;
        xfemele.num=num;
        xfemele.GetNode(0)=tet->get_fem_noeud(0)->get_id();
        xfemele.GetNode(1)=tet->get_fem_noeud(1)->get_id();
        xfemele.GetNode(2)=tet->get_fem_noeud(2)->get_id();
        xfemele.GetNode(3)=tet->get_fem_noeud(3)->get_id();
        int num=0;
        IBrep::iterator_fields itf=brep.begin_fields();
        for (int j=0;j<nbsol;j++) 
                {
                if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001)
                  {
                  xfemele.Fields()[num]=(*itf).num;
                  num++;
                  }
                itf++;
                }
        brep.AddElement(xfemele);
        i++;
        }


std::ofstream output(chemin.c_str());
output << brep << std::endl;
output.close();


return brep;
}


void TOIBREP::ajouter_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_ELEMENT3* tet,double val)
{
pair<double,FEM_ELEMENT3*> p(val,tet);
LISTE_FM_TRI::iterator it=lst.insert(p);
lstid[tet->get_id()]=it;
}

void TOIBREP::supprimer_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_ELEMENT3* tet)
{
LISTE_FM_TRI::iterator it2=lstid[tet->get_id()];
LISTE_FM_TRI_ID::iterator it3=lstid.find(tet->get_id());
lstid.erase(it3);
lst.erase(it2);
}

void TOIBREP::ajouter_liste(LISTE_FM& lst,FEM_ELEMENT3* tet)
{
lst.push_back(tet);
}

void TOIBREP::supprimer_liste(LISTE_FM& lst,FEM_ELEMENT3* tet)
{
LISTE_FM::iterator it=find(lst.begin(),lst.end(),tet);
if (it!=lst.end()) lst.erase(it);
}


double TOIBREP::resoudgradT(FEM_ELEMENT3* tet,int *signe)
{
double j[9];
double N[12];
double jN[12]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
double T[4];
double uv[2]={0.25,0.25};


tet->get_inverse_jacob(j,uv);
for (int i=0;i<3;i++)
        for (int k=0;k<4;k++)
                N[i*4+k]=tet->get_fonction_derive_interpolation(k+1,i+1,uv);
int premier=0;
double tmin=1e300;
double tmax=-1e300;
for (int i=0;i<4;i++)
        {
        if (i==tet->get_numero()) continue;
        T[i]=tet->get_fem_noeud(i)->get_solution();
        if (fabs(T[i])>0.000001)
                {
                if (premier==0)
                        if (T[i]>0) (*signe)=1; else (*signe)=-1;
                else if (T[i]*(*signe)<0) (*signe)=0;
                }
        T[i]=fabs(T[i]);
        if (tet->get_numero()!=i)
                {
                if (T[i]<tmin) tmin=T[i];
                if (T[i]>tmax) tmax=T[i];
                }
        premier=1;
        }
for (int i=0;i<3;i++)
        for (int k=0;k<4;k++)
                for (int l=0;l<3;l++)
                        jN[i*4+k]=jN[i*4+k]+j[i*3+l]*N[l*4+k];
double a=0.,b=0.,c=-1.;
for (int i=0;i<3;i++)
        {
        double coef=0.;
        double coefinc=0.;
        for (int l=0;l<4;l++)
                {
                if (tet->get_numero()!=l) coef=coef+jN[i*4+l]*T[l];
                else coefinc=coefinc+jN[i*4+l];
                }
        c=c+coef*coef;
        a=a+coefinc*coefinc;
        b=b+2*coef*coefinc;
        }
/*if (*signe==0)
        cout << "attention " <<endl;*/
double det=b*b-4.*a*c;
if (det<0.) det=0.; else det=sqrt(det);
double sol1=(-b-det)/2./a;
double sol2=(-b+det)/2./a;
double sol=sol1;
if (sol2>sol1) sol=sol2;
if (sol<tmin*0.99)
        sol=0.;
sol=sol*(*signe);
return sol;
}
Revision:	346
Committed:	Fri Jun 29 15:32:53 2012 UTC (12 years, 10 months ago) by francois
File size:	70449 byte(s)
Log Message:	un oubli dans les derniers commits.