ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur2d_outil.cpp
Revision: 632
Committed: Thu Jan 15 18:40:00 2015 UTC (10 years, 3 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur2d_outil.cpp
File size: 18680 byte(s)
Log Message:
Changement de l'espace de voisinage dans le mailleur 2D et 3D. On utilise un ntree (octree "anisotrope" qui ne se divise pas necessairement en 8). En 2D l'espace de voisinage est maintenant sur l'espace reelle au lieu d'être sur l'espace parametrique.
+Mise sous forme de parametres de certains choix stratégiques du mailleur

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 35 MG_SEGMENT* mgsegment=mg_maillage->ajouter_mg_segment(mgface,noeud1,noeud2,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 632 MG_NOEUD* noeud=mg_maillage->ajouter_mg_noeud(mgface,xyz[0],xyz[1],xyz[2],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 35 if (mgsegment1==NULL) mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,MAILLEUR_AUTO);
301     if (mgsegment2==NULL) mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,MAILLEUR_AUTO);
302     if (mgsegment3==NULL) mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,MAILLEUR_AUTO);
303 francois 54 M3D_TRIANGLE* mtriangle=new M3D_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,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