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 "fem_segment2.h"
#include "xfem_tetra4.h"
#include "xfem_triangle3.h"
#include "xfem_segment2.h"
#include "mg_maillage_outils.h"

class inter_ele_arete
{
public:
FEM_ELEMENT3* ele;
XFEM_ELEMENT3* xele;
MG_ARETE* are;
int status1;
int status2;
double xyz1[3];
double xyz2[3];
FEM_NOEUD* no1;
FEM_NOEUD* no2;
};


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;
}


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);
    }
      
}

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


void TOIBREP::oriente_tri(XFEM_ELEMENT2* tri,MG_FACE* face)
{
FEM_NOEUD* no1=tri->get_fem_noeud(0);
FEM_NOEUD* no2=tri->get_fem_noeud(1);
FEM_NOEUD* no3=tri->get_fem_noeud(2);
OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
OT_VECTEUR_3D normal=n1n2&n1n3;
double *xyz=no1->get_coord();
double uv[0];
face->inverser(uv,xyz);
double n[3];
face->calcul_normale_unitaire(uv,n);
OT_VECTEUR_3D dir(n);
normal.norme();
double ps=normal*dir;
if (ps<0) 
    {
    tri->change_noeud(1,no3);  
    tri->change_noeud(2,no2);  
    }
}

#ifdef ENABLE_IBREP
//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",MAGIC::ENTITE_SOLUTION::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",MAGIC::ENTITE_SOLUTION::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(1e300,i,j);
        lst_noeud.ajouter(noeud);
        i++;
        }
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base recherche boite");
octree_tetra.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_tetra.inserer(tet);
  for (int j=0;j<nbvraiface+lstaretetan.get_nb();j++)
      solution_ele->ecrire(0,i,j);
  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;
}
#endif


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++;
  }
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(1,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(1,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_tetra.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)==1) continue;
          int etat=outilfem.estdansletetra(tet,xyz[0],xyz[1],xyz[2]);
          if ((outilfem.compare_etat_tetra(etat,FEM_MAILLAGE_OUTILS::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=outilfem.estdansletetra(tet,xyztmp[0],xyztmp[1],xyztmp[2]);
                    if (outilfem.compare_etat_tetra(etattmp,FEM_MAILLAGE_OUTILS::SUR_FACE)) {etat=FEM_MAILLAGE_OUTILS::STRICTINTERIEUR;break;}
                    if (outilfem.compare_etat_tetra(etattmp,FEM_MAILLAGE_OUTILS::INTERIEUR)) {etat=FEM_MAILLAGE_OUTILS::STRICTINTERIEUR;break;}
                    if (outilfem.compare_etat_tetra(etattmp,FEM_MAILLAGE_OUTILS::STRICTINTERIEUR)) {etat=FEM_MAILLAGE_OUTILS::STRICTINTERIEUR;break;}
                    }
                }
          if (outilfem.compare_etat_tetra(etat,FEM_MAILLAGE_OUTILS::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,1);  //temp
                        tet->get_fem_noeud(k)->change_etat(1,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)==1) solution->ecrire(noeud->get_solution(),i,numsol);
        else solution->ecrire(1e300,i,numsol); 
        ++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)==1) solution_ele->ecrire(1.,i,numsol);
        ++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)
void TOIBREP::levelsetn(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lsttopo,MG_FACE *face,TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lstentitetopoface)
{
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"          Echantillonnage");
          affiche(mess);
          }
        
//solution->active_solution(numsol);
//solution_ele->active_solution(numsol);
vector<TOIBREP_POINT*> lst;
BOITE_2D boite=ot.get_boite_2D(face);
echantillonne_aretes(lst,face);
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(1,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(1,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,lstentitetopoface);
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(noeud->get_solution(),i,numsol);
        else solution->ecrire(1e300,i,numsol); 
        ++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(1.,i,numsol);
        ++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(1)!=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)==1) || (no1->get_etat(1)==-1))
            {
            nbactif++;
            if (no1->get_solution()>0.) nbpositif++;
            }
      if ((no2->get_etat(1)==1) || (no2->get_etat(1)==-1))
            {
            nbactif++;
            if (no2->get_solution()>0.) nbpositif++;
            }
      if ((no3->get_etat(1)==1) || (no3->get_etat(1)==-1))
            {
            nbactif++;
            if (no3->get_solution()>0.) nbpositif++;
            }
      if ((no4->get_etat(1)==1) || (no4->get_etat(1)==-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(1)==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];}
      if (no2->get_etat(1)==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];}
      if (no3->get_etat(1)==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];}
      if (no4->get_etat(1)==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(1,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(1)!=0)
      if (no2->get_etat(1)!=0)
      if (no3->get_etat(1)!=0)
      if (no4->get_etat(1)!=0)
      if (tet->get_etat(1)==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,1);
                   no1->change_etat(1,1);
                   no2->change_etat(1,1);
                   no3->change_etat(1,1);
                   no4->change_etat(1,1);
                   reactive=1;
                   couche_active=1;
                   
                 }
                
              }
           if (!reactive)
              {
              if (no1->get_etat(1)==2) no1->change_etat(1,-1);
              if (no2->get_etat(1)==2) no2->change_etat(1,-1);
              if (no3->get_etat(1)==2) no3->change_etat(1,-1);
              if (no4->get_etat(1)==2) no4->change_etat(1,-1);
              tet->change_etat(1,-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)!=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(1,0);
      }
}


void TOIBREP::remplir_trou(std::vector<TOIBREP_POINT*> *lst,MG_FACE* face,int sens,TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lstentitetopoface)
{
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(1)!=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)==1) || (no1->get_etat(1)==-1))
            {
            nbactif++;
            if (no1->get_solution()>0.) nbpositif++;
            }
      if ((no2->get_etat(1)==1) || (no2->get_etat(1)==-1))
            {
            nbactif++;
            if (no2->get_solution()>0.) nbpositif++;
            }
      if ((no3->get_etat(1)==1) || (no3->get_etat(1)==-1))
            {
            nbactif++;
            if (no3->get_solution()>0.) nbpositif++;
            }
      if ((no4->get_etat(1)==1) || (no4->get_etat(1)==-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(1)==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];}
            if (no2->get_etat(1)==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];}
            if (no3->get_etat(1)==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];}
            if (no4->get_etat(1)==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(1,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(1)!=0)
      if (no2->get_etat(1)!=0)
      if (no3->get_etat(1)!=0)
      if (no4->get_etat(1)!=0)
      if (tet->get_etat(1)==0)
          {
          castot++;
          double lstptdecoup[5*3];
          int reactive=0;
          for (int num=0;num<2;num++)
              for (int i=0;i<tet->get_nb_xfem(num);i++)
                {
                if (reactive) continue;
                MG_ELEMENT_TOPOLOGIQUE* mgt=tet->get_xfem(num,i)->get_lien_topologie();
                if (lstentitetopoface->existe(mgt)==1)
                    {
                    cas++;
                    tet->change_etat(1,1);
                    no1->change_etat(1,1);
                    no2->change_etat(1,1);
                    no3->change_etat(1,1);
                    no4->change_etat(1,1);
                    reactive=1;
                    couche_active=1;
                    }
                }
          int nb=0;
          decoupe_tetra_noeud(tet,lstptdecoup,&nb);
          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-10)
                  {
                  cas++;
                  tet->change_etat(1,1);
                  no1->change_etat(1,1);
                  no2->change_etat(1,1);
                  no3->change_etat(1,1);
                  no4->change_etat(1,1);
                  reactive=1;
                  couche_active=1;
                  }
                }
              }
           if (!reactive)
              {
              if (no1->get_etat(1)==2) no1->change_etat(1,-1);
              if (no2->get_etat(1)==2) no2->change_etat(1,-1);
              if (no3->get_etat(1)==2) no3->change_etat(1,-1);
              if (no4->get_etat(1)==2) no4->change_etat(1,-1);
              tet->change_etat(1,-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)!=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(1,0);
      }
if (castot!=0)
  {
  double res=cas*100.0/castot;
  char mess[255];
  sprintf(mess,"                  %lf%% en %d couches",res,nb_couche);
  affiche(mess);
  }
  else
  {
  char mess[255];
  sprintf(mess,"                  pas de trou");
  affiche(mess);
  }
}

int TOIBREP::decoupe_tetra(FEM_ELEMENT3* tet,MG_FACE* face)
{
  std::vector<FEM_NOEUD*> lstno;
  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());
        if (t<0.0000000001) lstno.push_back(no1);
        else if (t>1.-0.0000000001) lstno.push_back(no2);
        else 
          {
          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=inserer_noeud(face,x,y,z);
          no->change_solution(0.0);
          no->change_etat(0,INCONNU);
          lstno.push_back(no);
          }
        }
  if (no1->get_solution()*no3->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no3->get_solution());
        if (t<0.0000000001) lstno.push_back(no1);
        else if (t>1.-0.0000000001) lstno.push_back(no3);
        else 
          {
          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=inserer_noeud(face,x,y,z);
          no->change_solution(0.0);
          no->change_etat(0,INCONNU);
          lstno.push_back(no);
          }
        }
   if (no1->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no1->get_solution()/(no1->get_solution()-no4->get_solution());
        if (t<0.0000000001) lstno.push_back(no1);
        else if (t>1.-0.0000000001) lstno.push_back(no4);
        else 
          {
          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=inserer_noeud(face,x,y,z);
          no->change_solution(0.0);
          no->change_etat(0,INCONNU);
          lstno.push_back(no);
          }
        }
  if (no2->get_solution()*no3->get_solution()<-0.00000001)
        {
        double t=no2->get_solution()/(no2->get_solution()-no3->get_solution());
        if (t<0.0000000001) lstno.push_back(no2);
        else if (t>1.-0.0000000001) lstno.push_back(no3);
        else 
          {
          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=inserer_noeud(face,x,y,z);
          no->change_solution(0.0);
          no->change_etat(0,INCONNU);
          lstno.push_back(no);
          }
          
        }
    if (no2->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no2->get_solution()/(no2->get_solution()-no4->get_solution());
        if (t<0.0000000001) lstno.push_back(no2);
        else if (t>1.-0.0000000001) lstno.push_back(no4);
        else 
          {
          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=inserer_noeud(face,x,y,z);
          no->change_solution(0.0);
          no->change_etat(0,INCONNU);
          lstno.push_back(no);
          }
          
        }
    if (no3->get_solution()*no4->get_solution()<-0.00000001)
        {
        double t=no3->get_solution()/(no3->get_solution()-no4->get_solution());
        if (t<0.0000000001) lstno.push_back(no3);
        else if (t>1.-0.0000000001) lstno.push_back(no4);
        else 
          {
          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=inserer_noeud(face,x,y,z);
          no->change_solution(0.0);
          no->change_etat(0,INCONNU);
          lstno.push_back(no);
          }
          
          
        }
int nb=lstno.size();
if (nb==3)
      {
      FEM_NOEUD *lstnoeud[3];
      lstnoeud[0]=lstno[0];
      lstnoeud[1]=lstno[1];
      lstnoeud[2]=lstno[2];
      XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
      //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
      //oriente_tri(xtri,face);
      //mai->ajouter_xfem_element2(xtri);
//      XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
      xtri->change_etat(XINCONNU);
      tet->change_etat(0,FRONTIERE);
      return 1;
      }
if (nb==4)
      {
        if (test_du_point_milieu(lstno[0],lstno[1],tet)==1)
            {
            FEM_NOEUD *lstnoeud[3];
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[1];
            lstnoeud[2]=lstno[2];
            //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri,face);
            //mai->ajouter_xfem_element2(xtri);
            XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
            xtri->change_etat(XINCONNU);
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[1];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri2=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri2=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri2,face);
            //mai->ajouter_xfem_element2(xtri2);
            xtri2->change_etat(XINCONNU);
            tet->change_etat(0,FRONTIERE);
            }
         else if (test_du_point_milieu(lstno[0],lstno[2],tet)==1)
            {
            FEM_NOEUD *lstnoeud[3];
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[1];
            lstnoeud[2]=lstno[2];
            //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri,face);
            //mai->ajouter_xfem_element2(xtri);
            XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
            xtri->change_etat(XINCONNU);
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[2];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri2=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri2=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri2,face);
            //mai->ajouter_xfem_element2(xtri2);
            xtri2->change_etat(XINCONNU);
            tet->change_etat(0,FRONTIERE);
            }
        else if (test_du_point_milieu(lstno[0],lstno[3],tet)==1)
            {
            FEM_NOEUD *lstnoeud[3];
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[1];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri,face);
            //mai->ajouter_xfem_element2(xtri);
            xtri->change_etat(XINCONNU);
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[2];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri2=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri2=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri2,face);
            //mai->ajouter_xfem_element2(xtri2);
            xtri2->change_etat(XINCONNU);
            tet->change_etat(0,FRONTIERE);
            }
        else if (test_du_point_milieu(lstno[1],lstno[2],tet)==1)
            {
            FEM_NOEUD *lstnoeud[3];
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[1];
            lstnoeud[2]=lstno[2];
            XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri,face);
            //mai->ajouter_xfem_element2(xtri);
            xtri->change_etat(XINCONNU);
            lstnoeud[0]=lstno[1];
            lstnoeud[1]=lstno[2];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri2=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri2=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri2,face);
            //mai->ajouter_xfem_element2(xtri2);
            xtri2->change_etat(XINCONNU);
            tet->change_etat(0,FRONTIERE);
            }
         else if (test_du_point_milieu(lstno[1],lstno[3],tet)==1)
            {
            FEM_NOEUD *lstnoeud[3];
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[1];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri,face);
            //mai->ajouter_xfem_element2(xtri);
            xtri->change_etat(XINCONNU);
            lstnoeud[0]=lstno[1];
            lstnoeud[1]=lstno[2];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri2=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri2=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri2,face);
            //mai->ajouter_xfem_element2(xtri2);
            xtri2->change_etat(XINCONNU);
            tet->change_etat(0,FRONTIERE);
            }
         else
            {
            FEM_NOEUD *lstnoeud[3];
            lstnoeud[0]=lstno[0];
            lstnoeud[1]=lstno[2];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri,face);
            //mai->ajouter_xfem_element2(xtri);
            xtri->change_etat(XINCONNU);
            lstnoeud[0]=lstno[1];
            lstnoeud[1]=lstno[2];
            lstnoeud[2]=lstno[3];
            XFEM_TRIANGLE3* xtri2=inserer_xtriangle(tet,face,lstnoeud);
            //XFEM_TRIANGLE3* xtri2=new XFEM_TRIANGLE3(tet,face,lstnoeud);
            //oriente_tri(xtri2,face);
            //mai->ajouter_xfem_element2(xtri2);
            xtri2->change_etat(XINCONNU);
            tet->change_etat(0,FRONTIERE);
            }
      return 2;      
      }
}


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_tetra.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)==1) continue;
                        int etat=outilfem.estdansletetra(tet,xyz[0],xyz[1],xyz[2]);
                        if ((outilfem.compare_etat_tetra(etat,FEM_MAILLAGE_OUTILS::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=outilfem.estdansletetra(tet,xyztmp[0],xyztmp[1],xyztmp[2]);
                              if (outilfem.compare_etat_tetra(etattmp,FEM_MAILLAGE_OUTILS::SUR_FACE)) {etat=FEM_MAILLAGE_OUTILS::STRICTINTERIEUR;break;}
                              if (outilfem.compare_etat_tetra(etattmp,FEM_MAILLAGE_OUTILS::INTERIEUR)) {etat=FEM_MAILLAGE_OUTILS::STRICTINTERIEUR;break;}
                              if (outilfem.compare_etat_tetra(etattmp,FEM_MAILLAGE_OUTILS::STRICTINTERIEUR)) {etat=FEM_MAILLAGE_OUTILS::STRICTINTERIEUR;break;}
                              }
                          }
                        if (outilfem.compare_etat_tetra(etat,FEM_MAILLAGE_OUTILS::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 ((outilfem.compare_etat_tetra(etat,FEM_MAILLAGE_OUTILS::STRICTINTERIEUR)) && (pt->get_interieur()))
                                {
                                tet->change_etat(1,1);
                                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                                  {
                                  tet->get_fem_noeud(k)->change_etat(1,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(noeud->get_solution(),i,numsol); else sol->ecrire(1e300,i,numsol);
        ++i;
        }
}
#ifdef ENABLE_IBREP
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;
}
#endif


int TOIBREP::intersection_segment_triangle(FEM_NOEUD* nos1,FEM_NOEUD* nos2,FEM_NOEUD* no1,FEM_NOEUD *no2,FEM_NOEUD *no3,double *uv,double *xyz)
{
OT_VECTEUR_3D ori(no1->get_coord());
OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
OT_VECTEUR_3D vecdir(nos1->get_coord(),nos2->get_coord());
OT_VECTEUR_3D moinsvecdir=(-1.)*vecdir;
OT_VECTEUR_3D f(no1->get_coord(),nos1->get_coord());
OT_MATRICE_3D A(vec1,vec2,moinsvecdir);  
OT_MATRICE_3D Bu(f,vec2,moinsvecdir);  
OT_MATRICE_3D Bv(vec1,f,moinsvecdir);  
OT_MATRICE_3D Bt(vec1,vec2,f);  
double det=A.get_determinant();
if (fabs(det)<1e-10) return 0;
uv[0]=Bu.get_determinant()/det;
uv[1]=Bv.get_determinant()/det;
uv[2]=Bt.get_determinant()/det;
if (uv[2]<1e-10) return 0;
if (uv[2]>1.-1e-10) return 0;
if (uv[0]<-1e-10) return 0;
if (uv[0]>1+1e-10) return 0;
if (uv[1]<-1e-10) return 0;
if (uv[1]>1+1e-10) return 0;
if (uv[0]+uv[1]<-1e-10) return 0;
if (uv[0]+uv[1]>1+1e-10) return 0;
xyz[0]=(nos1->get_coord())[0]+uv[2]*vecdir.get_x();
xyz[1]=(nos1->get_coord())[1]+uv[2]*vecdir.get_y();
xyz[2]=(nos1->get_coord())[2]+uv[2]*vecdir.get_z();
return 1;
}


int TOIBREP::intersection_arete_triangle(MG_ARETE* are,FEM_NOEUD* no1,FEM_NOEUD *no2,FEM_NOEUD *no3,double *uv,double *xyz)
{
OT_VECTEUR_3D ori(no1->get_coord());
OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
uv[0]=0.;  
uv[1]=0.;  
double tdep=uv[2];
if (uv[2]>are->get_tmax()) return 0;
double titer=are->get_tmax()-tdep;
int ok=0;
int nbiter=0;
do 
  {
  nbiter++;
  if (nbiter>40) ok=2;
  double arexyz[3],darexyz[3];
  are->evaluer(uv[2],arexyz);
  are->deriver(uv[2],darexyz);
  double f1=ori.get_x()+uv[0]*vec1.get_x()+uv[1]*vec2.get_x()-arexyz[0];
  double f2=ori.get_y()+uv[0]*vec1.get_y()+uv[1]*vec2.get_y()-arexyz[1];
  double f3=ori.get_z()+uv[0]*vec1.get_z()+uv[1]*vec2.get_z()-arexyz[2];
  OT_VECTEUR_3D vecdare(darexyz);
  vecdare=(-1.)*vecdare;
  OT_VECTEUR_3D f(-f1,-f2,-f3);
  OT_MATRICE_3D A(vec1,vec2,vecdare);  
  OT_MATRICE_3D Bu(f,vec2,vecdare);  
  OT_MATRICE_3D Bv(vec1,f,vecdare);  
  OT_MATRICE_3D Bt(vec1,vec2,f);  
  double det=A.get_determinant();
  if (fabs(det)<1e-10) 
          {
            if (nbiter==1)
              {
              uv[2]=tdep+0.05*titer;
              continue;
              }
            else ok=2;
          }
  double du=Bu.get_determinant()/det;
  double dv=Bv.get_determinant()/det;
  double dt=Bt.get_determinant()/det;
  uv[0]=uv[0]+du;  
  uv[1]=uv[1]+dv;  
  uv[2]=uv[2]+dt; 
  if (are->get_courbe()->est_periodique())
    if (uv[2]<are->get_tmin()) uv[2]=uv[2]+are->get_courbe()->get_periode();
  if (fabs(du)<1e-11)
    if (fabs(dv)<1e-11)
      if (fabs(dt)<1e-11) ok=1;
  }
while (ok==0);
//if (ok==2) {uv[2]=tdep+0.1;return intersection_arete_triangle(are,no1,no2,no3,uv,xyz);}
if (ok==2) return 0;
if (uv[2]<tdep+1e-10) return 0;
if (uv[2]>are->get_tmax()-1e-10) return 0;
if (uv[0]<-1e-10) return 0;
if (uv[0]>1+1e-10) return 0;
if (uv[1]<-1e-10) return 0;
if (uv[1]>1+1e-10) return 0;
if (uv[0]+uv[1]<-1e-10) return 0;
if (uv[0]+uv[1]>1+1e-10) return 0;
are->evaluer(uv[2],xyz);
return 1;
}

void TOIBREP::decoupe_element2_par_element1et0(FEM_ELEMENT3* ele,TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lstentitetopo)
{
TPL_MAP_ENTITE<FEM_NOEUD*> lst;
int nb=ele->get_nb_xfem(0);
for (int i=0;i<nb;i++)
  if (lstentitetopo->existe(ele->get_xfem(0,i)->get_lien_topologie()))
        lst.ajouter(((XFEM_ELEMENT0*)ele->get_xfem(0,i))->get_fem_noeud(0));
nb=ele->get_nb_xfem(1);
for (int i=0;i<nb;i++)
 if (lstentitetopo->existe(ele->get_xfem(1,i)->get_lien_topologie()))
  {
  lst.ajouter(((XFEM_ELEMENT1*)ele->get_xfem(1,i))->get_fem_noeud(0));
  lst.ajouter(((XFEM_ELEMENT1*)ele->get_xfem(1,i))->get_fem_noeud(1));
  }
TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR it;
for (FEM_NOEUD* no=lst.get_premier(it);no!=NULL;no=lst.get_suivant(it))
  decoupe_element2(ele,no);
  
}

int TOIBREP::decoupe_element2(FEM_ELEMENT3* ele,FEM_NOEUD* no)
{
TPL_LISTE_ENTITE<XFEM_ELEMENT2*> lst;
int nb=ele->get_nb_xfem(2);
for (int i=0;i<nb;i++)
  lst.ajouter((XFEM_ELEMENT2*)ele->get_xfem(2,i));
for (int i=0;i<nb;i++)
  decoupe_element(lst.get(i),no);
  
}

int TOIBREP::decoupe_element(XFEM_ELEMENT2 *ele,FEM_NOEUD* no)
{
if (no->get_id()==598016)
  std::cout << std::endl;
if (ele->get_lien_topologie()->est_topologie_sousjacente(no->get_lien_topologie())==0) return 0;
int status=outilfem.projeteestdansletriangle(ele,no->get_x(),no->get_y(),no->get_z());
if (status==0) return 0;
int arete1=outilfem.compare_etat_triangle(status,FEM_MAILLAGE_OUTILS::ARETE1);
int arete2=outilfem.compare_etat_triangle(status,FEM_MAILLAGE_OUTILS::ARETE2);
int arete3=outilfem.compare_etat_triangle(status,FEM_MAILLAGE_OUTILS::ARETE3);
int aretetot=arete1+arete2+arete3;
if (aretetot==2) 
       return 0;
FEM_NOEUD *lstno[3];
if (arete1==0)
    {
    lstno[0]=ele->get_fem_noeud(0);
    lstno[1]=no;
    lstno[2]=ele->get_fem_noeud(1);
    /*
    XFEM_TRIANGLE3 *xtri=new XFEM_TRIANGLE3(ele->get_fem_element_maillage(),ele->get_lien_topologie(),lstno);
    mai->ajouter_xfem_element2(xtri);
    oriente_tri(xtri,(MG_FACE*)ele->get_lien_topologie());*/
    XFEM_TRIANGLE3 *xtri=inserer_xtriangle(ele->get_fem_element_maillage(),(MG_FACE*)ele->get_lien_topologie(),lstno);
    if (xtri!=NULL)
      {
      xtri->change_etat(ele->get_etat());
      ele->get_fem_element_maillage()->change_etat(0,FRONTIERE);  
      }
    }
if (arete2==0)
    {
    lstno[0]=ele->get_fem_noeud(1);
    lstno[1]=no;
    lstno[2]=ele->get_fem_noeud(2);
    //XFEM_TRIANGLE3 *xtri=new XFEM_TRIANGLE3(ele->get_fem_element_maillage(),ele->get_lien_topologie(),lstno);
    //mai->ajouter_xfem_element2(xtri);
    //oriente_tri(xtri,(MG_FACE*)ele->get_lien_topologie());
    //xtri->change_etat(ele->get_etat());
    //ele->get_fem_element_maillage()->change_etat(0,FRONTIERE);  
    XFEM_TRIANGLE3 *xtri=inserer_xtriangle(ele->get_fem_element_maillage(),(MG_FACE*)ele->get_lien_topologie(),lstno);
    if (xtri!=NULL)
      {
      xtri->change_etat(ele->get_etat());
      ele->get_fem_element_maillage()->change_etat(0,FRONTIERE);  
      }
    }
if (arete3==0)
    {
    lstno[0]=ele->get_fem_noeud(2);
    lstno[1]=no;
    lstno[2]=ele->get_fem_noeud(0);
    //XFEM_TRIANGLE3 *xtri=new XFEM_TRIANGLE3(ele->get_fem_element_maillage(),ele->get_lien_topologie(),lstno);
    //mai->ajouter_xfem_element2(xtri);
    //oriente_tri(xtri,(MG_FACE*)ele->get_lien_topologie());
    //xtri->change_etat(ele->get_etat());
    //ele->get_fem_element_maillage()->change_etat(0,FRONTIERE);  
    XFEM_TRIANGLE3 *xtri=inserer_xtriangle(ele->get_fem_element_maillage(),(MG_FACE*)ele->get_lien_topologie(),lstno);
    if (xtri!=NULL)
      {
      xtri->change_etat(ele->get_etat());
      ele->get_fem_element_maillage()->change_etat(0,FRONTIERE);  
      }
    }
mai->supprimer_xfem_element2id(ele->get_id());
return 1;
}


int TOIBREP::decoupe_element(int status,FEM_ELEMENT3 *ele,FEM_NOEUD* no,FEM_NOEUD* no1,FEM_NOEUD* no2,FEM_NOEUD* no3,FEM_NOEUD* no4)
{
int face1=outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::FACE1);
int face2=outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::FACE2);
int face3=outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::FACE3);
int face4=outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::FACE4);
int facetot=face1+face2+face3+face4;
if (facetot==3) 
       return 0;
FEM_NOEUD *lstno[4];
if (face3==0)
  {
  lstno[0]=no1;
  lstno[1]=no;
  lstno[2]=no3;
  lstno[3]=no4;
  XFEM_TETRA4* xtet1=new XFEM_TETRA4(ele,ele->get_lien_topologie(),lstno);
  mai->ajouter_xfem_element3(xtet1);
  xtet1->change_etat(XINCONNU);
  ele->change_etat(0,FRONTIERE);
  }
if (face2==0)
  {
  lstno[0]=no1;
  lstno[1]=no2;
  lstno[2]=no;
  lstno[3]=no4;
  XFEM_TETRA4* xtet2=new XFEM_TETRA4(ele,ele->get_lien_topologie(),lstno);
  mai->ajouter_xfem_element3(xtet2);
  xtet2->change_etat(XINCONNU);
  ele->change_etat(0,FRONTIERE);
  }
if (face1==0)
  { 
  lstno[0]=no1;
  lstno[1]=no2;
  lstno[2]=no3;
  lstno[3]=no;
  XFEM_TETRA4* xtet3=new XFEM_TETRA4(ele,ele->get_lien_topologie(),lstno);
  mai->ajouter_xfem_element3(xtet3);
  xtet3->change_etat(XINCONNU);
  ele->change_etat(0,FRONTIERE);
  }
if (face4==0)
  { 
  lstno[0]=no;
  lstno[1]=no2;
  lstno[2]=no3;
  lstno[3]=no4;
  XFEM_TETRA4* xtet4=new XFEM_TETRA4(ele,ele->get_lien_topologie(),lstno);
  mai->ajouter_xfem_element3(xtet4);
  xtet4->change_etat(XINCONNU);
  ele->change_etat(0,FRONTIERE);
  }
return 1;
}


void TOIBREP::decoupe_element(int status,FEM_NOEUD* no,FEM_ELEMENT3* ele)
{
double *xyz=no->get_coord();
if (ele->get_nb_xfem(3)==0)
        decoupe_element(status,ele,no,ele->get_fem_noeud(0),ele->get_fem_noeud(1),ele->get_fem_noeud(2),ele->get_fem_noeud(3));
else 
    {
    int nb=ele->get_nb_xfem(3);  
    for (int i=nb-1;i>-1;i--)
      {
      XFEM_ELEMENT3* xele3=(XFEM_ELEMENT3*)ele->get_xfem(3,i);
      int status=outilfem.estdansletetra(xele3,xyz[0],xyz[1],xyz[2]);  
      if (outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::INTERIEUR)==1)  
        {
        int res=decoupe_element(status,ele,no,xele3->get_fem_noeud(0),xele3->get_fem_noeud(1),xele3->get_fem_noeud(2),xele3->get_fem_noeud(3));
        if (res==1) mai->supprimer_xfem_element3id(xele3->get_id());  
        }
        
      }
      if (ele->verifie_validite_decoupage_xfem()==0)
                    std::cout << "pb " << ele->get_id() << " " << nb <<"," <<  ele->get_nb_xfem(3) << " xfem" << std::endl;
        
    }
}

