1 // MEFISTO2: a library to compute 2D triangulation from segmented boundaries
4 // Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
22 // File : aptrte.cxx le C++ de l'appel du trianguleur plan
24 // Author : Alain PERRONNET
25 // Date : 13 novembre 2006
29 #include "utilities.h"
41 areteideale()//( R3 xyz, R3 direction )
46 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
47 //dans la direction donnee
48 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
51 static double cpunew, cpuold=0;
57 tempscpu_( double & tempsec )
58 //Retourne le temps CPU utilise en secondes
60 tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
61 //MESSAGE( "temps cpu=" << tempsec );
69 deltacpu_( R & dtcpu )
70 //Retourne le temps CPU utilise en secondes depuis le precedent appel
73 dtcpu = R( cpunew - cpuold );
75 //MESSAGE( "delta temps cpu=" << dtcpu );
80 void aptrte( Z nutysu, R aretmx,
81 Z nblf, Z * nudslf, R2 * uvslf,
83 Z & nbst, R2 * & uvst,
86 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 // but : appel de la triangulation par un arbre-4 recouvrant
88 // ----- de triangles equilateraux
89 // le contour du domaine plan est defini par des lignes fermees
90 // la premiere ligne etant l'enveloppe de toutes les autres
91 // la fonction areteideale(s,d) donne la taille d'arete
92 // au point s dans la direction (actuellement inactive) d
93 // des lors toute arete issue d'un sommet s devrait avoir une longueur
94 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
97 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
98 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
102 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
103 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
104 // 1 il existe une fonction areteideale_(s,d)
105 // dont seules les 2 premieres composantes de uv sont actives
106 // ... autres options a definir ...
107 // aretmx : longueur maximale des aretes de la future triangulation
108 // nblf : nombre de lignes fermees de la surface
109 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
110 // nudslf(0)=0 pour permettre la difference sans test
111 // Attention le dernier sommet de chaque ligne est raccorde au premier
112 // tous les sommets et les points internes ont des coordonnees
113 // UV differentes <=> Pas de point double!
114 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
115 // nbpti : nombre de points internes futurs sommets de la triangulation
116 // uvpti : uv des points internes futurs sommets de la triangulation
120 // nbst : nombre de sommets de la triangulation finale
121 // uvst : coordonnees uv des nbst sommets de la triangulation
122 // nbt : nombre de triangles de la triangulation finale
123 // nust : 4 numeros dans uvst des sommets des nbt triangles
124 // s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
125 // ierr : 0 si pas d'erreur
127 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
128 // auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006
129 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
131 Z nbsttria=4; //Attention: 4 sommets stockes par triangle
132 //no st1, st2, st3, 0 (non quadrangle)
135 R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
136 Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
137 Z mxtrou = Max( 1024, nblf ); //nombre maximal de trous dans la surface
140 Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
141 Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
142 Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
143 Z *mnqueu=NULL, mxqueu;
145 Z *mnarcf=NULL, mxarcf;
154 R3 comxmi[2]; //coordonnees UV Min et Maximales
155 R aremin, aremax; //longueur minimale et maximale des aretes
156 R airemx; //aire maximale souhaitee d'un triangle
160 Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
161 Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
165 // initialisation du temps cpu
169 // quelques reservations de tableaux pour faire les calculs
170 // ========================================================
171 // declaration du tableau des coordonnees des sommets de la frontiere
172 // puis des sommets internes ajoutes
173 // majoration empirique du nombre de sommets de la triangulation
175 mxsomm = Max( 20000, 64*nbpti+i*i );
176 MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
177 MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx
178 << " mxsomm=" << mxsomm );
179 MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
182 //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
183 if( mnpxyd!=NULL ) delete [] mnpxyd;
184 mnpxyd = new R3[mxsomm];
185 if( mnpxyd==NULL ) goto ERREUR;
187 // le tableau mnsoar des aretes des triangles
188 // 1: sommet 1 dans pxyd,
189 // 2: sommet 2 dans pxyd,
190 // 3: numero de 1 a nblf de la ligne qui supporte l'arete
191 // 4: numero dans mnartr du triangle 1 partageant cette arete,
192 // 5: numero dans mnartr du triangle 2 partageant cette arete,
193 // 6: chainage des aretes frontalieres ou internes ou
194 // des aretes simples des etoiles de triangles,
195 // 7: chainage du hachage des aretes
196 // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
197 // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
198 // h(ns1,ns2) = min( ns1, ns2 )
199 if( mnsoar!=NULL ) delete [] mnsoar;
200 mxsoar = 3 * ( mxsomm + mxtrou );
201 mnsoar = new Z[mosoar*mxsoar];
202 if( mnsoar==NULL ) goto ERREUR;
203 //initialiser le tableau mnsoar pour le hachage des aretes
204 insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
206 // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
207 if( mnarst!=NULL ) delete [] mnarst;
208 mnarst = new Z[1+mxsomm];
209 if( mnarst==NULL ) goto ERREUR;
213 // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
214 // ou no du point si interne forc'e par l'utilisateur
215 // ou 0 si interne cree par le module
216 if( mnslig!=NULL ) delete [] mnslig;
217 mnslig = new Z[mxsomm];
218 if( mnslig==NULL ) goto ERREUR;
219 azeroi( mxsomm, mnslig );
221 // initialisation des aretes frontalieres de la triangulation future
222 // renumerotation des sommets des aretes des lignes pour la triangulation
223 // mise a l'echelle des coordonnees des sommets pour obtenir une
224 // meilleure precision lors des calculs + quelques verifications
225 // boucle sur les lignes fermees qui forment la frontiere
226 // ======================================================================
231 for (n=1; n<=nblf; n++)
233 //l'initialisation de la premiere arete de la ligne n dans la triangulation
234 //-------------------------------------------------------------------------
235 //le sommet ns0 est le numero de l'origine de la ligne
237 mnpxyd[ns0].x = uvslf[ns0].x;
238 mnpxyd[ns0].y = uvslf[ns0].y;
239 mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
240 // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
241 // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
243 //carre de la longueur de l'arete 1 de la ligne fermee n
244 d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
245 + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
246 aremin = Min( aremin, d );
247 aremax = Max( aremax, d );
249 //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
250 //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
251 //le numero des 2 sommets ns1 ns2 de la 1-ere arete
252 //Attention: les numeros ns debutent a 1 (ils ont >0)
253 // les tableaux c++ demarrent a zero!
254 // les tableaux fortran demarrent ou l'on veut!
259 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
260 mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
261 fasoar( ns1, ns2, moins1, moins1, n,
262 mosoar, mxsoar, n1soar, mnsoar, mnarst,
264 //pas de test sur ierr car pas de saturation possible a ce niveau
266 //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
269 //la nouvelle arete est la suivante de l'arete definie juste avant
271 mnsoar[mosoar * noar - mosoar + 5] = noar0;
273 //l'initialisation des aretes suivantes de la ligne dans la triangulation
274 //-----------------------------------------------------------------------
275 nbarli = nudslf[n] - nudslf[n-1]; //nombre d'aretes=sommets de la ligne n
276 for (i=2; i<=nbarli; i++)
278 ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
280 //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
283 //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
286 //l'arete precedente est dotee de sa suivante:celle cree ensuite
287 //les 2 coordonnees du sommet ns2 de la ligne
289 //debut ajout 5/10/2006 ................................................
290 nuds = Max( nuds, ns ); //le numero du dernier sommet traite
291 //fin ajout 5/10/2006 ................................................
292 mnpxyd[ns].x = uvslf[ns].x;
293 mnpxyd[ns].y = uvslf[ns].y;
294 mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
295 // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
296 // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
298 //carre de la longueur de l'arete
299 d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
300 + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
301 aremin = Min( aremin, d );
302 aremax = Max( aremax, d );
304 //debut ajout du 5/10/2006 .............................................
305 //la longueur de l'arete ns1-ns2
307 //longueur arete = Min ( aretmx, aretes incidentes )
308 mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d );
309 mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
310 //fin ajout du 5/10/2006 ...............................................
312 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
313 mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
315 //ajout de l'arete dans la liste
316 fasoar( ns1, ns2, moins1, moins1, n,
317 mosoar, mxsoar, n1soar, mnsoar,
318 mnarst, noar, ierr );
319 //pas de test sur ierr car pas de saturation possible a ce niveau
321 //chainage des aretes frontalieres en position 6 du tableau mnsoar
322 //la nouvelle arete est la suivante de l'arete definie juste avant
323 mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
326 //attention: la derniere arete de la ligne fermee enveloppe
327 // devient en fait la premiere arete de cette ligne
328 // dans le chainage des aretes de la frontiere!
330 if( ierr != 0 ) goto ERREUR;
332 aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
333 aremax = sqrt( aremax ); //longueur maximale d'une arete
335 //debut ajout 9/11/2006 ................................................
336 // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
338 // protection contre une arete max desiree trop grande ou trop petite
339 if( aretmx > aremax*2.05 ) aretmx = aremax;
341 // protection contre une arete max desiree trop petite
342 if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
343 aretmx =(aremin+aremax*2)/3.0;
345 if( aretmx < aremin && aremin > 0 )
348 //sauvegarde pour la fonction areteideale_
349 aretemaxface_ = aretmx;
351 //aire maximale souhaitee des triangles
352 airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral
354 for(i=0; i<=nuds; i++ )
355 mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
356 //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
357 //fin ajout 9/11/2006 .................................................
360 MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax );
361 MESSAGE("Triangulation: arete mx=" << aretmx
362 << " triangle aire mx=" << airemx );
364 //chainage des aretes frontalieres : la derniere arete frontaliere
365 mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
367 //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
368 //reservation du tableau des numeros des 3 aretes de chaque triangle
369 //mnartr( moartr, mxartr )
370 //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
371 // 3Triangles = 2 Aretes internes + Aretes frontalieres
372 // d'ou 3T/2 < AI + AF => T < 3T/2 - Sommets + 1-Trous
373 //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
374 if( mnartr!=NULL ) delete [] mnartr;
375 mxartr = 2 * ( mxsomm + mxtrou );
376 mnartr = new Z[moartr*mxartr];
377 if( mnartr==NULL ) goto ERREUR;
379 //Ajout des points internes
380 ns1 = nudslf[ nblf ];
381 for (i=0; i<nbpti; i++)
383 //les 2 coordonnees du point i de sommet nbs
384 mnpxyd[ns1].x = uvpti[i].x;
385 mnpxyd[ns1].y = uvpti[i].y;
386 mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
387 //le numero i du point interne
392 //nombre de sommets de la frontiere et internes
395 // creation de l'arbre-4 des te (tableau letree)
396 // ajout dans les te des sommets des lignes et des points internes imposes
397 // =======================================================================
398 // premiere estimation de mxtree
401 NEWTREE: //en cas de saturation de l'un des tableaux, on boucle
402 MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
403 if( mntree != NULL ) delete [] mntree;
405 mntree = new Z[motree*(1+mxtree)];
406 if( mntree==NULL ) goto ERREUR;
408 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
409 teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
415 //saturation de letree => sa taille est augmentee et relance
418 MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
424 MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
425 if( ierr != 0 ) goto ERREUR;
426 //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
428 // homogeneisation de l'arbre des te a un saut de taille au plus
429 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
430 // ===========================================================================
431 // reservation de la queue pour parcourir les te de l'arbre
432 if( mnqueu != NULL ) delete [] mnqueu;
434 mnqueu = new Z[mxqueu];
435 if( mnqueu==NULL) goto ERREUR;
437 tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
439 mntree, mxqueu, mnqueu,
444 MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
445 << d << " secondes");
448 //destruction du tableau auxiliaire et de l'arbre
453 MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
461 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
462 // et des points de la frontiere, des points internes imposes interieurs
463 // ==========================================================================
464 tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
465 mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
466 moartr, mxartr, n1artr, mnartr, mnarst,
469 // destruction de la queue et de l'arbre devenus inutiles
470 delete [] mnqueu; mnqueu=NULL;
471 delete [] mntree; mntree=NULL;
476 MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
478 // ierr =0 si pas d'erreur
479 // =1 si le tableau mnsoar est sature
480 // =2 si le tableau mnartr est sature
481 // =3 si aucun des triangles ne contient l'un des points internes
482 // =5 si saturation de la queue de parcours de l'arbre des te
483 if( ierr != 0 ) goto ERREUR;
485 //qualites de la triangulation actuelle
486 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
487 nbt, quamoy, quamin );
489 // boucle sur les aretes internes (non sur une ligne de la frontiere)
490 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
491 // ======================================================================
492 // formation du chainage 6 des aretes internes a echanger eventuellement
493 aisoar( mosoar, mxsoar, mnsoar, na );
494 tedela( mnpxyd, mnarst,
495 mosoar, mxsoar, n1soar, mnsoar, na,
496 moartr, mxartr, n1artr, mnartr, n );
498 MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
501 MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
502 << d << " secondes");
504 //qualites de la triangulation actuelle
505 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
506 nbt, quamoy, quamin );
508 // detection des aretes frontalieres initiales perdues
509 // triangulation frontale pour les restaurer
510 // ===================================================
512 if( mn1arcf != NULL ) delete [] mn1arcf;
513 if( mnarcf != NULL ) delete [] mnarcf;
514 if( mnarcf1 != NULL ) delete [] mnarcf1;
515 if( mnarcf2 != NULL ) delete [] mnarcf2;
516 mn1arcf = new Z[1+mxarcf];
517 if( mn1arcf == NULL ) goto ERREUR;
518 mnarcf = new Z[3*mxarcf];
519 if( mnarcf == NULL ) goto ERREUR;
520 mnarcf1 = new Z[mxarcf];
521 if( mnarcf1 == NULL ) goto ERREUR;
522 mnarcf2 = new Z[mxarcf];
523 if( mnarcf2 == NULL ) goto ERREUR;
525 terefr( nbarpi, mnpxyd,
526 mosoar, mxsoar, n1soar, mnsoar,
527 moartr, mxartr, n1artr, mnartr, mnarst,
528 mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
531 MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr );
534 MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
535 << d << " secondes");
537 if( ierr != 0 ) goto ERREUR;
539 //qualites de la triangulation actuelle
540 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
541 nbt, quamoy, quamin );
543 // fin de la triangulation avec respect des aretes initiales frontalieres
545 // suppression des triangles externes a la surface
546 // ===============================================
547 // recherche du dernier triangle utilise
548 mn = mxartr * moartr;
549 for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
552 if( mnartr[mn] != 0 ) break;
555 if( mntrsu != NULL ) delete [] mntrsu;
556 mntrsu = new Z[ndtri0];
557 if( mntrsu == NULL ) goto ERREUR;
559 if( mnlftr != NULL ) delete [] mnlftr;
560 mnlftr = new Z[nblf];
561 if( mnlftr == NULL ) goto ERREUR;
563 for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
566 tesuex( nblf, mnlftr,
567 ndtri0, nbsomm, mnpxyd, mnslig,
568 mosoar, mxsoar, mnsoar,
569 moartr, mxartr, n1artr, mnartr, mnarst,
572 delete [] mnlftr; mnlftr=NULL;
573 delete [] mntrsu; mntrsu=NULL;
577 MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
578 if( ierr != 0 ) goto ERREUR;
580 //qualites de la triangulation actuelle
581 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
582 nbt, quamoy, quamin );
584 // amelioration de la qualite de la triangulation par
585 // barycentrage des sommets internes a la triangulation
586 // suppression des aretes trop longues ou trop courtes
587 // modification de la topologie des groupes de triangles
588 // mise en delaunay de la triangulation
589 // =====================================================
590 mnarcf3 = new Z[mxarcf];
591 if( mnarcf3 == NULL )
593 MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
596 teamqt( nutysu, aretmx, airemx,
597 mnarst, mosoar, mxsoar, n1soar, mnsoar,
598 moartr, mxartr, n1artr, mnartr,
599 mxarcf, mnarcf2, mnarcf3,
600 mn1arcf, mnarcf, mnarcf1,
601 nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
603 if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
604 if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
605 if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
606 if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
607 if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
611 MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
612 if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
613 if( ierr != 0 ) goto ERREUR;
615 //qualites de la triangulation finale
616 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
617 nbt, quamoy, quamin );
619 // renumerotation des sommets internes: mnarst(i)=numero final du sommet
620 // ===================================
621 for (i=0; i<=nbsomm; i++)
624 for (nt=1; nt<=mxartr; nt++)
626 if( mnartr[nt*moartr-moartr] != 0 )
628 //le numero des 3 sommets du triangle nt
629 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
630 //les 3 sommets du triangle sont actifs
631 mnarst[ nosotr[0] ] = 1;
632 mnarst[ nosotr[1] ] = 1;
633 mnarst[ nosotr[2] ] = 1;
637 for (i=1; i<=nbsomm; i++)
643 // generation du tableau uvst de la surface triangulee
644 // ---------------------------------------------------
645 if( uvst != NULL ) delete [] uvst;
647 if( uvst == NULL ) goto ERREUR;
650 for (i=0; i<nbsomm; i++ )
655 uvst[nbst].x = mnpxyd[i].x;
656 uvst[nbst].y = mnpxyd[i].y;
658 //si le sommet est un point ou appartient a une ligne
659 //ses coordonnees initiales sont restaurees
666 //retour aux coordonnees initiales dans uvslf
668 n = n - 1000000 * l + nudslf[l-1] - 1;
669 uvst[nbst].x = uvslf[n].x;
670 uvst[nbst].y = uvslf[n].y;
674 //point utilisateur n interne impose
675 //retour aux coordonnees initiales dans uvpti
676 uvst[nbst].x = uvpti[n-1].x;
677 uvst[nbst].y = uvpti[n-1].y;
684 // generation du tableau 'nsef' de la surface triangulee
685 // -----------------------------------------------------
686 // boucle sur les triangles occupes (internes et externes)
687 if( nust != NULL ) delete [] nust;
688 nust = new Z[nbsttria*nbt];
689 if( nust == NULL ) goto ERREUR;
691 for (i=1; i<=mxartr; i++)
693 //le triangle i de mnartr
694 if( mnartr[i*moartr-moartr] != 0 )
696 //le triangle i est interne => nosotr numero de ses 3 sommets
697 nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
698 nust[nbt++] = mnarst[ nosotr[0] ];
699 nust[nbt++] = mnarst[ nosotr[1] ];
700 nust[nbt++] = mnarst[ nosotr[2] ];
704 nbt /= nbsttria; //le nombre final de triangles de la surface
705 MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
706 << nbt << " triangles" );
709 MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
711 // destruction des tableaux auxiliaires
712 // ------------------------------------
714 if( mnarst != NULL ) delete [] mnarst;
715 if( mnartr != NULL ) delete [] mnartr;
716 if( mnslig != NULL ) delete [] mnslig;
717 if( mnsoar != NULL ) delete [] mnsoar;
718 if( mnpxyd != NULL ) delete [] mnpxyd;
719 if( mntree != NULL ) delete [] mntree;
720 if( mnqueu != NULL ) delete [] mnqueu;
721 if( mntrsu != NULL ) delete [] mntrsu;
722 if( mnlftr != NULL ) delete [] mnlftr;
723 if( mn1arcf != NULL ) delete [] mn1arcf;
724 if( mnarcf != NULL ) delete [] mnarcf;
725 if( mnarcf1 != NULL ) delete [] mnarcf1;
726 if( mnarcf2 != NULL ) delete [] mnarcf2;
727 if( mnarcf3 != NULL ) delete [] mnarcf3;
731 if( ierr == 51 || ierr == 52 )
733 //saturation des sommets => redepart avec 2 fois plus de sommets
740 MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
741 if( ierr == 0 ) ierr=1;
749 qualitetrte( R3 *mnpxyd,
750 Z & mosoar, Z & mxsoar, Z *mnsoar,
751 Z & moartr, Z & mxartr, Z *mnartr,
752 Z & nbtria, R & quamoy, R & quamin )
753 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
754 // but : calculer la qualite moyenne et minimale de la triangulation
755 // ----- actuelle definie par les tableaux mnsoar et mnartr
758 // mnpxyd : tableau des coordonnees 2d des points
759 // par point : x y distance_souhaitee
760 // mosoar : nombre maximal d'entiers par arete et
761 // indice dans mnsoar de l'arete suivante dans le hachage
762 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
763 // attention: mxsoar>3*mxsomm obligatoire!
764 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
765 // chainage des aretes frontalieres, chainage du hachage des aretes
766 // hachage des aretes = mnsoar(1)+mnsoar(2)*2
767 // avec mxsoar>=3*mxsomm
768 // une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
769 // mnsoar(2,arete vide)=l'arete vide qui precede
770 // mnsoar(3,arete vide)=l'arete vide qui suit
771 // moartr : nombre maximal d'entiers par arete du tableau mnartr
772 // mxartr : nombre maximal de triangles declarables
773 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
774 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
777 // nbtria : nombre de triangles internes au domaine
778 // quamoy : qualite moyenne des triangles actuels
779 // quamin : qualite minimale des triangles actuels
780 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
783 Z nosotr[3], mn, nbtrianeg, nt, ntqmin;
793 for ( nt=1; nt<=mxartr; nt++ )
798 //un triangle occupe de plus
801 //le numero des 3 sommets du triangle nt
802 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
804 //la qualite du triangle ns1 ns2 ns3
805 qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
811 //la qualite minimale
812 if( qualite < quamin )
818 //aire signee du triangle nt
819 d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
822 //un triangle d'aire negative de plus
824 MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
825 << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
826 << " a une aire " << d <<"<=0");
829 //aire des triangles actuels
836 MESSAGE("Qualite moyenne=" << quamoy
837 << " Qualite minimale=" << quamin
838 << " des " << nbtria << " triangles de surface plane totale="
843 //le numero des 3 sommets du triangle ntqmin de qualite minimale
844 nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr, nosotr );
845 MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
846 <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
847 for (int i=0;i<3;i++)
848 MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
849 <<" y="<< mnpxyd[nosotr[i]-1].y);
853 MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );