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"


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

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

TOIBREP::~TOIBREP()
{
}


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

int TOIBREP::compare_etat(int etat,int valeur)
{
int ope=etat&valeur;
if (ope==valeur) 
  return 1;
return 0;
}


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


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

//void TOIBREP::importer(std::string nomfichier,class MagXchange* data,std::string nomfichierout)
IBrep TOIBREP::importer(std::string nomfichier,std::string nomfichieribrep,MG_GROUPE_TOPOLOGIQUE* mggt)
{
affiche((char*)"Traitement de base");
TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lst;
int nbface=geo->get_nb_mg_face();
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);
                }
        }
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_tetra=mai->get_nb_fem_tetra();
std::string nomfichier2=nomfichier+"_dis.sol";
FEM_SOLUTION* solution=new FEM_SOLUTION(mai,nbvraiface,(char*)nomfichier2.c_str(),nb_noeud,"LS",ENTITE_NOEUD);
int i=0;
gest->ajouter_fem_solution(solution);
nomfichier2=nomfichier+"_ele.sol";
FEM_SOLUTION* solution_ele=new FEM_SOLUTION(mai,nbvraiface,(char*)nomfichier2.c_str(),nb_tetra,"A",ENTITE_TETRA);
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;
        }

//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;j++)
                solution->ecrire(i,j,1e300);
        lst_noeud.ajouter(noeud);
        i++;
        }
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base recherche boite");
octree.initialiser(&lst_noeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base ini octree");
LISTE_FEM_TETRA::iterator it3;
i=0;
for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3))
  {
  octree.inserer(tet);
  for (int j=0;j<nbvraiface;j++)
      solution_ele->ecrire(i,j,0);
  i++;
  }
if (compteur!=NULL) compteur->ajouter_etape("Traitement de base");
for (int i=0;i<nbface;i++)
        {
        if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"      Face %d",i);
          affiche(mess);
          }
        MG_FACE* face=geo->get_mg_face(i);
        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);
        }
affiche((char*)"Exportation IBREP");
IBrep brep=exporter_IBrep(nomfichieribrep,solution,solution_ele,mggt);
if (compteur!=NULL) compteur->ajouter_etape("Exportation IBREP");
return brep;
}

void TOIBREP::levelset0(FEM_SOLUTION *solution,int numsol)
{
/*solution->active_solution(numsol);
int nbface=geo->get_nb_mg_face();
vector<TOIBREP_POINT*> lst;
for (int iface=0;iface<nbface;iface++)
        {
        MG_FACE* face=geo->get_mg_face(iface);
        int nbpoint=lst.size();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
        for (MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get_premier(it);ele!=NULL;ele=face->get_lien_maillage()->get_suivant(it))
                {
                MG_TRIANGLE *tri=(Mvoid TOIBREP_POINT::get_coord2(double *uv)
{
uv[0]=u;
uv[1]=v;
}G_TRIANGLE*)ele;
                TOIBREP_POINT *pt1=new TOIBREP_POINT(tri->get_noeud1()->get_x(),tri->get_noeud1()->get_y(),tri->get_noeud1()->get_z(),tri->get_noeud1()->get_lien_topologie());
                TOIBREP_POINT *pt2=new TOIBREP_POINT(tri->get_noeud2()->get_x(),tri->get_noeud2()->get_y(),tri->get_noeud2()->get_z(),tri->get_noeud2()->get_lien_topologie());
                TOIBREP_POINT *pt3=new TOIBREP_POINT(tri->get_noeud3()->get_x(),tri->get_noeud3()->get_y(),tri->get_noeud3()->get_z(),tri->get_noeud3()->get_lien_topologie());
                lst.push_back(pt1);
                lst.push_back(pt2);
                lst.push_back(pt3);
                }
        calcullevelsetpremierepasse(solution,numsol,face,8,&lst,nbpoint,lst.size());
        }
calcullevelsetdeuxiemepasse(solution,numsol,&lst);*/
}
void TOIBREP::echantillonne_sommet(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_arete(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face)
{
echantillonne_sommet(lst,face);
int nbbou=face->get_nb_mg_boucle();
for (int i=0;i<nbbou;i++)
    {
    MG_BOUCLE* boucle=face->get_mg_boucle(i);
    int nb_arete=boucle->get_nb_mg_coarete();
    for (int j=0;j<nb_arete;j++)
      {
      MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete();
      double tmin=are->get_tmin();
      double tmax=are->get_tmax();
      double N=NPAS;
      for (int k=0;k<N+1;k++)
          {
          double t=tmin+k*1.0/N*(tmax-tmin);
          double xyz[3];
          are->evaluer(t,xyz);
          TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
          lst.push_back(pt);
//        octree.inserer(pt);
          }
     /* double t=tmin+0.2/NPAS*(tmax-tmin);
      double xyz[3];
      are->evaluer(t,xyz);
      TOIBREP_POINT *pt1=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
      lst.push_back(pt1);
      t=tmax-0.2/NPAS*(tmax-tmin);
      are->evaluer(t,xyz);
      TOIBREP_POINT *pt2=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
      lst.push_back(pt2);*/
      
      
      }
    }
}

void TOIBREP::levelsetn(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lsttopo,vector<MG_FACE*> *lstface,class FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele,int numsol)
{
if (affichageactif)
          { 
          char mess[300];
          sprintf(mess,"          Echantillonnage");
          affiche(mess);
          }
        
solution->active_solution(numsol);
solution_ele->active_solution(numsol);
int nbface=lstface->size();
if (nbface==0) return;
vector<TOIBREP_POINT*> lst;
BOITE_2D boite=ot.get_boite_2D((*lstface)[0]);
echantillonne_arete(lst,(*lstface)[0]);
for (int iface=1;iface<nbface;iface++)
        {
        MG_FACE* face=(*lstface)[iface];
        BOITE_2D boiteface=ot.get_boite_2D(face);
        boite=boite+boiteface;
        echantillonne_arete(lst,(*lstface)[iface]);
        }
MG_FACE* face=(*lstface)[0];
int trouve=0;
int orientation;
int j=0;
do 
        {
        MG_COFACE* coface=face->get_mg_coface(j);
        MG_VOLUME* vol=coface->get_coquille()->get_mg_volume();
        if (lsttopo->existe(vol)) 
                {
                orientation=coface->get_orientation();
                trouve=1;
                }
        j++;
        }
while (trouve==0);
LISTE_FEM_TETRA::iterator ittetra;
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra))
        tet->change_etat(0);