void TOIBREP::decoupe_element(XFEM_ELEMENT1* xseg,inter_ele_arete *ele_inter)
{
std::map<double,FEM_NOEUD*> lstnoeud;
for (int i=0;i<ele_inter->ele->get_nb_xfem(3);i++)
      {
      XFEM_TETRA4* xtet=(XFEM_TETRA4*)ele_inter->ele->get_xfem(3,i);
      double xyz1[3],xyz2[3],xyz3[3],xyz4[3],uv1[3],uv2[3],uv3[3],uv4[3];
      int inter_face1=intersection_segment_triangle(xseg->get_fem_noeud(0),xseg->get_fem_noeud(1),xtet->get_fem_noeud(0),xtet->get_fem_noeud(1),xtet->get_fem_noeud(2),uv1,xyz1);
      int inter_face2=intersection_segment_triangle(xseg->get_fem_noeud(0),xseg->get_fem_noeud(1),xtet->get_fem_noeud(0),xtet->get_fem_noeud(1),xtet->get_fem_noeud(3),uv2,xyz2);
      int inter_face3=intersection_segment_triangle(xseg->get_fem_noeud(0),xseg->get_fem_noeud(1),xtet->get_fem_noeud(0),xtet->get_fem_noeud(2),xtet->get_fem_noeud(3),uv3,xyz3);
      int inter_face4=intersection_segment_triangle(xseg->get_fem_noeud(0),xseg->get_fem_noeud(1),xtet->get_fem_noeud(1),xtet->get_fem_noeud(2),xtet->get_fem_noeud(3),uv4,xyz4);
      if (inter_face1==1)
        {
        std::map<double,FEM_NOEUD*>::iterator it=lstnoeud.lower_bound(uv1[2]-1e-6);
        double l=0.;
        if (it!=lstnoeud.end())
            l=(uv1[2]-it->first)*(uv1[2]-it->first);
        if ((l>1e-10) || (it==lstnoeud.end()))
                {
                FEM_NOEUD* no=inserer_noeud(ele_inter->are,xyz1[0],xyz1[1],xyz1[2]);
                no->change_etat(0,FRONTIERE);
                std::pair<double,FEM_NOEUD*> tmp(uv1[2],no);
                lstnoeud.insert(tmp);
                }
        }
      if (inter_face2==1)
        {
        std::map<double,FEM_NOEUD*>::iterator it=lstnoeud.lower_bound(uv2[2]-1e-6);
        double l=0.;
        if (it!=lstnoeud.end())
            l=(uv2[2]-it->first)*(uv2[2]-it->first);
        if ((l>1e-10) || (it==lstnoeud.end()))
                {
                FEM_NOEUD* no=inserer_noeud(ele_inter->are,xyz2[0],xyz2[1],xyz2[2]);
                no->change_etat(0,FRONTIERE);
                std::pair<double,FEM_NOEUD*> tmp(uv2[2],no);
                lstnoeud.insert(tmp);
                }
        }  
      if (inter_face1==3)
        {
        std::map<double,FEM_NOEUD*>::iterator it=lstnoeud.lower_bound(uv3[2]-1e-6);
        double l=0.;
        if (it!=lstnoeud.end())
            l=(uv3[2]-it->first)*(uv3[2]-it->first);
        if ((l>1e-10) || (it==lstnoeud.end()))
                {
                FEM_NOEUD* no=inserer_noeud(ele_inter->are,xyz3[0],xyz3[1],xyz3[2]);
                no->change_etat(0,FRONTIERE);
                std::pair<double,FEM_NOEUD*> tmp(uv3[2],no);
                lstnoeud.insert(tmp);
                }
        }
        if (inter_face4==1)
        {
        std::map<double,FEM_NOEUD*>::iterator it=lstnoeud.lower_bound(uv4[2]-1e-6);
        double l=0.;
        if (it!=lstnoeud.end())
            l=(uv4[2]-it->first)*(uv4[2]-it->first);
        if ((l>1e-10) || (it==lstnoeud.end()))
                {
                FEM_NOEUD* no=inserer_noeud(ele_inter->are,xyz4[0],xyz4[1],xyz4[2]);
                no->change_etat(0,FRONTIERE);
                std::pair<double,FEM_NOEUD*> tmp(uv4[2],no);
                lstnoeud.insert(tmp);
                }
        }
      }
      int nbnoeudcree=lstnoeud.size();
      for (std::map<double,FEM_NOEUD*>::iterator it=lstnoeud.begin();it!=lstnoeud.end();it++)
        decoupe_element(-1,it->second,ele_inter->ele);
      std::pair<double,FEM_NOEUD*> tmp(0.,xseg->get_fem_noeud(0));
      lstnoeud.insert(tmp);
      std::pair<double,FEM_NOEUD*> tmp2(1.,xseg->get_fem_noeud(1));
      lstnoeud.insert(tmp2);
      std::map<double,FEM_NOEUD*>::iterator it=lstnoeud.begin();
      std::map<double,FEM_NOEUD*>::iterator it2=lstnoeud.begin();
      it2++;
      for (int i=0;i<nbnoeudcree+1;i++)
              {
              FEM_NOEUD* tab[2];
              tab[0]=it->second;
              tab[1]=it2->second;
              XFEM_SEGMENT2* sxseg=new XFEM_SEGMENT2(ele_inter->ele,xseg->get_lien_topologie(),tab);
              mai->ajouter_xfem_element1(sxseg);
              sxseg->change_etat(XFRONTIERE);
              it++;
              it2++;
             OT_VECTEUR_3D vec(tab[0]->get_coord(),tab[1]->get_coord());
             if (vec.get_longueur()<1e-10)
               std::cout << " toto " << std::endl; 
                
              }
      mai->supprimer_xfem_element1id(xseg->get_id());
      
        
/*if (ele_inter->ele->get_nb_xfem(3)==0)
    {
    decoupe_element(ele_inter->status1,ele_inter->ele,ele_inter->no1,ele_inter->ele->get_fem_noeud(0),ele_inter->ele->get_fem_noeud(1),ele_inter->ele->get_fem_noeud(2),ele_inter->ele->get_fem_noeud(3)); 
    if (ele_inter->ele->get_nb_xfem(3)==0) decoupe_element(ele_inter->status2,ele_inter->ele,ele_inter->no2,ele_inter->ele->get_fem_noeud(0),ele_inter->ele->get_fem_noeud(1),ele_inter->ele->get_fem_noeud(2),ele_inter->ele->get_fem_noeud(3)); 
    for (int i=0;i<ele_inter->ele->get_nb_xfem(3);i++)
      {
      XFEM_ELEMENT3* xele;//=ele_inter->ele->get_xfem(i);
      int status=outilfem.estdansletetra(xele,ele_inter->xyz2[0],ele_inter->xyz2[1],ele_inter->xyz2[2]);  
      if (outilfem.compare_etat(status,FEM_MAILLAGE_OUTILS::INTERIEUR)==1)  
        {
        decoupe_element(ele_inter->status2,ele_inter->ele,ele_inter->no2,xele->get_fem_noeud(0),xele->get_fem_noeud(1),xele->get_fem_noeud(2),xele->get_fem_noeud(3)); 
        mai->supprimer_xfem_element3id(xele->get_id());  
        }
        
      }
      
    }
 else
    {
    int inter;
    double t=0.;
    FEM_NOEUD* no1=ele_inter->no1;
    double xyz1[3]={ele_inter->xyz1[0],ele_inter->xyz1[1],ele_inter->xyz1[2]};
    do 
          {
          inter_ele_arete xele_inter;
          inter=decoupe_segment_xtetra(ele_inter->ele,xseg,xyz1,ele_inter->xyz2,t,xele_inter);
          if (inter==1)
              {
              FEM_NOEUD* no2=new FEM_NOEUD(ele_inter->ele->get_mg_element_maillage()->get_lien_topologie(),xele_inter.xyz2[0],xele_inter.xyz2[1],xele_inter.xyz2[2]);
              mai->ajouter_fem_noeud(no2);
              FEM_NOEUD* tab[2];
              tab[0]=no1;
              tab[1]=no2;
              XFEM_SEGMENT2* xseg;//=new XFEM_SEGMENT2((FEM_SEGMENT2*)seg,tab);
              mai->ajouter_xfem_element1(xseg);
              xele_inter.no1=no1;
              xele_inter.no2=no2;
              //decoupe_element(xseg,&xele_inter);
              xyz1[0]=xele_inter.xyz2[0];
              xyz1[1]=xele_inter.xyz2[1];
              xyz1[2]=xele_inter.xyz2[2];
              no1=no2;
              mai->supprimer_xfem_element3id(xele_inter.xele->get_id());
              }
            
          }
        while (inter!=0);
        FEM_NOEUD* tab[2];
        tab[0]=no1;
        tab[1]=ele_inter->no2;
        //FEM_SEGMENT2* seg=new FEM_SEGMENT2(are,tab);
        //mai->ajouter_fem_element1(seg);
        //mai->supprimer_xfem_element3id(xele->get_id());
    }*/
}
 void  TOIBREP::decoupe_xtri(FEM_ELEMENT3 *tet,int nb_xtri)
 {
 if (tet->get_id()==433021)
   std::cout <<"ttoot";
 if (nb_xtri!=2) return;
 int nbseg=5;
 int nbno=4;
 FEM_NOEUD** tabnoeud=new FEM_NOEUD*[nbno];
 FEM_NOEUD* tab1[3];
 FEM_NOEUD* tab2[3];
 XFEM_ELEMENT2* xtri1=(XFEM_ELEMENT2*)tet->get_xfem(2,tet->get_nb_xfem(2)-1);
 tab1[0]=xtri1->get_fem_noeud(0);
 tab1[1]=xtri1->get_fem_noeud(1);
 tab1[2]=xtri1->get_fem_noeud(2);
 XFEM_ELEMENT2* xtri2=(XFEM_ELEMENT2*)tet->get_xfem(2,tet->get_nb_xfem(2)-2);
  tab2[0]=xtri2->get_fem_noeud(0);
  tab2[1]=xtri2->get_fem_noeud(1);
  tab2[2]=xtri2->get_fem_noeud(2);
  if (tab1[0]!=tab2[0])
    if (tab1[0]!=tab2[1])
      if (tab1[0]!=tab2[2])
            {
            tabnoeud[0]=tab1[0]; 
            tabnoeud[1]=tab1[1]; 
            tabnoeud[2]=tab1[2]; 
            }
  if (tab1[1]!=tab2[0])
    if (tab1[1]!=tab2[1])
      if (tab1[1]!=tab2[2])
            {
            tabnoeud[0]=tab1[1]; 
            tabnoeud[1]=tab1[2]; 
            tabnoeud[2]=tab1[0]; 
            }
  if (tab1[2]!=tab2[0])
    if (tab1[2]!=tab2[1])
      if (tab1[2]!=tab2[2])
            {
            tabnoeud[0]=tab1[2]; 
            tabnoeud[1]=tab1[0]; 
            tabnoeud[2]=tab1[1]; 
            }
  if (tab2[2]!=tab1[0])
    if (tab2[2]!=tab1[1])
      if (tab2[2]!=tab1[2])
          tabnoeud[3]=tab2[2];
  if (tab2[1]!=tab1[0])
    if (tab2[1]!=tab1[1])
      if (tab2[1]!=tab1[2])
          tabnoeud[3]=tab2[1];
if (tab2[0]!=tab1[0])
    if (tab2[0]!=tab1[1])
      if (tab2[0]!=tab1[2])
          tabnoeud[3]=tab2[0];
 TPL_LISTE_ENTITE<XFEM_ELEMENT1*> listeasupprimer;
 int nbxfem1=tet->get_nb_xfem(1);
 for (int i=0;i<nbxfem1;i++)
    {
    XFEM_ELEMENT1* xfem1=(XFEM_ELEMENT1*)tet->get_xfem(1,i);
    FEM_NOEUD* no1=xfem1->get_fem_noeud(0);
    FEM_NOEUD* no2=xfem1->get_fem_noeud(1);
    if ((no1->get_id()==598025) && (no2->get_id()==598034))
            std::cout <<"ttoot";
    if ((no2->get_id()==598025) && (no1->get_id()==598034))
            std::cout <<"ttoot";
    double xyz[3];
    int inter=inter_segment_segment_plan(xtri1,tabnoeud[1],tabnoeud[2],no1,no2,xyz);
    if (xtri1->get_lien_topologie()->get_id()==789)
     if (tet->get_nb_xfem(1)==1) 
      //if (tet->get_nb_xfem(2)==2) 
           if (inter==0)
             std::cout << tet->get_id();
    if (inter==1)
            {
            FEM_NOEUD* no=inserer_noeud(xfem1->get_lien_topologie(),xyz[0],xyz[1],xyz[2]);
            no->change_etat(0,INCONNU);
            FEM_NOEUD* tab[2];
            tab[0]=no1;
            tab[1]=no;
            XFEM_SEGMENT2* sxseg1=new XFEM_SEGMENT2(tet,xfem1->get_lien_topologie(),tab);
            mai->ajouter_xfem_element1(sxseg1);
            sxseg1->change_etat(XFRONTIERE);
            tab[0]=no;
            tab[1]=no2;
            XFEM_SEGMENT2* sxseg2=new XFEM_SEGMENT2(tet,xfem1->get_lien_topologie(),tab);
            mai->ajouter_xfem_element1(sxseg2);
            sxseg2->change_etat(XFRONTIERE);
            listeasupprimer.ajouter(xfem1);
            }
    }
 for (int i=0;i<listeasupprimer.get_nb();i++)
   mai->supprimer_xfem_element1id(listeasupprimer.get(i)->get_id());
 delete [] tabnoeud;
 }
  
void TOIBREP::recherche_interieur_face(FEM_ELEMENT3* tet,MG_FACE* face,TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lstentitetopo)
{
int nb=tet->get_nb_xfem(1);
int traite=0;
for (int i=0;i<nb;i++)
    {
    XFEM_ELEMENT1* xfem=(XFEM_ELEMENT1*)tet->get_xfem(1,i);
    if (lstentitetopo->existe(xfem->get_lien_topologie())==1)
            {
            traite=1;
            MG_ARETE* are=(MG_ARETE*)xfem->get_lien_topologie();
            MG_COARETE* coare;
            int icor=0;
            do
                  {
                    coare=are->get_mg_coarete(icor);
                    icor++;
                  }
            while (coare->get_boucle()->get_mg_face()!=face);
            double ori=coare->get_orientation();
            double *xyz1=xfem->get_fem_noeud(0)->get_coord();
            double *xyz2=xfem->get_fem_noeud(1)->get_coord();
            double t;
            are->inverser(t,xyz1);
            double dxyz[3];
            are->deriver(t,dxyz);
            OT_VECTEUR_3D dir(dxyz);
            dir=ori*dir;
            dir.norme();
            OT_VECTEUR_3D dircor(xyz1,xyz2);
            dircor.norme();
            if (dir*dircor<0) dircor=(-1.)*dircor;
            double uv[2];
            face->inverser(uv,xyz1);
            double vecnormal[3];
            face->calcul_normale_unitaire(uv,vecnormal);
            OT_VECTEUR_3D nor(vecnormal);
            nor=face->get_mg_coface(0)->get_orientation()*nor;
            OT_VECTEUR_3D vecext=nor&dircor;
            int nbxtri=tet->get_nb_xfem(2);
            for (int j=0;j<nbxtri;j++)
                  {
                  XFEM_ELEMENT2* xfem2=(XFEM_ELEMENT2*)tet->get_xfem(2,j);
                  if (xfem2->get_lien_topologie()==face) 
                  if ((xfem->get_fem_noeud(0)==xfem2->get_fem_noeud(0)) || (xfem->get_fem_noeud(0)==xfem2->get_fem_noeud(1)) || (xfem->get_fem_noeud(0)==xfem2->get_fem_noeud(2)) )
                  if ((xfem->get_fem_noeud(1)==xfem2->get_fem_noeud(0)) || (xfem->get_fem_noeud(1)==xfem2->get_fem_noeud(1)) || (xfem->get_fem_noeud(1)==xfem2->get_fem_noeud(2)) )
                        {
                        double x=xfem2->get_fem_noeud(0)->get_x()+xfem2->get_fem_noeud(1)->get_x()+xfem2->get_fem_noeud(2)->get_x();  
                        double y=xfem2->get_fem_noeud(0)->get_y()+xfem2->get_fem_noeud(1)->get_y()+xfem2->get_fem_noeud(2)->get_y();  
                        double z=xfem2->get_fem_noeud(0)->get_z()+xfem2->get_fem_noeud(1)->get_z()+xfem2->get_fem_noeud(2)->get_z();
                        OT_VECTEUR_3D dirxtri(x/3.-xyz1[0],y/3.-xyz1[1],z/3.-xyz1[2]);
                        dirxtri.norme();
                        if (dirxtri*vecext>0.) 
                            {
                              xfem2->change_etat(XINTERIEUR); 
                              if (xfem2->get_fem_noeud(0)->get_etat(0)==INCONNU) xfem2->get_fem_noeud(0)->change_etat(0,FRONTIERE);
                              if (xfem2->get_fem_noeud(1)->get_etat(0)==INCONNU) xfem2->get_fem_noeud(1)->change_etat(0,FRONTIERE);
                              if (xfem2->get_fem_noeud(2)->get_etat(0)==INCONNU) xfem2->get_fem_noeud(2)->change_etat(0,FRONTIERE);
                            }
                            else 
                            {
                            xfem2->change_etat(XEXTERIEUR);
                            if (xfem2->get_fem_noeud(0)->get_etat(0)==INCONNU) xfem2->get_fem_noeud(0)->change_etat(0,EXTERIEUR);
                            if (xfem2->get_fem_noeud(1)->get_etat(0)==INCONNU) xfem2->get_fem_noeud(1)->change_etat(0,EXTERIEUR);
                            if (xfem2->get_fem_noeud(2)->get_etat(0)==INCONNU) xfem2->get_fem_noeud(2)->change_etat(0,EXTERIEUR);
                             }
                        }
                  }
              
            }
    }/*
if (traite>0)
{
int changement=0;
do
{
int nb2=tet->get_nb_xfem(2);
changement=0;
int casatraiter=0;
for (int i=0;i<nb2;i++)
    {
    XFEM_ELEMENT2* xfem=(XFEM_ELEMENT2*)tet->get_xfem(2,i);
    if (xfem->get_lien_topologie()==face)
    if (xfem->get_etat()==XINCONNU)
          {
          casatraiter=1;
          if ((xfem->get_fem_noeud(0)->get_etat(0)==EXTERIEUR) || (xfem->get_fem_noeud(1)->get_etat(0)==EXTERIEUR) || (xfem->get_fem_noeud(2)->get_etat(0)==EXTERIEUR))  
              {
              xfem->change_etat(XEXTERIEUR);
              if (xfem->get_fem_noeud(0)->get_etat(0)==INCONNU) xfem->get_fem_noeud(0)->change_etat(0,EXTERIEUR);
              if (xfem->get_fem_noeud(1)->get_etat(0)==INCONNU) xfem->get_fem_noeud(1)->change_etat(0,EXTERIEUR);
              if (xfem->get_fem_noeud(2)->get_etat(0)==INCONNU) xfem->get_fem_noeud(2)->change_etat(0,EXTERIEUR);
              changement=1;
              }
          else
              if (((xfem->get_fem_noeud(0)->get_etat(0)==FRONTIERE) && (xfem->get_fem_noeud(0)->get_lien_topologie()->get_dimension()==2)) || ((xfem->get_fem_noeud(1)->get_etat(0)==FRONTIERE) && (xfem->get_fem_noeud(1)->get_lien_topologie()->get_dimension()==2)) || ((xfem->get_fem_noeud(2)->get_etat(0)==FRONTIERE) && (xfem->get_fem_noeud(2)->get_lien_topologie()->get_dimension()==2)))  
              {
              xfem->change_etat(XINTERIEUR);
              if (xfem->get_fem_noeud(0)->get_etat(0)==INCONNU) xfem->get_fem_noeud(0)->change_etat(0,FRONTIERE);
              if (xfem->get_fem_noeud(1)->get_etat(0)==INCONNU) xfem->get_fem_noeud(1)->change_etat(0,FRONTIERE);
              if (xfem->get_fem_noeud(2)->get_etat(0)==INCONNU) xfem->get_fem_noeud(2)->change_etat(0,FRONTIERE);
              changement=1;
              }
          }
    }
    if (casatraiter==0) changement=1;
}
while (changement==0);
}*/
}            
XFEM_TRIANGLE3* TOIBREP::inserer_xtriangle(FEM_ELEMENT_MAILLAGE* tet,MG_FACE* face,FEM_NOEUD** tabnoeud)
{
      OT_VECTEUR_3D vec1(tabnoeud[0]->get_coord(),tabnoeud[1]->get_coord());
      OT_VECTEUR_3D vec2(tabnoeud[0]->get_coord(),tabnoeud[2]->get_coord());
      OT_VECTEUR_3D vec3=vec1&vec2;
      double eps=0.0000000001*longueur_caracteristique*longueur_caracteristique;
      double aire=vec3.get_longueur()*0.5;
      if (aire<eps) 
        return NULL;
      XFEM_TRIANGLE3* xtri=new XFEM_TRIANGLE3(tet,face,tabnoeud);
      oriente_tri(xtri,face);
      mai->ajouter_xfem_element2(xtri);  
      if (xtri->get_id()==616861)
            std::cout<<""<<std::endl;
      if (xtri->get_id()==616862)
            std::cout<<""<<std::endl;
      if (xtri->get_id()==616855)
            std::cout<<""<<std::endl;
      return xtri;
}


FEM_NOEUD* TOIBREP::inserer_noeud(MG_ELEMENT_TOPOLOGIQUE* topo,double x,double y, double z)
{
TPL_MAP_ENTITE<FEM_NOEUD*> lstnoeud;
octree_noeud.rechercher(x,y,z,0.000000001,lstnoeud);
TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR it;
for (FEM_NOEUD* no=lstnoeud.get_premier(it);no!=NULL;no=lstnoeud.get_suivant(it))
      {
      double *xyz=no->get_coord();
      OT_VECTEUR_3D vec(xyz[0]-x,xyz[1]-y,xyz[2]-z);
      if (vec.get_longueur()<0.000000001*longueur_caracteristique)
            {
            if (no->get_lien_topologie()==NULL)
              no->change_topologie_null(topo);
            return no;
            }
      }
FEM_NOEUD* no=new FEM_NOEUD(topo,x,y,z);
mai->ajouter_fem_noeud(no);
octree_noeud.inserer(no);
return no;
}
  
  
void TOIBREP::importer_et_decouper(MG_GROUPE_TOPOLOGIQUE* mggt)
{
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;
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];
        lst_noeud.ajouter(noeud);
        }
BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
boite.change_grosseur(1.1);
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base recherche boite");
longueur_caracteristique=boite.get_rayon()*2.;
octree_tetra.initialiser(&lst_noeud,1,boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax());
octree_noeud.initialiser(&octree_tetra);
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base ini octree");
LISTE_FEM_ELEMENT3::iterator it3;
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
  {
  octree_tetra.inserer(tet);
  tet->change_etat(0,INCONNU);
  }
for (FEM_NOEUD* noeud=mai->get_premier_noeud(it2);noeud!=NULL;noeud=mai->get_suivant_noeud(it2))
  octree_noeud.inserer(noeud);
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);
                }
        }
        
TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*>::ITERATEUR it;
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"      Predecoupage sommets");
          affiche(mess);
          }    
for (MG_ELEMENT_TOPOLOGIQUE* eletopo=lst.get_premier(it);eletopo!=NULL;eletopo=lst.get_suivant(it))
  if (eletopo->get_dimension()==0)
      {
        MG_SOMMET* som=(MG_SOMMET*)eletopo;
        double xyz[3];
        som->get_point()->evaluer(xyz);
        TPL_MAP_ENTITE<FEM_ELEMENT3*> liste;
        octree_tetra.rechercher(xyz[0],xyz[1],xyz[2],1e-1,liste);
        TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR itl;
        for (FEM_ELEMENT3* ele=liste.get_premier(itl);ele!=NULL;ele=liste.get_suivant(itl))
            {
            int status=outilfem.estdansletetra(ele,xyz[0],xyz[1],xyz[2]);  
            if (outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::INTERIEUR)==1)  
              { 
              FEM_NOEUD* no;
              if (som->get_lien_fem_maillage()->get_nb()!=0) no=(FEM_NOEUD*)som->get_lien_fem_maillage()->get(0);
              else 
                {
                no=inserer_noeud(eletopo,xyz[0],xyz[1],xyz[2]);
                no->change_etat(0,FRONTIERE);
                }
              XFEM_ELEMENT0 *xele0=new XFEM_ELEMENT0(ele,eletopo,&no);
              mai->ajouter_xfem_element0(xele0);
              xele0->change_etat(XFRONTIERE);
              decoupe_element(status,no,ele);
              }
        
        
            }
        
      }
if (compteur!=NULL) 
          {
          double tps=compteur->ajouter_etape("Predecoupage sommets");
          if (affichageactif)
            { 
            char mess[300];
            sprintf(mess,"           ---> CPU : %lf",tps);
            affiche(mess);
            }    
          }
          testtopo();
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"      Predecoupage aretes");
          affiche(mess);
          }
for (MG_ELEMENT_TOPOLOGIQUE* eletopo=lst.get_premier(it);eletopo!=NULL;eletopo=lst.get_suivant(it))
  if (eletopo->get_dimension()==1)
      {
        MG_ARETE* are=(MG_ARETE*)eletopo;
        MG_SOMMET* som1=are->get_cosommet1()->get_sommet();
        MG_SOMMET* som2=are->get_cosommet2()->get_sommet();
        double xyz1[3],xyzcible[3];
        inter_ele_arete ele_inter;
        ele_inter.are=are;
        som1->get_point()->evaluer(xyz1);
        som2->get_point()->evaluer(xyzcible);
        double t=are->get_tmin();
        FEM_NOEUD* no1=(FEM_NOEUD*)som1->get_lien_fem_maillage()->get(0);
        FEM_NOEUD* nofin=(FEM_NOEUD*)som2->get_lien_fem_maillage()->get(0);
        int inter;
        do 
          {
          inter=decoupe_arete_tetra(are,xyz1,xyzcible,t,ele_inter);
          if (inter==1)
              {
              FEM_NOEUD* no2=inserer_noeud(are,ele_inter.xyz2[0],ele_inter.xyz2[1],ele_inter.xyz2[2]);
              no2->change_etat(0,FRONTIERE);
              FEM_NOEUD* tab[2];
              tab[0]=no1;
              tab[1]=no2;
              XFEM_SEGMENT2* xseg=new XFEM_SEGMENT2(ele_inter.ele,are,tab);
              mai->ajouter_xfem_element1(xseg);
              xseg->change_etat(XFRONTIERE);
              decoupe_element(ele_inter.status1,no1,ele_inter.ele);
              decoupe_element(ele_inter.status2,no2,ele_inter.ele);
              ele_inter.no1=no1;
              ele_inter.no2=no2;
              decoupe_element(xseg,&ele_inter);
              xyz1[0]=ele_inter.xyz2[0];
              xyz1[1]=ele_inter.xyz2[1];
              xyz1[2]=ele_inter.xyz2[2];
              no1=no2;
              }
            
          }
        while (inter!=0);
        FEM_NOEUD* tab[2];
        tab[0]=no1;
        tab[1]=nofin;
        XFEM_SEGMENT2* xseg=new XFEM_SEGMENT2(ele_inter.ele,are,tab);
         xseg->change_etat(XFRONTIERE);
         mai->ajouter_xfem_element1(xseg);
        ele_inter.no1=no1;
        ele_inter.no2=nofin;
        decoupe_element(-1,no1,ele_inter.ele);
        decoupe_element(-1,nofin,ele_inter.ele);
        decoupe_element(xseg,&ele_inter);
                    
      }
 if (compteur!=NULL) 
          {
          double tps=compteur->ajouter_etape("Predecoupage aretes");
          if (affichageactif)
            { 
            char mess[300];
            sprintf(mess,"           ---> CPU : %lf",tps);
            affiche(mess);
            }    
          }
testtopo();
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"      Decoupage XFEM");
          affiche(mess);
          }  
int numdebuttps;
if (compteur!=NULL) 
          {
          std::string nom;
          numdebuttps=compteur->get_nb_etape()-1;
          }
int compteur_dim2=0;
for (MG_ELEMENT_TOPOLOGIQUE* eletopo=lst.get_premier(it);eletopo!=NULL;eletopo=lst.get_suivant(it))
  if ((eletopo->get_dimension()==2) && (compteur_dim2++==2))
      {
      for (FEM_NOEUD* noeud=mai->get_premier_noeud(it2);noeud!=NULL;noeud=mai->get_suivant_noeud(it2))
            noeud->change_solution(1e300);
      for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
            tet->change_solution(0.);
        MG_FACE* face=(MG_FACE*)eletopo;
      if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"        Face id %lu",face->get_id());
          affiche(mess);
          }
        //if (!lst.existe(face)) continue;
        TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lstentitetopo;
        face->get_topologie_sousjacente(&lstentitetopo);
        levelsetn(&lst,face,&lstentitetopo);
        for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
            if (tet->get_etat(1)==1)
            {
            int nb_xtri=decoupe_tetra(tet,face);
            decoupe_xtri(tet,nb_xtri);
            decoupe_element2_par_element1et0(tet,&lstentitetopo);
            //recherche_interieur_face(tet,face,&lstentitetopo);
            }
            else
            {
            for (int num=0;num<2;num++)
              for (int i=0;i<tet->get_nb_xfem(num);i++)
                {
                MG_ELEMENT_TOPOLOGIQUE* mgt=tet->get_xfem(num,i)->get_lien_topologie();
                if (lstentitetopo.existe(mgt)==1)
                    std::cout << "probleme tetra "<< num << ":" << tet->get_id() << " " << tet->get_etat(1) <<" = " << tet->get_fem_noeud(0)->get_etat(1)<< "-"<< tet->get_fem_noeud(1)->get_etat(1)<< "-"<< tet->get_fem_noeud(2)->get_etat(1)<< "-"<< tet->get_fem_noeud(3)->get_etat(1)  << std::endl;
                
                }
            }
        
      }
      
if (compteur!=NULL) 
          {
          std::string nom;
          double tps;
          compteur->get_etape(numdebuttps,compteur->get_nb_etape()-1,nom,tps);
          if (affichageactif)
            { 
            char mess[300];
            sprintf(mess,"           ---> CPU : %lf",tps);
            affiche(mess);
            }    
            
          } 
 
 
