ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/toIbrep/src/toibrep.cpp
Revision: 280
Committed: Thu Jun 30 19:51:22 2011 UTC (13 years, 10 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

File Contents

# User Rev Content
1 bechet 105 #include "gestionversion.h"
2 francois 259 #include "toibrep.h"
3 francois 104 #include "mg_file.h"
4     #include "ot_mathematique.h"
5     #include "tpl_map_entite.h"
6 francois 222 #include "tpl_relation_entite.h"
7 francois 104 #include "tpl_octree.h"
8 francois 259 #include "toibrep_point.h"
9 francois 106 #include "vct.h"
10 francois 104 #include <math.h>
11 francois 259 #include <fstream>
12 francois 262 #include <string>
13 francois 259 #include "IBrep.h"
14 francois 276 #include "ot_cpu.h"
15 francois 104
16    
17 francois 222
18 francois 276 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)
19 francois 104 {
20     }
21    
22 francois 276 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)
23     {
24     }
25     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)
26     {
27     }
28     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)
29     {
30     }
31 francois 104
32 francois 259 TOIBREP::~TOIBREP()
33 francois 104 {
34     }
35    
36    
37    
38 francois 276 void TOIBREP::active_affichage(void (*fonc)(char*))
39     {
40     affiche=fonc;
41     affichageactif=1;
42     }
43    
44     int TOIBREP::compare_etat(int etat,int valeur)
45     {
46     int ope=etat&valeur;
47     if (ope==valeur)
48     return 1;
49     return 0;
50     }
51    
52    
53    
54 francois 259 int TOIBREP::estdansletetra(FEM_TETRA *tet,double x,double y, double z)
55 francois 104 {
56 francois 222 double *xyz1,*xyz2,*xyz3,*xyz4;
57     if (tet->get_nb_fem_noeud()==4)
58     {
59     xyz1=tet->get_fem_noeud(0)->get_coord();
60     xyz2=tet->get_fem_noeud(1)->get_coord();
61     xyz3=tet->get_fem_noeud(2)->get_coord();
62     xyz4=tet->get_fem_noeud(3)->get_coord();
63     }
64     if (tet->get_nb_fem_noeud()==10)
65     {
66     xyz1=tet->get_fem_noeud(0)->get_coord();
67     xyz2=tet->get_fem_noeud(2)->get_coord();
68     xyz3=tet->get_fem_noeud(4)->get_coord();
69     xyz4=tet->get_fem_noeud(9)->get_coord();
70     }
71 francois 104 OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]);
72     OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]);
73     OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]);
74     OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]);
75     OT_MATRICE_3D mat(v1,v2,v3);
76     OT_MATRICE_3D mat1(v4,v2,v3);
77     OT_MATRICE_3D mat2(v1,v4,v3);
78     OT_MATRICE_3D mat3(v1,v2,v4);
79     double det=mat.get_determinant();
80     double xsi=mat1.get_determinant()/det;
81     double eta=mat2.get_determinant()/det;
82     double dseta=mat3.get_determinant()/det;
83 francois 276 int reponse1=1;
84     int reponse2=1;
85     double eps=0.0001;
86     if (xsi<-eps) reponse1=0;
87     if (eta<-eps) reponse1=0;
88     if (dseta<-eps) reponse1=0;
89     if (xsi+eta+dseta>1.+eps) reponse1=0;
90     if (xsi<eps) reponse2=0;
91     if (eta<eps) reponse2=0;
92     if (dseta<eps) reponse2=0;
93     if (xsi+eta+dseta>1.-eps) reponse2=0;
94     return reponse1+2*reponse2;
95 francois 104 }
96    
97 francois 276
98    
99 francois 259 double TOIBREP::calculdist(double *n,double x,double y,double z,FEM_NOEUD* noeud)
100 francois 104 {
101     double* xyz=noeud->get_coord();
102     double dist=sqrt((xyz[0]-x)*(xyz[0]-x)+(xyz[1]-y)*(xyz[1]-y)+(xyz[2]-z)*(xyz[2]-z));
103     OT_VECTEUR_3D vec1(n[0],n[1],n[2]);
104     OT_VECTEUR_3D vec2(xyz[0]-x,xyz[1]-y,xyz[2]-z);
105     double ps=vec1*vec2;
106     if (ps<0.) dist=-dist;
107     return dist;
108     }
109    
110 francois 259 //void TOIBREP::importer(std::string nomfichier,class MagXchange* data,std::string nomfichierout)
111 francois 276 IBrep TOIBREP::importer(std::string nomfichier,std::string nomfichieribrep,MG_GROUPE_TOPOLOGIQUE* mggt)
112 francois 104 {
113 francois 276 affiche((char*)"Traitement de base");
114 francois 222 TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lst;
115 francois 106 int nbface=geo->get_nb_mg_face();
116 francois 222 if (mggt!=NULL)
117     {
118     int nb=mggt->get_nb();
119     for (int i=0;i<nb;i++)
120     {
121     lst.ajouter(mggt->get(i));
122     mggt->get(i)->get_topologie_sousjacente(&lst);
123     }
124     }
125     else
126     {
127     int nbvol=geo->get_nb_mg_volume();
128     for (int i=0;i<nbvol;i++)
129     {
130     lst.ajouter(geo->get_mg_volume(i));
131     geo->get_mg_volume(i)->get_topologie_sousjacente(&lst);
132     }
133     }
134     int nbvraiface=0;
135     TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*>::ITERATEUR it;
136     for (MG_ELEMENT_TOPOLOGIQUE* ele=lst.get_premier(it);ele!=NULL;ele=lst.get_suivant(it))
137     if (ele->get_dimension()==2) nbvraiface++;
138 francois 276 //if (compteur!=NULL) compteur->ajouter_etape("Traitement de base selection");
139 francois 222 int nb_noeud=mai->get_nb_fem_noeud();
140 francois 276 int nb_tetra=mai->get_nb_fem_tetra();
141     std::string nomfichier2=nomfichier+"_dis.sol";
142     FEM_SOLUTION* solution=new FEM_SOLUTION(mai,nbvraiface,(char*)nomfichier2.c_str(),nb_noeud,"LS",ENTITE_NOEUD);
143 francois 222 int i=0;
144 francois 276 gest->ajouter_fem_solution(solution);
145     nomfichier2=nomfichier+"_ele.sol";
146     FEM_SOLUTION* solution_ele=new FEM_SOLUTION(mai,nbvraiface,(char*)nomfichier2.c_str(),nb_tetra,"A",ENTITE_TETRA);
147     gest->ajouter_fem_solution(solution_ele);
148 francois 222 for (MG_ELEMENT_TOPOLOGIQUE* ele=lst.get_premier(it);ele!=NULL;ele=lst.get_suivant(it))
149 francois 106 {
150 francois 222 if (ele->get_dimension()!=2) continue;
151 francois 106 char mess[255];
152 francois 262 sprintf(mess,"F %lu",ele->get_id());
153 francois 276 solution->change_legende(i,mess);
154     solution_ele->change_legende(i,mess);
155 francois 222 ++i;
156 francois 106 }
157 francois 276
158     //if (compteur!=NULL) compteur->ajouter_etape("Traitement de base creation solution");
159 francois 104 double xmin=1e300,ymin=1e300,zmin=1e308;
160     double xmax=-1e300,ymax=-1e300,zmax=-1e308;
161 francois 222 TPL_MAP_ENTITE<FEM_NOEUD*> lst_noeud;
162     LISTE_FEM_NOEUD::iterator it2;
163     i=0;
164     for (FEM_NOEUD* noeud=mai->get_premier_noeud(it2);noeud!=NULL;noeud=mai->get_suivant_noeud(it2))
165 francois 104 {
166     double* xyz=noeud->get_coord();
167     if (xyz[0]<xmin) xmin=xyz[0];
168     if (xyz[1]<ymin) ymin=xyz[1];
169     if (xyz[2]<zmin) zmin=xyz[2];
170     if (xyz[0]>xmax) xmax=xyz[0];
171     if (xyz[1]>ymax) ymax=xyz[1];
172     if (xyz[2]>zmax) zmax=xyz[2];
173 francois 276 for (int j=0;j<nbvraiface;j++)
174 francois 106 solution->ecrire(i,j,1e300);
175 francois 104 lst_noeud.ajouter(noeud);
176 francois 222 i++;
177 francois 104 }
178 francois 276 //if (compteur!=NULL) compteur->ajouter_etape("Traitement de base recherche boite");
179 francois 104 octree.initialiser(&lst_noeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
180 francois 276 //if (compteur!=NULL) compteur->ajouter_etape("Traitement de base ini octree");
181 francois 222 LISTE_FEM_TETRA::iterator it3;
182 francois 276 i=0;
183 francois 222 for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3))
184 francois 276 {
185     octree.inserer(tet);
186     for (int j=0;j<nbvraiface;j++)
187     solution_ele->ecrire(i,j,0);
188     i++;
189     }
190     if (compteur!=NULL) compteur->ajouter_etape("Traitement de base");
191 francois 106 for (int i=0;i<nbface;i++)
192     {
193 francois 276 if (affichageactif)
194     {
195     char mess[300];
196     sprintf(mess," Face %d",i);
197     affiche(mess);
198     }
199 francois 222 MG_FACE* face=geo->get_mg_face(i);
200     if (!lst.existe(face)) continue;
201 francois 106 vector<MG_FACE*> lstface;
202 francois 222 lstface.push_back(face);
203 francois 276 /// ici topologie virtuelle
204    
205    
206 francois 222 /*TPL_LISTE_ENTITE<double> lst;
207 francois 106 int type=face->get_surface()->get_type_geometrique(lst);
208     int idem=0;
209     for (int j=0;j<i;j++)
210     {
211     TPL_LISTE_ENTITE<double> lst2;
212     MG_FACE* face2=geo->get_mg_face(j);
213     int type2=face2->get_surface()->get_type_geometrique(lst2);
214     if (type==type2)
215     if (lst.get_nb()==lst2.get_nb())
216     {
217     int diff=0;
218     for (int i=0;i<lst.get_nb();i++)
219     if (fabs(lst.get(i)-lst2.get(i))>0.000001) diff=1;
220     if (diff==0) idem=1;
221     }
222     }
223     if (!idem) lstface.push_back(geo->get_mg_face(i));
224     for (int j=i+1;j<nbface;j++)
225     {
226     TPL_LISTE_ENTITE<double> lst2;
227     MG_FACE* face2=geo->get_mg_face(j);
228     int type2=face2->get_surface()->get_type_geometrique(lst2);
229     if (type==type2)
230     if (lst.get_nb()==lst2.get_nb())
231     {
232     int diff=0;
233     for (int i=0;i<lst.get_nb();i++)
234     if (fabs(lst.get(i)-lst2.get(i))>0.000001) diff=1;
235     if (diff==0) lstface.push_back(face2);
236     }
237 francois 222 }*/
238 francois 276 levelsetn(&lst,&lstface,solution,solution_ele,i);
239 francois 106 }
240 francois 276 affiche((char*)"Exportation IBREP");
241     IBrep brep=exporter_IBrep(nomfichieribrep,solution,solution_ele,mggt);
242     if (compteur!=NULL) compteur->ajouter_etape("Exportation IBREP");
243     return brep;
244 francois 106 }
245    
246 francois 259 void TOIBREP::levelset0(FEM_SOLUTION *solution,int numsol)
247 francois 106 {
248 francois 276 /*solution->active_solution(numsol);
249 francois 104 int nbface=geo->get_nb_mg_face();
250 francois 259 vector<TOIBREP_POINT*> lst;
251 francois 104 for (int iface=0;iface<nbface;iface++)
252     {
253 francois 106 MG_FACE* face=geo->get_mg_face(iface);
254     int nbpoint=lst.size();
255 francois 104 TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
256     for (MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get_premier(it);ele!=NULL;ele=face->get_lien_maillage()->get_suivant(it))
257     {
258 francois 276 MG_TRIANGLE *tri=(Mvoid TOIBREP_POINT::get_coord2(double *uv)
259     {
260     uv[0]=u;
261     uv[1]=v;
262     }G_TRIANGLE*)ele;
263 francois 259 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());
264     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());
265     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());
266 francois 106 lst.push_back(pt1);
267     lst.push_back(pt2);
268     lst.push_back(pt3);
269 francois 104 }
270 francois 222 calcullevelsetpremierepasse(solution,numsol,face,8,&lst,nbpoint,lst.size());
271 francois 106 }
272 francois 276 calcullevelsetdeuxiemepasse(solution,numsol,&lst);*/
273 francois 106 }
274 francois 276 void TOIBREP::echantillonne_sommet(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face)
275 francois 106 {
276 francois 276 TPL_MAP_ENTITE<MG_SOMMET*> lstsom;
277     int nbbou=face->get_nb_mg_boucle();
278     for (int i=0;i<nbbou;i++)
279     {
280     MG_BOUCLE* boucle=face->get_mg_boucle(i);
281     int nb_arete=boucle->get_nb_mg_coarete();
282     for (int j=0;j<nb_arete;j++)
283     {
284     MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete();
285     lstsom.ajouter(are->get_cosommet1()->get_sommet());
286     lstsom.ajouter(are->get_cosommet2()->get_sommet());
287     }
288     }
289     TPL_MAP_ENTITE<MG_SOMMET*>::ITERATEUR it;
290     for (MG_SOMMET* som=lstsom.get_premier(it);som!=NULL;som=lstsom.get_suivant(it))
291     {
292     double xyz[3];
293     som->get_point()->evaluer(xyz);
294     double uv[2];
295     face->inverser(uv,xyz);
296     int N=NPAS;
297     double rayon=1e-2;
298     for (int k=0;k<N+1;k++)
299     {
300     double uv2[2];
301     uv2[0]=uv[0]+rayon*cos(k*1.0/N*2.*M_PI);
302     uv2[1]=uv[1]+rayon*sin(k*1.0/N*2.*M_PI);
303     if ((face->valide_parametre_u(uv2[0])) && (face->valide_parametre_v(uv2[1])))
304     {
305     double distance=ot.calcule_distance_contour_face_uv(uv2,face);
306     if (distance>0.)
307     {
308     double xyz[3];
309     face->evaluer(uv2,xyz);
310     TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],1,face);
311     lst.push_back(pt);
312     }
313     }
314     }
315    
316     }
317     }
318     void TOIBREP::echantillonne_arete(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face)
319     {
320     echantillonne_sommet(lst,face);
321     int nbbou=face->get_nb_mg_boucle();
322     for (int i=0;i<nbbou;i++)
323     {
324     MG_BOUCLE* boucle=face->get_mg_boucle(i);
325     int nb_arete=boucle->get_nb_mg_coarete();
326     for (int j=0;j<nb_arete;j++)
327     {
328     MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete();
329     double tmin=are->get_tmin();
330     double tmax=are->get_tmax();
331     double N=NPAS;
332     for (int k=0;k<N+1;k++)
333     {
334     double t=tmin+k*1.0/N*(tmax-tmin);
335     double xyz[3];
336     are->evaluer(t,xyz);
337     TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
338     lst.push_back(pt);
339     // octree.inserer(pt);
340     }
341     /* double t=tmin+0.2/NPAS*(tmax-tmin);
342     double xyz[3];
343     are->evaluer(t,xyz);
344     TOIBREP_POINT *pt1=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
345     lst.push_back(pt1);
346     t=tmax-0.2/NPAS*(tmax-tmin);
347     are->evaluer(t,xyz);
348     TOIBREP_POINT *pt2=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face);
349     lst.push_back(pt2);*/
350    
351    
352    
353     }
354     }
355     }
356    
357     void TOIBREP::levelsetn(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lsttopo,vector<MG_FACE*> *lstface,class FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele,int numsol)
358     {
359     if (affichageactif)
360     {
361     char mess[300];
362     sprintf(mess," Echantillonnage");
363     affiche(mess);
364     }
365    
366     solution->active_solution(numsol);
367     solution_ele->active_solution(numsol);
368     int nbface=lstface->size();
369     if (nbface==0) return;
370     vector<TOIBREP_POINT*> lst;
371     BOITE_2D boite=ot.get_boite_2D((*lstface)[0]);
372     echantillonne_arete(lst,(*lstface)[0]);
373     for (int iface=1;iface<nbface;iface++)
374     {
375     MG_FACE* face=(*lstface)[iface];
376     BOITE_2D boiteface=ot.get_boite_2D(face);
377     boite=boite+boiteface;
378     echantillonne_arete(lst,(*lstface)[iface]);
379     }
380     MG_FACE* face=(*lstface)[0];
381     int trouve=0;
382     int orientation;
383     int j=0;
384     do
385     {
386     MG_COFACE* coface=face->get_mg_coface(j);
387     MG_VOLUME* vol=coface->get_coquille()->get_mg_volume();
388     if (lsttopo->existe(vol))
389 francois 104 {
390 francois 276 orientation=coface->get_orientation();
391     trouve=1;
392     }
393     j++;
394     }
395     while (trouve==0);
396     LISTE_FEM_TETRA::iterator ittetra;
397     for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra))
398     tet->change_etat(0);
399     LISTE_FEM_NOEUD::iterator itnoeud;
400     for (FEM_NOEUD* no=mai->get_premier_noeud(itnoeud);no!=NULL;no=mai->get_suivant_noeud(itnoeud))
401     no->change_etat(0);
402     double umin=boite.get_xmin();
403     double umax=boite.get_xmax();
404     double vmin=boite.get_ymin();
405     double vmax=boite.get_ymax();
406     if (face->get_surface()->est_periodique_u())
407     {
408     umin=0;
409     umax=face->get_surface()->get_periode_u();
410     }
411     if (face->get_surface()->est_periodique_v())
412     {
413     vmin=0;
414     vmax=face->get_surface()->get_periode_v();
415     }
416     for (int i=-5;i<NPAS+5;i++)
417     for (int j=-5;j<NPAS+5;j++)
418     {
419     double uv[2];
420     uv[0]=umin+i*1.0/NPAS*(umax-umin);
421     uv[1]=vmin+j*1.0/NPAS*(vmax-vmin);
422     double xyz[3];
423     if ((face->valide_parametre_u(uv[0])) && (face->valide_parametre_v(uv[1])))
424     {
425     double distance=ot.calcule_distance_contour_face_uv(uv,face);
426     int inte=0;
427     if (distance>0.) inte=1.; else continue;
428     face->evaluer(uv,xyz);
429     TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],inte,face);
430     lst.push_back(pt);
431     // octree.inserer(pt);
432     }
433     }
434     // debut temp
435     /*
436     MG_GESTIONNAIRE gestemp;
437     MG_MAILLAGE *mai=new MG_MAILLAGE(NULL);
438     gestemp.ajouter_mg_maillage(mai);
439     for (int i=0;i<lst.size();i++)
440     {
441     double xyz[3];
442     lst[i]->get_coord3(xyz);
443     MG_NOEUD* no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],TRIANGULATION);
444     mai->ajouter_mg_noeud(no);
445     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);}
446     }
447    
448     gestemp.enregistrer("debug.magic");
449     */
450     // fin temp
451     if (compteur!=NULL) compteur->ajouter_etape("Echantillonnage");
452     if (affichageactif)
453     {
454     char mess[300];
455     sprintf(mess," Preparation newton");
456     affiche(mess);
457     }
458     calcullevelsetpremierepasse(face,orientation,&lst,0,lst.size());
459     if (compteur!=NULL) compteur->ajouter_etape("Preparation newton");
460    
461     if (affichageactif)
462     {
463     char mess[300];
464     sprintf(mess," Newton");
465     affiche(mess);
466     }
467     calcullevelsetdeuxiemepasse(&lst);
468     if (compteur!=NULL) compteur->ajouter_etape("Newton");
469     if (affichageactif)
470     {
471     char mess[300];
472     sprintf(mess," Remplissage des trous");
473     affiche(mess);
474     }
475     remplir_trou(&lst,face,orientation);
476     if (compteur!=NULL) compteur->ajouter_etape("Remplissage des trous");
477     //etendrelevelset(solution,numsol);
478     int i=0;
479     for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
480     {
481 francois 280 if (noeud->get_etat()==1) solution->ecrire(i,numsol,noeud->get_solution());
482 francois 276 else solution->ecrire(i,numsol,1e300);
483     ++i;
484     }
485     i=0;
486     for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra))
487     {
488 francois 280 if (tet->get_etat()==1) solution_ele->ecrire(i,numsol,1.);
489 francois 276 ++i;
490     }
491     int nbpt=lst.size();
492     for (int i=0;i<nbpt;i++) delete lst[i];
493     TOIBREP_POINT::remisecompteurid();
494    
495    
496    
497     }
498    
499     void TOIBREP::remplir_trou(std::vector<TOIBREP_POINT*> *lst,MG_FACE* face,int sens)
500     {
501     int castot=0;
502     int cas=0;
503 francois 280 int couche_active;
504     int nb_couche=0;
505     do
506     {
507     nb_couche++;
508 francois 276 LISTE_FEM_TETRA::iterator ittetra;
509 francois 280 std::vector<FEM_TETRA*> couche;
510     couche_active=0;
511     for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra))
512 francois 276 {
513 francois 280 if (tet->get_etat()!=0) continue;
514 francois 276 FEM_NOEUD *no1=tet->get_fem_noeud(0);
515     FEM_NOEUD *no2=tet->get_fem_noeud(1);
516     FEM_NOEUD *no3=tet->get_fem_noeud(2);
517     FEM_NOEUD *no4=tet->get_fem_noeud(3);
518     int nbactif=0;
519     int nbpositif=0;
520 francois 280 if ((no1->get_etat()==1) || (no1->get_etat()==-1))
521 francois 276 {
522     nbactif++;
523     if (no1->get_solution()>0.) nbpositif++;
524     }
525 francois 280 if ((no2->get_etat()==1) || (no2->get_etat()==-1))
526 francois 276 {
527     nbactif++;
528     if (no2->get_solution()>0.) nbpositif++;
529     }
530 francois 280 if ((no3->get_etat()==1) || (no3->get_etat()==-1))
531 francois 276 {
532     nbactif++;
533     if (no3->get_solution()>0.) nbpositif++;
534     }
535 francois 280 if ((no4->get_etat()==1) || (no4->get_etat()==-1))
536 francois 276 {
537     nbactif++;
538     if (no4->get_solution()>0.) nbpositif++;
539     }
540     if ((nbactif==3) && (nbpositif<3) && (nbpositif>0))
541     {
542     couche.push_back(tet);
543     }
544     }
545     for (int i=0;i<couche.size();i++)
546     {
547     FEM_TETRA* tet=couche[i];
548     FEM_NOEUD *no1=tet->get_fem_noeud(0);
549     FEM_NOEUD *no2=tet->get_fem_noeud(1);
550     FEM_NOEUD *no3=tet->get_fem_noeud(2);
551     FEM_NOEUD *no4=tet->get_fem_noeud(3);
552     TOIBREP_POINT* pt;
553     FEM_NOEUD* noeudref;
554 francois 280 if (no1->get_etat()==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];}
555     if (no2->get_etat()==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];}
556     if (no3->get_etat()==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];}
557     if (no4->get_etat()==0) {noeudref=no4;pt=(*lst)[no1->get_numero()];}
558 francois 276 MG_FACE* face=(MG_FACE*)pt->get_mg_element_topologique();
559 francois 280 double uv[2];
560 francois 276 double xyz[3];
561 francois 280 double dist=calcul_distance(noeudref,face,pt,xyz,uv);
562     double normal[3];
563 francois 276 face->calcul_normale_unitaire(uv,normal);
564     normal[0]=normal[0]*sens*(-1.);
565     normal[1]=normal[1]*sens*(-1.);
566     normal[2]=normal[2]*sens*(-1.);
567     double dist2=calculdist(normal,xyz[0],xyz[1],xyz[2],noeudref);
568     if (dist2<0.) dist=dist*(-1);
569     noeudref->change_solution(dist);
570 francois 280 noeudref->change_etat(2);
571     TOIBREP_POINT *nvpt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],1,face);
572     lst->push_back(nvpt);
573     noeudref->change_numero(nvpt->get_id());
574 francois 276
575     }
576     for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra))
577     {
578     FEM_NOEUD *no1=tet->get_fem_noeud(0);
579     FEM_NOEUD *no2=tet->get_fem_noeud(1);
580     FEM_NOEUD *no3=tet->get_fem_noeud(2);
581     FEM_NOEUD *no4=tet->get_fem_noeud(3);
582 francois 280 if (no1->get_etat()!=0)
583     if (no2->get_etat()!=0)
584     if (no3->get_etat()!=0)
585     if (no4->get_etat()!=0)
586     if (tet->get_etat()==0)
587 francois 276 {
588     castot++;
589 francois 280 double lstptdecoup[5*3];
590 francois 276 int nb=0;
591 francois 280 decoupe_tetra(tet,lstptdecoup,&nb);
592     int reactive=0;
593 francois 276 for (int i=0;i<nb;i++)
594     {
595 francois 280 if (reactive) continue;
596 francois 276 double uv[2];
597 francois 280 face->inverser(uv,lstptdecoup+3*i);
598 francois 276 if ((face->valide_parametre_u(uv[0])) && (face->valide_parametre_v(uv[1])))
599     {
600     double distance=ot.calcule_distance_contour_face_uv(uv,face);
601     if (distance>1e-6)
602     {
603     cas++;
604     tet->change_etat(1);
605 francois 280 no1->change_etat(1);
606     no2->change_etat(1);
607     no3->change_etat(1);
608     no4->change_etat(1);
609     reactive=1;
610     couche_active=1;
611 francois 276 }
612     }
613     }
614 francois 280 if (!reactive)
615     {
616     if (no1->get_etat()==2) no1->change_etat(-1);
617     if (no2->get_etat()==2) no2->change_etat(-1);
618     if (no3->get_etat()==2) no3->change_etat(-1);
619     if (no4->get_etat()==2) no4->change_etat(-1);
620     tet->change_etat(-1);
621     }
622 francois 276 }
623     }
624 francois 280 }
625     while (couche_active==1);
626 francois 276 double res=cas*100.0/castot;
627     char mess[255];
628 francois 280 sprintf(mess," %lf%% en %d couches",res,nb_couche);
629 francois 276 affiche(mess);
630     }
631    
632     void TOIBREP::decoupe_tetra(FEM_TETRA* tet,double* lst,int *nb)
633     {
634     FEM_NOEUD* no1=tet->get_fem_noeud(0);
635     FEM_NOEUD* no2=tet->get_fem_noeud(1);
636     FEM_NOEUD* no3=tet->get_fem_noeud(2);
637     FEM_NOEUD* no4=tet->get_fem_noeud(3);
638     if (no1->get_solution()*no2->get_solution()<-0.00000001)
639     {
640     double t=no1->get_solution()/(no1->get_solution()-no2->get_solution());
641     double x=no1->get_x()+t*(no2->get_x()-no1->get_x());
642     double y=no1->get_y()+t*(no2->get_y()-no1->get_y());
643     double z=no1->get_z()+t*(no2->get_z()-no1->get_z());
644     lst[3*(*nb)]=x;
645     lst[3*(*nb)+1]=y;
646     lst[3*(*nb)+2]=z;
647     (*nb)++;
648     }
649     if (no1->get_solution()*no3->get_solution()<-0.00000001)
650     {
651     double t=no1->get_solution()/(no1->get_solution()-no3->get_solution());
652     double x=no1->get_x()+t*(no3->get_x()-no1->get_x());
653     double y=no1->get_y()+t*(no3->get_y()-no1->get_y());
654     double z=no1->get_z()+t*(no3->get_z()-no1->get_z());
655     lst[3*(*nb)]=x;
656     lst[3*(*nb)+1]=y;
657     lst[3*(*nb)+2]=z;
658     (*nb)++;
659     }
660     if (no1->get_solution()*no4->get_solution()<-0.00000001)
661     {
662     double t=no1->get_solution()/(no1->get_solution()-no4->get_solution());
663     double x=no1->get_x()+t*(no4->get_x()-no1->get_x());
664     double y=no1->get_y()+t*(no4->get_y()-no1->get_y());
665     double z=no1->get_z()+t*(no4->get_z()-no1->get_z());
666     lst[3*(*nb)]=x;
667     lst[3*(*nb)+1]=y;
668     lst[3*(*nb)+2]=z;
669     (*nb)++;
670     }
671     if (no2->get_solution()*no3->get_solution()<-0.00000001)
672     {
673     double t=no2->get_solution()/(no2->get_solution()-no3->get_solution());
674     double x=no2->get_x()+t*(no3->get_x()-no2->get_x());
675     double y=no2->get_y()+t*(no3->get_y()-no2->get_y());
676     double z=no2->get_z()+t*(no3->get_z()-no2->get_z());
677     lst[3*(*nb)]=x;
678     lst[3*(*nb)+1]=y;
679     lst[3*(*nb)+2]=z;
680     (*nb)++;
681     }
682     if (no2->get_solution()*no4->get_solution()<-0.00000001)
683     {
684     double t=no2->get_solution()/(no2->get_solution()-no4->get_solution());
685     double x=no2->get_x()+t*(no4->get_x()-no2->get_x());
686     double y=no2->get_y()+t*(no4->get_y()-no2->get_y());
687     double z=no2->get_z()+t*(no4->get_z()-no2->get_z());
688     lst[3*(*nb)]=x;
689     lst[3*(*nb)+1]=y;
690     lst[3*(*nb)+2]=z;
691     (*nb)++;
692     }
693     if (no3->get_solution()*no4->get_solution()<-0.00000001)
694     {
695     double t=no3->get_solution()/(no3->get_solution()-no4->get_solution());
696     double x=no3->get_x()+t*(no4->get_x()-no3->get_x());
697     double y=no3->get_y()+t*(no4->get_y()-no3->get_y());
698     double z=no3->get_z()+t*(no4->get_z()-no3->get_z());
699     lst[3*(*nb)]=x;
700     lst[3*(*nb)+1]=y;
701     lst[3*(*nb)+2]=z;
702     (*nb)++;
703     }
704     }
705    
706     void TOIBREP::calcullevelsetpremierepasse(MG_FACE* face,int sens,vector<TOIBREP_POINT*> *lst,int n1,int n2)
707     {
708     /*LISTE_FEM_TETRA::iterator it3;
709     for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3))
710     {
711     double *xyz1,*xyz2,*xyz3,*xyz4;
712     if (tet->get_nb_fem_noeud()==4)
713     {
714     xyz1=tet->get_fem_noeud(0)->get_coord();
715     xyz2=tet->get_fem_noeud(1)->get_coord();
716     xyz3=tet->get_fem_noeud(2)->get_coord();
717     xyz4=tet->get_fem_noeud(3)->get_coord();
718     }
719     if (tet->get_nb_fem_noeud()==10)
720     {
721     xyz1=tet->get_fem_noeud(0)->get_coord();
722     xyz2=tet->get_fem_noeud(2)->get_coord();
723     xyz3=tet->get_fem_noeud(4)->get_coord();
724     xyz4=tet->get_fem_noeud(9)->get_coord();
725     }
726     double xcentre=0.25*(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0]);
727     double ycentre=0.25*(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1]);
728     double zcentre=0.25*(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2]);
729     double rayon1=(xyz1[0]-xcentre)*(xyz1[0]-xcentre)+(xyz1[1]-ycentre)*(xyz1[1]-ycentre)+(xyz1[2]-zcentre)*(xyz1[2]-zcentre);
730     double rayon2=(xyz2[0]-xcentre)*(xyz2[0]-xcentre)+(xyz2[1]-ycentre)*(xyz2[1]-ycentre)+(xyz2[2]-zcentre)*(xyz2[2]-zcentre);
731     double rayon3=(xyz3[0]-xcentre)*(xyz3[0]-xcentre)+(xyz3[1]-ycentre)*(xyz3[1]-ycentre)+(xyz3[2]-zcentre)*(xyz3[2]-zcentre);
732     double rayon4=(xyz4[0]-xcentre)*(xyz4[0]-xcentre)+(xyz4[1]-ycentre)*(xyz4[1]-ycentre)+(xyz4[2]-zcentre)*(xyz4[2]-zcentre);
733     double rayon=std::max(rayonETRA_ACTIF1,std::max(rayon2,std::max(rayon3,rayon4)));
734     rayon=sqrt(rayon);
735     TPL_MAP_ENTITE<TOIBREP_POINT*> liste;
736     octree.rechercher(xcentre,ycentre,zcentre,rayon,liste);
737     TPL_MAP_ENTITE<TOIBREP_POINT*>::ITERATEUR it;
738     for (TOIBREP_POINT* pt=liste.get_premier(it);pt!=NULL;pt=liste.get_suivant(it))
739     {
740     double xyzpt[3];
741     pt->get_coord3(xyzpt);
742     int etat=estdansletetra(tet,xyzpt[0],xyzpt[1],xyzpt[2]);
743     if (compare_etat(etat,INTERIEUR))
744     {
745     for (int k=0;k<tet->get_nb_fem_noeud();k++)
746     {
747     double normal[3];
748     double uv[2];
749     face->inverser(uv,xyzpt);
750     face->calcul_normale_unitaire(uv,normal);
751     normal[0]=normal[0]*sens*(-1.);
752     normal[1]=normal[1]*sens*(-1.);
753     normal[2]=normal[2]*sens*(-1.);
754     double dist=calculdist(normal,xyzpt[0],xyzpt[1],xyzpt[2],tet->get_fem_noeud(k));
755     if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution()))
756     {
757     tet->get_fem_noeud(k)->change_solution(dist);
758     tet->get_fem_noeud(k)->change_numero(pt->get_id());
759     pt->change_coord2(uv);
760     }
761     }
762     }
763     if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur()))
764     {
765     tet->change_etat(1);
766     for (int k=0;k<tet->get_nb_fem_noeud();k++)
767     tet->get_fem_noeud(k)->change_etat(1);
768     }
769     if (tet->get_etat()==1) break;
770     }
771    
772     }
773     octree.vide();
774     }*/
775     /*LISTE_FEM_TETRA::iterator it3;
776     for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3))
777     octree.inserer(tet);*/
778     for (int i=n1;i<n2;i++)
779     {
780 francois 259 TOIBREP_POINT* pt=(*lst)[i];
781 francois 104 double uv[2];
782 francois 106 double xyz[3];
783     pt->get_coord3(xyz);
784 francois 222 double normal[3];
785 francois 104 face->inverser(uv,xyz);
786     face->calcul_normale_unitaire(uv,normal);
787 francois 222 normal[0]=normal[0]*sens*(-1.);
788     normal[1]=normal[1]*sens*(-1.);
789     normal[2]=normal[2]*sens*(-1.);
790     TPL_MAP_ENTITE<FEM_TETRA*> liste;
791 francois 259 octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste);
792 francois 222 TPL_MAP_ENTITE<FEM_TETRA*>::ITERATEUR it;
793     for (FEM_TETRA* tet=liste.get_premier(it);tet!=NULL;tet=liste.get_suivant(it))
794 francois 104 {
795 francois 276 if (tet->get_etat()==1) continue;
796     int etat=estdansletetra(tet,xyz[0],xyz[1],xyz[2]);
797     if (compare_etat(etat,INTERIEUR))
798 francois 104 {
799 francois 222 for (int k=0;k<tet->get_nb_fem_noeud();k++)
800 francois 106 {
801 francois 222 double dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_fem_noeud(k));
802     if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution()))
803     {
804     tet->get_fem_noeud(k)->change_solution(dist);
805     tet->get_fem_noeud(k)->change_numero(pt->get_id());
806     pt->change_coord2(uv);
807     }
808 francois 106 }
809 francois 104 }
810 francois 276 if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur()))
811     {
812     tet->change_etat(1);
813     for (int k=0;k<tet->get_nb_fem_noeud();k++)
814     {
815     tet->get_fem_noeud(k)->change_etat(1);
816     //octree.supprimer(tet);
817     }
818     }
819 francois 104 }
820     }
821 francois 276 //octree.vide();
822    
823    
824    
825    
826 francois 106 }
827    
828 francois 276 void TOIBREP::calcullevelsetdeuxiemepasse(vector<TOIBREP_POINT*> *lst)
829 francois 106 {
830 francois 276 LISTE_FEM_NOEUD::iterator itno;
831     for (FEM_NOEUD* no=mai->get_premier_noeud(itno);no!=NULL;no=mai->get_suivant_noeud(itno))
832     if (no->get_solution()<1e250)
833 francois 106 {
834 francois 276 TOIBREP_POINT* pt=(*lst)[no->get_numero()];
835 francois 106 int dim=pt->get_mg_element_topologique()->get_dimension();
836     if (dim==2)
837     {
838     MG_FACE* face=(MG_FACE*)pt->get_mg_element_topologique();
839 francois 280 double xyz[3],uv[2];
840     double dist=calcul_distance(no,face,pt,xyz,uv);
841 francois 106 int signe=1;
842 francois 276 if (no->get_solution()<0) signe=-1;
843     no->change_solution(signe*dist);
844 francois 106 }
845     if (dim==1)
846     {
847     MG_ARETE* arete=(MG_ARETE*)pt->get_mg_element_topologique();
848 francois 276 double dist=calcul_distance(no,arete,pt);
849     no->change_solution(dist);
850 francois 106 }
851     if (dim==0)
852     {
853     MG_SOMMET* sommet=(MG_SOMMET*)pt->get_mg_element_topologique();
854     MG_ARETE* arete=sommet->get_mg_cosommet(0)->get_arete();
855 francois 276 double dist=calcul_distance(no,arete,pt);
856     no->change_solution(dist);
857 francois 106 }
858    
859     }
860 francois 104
861 francois 276
862    
863    
864 francois 106 }
865    
866 francois 259 double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_ARETE* are,TOIBREP_POINT* pt,double precision)
867 francois 106 {
868     double tii,eps;
869     are->inverser(tii,noeud->get_coord());
870     int compteur=0;
871     OT_VECTEUR_3D Pt(noeud->get_x(),noeud->get_y(),noeud->get_z());
872     do
873     {
874     compteur++;
875     double ti=tii;
876     double xyz[3],dxyz[3],ddxyz[3];
877     are->deriver_seconde(ti,ddxyz,dxyz,xyz);
878     OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]);
879     OT_VECTEUR_3D Ct_deriver(dxyz[0],dxyz[1],dxyz[2]);
880     OT_VECTEUR_3D Ct_deriver_seconde(ddxyz[0],ddxyz[1],ddxyz[2]);
881     OT_VECTEUR_3D Distance = Ct-Pt;
882     tii=ti-Ct_deriver*Distance/(Ct_deriver_seconde*Distance+Ct_deriver.get_longueur2());
883     eps=fabs(tii-ti);
884     if (compteur>500) return 1e300;
885     if (tii<are->get_tmin())
886     {
887     tii=are->get_tmin();
888     eps=0.;
889     }
890     if (tii>are->get_tmax())
891     {
892     tii=are->get_tmax();
893     eps=0.;
894     }
895     }
896     while (eps>precision);
897     double xyz[3],dxyz[3];
898     are->evaluer(tii,xyz);
899     OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]);
900     double distance=(Ct-Pt).get_longueur();
901     MG_FACE* face1=are->get_mg_coarete(0)->get_boucle()->get_mg_face();
902     MG_FACE* face2=are->get_mg_coarete(1)->get_boucle()->get_mg_face();
903     are->deriver(tii,dxyz);
904     OT_VECTEUR_3D x1(dxyz[0],dxyz[1],dxyz[2]);
905     OT_VECTEUR_3D x2(dxyz[0],dxyz[1],dxyz[2]);
906     x1=are->get_mg_coarete(0)->get_orientation()*x1;
907     x2=are->get_mg_coarete(1)->get_orientation()*x2;
908     x1.norme();
909     x2.norme();
910     double uv1[2],uv2[2];
911     double normal1[3],normal2[3];
912     face1->inverser(uv1,xyz);
913     face2->inverser(uv2,xyz);
914     face1->calcul_normale_unitaire(uv1,normal1);
915     face2->calcul_normale_unitaire(uv2,normal2);
916     OT_VECTEUR_3D z1(normal1[0],normal1[1],normal1[2]);
917     OT_VECTEUR_3D z2(normal2[0],normal2[1],normal2[2]);
918     z1=face1->get_mg_coface(0)->get_orientation()*z1;
919     z2=face2->get_mg_coface(0)->get_orientation()*z2;
920     OT_VECTEUR_3D y1=z1&x1;
921     OT_VECTEUR_3D y2=z2&x2;
922     double test=(z1&z2)*x1;
923     int signe=-1;
924     OT_VECTEUR_3D dirpt=Pt-Ct;
925     if (test>0)
926     {if ((z1*dirpt>0.) || (z2*dirpt>0.)) signe=1;}
927     else
928     {if ((z1*dirpt>0.) && (z2*dirpt>0.)) signe=1;}
929     return signe*distance;
930     }
931    
932 francois 280 double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_FACE* face,TOIBREP_POINT* pt,double *xyz,double *uv,double precision)
933 francois 106 {
934     double uvii[2],eps;
935     pt->get_coord2(uvii);
936     int compteur=0;
937     OT_VECTEUR_3D Pt(noeud->get_x(),noeud->get_y(),noeud->get_z());
938     double delta_u,delta_v;
939     do
940     {
941     compteur++;
942     double uvi[2];
943     uvi[0]=uvii[0];
944     uvi[1]=uvii[1];
945     double xyzduu[3],xyzdvv[3],xyzduv[3],xyzdu[3],xyzdv[3],xyz[3];
946     face->deriver_seconde(uvi,xyzduu,xyzduv,xyzdvv,xyz,xyzdu,xyzdv);
947     OT_VECTEUR_3D S(xyz[0],xyz[1],xyz[2]);
948     OT_VECTEUR_3D Su(xyzdu[0],xyzdu[1],xyzdu[2]);
949     OT_VECTEUR_3D Sv(xyzdv[0],xyzdv[1],xyzdv[2]);
950     OT_VECTEUR_3D Suu(xyzduu[0],xyzduu[1],xyzduu[2]);
951     OT_VECTEUR_3D Suv(xyzduv[0],xyzduv[1],xyzduv[2]);
952     OT_VECTEUR_3D Svv(xyzdvv[0],xyzdvv[1],xyzdvv[2]);
953     OT_VECTEUR_3D Distance = S-Pt;
954     double a[4],b[2];
955     a[0]=Su.get_longueur2()+Distance*Suu;
956     a[1]=Su*Sv+Distance*Suv;
957     a[2]=Su*Sv+Distance*Suv;
958     a[3]=Sv.get_longueur2()+Distance*Svv;
959     b[0]=Distance*Su;b[0]=-b[0];
960     b[1]=Distance*Sv;b[1]=-b[1];
961     double deltau,deltav;
962     double denominateur_delta=(a[0]*a[3]-a[2]*a[1]);
963     if (a[0]<1E-12)
964     deltau=0;
965     else delta_u=(b[0]*a[3]-b[1]*a[1])/denominateur_delta;
966     if (a[3]<1E-12)
967     deltav=0;
968     else delta_v=(a[0]*b[1]-a[2]*b[0])/denominateur_delta;
969     /*if (fabs(denominateur_delta) < ( (fabs(a[0])+fabs(a[1])+fabs(a[2])+fabs(a[3]))*1e-12 ) )
970     return 1e300;*/
971     uvii[0]=uvi[0]+delta_u;
972     uvii[1]=uvi[1]+delta_v;
973     if (face->get_surface()->est_periodique_u()==1)
974     {
975     if(uvii[0]<0.) uvii[0]=face->get_surface()->get_periode_u()-uvii[0];
976     if(uvii[0]>face->get_surface()->get_periode_u()) uvii[0]=uvii[0]-face->get_surface()->get_periode_u();
977     }
978     if (face->get_surface()->est_periodique_v()==1)
979     {
980     if(uvii[1]<0.) uvii[0]=face->get_surface()->get_periode_v()-uvii[1];
981     if(uvii[1]>face->get_surface()->get_periode_v()) uvii[1]=uvii[1]-face->get_surface()->get_periode_v();
982     }
983     delta_u=uvii[0]-uvi[0];
984     delta_v=uvii[1]-uvi[1];
985     if (compteur>500) return 1e300;
986     }
987    
988     while ((fabs(delta_u)>precision)||(fabs(delta_v)>precision));
989     face->evaluer(uvii,xyz);
990     OT_VECTEUR_3D S(xyz[0],xyz[1],xyz[2]);
991     double distance=(S-Pt).get_longueur();
992 francois 280 uv[0]=uvii[0];
993     uv[1]=uvii[1];
994 francois 106 return distance;
995     }
996    
997    
998 francois 276
999 francois 259 void TOIBREP::etendrelevelset(FEM_SOLUTION* sol,int numsol)
1000 francois 222 {
1001     sol->active_solution(numsol);
1002    
1003    
1004     LISTE_FM know;
1005     LISTE_FM trial;
1006     LISTE_FM far;
1007     LISTE_FM exterieur;
1008     LISTE_FM_TRI trialtri;
1009     LISTE_FM_TRI_ID trialtriid;
1010    
1011    
1012     LISTE_FEM_TETRA::iterator ittet;
1013     for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
1014     {
1015     tet->change_solution(0.);
1016     int numsol=0;
1017     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);
1018     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);
1019     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);
1020     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);
1021     if (numsol==4)
1022     ajouter_liste(know,tet);
1023    
1024     else if (numsol==3)
1025     {
1026     if (tet->get_fem_noeud(0)->get_numero()==0) tet->change_numero(0);
1027     else if (tet->get_fem_noeud(1)->get_numero()==0) tet->change_numero(1);
1028     else if (tet->get_fem_noeud(2)->get_numero()==0) tet->change_numero(2);
1029     else tet->change_numero(3);
1030     ajouter_liste(trial,tet);
1031     }
1032     else
1033     ajouter_liste(far,tet);
1034    
1035     }
1036     for (LISTE_FM::iterator i=trial.begin();i!=trial.end();i++)
1037     {
1038     FEM_TETRA* tet=*i;
1039     int signe;
1040     double sol=resoudgradT(tet,&signe);
1041     if (fabs(sol)>0.00000001)
1042     {
1043     if (fabs(sol)<fabs(tet->get_fem_noeud(tet->get_numero())->get_solution()))
1044     {
1045     tet->get_fem_noeud(tet->get_numero())->change_solution(sol);
1046     ajouter_liste(trialtri,trialtriid,tet,fabs(tet->get_fem_noeud(tet->get_numero())->get_solution()));
1047     }
1048     }
1049     else
1050     {
1051     ajouter_liste(exterieur,tet);
1052     tet->change_solution(1e300);
1053     }
1054     }
1055     int fin=0;
1056     LISTE_FM_TRI::iterator itfin=trialtri.end();
1057     itfin--;
1058     double longref=(*itfin).first;
1059     do
1060     {
1061     LISTE_FM_TRI::iterator it=trialtri.begin();
1062     FEM_TETRA* tet=(*it).second;
1063     double longcourant=(*it).first;
1064     supprimer_liste(trialtri,trialtriid,tet);
1065     ajouter_liste(know,tet);
1066     FEM_NOEUD* noeud=tet->get_fem_noeud(tet->get_numero());
1067     noeud->change_numero(1);
1068     if (noeud->get_solution()>20000)
1069     cout << "BUGGGGGGG" <<endl;
1070     int nbtetra=noeud->get_lien_tetra()->get_nb();
1071     for (int i=0;i<nbtetra;i++)
1072     {
1073     FEM_TETRA* tet2=noeud->get_lien_tetra()->get(i);
1074     if (tet2==tet) continue;
1075     LISTE_FM_TRI_ID::iterator it=trialtriid.find(tet2->get_id());
1076     if (it!=trialtriid.end())
1077     {
1078     int signe;
1079     double sol=resoudgradT(tet2,&signe);
1080     double solution=tet2->get_fem_noeud(tet2->get_numero())->get_solution();
1081     if (fabs(sol)>0.00000001)
1082     if (!((solution<1e200)&&(sol*solution<0)))
1083     if (fabs(sol)<fabs(solution))
1084     {
1085     supprimer_liste(trialtri,trialtriid,tet2);
1086     ajouter_liste(trialtri,trialtriid,tet2,fabs(sol));
1087     tet2->get_fem_noeud(tet2->get_numero())->change_solution(sol);
1088     }
1089     }
1090     LISTE_FM::iterator it2=find(far.begin(),far.end(),tet2);
1091     if (it2!=far.end())
1092     {
1093     int numsol=0;
1094     if (tet2->get_fem_noeud(0)->get_numero()==1) numsol++;
1095     if (tet2->get_fem_noeud(1)->get_numero()==1) numsol++;
1096     if (tet2->get_fem_noeud(2)->get_numero()==1) numsol++;
1097     if (tet2->get_fem_noeud(3)->get_numero()==1) numsol++;
1098     //if (numsol==4)
1099     //cout << " BUG " <<endl;
1100     if (numsol==3)
1101     {
1102     if (tet2->get_fem_noeud(0)->get_numero()==0) tet2->change_numero(0);
1103     else if (tet2->get_fem_noeud(1)->get_numero()==0) tet2->change_numero(1);
1104     else if (tet2->get_fem_noeud(2)->get_numero()==0) tet2->change_numero(2);
1105     else tet2->change_numero(3);
1106     int signe;
1107     double sol=resoudgradT(tet2,&signe);
1108     double ancsol=tet2->get_fem_noeud(tet2->get_numero())->get_solution();
1109     if (fabs(sol)>0.00000001)
1110     {
1111     if (!((ancsol<1e200) && (ancsol*sol<0.) ))
1112     {
1113     tet2->get_fem_noeud(tet2->get_numero())->change_solution(sol);
1114     supprimer_liste(far,tet2);
1115     ajouter_liste(trialtri,trialtriid,tet2,fabs(sol));
1116     }
1117     }
1118     else
1119     {
1120     tet2->change_solution(1e300);
1121     ajouter_liste(exterieur,tet2);
1122     }
1123     }
1124     }
1125     }
1126     if (trialtri.size()==0) fin=1;
1127     if (exterieur.size()>0)
1128     if (fin==0)
1129     if (longcourant>longref)
1130     {
1131     int nombre=exterieur.size();
1132     for (int i=0;i<nombre;i++)
1133     {
1134     FEM_TETRA* tet2=(*(exterieur.begin()));
1135     supprimer_liste(exterieur,tet2);
1136     ajouter_liste(know,tet2);
1137     for (int nd=0;nd<4;nd++)
1138     {
1139     FEM_NOEUD* noeud=tet2->get_fem_noeud(nd);
1140     int nbtetra=noeud->get_lien_tetra()->get_nb();
1141     for (int i=0;i<nbtetra;i++)
1142     {
1143     FEM_TETRA* tet3=noeud->get_lien_tetra()->get(i);
1144     if (tet2==tet3) continue;
1145     LISTE_FM::iterator it2=find(far.begin(),far.end(),tet3);
1146     if (it2!=far.end())
1147     {
1148     supprimer_liste(far,tet3);
1149     ajouter_liste(exterieur,tet3);
1150     tet3->change_solution(1e300);
1151     }
1152     }
1153     }
1154     }
1155     LISTE_FM_TRI::iterator itfin=trialtri.end();
1156     itfin--;
1157     longref=(*itfin).first;
1158     }
1159     }
1160     while (fin==0);
1161     LISTE_FEM_NOEUD::iterator itnoeud;
1162     for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
1163     noeud->change_numero(0);
1164     for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
1165     {
1166     if (tet->get_solution()>1e200) continue;
1167     tet->get_fem_noeud(0)->change_numero(1);
1168     tet->get_fem_noeud(1)->change_numero(1);
1169     tet->get_fem_noeud(2)->change_numero(1);
1170     tet->get_fem_noeud(3)->change_numero(1);
1171     }
1172     int i=0;
1173     for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
1174     {
1175     if (noeud->get_numero()==1) sol->ecrire(i,numsol,noeud->get_solution()); else sol->ecrire(i,numsol,1e300);
1176     ++i;
1177     }
1178     }
1179 francois 276 IBrep TOIBREP::exporter_IBrep(string chemin,FEM_SOLUTION* solution,FEM_SOLUTION *solution_ele,MG_GROUPE_TOPOLOGIQUE* mggt)
1180 francois 259 {
1181     TPL_MAP_ENTITE<MG_VOLUME*> lst;
1182     if (mggt==NULL)
1183     {
1184     int nbvol=geo->get_nb_mg_volume();
1185     for (int i=0;i<nbvol;i++)
1186     lst.ajouter(geo->get_mg_volume(i));
1187     }
1188     else
1189     {
1190     int nbmggt=mggt->get_nb();
1191     for (int i=0;i<nbmggt;i++)
1192     if (mggt->get(i)->get_dimension()==3)
1193     lst.ajouter((MG_VOLUME*)mggt->get(i));
1194     }
1195     IBrep brep(100);
1196     int nbvolexp=lst.get_nb();
1197     for (int i=0;i<nbvolexp;i++)
1198     {
1199     MG_VOLUME* vol=lst.get(i);
1200     int nbcoq=vol->get_nb_mg_coquille();
1201     IVolumeN ivoltmp(vol->get_id());
1202 francois 263 pair<IVolume*,bool> pair_ivol=brep.AddVolume(ivoltmp);
1203     IVolume* ivol=pair_ivol.first;
1204 francois 259 for (int j=0;j<nbcoq;j++)
1205     {
1206     MG_COQUILLE* coq=vol->get_mg_coquille(j);
1207     int nbface=coq->get_nb_mg_coface();
1208     IShell ishell(nbface);
1209     for (int k=0;k<nbface;k++)
1210     {
1211     MG_COFACE* coface=coq->get_mg_coface(k);
1212     ishell[k]=coface->get_id();
1213     MG_FACE* face=coface->get_face();
1214     int sens=coface->get_orientation();
1215     ICoFaceN icoface(coface->get_id(),face->get_id(),sens);
1216     brep.AddCoFace(icoface);
1217     IFace* iface=brep.GetFace(face->get_id());
1218     if (iface==NULL)
1219     {
1220     IFaceN ifacenew(face->get_id());
1221 francois 263 pair<IFace*,bool> pair_iface=brep.AddFace(ifacenew);
1222     iface=pair_iface.first;
1223 francois 259 int nbbou=face->get_nb_mg_boucle();
1224     for (int l=0;l<nbbou;l++)
1225     {
1226     MG_BOUCLE* bou=face->get_mg_boucle(l);
1227     int nbare=bou->get_nb_mg_coarete();
1228     ILoop iloop(nbare);
1229     for (int m=0;m<nbare;m++)
1230     {
1231     MG_COARETE* coare=bou->get_mg_coarete(m);
1232     iloop[m]=coare->get_id();
1233     MG_ARETE* are=coare->get_arete();
1234     int sens=coare->get_orientation();
1235     ICoEdgeN icoedge(coare->get_id(),are->get_id(),sens,face->get_id());
1236     brep.AddCoEdge(icoedge);
1237     IEdge* iedge=brep.GetEdge(are->get_id());
1238     IVertex *ivertex1,*ivertex2;
1239     ICoVertex *icover1,*icover2;
1240     if (iedge==NULL)
1241     {
1242     MG_COSOMMET* cover1=are->get_cosommet1();
1243     MG_COSOMMET* cover2=are->get_cosommet2();
1244     MG_SOMMET* ver1=cover1->get_sommet();
1245     MG_SOMMET* ver2=cover2->get_sommet();
1246     ICoVertexN icovertmp1(cover1->get_id(),are->get_id(),ver1->get_id(),cover1->get_t());
1247     ICoVertexN icovertmp2(cover2->get_id(),are->get_id(),ver2->get_id(),cover2->get_t());
1248 francois 263 pair<ICoVertex*,bool> pair_icover1=brep.AddCoVertex(icovertmp1);
1249     pair<ICoVertex*,bool> pair_icover2=brep.AddCoVertex(icovertmp2);
1250     icover1=pair_icover1.first;
1251     icover2=pair_icover2.first;
1252 francois 259 IEdgeN iedgetmp(are->get_id(),cover1->get_id(),cover2->get_id());
1253     ivertex1=brep.GetVertex(ver1->get_id());
1254     ivertex2=brep.GetVertex(ver2->get_id());
1255     if (ivertex1==NULL)
1256     {
1257     MG_POINT* pt=ver1->get_point();
1258     double xyz[3];
1259     pt->evaluer(xyz);
1260     IVertexN ivertex(ver1->get_id(),xyz);
1261 francois 263 pair<IVertex*,bool> pair_ivertex1=brep.AddVertex(ivertex);
1262     ivertex1=pair_ivertex1.first;
1263 francois 259 }
1264     if (ivertex2==NULL)
1265     {
1266     MG_POINT* pt=ver2->get_point();
1267     double xyz[3];
1268     pt->evaluer(xyz);
1269     IVertexN ivertex(ver2->get_id(),xyz);
1270 francois 263 pair<IVertex*,bool> pair_ivertex2=brep.AddVertex(ivertex);
1271     ivertex2=pair_ivertex2.first;
1272 francois 259 }
1273     ivertex1->AddCoVertex(cover1->get_id());
1274     ivertex2->AddCoVertex(cover2->get_id());
1275 francois 263 pair<IEdge*,bool> pair_iedge=brep.AddEdge(iedgetmp);
1276     iedge=pair_iedge.first;
1277 francois 259 }
1278     iedge->AddCoEdge(coare->get_id(),brep);
1279     if (sens==1)
1280     {
1281     icover1=brep.GetCoVertex(iedge->FromCoVertex());
1282     ivertex1=brep.GetVertex(icover1->Vertex());
1283     ivertex1->AddFace(face->get_id());
1284     }
1285     else
1286     {
1287     icover2=brep.GetCoVertex(iedge->ToCoVertex());
1288     ivertex2=brep.GetVertex(icover2->Vertex());
1289     ivertex2->AddFace(face->get_id());
1290     }
1291     }
1292     iface->AddLoop(iloop);
1293     }
1294     }
1295     iface->AddCoFace(coface->get_id());
1296     }
1297     ivol->AddShell(ishell);
1298     }
1299     }
1300     int nbsol=solution->get_nb_champ();
1301 francois 262 unsigned long *correspondface=new unsigned long[nbsol];
1302     for (int i=0;i<nbsol;i++)
1303     {
1304     std::string nom=solution->get_legende(i);
1305     char nom2[2];
1306     unsigned long id;
1307 francois 272 sscanf(nom.c_str(),"%s %lu",nom2,&id);
1308 francois 262 correspondface[i]=id;
1309     }
1310 francois 259 LISTE_FEM_NOEUD::iterator itnoeud;
1311     int i=0;
1312     for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud))
1313     {
1314     int nbsolnoeud=0;
1315     for (int j=0;j<nbsol;j++)
1316     {
1317     if (solution->lire(i,j)<1e200) nbsolnoeud++;
1318     }
1319     double *xyz=noeud->get_coord();
1320     INodeN inoeudtmp(noeud->get_id(),xyz[0],xyz[1],xyz[2],nbsolnoeud);
1321 francois 263 pair<INode*,bool> pair_inoeud=brep.AddNode(inoeudtmp);
1322     INode* inoeud=pair_inoeud.first;
1323 francois 259 int num=0;
1324     for (int j=0;j<nbsol;j++)
1325     {
1326     if (solution->lire(i,j)<1e200)
1327     {
1328 francois 262 inoeud->Fields()[num]=correspondface[j];
1329 francois 259 inoeud->Values()[num]=solution->lire(i,j);
1330     num++;
1331     }
1332     }
1333     i++;
1334     }
1335     LISTE_FEM_TETRA::iterator ittet;
1336     int num=0;
1337 francois 276 i=0;
1338 francois 259 for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
1339     {
1340 francois 276 int nbsolele=0;
1341     for (int j=0;j<nbsol;j++)
1342     {
1343     if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001) nbsolele++;
1344     }
1345     IElementN xfemele(IElement::TETRAHEDRON,nbsolele);
1346 francois 259 num++;
1347     xfemele.num=num;
1348     xfemele.GetNode(0)=tet->get_fem_noeud(0)->get_id();
1349     xfemele.GetNode(1)=tet->get_fem_noeud(1)->get_id();
1350     xfemele.GetNode(2)=tet->get_fem_noeud(2)->get_id();
1351     xfemele.GetNode(3)=tet->get_fem_noeud(3)->get_id();
1352 francois 276 int num=0;
1353     for (int j=0;j<nbsol;j++)
1354     {
1355     if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001)
1356     {
1357     xfemele.Fields()[num]=correspondface[j];
1358     num++;
1359     }
1360     }
1361 francois 259 brep.AddElement(xfemele);
1362 francois 276 i++;
1363 francois 259 }
1364 francois 222
1365 francois 276 delete [] correspondface;
1366 francois 259
1367    
1368    
1369     std::ofstream output(chemin.c_str());
1370     output << brep << std::endl;
1371     output.close();
1372    
1373 francois 262
1374    
1375    
1376 francois 276 return brep;
1377 francois 259 }
1378    
1379    
1380    
1381     void TOIBREP::ajouter_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_TETRA* tet,double val)
1382 francois 222 {
1383     pair<double,FEM_TETRA*> p(val,tet);
1384     LISTE_FM_TRI::iterator it=lst.insert(p);
1385     lstid[tet->get_id()]=it;
1386     }
1387    
1388 francois 259 void TOIBREP::supprimer_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_TETRA* tet)
1389 francois 222 {
1390     LISTE_FM_TRI::iterator it2=lstid[tet->get_id()];
1391     LISTE_FM_TRI_ID::iterator it3=lstid.find(tet->get_id());
1392     lstid.erase(it3);
1393     lst.erase(it2);
1394     }
1395    
1396 francois 259 void TOIBREP::ajouter_liste(LISTE_FM& lst,FEM_TETRA* tet)
1397 francois 222 {
1398     lst.push_back(tet);
1399     }
1400    
1401 francois 259 void TOIBREP::supprimer_liste(LISTE_FM& lst,FEM_TETRA* tet)
1402 francois 222 {
1403     LISTE_FM::iterator it=find(lst.begin(),lst.end(),tet);
1404     if (it!=lst.end()) lst.erase(it);
1405     }
1406    
1407    
1408    
1409 francois 259 double TOIBREP::resoudgradT(FEM_TETRA* tet,int *signe)
1410 francois 222 {
1411     double j[9];
1412     double N[12];
1413     double jN[12]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
1414     double T[4];
1415     double uv[2]={0.25,0.25};
1416    
1417    
1418     tet->get_inverse_jacob(j,uv);
1419     for (int i=0;i<3;i++)
1420     for (int k=0;k<4;k++)
1421     N[i*4+k]=tet->get_fonction_derive_interpolation(k+1,i+1,uv);
1422     int premier=0;
1423     double tmin=1e300;
1424     double tmax=-1e300;
1425     for (int i=0;i<4;i++)
1426     {
1427     if (i==tet->get_numero()) continue;
1428     T[i]=tet->get_fem_noeud(i)->get_solution();
1429     if (fabs(T[i])>0.000001)
1430     {
1431     if (premier==0)
1432     if (T[i]>0) (*signe)=1; else (*signe)=-1;
1433     else if (T[i]*(*signe)<0) (*signe)=0;
1434     }
1435     T[i]=fabs(T[i]);
1436     if (tet->get_numero()!=i)
1437     {
1438     if (T[i]<tmin) tmin=T[i];
1439     if (T[i]>tmax) tmax=T[i];
1440     }
1441     premier=1;
1442     }
1443     for (int i=0;i<3;i++)
1444     for (int k=0;k<4;k++)
1445     for (int l=0;l<3;l++)
1446     jN[i*4+k]=jN[i*4+k]+j[i*3+l]*N[l*4+k];
1447     double a=0.,b=0.,c=-1.;
1448     for (int i=0;i<3;i++)
1449     {
1450     double coef=0.;
1451     double coefinc=0.;
1452     for (int l=0;l<4;l++)
1453     {
1454     if (tet->get_numero()!=l) coef=coef+jN[i*4+l]*T[l];
1455     else coefinc=coefinc+jN[i*4+l];
1456     }
1457     c=c+coef*coef;
1458     a=a+coefinc*coefinc;
1459     b=b+2*coef*coefinc;
1460     }
1461     /*if (*signe==0)
1462     cout << "attention " <<endl;*/
1463     double det=b*b-4.*a*c;
1464     if (det<0.) det=0.; else det=sqrt(det);
1465     double sol1=(-b-det)/2./a;
1466     double sol2=(-b+det)/2./a;
1467     double sol=sol1;
1468     if (sol2>sol1) sol=sol2;
1469     if (sol<tmin*0.99)
1470     sol=0.;
1471     sol=sol*(*signe);
1472     return sol;
1473     }