ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur2d_outil.cpp
Revision: 469
Committed: Wed Nov 27 20:02:22 2013 UTC (11 years, 5 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur2d_outil.cpp
File size: 26532 byte(s)
Log Message:
ajout d'un generateur de carte de taille constante pour les nouvelles cartes et decoupage de la librairie mailleur en 2 pour cause d'appel cyclique

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     #include "ot_boite_2D.h"
30     #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    
93     if (type_validation==TOUS_SEGMENT)
94     {
95     TPL_MAP_ENTITE<MG_SEGMENT*> liste_trouvee;
96     quadtree_de_segment->rechercher(uii,vii,longueur/2.,liste_trouvee);
97     for (int j=0;j<liste_trouvee.get_nb();j++)
98     {
99     MG_SEGMENT* mgsegment=liste_trouvee.get(j);
100     if ( ((noeud1==mgsegment->get_noeud1()) && (noeud2==mgsegment->get_noeud2())) || ((noeud1==mgsegment->get_noeud2()) && (noeud2==mgsegment->get_noeud1())))
101     return(0);
102     int res=intersection_segment_segment(noeud1,noeud2,mgsegment->get_noeud1(),mgsegment->get_noeud2());
103     if (res!=0)
104     return(0);
105     }
106     }
107     if (type_validation==TOUS_FRONT)
108     {
109     TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
110     quadtree_de_front->rechercher(uii,vii,longueur/2.,liste_trouvee);
111     for (int j=0;j<liste_trouvee.get_nb();j++)
112     {
113     MG_FRONT_2D* ft=liste_trouvee.get(j);
114     MG_SEGMENT* mgsegment=ft->get_segment();
115     if ( ((noeud1==mgsegment->get_noeud1()) && (noeud2==mgsegment->get_noeud2())) || ((noeud1==mgsegment->get_noeud2()) && (noeud2==mgsegment->get_noeud1())))
116     return(0);
117     int res=intersection_segment_segment(noeud1,noeud2,mgsegment->get_noeud1(),mgsegment->get_noeud2());
118     if (res!=0)
119     return(0);
120     }
121     }
122    
123 francois 35 MG_SEGMENT* mgsegment=mg_maillage->ajouter_mg_segment(mgface,noeud1,noeud2,MAILLEUR_AUTO);
124 5 quadtree_de_segment->inserer(mgsegment);
125     *nv_segment=mgsegment;
126     return(1);
127     }
128    
129    
130     void MAILLEUR2D::supprime_segment(MG_SEGMENT* mgsegment)
131     {
132     quadtree_de_segment->supprimer(mgsegment);
133     mg_maillage->supprimer_mg_segmentid(mgsegment->get_id());
134     }
135    
136    
137     int MAILLEUR2D::genere_noeud(MG_FACE* mgface,MG_FRONT_2D* front,MG_FRONT_2D **front_rencontre,MG_NOEUD **noeud_rencontre)
138     {
139     OT_VECTEUR_3D w(0.,0.,1.);
140     MG_NOEUD* noeud1=front->get_noeud1();
141     MG_NOEUD* noeud2=front->get_noeud2();
142     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
143     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
144     double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
145     double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
146     double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
147     double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
148    
149     double longueur_segment=calcule_longueur_segment_metrique(mgface,front->get_segment());
150     double longueur_desiree=0.8660254037844386*(MAILLEUR2D::priorite_metrique+longueur_segment-longueur_segment*MAILLEUR2D::priorite_metrique);
151     longueur_desiree=longueur_desiree/front->get_ifail();
152    
153    
154     double longueur_calculee=0.;
155     double longueur_calculee_pas_precedent;
156     double umilieu=0.5*(u1+u2);
157     double vmilieu=0.5*(v1+v2);
158     double param_milieu[2]={umilieu-du,vmilieu-dv};
159     double E,F,G;
160    
161     // je pense qu il faudrait changer ici
162     // perpendiculaire dans l espace parametrique ou perpendiculaire dans la metrique ou dans la combinaison des deux
163     mgface->get_EFG(param_milieu,E,F,G);
164     double deno=E*(u2-u1)*(u2-u1)+G*(v2-v1)*(v2-v1)+2*F*(u2-u1)*(v2-v1);
165     double nume=(v2-v1)*(u2-u1)*(E-G)+F*((v2-v1)*(v2-v1)-(u2-u1)*(u2-u1));
166     double teta=atan(nume/deno);
167    
168     /*double tenseur_metrique[9];
169     double xyzdu[3];
170     double xyzdv[3];
171     metrique->evaluer(param_milieu,tenseur_metrique);
172     mgface->deriver(param_milieu,xyzdu,xyzdv);
173    
174     E = (xyzdu[0]*tenseur_metrique[0]+xyzdu[1]*tenseur_metrique[1]+xyzdu[2]*tenseur_metrique[2])*xyzdu[0]+(xyzdu[0]*tenseur_metrique[1]+xyzdu[1]*tenseur_metrique[4]+xyzdu[2]*tenseur_metrique[5])*xyzdu[1]+(xyzdu[0]*tenseur_metrique[2]+xyzdu[1]*tenseur_metrique[5]+xyzdu[2]*tenseur_metrique[8])*xyzdu[2];
175     F = (xyzdv[0]*tenseur_metrique[0]+xyzdv[1]*tenseur_metrique[1]+xyzdv[2]*tenseur_metrique[2])*xyzdu[0]+(xyzdv[0]*tenseur_metrique[1]+xyzdv[1]*tenseur_metrique[4]+xyzdv[2]*tenseur_metrique[5])*xyzdu[1]+(xyzdv[0]*tenseur_metrique[2]+xyzdv[1]*tenseur_metrique[5]+xyzdv[2]*tenseur_metrique[8])*xyzdu[2];
176     G = (xyzdv[0]*tenseur_metrique[0]+xyzdv[1]*tenseur_metrique[1]+xyzdv[2]*tenseur_metrique[2])*xyzdv[0]+(xyzdv[0]*tenseur_metrique[1]+xyzdv[1]*tenseur_metrique[4]+xyzdv[2]*tenseur_metrique[5])*xyzdv[1]+(xyzdv[0]*tenseur_metrique[2]+xyzdv[1]*tenseur_metrique[5]+xyzdv[2]*tenseur_metrique[8])*xyzdv[2];
177    
178     double deno=E*(u2-u1)*(u2-u1)+G*(v2-v1)*(v2-v1)+2*F*(u2-u1)*(v2-v1);
179     double nume=(v2-v1)*(u2-u1)*(E-G)+F*((v2-v1)*(v2-v1)-(u2-u1)*(u2-u1));
180     double teta=atan(nume/deno);
181    
182     */
183    
184     OT_VECTEUR_3D vecteur_du_front(u2-u1,v2-v1,0.);
185     OT_VECTEUR_3D vecteur_normal_front=w&vecteur_du_front;
186     double longueur_parametrique=vecteur_du_front.get_longueur();
187     vecteur_du_front.norme();
188     vecteur_normal_front.norme();
189     double dt=0.03125;
190     double ti;
191     double tii=0;
192     double uii;
193     double vii;
194     while (longueur_calculee<longueur_desiree)
195     {
196     refresh();
197     ti=tii;
198     dt=0.03125;
199     double param_tii[2];
200     int valide_u=0;
201     int valide_v=0;
202     int valide_xyz=0;
203     while ((valide_u & valide_v & valide_xyz)!=1)
204     {
205     tii=ti+longueur_parametrique*dt;
206     uii=umilieu+tii*(vecteur_du_front.get_x()*sin(teta)+vecteur_normal_front.get_x()*cos(teta));
207     vii=vmilieu+tii*(vecteur_du_front.get_y()*sin(teta)+vecteur_normal_front.get_y()*cos(teta));
208     param_tii[0]=uii-du;
209     param_tii[1]=vii-dv;
210     if (mgface->valide_parametre_u(param_tii[0])) valide_u=1;
211     if (mgface->valide_parametre_v(param_tii[1])) valide_v=1;
212     if ((valide_u & valide_v)!=1) dt=dt/2.;
213     if ((valide_u & valide_v)==1)
214     {
215     double xyz[3];
216     mgface->evaluer(param_tii,xyz);
217     valide_xyz=metrique->valide_parametre(xyz);
218     if ((valide_xyz)!=1) dt=dt/2.;
219     }
220     if (dt<0.0000001) return ERREUR;
221     }
222     longueur_calculee_pas_precedent=longueur_calculee;
223     double t=0.7886751345*ti+0.2113248654*tii;
224     double ut=vecteur_du_front.get_x()*sin(teta)+vecteur_normal_front.get_x()*cos(teta);
225     double vt=vecteur_du_front.get_y()*sin(teta)+vecteur_normal_front.get_y()*cos(teta);
226     double u=umilieu+t*ut;
227     double v=vmilieu+t*vt;;
228     double param_integration1[2]={u-du,v-dv};
229     double tenseur_metrique[9];
230     double xyzdu[3];
231     double xyzdv[3];
232 francois 446 double xyz[3];
233 5 mgface->evaluer(param_integration1,xyz);
234     metrique->evaluer(xyz,tenseur_metrique);
235     mgface->deriver(param_integration1,xyzdu,xyzdv);
236     double t3 = xyzdu[0]*ut+xyzdv[0]*vt;
237     double t7 = xyzdu[1]*ut+xyzdv[1]*vt;
238     double t11 = xyzdu[2]*ut+xyzdv[2]*vt;
239     double t13 = t3*tenseur_metrique[0]+t7*tenseur_metrique[3]+t11*tenseur_metrique[6];
240     double t18 = t3*tenseur_metrique[1]+t7*tenseur_metrique[4]+t11*tenseur_metrique[7];
241     double t23 = t3*tenseur_metrique[2]+t7*tenseur_metrique[5]+t11*tenseur_metrique[8];
242     double facteur = (t13*xyzdu[0]+t18*xyzdu[1]+t23*xyzdu[2])*ut+(t13*xyzdv[0]+t18*xyzdv[1]+t23*xyzdv[2])*vt;
243     longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(facteur);
244     t=0.7886751345*tii+0.2113248654*ti;
245     u=umilieu+t*ut;
246     v=vmilieu+t*vt;;
247     double param_integration2[2]={u-du,v-dv};
248     mgface->evaluer(param_integration2,xyz);
249     metrique->evaluer(xyz,tenseur_metrique);
250     mgface->deriver(param_integration2,xyzdu,xyzdv);
251     t3 = xyzdu[0]*ut+xyzdv[0]*vt;
252     t7 = xyzdu[1]*ut+xyzdv[1]*vt;
253     t11 = xyzdu[2]*ut+xyzdv[2]*vt;
254     t13 = t3*tenseur_metrique[0]+t7*tenseur_metrique[3]+t11*tenseur_metrique[6];
255     t18 = t3*tenseur_metrique[1]+t7*tenseur_metrique[4]+t11*tenseur_metrique[7];
256     t23 = t3*tenseur_metrique[2]+t7*tenseur_metrique[5]+t11*tenseur_metrique[8];
257     facteur = (t13*xyzdu[0]+t18*xyzdu[1]+t23*xyzdu[2])*ut+(t13*xyzdv[0]+t18*xyzdv[1]+t23*xyzdv[2])*vt;
258     longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(facteur);
259     }
260     double ut=vecteur_du_front.get_x()*sin(teta)+vecteur_normal_front.get_x()*cos(teta);
261     double vt=vecteur_du_front.get_y()*sin(teta)+vecteur_normal_front.get_y()*cos(teta);
262     double t=ti+(tii-ti)*(longueur_desiree-longueur_calculee_pas_precedent)/(longueur_calculee-longueur_calculee_pas_precedent);
263     double udecale=umilieu+t*ut;
264     double vdecale=vmilieu+t*vt;
265     double u=decalage->decalage_parametre_u(udecale,-du);
266     double v=decalage->decalage_parametre_v(vdecale,-dv);
267     double param_noeud_cree[2]={u,v};
268    
269    
270     // recherhce des element du front proche de u,v;
271     TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
272    
273     double umin=std::min(u1,u2);
274     umin=std::min(umin,udecale);
275     double umax=std::max(u1,u2);
276     umax=std::max(umax,udecale);
277     double vmin=std::min(v1,v2);
278     vmin=std::min(vmin,vdecale);
279     double vmax=std::max(v1,v2);
280     vmax=std::max(vmax,vdecale);
281     BOITE_2D boite_recherche(umin,vmin,umax,vmax);
282     double longueur_recherche=boite_recherche.get_rayon()*2.;
283    
284     quadtree_de_front->rechercher(u,v,longueur_recherche,liste_trouvee);
285     int nb_entite=liste_trouvee.get_nb();
286    
287     double distance_reference = -1.;
288     double angle_reference=2.*M_PI;
289     unsigned long id_noeud_reference = 0;
290     MG_FRONT_2D* front_reference=NULL;
291    
292     for (int i=0;i<nb_entite;i++)
293     {
294     MG_FRONT_2D* front_courant=liste_trouvee.get(i);
295     MG_NOEUD* noeud_front1=front_courant->get_noeud1();
296     MG_NOEUD* noeud_front2=front_courant->get_noeud2();
297     double unoeudfront1=decalage->decalage_parametre_u(noeud_front1->get_u(),du);
298     double vnoeudfront1=decalage->decalage_parametre_v(noeud_front1->get_v(),dv);
299     double unoeudfront2=decalage->decalage_parametre_u(noeud_front2->get_u(),du);
300     double vnoeudfront2=decalage->decalage_parametre_v(noeud_front2->get_v(),dv);
301     double distance_noeudfront1=calcule_distance_metrique(mgface,udecale,vdecale,unoeudfront1,vnoeudfront1,du,dv);
302     if ((distance_noeudfront1<distance_reference) || (distance_reference<0.) || (id_noeud_reference==noeud_front1->get_id()))
303     if (noeud_front1->get_id()!=noeud1->get_id())
304     if (noeud_front1->get_id()!=noeud2->get_id())
305     if (id_noeud_reference==noeud_front1->get_id())
306     {
307     refresh();
308     OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
309     OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
310     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
311     vecteur_baseu.norme();
312     vecteur_basev.norme();
313     vecteur_front.norme();
314     double cosangle=vecteur_baseu*vecteur_front;
315     double sinangle=vecteur_basev*vecteur_front;
316     sinangle=-sinangle;
317     if (cosangle>1.) cosangle=1.;
318     if (cosangle<-1.) cosangle=(-1.);
319     double angle=acos(cosangle);
320     if (sinangle<(-0.0001)) angle=(-angle);
321     if (angle<0.) angle=angle+2.*M_PI;
322     if (angle<angle_reference)
323     {
324     angle_reference=angle;
325     front_reference=front_courant;
326     }
327     }
328     else
329     {
330     refresh();
331     OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
332     OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
333     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
334     vecteur_baseu.norme();
335     vecteur_basev.norme();
336     vecteur_front.norme();
337     double cosangle=vecteur_baseu*vecteur_front;
338     double sinangle=vecteur_basev*vecteur_front;
339     sinangle=-sinangle;
340     if (cosangle>1.) cosangle=1.;
341     if (cosangle<-1.) cosangle=(-1.);
342     double angle=acos(cosangle);
343     if (sinangle<(-0.0001)) angle=(-angle);
344     if (angle<0.) angle=angle+2.*M_PI;
345     distance_reference=distance_noeudfront1;
346     angle_reference=angle;
347     id_noeud_reference=noeud_front1->get_id();
348     front_reference=front_courant;
349     }
350     }
351    
352     if ((front_reference!=NULL) && (distance_reference<longueur_desiree*2./sqrt(3)))
353     {
354     MG_NOEUD* noeud=mg_maillage->get_mg_noeudid(id_noeud_reference);
355     double uref=decalage->decalage_parametre_u(noeud->get_u(),du);
356     double vref=decalage->decalage_parametre_v(noeud->get_v(),dv);
357     double distance=calcule_distance_metrique(mgface,umilieu,vmilieu,uref,vref,du,dv);
358     if (distance<1.5*longueur_desiree)
359     {
360     (*front_rencontre)=front_reference;
361     (*noeud_rencontre)=noeud;
362     return FRONT_RENCONTRE;
363     }
364     }
365     double coord[3];
366     mgface->evaluer(param_noeud_cree,coord);
367 francois 35 MG_NOEUD* noeud=mg_maillage->ajouter_mg_noeud(mgface,coord[0],coord[1],coord[2],MAILLEUR_AUTO);
368 5 noeud->change_u(param_noeud_cree[0]);
369     noeud->change_v(param_noeud_cree[1]);
370     (*noeud_rencontre)=noeud;
371     (*front_rencontre)=NULL;
372     return NOEUD_CREE;
373     }
374    
375    
376     double MAILLEUR2D::calcule_longueur_segment_metrique(MG_FACE* mgface,MG_SEGMENT* mgsegment)
377     {
378     MG_NOEUD* noeud1=mgsegment->get_noeud1();
379     MG_NOEUD* noeud2=mgsegment->get_noeud2();
380    
381     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
382     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
383     double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
384     double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
385     double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
386     double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
387     return calcule_distance_metrique(mgface,u1,v1,u2,v2,du,dv);
388     }
389    
390    
391     double MAILLEUR2D::calcule_distance_metrique(MG_FACE* mgface,double u1,double v1,double u2,double v2,double du,double dv)
392     {
393     double longueur_calculee=0.;
394     double param_debut[2]={u1-du,v1-dv};
395     double E,F,G;
396     mgface->get_EFG(param_debut,E,F,G);
397    
398     double dt=0.03125;
399     double ti;
400     double tii=0;
401    
402     while (tii<1.)
403     {
404     refresh();
405     ti=tii;
406     tii=ti+dt;
407     if (tii>1.) tii=1.;
408     double t=0.7886751345*ti+0.2113248654*tii;
409     double u=u1+t*(u2-u1);
410     double v=v1+t*(v2-v1);
411     double ut=u2-u1;
412     double vt=v2-v1;
413     double param_integration1[2]={u-du,v-dv};
414     double tenseur_metrique[9];
415     double xyzdu[3];
416     double xyzdv[3];
417     double xyz[3];
418     mgface->evaluer(param_integration1,xyz);
419 francois 331 if (metrique->valide_parametre(xyz)) metrique->evaluer(xyz,tenseur_metrique);
420     else return 1e308;
421 francois 446 mgface->deriver(param_integration1,xyzdu,xyzdv);
422 5 double t3 = xyzdu[0]*ut+xyzdv[0]*vt;
423     double t7 = xyzdu[1]*ut+xyzdv[1]*vt;
424     double t11 = xyzdu[2]*ut+xyzdv[2]*vt;
425     double t13 = t3*tenseur_metrique[0]+t7*tenseur_metrique[3]+t11*tenseur_metrique[6];
426     double t18 = t3*tenseur_metrique[1]+t7*tenseur_metrique[4]+t11*tenseur_metrique[7];
427     double t23 = t3*tenseur_metrique[2]+t7*tenseur_metrique[5]+t11*tenseur_metrique[8];
428     double facteur = (t13*xyzdu[0]+t18*xyzdu[1]+t23*xyzdu[2])*ut+(t13*xyzdv[0]+t18*xyzdv[1]+t23*xyzdv[2])*vt;
429     longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(facteur);
430     t=0.7886751345*tii+0.2113248654*ti;
431     u=u1+t*(u2-u1);
432     v=v1+t*(v2-v1);
433     double param_integration2[2]={u-du,v-dv};
434     mgface->evaluer(param_integration2,xyz);
435 francois 331 if (metrique->valide_parametre(xyz)) metrique->evaluer(xyz,tenseur_metrique);
436     else return 1e308;
437 5 mgface->deriver(param_integration2,xyzdu,xyzdv);
438     t3 = xyzdu[0]*ut+xyzdv[0]*vt;
439     t7 = xyzdu[1]*ut+xyzdv[1]*vt;
440     t11 = xyzdu[2]*ut+xyzdv[2]*vt;
441     t13 = t3*tenseur_metrique[0]+t7*tenseur_metrique[3]+t11*tenseur_metrique[6];
442     t18 = t3*tenseur_metrique[1]+t7*tenseur_metrique[4]+t11*tenseur_metrique[7];
443     t23 = t3*tenseur_metrique[2]+t7*tenseur_metrique[5]+t11*tenseur_metrique[8];
444     facteur = (t13*xyzdu[0]+t18*xyzdu[1]+t23*xyzdu[2])*ut+(t13*xyzdv[0]+t18*xyzdv[1]+t23*xyzdv[2])*vt;
445     longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(facteur);
446     }
447     return longueur_calculee;
448     }
449    
450    
451    
452    
453    
454    
455     MG_TRIANGLE* MAILLEUR2D::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3)
456     {
457     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
458     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
459     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
460 francois 35 if (mgsegment1==NULL) mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,MAILLEUR_AUTO);
461     if (mgsegment2==NULL) mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,MAILLEUR_AUTO);
462     if (mgsegment3==NULL) mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,MAILLEUR_AUTO);
463 francois 54 M3D_TRIANGLE* mtriangle=new M3D_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAILLEUR_AUTO);
464 5 mg_maillage->ajouter_mg_triangle(mtriangle);
465     double qual=OPERATEUR::qualite_triangle(mgnoeud1->get_coord(),mgnoeud2->get_coord(),mgnoeud3->get_coord());
466     mtriangle->change_qualite(qual);
467     std::pair<const double,M3D_TRIANGLE*> tmp(mtriangle->get_qualite(),mtriangle);
468     return mtriangle;
469     }
470    
471     int MAILLEUR2D::triangle_est_dans_bon_sens(MG_FACE* face,MG_NOEUD* noeud1,MG_NOEUD *noeud2,MG_NOEUD *noeud3)
472     {
473     double* xyz1=noeud1->get_coord();
474     double* xyz2=noeud2->get_coord();
475     double* xyz3=noeud3->get_coord();
476     OT_VECTEUR_3D n1n3(xyz1,xyz3);
477     OT_VECTEUR_3D n1n2(xyz1,xyz2);
478     double uv[2];
479     face->inverser(uv,xyz1);
480     double n[3];
481     face->calcul_normale_unitaire(uv,n);
482     OT_VECTEUR_3D nmat=n&n1n3;
483     nmat.norme();
484     n1n2.norme();
485     if (nmat*n1n2<0) return 0;
486     return 1;
487     }
488    
489    
490    
491    
492    
493    
494    
495     // PARTIE DE CODE RECUPERER DE MON PROJET DE DEA
496     // TRES ILLISIBLE MAIS FONCTIONNEL
497     // POIL AU MARRON
498    
499     #define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
500     #define EGAL(x,y,eps) (float)fabs((double)(x-y))<eps
501     #define DETER(a,b,c,d) (a*d-b*c)
502    
503     int MAILLEUR2D::intersection_segment_segment(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
504     {
505    
506     double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
507     double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
508     double ua=decalage->decalage_parametre_u(noeud1->get_u(),du);
509     double va=decalage->decalage_parametre_v(noeud1->get_v(),dv);
510     double ub=decalage->decalage_parametre_u(noeud2->get_u(),du);
511     double vb=decalage->decalage_parametre_v(noeud2->get_v(),dv);
512     double um=decalage->decalage_parametre_u(noeud3->get_u(),du);
513     double vm=decalage->decalage_parametre_v(noeud3->get_v(),dv);
514     double un=decalage->decalage_parametre_u(noeud4->get_u(),du);
515     double vn=decalage->decalage_parametre_v(noeud4->get_v(),dv);
516     double ab[3];
517     double nm[3];
518     double am[3];
519     ab[0]=ub-ua;
520     ab[1]=vb-va;
521     ab[2]=0.;
522     nm[0]=um-un;
523     nm[1]=vm-vn;
524     nm[2]=0.;
525     am[0]=um-ua;
526     am[1]=vm-va;
527     am[2]=0.;
528     int equation[4];
529     equation[0]=1; /* etat de l'equation 0 */
530     equation[1]=1;
531     equation[2]=1;
532     equation[3]=3; /* cette variable comporte le bilan du nombre d'equation */
533     double eps2=PSCA(ab,ab);
534     double eps=sqrt(eps2);
535     eps=eps*0.0001;
536     eps2=eps2*0.0001;
537     /* recherche du nombre d'equation -> inter franche ou para ou confondu */
538     if ( (EGAL(ab[0],0,eps)) && (EGAL(nm[0],0,eps)) )
539     if (EGAL(am[0],0,eps)) equation[0]=0; else return(0);
540     if ( (EGAL(ab[1],0,eps)) && (EGAL(nm[1],0,eps)) )
541     if (EGAL(am[1],0,eps)) equation[1]=0; else return(0);
542     if ( (EGAL(ab[2],0,eps)) && (EGAL(nm[2],0,eps)) )
543     if (EGAL(am[2],0,eps)) equation[2]=0; else return(0);
544     equation[3]=equation[0]+equation[1]+equation[2];
545     if (equation[3]==3)
546     {
547     double det=DETER(ab[0],nm[0],ab[1],nm[1]);
548     if (fabs(det)>eps2)
549     {
550     det=1/det;
551     double sol1=det*DETER(am[0],nm[0],am[1],nm[1]);
552     double sol2=det*DETER(ab[0],am[0],ab[1],am[1]);
553     if ( (float)fabs((double)(sol1*ab[2]-sol2*nm[2]-am[2]))>eps2) return(0);
554     return(examine_solution(sol1,sol2,1));
555     }
556     else
557     {
558     equation[0]=0;
559     equation[3]=2;
560     /* on verifie la compatibilite des deux equations dont le det est nul*/
561     double tmp;
562     if (ab[0]!=0) tmp=ab[1]*am[0]/ab[0]; else tmp=nm[1]*am[0]/nm[0];
563     if (!(EGAL(tmp,am[1],eps))) return(0);
564     }
565     }
566     if (equation[3]==2)
567     {
568     /* on repere les equations qui existent */
569     int ne1;
570     int ne2;
571     if (equation[0]!=0)
572     {
573     ne1=0;
574     if (equation[1]!=0) ne2=1; else ne2=2;
575     }
576     else
577     {
578     ne1=1;
579     ne2=2;
580     }
581    
582     double det=DETER(ab[ne1],nm[ne1],ab[ne2],nm[ne2]);
583     if (fabs(det)>eps2)
584     {
585     det=1/det;
586     double sol1=det*DETER(am[ne1],nm[ne1],am[ne2],nm[ne2]);
587     double sol2=det*DETER(ab[ne1],am[ne1],ab[ne2],am[ne2]);
588     return(examine_solution(sol1,sol2,1));
589     }
590     else
591     {
592     equation[ne1]=0;
593     equation[3]=1;
594     /* on verifie la compatibilite des deux equations dont le det est nul */
595     double tmp;
596     if (ab[ne1]!=0) tmp=ab[ne2]*am[ne1]/ab[ne1]; else tmp=nm[ne2]*am[ne1]/nm[ne1];
597     if (!(EGAL(tmp,am[ne2],eps))) return(0);
598     }
599    
600     }
601     if (equation[3]==1)
602     {
603     /* on repere l' equation qui existe */
604     int ne1;
605     if (equation[0]!=0) ne1=0; else
606     if (equation[1]!=0) ne1=1; else ne1=2;
607     double an[3];
608     an[0]=un-ua;
609     an[1]=vn-va;
610     an[2]=0.;
611     double tmp=1/ab[ne1];
612     double sol1=am[ne1]*tmp;
613     double sol2=an[ne1]*tmp;
614     return(examine_solution(sol1,sol2,2));
615     }
616     return(0);
617     }
618    
619    
620    
621     int MAILLEUR2D::examine_solution(double sol1,double sol2,int type)
622     {
623     double epsilon=0.0001;
624    
625     if (type==1)
626     {
627     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
628     if ( ( (EGAL(sol1,0,epsilon)) || (EGAL(sol1,1,epsilon))) && ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) ) return 1;
629     if ( ( (EGAL(sol2,0,epsilon)) || (EGAL(sol2,1,epsilon))) && ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) ) return 1;
630     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>(-0.1-epsilon)) && ((sol2)<(1.1-epsilon)) ) return 1;
631     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) && (sol1>(-0.1-epsilon)) && ((sol1)<(1.1-epsilon)) ) return 1;
632    
633     }
634     if (type==2)
635     {
636     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) return 1;
637     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
638     if ( ((sol1)>(1+epsilon)) && ((-sol2)>epsilon) ) return 1;
639     if ( ((sol2)>(1+epsilon)) && ((-sol1)>epsilon) ) return 1;
640     }
641     return 0;
642     }
643     #undef EGAL
644     #undef PSCA
645     #undef DETER
646     ///FIN DU CODE DE DEA