1 // SMESH MEFISTO2 : algorithm for meshing
3 // Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris
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.
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.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
24 // Author : Alain PERRONNET
26 // Date : 13 novembre 2006
31 #include <climits> // limites min max int long real ...
33 #include <unistd.h> // gethostname, ...
37 #include <iostream> // pour cout cin ...
38 #include <iomanip> // pour le format des io setw, stx, setfill, ...
40 #include <string.h> // pour les fonctions sur les chaines de caracteres
43 #include <math.h> // pour les fonctions mathematiques
46 #include <sys/types.h>
52 #if defined MEFISTO2D_EXPORTS
53 #define MEFISTO2D_EXPORT __declspec( dllexport )
55 #define MEFISTO2D_EXPORT __declspec( dllimport )
58 #define MEFISTO2D_EXPORT
63 void aptrte( Z nutysu, R aretmx,
64 Z nblf, Z *nudslf, R2 *uvslf,
66 Z & nbst, R2 * & uvst, Z & nbt, Z * & nust,
68 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
69 // but : appel de la triangulation par un arbre-4 recouvrant
70 // ----- de triangles equilateraux
71 // le contour du domaine plan est defini par des lignes fermees
72 // la premiere ligne etant l'enveloppe de toutes les autres
73 // la fonction areteideale_(s,d) donne la taille d'arete
74 // au point s dans la direction d (direction inactive pour l'instant)
75 // des lors toute arete issue d'un sommet s devrait avoir une longueur
76 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
79 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
80 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
84 // nutysu : numero de traitement de areteideale_() selon le type de surface
85 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
86 // 1 il existe une fonction areteideale_(s,d)
87 // dont seules les 2 premieres composantes de uv sont actives
88 // ... autres options a definir ...
89 // aretmx : longueur maximale des aretes de la future triangulation
90 // nblf : nombre de lignes fermees de la surface
91 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
92 // nudslf(0)=0 pour permettre la difference sans test
93 // Attention le dernier sommet de chaque ligne est raccorde au premier
94 // tous les sommets et les points internes ont des coordonnees
95 // UV differentes <=> Pas de point double!
96 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
97 // nbpti : nombre de points internes futurs sommets de la triangulation
98 // uvpti : uv des points internes futurs sommets de la triangulation
102 // nbst : nombre de sommets de la triangulation finale
103 // uvst : coordonnees uv des nbst sommets de la triangulation
104 // nbt : nombre de triangles de la triangulation finale
105 // nust : 3 numeros dans uvst des sommets des nbt triangles
106 // ierr : 0 si pas d'erreur
108 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
109 // auteur : Alain Perronnet Analyse Numerique Paris UPMC decembre 2001
110 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113 #define tempscpu TEMPSCPU
114 #define deltacpu DELTACPU
115 #define insoar INSOAR
116 #define azeroi AZEROI
117 #define fasoar FASOAR
118 #define teajte TEAJTE
119 #define tehote TEHOTE
120 #define tetrte TETRTE
121 #define aisoar AISOAR
122 #define tedela TEDELA
123 #define terefr TEREFR
124 #define tesuex TESUEX
125 #define teamqt TEAMQT
126 #define nusotr NUSOTR
127 #define qutr2d QUTR2D
128 #define surtd2 SURTD2
129 #define qualitetrte QUALITETRTE
131 #define areteideale ARETEIDEALE
134 #define tempscpu tempscpu_
135 #define deltacpu deltacpu_
136 #define insoar insoar_
137 #define azeroi azeroi_
138 #define fasoar fasoar_
139 #define teajte teajte_
140 #define tehote tehote_
141 #define tetrte tetrte_
142 #define aisoar aisoar_
143 #define tedela tedela_
144 #define terefr terefr_
145 #define tesuex tesuex_
146 #define teamqt teamqt_
147 #define nusotr nusotr_
148 #define qutr2d qutr2d_
149 #define surtd2 surtd2_
150 #define qualitetrte qualitetrte_
152 #define areteideale areteideale_
164 qualitetrte( R3 *mnpxyd,
165 Z & mosoar, Z & mxsoar, Z *mnsoar,
166 Z & moartr, Z & mxartr, Z *mnartr,
167 Z & nbtria, R & quamoy, R & quamin ); }
168 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
169 // but : calculer la qualite moyenne et minimale de la triangulation
170 // ----- actuelle definie par les tableaux nosoar et noartr
173 // mnpxyd : tableau des coordonnees 2d des points
174 // par point : x y distance_souhaitee
175 // mosoar : nombre maximal d'entiers par arete et
176 // indice dans nosoar de l'arete suivante dans le hachage
177 // mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
178 // attention: mxsoar>3*mxsomm obligatoire!
179 // nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
180 // chainage des aretes frontalieres, chainage du hachage des aretes
181 // hachage des aretes = nosoar(1)+nosoar(2)*2
182 // avec mxsoar>=3*mxsomm
183 // une arete i de nosoar est vide <=> nosoar(1,i)=0 et
184 // nosoar(2,arete vide)=l'arete vide qui precede
185 // nosoar(3,arete vide)=l'arete vide qui suit
186 // moartr : nombre maximal d'entiers par arete du tableau noartr
187 // mxartr : nombre maximal de triangles declarables
188 // noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
189 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
192 // nbtria : nombre de triangles internes au domaine
193 // quamoy : qualite moyenne des triangles actuels
194 // quamin : qualite minimale des triangles actuels
195 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
204 tempscpu( double & tempsec );
207 //Retourne le temps CPU utilise en secondes
216 deltacpu( R & dtcpu );
219 //Retourne le temps CPU utilise en secondes depuis le precedent appel
221 //initialiser le tableau mnsoar pour le hachage des aretes
229 insoar( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );
232 //mettre a zero les nb entiers de tab
240 azeroi( Z & nb, Z * tab );
250 fasoar( Z & ns1, Z & ns2, Z & nt1, Z & nt2, Z & nolign,
251 Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z * mnarst,
252 Z & noar, Z & ierr );
254 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
255 // but : former l'arete de sommet ns1-ns2 dans le hachage du tableau
256 // ----- nosoar des aretes de la triangulation
259 // ns1 ns2: numero pxyd des 2 sommets de l'arete
260 // nt1 : numero du triangle auquel appartient l'arete
261 // nt1=-1 si numero inconnu
262 // nt2 : numero de l'eventuel second triangle de l'arete si connu
263 // nt2=-1 si numero inconnu
264 // nolign : numero de la ligne fermee de l'arete
265 // =0 si l'arete n'est une arete de ligne
266 // ce numero est ajoute seulement si l'arete est creee
267 // mosoar : nombre maximal d'entiers par arete du tableau nosoar
268 // mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
271 // n1soar : numero de la premiere arete vide dans le tableau nosoar
272 // une arete i de nosoar est vide <=> nosoar(1,i)=0
273 // chainage des aretes vides amont et aval
274 // l'arete vide qui precede=nosoar(4,i)
275 // l'arete vide qui suit =nosoar(5,i)
276 // nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
277 // chainage momentan'e d'aretes, chainage du hachage des aretes
278 // hachage des aretes = min( nosoar(1), nosoar(2) )
279 // noarst : noarst(np) numero d'une arete du sommet np
281 // ierr : si < 0 en entree pas d'affichage en cas d'erreur du type
282 // "arete appartenant a plus de 2 triangles et a creer!"
283 // si >=0 en entree affichage de ce type d'erreur
286 // noar : >0 numero de l'arete retrouvee ou ajoutee
287 // ierr : =0 si pas d'erreur
288 // =1 si le tableau nosoar est sature
289 // =2 si arete a creer et appartenant a 2 triangles distincts
290 // des triangles nt1 et nt2
291 // =3 si arete appartenant a 2 triangles distincts
292 // differents des triangles nt1 et nt2
293 // =4 si arete appartenant a 2 triangles distincts
294 // dont le second n'est pas le triangle nt2
295 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
297 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
305 teajte( Z & mxsomm, Z & nbsomm, R3 * mnpxyd, R3 * comxmi,
306 R & aretmx, Z & mxtree, Z * letree,
317 tehote( Z & nutysu, Z & nbarpi, Z & mxsomm, Z & nbsomm, R3 * mnpxyd,
318 R3 * comxmi, R & aretmx,
319 Z * letree, Z & mxqueu, Z * mnqueu,
322 // homogeneisation de l'arbre des te a un saut de taille au plus
323 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
332 tetrte( R3 * comxmi, R & aretmx, Z & nbarpi, Z & mxsomm, R3 * mnpxyd,
333 Z & mxqueu, Z * mnqueu, Z * mntree,
334 Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
335 Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
338 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
339 // et des points de la frontiere, des points internes imposes interieurs
348 aisoar( Z & mosoar, Z & mxsoar, Z * mnsoar, Z & na );
350 // formation du chainage 6 des aretes internes a echanger eventuellement
359 tedela( R3 * mnpxyd, Z * mnarst,
360 Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & na,
361 Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z & n );
363 // boucle sur les aretes internes (non sur une ligne de la frontiere)
364 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
373 terefr( Z & nbarpi, R3 * mnpxyd,
374 Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
375 Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
376 Z & mxarcf, Z * mnarc1, Z * mnarc2,
377 Z * mnarc3, Z * mnarc4,
380 // detection des aretes frontalieres initiales perdues
381 // triangulation frontale pour les restaurer
390 tesuex( Z & nblf, Z * nulftr,
391 Z & ndtri0, Z & nbsomm, R3 * mnpxyd, Z * mnslig,
392 Z & mosoar, Z & mxsoar, Z * mnsoar,
393 Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
394 Z & nbtria, Z * mntrsu, Z & ierr );
396 // suppression des triangles externes a la surface
405 teamqt( Z & nutysu, R & aretmx, R & airemx,
406 Z * mnarst, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
407 Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr,
408 Z & mxarcf, Z * mntrcf, Z * mnstbo,
409 Z * n1arcf, Z * mnarcf, Z * mnarc1,
410 Z & nbarpi, Z & nbsomm, Z & mxsomm,
411 R3 * mnpxyd, Z * mnslig,
414 // amelioration de la qualite de la triangulation par
415 // barycentrage des sommets internes a la triangulation
416 // suppression des aretes trop longues ou trop courtes
417 // modification de la topologie des groupes de triangles
418 // mise en delaunay de la triangulation
427 nusotr( Z & nt, Z & mosoar, Z * mnsoar, Z & moartr, Z * mnartr,Z * nosotr );
429 //retrouver les numero des 3 sommets du triangle nt
438 qutr2d( R3 & p1, R3 & p2, R3 & p3, R & qualite );
440 //calculer la qualite d'un triangle de R2 de sommets p1, p2, p3
449 surtd2( R3 & p1, R3 & p2, R3 & p3 );
451 //calcul de la surface d'un triangle defini par 3 points de r**2