ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur2d_outil.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (10 months, 4 weeks ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur2d_outil.cpp
File size: 19235 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

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