for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
  {
  double vol;
  int res=tet->verifie_validite_decoupage_xfem(&vol);
  int valide=res&1;
  int volnul=(res<<1)&1;
  if (valide==0) 
        {
        std::cout << "Erreur decoupage tetra de id " << tet->get_id() << " volume residuel " << vol << std::endl;   
        }
  if (volnul==1) 
        {
        std::cout << "Erreur decoupage tetra de id " << tet->get_id() << " volume nul d'un xfem" << std::endl;   
        }
  int nb1=tet->get_nb_xfem(1);
  for (int i=0;i<nb1;i++)
    {
    XFEM_ELEMENT1* xseg=(XFEM_ELEMENT1*)tet->get_xfem(1,i);
    FEM_NOEUD* no1=xseg->get_fem_noeud(0);
    FEM_NOEUD* no2=xseg->get_fem_noeud(1);
    int nb3=tet->get_nb_xfem(3);
    int trouve=0;
    for (int j=0;j<nb3;j++)
      {
      XFEM_ELEMENT3* xtettmp=(XFEM_ELEMENT3*)tet->get_xfem(3,j);
      FEM_NOEUD* not1=xtettmp->get_fem_noeud(0);
      FEM_NOEUD* not2=xtettmp->get_fem_noeud(1);
      FEM_NOEUD* not3=xtettmp->get_fem_noeud(2);
      FEM_NOEUD* not4=xtettmp->get_fem_noeud(3);
      if ((no1==not1) || (no1==not2) ||  (no1==not3) || (no1==not4))
      if ((no2==not1) || (no2==not2) ||  (no2==not3) || (no2==not4))
              trouve=1;
        
      }
    if (trouve==0) 
      std::cout << "Erreur decoupage tetra de id " << tet->get_id() << " xfem1 pas coherent " << std::endl;
    }
    if (tet->get_id()==-407234)
            {
            MG_MAILLAGE* mgmai=mai->get_mg_maillage();
            FEM_MAILLAGE *femtmp=new FEM_MAILLAGE(mgmai->get_mg_geometrie(),mgmai,1);  
            mgmai->get_gestionnaire()->ajouter_fem_maillage(femtmp);
            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);
            FEM_NOEUD* notmp1=new FEM_NOEUD(no1->get_lien_topologie(),no1->get_x(),no1->get_y(),no1->get_z());
            FEM_NOEUD* notmp2=new FEM_NOEUD(no2->get_lien_topologie(),no2->get_x(),no2->get_y(),no2->get_z());
            FEM_NOEUD* notmp3=new FEM_NOEUD(no3->get_lien_topologie(),no3->get_x(),no3->get_y(),no3->get_z());
            FEM_NOEUD* notmp4=new FEM_NOEUD(no4->get_lien_topologie(),no4->get_x(),no4->get_y(),no4->get_z());
            femtmp->ajouter_fem_noeud(notmp1);
            femtmp->ajouter_fem_noeud(notmp2);
            femtmp->ajouter_fem_noeud(notmp3);
            femtmp->ajouter_fem_noeud(notmp4);
            FEM_NOEUD* tab[4];
            tab[0]=notmp1;
            tab[1]=notmp2;
            tab[2]=notmp3;
            tab[3]=notmp4;
            FEM_TETRA4* tettmp=new FEM_TETRA4(tet->get_lien_topologie(),tet->get_mg_element_maillage(),tab);
            femtmp->ajouter_fem_element3(tettmp);
            int nb1=tet->get_nb_xfem(1);
            for (int j=0;j<nb1;j++)
                        {
                        XFEM_ELEMENT1* xtettmp=(XFEM_ELEMENT1*)tet->get_xfem(1,j);
                        FEM_NOEUD* not1=xtettmp->get_fem_noeud(0);
                        FEM_NOEUD* not2=xtettmp->get_fem_noeud(1);
                        FEM_NOEUD* notmp1=new FEM_NOEUD(not1->get_lien_topologie(),not1->get_x(),not1->get_y(),not1->get_z());
                        FEM_NOEUD* notmp2=new FEM_NOEUD(not2->get_lien_topologie(),not2->get_x(),not2->get_y(),not2->get_z());
                        femtmp->ajouter_fem_noeud(notmp1);
                        femtmp->ajouter_fem_noeud(notmp2);
                        FEM_NOEUD* tab[4];
                        tab[0]=notmp1;
                        tab[1]=notmp2;
                        XFEM_SEGMENT2* xtet2=new XFEM_SEGMENT2(tettmp,tet->get_lien_topologie(),tab);
                        femtmp->ajouter_xfem_element1(xtet2);
                        OT_VECTEUR_3D vec(notmp1->get_coord(),notmp2->get_coord());
                        std::cout << vec.get_longueur() << std::endl;
                        }
            int nb3=tet->get_nb_xfem(3);
            for (int j=0;j<nb3;j++)
                        {
                        XFEM_ELEMENT3* xtettmp=(XFEM_ELEMENT3*)tet->get_xfem(3,j);
                        FEM_NOEUD* not1=xtettmp->get_fem_noeud(0);
                        FEM_NOEUD* not2=xtettmp->get_fem_noeud(1);
                        FEM_NOEUD* not3=xtettmp->get_fem_noeud(2);
                        FEM_NOEUD* not4=xtettmp->get_fem_noeud(3);
                        FEM_NOEUD* notmp1=new FEM_NOEUD(not1->get_lien_topologie(),not1->get_x(),not1->get_y(),not1->get_z());
                        FEM_NOEUD* notmp2=new FEM_NOEUD(not2->get_lien_topologie(),not2->get_x(),not2->get_y(),not2->get_z());
                        FEM_NOEUD* notmp3=new FEM_NOEUD(not3->get_lien_topologie(),not3->get_x(),not3->get_y(),not3->get_z());
                        FEM_NOEUD* notmp4=new FEM_NOEUD(not4->get_lien_topologie(),not4->get_x(),not4->get_y(),not4->get_z());
                        femtmp->ajouter_fem_noeud(notmp1);
                        femtmp->ajouter_fem_noeud(notmp2);
                        femtmp->ajouter_fem_noeud(notmp3);
                        femtmp->ajouter_fem_noeud(notmp4);
                        FEM_NOEUD* tab[4];
                        tab[0]=notmp1;
                        tab[1]=notmp2;
                        tab[2]=notmp3;
                        tab[3]=notmp4;
                        XFEM_TETRA4* xtet2=new XFEM_TETRA4(tettmp,tet->get_lien_topologie(),tab);
                        femtmp->ajouter_xfem_element3(xtet2);
                        }
            
      
    }
    
  }
 

}


int TOIBREP::decoupe_segment_xtetra(FEM_ELEMENT3* ele,XFEM_ELEMENT1* seg,double *xyz,double *xyzcible,double &t,inter_ele_arete& xele_inter)
{
int code_retour=0;
int nb=ele->get_nb_xfem(3);
for (int i=0;i<nb;i++)
      {
      XFEM_ELEMENT3* xele=(XFEM_ELEMENT3*)ele->get_xfem(3,i);
      int status=outilfem.estdansletetra(xele,xyz[0],xyz[1],xyz[2]);
      if (outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::INTERIEUR)==1) 
              {
              int status2=outilfem.estdansletetra(xele,xyzcible[0],xyzcible[1],xyzcible[2]);  
              if (outilfem.compare_etat_tetra(status2,FEM_MAILLAGE_OUTILS::INTERIEUR)==1)  return 0;
              double xyz1[3],uv1[3],xyz4[3],uv4[3],xyz2[3],uv2[3],xyz3[3],uv3[3];
              uv1[2]=t;
              uv2[2]=t;
              uv3[2]=t;
              uv4[2]=t;
              int inter_face1=intersection_segment_triangle(seg->get_fem_noeud(0),seg->get_fem_noeud(1),xele->get_fem_noeud(0),xele->get_fem_noeud(1),xele->get_fem_noeud(2),xyz1,uv1);
              int inter_face2=intersection_segment_triangle(seg->get_fem_noeud(0),seg->get_fem_noeud(1),xele->get_fem_noeud(0),xele->get_fem_noeud(1),xele->get_fem_noeud(3),xyz2,uv2);
              int inter_face3=intersection_segment_triangle(seg->get_fem_noeud(0),seg->get_fem_noeud(1),xele->get_fem_noeud(0),xele->get_fem_noeud(2),xele->get_fem_noeud(3),xyz3,uv3);
              int inter_face4=intersection_segment_triangle(seg->get_fem_noeud(0),seg->get_fem_noeud(1),xele->get_fem_noeud(1),xele->get_fem_noeud(2),xele->get_fem_noeud(3),xyz4,uv4);
              xele_inter.status2=1;
              if (inter_face1==1)
                  {
                  t=uv1[2];
                  xele_inter.xele=xele;
                  xele_inter.xyz1[0]=xyz[0];
                  xele_inter.xyz1[1]=xyz[1];
                  xele_inter.xyz1[2]=xyz[2];
                  xele_inter.xyz2[0]=xyz1[0];
                  xele_inter.xyz2[1]=xyz1[1];
                  xele_inter.xyz2[2]=xyz1[2];
                  xele_inter.status1=status;
                  xele_inter.status2=xele_inter.status2+4;
                  code_retour=1;
                  }
              if (inter_face2==1)
                  {
                  t=uv2[1];
                  xele_inter.xele=xele;
                  xele_inter.xyz1[0]=xyz[0];
                  xele_inter.xyz1[1]=xyz[1];
                  xele_inter.xyz1[2]=xyz[2];
                  xele_inter.xyz2[0]=xyz2[0];
                  xele_inter.xyz2[1]=xyz2[1];
                  xele_inter.xyz2[2]=xyz2[2];
                  xele_inter.status1=status;
                  xele_inter.status2=xele_inter.status2+8;
                  code_retour=1;
                  }
                if (inter_face3==1)
                  {
                  t=uv3[2];
                  xele_inter.xele=xele;
                  xele_inter.xyz1[0]=xyz[0];
                  xele_inter.xyz1[1]=xyz[1];
                  xele_inter.xyz1[2]=xyz[2];
                  xele_inter.xyz2[0]=xyz3[0];
                  xele_inter.xyz2[1]=xyz3[1];
                  xele_inter.xyz2[2]=xyz3[2];
                  xele_inter.status1=status;
                  xele_inter.status2=xele_inter.status2+16;
                  code_retour=1;
                  }
              if (inter_face4==1)
                  {
                  t=uv4[2];
                  xele_inter.xele=xele;
                  xele_inter.xyz1[0]=xyz[0];
                  xele_inter.xyz1[1]=xyz[1];
                  xele_inter.xyz1[2]=xyz[2];
                  xele_inter.xyz2[0]=xyz4[0];
                  xele_inter.xyz2[1]=xyz4[1];
                  xele_inter.xyz2[2]=xyz4[2];
                  xele_inter.status1=status;
                  xele_inter.status2=xele_inter.status2+32;
                  code_retour=1;
                  }
                
              }
      }
}


int TOIBREP::decoupe_arete_tetra(MG_ARETE* are,double *xyz,double *xyzcible,double &t,inter_ele_arete& ele_inter)
{
int code_retour=0;
TPL_MAP_ENTITE<FEM_ELEMENT3*> liste;
octree_tetra.rechercher(xyz[0],xyz[1],xyz[2],1e-10,liste);
TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR itl;
for (FEM_ELEMENT3* ele=liste.get_premier(itl);ele!=NULL;ele=liste.get_suivant(itl))
    {
    int status=outilfem.estdansletetra(ele,xyz[0],xyz[1],xyz[2]);  
    if (outilfem.compare_etat_tetra(status,FEM_MAILLAGE_OUTILS::INTERIEUR)==1)  
          {
          int status2=outilfem.estdansletetra(ele,xyzcible[0],xyzcible[1],xyzcible[2]);  
          if (outilfem.compare_etat_tetra(status2,FEM_MAILLAGE_OUTILS::INTERIEUR)==1)  {ele_inter.ele=ele;return 0;}
          double uv1[3],xyz1[3],uv2[3],xyz2[3],uv3[3],xyz3[3],uv4[3],xyz4[3];
          int inter_face1,inter_face2,inter_face3,inter_face4;
          int boucle=0;
          double titer=are->get_tmax()-t;
          do
          {
          if (t+boucle*0.05*titer>are->get_tmax()) break;
          uv1[2]=t+boucle*0.05*titer;
          uv2[2]=t+boucle*0.05*titer;
          uv3[2]=t+boucle*0.05*titer;
          uv4[2]=t+boucle*0.05*titer;
          inter_face1=intersection_arete_triangle(are,ele->get_fem_noeud(0),ele->get_fem_noeud(1),ele->get_fem_noeud(2),uv1,xyz1);
          inter_face2=intersection_arete_triangle(are,ele->get_fem_noeud(0),ele->get_fem_noeud(1),ele->get_fem_noeud(3),uv2,xyz2);
          inter_face3=intersection_arete_triangle(are,ele->get_fem_noeud(0),ele->get_fem_noeud(2),ele->get_fem_noeud(3),uv3,xyz3);
          inter_face4=intersection_arete_triangle(are,ele->get_fem_noeud(1),ele->get_fem_noeud(2),ele->get_fem_noeud(3),uv4,xyz4);
          boucle++;
          }
          while (inter_face1+inter_face2+inter_face3+inter_face4==0);
          ele_inter.status2=1;
          if (inter_face1==1)
              {
              OT_VECTEUR_3D vec(xyz,xyz1);
              if (vec.get_longueur()<1e-10) inter_face1=0;
              else 
                  {
                   t=uv1[2];
                   ele_inter.ele=ele;
                   ele_inter.status2=ele_inter.status2+4;
                   ele_inter.status1=status;
                   ele_inter.xyz1[0]=xyz[0];
                   ele_inter.xyz1[1]=xyz[1];
                   ele_inter.xyz1[2]=xyz[2];
                   ele_inter.xyz2[0]=xyz1[0];
                   ele_inter.xyz2[1]=xyz1[1];
                   ele_inter.xyz2[2]=xyz1[2];
                   code_retour=1;
                   
                  }
              }
          if (inter_face2==1)
              {
              OT_VECTEUR_3D vec(xyz,xyz2);
              if (vec.get_longueur()<1e-10) inter_face1=0;
               else 
                  {
                   t=uv2[2];
                   ele_inter.ele=ele;
                   ele_inter.status2=ele_inter.status2+8;
                   ele_inter.status1=status;
                   ele_inter.xyz1[0]=xyz[0];
                   ele_inter.xyz1[1]=xyz[1];
                   ele_inter.xyz1[2]=xyz[2];
                   ele_inter.xyz2[0]=xyz2[0];
                   ele_inter.xyz2[1]=xyz2[1];
                   ele_inter.xyz2[2]=xyz2[2];
                   code_retour=1; 
                  }
                
              }
          if (inter_face3==1)
              {
              OT_VECTEUR_3D vec(xyz,xyz3);
              if (vec.get_longueur()<1e-10) inter_face1=0;
               else 
                  {
                   t=uv3[2];
                   ele_inter.ele=ele;
                   ele_inter.status2=ele_inter.status2+16;
                   ele_inter.status1=status;
                   ele_inter.xyz1[0]=xyz[0];
                   ele_inter.xyz1[1]=xyz[1];
                   ele_inter.xyz1[2]=xyz[2];
                   ele_inter.xyz2[0]=xyz3[0];
                   ele_inter.xyz2[1]=xyz3[1];
                   ele_inter.xyz2[2]=xyz3[2];
                   code_retour=1; 
                  }
                
              }
          if (inter_face4==1)
              {
              OT_VECTEUR_3D vec(xyz,xyz4);
              if (vec.get_longueur()<1e-10) inter_face1=0;
               else 
                  {
                   t=uv4[2]; 
                   ele_inter.ele=ele;
                   ele_inter.status2=ele_inter.status2+32;
                   ele_inter.status1=status;
                   ele_inter.xyz1[0]=xyz[0];
                   ele_inter.xyz1[1]=xyz[1];
                   ele_inter.xyz1[2]=xyz[2];
                   ele_inter.xyz2[0]=xyz4[0];
                   ele_inter.xyz2[1]=xyz4[1];
                   ele_inter.xyz2[2]=xyz4[2];
                   code_retour=1;
                   
                  }
                
              }    
          }
  }
return code_retour;
}


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;
}


int TOIBREP::inter_segment_segment_plan(XFEM_ELEMENT2 *xfem2,FEM_NOEUD *no1,FEM_NOEUD *no2,FEM_NOEUD *nos1,FEM_NOEUD *nos2,double *xyz)
{
OT_VECTEUR_3D n1n2(xfem2->get_fem_noeud(0)->get_coord(),xfem2->get_fem_noeud(1)->get_coord());
OT_VECTEUR_3D n1n3(xfem2->get_fem_noeud(0)->get_coord(),xfem2->get_fem_noeud(2)->get_coord());
OT_VECTEUR_3D n=n1n2&n1n3;
n.norme();
double d=-(n.get_x()*xfem2->get_fem_noeud(0)->get_x()+n.get_y()*xfem2->get_fem_noeud(0)->get_y()+n.get_z()*xfem2->get_fem_noeud(0)->get_z());
double t1=-(n.get_x()*nos1->get_x()+n.get_y()*nos1->get_y()+n.get_z()*nos1->get_z()+d);
double t2=-(n.get_x()*nos2->get_x()+n.get_y()*nos2->get_y()+n.get_z()*nos2->get_z()+d);
double xyzs1p[3];
xyzs1p[0]=nos1->get_x()+t2*n.get_x();
xyzs1p[1]=nos1->get_y()+t2*n.get_y();
xyzs1p[2]=nos1->get_z()+t2*n.get_z();
double xyzs2p[3];
xyzs2p[0]=nos2->get_x()+t2*n.get_x();
xyzs2p[1]=nos2->get_y()+t2*n.get_y();
xyzs2p[2]=nos2->get_z()+t2*n.get_z();
double *xyz1=no1->get_coord();
double *xyz2=no2->get_coord();
double ab[3],cd[3],ac[3];
ab[0]=xyz2[0]-xyz1[0];
ab[1]=xyz2[1]-xyz1[1];
ab[2]=xyz2[2]-xyz1[2];
cd[0]=xyzs2p[0]-xyzs1p[0];
cd[1]=xyzs2p[1]-xyzs1p[1];
cd[2]=xyzs2p[2]-xyzs1p[2];
ac[0]=xyzs1p[0]-xyz1[0];
ac[1]=xyzs1p[1]-xyz1[1];
ac[2]=xyzs1p[2]-xyz1[2];
int n1,n2;
double det1=-ab[0]*cd[1]+ab[1]*cd[0];
double det2=-ab[0]*cd[2]+ab[2]*cd[0];
double det3=-ab[1]*cd[2]+ab[2]*cd[1];
double det;
if ( (fabs(det1)>=fabs(det2)) && (fabs(det1)>=fabs(det3)) ) 
    {
    n1=0;
    n2=1;
    det=det1;
    }
else
    {
    if ( (fabs(det2)>=fabs(det1)) && (fabs(det2)>=fabs(det3)) ) 
        {
        n1=0;
        n2=2;
       det=det2;
        }
    else
        {
        if ( (fabs(det3)>=fabs(det2)) && (fabs(det3)>=fabs(det1)) ) 
          {
          n1=1;
          n2=2;
          det=det3;
          }
        else return 0;  
        }
    }
double tt=ab[n1]*ac[n2]-ab[n2]*ac[n1];
tt=tt/det;
if (tt<0.0000000001) return 0;
if (tt>1.-0.0000000001) return 0;
xyz[0]=nos1->get_x()+tt*(nos2->get_x()-nos1->get_x());
xyz[1]=nos1->get_y()+tt*(nos2->get_y()-nos1->get_y());
xyz[2]=nos1->get_z()+tt*(nos2->get_z()-nos1->get_z());
return 1;
}


void TOIBREP::testtopo(void)
{
LISTE_XFEM_ELEMENT0::iterator it0;
LISTE_XFEM_ELEMENT1::iterator it1;
LISTE_XFEM_ELEMENT2::iterator it2;
LISTE_XFEM_ELEMENT3::iterator it3;
for (XFEM_ELEMENT0 *ele=mai->get_premier_xelement0(it0);ele!=NULL;ele=mai->get_suivant_xelement0(it0))
  {
  if (ele->get_fem_noeud(0)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 0 " << ele->get_fem_noeud(0)->get_id() << " non attache à la topologie" << std::endl; 
  }
for (XFEM_ELEMENT1 *ele=mai->get_premier_xelement1(it1);ele!=NULL;ele=mai->get_suivant_xelement1(it1))
  {
  if (ele->get_lien_topologie()==NULL) std::cout << "Xele1 " << ele->get_id() << " non attache à la topologie " << std::endl;
  if (ele->get_fem_noeud(0)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 1 " << ele->get_fem_noeud(0)->get_id() << " non attache à la topologie. " << ele->get_fem_noeud(0)->get_lien_element3()->get_nb() << " tetras lies"   << std::endl; 
  if (ele->get_fem_noeud(1)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 1 " << ele->get_fem_noeud(1)->get_id() << " non attache à la topologie. " << ele->get_fem_noeud(0)->get_lien_element3()->get_nb() << " tetras lies"  << std::endl; 
  }
for (XFEM_ELEMENT2 *ele=mai->get_premier_xelement2(it2);ele!=NULL;ele=mai->get_suivant_xelement2(it2))
  {
  if (ele->get_lien_topologie()==NULL) std::cout << "Xele2 " << ele->get_id() << " non attache à la topologie " << std::endl;
  if (ele->get_fem_noeud(0)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 2 " << ele->get_fem_noeud(0)->get_id() << " non attache à la topologie" << std::endl; 
  if (ele->get_fem_noeud(1)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 2 " << ele->get_fem_noeud(1)->get_id() << " non attache à la topologie" << std::endl; 
  }
/*for (XFEM_ELEMENT3 *ele=mai->get_premier_xelement3(it3);ele!=NULL;ele=mai->get_suivant_xelement3(it3))
  {
  if (ele->get_fem_noeud(0)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 3" << ele->get_fem_noeud(0)->get_id() << " non attache à la topologie" << std::endl; 
  if (ele->get_fem_noeud(1)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 3" << ele->get_fem_noeud(1)->get_id() << " non attache à la topologie" << std::endl; 
  if (ele->get_fem_noeud(2)->get_lien_topologie()==NULL) std::cout << "Noeud de xele 3" << ele->get_fem_noeud(2)->get_id() << " non attache à la topologie" << std::endl; 
  }*/


}

Revision:	488
Committed:	Mon Feb 10 18:09:02 2014 UTC (11 years, 3 months ago) by francois
File size:	122668 byte(s)
Log Message:	resolution de mise en module independant du code d'eric