ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur2d_outil.cpp
Revision: 1150
Committed: Tue May 28 12:57:28 2024 UTC (11 months, 2 weeks ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur2d_outil.cpp
File size: 19376 byte(s)
Log Message:
deplacement des contantes 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 francois 1150 if (type_validation==MAGIC::MAILLEURFRONTALETAT::TOUS_SEGMENT)
103 5 {
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 francois 1150 if (type_validation==MAGIC::MAILLEURFRONTALETAT::TOUS_FRONT)
117 5 {
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 1150 if (res==0) return MAGIC::MAILLEURFRONTALETAT::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 francois 1095 double du2=decalage->calcul_decalage_parametre_u(noeud_front1->get_u());
213     double dv2=decalage->calcul_decalage_parametre_v(noeud_front1->get_v());
214     double u1d2=decalage->decalage_parametre_u(noeud1->get_u(),du2);
215     double v1d2=decalage->decalage_parametre_v(noeud1->get_v(),dv2);
216     double udecale2=decalage->decalage_parametre_u(u,du2);
217     double vdecale2=decalage->decalage_parametre_v(v,dv2);
218     double unoeudfront1=decalage->decalage_parametre_u(noeud_front1->get_u(),du2);
219     double vnoeudfront1=decalage->decalage_parametre_v(noeud_front1->get_v(),dv2);
220     double unoeudfront2=decalage->decalage_parametre_u(noeud_front2->get_u(),du2);
221     double vnoeudfront2=decalage->decalage_parametre_v(noeud_front2->get_v(),dv2);
222     double distance_noeudfront1=metrique->calcule_distance_metrique(mgface,udecale2,vdecale2,unoeudfront1,vnoeudfront1,du2,dv2,pas);
223 5 if ((distance_noeudfront1<distance_reference) || (distance_reference<0.) || (id_noeud_reference==noeud_front1->get_id()))
224     if (noeud_front1->get_id()!=noeud1->get_id())
225     if (noeud_front1->get_id()!=noeud2->get_id())
226     if (id_noeud_reference==noeud_front1->get_id())
227     {
228     refresh();
229 francois 1095 OT_VECTEUR_3D vecteur_baseu(u1d2-unoeudfront1,v1d2-vnoeudfront1,0.);
230 5 OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
231     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
232     vecteur_baseu.norme();
233     vecteur_basev.norme();
234     vecteur_front.norme();
235     double cosangle=vecteur_baseu*vecteur_front;
236     double sinangle=vecteur_basev*vecteur_front;
237     sinangle=-sinangle;
238     if (cosangle>1.) cosangle=1.;
239     if (cosangle<-1.) cosangle=(-1.);
240     double angle=acos(cosangle);
241     if (sinangle<(-0.0001)) angle=(-angle);
242     if (angle<0.) angle=angle+2.*M_PI;
243     if (angle<angle_reference)
244     {
245     angle_reference=angle;
246     front_reference=front_courant;
247     }
248     }
249     else
250     {
251     refresh();
252 francois 1095 OT_VECTEUR_3D vecteur_baseu(u1d2-unoeudfront1,v1d2-vnoeudfront1,0.);
253 5 OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
254     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
255     vecteur_baseu.norme();
256     vecteur_basev.norme();
257     vecteur_front.norme();
258     double cosangle=vecteur_baseu*vecteur_front;
259     double sinangle=vecteur_basev*vecteur_front;
260     sinangle=-sinangle;
261     if (cosangle>1.) cosangle=1.;
262     if (cosangle<-1.) cosangle=(-1.);
263     double angle=acos(cosangle);
264     if (sinangle<(-0.0001)) angle=(-angle);
265     if (angle<0.) angle=angle+2.*M_PI;
266     distance_reference=distance_noeudfront1;
267     angle_reference=angle;
268     id_noeud_reference=noeud_front1->get_id();
269     front_reference=front_courant;
270     }
271     }
272    
273     if ((front_reference!=NULL) && (distance_reference<longueur_desiree*2./sqrt(3)))
274     {
275     MG_NOEUD* noeud=mg_maillage->get_mg_noeudid(id_noeud_reference);
276     double uref=decalage->decalage_parametre_u(noeud->get_u(),du);
277     double vref=decalage->decalage_parametre_v(noeud->get_v(),dv);
278 francois 532 double distance=metrique->calcule_distance_metrique(mgface,umilieu,vmilieu,uref,vref,du,dv,pas);
279 5 if (distance<1.5*longueur_desiree)
280     {
281     (*front_rencontre)=front_reference;
282     (*noeud_rencontre)=noeud;
283 francois 1150 return MAGIC::MAILLEURFRONTALETAT::FRONT_RENCONTRE;
284 5 }
285     }
286 francois 791 MG_NOEUD* noeud=mg_maillage->ajouter_mg_noeud(mgface,xyz[0],xyz[1],xyz[2],MAGIC::ORIGINE::MAILLEUR_AUTO);
287 5 noeud->change_u(param_noeud_cree[0]);
288     noeud->change_v(param_noeud_cree[1]);
289     (*noeud_rencontre)=noeud;
290     (*front_rencontre)=NULL;
291 francois 1150 return MAGIC::MAILLEURFRONTALETAT::NOEUD_CREE;
292 5 }
293    
294    
295    
296    
297    
298    
299    
300    
301    
302    
303     MG_TRIANGLE* MAILLEUR2D::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3)
304     {
305     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
306     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
307     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
308 francois 791 if (mgsegment1==NULL) mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
309     if (mgsegment2==NULL) mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
310     if (mgsegment3==NULL) mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,MAGIC::ORIGINE::MAILLEUR_AUTO);
311     M3D_TRIANGLE* mtriangle=new M3D_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::MAILLEUR_AUTO);
312 5 mg_maillage->ajouter_mg_triangle(mtriangle);
313     double qual=OPERATEUR::qualite_triangle(mgnoeud1->get_coord(),mgnoeud2->get_coord(),mgnoeud3->get_coord());
314     mtriangle->change_qualite(qual);
315     std::pair<const double,M3D_TRIANGLE*> tmp(mtriangle->get_qualite(),mtriangle);
316     return mtriangle;
317     }
318    
319     int MAILLEUR2D::triangle_est_dans_bon_sens(MG_FACE* face,MG_NOEUD* noeud1,MG_NOEUD *noeud2,MG_NOEUD *noeud3)
320     {
321     double* xyz1=noeud1->get_coord();
322     double* xyz2=noeud2->get_coord();
323     double* xyz3=noeud3->get_coord();
324     OT_VECTEUR_3D n1n3(xyz1,xyz3);
325     OT_VECTEUR_3D n1n2(xyz1,xyz2);
326     double uv[2];
327     face->inverser(uv,xyz1);
328     double n[3];
329     face->calcul_normale_unitaire(uv,n);
330     OT_VECTEUR_3D nmat=n&n1n3;
331     nmat.norme();
332     n1n2.norme();
333     if (nmat*n1n2<0) return 0;
334     return 1;
335     }
336    
337    
338    
339    
340    
341    
342    
343     // PARTIE DE CODE RECUPERER DE MON PROJET DE DEA
344     // TRES ILLISIBLE MAIS FONCTIONNEL
345     // POIL AU MARRON
346    
347     #define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
348     #define EGAL(x,y,eps) (float)fabs((double)(x-y))<eps
349     #define DETER(a,b,c,d) (a*d-b*c)
350    
351     int MAILLEUR2D::intersection_segment_segment(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
352     {
353    
354     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
355     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
356     double ua=decalage->decalage_parametre_u(noeud1->get_u(),du);
357     double va=decalage->decalage_parametre_v(noeud1->get_v(),dv);
358     double ub=decalage->decalage_parametre_u(noeud2->get_u(),du);
359     double vb=decalage->decalage_parametre_v(noeud2->get_v(),dv);
360     double um=decalage->decalage_parametre_u(noeud3->get_u(),du);
361     double vm=decalage->decalage_parametre_v(noeud3->get_v(),dv);
362 francois 536 double un=decalage->decalage_parametre_u(noeud4->get_u(),um,du);
363     double vn=decalage->decalage_parametre_v(noeud4->get_v(),vm,dv);
364    
365    
366 5 double ab[3];
367     double nm[3];
368     double am[3];
369     ab[0]=ub-ua;
370     ab[1]=vb-va;
371     ab[2]=0.;
372     nm[0]=um-un;
373     nm[1]=vm-vn;
374     nm[2]=0.;
375     am[0]=um-ua;
376     am[1]=vm-va;
377     am[2]=0.;
378     int equation[4];
379     equation[0]=1; /* etat de l'equation 0 */
380     equation[1]=1;
381     equation[2]=1;
382     equation[3]=3; /* cette variable comporte le bilan du nombre d'equation */
383     double eps2=PSCA(ab,ab);
384     double eps=sqrt(eps2);
385     eps=eps*0.0001;
386     eps2=eps2*0.0001;
387     /* recherche du nombre d'equation -> inter franche ou para ou confondu */
388     if ( (EGAL(ab[0],0,eps)) && (EGAL(nm[0],0,eps)) )
389     if (EGAL(am[0],0,eps)) equation[0]=0; else return(0);
390     if ( (EGAL(ab[1],0,eps)) && (EGAL(nm[1],0,eps)) )
391     if (EGAL(am[1],0,eps)) equation[1]=0; else return(0);
392     if ( (EGAL(ab[2],0,eps)) && (EGAL(nm[2],0,eps)) )
393     if (EGAL(am[2],0,eps)) equation[2]=0; else return(0);
394     equation[3]=equation[0]+equation[1]+equation[2];
395     if (equation[3]==3)
396     {
397     double det=DETER(ab[0],nm[0],ab[1],nm[1]);
398     if (fabs(det)>eps2)
399     {
400     det=1/det;
401     double sol1=det*DETER(am[0],nm[0],am[1],nm[1]);
402     double sol2=det*DETER(ab[0],am[0],ab[1],am[1]);
403     if ( (float)fabs((double)(sol1*ab[2]-sol2*nm[2]-am[2]))>eps2) return(0);
404     return(examine_solution(sol1,sol2,1));
405     }
406     else
407     {
408     equation[0]=0;
409     equation[3]=2;
410     /* on verifie la compatibilite des deux equations dont le det est nul*/
411     double tmp;
412     if (ab[0]!=0) tmp=ab[1]*am[0]/ab[0]; else tmp=nm[1]*am[0]/nm[0];
413     if (!(EGAL(tmp,am[1],eps))) return(0);
414     }
415     }
416     if (equation[3]==2)
417     {
418     /* on repere les equations qui existent */
419     int ne1;
420     int ne2;
421     if (equation[0]!=0)
422     {
423     ne1=0;
424     if (equation[1]!=0) ne2=1; else ne2=2;
425     }
426     else
427     {
428     ne1=1;
429     ne2=2;
430     }
431    
432     double det=DETER(ab[ne1],nm[ne1],ab[ne2],nm[ne2]);
433     if (fabs(det)>eps2)
434     {
435     det=1/det;
436     double sol1=det*DETER(am[ne1],nm[ne1],am[ne2],nm[ne2]);
437     double sol2=det*DETER(ab[ne1],am[ne1],ab[ne2],am[ne2]);
438     return(examine_solution(sol1,sol2,1));
439     }
440     else
441     {
442     equation[ne1]=0;
443     equation[3]=1;
444     /* on verifie la compatibilite des deux equations dont le det est nul */
445     double tmp;
446     if (ab[ne1]!=0) tmp=ab[ne2]*am[ne1]/ab[ne1]; else tmp=nm[ne2]*am[ne1]/nm[ne1];
447     if (!(EGAL(tmp,am[ne2],eps))) return(0);
448     }
449    
450     }
451     if (equation[3]==1)
452     {
453     /* on repere l' equation qui existe */
454     int ne1;
455     if (equation[0]!=0) ne1=0; else
456     if (equation[1]!=0) ne1=1; else ne1=2;
457     double an[3];
458     an[0]=un-ua;
459     an[1]=vn-va;
460     an[2]=0.;
461     double tmp=1/ab[ne1];
462     double sol1=am[ne1]*tmp;
463     double sol2=an[ne1]*tmp;
464     return(examine_solution(sol1,sol2,2));
465     }
466     return(0);
467     }
468    
469    
470    
471     int MAILLEUR2D::examine_solution(double sol1,double sol2,int type)
472     {
473     double epsilon=0.0001;
474    
475     if (type==1)
476     {
477     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
478     if ( ( (EGAL(sol1,0,epsilon)) || (EGAL(sol1,1,epsilon))) && ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) ) return 1;
479     if ( ( (EGAL(sol2,0,epsilon)) || (EGAL(sol2,1,epsilon))) && ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) ) return 1;
480     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>(-0.1-epsilon)) && ((sol2)<(1.1-epsilon)) ) return 1;
481     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) && (sol1>(-0.1-epsilon)) && ((sol1)<(1.1-epsilon)) ) return 1;
482    
483     }
484     if (type==2)
485     {
486     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) return 1;
487     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
488     if ( ((sol1)>(1+epsilon)) && ((-sol2)>epsilon) ) return 1;
489     if ( ((sol2)>(1+epsilon)) && ((-sol1)>epsilon) ) return 1;
490     }
491     return 0;
492     }
493     #undef EGAL
494     #undef PSCA
495     #undef DETER
496     ///FIN DU CODE DE DEA