ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur2d_outil.cpp
Revision: 820
Committed: Fri Aug 19 19:45:21 2016 UTC (8 years, 8 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur2d_outil.cpp
File size: 18822 byte(s)
Log Message:
prise en compte des poles dans le maillage des spheres pour le mailleur 2D. 20 ans après ca marche un peu pres.

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4 francois 331 // Jean Christophe Cuilli�re et Vincent FRANCOIS
5     // D�partement de G�nie M�canique - UQTR
6 5 //------------------------------------------------------------
7 francois 331 // Le projet MAGIC est un projet de recherche du d�partement
8     // de g�nie m�canique de l'Universit� du Qu�bec �
9     // Trois Rivi�res
10     // Les librairies ne peuvent �tre utilis�es sans l'accord
11 5 // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur2d_outil.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20 francois 331 // Version du 02/03/2006 � 11H23
21 5 //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26     #include <math.h>
27     #include "mailleur2d.h"
28     #include "ot_mathematique.h"
29 francois 481 #include "ot_boite_2d.h"
30 5 #include "tpl_map_entite.h"
31     #include "m3d_triangle.h"
32    
33    
34    
35    
36    
37    
38     int MAILLEUR2D::noeud_est_dans_triangle(MG_NOEUD* noeud,MG_TRIANGLE *triangle)
39     {
40     return noeud_est_dans_triangle(noeud,triangle->get_noeud1(),triangle->get_noeud2(),triangle->get_noeud3());
41     }
42    
43    
44     int MAILLEUR2D::noeud_est_dans_triangle(MG_NOEUD* noeud,MG_NOEUD *noeud1,MG_NOEUD *noeud2,MG_NOEUD *noeud3)
45     {
46     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
47     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
48    
49     double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
50     double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
51     double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
52     double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
53     double u3=decalage->decalage_parametre_u(noeud3->get_u(),du);
54     double v3=decalage->decalage_parametre_v(noeud3->get_v(),dv);
55     double u=decalage->decalage_parametre_u(noeud->get_u(),du);
56     double v=decalage->decalage_parametre_v(noeud->get_v(),dv);
57     double eps=0.0001;
58     double delta=(u2-u1)*(v3-v1)-(v2-v1)*(u3-u1);
59     double precision=std::max(fabs(u1),fabs(v1));
60     precision=std::max(precision,fabs(u2));
61     precision=std::max(precision,fabs(v2));
62     precision=std::max(precision,fabs(u3));
63     precision=std::max(precision,fabs(v3));
64     precision=std::max(precision,fabs(u));
65     precision=std::max(precision,fabs(v));
66     if (OPERATEUR::egal(delta,0.0,precision*eps)) return(0);
67     double xsi=1.0/delta*((v3-v1)*(u-u1)-(u3-u1)*(v-v1));
68     double eta=1.0/delta*((u2-u1)*(v-v1)-(v2-v1)*(u-u1));
69     if (!((eta > eps) && (xsi > eps) && ((eta+xsi) < 1.0-eps))) return(0);
70     return (1);
71     }
72    
73    
74    
75    
76    
77    
78     int MAILLEUR2D::insere_segment(MG_FACE *mgface,MG_SEGMENT **nv_segment,MG_NOEUD* noeud1,MG_NOEUD* noeud2,int type_validation)
79     {
80    
81     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
82     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
83     double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
84     double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
85     double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
86     double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
87     double ui=0.5*(u1+u2);
88     double vi=0.5*(v1+v2);
89     double uii=decalage->decalage_parametre_u(ui,-du);
90     double vii=decalage->decalage_parametre_v(vi,-dv);
91     double longueur=sqrt((u2-u1)*(u2-u1)+(v2-v1)*(v2-v1));
92 francois 632 double xyz[3],uv[2]={uii,vii};
93     double xyz1[3],uv1[2]={noeud1->get_u(),noeud1->get_v()};
94     double xyz2[3],uv2[2]={noeud2->get_u(),noeud2->get_v()};
95     mgface->evaluer(uv1,xyz1);
96     mgface->evaluer(uv2,xyz2);
97     OT_VECTEUR_3D vec(xyz1,xyz2);
98     longueur=vec.get_longueur();
99     xyz[0]=0.5*(xyz1[0]+xyz2[0]);
100     xyz[1]=0.5*(xyz1[1]+xyz2[1]);
101     xyz[2]=0.5*(xyz1[2]+xyz2[2]);
102 5 if (type_validation==TOUS_SEGMENT)
103     {
104     TPL_MAP_ENTITE<MG_SEGMENT*> liste_trouvee;
105 francois 632 ntree_de_segment->rechercher(xyz[0],xyz[1],xyz[2],longueur/2.,liste_trouvee);
106 5 for (int j=0;j<liste_trouvee.get_nb();j++)
107     {
108     MG_SEGMENT* mgsegment=liste_trouvee.get(j);
109     if ( ((noeud1==mgsegment->get_noeud1()) && (noeud2==mgsegment->get_noeud2())) || ((noeud1==mgsegment->get_noeud2()) && (noeud2==mgsegment->get_noeud1())))
110     return(0);
111     int res=intersection_segment_segment(noeud1,noeud2,mgsegment->get_noeud1(),mgsegment->get_noeud2());
112     if (res!=0)
113     return(0);
114     }
115     }
116     if (type_validation==TOUS_FRONT)
117     {
118     TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
119 francois 632 ntree_de_front->rechercher(xyz[0],xyz[1],xyz[2],longueur/2.,liste_trouvee);
120 5 for (int j=0;j<liste_trouvee.get_nb();j++)
121     {
122     MG_FRONT_2D* ft=liste_trouvee.get(j);
123     MG_SEGMENT* mgsegment=ft->get_segment();
124     if ( ((noeud1==mgsegment->get_noeud1()) && (noeud2==mgsegment->get_noeud2())) || ((noeud1==mgsegment->get_noeud2()) && (noeud2==mgsegment->get_noeud1())))
125     return(0);
126     int res=intersection_segment_segment(noeud1,noeud2,mgsegment->get_noeud1(),mgsegment->get_noeud2());
127     if (res!=0)
128     return(0);
129     }
130     }
131 francois 820 if ((*nv_segment)==NULL)
132     {
133     MG_SEGMENT* mgsegment=mg_maillage->ajouter_mg_segment(mgface,noeud1,noeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
134     *nv_segment=mgsegment;
135     }
136     ntree_de_segment->inserer(*nv_segment);
137 5 return(1);
138     }
139    
140    
141     void MAILLEUR2D::supprime_segment(MG_SEGMENT* mgsegment)
142     {
143 francois 632 ntree_de_segment->supprimer(mgsegment);
144 5 mg_maillage->supprimer_mg_segmentid(mgsegment->get_id());
145     }
146    
147    
148     int MAILLEUR2D::genere_noeud(MG_FACE* mgface,MG_FRONT_2D* front,MG_FRONT_2D **front_rencontre,MG_NOEUD **noeud_rencontre)
149     {
150     OT_VECTEUR_3D w(0.,0.,1.);
151     MG_NOEUD* noeud1=front->get_noeud1();
152     MG_NOEUD* noeud2=front->get_noeud2();
153     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
154     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
155     double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
156     double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
157     double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
158     double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
159    
160 francois 532 double longueur_segment=metrique->calcule_longueur_segment_metrique(mgface,front->get_segment(),decalage,pas);
161 5 double longueur_desiree=0.8660254037844386*(MAILLEUR2D::priorite_metrique+longueur_segment-longueur_segment*MAILLEUR2D::priorite_metrique);
162     longueur_desiree=longueur_desiree/front->get_ifail();
163    
164    
165     double umilieu=0.5*(u1+u2);
166     double vmilieu=0.5*(v1+v2);
167     double param_milieu[2]={umilieu-du,vmilieu-dv};
168     double E,F,G;
169    
170     mgface->get_EFG(param_milieu,E,F,G);
171     double deno=E*(u2-u1)*(u2-u1)+G*(v2-v1)*(v2-v1)+2*F*(u2-u1)*(v2-v1);
172     double nume=(v2-v1)*(u2-u1)*(E-G)+F*((v2-v1)*(v2-v1)-(u2-u1)*(u2-u1));
173     double teta=atan(nume/deno);
174    
175    
176 francois 532 double udecale;
177     double vdecale;
178 francois 820 int res=metrique->ajuste_distance_ortho_metrique(mgface,u1,v1,u2,v2,udecale,vdecale,longueur_desiree,du,dv,teta,pas);
179 francois 632 if (res==0) return ERREUR;
180 5 double u=decalage->decalage_parametre_u(udecale,-du);
181     double v=decalage->decalage_parametre_v(vdecale,-dv);
182     double param_noeud_cree[2]={u,v};
183    
184    
185     // recherhce des element du front proche de u,v;
186     TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
187 francois 632 double xyz[3],uv[2]={u,v};
188     mgface->evaluer(uv,xyz);
189     double xmin=std::min(noeud1->get_x(),noeud2->get_x());xmin=std::min(xmin,xyz[0]);
190     double xmax=std::max(noeud1->get_x(),noeud2->get_x());xmax=std::max(xmax,xyz[0]);
191     double ymin=std::min(noeud1->get_y(),noeud2->get_y());ymin=std::min(ymin,xyz[1]);
192     double ymax=std::max(noeud1->get_y(),noeud2->get_y());ymax=std::max(ymax,xyz[1]);
193     double zmin=std::min(noeud1->get_z(),noeud2->get_z());zmin=std::min(zmin,xyz[2]);
194     double zmax=std::max(noeud1->get_z(),noeud2->get_z());zmax=std::max(zmax,xyz[2]);
195     BOITE_3D boite_recherche(xmin,ymin,zmin,xmax,ymax,zmax);
196 5 double longueur_recherche=boite_recherche.get_rayon()*2.;
197    
198 francois 632
199    
200     ntree_de_front->rechercher(xyz[0],xyz[1],xyz[2],longueur_recherche,liste_trouvee);
201 5 int nb_entite=liste_trouvee.get_nb();
202     double distance_reference = -1.;
203     double angle_reference=2.*M_PI;
204     unsigned long id_noeud_reference = 0;
205     MG_FRONT_2D* front_reference=NULL;
206    
207     for (int i=0;i<nb_entite;i++)
208     {
209     MG_FRONT_2D* front_courant=liste_trouvee.get(i);
210     MG_NOEUD* noeud_front1=front_courant->get_noeud1();
211     MG_NOEUD* noeud_front2=front_courant->get_noeud2();
212     double unoeudfront1=decalage->decalage_parametre_u(noeud_front1->get_u(),du);
213     double vnoeudfront1=decalage->decalage_parametre_v(noeud_front1->get_v(),dv);
214     double unoeudfront2=decalage->decalage_parametre_u(noeud_front2->get_u(),du);
215     double vnoeudfront2=decalage->decalage_parametre_v(noeud_front2->get_v(),dv);
216 francois 532 double distance_noeudfront1=metrique->calcule_distance_metrique(mgface,udecale,vdecale,unoeudfront1,vnoeudfront1,du,dv,pas);
217 5 if ((distance_noeudfront1<distance_reference) || (distance_reference<0.) || (id_noeud_reference==noeud_front1->get_id()))
218     if (noeud_front1->get_id()!=noeud1->get_id())
219     if (noeud_front1->get_id()!=noeud2->get_id())
220     if (id_noeud_reference==noeud_front1->get_id())
221     {
222     refresh();
223     OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
224     OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
225     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
226     vecteur_baseu.norme();
227     vecteur_basev.norme();
228     vecteur_front.norme();
229     double cosangle=vecteur_baseu*vecteur_front;
230     double sinangle=vecteur_basev*vecteur_front;
231     sinangle=-sinangle;
232     if (cosangle>1.) cosangle=1.;
233     if (cosangle<-1.) cosangle=(-1.);
234     double angle=acos(cosangle);
235     if (sinangle<(-0.0001)) angle=(-angle);
236     if (angle<0.) angle=angle+2.*M_PI;
237     if (angle<angle_reference)
238     {
239     angle_reference=angle;
240     front_reference=front_courant;
241     }
242     }
243     else
244     {
245     refresh();
246     OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
247     OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
248     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
249     vecteur_baseu.norme();
250     vecteur_basev.norme();
251     vecteur_front.norme();
252     double cosangle=vecteur_baseu*vecteur_front;
253     double sinangle=vecteur_basev*vecteur_front;
254     sinangle=-sinangle;
255     if (cosangle>1.) cosangle=1.;
256     if (cosangle<-1.) cosangle=(-1.);
257     double angle=acos(cosangle);
258     if (sinangle<(-0.0001)) angle=(-angle);
259     if (angle<0.) angle=angle+2.*M_PI;
260     distance_reference=distance_noeudfront1;
261     angle_reference=angle;
262     id_noeud_reference=noeud_front1->get_id();
263     front_reference=front_courant;
264     }
265     }
266    
267     if ((front_reference!=NULL) && (distance_reference<longueur_desiree*2./sqrt(3)))
268     {
269     MG_NOEUD* noeud=mg_maillage->get_mg_noeudid(id_noeud_reference);
270     double uref=decalage->decalage_parametre_u(noeud->get_u(),du);
271     double vref=decalage->decalage_parametre_v(noeud->get_v(),dv);
272 francois 532 double distance=metrique->calcule_distance_metrique(mgface,umilieu,vmilieu,uref,vref,du,dv,pas);
273 5 if (distance<1.5*longueur_desiree)
274     {
275     (*front_rencontre)=front_reference;
276     (*noeud_rencontre)=noeud;
277     return FRONT_RENCONTRE;
278     }
279     }
280 francois 791 MG_NOEUD* noeud=mg_maillage->ajouter_mg_noeud(mgface,xyz[0],xyz[1],xyz[2],MAGIC::ORIGINE::MAILLEUR_AUTO);
281 5 noeud->change_u(param_noeud_cree[0]);
282     noeud->change_v(param_noeud_cree[1]);
283     (*noeud_rencontre)=noeud;
284     (*front_rencontre)=NULL;
285     return NOEUD_CREE;
286     }
287    
288    
289    
290    
291    
292    
293    
294    
295    
296    
297     MG_TRIANGLE* MAILLEUR2D::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3)
298     {
299     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
300     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
301     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
302 francois 791 if (mgsegment1==NULL) mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
303     if (mgsegment2==NULL) mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
304     if (mgsegment3==NULL) mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,MAGIC::ORIGINE::MAILLEUR_AUTO);
305     M3D_TRIANGLE* mtriangle=new M3D_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::MAILLEUR_AUTO);
306 5 mg_maillage->ajouter_mg_triangle(mtriangle);
307     double qual=OPERATEUR::qualite_triangle(mgnoeud1->get_coord(),mgnoeud2->get_coord(),mgnoeud3->get_coord());
308     mtriangle->change_qualite(qual);
309     std::pair<const double,M3D_TRIANGLE*> tmp(mtriangle->get_qualite(),mtriangle);
310     return mtriangle;
311     }
312    
313     int MAILLEUR2D::triangle_est_dans_bon_sens(MG_FACE* face,MG_NOEUD* noeud1,MG_NOEUD *noeud2,MG_NOEUD *noeud3)
314     {
315     double* xyz1=noeud1->get_coord();
316     double* xyz2=noeud2->get_coord();
317     double* xyz3=noeud3->get_coord();
318     OT_VECTEUR_3D n1n3(xyz1,xyz3);
319     OT_VECTEUR_3D n1n2(xyz1,xyz2);
320     double uv[2];
321     face->inverser(uv,xyz1);
322     double n[3];
323     face->calcul_normale_unitaire(uv,n);
324     OT_VECTEUR_3D nmat=n&n1n3;
325     nmat.norme();
326     n1n2.norme();
327     if (nmat*n1n2<0) return 0;
328     return 1;
329     }
330    
331    
332    
333    
334    
335    
336    
337     // PARTIE DE CODE RECUPERER DE MON PROJET DE DEA
338     // TRES ILLISIBLE MAIS FONCTIONNEL
339     // POIL AU MARRON
340    
341     #define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
342     #define EGAL(x,y,eps) (float)fabs((double)(x-y))<eps
343     #define DETER(a,b,c,d) (a*d-b*c)
344    
345     int MAILLEUR2D::intersection_segment_segment(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
346     {
347    
348     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
349     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
350     double ua=decalage->decalage_parametre_u(noeud1->get_u(),du);
351     double va=decalage->decalage_parametre_v(noeud1->get_v(),dv);
352     double ub=decalage->decalage_parametre_u(noeud2->get_u(),du);
353     double vb=decalage->decalage_parametre_v(noeud2->get_v(),dv);
354     double um=decalage->decalage_parametre_u(noeud3->get_u(),du);
355     double vm=decalage->decalage_parametre_v(noeud3->get_v(),dv);
356 francois 536 double un=decalage->decalage_parametre_u(noeud4->get_u(),um,du);
357     double vn=decalage->decalage_parametre_v(noeud4->get_v(),vm,dv);
358    
359    
360 5 double ab[3];
361     double nm[3];
362     double am[3];
363     ab[0]=ub-ua;
364     ab[1]=vb-va;
365     ab[2]=0.;
366     nm[0]=um-un;
367     nm[1]=vm-vn;
368     nm[2]=0.;
369     am[0]=um-ua;
370     am[1]=vm-va;
371     am[2]=0.;
372     int equation[4];
373     equation[0]=1; /* etat de l'equation 0 */
374     equation[1]=1;
375     equation[2]=1;
376     equation[3]=3; /* cette variable comporte le bilan du nombre d'equation */
377     double eps2=PSCA(ab,ab);
378     double eps=sqrt(eps2);
379     eps=eps*0.0001;
380     eps2=eps2*0.0001;
381     /* recherche du nombre d'equation -> inter franche ou para ou confondu */
382     if ( (EGAL(ab[0],0,eps)) && (EGAL(nm[0],0,eps)) )
383     if (EGAL(am[0],0,eps)) equation[0]=0; else return(0);
384     if ( (EGAL(ab[1],0,eps)) && (EGAL(nm[1],0,eps)) )
385     if (EGAL(am[1],0,eps)) equation[1]=0; else return(0);
386     if ( (EGAL(ab[2],0,eps)) && (EGAL(nm[2],0,eps)) )
387     if (EGAL(am[2],0,eps)) equation[2]=0; else return(0);
388     equation[3]=equation[0]+equation[1]+equation[2];
389     if (equation[3]==3)
390     {
391     double det=DETER(ab[0],nm[0],ab[1],nm[1]);
392     if (fabs(det)>eps2)
393     {
394     det=1/det;
395     double sol1=det*DETER(am[0],nm[0],am[1],nm[1]);
396     double sol2=det*DETER(ab[0],am[0],ab[1],am[1]);
397     if ( (float)fabs((double)(sol1*ab[2]-sol2*nm[2]-am[2]))>eps2) return(0);
398     return(examine_solution(sol1,sol2,1));
399     }
400     else
401     {
402     equation[0]=0;
403     equation[3]=2;
404     /* on verifie la compatibilite des deux equations dont le det est nul*/
405     double tmp;
406     if (ab[0]!=0) tmp=ab[1]*am[0]/ab[0]; else tmp=nm[1]*am[0]/nm[0];
407     if (!(EGAL(tmp,am[1],eps))) return(0);
408     }
409     }
410     if (equation[3]==2)
411     {
412     /* on repere les equations qui existent */
413     int ne1;
414     int ne2;
415     if (equation[0]!=0)
416     {
417     ne1=0;
418     if (equation[1]!=0) ne2=1; else ne2=2;
419     }
420     else
421     {
422     ne1=1;
423     ne2=2;
424     }
425    
426     double det=DETER(ab[ne1],nm[ne1],ab[ne2],nm[ne2]);
427     if (fabs(det)>eps2)
428     {
429     det=1/det;
430     double sol1=det*DETER(am[ne1],nm[ne1],am[ne2],nm[ne2]);
431     double sol2=det*DETER(ab[ne1],am[ne1],ab[ne2],am[ne2]);
432     return(examine_solution(sol1,sol2,1));
433     }
434     else
435     {
436     equation[ne1]=0;
437     equation[3]=1;
438     /* on verifie la compatibilite des deux equations dont le det est nul */
439     double tmp;
440     if (ab[ne1]!=0) tmp=ab[ne2]*am[ne1]/ab[ne1]; else tmp=nm[ne2]*am[ne1]/nm[ne1];
441     if (!(EGAL(tmp,am[ne2],eps))) return(0);
442     }
443    
444     }
445     if (equation[3]==1)
446     {
447     /* on repere l' equation qui existe */
448     int ne1;
449     if (equation[0]!=0) ne1=0; else
450     if (equation[1]!=0) ne1=1; else ne1=2;
451     double an[3];
452     an[0]=un-ua;
453     an[1]=vn-va;
454     an[2]=0.;
455     double tmp=1/ab[ne1];
456     double sol1=am[ne1]*tmp;
457     double sol2=an[ne1]*tmp;
458     return(examine_solution(sol1,sol2,2));
459     }
460     return(0);
461     }
462    
463    
464    
465     int MAILLEUR2D::examine_solution(double sol1,double sol2,int type)
466     {
467     double epsilon=0.0001;
468    
469     if (type==1)
470     {
471     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
472     if ( ( (EGAL(sol1,0,epsilon)) || (EGAL(sol1,1,epsilon))) && ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) ) return 1;
473     if ( ( (EGAL(sol2,0,epsilon)) || (EGAL(sol2,1,epsilon))) && ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) ) return 1;
474     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>(-0.1-epsilon)) && ((sol2)<(1.1-epsilon)) ) return 1;
475     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) && (sol1>(-0.1-epsilon)) && ((sol1)<(1.1-epsilon)) ) return 1;
476    
477     }
478     if (type==2)
479     {
480     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) return 1;
481     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
482     if ( ((sol1)>(1+epsilon)) && ((-sol2)>epsilon) ) return 1;
483     if ( ((sol2)>(1+epsilon)) && ((-sol1)>epsilon) ) return 1;
484     }
485     return 0;
486     }
487     #undef EGAL
488     #undef PSCA
489     #undef DETER
490     ///FIN DU CODE DE DEA