LISTE_FEM_NOEUD::iterator itnoeud;
for (FEM_NOEUD* no=mai->get_premier_noeud(itnoeud);no!=NULL;no=mai->get_suivant_noeud(itnoeud))
        no->change_etat(0);
double umin=boite.get_xmin();
double umax=boite.get_xmax();
double vmin=boite.get_ymin();
double vmax=boite.get_ymax();
if (face->get_surface()->est_periodique_u())
        {
        umin=0;
        umax=face->get_surface()->get_periode_u();
        }
if (face->get_surface()->est_periodique_v())
        {
        vmin=0;
        vmax=face->get_surface()->get_periode_v();
        }
for (int i=-5;i<NPAS+5;i++)
        for (int j=-5;j<NPAS+5;j++)
                {
                double uv[2];
                uv[0]=umin+i*1.0/NPAS*(umax-umin);
                uv[1]=vmin+j*1.0/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);
//                      octree.inserer(pt);
                        }
                }
// debut temp
/*
MG_GESTIONNAIRE gestemp;
MG_MAILLAGE *mai=new MG_MAILLAGE(NULL);
gestemp.ajouter_mg_maillage(mai);
for (int i=0;i<lst.size();i++)
{
double xyz[3];
lst[i]->get_coord3(xyz);
MG_NOEUD* no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],TRIANGULATION);
mai->ajouter_mg_noeud(no);
if (i!=0) {MG_SEGMENT* seg=new MG_SEGMENT(NULL,no,mai->get_mg_noeud(mai->get_nb_mg_noeud()-2),TRIANGULATION);mai->ajouter_mg_segment(seg);}
}

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

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


}

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

void TOIBREP::decoupe_tetra(FEM_TETRA* 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_TETRA::iterator it3;
for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(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();
}*/
/*LISTE_FEM_TETRA::iterator it3;
for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(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_TETRA*> liste;
                octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste);
                TPL_MAP_ENTITE<FEM_TETRA*>::ITERATEUR it;
                for (FEM_TETRA* tet=liste.get_premier(it);tet!=NULL;tet=liste.get_suivant(it))
                        {
                        if (tet->get_etat()==1) continue;
                        int etat=estdansletetra(tet,xyz[0],xyz[1],xyz[2]);
                        if (compare_etat(etat,INTERIEUR))
                                {
                                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                                        {
                                        double dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_fem_noeud(k));
                                        if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution()))
                                                {
                                                tet->get_fem_noeud(k)->change_solution(dist);
                                                tet->get_fem_noeud(k)->change_numero(pt->get_id());
                                                pt->change_coord2(uv);
                                                }
                                        }
                                }
                        if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur()))
                                {
                                tet->change_etat(1);
                                for (int k=0;k<tet->get_nb_fem_noeud();k++)
                                  {
                                  tet->get_fem_noeud(k)->change_etat(1);
                                  //octree.supprimer(tet);
                                  }
                                }
                        }
                }
//octree.vide();


}       

void TOIBREP::calcullevelsetdeuxiemepasse(vector<TOIBREP_POINT*> *lst)
{
LISTE_FEM_NOEUD::iterator itno;
for (FEM_NOEUD* no=mai->get_premier_noeud(itno);no!=NULL;no=mai->get_suivant_noeud(itno))
        if (no->get_solution()<1e250)
                {
                TOIBREP_POINT* pt=(*lst)[no->get_numero()];
                int dim=pt->get_mg_element_topologique()->get_dimension();
                if (dim==2)
                        {
                        MG_FACE* face=(MG_FACE*)pt->get_mg_element_topologique();
                        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);
                        }
                if (dim==1)
                        {
                        MG_ARETE* arete=(MG_ARETE*)pt->get_mg_element_topologique();
                        double dist=calcul_distance(no,arete,pt);
                        no->change_solution(dist);
                        }
                if (dim==0)
                        {
                        MG_SOMMET* sommet=(MG_SOMMET*)pt->get_mg_element_topologique();
                        MG_ARETE* arete=sommet->get_mg_cosommet(0)->get_arete();
                        double dist=calcul_distance(no,arete,pt);
                        no->change_solution(dist);
                        }
                 
                }


}

double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_ARETE* are,TOIBREP_POINT* pt,double precision)
{
double tii,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);
double xyz[3],dxyz[3];
are->evaluer(tii,xyz);
OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]);
double distance=(Ct-Pt).get_longueur();
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();
are->deriver(tii,dxyz);
OT_VECTEUR_3D x1(dxyz[0],dxyz[1],dxyz[2]);
OT_VECTEUR_3D x2(dxyz[0],dxyz[1],dxyz[2]);
x1=are->get_mg_coarete(0)->get_orientation()*x1;
x2=are->get_mg_coarete(1)->get_orientation()*x2;
x1.norme();
x2.norme();
double uv1[2],uv2[2];
double normal1[3],normal2[3];
face1->inverser(uv1,xyz);
face2->inverser(uv2,xyz);
face1->calcul_normale_unitaire(uv1,normal1);
face2->calcul_normale_unitaire(uv2,normal2);
OT_VECTEUR_3D z1(normal1[0],normal1[1],normal1[2]);
OT_VECTEUR_3D z2(normal2[0],normal2[1],normal2[2]);
z1=face1->get_mg_coface(0)->get_orientation()*z1;
z2=face2->get_mg_coface(0)->get_orientation()*z2;
OT_VECTEUR_3D y1=z1&x1;
OT_VECTEUR_3D y2=z2&x2;
double test=(z1&z2)*x1;
int signe=-1;
OT_VECTEUR_3D dirpt=Pt-Ct;
if (test>0)
        {if ((z1*dirpt>0.) || (z2*dirpt>0.)) signe=1;}
else
        {if ((z1*dirpt>0.) && (z2*dirpt>0.)) signe=1;}
return signe*distance;
}

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_TETRA::iterator ittet;
for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(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_TETRA* 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_TETRA* 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_tetra()->get_nb();
        for (int i=0;i<nbtetra;i++)
                {
                FEM_TETRA* tet2=noeud->get_lien_tetra()->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_TETRA* 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_tetra()->get_nb();
                                for (int i=0;i<nbtetra;i++)
                                        {
                                        FEM_TETRA* tet3=noeud->get_lien_tetra()->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_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
        {
        if (tet->get_solution()>1e200) continue;
        tet->get_fem_noeud(0)->change_numero(1);
        tet->get_fem_noeud(1)->change_numero(1);
        tet->get_fem_noeud(2)->change_numero(1);
        tet->get_fem_noeud(3)->change_numero(1);
        }
int i=0;        
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
        {
        if (noeud->get_numero()==1) sol->ecrire(i,numsol,noeud->get_solution()); else sol->ecrire(i,numsol,1e300);
        ++i;
        }
}
IBrep TOIBREP::exporter_IBrep(string chemin,FEM_SOLUTION* solution,FEM_SOLUTION *solution_ele,MG_GROUPE_TOPOLOGIQUE* mggt)
{
TPL_MAP_ENTITE<MG_VOLUME*> lst;
if (mggt==NULL)
        {
        int nbvol=geo->get_nb_mg_volume();
        for (int i=0;i<nbvol;i++)
                lst.ajouter(geo->get_mg_volume(i));
        }
else 
        {
        int nbmggt=mggt->get_nb();
        for (int i=0;i<nbmggt;i++)
                if (mggt->get(i)->get_dimension()==3) 
                        lst.ajouter((MG_VOLUME*)mggt->get(i));
        }
IBrep brep(100);
int nbvolexp=lst.get_nb();
for (int i=0;i<nbvolexp;i++)
        {
        MG_VOLUME* vol=lst.get(i);
        int nbcoq=vol->get_nb_mg_coquille();
        IVolumeN ivoltmp(vol->get_id());
        pair<IVolume*,bool> pair_ivol=brep.AddVolume(ivoltmp);
        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);
                                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);
                                                                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);
                                                                ivertex2=pair_ivertex2.first;
                                                                }
                                                        ivertex1->AddCoVertex(cover1->get_id());
                                                        ivertex2->AddCoVertex(cover2->get_id());
                                                        pair<IEdge*,bool> pair_iedge=brep.AddEdge(iedgetmp);
                                                        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);
        correspondface[i]=id;
        }
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;
        for (int j=0;j<nbsol;j++) 
                {
                if (solution->lire(i,j)<1e200)
                        {
                        inoeud->Fields()[num]=correspondface[j];
                        inoeud->Values()[num]=solution->lire(i,j);
                        num++;
                        } 
                }
        i++;
        }
LISTE_FEM_TETRA::iterator ittet;
int num=0;
i=0;
for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(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;
        for (int j=0;j<nbsol;j++) 
                {
                if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001)
                  {
                  xfemele.Fields()[num]=correspondface[j];
                  num++;
                  }
                }
        brep.AddElement(xfemele);
        i++;
        }

delete [] correspondface;


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


return brep;
}


void TOIBREP::ajouter_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_TETRA* tet,double val)
{
pair<double,FEM_TETRA*> 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_TETRA* 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_TETRA* tet)
{
lst.push_back(tet);
}

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


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


tet->get_inverse_jacob(j,uv);
for (int i=0;i<3;i++)
        for (int k=0;k<4;k++)
                N[i*4+k]=tet->get_fonction_derive_interpolation(k+1,i+1,uv);
int premier=0;
double tmin=1e300;
double tmax=-1e300;
for (int i=0;i<4;i++)
        {
        if (i==tet->get_numero()) continue;
        T[i]=tet->get_fem_noeud(i)->get_solution();
        if (fabs(T[i])>0.000001)
                {
                if (premier==0)
                        if (T[i]>0) (*signe)=1; else (*signe)=-1;
                else if (T[i]*(*signe)<0) (*signe)=0;
                }
        T[i]=fabs(T[i]);
        if (tet->get_numero()!=i)
                {
                if (T[i]<tmin) tmin=T[i];
                if (T[i]>tmax) tmax=T[i];
                }
        premier=1;
        }
for (int i=0;i<3;i++)
        for (int k=0;k<4;k++)
                for (int l=0;l<3;l++)
                        jN[i*4+k]=jN[i*4+k]+j[i*3+l]*N[l*4+k];
double a=0.,b=0.,c=-1.;
for (int i=0;i<3;i++)
        {
        double coef=0.;
        double coefinc=0.;
        for (int l=0;l<4;l++)
                {
                if (tet->get_numero()!=l) coef=coef+jN[i*4+l]*T[l];
                else coefinc=coefinc+jN[i*4+l];
                }
        c=c+coef*coef;
        a=a+coefinc*coefinc;
        b=b+2*coef*coefinc;
        }
/*if (*signe==0)
        cout << "attention " <<endl;*/
double det=b*b-4.*a*c;
if (det<0.) det=0.; else det=sqrt(det);
double sol1=(-b-det)/2./a;
double sol2=(-b+det)/2./a;
double sol=sol1;
if (sol2>sol1) sol=sol2;
if (sol<tmin*0.99)
        sol=0.;
sol=sol*(*signe);
return sol;
}
Revision:	280
Committed:	Thu Jun 30 19:51:22 2011 UTC (14 years, 2 months ago) by francois
File size:	44691 byte(s)
Log Message:	progression de couche dans le bouchage de trou lors de la construction des level-sets dans toIbrep