1 // MEFISTO2: a library to compute 2D triangulation from segmented boundaries
3 // Copyright (C) 2006-2019 CEA/DEN, EDF R&D, OPEN CASCADE
5 // This library is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU Lesser General Public
7 // License as published by the Free Software Foundation; either
8 // version 2.1 of the License, or (at your option) any later version.
10 // This library is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // Lesser General Public License for more details.
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with this library; if not, write to the Free Software
17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 // File : aptrte.cxx le C++ de l'appel du trianguleur plan
23 // Author : Alain PERRONNET
24 // Date : 13 novembre 2006
28 #include "utilities.h"
43 areteideale()//( R3 xyz, R3 direction )
48 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
49 //dans la direction donnee
50 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
53 static double cpunew, cpuold=0;
62 tempscpu_( double & tempsec )
63 //Retourne le temps CPU utilise en secondes
65 tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
66 //MESSAGE( "temps cpu=" << tempsec );
77 deltacpu_( R & dtcpu )
78 //Retourne le temps CPU utilise en secondes depuis le precedent appel
81 dtcpu = R( cpunew - cpuold );
83 //MESSAGE( "delta temps cpu=" << dtcpu );
88 void aptrte( Z nutysu, R aretmx,
89 Z nblf, Z * nudslf, R2 * uvslf,
91 Z & nbst, R2 * & uvst,
94 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
95 // but : appel de la triangulation par un arbre-4 recouvrant
96 // ----- de triangles equilateraux
97 // le contour du domaine plan est defini par des lignes fermees
98 // la premiere ligne etant l'enveloppe de toutes les autres
99 // la fonction areteideale(s,d) donne la taille d'arete
100 // au point s dans la direction (actuellement inactive) d
101 // des lors toute arete issue d'un sommet s devrait avoir une longueur
102 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
105 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
106 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
110 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
111 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
112 // 1 il existe une fonction areteideale_(s,d)
113 // dont seules les 2 premieres composantes de uv sont actives
114 // ... autres options a definir ...
115 // aretmx : longueur maximale des aretes de la future triangulation
116 // nblf : nombre de lignes fermees de la surface
117 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
118 // nudslf(0)=0 pour permettre la difference sans test
119 // Attention le dernier sommet de chaque ligne est raccorde au premier
120 // tous les sommets et les points internes ont des coordonnees
121 // UV differentes <=> Pas de point double!
122 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
123 // nbpti : nombre de points internes futurs sommets de la triangulation
124 // uvpti : uv des points internes futurs sommets de la triangulation
128 // nbst : nombre de sommets de la triangulation finale
129 // uvst : coordonnees uv des nbst sommets de la triangulation
130 // nbt : nombre de triangles de la triangulation finale
131 // nust : 4 numeros dans uvst des sommets des nbt triangles
132 // s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
133 // ierr : 0 si pas d'erreur
135 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
136 // auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006
137 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139 Z nbsttria=4; //Attention: 4 sommets stockes par triangle
140 //no st1, st2, st3, 0 (non quadrangle)
143 // R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
144 Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
145 Z mxtrou = Max( 1024, nblf ); //nombre maximal de trous dans la surface
148 Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
149 Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
150 Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
151 Z *mnqueu=NULL, mxqueu;
153 Z *mnarcf=NULL, mxarcf;
162 R3 comxmi[2]; //coordonnees UV Min et Maximales
163 R aremin, aremax; //longueur minimale et maximale des aretes
164 R airemx; //aire maximale souhaitee d'un triangle
168 Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
169 Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
173 // initialisation du temps cpu
177 // quelques reservations de tableaux pour faire les calculs
178 // ========================================================
179 // declaration du tableau des coordonnees des sommets de la frontiere
180 // puis des sommets internes ajoutes
181 // majoration empirique du nombre de sommets de la triangulation
183 mxsomm = Max( 20000, 64*nbpti+i*i );
184 // MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
185 // MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx
186 // << " mxsomm=" << mxsomm );
187 // MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
190 //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
191 if( mnpxyd!=NULL ) delete [] mnpxyd;
192 mnpxyd = new R3[mxsomm];
193 if( mnpxyd==NULL ) goto ERREUR;
195 // le tableau mnsoar des aretes des triangles
196 // 1: sommet 1 dans pxyd,
197 // 2: sommet 2 dans pxyd,
198 // 3: numero de 1 a nblf de la ligne qui supporte l'arete
199 // 4: numero dans mnartr du triangle 1 partageant cette arete,
200 // 5: numero dans mnartr du triangle 2 partageant cette arete,
201 // 6: chainage des aretes frontalieres ou internes ou
202 // des aretes simples des etoiles de triangles,
203 // 7: chainage du hachage des aretes
204 // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
205 // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
206 // h(ns1,ns2) = min( ns1, ns2 )
207 if( mnsoar!=NULL ) delete [] mnsoar;
208 mxsoar = 3 * ( mxsomm + mxtrou );
209 mnsoar = new Z[mosoar*mxsoar];
210 if( mnsoar==NULL ) goto ERREUR;
211 //initialiser le tableau mnsoar pour le hachage des aretes
212 insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
214 // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
215 if( mnarst!=NULL ) delete [] mnarst;
216 mnarst = new Z[1+mxsomm];
217 if( mnarst==NULL ) goto ERREUR;
221 // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
222 // ou no du point si interne forc'e par l'utilisateur
223 // ou 0 si interne cree par le module
224 if( mnslig!=NULL ) delete [] mnslig;
225 mnslig = new Z[mxsomm];
226 if( mnslig==NULL ) goto ERREUR;
227 azeroi( mxsomm, mnslig );
229 // initialisation des aretes frontalieres de la triangulation future
230 // renumerotation des sommets des aretes des lignes pour la triangulation
231 // mise a l'echelle des coordonnees des sommets pour obtenir une
232 // meilleure precision lors des calculs + quelques verifications
233 // boucle sur les lignes fermees qui forment la frontiere
234 // ======================================================================
239 for (n=1; n<=nblf; n++)
241 //l'initialisation de la premiere arete de la ligne n dans la triangulation
242 //-------------------------------------------------------------------------
243 //le sommet ns0 est le numero de l'origine de la ligne
245 mnpxyd[ns0].x = uvslf[ns0].x;
246 mnpxyd[ns0].y = uvslf[ns0].y;
247 mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
248 // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
249 // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
251 //carre de la longueur de l'arete 1 de la ligne fermee n
252 d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
253 + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
254 aremin = Min( aremin, d );
255 aremax = Max( aremax, d );
257 //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
258 //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
259 //le numero des 2 sommets ns1 ns2 de la 1-ere arete
260 //Attention: les numeros ns debutent a 1 (ils ont >0)
261 // les tableaux c++ demarrent a zero!
262 // les tableaux fortran demarrent ou l'on veut!
267 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
268 mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
269 fasoar( ns1, ns2, moins1, moins1, n,
270 mosoar, mxsoar, n1soar, mnsoar, mnarst,
272 //pas de test sur ierr car pas de saturation possible a ce niveau
274 //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
277 //la nouvelle arete est la suivante de l'arete definie juste avant
279 mnsoar[mosoar * noar - mosoar + 5] = noar0;
281 //l'initialisation des aretes suivantes de la ligne dans la triangulation
282 //-----------------------------------------------------------------------
283 nbarli = nudslf[n] - nudslf[n-1]; //nombre d'aretes=sommets de la ligne n
284 for (i=2; i<=nbarli; i++)
286 ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
288 //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
291 //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
294 //l'arete precedente est dotee de sa suivante:celle cree ensuite
295 //les 2 coordonnees du sommet ns2 de la ligne
297 //debut ajout 5/10/2006 ................................................
298 nuds = Max( nuds, ns ); //le numero du dernier sommet traite
299 //fin ajout 5/10/2006 ................................................
300 mnpxyd[ns].x = uvslf[ns].x;
301 mnpxyd[ns].y = uvslf[ns].y;
302 mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
303 // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
304 // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
306 //carre de la longueur de l'arete
307 d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
308 + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
309 aremin = Min( aremin, d );
310 aremax = Max( aremax, d );
312 //debut ajout du 5/10/2006 .............................................
313 //la longueur de l'arete ns1-ns2
315 //longueur arete = Min ( aretmx, aretes incidentes )
316 mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d );
317 mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
318 //fin ajout du 5/10/2006 ...............................................
320 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
321 mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
323 //ajout de l'arete dans la liste
324 fasoar( ns1, ns2, moins1, moins1, n,
325 mosoar, mxsoar, n1soar, mnsoar,
326 mnarst, noar, ierr );
327 //pas de test sur ierr car pas de saturation possible a ce niveau
329 //chainage des aretes frontalieres en position 6 du tableau mnsoar
330 //la nouvelle arete est la suivante de l'arete definie juste avant
331 mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
334 //attention: la derniere arete de la ligne fermee enveloppe
335 // devient en fait la premiere arete de cette ligne
336 // dans le chainage des aretes de la frontiere!
338 if( ierr != 0 ) goto ERREUR;
340 aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
341 aremax = sqrt( aremax ); //longueur maximale d'une arete
343 //debut ajout 9/11/2006 ................................................
344 // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
346 // protection contre une arete max desiree trop grande ou trop petite
347 if( aretmx > aremax*2.05 ) aretmx = aremax;
349 // protection contre une arete max desiree trop petite
350 if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
351 aretmx =(aremin+aremax*2)/3.0;
353 if( aretmx < aremin && aremin > 0 )
356 //sauvegarde pour la fonction areteideale_
357 aretemaxface_ = aretmx;
359 //aire maximale souhaitee des triangles
360 airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral
362 for(i=0; i<=nuds; i++ )
363 mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
364 //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
365 //fin ajout 9/11/2006 .................................................
368 // MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax );
369 // MESSAGE("Triangulation: arete mx=" << aretmx
370 // << " triangle aire mx=" << airemx );
372 //chainage des aretes frontalieres : la derniere arete frontaliere
373 mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
375 //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
376 //reservation du tableau des numeros des 3 aretes de chaque triangle
377 //mnartr( moartr, mxartr )
378 //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
379 // 3Triangles = 2 Aretes internes + Aretes frontalieres
380 // d'ou 3T/2 < AI + AF => T < 3T/2 - Sommets + 1-Trous
381 //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
382 if( mnartr!=NULL ) delete [] mnartr;
383 mxartr = 2 * ( mxsomm + mxtrou );
384 mnartr = new Z[moartr*mxartr];
385 if( mnartr==NULL ) goto ERREUR;
387 //Ajout des points internes
388 ns1 = nudslf[ nblf ];
389 for (i=0; i<nbpti; i++)
391 //les 2 coordonnees du point i de sommet nbs
392 mnpxyd[ns1].x = uvpti[i].x;
393 mnpxyd[ns1].y = uvpti[i].y;
394 mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
395 //le numero i du point interne
400 //nombre de sommets de la frontiere et internes
403 // creation de l'arbre-4 des te (tableau letree)
404 // ajout dans les te des sommets des lignes et des points internes imposes
405 // =======================================================================
406 // premiere estimation de mxtree
409 NEWTREE: //en cas de saturation de l'un des tableaux, on boucle
410 //MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
411 if( mntree != NULL ) delete [] mntree;
413 mntree = new Z[motree*(1+mxtree)];
414 if( mntree==NULL ) goto ERREUR;
416 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
417 comxmi[0].x = comxmi[1].x = uvslf[0].x;
418 comxmi[0].y = comxmi[1].y = uvslf[0].y;
419 teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
425 //saturation de letree => sa taille est augmentee et relance
428 //MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
434 //MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
435 if( ierr != 0 ) goto ERREUR;
436 //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
438 // homogeneisation de l'arbre des te a un saut de taille au plus
439 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
440 // ===========================================================================
441 // reservation de la queue pour parcourir les te de l'arbre
442 if( mnqueu != NULL ) delete [] mnqueu;
444 mnqueu = new Z[mxqueu];
445 if( mnqueu==NULL) goto ERREUR;
447 tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
449 mntree, mxqueu, mnqueu,
454 //MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
455 // << d << " secondes");
458 //destruction du tableau auxiliaire et de l'arbre
463 //MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
471 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
472 // et des points de la frontiere, des points internes imposes interieurs
473 // ==========================================================================
474 tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
475 mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
476 moartr, mxartr, n1artr, mnartr, mnarst,
479 // destruction de la queue et de l'arbre devenus inutiles
480 delete [] mnqueu; mnqueu=NULL;
481 delete [] mntree; mntree=NULL;
486 //MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
488 // ierr =0 si pas d'erreur
489 // =1 si le tableau mnsoar est sature
490 // =2 si le tableau mnartr est sature
491 // =3 si aucun des triangles ne contient l'un des points internes
492 // =5 si saturation de la queue de parcours de l'arbre des te
493 if( ierr != 0 ) goto ERREUR;
495 //qualites de la triangulation actuelle
496 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
497 nbt, quamoy, quamin );
499 // boucle sur les aretes internes (non sur une ligne de la frontiere)
500 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
501 // ======================================================================
502 // formation du chainage 6 des aretes internes a echanger eventuellement
503 aisoar( mosoar, mxsoar, mnsoar, na );
504 tedela( mnpxyd, mnarst,
505 mosoar, mxsoar, n1soar, mnsoar, na,
506 moartr, mxartr, n1artr, mnartr, n );
508 //MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
511 // MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
512 // << d << " secondes");
514 //qualites de la triangulation actuelle
515 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
516 nbt, quamoy, quamin );
518 // detection des aretes frontalieres initiales perdues
519 // triangulation frontale pour les restaurer
520 // ===================================================
522 if( mn1arcf != NULL ) delete [] mn1arcf;
523 if( mnarcf != NULL ) delete [] mnarcf;
524 if( mnarcf1 != NULL ) delete [] mnarcf1;
525 if( mnarcf2 != NULL ) delete [] mnarcf2;
526 mn1arcf = new Z[1+mxarcf];
527 if( mn1arcf == NULL ) goto ERREUR;
528 mnarcf = new Z[3*mxarcf];
529 if( mnarcf == NULL ) goto ERREUR;
530 mnarcf1 = new Z[mxarcf];
531 if( mnarcf1 == NULL ) goto ERREUR;
532 mnarcf2 = new Z[mxarcf];
533 if( mnarcf2 == NULL ) goto ERREUR;
535 terefr( nbarpi, mnpxyd,
536 mosoar, mxsoar, n1soar, mnsoar,
537 moartr, mxartr, n1artr, mnartr, mnarst,
538 mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
541 //MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr );
544 //MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
545 // << d << " secondes");
547 if( ierr != 0 ) goto ERREUR;
549 //qualites de la triangulation actuelle
550 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
551 nbt, quamoy, quamin );
553 // fin de la triangulation avec respect des aretes initiales frontalieres
555 // suppression des triangles externes a la surface
556 // ===============================================
557 // recherche du dernier triangle utilise
558 mn = mxartr * moartr;
559 for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
562 if( mnartr[mn] != 0 ) break;
565 if( mntrsu != NULL ) delete [] mntrsu;
566 mntrsu = new Z[ndtri0];
567 if( mntrsu == NULL ) goto ERREUR;
569 if( mnlftr != NULL ) delete [] mnlftr;
570 mnlftr = new Z[nblf];
571 if( mnlftr == NULL ) goto ERREUR;
573 for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
576 tesuex( nblf, mnlftr,
577 ndtri0, nbsomm, mnpxyd, mnslig,
578 mosoar, mxsoar, mnsoar,
579 moartr, mxartr, n1artr, mnartr, mnarst,
582 delete [] mnlftr; mnlftr=NULL;
583 delete [] mntrsu; mntrsu=NULL;
587 //MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
588 if( ierr != 0 ) goto ERREUR;
590 //qualites de la triangulation actuelle
591 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
592 nbt, quamoy, quamin );
594 // amelioration de la qualite de la triangulation par
595 // barycentrage des sommets internes a la triangulation
596 // suppression des aretes trop longues ou trop courtes
597 // modification de la topologie des groupes de triangles
598 // mise en delaunay de la triangulation
599 // =====================================================
600 mnarcf3 = new Z[mxarcf];
601 if( mnarcf3 == NULL )
603 MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
606 teamqt( nutysu, aretmx, airemx,
607 mnarst, mosoar, mxsoar, n1soar, mnsoar,
608 moartr, mxartr, n1artr, mnartr,
609 mxarcf, mnarcf2, mnarcf3,
610 mn1arcf, mnarcf, mnarcf1,
611 nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
613 if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
614 if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
615 if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
616 if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
617 if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
621 //MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
622 if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
623 if( ierr != 0 ) goto ERREUR;
625 //qualites de la triangulation finale
626 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
627 nbt, quamoy, quamin );
629 // renumerotation des sommets internes: mnarst(i)=numero final du sommet
630 // ===================================
631 for (i=0; i<=nbsomm; i++)
634 for (nt=1; nt<=mxartr; nt++)
636 if( mnartr[nt*moartr-moartr] != 0 )
638 //le numero des 3 sommets du triangle nt
639 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
640 //les 3 sommets du triangle sont actifs
641 mnarst[ nosotr[0] ] = 1;
642 mnarst[ nosotr[1] ] = 1;
643 mnarst[ nosotr[2] ] = 1;
647 for (i=1; i<=nbsomm; i++)
653 // generation du tableau uvst de la surface triangulee
654 // ---------------------------------------------------
655 if( uvst != NULL ) delete [] uvst;
657 if( uvst == NULL ) goto ERREUR;
660 for (i=0; i<nbsomm; i++ )
665 uvst[nbst].x = mnpxyd[i].x;
666 uvst[nbst].y = mnpxyd[i].y;
668 //si le sommet est un point ou appartient a une ligne
669 //ses coordonnees initiales sont restaurees
676 //retour aux coordonnees initiales dans uvslf
678 n = n - 1000000 * l + nudslf[l-1] - 1;
679 uvst[nbst].x = uvslf[n].x;
680 uvst[nbst].y = uvslf[n].y;
684 //point utilisateur n interne impose
685 //retour aux coordonnees initiales dans uvpti
686 uvst[nbst].x = uvpti[n-1].x;
687 uvst[nbst].y = uvpti[n-1].y;
694 // generation du tableau 'nsef' de la surface triangulee
695 // -----------------------------------------------------
696 // boucle sur les triangles occupes (internes et externes)
697 if( nust != NULL ) delete [] nust;
698 nust = new Z[nbsttria*nbt];
699 if( nust == NULL ) goto ERREUR;
701 for (i=1; i<=mxartr; i++)
703 //le triangle i de mnartr
704 if( mnartr[i*moartr-moartr] != 0 )
706 //le triangle i est interne => nosotr numero de ses 3 sommets
707 nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
708 nust[nbt++] = mnarst[ nosotr[0] ];
709 nust[nbt++] = mnarst[ nosotr[1] ];
710 nust[nbt++] = mnarst[ nosotr[2] ];
714 nbt /= nbsttria; //le nombre final de triangles de la surface
715 // MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
716 // << nbt << " triangles" );
719 // MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
721 // destruction des tableaux auxiliaires
722 // ------------------------------------
724 if( mnarst != NULL ) delete [] mnarst;
725 if( mnartr != NULL ) delete [] mnartr;
726 if( mnslig != NULL ) delete [] mnslig;
727 if( mnsoar != NULL ) delete [] mnsoar;
728 if( mnpxyd != NULL ) delete [] mnpxyd;
729 if( mntree != NULL ) delete [] mntree;
730 if( mnqueu != NULL ) delete [] mnqueu;
731 if( mntrsu != NULL ) delete [] mntrsu;
732 if( mnlftr != NULL ) delete [] mnlftr;
733 if( mn1arcf != NULL ) delete [] mn1arcf;
734 if( mnarcf != NULL ) delete [] mnarcf;
735 if( mnarcf1 != NULL ) delete [] mnarcf1;
736 if( mnarcf2 != NULL ) delete [] mnarcf2;
737 if( mnarcf3 != NULL ) delete [] mnarcf3;
741 if( ierr == 51 || ierr == 52 )
743 //saturation des sommets => redepart avec 2 fois plus de sommets
750 MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
751 if( ierr == 0 ) ierr=1;
762 qualitetrte( R3 *mnpxyd,
763 Z & mosoar, Z & mxsoar, Z *mnsoar,
764 Z & moartr, Z & mxartr, Z *mnartr,
765 Z & nbtria, R & quamoy, R & quamin )
766 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
767 // but : calculer la qualite moyenne et minimale de la triangulation
768 // ----- actuelle definie par les tableaux mnsoar et mnartr
771 // mnpxyd : tableau des coordonnees 2d des points
772 // par point : x y distance_souhaitee
773 // mosoar : nombre maximal d'entiers par arete et
774 // indice dans mnsoar de l'arete suivante dans le hachage
775 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
776 // attention: mxsoar>3*mxsomm obligatoire!
777 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
778 // chainage des aretes frontalieres, chainage du hachage des aretes
779 // hachage des aretes = mnsoar(1)+mnsoar(2)*2
780 // avec mxsoar>=3*mxsomm
781 // une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
782 // mnsoar(2,arete vide)=l'arete vide qui precede
783 // mnsoar(3,arete vide)=l'arete vide qui suit
784 // moartr : nombre maximal d'entiers par arete du tableau mnartr
785 // mxartr : nombre maximal de triangles declarables
786 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
787 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
790 // nbtria : nombre de triangles internes au domaine
791 // quamoy : qualite moyenne des triangles actuels
792 // quamin : qualite minimale des triangles actuels
793 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
796 Z nosotr[3], mn, nbtrianeg, nt, ntqmin;
806 for ( nt=1; nt<=mxartr; nt++ )
811 //un triangle occupe de plus
814 //le numero des 3 sommets du triangle nt
815 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
817 //la qualite du triangle ns1 ns2 ns3
818 qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
824 //la qualite minimale
825 if( qualite < quamin )
831 //aire signee du triangle nt
832 d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
835 //un triangle d'aire negative de plus
837 MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
838 << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
839 << " a une aire " << d <<"<=0");
842 //aire des triangles actuels
849 // MESSAGE("Qualite moyenne=" << quamoy
850 // << " Qualite minimale=" << quamin
851 // << " des " << nbtria << " triangles de surface plane totale="
856 //le numero des 3 sommets du triangle ntqmin de qualite minimale
857 nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr, nosotr );
858 // MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
859 // <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
860 // for (int i=0;i<3;i++)
861 // MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
862 // <<" y="<< mnpxyd[nosotr[i]-1].y);
866 MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );