9 R areteideale_( R3 xyz, R3 direction )
14 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
15 //dans la direction donnee
16 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
19 static double cpunew, cpuold=0;
21 void tempscpu_( double & tempsec )
22 //Retourne le temps CPU utilise en secondes
24 tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
25 //MESSAGE( "temps cpu=" << tempsec );
29 void deltacpu_( R & dtcpu )
30 //Retourne le temps CPU utilise en secondes depuis le precedent appel
33 dtcpu = R( cpunew - cpuold );
35 //MESSAGE( "delta temps cpu=" << dtcpu );
40 void aptrte( Z nutysu, R aretmx,
41 Z nblf, Z * nudslf, R2 * uvslf,
43 Z & nbst, R2 * & uvst, Z & nbt, Z * & nust,
45 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
46 // but : appel de la triangulation par un arbre-4 recouvrant
47 // ----- de triangles equilateraux
48 // le contour du domaine plan est defini par des lignes fermees
49 // la premiere ligne etant l'enveloppe de toutes les autres
50 // la fonction areteideale(s,d) donne la taille d'arete
51 // au point s dans la direction (actuellement inactive) d
52 // des lors toute arete issue d'un sommet s devrait avoir une longueur
53 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
56 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
57 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
61 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
62 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
63 // 1 il existe une fonction areteideale_(s,d)
64 // dont seules les 2 premieres composantes de uv sont actives
65 // ... autres options a definir ...
66 // aretmx : longueur maximale des aretes de la future triangulation
67 // nblf : nombre de lignes fermees de la surface
68 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
69 // nudslf(0)=0 pour permettre la difference sans test
70 // Attention le dernier sommet de chaque ligne est raccorde au premier
71 // tous les sommets et les points internes ont des coordonnees
72 // UV differentes <=> Pas de point double!
73 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
74 // nbpti : nombre de points internes futurs sommets de la triangulation
75 // uvpti : uv des points internes futurs sommets de la triangulation
79 // nbst : nombre de sommets de la triangulation finale
80 // uvst : coordonnees uv des nbst sommets de la triangulation
81 // nbt : nombre de triangles de la triangulation finale
82 // nust : 4 numeros dans uvst des sommets des nbt triangles
83 // s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
84 // ierr : 0 si pas d'erreur
86 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 // auteur : Alain Perronnet Analyse Numerique Paris UPMC decembre 2001
88 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
91 R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
92 Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
93 Z mxtrou = Max( 1024, nblf ); //nombre maximal de trous dans la surface
96 Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
97 Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
98 Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
99 Z *mnqueu=NULL, mxqueu;
101 Z *mnarcf=NULL, mxarcf;
111 R3 comxmi[2]; //coordonnees UV Min et Maximales
112 R aremin, aremax; //longueur minimale et maximale des aretes
116 Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
117 Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
120 aretemaxface_ = aretmx;
122 // initialisation du temps cpu
126 // quelques reservations de tableaux pour faire les calculs
127 // ========================================================
128 // le tableau pointeur sur la premiere arete de chaque ligne fermee
129 if( mndalf!=NULL ) delete [] mndalf;
130 mndalf = new Z[1+nblf];
131 if( mndalf==NULL ) goto ERREUR;
134 // declaration du tableau des coordonnees des sommets de la frontiere
135 // puis des sommets internes ajoutes
136 // majoration empirique du nombre de sommets de la triangulation
138 mxsomm = Max( 20000, 64*nbpti+i*i );
139 MESSAGE( "APTRTE: Depart de la triangulation avec " );
140 MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx << " mxsomm=" << mxsomm );
143 //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
144 if( mnpxyd!=NULL ) delete [] mnpxyd;
145 mnpxyd = new R3[mxsomm];
146 if( mnpxyd==NULL ) goto ERREUR;
148 // le tableau mnsoar des aretes des triangles
149 // 1: sommet 1 dans pxyd,
150 // 2: sommet 2 dans pxyd,
151 // 3: numero de 1 a nblf de la ligne qui supporte l'arete
152 // 4: numero dans mnartr du triangle 1 partageant cette arete,
153 // 5: numero dans mnartr du triangle 2 partageant cette arete,
154 // 6: chainage des aretes frontalieres ou internes ou
155 // des aretes simples des etoiles de triangles,
156 // 7: chainage du hachage des aretes
157 // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
158 // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
159 // h(ns1,ns2) = min( ns1, ns2 )
160 if( mnsoar!=NULL ) delete [] mnsoar;
161 mxsoar = 3 * ( mxsomm + mxtrou );
162 mnsoar = new Z[mosoar*mxsoar];
163 if( mnsoar==NULL ) goto ERREUR;
164 //initialiser le tableau mnsoar pour le hachage des aretes
165 insoar_( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
167 // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
168 if( mnarst!=NULL ) delete [] mnarst;
169 mnarst = new Z[1+mxsomm];
170 if( mnarst==NULL ) goto ERREUR;
172 azeroi_( n, mnarst );
174 // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
175 // ou no du point si interne forc'e par l'utilisateur
176 // ou 0 si interne cree par le module
177 if( mnslig!=NULL ) delete [] mnslig;
178 mnslig = new Z[mxsomm];
179 if( mnslig==NULL ) goto ERREUR;
180 azeroi_( mxsomm, mnslig );
182 // initialisation des aretes frontalieres de la triangulation future
183 // renumerotation des sommets des aretes des lignes pour la triangulation
184 // mise a l'echelle des coordonnees des sommets pour obtenir une
185 // meilleure precision lors des calculs + quelques verifications
186 // boucle sur les lignes fermees qui forment la frontiere
187 // ======================================================================
192 for (n=1; n<=nblf; n++)
194 //l'initialisation de la premiere arete de la ligne n dans la triangulation
195 //-------------------------------------------------------------------------
196 //le sommet ns0 est le numero de l'origine de la ligne
198 mnpxyd[ns0].x = uvslf[ns0].x;
199 mnpxyd[ns0].y = uvslf[ns0].y;
200 mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction );
201 // cout << "Sommet " << ns0 << ": " << mnpxyd[ns0].x
202 // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z << endl;
204 //carre de la longueur de l'arete 1 de la ligne fermee n
205 d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
206 + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
207 aremin = Min( aremin, d );
208 aremax = Max( aremax, d );
210 //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
211 //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
212 //le numero des 2 sommets ns1 ns2 de la 1-ere arete
213 //Attention: les numeros ns debutent a 1 (ils ont >0)
214 // les tableaux c++ demarrent a zero!
215 // les tableaux fortran demarrent ou l'on veut!
220 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
221 mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
222 fasoar_( ns1, ns2, moins1, moins1, n,
223 mosoar, mxsoar, n1soar, mnsoar, mnarst,
225 //pas de test sur ierr car pas de saturation possible a ce niveau
227 //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
230 //la nouvelle arete est la suivante de l'arete definie juste avant
232 mnsoar[mosoar * noar - mosoar + 5] = noar0;
234 //l'initialisation des aretes suivantes de la ligne dans la triangulation
235 //-----------------------------------------------------------------------
236 nbarli = nudslf[n] - nudslf[n-1]; //nombre d'aretes=sommets de la ligne n
237 for (i=2; i<=nbarli; i++)
239 ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
241 //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
244 //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
247 //l'arete precedente est dotee de sa suivante:celle cree ensuite
248 //les 2 coordonnees du sommet ns2 de la ligne
250 mnpxyd[ns].x = uvslf[ns].x;
251 mnpxyd[ns].y = uvslf[ns].y;
252 mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction );
253 // cout << "Sommet " << ns << ": " << mnpxyd[ns].x
254 // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z << endl;
256 //carre de la longueur de l'arete
257 d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
258 + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
259 aremin = Min( aremin, d );
260 aremax = Max( aremax, d );
262 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
263 mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
265 //ajout de l'arete dans la liste
266 fasoar_( ns1, ns2, moins1, moins1, n,
267 mosoar, mxsoar, n1soar, mnsoar,
268 mnarst, noar, ierr );
269 //pas de test sur ierr car pas de saturation possible a ce niveau
271 //chainage des aretes frontalieres en position 6 du tableau mnsoar
272 //la nouvelle arete est la suivante de l'arete definie juste avant
273 mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
276 //attention: la derniere arete de la ligne fermee enveloppe
277 // devient en fait la premiere arete de cette ligne
278 // dans le chainage des aretes de la frontiere!
280 if( ierr != 0 ) goto ERREUR;
282 aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
283 aremax = sqrt( aremax ); //longueur maximale d'une arete
285 aretmx = Min( aretmx, aremax ); //pour homogeneiser
286 cout << "nutysu=" << nutysu << " aretmx=" << aretmx
287 << " arete min=" << aremin << " arete max=" << aremax << endl;
289 //chainage des aretes frontalieres : la derniere arete frontaliere
290 mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
292 //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
293 //reservation du tableau des numeros des 3 aretes de chaque triangle
294 //mnartr( moartr, mxartr )
295 //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
296 // 3Triangles = 2 Aretes internes + Aretes frontalieres
297 // d'ou 3T/2 < AI + AF => T < 3T/2 - Sommets + 1-Trous
298 //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
299 if( mnartr!=NULL ) delete [] mnartr;
300 mxartr = 2 * ( mxsomm + mxtrou );
301 mnartr = new Z[moartr*mxartr];
302 if( mnartr==NULL ) goto ERREUR;
304 //Ajout des points internes
305 ns1 = nudslf[ nblf ];
306 for (i=0; i<nbpti; i++)
308 //les 2 coordonnees du point i de sommet nbs
309 mnpxyd[ns1].x = uvpti[i].x;
310 mnpxyd[ns1].y = uvpti[i].y;
311 mnpxyd[ns1].z = areteideale_( mnpxyd[ns1], direction );
312 //le numero i du point interne
317 //nombre de sommets de la frontiere et internes
320 // creation de l'arbre-4 des te (tableau letree)
321 // ajout dans les te des sommets des lignes et des points internes imposes
322 // =======================================================================
323 // premiere estimation de mxtree
326 NEWTREE: //en cas de saturation de l'un des tableaux, on boucle
327 MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
328 if( mntree != NULL ) delete [] mntree;
330 mntree = new Z[motree*(1+mxtree)];
331 if( mntree==NULL ) goto ERREUR;
333 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
334 teajte_( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
340 //saturation de letree => sa taille est augmentee et relance
343 MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
349 MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
350 if( ierr != 0 ) goto ERREUR;
351 //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
353 // homogeneisation de l'arbre des te a un saut de taille au plus
354 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
355 // ===========================================================================
356 // reservation de la queue pour parcourir les te de l'arbre
357 if( mnqueu != NULL ) delete [] mnqueu;
359 mnqueu = new Z[mxqueu];
360 if( mnqueu==NULL) goto ERREUR;
362 tehote_( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
364 mntree, mxqueu, mnqueu,
369 cout << "Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
370 << d << " secondes" << endl;
373 //destruction du tableau auxiliaire et de l'arbre
378 MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
386 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
387 // et des points de la frontiere, des points internes imposes interieurs
388 // ==========================================================================
389 tetrte_( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
390 mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
391 moartr, mxartr, n1artr, mnartr, mnarst,
394 // destruction de la queue et de l'arbre devenus inutiles
395 delete [] mnqueu; mnqueu=NULL;
396 delete [] mntree; mntree=NULL;
401 MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
403 // ierr =0 si pas d'erreur
404 // =1 si le tableau mnsoar est sature
405 // =2 si le tableau mnartr est sature
406 // =3 si aucun des triangles ne contient l'un des points internes
407 // =5 si saturation de la queue de parcours de l'arbre des te
408 if( ierr != 0 ) goto ERREUR;
410 //qualites de la triangulation actuelle
411 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
412 nbt, quamoy, quamin );
414 // boucle sur les aretes internes (non sur une ligne de la frontiere)
415 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
416 // ======================================================================
417 // formation du chainage 6 des aretes internes a echanger eventuellement
418 aisoar_( mosoar, mxsoar, mnsoar, na );
419 tedela_( mnpxyd, mnarst,
420 mosoar, mxsoar, n1soar, mnsoar, na,
421 moartr, mxartr, n1artr, mnartr, n );
423 MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
426 cout << "Temps de la triangulation Delaunay par echange des diagonales="
427 << d << " secondes" << endl;
429 //qualites de la triangulation actuelle
430 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
431 nbt, quamoy, quamin );
433 // detection des aretes frontalieres initiales perdues
434 // triangulation frontale pour les restaurer
435 // ===================================================
437 if( mn1arcf != NULL ) delete [] mn1arcf;
438 if( mnarcf != NULL ) delete [] mnarcf;
439 if( mnarcf1 != NULL ) delete [] mnarcf1;
440 if( mnarcf2 != NULL ) delete [] mnarcf2;
441 mn1arcf = new Z[1+mxarcf];
442 if( mn1arcf == NULL ) goto ERREUR;
443 mnarcf = new Z[3*mxarcf];
444 if( mnarcf == NULL ) goto ERREUR;
445 mnarcf1 = new Z[mxarcf];
446 if( mnarcf1 == NULL ) goto ERREUR;
447 mnarcf2 = new Z[mxarcf];
448 if( mnarcf2 == NULL ) goto ERREUR;
450 terefr_( nbarpi, mnpxyd,
451 mosoar, mxsoar, n1soar, mnsoar,
452 moartr, n1artr, mnartr, mnarst,
453 mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
456 MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" );
459 cout << "Temps de la recuperation des aretes perdues de la frontiere="
460 << d << " secondes" << endl;
462 if( ierr != 0 ) goto ERREUR;
464 //qualites de la triangulation actuelle
465 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
466 nbt, quamoy, quamin );
468 // fin de la triangulation avec respect des aretes initiales frontalieres
470 // suppression des triangles externes a la surface
471 // ===============================================
472 // recherche du dernier triangle utilise
473 mn = mxartr * moartr;
474 for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
477 if( mnartr[mn] != 0 ) break;
480 if( mntrsu != NULL ) delete [] mntrsu;
481 mntrsu = new Z[ndtri0];
482 if( mntrsu == NULL ) goto ERREUR;
484 if( mnlftr != NULL ) delete [] mnlftr;
485 mnlftr = new Z[nblf];
486 if( mnlftr == NULL ) goto ERREUR;
488 for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
491 tesuex_( nblf, mnlftr,
492 ndtri0, nbsomm, mnpxyd, mnslig,
493 mosoar, mxsoar, mnsoar,
494 moartr, mxartr, n1artr, mnartr, mnarst,
497 delete [] mnlftr; mnlftr=NULL;
498 delete [] mntrsu; mntrsu=NULL;
502 MESSAGE( "Temps de la suppression des triangles externes=" << d );
503 if( ierr != 0 ) goto ERREUR;
505 //qualites de la triangulation actuelle
506 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
507 nbt, quamoy, quamin );
509 // amelioration de la qualite de la triangulation par
510 // barycentrage des sommets internes a la triangulation
511 // suppression des aretes trop longues ou trop courtes
512 // modification de la topologie des groupes de triangles
513 // mise en delaunay de la triangulation
514 // =====================================================
515 mnarcf3 = new Z[mxarcf];
516 if( mnarcf3 == NULL ) goto ERREUR;
519 mnarst, mosoar, mxsoar, n1soar, mnsoar,
520 moartr, mxartr, n1artr, mnartr,
521 mxarcf, mnarcf2, mnarcf3,
522 mn1arcf, mnarcf, mnarcf1,
523 comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
525 if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
526 if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
527 if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
528 if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
529 if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
533 MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
534 if( ierr != 0 ) goto ERREUR;
536 //qualites de la triangulation finale
537 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
538 nbt, quamoy, quamin );
540 // renumerotation des sommets internes: mnarst(i)=numero final du sommet
541 // ===================================
542 for (i=0; i<=nbsomm; i++)
545 for (nt=1; nt<=mxartr; nt++)
547 if( mnartr[nt*moartr-moartr] != 0 )
549 //le numero des 3 sommets du triangle nt
550 nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
551 //les 3 sommets du triangle sont actifs
552 mnarst[ nosotr[0] ] = 1;
553 mnarst[ nosotr[1] ] = 1;
554 mnarst[ nosotr[2] ] = 1;
558 for (i=1; i<=nbsomm; i++)
564 // generation du tableau uvst de la surface triangulee
565 // ---------------------------------------------------
566 if( uvst != NULL ) delete [] uvst;
568 if( uvst == NULL ) goto ERREUR;
571 for (i=0; i<nbsomm; i++ )
576 uvst[nbst].x = mnpxyd[i].x;
577 uvst[nbst].y = mnpxyd[i].y;
579 //si le sommet est un point ou appartient a une ligne
580 //ses coordonnees initiales sont restaurees
587 //retour aux coordonnees initiales dans uvslf
589 n = n - 1000000 * l + nudslf[l-1] - 1;
590 uvst[nbst].x = uvslf[n].x;
591 uvst[nbst].y = uvslf[n].y;
595 //point utilisateur n interne impose
596 //retour aux coordonnees initiales dans uvpti
597 uvst[nbst].x = uvpti[n-1].x;
598 uvst[nbst].y = uvpti[n-1].y;
605 // generation du tableau 'nsef' de la surface triangulee
606 // -----------------------------------------------------
607 // boucle sur les triangles occupes (internes et externes)
608 if( nust != NULL ) delete [] nust;
610 if( nust == NULL ) goto ERREUR;
612 for (i=1; i<=mxartr; i++)
614 //le triangle i de mnartr
615 if( mnartr[i*moartr-moartr] != 0 )
617 //le triangle i est interne => nosotr numero de ses 3 sommets
618 nusotr_( i, mosoar, mnsoar, moartr, mnartr, nosotr );
619 nust[nbt++] = mnarst[ nosotr[0] ];
620 nust[nbt++] = mnarst[ nosotr[1] ];
621 nust[nbt++] = mnarst[ nosotr[2] ];
625 nbt /= 4; //le nombre final de triangles de la surface
626 cout << "Nombre de sommets=" << nbst
627 << " Nombre de triangles=" << nbt << endl;
631 MESSAGE( "Temps total de la triangulation=" << tcpu << " secondes" );
633 // destruction des tableaux auxiliaires
634 // ------------------------------------
636 if( mnarst != NULL ) delete [] mnarst;
637 if( mnartr != NULL ) delete [] mnartr;
638 if( mnslig != NULL ) delete [] mnslig;
639 if( mnsoar != NULL ) delete [] mnsoar;
640 if( mnpxyd != NULL ) delete [] mnpxyd;
641 if( mndalf != NULL ) delete [] mndalf;
642 if( mntree != NULL ) delete [] mntree;
643 if( mnqueu != NULL ) delete [] mnqueu;
644 if( mntrsu != NULL ) delete [] mntrsu;
645 if( mnlftr != NULL ) delete [] mnlftr;
646 if( mn1arcf != NULL ) delete [] mn1arcf;
647 if( mnarcf != NULL ) delete [] mnarcf;
648 if( mnarcf1 != NULL ) delete [] mnarcf1;
649 if( mnarcf2 != NULL ) delete [] mnarcf2;
650 if( mnarcf3 != NULL ) delete [] mnarcf3;
654 if( ierr == 51 || ierr == 52 )
656 //saturation des sommets => redepart avec 2 fois plus de sommets
663 MESSAGE( "Triangulation non realisee " << ierr );
664 if( ierr == 0 ) ierr=1;
670 void qualitetrte( R3 *mnpxyd,
671 Z & mosoar, Z & mxsoar, Z *mnsoar,
672 Z & moartr, Z & mxartr, Z *mnartr,
673 Z & nbtria, R & quamoy, R & quamin )
674 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
675 // but : calculer la qualite moyenne et minimale de la triangulation
676 // ----- actuelle definie par les tableaux mnsoar et mnartr
679 // mnpxyd : tableau des coordonnees 2d des points
680 // par point : x y distance_souhaitee
681 // mosoar : nombre maximal d'entiers par arete et
682 // indice dans mnsoar de l'arete suivante dans le hachage
683 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
684 // attention: mxsoar>3*mxsomm obligatoire!
685 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
686 // chainage des aretes frontalieres, chainage du hachage des aretes
687 // hachage des aretes = mnsoar(1)+mnsoar(2)*2
688 // avec mxsoar>=3*mxsomm
689 // une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
690 // mnsoar(2,arete vide)=l'arete vide qui precede
691 // mnsoar(3,arete vide)=l'arete vide qui suit
692 // moartr : nombre maximal d'entiers par arete du tableau mnartr
693 // mxartr : nombre maximal de triangles declarables
694 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
695 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
698 // nbtria : nombre de triangles internes au domaine
699 // quamoy : qualite moyenne des triangles actuels
700 // quamin : qualite minimale des triangles actuels
701 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
704 Z nosotr[3], mn, nbtrianeg, nt;
713 for ( nt=1; nt<=mxartr; nt++ )
718 //un triangle occupe de plus
721 //le numero des 3 sommets du triangle nt
722 nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
724 //la qualite du triangle ns1 ns2 ns3
725 qutr2d_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
731 //la qualite minimale
732 quamin = Min( quamin, qualite );
734 //aire signee du triangle nt
735 d = surtd2_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
738 //un triangle d'aire negative de plus
740 cout << "ATTENTION: le triangle " << nt << " de sommets:"
741 << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
742 << " a une aire " << d <<"<=0" << endl;
745 //aire des triangles actuels
752 cout << "Qualite moyenne=" << quamoy
753 << " Qualite minimale=" << quamin
754 << " des " << nbtria << " triangles de surface totale="
758 MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );