ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur2d_outil.cpp
Revision: 791
Committed: Fri Mar 18 21:23:00 2016 UTC (9 years, 1 month ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur2d_outil.cpp
File size: 18776 byte(s)
Log Message:
Transferts de la definition des origines dans le namespace MAGiC

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    
132 francois 791 MG_SEGMENT* mgsegment=mg_maillage->ajouter_mg_segment(mgface,noeud1,noeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
133 francois 632 ntree_de_segment->inserer(mgsegment);
134 5 *nv_segment=mgsegment;
135     return(1);
136     }
137    
138    
139     void MAILLEUR2D::supprime_segment(MG_SEGMENT* mgsegment)
140     {
141 francois 632 ntree_de_segment->supprimer(mgsegment);
142 5 mg_maillage->supprimer_mg_segmentid(mgsegment->get_id());
143     }
144    
145    
146     int MAILLEUR2D::genere_noeud(MG_FACE* mgface,MG_FRONT_2D* front,MG_FRONT_2D **front_rencontre,MG_NOEUD **noeud_rencontre)
147     {
148     OT_VECTEUR_3D w(0.,0.,1.);
149     MG_NOEUD* noeud1=front->get_noeud1();
150     MG_NOEUD* noeud2=front->get_noeud2();
151     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
152     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
153     double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
154     double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
155     double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
156     double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
157    
158 francois 532 double longueur_segment=metrique->calcule_longueur_segment_metrique(mgface,front->get_segment(),decalage,pas);
159 5 double longueur_desiree=0.8660254037844386*(MAILLEUR2D::priorite_metrique+longueur_segment-longueur_segment*MAILLEUR2D::priorite_metrique);
160     longueur_desiree=longueur_desiree/front->get_ifail();
161    
162    
163     double umilieu=0.5*(u1+u2);
164     double vmilieu=0.5*(v1+v2);
165     double param_milieu[2]={umilieu-du,vmilieu-dv};
166     double E,F,G;
167    
168     mgface->get_EFG(param_milieu,E,F,G);
169     double deno=E*(u2-u1)*(u2-u1)+G*(v2-v1)*(v2-v1)+2*F*(u2-u1)*(v2-v1);
170     double nume=(v2-v1)*(u2-u1)*(E-G)+F*((v2-v1)*(v2-v1)-(u2-u1)*(u2-u1));
171     double teta=atan(nume/deno);
172    
173    
174 francois 532 double udecale;
175     double vdecale;
176 francois 632 int res=metrique->ajuste_distance_metrique(mgface,u1,v1,u2,v2,udecale,vdecale,longueur_desiree,du,dv,teta,pas);
177     if (res==0) return ERREUR;
178 5 double u=decalage->decalage_parametre_u(udecale,-du);
179     double v=decalage->decalage_parametre_v(vdecale,-dv);
180     double param_noeud_cree[2]={u,v};
181    
182    
183     // recherhce des element du front proche de u,v;
184     TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
185 francois 632 double xyz[3],uv[2]={u,v};
186     mgface->evaluer(uv,xyz);
187     double xmin=std::min(noeud1->get_x(),noeud2->get_x());xmin=std::min(xmin,xyz[0]);
188     double xmax=std::max(noeud1->get_x(),noeud2->get_x());xmax=std::max(xmax,xyz[0]);
189     double ymin=std::min(noeud1->get_y(),noeud2->get_y());ymin=std::min(ymin,xyz[1]);
190     double ymax=std::max(noeud1->get_y(),noeud2->get_y());ymax=std::max(ymax,xyz[1]);
191     double zmin=std::min(noeud1->get_z(),noeud2->get_z());zmin=std::min(zmin,xyz[2]);
192     double zmax=std::max(noeud1->get_z(),noeud2->get_z());zmax=std::max(zmax,xyz[2]);
193     BOITE_3D boite_recherche(xmin,ymin,zmin,xmax,ymax,zmax);
194 5 double longueur_recherche=boite_recherche.get_rayon()*2.;
195    
196 francois 632
197    
198     ntree_de_front->rechercher(xyz[0],xyz[1],xyz[2],longueur_recherche,liste_trouvee);
199 5 int nb_entite=liste_trouvee.get_nb();
200     double distance_reference = -1.;
201     double angle_reference=2.*M_PI;
202     unsigned long id_noeud_reference = 0;
203     MG_FRONT_2D* front_reference=NULL;
204    
205     for (int i=0;i<nb_entite;i++)
206     {
207     MG_FRONT_2D* front_courant=liste_trouvee.get(i);
208     MG_NOEUD* noeud_front1=front_courant->get_noeud1();
209     MG_NOEUD* noeud_front2=front_courant->get_noeud2();
210     double unoeudfront1=decalage->decalage_parametre_u(noeud_front1->get_u(),du);
211     double vnoeudfront1=decalage->decalage_parametre_v(noeud_front1->get_v(),dv);
212     double unoeudfront2=decalage->decalage_parametre_u(noeud_front2->get_u(),du);
213     double vnoeudfront2=decalage->decalage_parametre_v(noeud_front2->get_v(),dv);
214 francois 532 double distance_noeudfront1=metrique->calcule_distance_metrique(mgface,udecale,vdecale,unoeudfront1,vnoeudfront1,du,dv,pas);
215 5 if ((distance_noeudfront1<distance_reference) || (distance_reference<0.) || (id_noeud_reference==noeud_front1->get_id()))
216     if (noeud_front1->get_id()!=noeud1->get_id())
217     if (noeud_front1->get_id()!=noeud2->get_id())
218     if (id_noeud_reference==noeud_front1->get_id())
219     {
220     refresh();
221     OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
222     OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
223     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
224     vecteur_baseu.norme();
225     vecteur_basev.norme();
226     vecteur_front.norme();
227     double cosangle=vecteur_baseu*vecteur_front;
228     double sinangle=vecteur_basev*vecteur_front;
229     sinangle=-sinangle;
230     if (cosangle>1.) cosangle=1.;
231     if (cosangle<-1.) cosangle=(-1.);
232     double angle=acos(cosangle);
233     if (sinangle<(-0.0001)) angle=(-angle);
234     if (angle<0.) angle=angle+2.*M_PI;
235     if (angle<angle_reference)
236     {
237     angle_reference=angle;
238     front_reference=front_courant;
239     }
240     }
241     else
242     {
243     refresh();
244     OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
245     OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
246     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
247     vecteur_baseu.norme();
248     vecteur_basev.norme();
249     vecteur_front.norme();
250     double cosangle=vecteur_baseu*vecteur_front;
251     double sinangle=vecteur_basev*vecteur_front;
252     sinangle=-sinangle;
253     if (cosangle>1.) cosangle=1.;
254     if (cosangle<-1.) cosangle=(-1.);
255     double angle=acos(cosangle);
256     if (sinangle<(-0.0001)) angle=(-angle);
257     if (angle<0.) angle=angle+2.*M_PI;
258     distance_reference=distance_noeudfront1;
259     angle_reference=angle;
260     id_noeud_reference=noeud_front1->get_id();
261     front_reference=front_courant;
262     }
263     }
264    
265     if ((front_reference!=NULL) && (distance_reference<longueur_desiree*2./sqrt(3)))
266     {
267     MG_NOEUD* noeud=mg_maillage->get_mg_noeudid(id_noeud_reference);
268     double uref=decalage->decalage_parametre_u(noeud->get_u(),du);
269     double vref=decalage->decalage_parametre_v(noeud->get_v(),dv);
270 francois 532 double distance=metrique->calcule_distance_metrique(mgface,umilieu,vmilieu,uref,vref,du,dv,pas);
271 5 if (distance<1.5*longueur_desiree)
272     {
273     (*front_rencontre)=front_reference;
274     (*noeud_rencontre)=noeud;
275     return FRONT_RENCONTRE;
276     }
277     }
278 francois 791 MG_NOEUD* noeud=mg_maillage->ajouter_mg_noeud(mgface,xyz[0],xyz[1],xyz[2],MAGIC::ORIGINE::MAILLEUR_AUTO);
279 5 noeud->change_u(param_noeud_cree[0]);
280     noeud->change_v(param_noeud_cree[1]);
281     (*noeud_rencontre)=noeud;
282     (*front_rencontre)=NULL;
283     return NOEUD_CREE;
284     }
285    
286    
287    
288    
289    
290    
291    
292    
293    
294    
295     MG_TRIANGLE* MAILLEUR2D::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3)
296     {
297     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
298     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
299     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
300 francois 791 if (mgsegment1==NULL) mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
301     if (mgsegment2==NULL) mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
302     if (mgsegment3==NULL) mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,MAGIC::ORIGINE::MAILLEUR_AUTO);
303     M3D_TRIANGLE* mtriangle=new M3D_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::MAILLEUR_AUTO);
304 5 mg_maillage->ajouter_mg_triangle(mtriangle);
305     double qual=OPERATEUR::qualite_triangle(mgnoeud1->get_coord(),mgnoeud2->get_coord(),mgnoeud3->get_coord());
306     mtriangle->change_qualite(qual);
307     std::pair<const double,M3D_TRIANGLE*> tmp(mtriangle->get_qualite(),mtriangle);
308     return mtriangle;
309     }
310    
311     int MAILLEUR2D::triangle_est_dans_bon_sens(MG_FACE* face,MG_NOEUD* noeud1,MG_NOEUD *noeud2,MG_NOEUD *noeud3)
312     {
313     double* xyz1=noeud1->get_coord();
314     double* xyz2=noeud2->get_coord();
315     double* xyz3=noeud3->get_coord();
316     OT_VECTEUR_3D n1n3(xyz1,xyz3);
317     OT_VECTEUR_3D n1n2(xyz1,xyz2);
318     double uv[2];
319     face->inverser(uv,xyz1);
320     double n[3];
321     face->calcul_normale_unitaire(uv,n);
322     OT_VECTEUR_3D nmat=n&n1n3;
323     nmat.norme();
324     n1n2.norme();
325     if (nmat*n1n2<0) return 0;
326     return 1;
327     }
328    
329    
330    
331    
332    
333    
334    
335     // PARTIE DE CODE RECUPERER DE MON PROJET DE DEA
336     // TRES ILLISIBLE MAIS FONCTIONNEL
337     // POIL AU MARRON
338    
339     #define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
340     #define EGAL(x,y,eps) (float)fabs((double)(x-y))<eps
341     #define DETER(a,b,c,d) (a*d-b*c)
342    
343     int MAILLEUR2D::intersection_segment_segment(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
344     {
345    
346     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
347     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
348     double ua=decalage->decalage_parametre_u(noeud1->get_u(),du);
349     double va=decalage->decalage_parametre_v(noeud1->get_v(),dv);
350     double ub=decalage->decalage_parametre_u(noeud2->get_u(),du);
351     double vb=decalage->decalage_parametre_v(noeud2->get_v(),dv);
352     double um=decalage->decalage_parametre_u(noeud3->get_u(),du);
353     double vm=decalage->decalage_parametre_v(noeud3->get_v(),dv);
354 francois 536 double un=decalage->decalage_parametre_u(noeud4->get_u(),um,du);
355     double vn=decalage->decalage_parametre_v(noeud4->get_v(),vm,dv);
356    
357    
358 5 double ab[3];
359     double nm[3];
360     double am[3];
361     ab[0]=ub-ua;
362     ab[1]=vb-va;
363     ab[2]=0.;
364     nm[0]=um-un;
365     nm[1]=vm-vn;
366     nm[2]=0.;
367     am[0]=um-ua;
368     am[1]=vm-va;
369     am[2]=0.;
370     int equation[4];
371     equation[0]=1; /* etat de l'equation 0 */
372     equation[1]=1;
373     equation[2]=1;
374     equation[3]=3; /* cette variable comporte le bilan du nombre d'equation */
375     double eps2=PSCA(ab,ab);
376     double eps=sqrt(eps2);
377     eps=eps*0.0001;
378     eps2=eps2*0.0001;
379     /* recherche du nombre d'equation -> inter franche ou para ou confondu */
380     if ( (EGAL(ab[0],0,eps)) && (EGAL(nm[0],0,eps)) )
381     if (EGAL(am[0],0,eps)) equation[0]=0; else return(0);
382     if ( (EGAL(ab[1],0,eps)) && (EGAL(nm[1],0,eps)) )
383     if (EGAL(am[1],0,eps)) equation[1]=0; else return(0);
384     if ( (EGAL(ab[2],0,eps)) && (EGAL(nm[2],0,eps)) )
385     if (EGAL(am[2],0,eps)) equation[2]=0; else return(0);
386     equation[3]=equation[0]+equation[1]+equation[2];
387     if (equation[3]==3)
388     {
389     double det=DETER(ab[0],nm[0],ab[1],nm[1]);
390     if (fabs(det)>eps2)
391     {
392     det=1/det;
393     double sol1=det*DETER(am[0],nm[0],am[1],nm[1]);
394     double sol2=det*DETER(ab[0],am[0],ab[1],am[1]);
395     if ( (float)fabs((double)(sol1*ab[2]-sol2*nm[2]-am[2]))>eps2) return(0);
396     return(examine_solution(sol1,sol2,1));
397     }
398     else
399     {
400     equation[0]=0;
401     equation[3]=2;
402     /* on verifie la compatibilite des deux equations dont le det est nul*/
403     double tmp;
404     if (ab[0]!=0) tmp=ab[1]*am[0]/ab[0]; else tmp=nm[1]*am[0]/nm[0];
405     if (!(EGAL(tmp,am[1],eps))) return(0);
406     }
407     }
408     if (equation[3]==2)
409     {
410     /* on repere les equations qui existent */
411     int ne1;
412     int ne2;
413     if (equation[0]!=0)
414     {
415     ne1=0;
416     if (equation[1]!=0) ne2=1; else ne2=2;
417     }
418     else
419     {
420     ne1=1;
421     ne2=2;
422     }
423    
424     double det=DETER(ab[ne1],nm[ne1],ab[ne2],nm[ne2]);
425     if (fabs(det)>eps2)
426     {
427     det=1/det;
428     double sol1=det*DETER(am[ne1],nm[ne1],am[ne2],nm[ne2]);
429     double sol2=det*DETER(ab[ne1],am[ne1],ab[ne2],am[ne2]);
430     return(examine_solution(sol1,sol2,1));
431     }
432     else
433     {
434     equation[ne1]=0;
435     equation[3]=1;
436     /* on verifie la compatibilite des deux equations dont le det est nul */
437     double tmp;
438     if (ab[ne1]!=0) tmp=ab[ne2]*am[ne1]/ab[ne1]; else tmp=nm[ne2]*am[ne1]/nm[ne1];
439     if (!(EGAL(tmp,am[ne2],eps))) return(0);
440     }
441    
442     }
443     if (equation[3]==1)
444     {
445     /* on repere l' equation qui existe */
446     int ne1;
447     if (equation[0]!=0) ne1=0; else
448     if (equation[1]!=0) ne1=1; else ne1=2;
449     double an[3];
450     an[0]=un-ua;
451     an[1]=vn-va;
452     an[2]=0.;
453     double tmp=1/ab[ne1];
454     double sol1=am[ne1]*tmp;
455     double sol2=an[ne1]*tmp;
456     return(examine_solution(sol1,sol2,2));
457     }
458     return(0);
459     }
460    
461    
462    
463     int MAILLEUR2D::examine_solution(double sol1,double sol2,int type)
464     {
465     double epsilon=0.0001;
466    
467     if (type==1)
468     {
469     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
470     if ( ( (EGAL(sol1,0,epsilon)) || (EGAL(sol1,1,epsilon))) && ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) ) return 1;
471     if ( ( (EGAL(sol2,0,epsilon)) || (EGAL(sol2,1,epsilon))) && ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) ) return 1;
472     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>(-0.1-epsilon)) && ((sol2)<(1.1-epsilon)) ) return 1;
473     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) && (sol1>(-0.1-epsilon)) && ((sol1)<(1.1-epsilon)) ) return 1;
474    
475     }
476     if (type==2)
477     {
478     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) return 1;
479     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
480     if ( ((sol1)>(1+epsilon)) && ((-sol2)>epsilon) ) return 1;
481     if ( ((sol2)>(1+epsilon)) && ((-sol1)>epsilon) ) return 1;
482     }
483     return 0;
484     }
485     #undef EGAL
486     #undef PSCA
487     #undef DETER
488     ///FIN DU CODE DE DEA