From 1db8a57bb11dbde1bb12e338257691936f4dfcb3 Mon Sep 17 00:00:00 2001 From: abd Date: Wed, 26 Apr 2006 10:41:24 +0000 Subject: [PATCH] Compile new version of MEFISTO with DF6.0 --- src/MEFISTO2/aptrte.cxx | 254 +++++++++++++++------------------------- src/MEFISTO2/aptrte.h | 231 +++++++++++++++++------------------- src/MEFISTO2/trte.f | 227 +++++++++++++++++++++++++---------- 3 files changed, 362 insertions(+), 350 deletions(-) diff --git a/src/MEFISTO2/aptrte.cxx b/src/MEFISTO2/aptrte.cxx index 8562631a6..d2b4b3693 100755 --- a/src/MEFISTO2/aptrte.cxx +++ b/src/MEFISTO2/aptrte.cxx @@ -1,6 +1,6 @@ -// MEFISTO : library to compute 2D triangulation from segmented boundaries +// MEFISTO2: a library to compute 2D triangulation from segmented boundaries // -// Copyright (C) 2003 Laboratoire J.-L. Lions UPMC Paris +// Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -18,23 +18,26 @@ // // See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr // -// -// File : aptrte.cxx +// File : aptrte.cxx le C++ de l'appel du trianguleur plan // Module : SMESH -// Author: Alain PERRONNET +// Author : Alain PERRONNET +// Date : 16 mars 2006 #include "Rn.h" #include "aptrte.h" -#ifndef WIN32 #include "utilities.h" using namespace std; -#endif extern "C" { R aretemaxface_; - R areteideale_( R3 xyz, R3 direction ) + MEFISTO2D_EXPORT + R + #ifdef WIN32 + __stdcall + #endif + areteideale()//( R3 xyz, R3 direction ) { return aretemaxface_; } @@ -50,7 +53,7 @@ void tempscpu_( double & tempsec ) //Retourne le temps CPU utilise en secondes { tempsec = ( (double) clock() ) / CLOCKS_PER_SEC; - // MESSAGEE( "temps cpu=" << tempsec ); + //MESSAGE( "temps cpu=" << tempsec ); } @@ -60,15 +63,15 @@ void deltacpu_( R & dtcpu ) tempscpu_( cpunew ); dtcpu = R( cpunew - cpuold ); cpuold = cpunew; - // MESSAGEE( "delta temps cpu=" << dtcpu ); + //MESSAGE( "delta temps cpu=" << dtcpu ); return; } - -void aptrte( Z nutysu, R aretmx, - Z nblf, Z * nudslf, R2 * uvslf, - Z nbpti, R2 *uvpti, - Z & nbst, R2 * & uvst, Z & nbt, Z * & nust, +void aptrte( Z nutysu, R aretmx, + Z nblf, Z * nudslf, R2 * uvslf, + Z nbpti, R2 * uvpti, + Z & nbst, R2 * & uvst, + Z & nbt, Z * & nust, Z & ierr ) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // but : appel de la triangulation par un arbre-4 recouvrant @@ -112,9 +115,12 @@ void aptrte( Z nutysu, R aretmx, // ierr : 0 si pas d'erreur // > 0 sinon //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -// auteur : Alain Perronnet Analyse Numerique Paris UPMC decembre 2001 +// auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ { + Z nbsttria=4; //Attention: 4 sommets stockes par triangle + //no st1, st2, st3, 0 (non quadrangle) + R d, tcpu=0; R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici! Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees @@ -123,7 +129,7 @@ void aptrte( Z nutysu, R aretmx, R3 *mnpxyd=NULL; Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles - Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE + Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE Z *mnqueu=NULL, mxqueu; Z *mn1arcf=NULL; Z *mnarcf=NULL, mxarcf; @@ -131,7 +137,6 @@ void aptrte( Z nutysu, R aretmx, Z *mnarcf2=NULL; Z *mnarcf3=NULL; Z *mntrsu=NULL; - Z *mndalf=NULL; Z *mnslig=NULL; Z *mnarst=NULL; Z *mnlftr=NULL; @@ -144,6 +149,7 @@ void aptrte( Z nutysu, R aretmx, Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt; Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn; Z moins1=-1; + R dist; aretemaxface_ = aretmx; @@ -153,19 +159,14 @@ void aptrte( Z nutysu, R aretmx, // quelques reservations de tableaux pour faire les calculs // ======================================================== - // le tableau pointeur sur la premiere arete de chaque ligne fermee - if( mndalf!=NULL ) delete [] mndalf; - mndalf = new Z[1+nblf]; - if( mndalf==NULL ) goto ERREUR; - mndalf[0]=0; - // declaration du tableau des coordonnees des sommets de la frontiere // puis des sommets internes ajoutes // majoration empirique du nombre de sommets de la triangulation - i = 4*nbarfr/10; + i = 4*nbarfr/10; mxsomm = Max( 20000, 64*nbpti+i*i ); - // MESSAGEE( "APTRTE: Depart de la triangulation avec " ); - // MESSAGEE( "nutysu=" << nutysu << " aretmx=" << aretmx << " mxsomm=" << mxsomm ); + MESSAGE( "APTRTE: Debut de la triangulation plane avec " ); + MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx << " mxsomm=" << mxsomm ); + MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes"); NEWDEPART: //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets @@ -190,23 +191,14 @@ void aptrte( Z nutysu, R aretmx, mnsoar = new Z[mosoar*mxsoar]; if( mnsoar==NULL ) goto ERREUR; //initialiser le tableau mnsoar pour le hachage des aretes -#ifdef DFORTRAN - INSOAR( mxsomm, mosoar, mxsoar, n1soar, mnsoar ); -#else - insoar_( mxsomm, mosoar, mxsoar, n1soar, mnsoar ); -#endif + insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar ); // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets if( mnarst!=NULL ) delete [] mnarst; mnarst = new Z[1+mxsomm]; if( mnarst==NULL ) goto ERREUR; n = 1+mxsomm; - -#ifdef DFORTRAN - AZEROI( n, mnarst ); -#else - azeroi_( n, mnarst ); -#endif + azeroi( n, mnarst ); // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier // ou no du point si interne forc'e par l'utilisateur @@ -214,11 +206,7 @@ void aptrte( Z nutysu, R aretmx, if( mnslig!=NULL ) delete [] mnslig; mnslig = new Z[mxsomm]; if( mnslig==NULL ) goto ERREUR; -#ifdef DFORTRAN - AZEROI( mxsomm, mnslig ); -#else - azeroi_( mxsomm, mnslig ); -#endif + azeroi( mxsomm, mnslig ); // initialisation des aretes frontalieres de la triangulation future // renumerotation des sommets des aretes des lignes pour la triangulation @@ -238,13 +226,13 @@ void aptrte( Z nutysu, R aretmx, ns0 = nudslf[n-1]; mnpxyd[ns0].x = uvslf[ns0].x; mnpxyd[ns0].y = uvslf[ns0].y; - mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction ); + mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction ); // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z); //carre de la longueur de l'arete 1 de la ligne fermee n - d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ); - d = d + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ; + d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ) + + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ; aremin = Min( aremin, d ); aremax = Max( aremax, d ); @@ -260,17 +248,13 @@ void aptrte( Z nutysu, R aretmx, //le numero n de la ligne du sommet et son numero ns1 dans la ligne mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1]; -#ifdef DFORTRAN - FASOAR( ns1, ns2, moins1, moins1, n, -#else - fasoar_( ns1, ns2, moins1, moins1, n, -#endif + fasoar( ns1, ns2, moins1, moins1, n, mosoar, mxsoar, n1soar, mnsoar, mnarst, noar0, ierr ); //pas de test sur ierr car pas de saturation possible a ce niveau //le pointeur dans le hachage sur la premiere arete de la ligne fermee n - mndalf[n] = noar0; + //mndalf[n] = noar0; //la nouvelle arete est la suivante de l'arete definie juste avant if( noar > 0 ) @@ -294,13 +278,13 @@ void aptrte( Z nutysu, R aretmx, ns = ns1 - 1; mnpxyd[ns].x = uvslf[ns].x; mnpxyd[ns].y = uvslf[ns].y; - mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction ); + mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction ); // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z); //carre de la longueur de l'arete - d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2); - d = d + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2); + d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2) + + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2); aremin = Min( aremin, d ); aremax = Max( aremax, d ); @@ -308,11 +292,7 @@ void aptrte( Z nutysu, R aretmx, mnslig[ns] = 1000000 * n + ns1-nudslf[n-1]; //ajout de l'arete dans la liste -#ifdef DFORTRAN - FASOAR( ns1, ns2, moins1, moins1, n, -#else - fasoar_( ns1, ns2, moins1, moins1, n, -#endif + fasoar( ns1, ns2, moins1, moins1, n, mosoar, mxsoar, n1soar, mnsoar, mnarst, noar, ierr ); //pas de test sur ierr car pas de saturation possible a ce niveau @@ -332,8 +312,8 @@ void aptrte( Z nutysu, R aretmx, aremax = sqrt( aremax ); //longueur maximale d'une arete aretmx = Min( aretmx, aremax ); //pour homogeneiser - // MESSAGEE("nutysu=" << nutysu << " aretmx=" << aretmx - // << " arete min=" << aremin << " arete max=" << aremax); + MESSAGE("nutysu=" << nutysu << " aretmx=" << aretmx + << " arete min=" << aremin << " arete max=" << aremax); //chainage des aretes frontalieres : la derniere arete frontaliere mnsoar[ mosoar * noar - mosoar + 5 ] = 0; @@ -357,7 +337,7 @@ void aptrte( Z nutysu, R aretmx, //les 2 coordonnees du point i de sommet nbs mnpxyd[ns1].x = uvpti[i].x; mnpxyd[ns1].y = uvpti[i].y; - mnpxyd[ns1].z = areteideale_( mnpxyd[ns1], direction ); + mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction ); //le numero i du point interne mnslig[ns1] = i+1; ns1++; @@ -373,18 +353,14 @@ void aptrte( Z nutysu, R aretmx, mxtree = 2 * mxsomm; NEWTREE: //en cas de saturation de l'un des tableaux, on boucle - // MESSAGEE( "Debut triangulation avec mxsomm=" << mxsomm ); + MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm ); if( mntree != NULL ) delete [] mntree; nbsomm = nbarpi; mntree = new Z[motree*(1+mxtree)]; if( mntree==NULL ) goto ERREUR; //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm -#ifdef DFORTRAN - TEAJTE( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr ); -#else - teajte_( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr ); -#endif + teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr ); comxmi[0].z=0; comxmi[1].z=0; @@ -393,13 +369,13 @@ void aptrte( Z nutysu, R aretmx, //saturation de letree => sa taille est augmentee et relance mxtree = mxtree * 2; ierr = 0; - // MESSAGEE( "Nouvelle valeur de mxtree=" << mxtree ); + MESSAGE( "Nouvelle valeur de mxtree=" << mxtree ); goto NEWTREE; } deltacpu_( d ); tcpu += d; - // MESSAGEE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" ); + MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" ); if( ierr != 0 ) goto ERREUR; //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes @@ -412,19 +388,15 @@ void aptrte( Z nutysu, R aretmx, mnqueu = new Z[mxqueu]; if( mnqueu==NULL) goto ERREUR; -#ifdef DFORTRAN - TEHOTE( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd, -#else - tehote_( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd, -#endif + tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mntree, mxqueu, mnqueu, ierr ); deltacpu_( d ); tcpu += d; - // MESSAGEE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE=" - // << d << " secondes"); + MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE=" + << d << " secondes"); if( ierr != 0 ) { //destruction du tableau auxiliaire et de l'arbre @@ -432,7 +404,7 @@ void aptrte( Z nutysu, R aretmx, { //letree sature mxtree = mxtree * 2; - // MESSAGEE( "Redemarrage avec la valeur de mxtree=" << mxtree ); + MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree ); ierr = 0; goto NEWTREE; } @@ -443,11 +415,7 @@ void aptrte( Z nutysu, R aretmx, // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets // et des points de la frontiere, des points internes imposes interieurs // ========================================================================== -#ifdef DFORTRAN - TETRTE( comxmi, aretmx, nbarpi, mxsomm, mnpxyd, -#else - tetrte_( comxmi, aretmx, nbarpi, mxsomm, mnpxyd, -#endif + tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd, mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar, moartr, mxartr, n1artr, mnartr, mnarst, ierr ); @@ -459,7 +427,7 @@ void aptrte( Z nutysu, R aretmx, //Temps calcul deltacpu_( d ); tcpu += d; - // MESSAGEE( "Temps de la triangulation des TE=" << d << " secondes" ); + MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" ); // ierr =0 si pas d'erreur // =1 si le tableau mnsoar est sature @@ -476,22 +444,16 @@ void aptrte( Z nutysu, R aretmx, // avec echange des 2 diagonales afin de rendre la triangulation delaunay // ====================================================================== // formation du chainage 6 des aretes internes a echanger eventuellement -#ifdef DFORTRAN - AISOAR( mosoar, mxsoar, mnsoar, na ); - TEDELA( mnpxyd, mnarst, -#else - aisoar_( mosoar, mxsoar, mnsoar, na ); - tedela_( mnpxyd, mnarst, -#endif - + aisoar( mosoar, mxsoar, mnsoar, na ); + tedela( mnpxyd, mnarst, mosoar, mxsoar, n1soar, mnsoar, na, moartr, mxartr, n1artr, mnartr, n ); - // MESSAGEE( "Nombre d'echanges des diagonales de 2 triangles=" << n ); + MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n ); deltacpu_( d ); tcpu += d; - // MESSAGEE("Temps de la triangulation Delaunay par echange des diagonales=" - // << d << " secondes"); + MESSAGE("Temps de la triangulation Delaunay par echange des diagonales=" + << d << " secondes"); //qualites de la triangulation actuelle qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, @@ -514,21 +476,17 @@ void aptrte( Z nutysu, R aretmx, mnarcf2 = new Z[mxarcf]; if( mnarcf2 == NULL ) goto ERREUR; -#ifdef DFORTRAN - TEREFR( nbarpi, mnpxyd, -#else - terefr_( nbarpi, mnpxyd, -#endif + terefr( nbarpi, mnpxyd, mosoar, mxsoar, n1soar, mnsoar, moartr, n1artr, mnartr, mnarst, mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2, n, ierr ); - // MESSAGEE( "Restauration de " << n << " aretes perdues de la frontiere" ); + MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" ); deltacpu_( d ); tcpu += d; - // MESSAGEE("Temps de la recuperation des aretes perdues de la frontiere=" - // << d << " secondes"); + MESSAGE("Temps de la recuperation des aretes perdues de la frontiere=" + << d << " secondes"); if( ierr != 0 ) goto ERREUR; @@ -559,11 +517,7 @@ void aptrte( Z nutysu, R aretmx, for (n=0; n0 si une erreur est survenue c 11 algorithme defaillant @@ -2958,8 +2960,8 @@ c equivalence (nosotr(1),ns1), (nosotr(2),ns2), % (nosotr(3),ns3) c -c le nombre de sommets de te supprimes - nbstsu = 0 +cccc le nombre de sommets de te supprimes +ccc nbstsu = 0 c c initialisation du chainage des aretes des cf => 0 arete de cf do 10 narete=1,mxsoar @@ -2990,7 +2992,8 @@ c recherche des triangles de sommet ns c ils doivent former un contour ferme de type etoile call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr, % mxarcf, nbtrcf, notrcf ) - if( nbtrcf .le. 0 ) then + if( nbtrcf .eq. 0 ) goto 100 + if( nbtrcf .lt. 0 ) then c erreur: impossible de trouver tous les triangles de sommet ns c seule une partie est a priori retrouvee nbtrcf = -nbtrcf @@ -3050,8 +3053,9 @@ c ========================================== % mxarcf, n1arcf, noarcf, % larmin, notrcf, liarcf, ierr ) if( ierr .eq. 0 ) then -c un sommet de te supprime de plus - nbstsu = nbstsu + 1 +cccc un sommet de te supprime de plus +ccc nbstsu = nbstsu + 1 + goto 100 else if( ierr .lt. 0 ) then c le sommet nste est externe donc non supprime c ou bien le sommet nste est le centre d'un cf dont toutes @@ -3064,18 +3068,21 @@ c erreur motivant un arret de la triangulation return endif c -c boucle jusqu'a obtenir une qualite suffisante -c si triangulation tres irreguliere => -c destruction de beaucoup de points internes -c les 2 variables suivantes brident ces destructions massives +cccc boucle jusqu'a obtenir une qualite suffisante +cccc si triangulation tres irreguliere => +cccc destruction de beaucoup de points internes +cccc les 2 variables suivantes brident ces destructions massives ccc nbsuns = nbsuns + 1 - quaopt = quaopt * 0.8 -ccc if( nbsuns .le. 5 ) goto 15 - goto 15 +ccc quaopt = quaopt * 0.8 +ccc if( nbsuns .lt. 5 ) goto 15 endif endif c 100 continue +c +c write(imprim,*)'retrait de',nbstsu, +c % ' sommets de te trop proches de la frontiere' + return end @@ -3370,7 +3377,10 @@ c l'arete suivante % mxtrcf, n1arcf, noarcf, % larmin, notrcf, nostbo, % ierr ) - if( ierr .lt. 0 ) then + if( ierr .eq. -543 ) then + ierr = 0 + goto 1000 + else if( ierr .lt. 0 ) then c le sommet ns est externe donc non supprime c ou bien le sommet ns est le centre d'un cf dont toutes c les aretes simples sont frontalieres @@ -3704,12 +3714,13 @@ c n1arcf : tableau (0:mxtrcf) auxiliaire d'entiers c noarcf : tableau (3,mxtrcf) auxiliaire d'entiers c larmin : tableau ( mxtrcf ) auxiliaire d'entiers c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -c auteur : alain perronnet analyse numerique paris upmc mai 1997 +c auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006 c....................................................................012 parameter (lchain=6) common / unites / lecteu, imprim, nunite(30) double precision pxyd(3,*) double precision ponder, ponde1, xbar, ybar, x, y, d, dmin, dmax + double precision surtd2 double precision d2d3(3,3) real origin(3), xyz(3) integer noartr(moartr,*), @@ -3745,7 +3756,8 @@ c les compteurs de passage sur les differents cas nbst8 = 0 c c coefficient de ponderation croissant avec les iterations - ponder = min( 1d0, ( 50 + (50*iter)/nbitaq ) * 0.01d0 ) + ponder = min( 1d0, 0.1d0 + iter * 0.9d0 / nbitaq ) +ccc 9 mars 2006 ponder = min( 1d0, ( 50 + (50*iter)/nbitaq ) * 0.01d0 ) ponde1 = 1d0 - ponder c c l'ordre du parcours dans le sens croissant ou decroissant @@ -3774,7 +3786,7 @@ c recherche des triangles de sommet ns c ils doivent former un contour ferme de type etoile call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr, % mxtrcf, nbtrcf, notrcf ) - if( nbtrcf .le. 0 ) goto 1000 + if( nbtrcf .le. 2 ) goto 1000 c c boucle sur les triangles qui forment une boule autour du sommet ns nbstbo = 0 @@ -3840,9 +3852,14 @@ c pas de modification de la topologie lors de la derniere iteration c ================================================================= if( iter .ge. nbitaq ) goto 200 c -c si la boule de ns contient 3 ou 4 triangles le sommet ns est detruit -c ==================================================================== - if( nbtrcf .le. 4 ) then +c si la boule de ns contient au plus 3 triangles +c => pas de changement de topologie +c ============================================== + if( nbtrcf .le. 3 ) goto 200 +c +c si la boule de ns contient 4 triangles le sommet ns est detruit +c =============================================================== + if( nbtrcf .eq. 4 ) then c c remise a -1 du chainage des aretes peripheriques de la boule ns noar = noar0 @@ -3861,7 +3878,10 @@ c l'arete suivante % mxtrcf, n1arcf, noarcf, % larmin, notrcf, nostbo, % ierr ) - if( ierr .lt. 0 ) then + if( ierr .eq. -543 ) then + ierr = 0 + goto 1000 + else if( ierr .lt. 0 ) then c le sommet ns est externe donc non supprime c ou bien le sommet ns est le centre d'un cf dont toutes c les aretes simples sont frontalieres @@ -3873,6 +3893,8 @@ c dans les 2 cas le sommet ns n'est pas supprime nbstsu = nbstsu + 1 else c erreur irrecuperable + write(imprim,*) + % 'teamqs: erreur1 irrecuperable en sortie te1stm' goto 9999 endif goto 1000 @@ -3910,6 +3932,9 @@ c l'arete suivante endif c c le point ns1 devient le milieu de l'arete ns-ns1 + x = pxyd(1,ns1) + y = pxyd(2,ns1) + d = pxyd(3,ns1) do 75 j=1,3 pxyd(j,ns1) = (pxyd(j,ns) + pxyd(j,ns1)) * 0.5d0 75 continue @@ -3931,16 +3956,22 @@ c suppression du point ns et mise en delaunay if( ierr .lt. 0 ) then c le sommet ns est externe donc non supprime c ou bien le sommet ns est le centre d'un cf dont toutes -c les aretes simples sont frontalieres -c dans les 2 cas le sommet ns n'est pas supprime +c les aretes simples sont frontalieres ou erreur +c dans les 3 cas le sommet ns n'est pas supprime +c restauration du sommet ns1 a son ancienne place + pxyd(1,ns1) = x + pxyd(2,ns1) = y + pxyd(3,ns1) = d ierr = 0 - goto 200 + goto 1000 else if( ierr .eq. 0 ) then nbstsu = nbstsu + 1 nbst5 = nbst5 + 1 goto 1000 else c erreur irrecuperable + write(imprim,*) + % 'teamqs: erreur2 irrecuperable en sortie te1stm' goto 9999 endif endif @@ -3983,9 +4014,9 @@ c le numero pxyd de ses 3 sommets c c ajout du nouveau barycentre if( nbsomm .ge. mxsomm ) then - write(imprim,*) 'saturation du tableau pxyd' + write(imprim,*) 'teamqs: saturation du tableau pxyd' c abandon de l'amelioration - goto 1100 + goto 9999 endif nbsomm = nbsomm + 1 do 120 i=1,3 @@ -4024,7 +4055,11 @@ c protection a ne pas modifier sinon erreur! % moartr, mxartr, n1artr, noartr, % noarst, % nosotr, ierr ) - if( ierr .ne. 0 ) goto 9999 + if( ierr .ne. 0 ) then + write(imprim,*) + % 'teamqs: erreur irrecuperable en sortie tr3str' + goto 9999 + endif 140 continue c nbst8 = nbst8 + 1 @@ -4041,10 +4076,56 @@ c simples de la boule du sommet ns 200 xbar = xbar / nbstbo ybar = ybar / nbstbo c -c ponderation pour eviter les degenerescenses +C DEBUT AJOUT 21/MAI/2005 +C PONDERATION POUR EVITER LES DEGENERESCENSES AVEC PROTECTION +C SI UN TRIANGLE DE SOMMET NS A UNE AIRE NEGATIVE APRES BARYCENTRAGE +C ALORS LE SOMMET NS N'EST PAS BOUGE +c +c protection des XY du point initial + xxx = pxyd(1,ns) + yyy = pxyd(2,ns) +c pxyd(1,ns) = ponde1 * pxyd(1,ns) + ponder * xbar pxyd(2,ns) = ponde1 * pxyd(2,ns) + ponder * ybar c +ccc write(imprim,*)'teamqs 200: ns=',ns,' ancien =',xxx,yyy +ccc write(imprim,*)'teamqs 200: ns=',ns,' nouveau=',pxyd(1,ns),pxyd(2,ns) +c + do 240 i=1,nbtrcf +c le numero de l'arete du triangle nt ne contenant pas le sommet ns + nt = notrcf(i) + do 220 na=1,3 +c le numero de l'arete na dans le tableau nosoar + noar = abs( noartr(na,nt) ) + if( nosoar(1,noar) .ne. ns .and. + % nosoar(2,noar) .ne. ns ) then + if( noartr(na,nt) .ge. 0 ) then + ns2 = nosoar(1,noar) + ns3 = nosoar(2,noar) + else + ns3 = nosoar(1,noar) + ns2 = nosoar(2,noar) + endif + goto 225 + endif + 220 continue + +c aire signee du triangle nt + 225 d = surtd2( pxyd(1,ns), pxyd(1,ns2), pxyd(1,ns3) ) + if( d .le. 0d0 ) then +ccc write(imprim,*),'iter=',iter, +ccc % ' Barycentrage au point ns=',ns, +ccc % ' XB=',pxyd(1,ns),' YB=',pxyd(2,ns), +ccc % ' => triangle avec AIRE<0 => Pt REMIS en X =',xxx, +ccc % ' Y =',yyy + pxyd(1,ns) = xxx + pxyd(2,ns) = yyy + goto 1000 + endif + 240 continue +C +C FIN AJOUT 21/MAI/2005 +c c les aretes chainees de la boule sont rendues delaunay 300 call tedela( pxyd, noarst, % mosoar, mxsoar, n1soar, nosoar, noar0, @@ -4052,11 +4133,9 @@ c les aretes chainees de la boule sont rendues delaunay c 1000 continue c -c trace de la triangulation actuelle et calcul de la qualite - 1100 continue -c -ccc write(imprim,11000) nbst4, nbst5, nbst8 -ccc11000 format( i7,' sommets de 4t', +ccc write(imprim,11000) iter, nbitaq, nbst4, nbst5, nbst8 +ccc11000 format( 'teamqs iter=',i2,' max iter=',i2,':', +ccc % i7,' sommets de 4t', ccc % i7,' sommets 5t+5t', ccc % i7,' sommets >7t' ) c @@ -4083,7 +4162,7 @@ c % comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign, % ierr ) c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -c but : amelioration de la qualite de la triangulation issue de teabr4 +c but : amelioration de la qualite de la triangulation c ----- c c entrees: @@ -4153,10 +4232,8 @@ c ============================================================== % mosoar, mxsoar, n1soar, nosoar, % moartr, mxartr, n1artr, noartr, % mxarcf, n1arcf, noarcf, larmin, notrcf, nostbo, - % nbstsu, ierr ) + % ierr ) if( ierr .ne. 0 ) goto 9999 -c write(imprim,*) 'retrait de',nbstsu, -c % ' sommets de te trop proches de la frontiere' c c ajustage des tailles moyennes des aretes avec ampli=1.34d0 entre c ampli/2 x taille_souhaitee et ampli x taille_souhaitee @@ -4171,18 +4248,18 @@ c ================================================================ % ierr ) if( ierr .ne. 0 ) goto 9999 c -c modification de la topologie autour des sommets frontaliers -c pour avoir un nombre de triangles egal a l'angle/60 degres -c et mise en triangulation delaunay locale -c =========================================================== - call teamsf( nutysu, - % noarst, mosoar, mxsoar, n1soar, nosoar, - % moartr, mxartr, n1artr, noartr, - % mxarcf, notrcf, nostbo, - % n1arcf, noarcf, larmin, - % comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign, - % ierr ) - if( ierr .ne. 0 ) goto 9999 +cccc modification de la topologie autour des sommets frontaliers +cccc pour avoir un nombre de triangles egal a l'angle/60 degres +cccc et mise en triangulation delaunay locale +cccc =========================================================== +ccc call teamsf( nutysu, +ccc % noarst, mosoar, mxsoar, n1soar, nosoar, +ccc % moartr, mxartr, n1artr, noartr, +ccc % mxarcf, notrcf, nostbo, +ccc % n1arcf, noarcf, larmin, +ccc % comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign, +ccc % ierr ) +ccc if( ierr .ne. 0 ) goto 9999 c c quelques iterations de barycentrage des points internes c modification de la topologie pour avoir 4 ou 5 ou 6 triangles @@ -5540,7 +5617,7 @@ c dans les 2 cas => retour sans modifs c >0 si une erreur est survenue c =11 algorithme defaillant c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -c auteur : alain perronnet analyse numerique paris upmc mars 1997 +c auteur : alain perronnet analyse numerique paris upmc mars 2006 c....................................................................012 parameter ( lchain=6, quamal=0.3) common / unites / lecteu,imprim,intera,nunite(29) @@ -5601,7 +5678,11 @@ c forme a partir des aretes des triangles de l'etoile du sommet nsasup % moartr, n1artr, noartr, % nbarcf, n1arcf, noarcf, % ierr ) - if( ierr .ne. 0 ) return + if( ierr .ne. 0 ) then +c modification de ierr pour continuer le calcul + ierr = -543 + return + endif c c ici le sommet nsasup appartient a aucune arete noarst( nsasup ) = 0 @@ -5655,7 +5736,7 @@ c mise en delaunay des aretes chainees call tedela( pxyd, noarst, % mosoar, mxsoar, n1soar, nosoar, liarcf(1), % moartr, mxartr, n1artr, noartr, modifs ) -ccc write(imprim,*) 'nombre echanges diagonales =',modifs +ccc write(imprim,*) 'nombre echanges diagonales =',modifs return end @@ -7247,11 +7328,12 @@ c n1aeoc : numero dans nosoar de la premiere arete simple de l'etoile c c sortie : c -------- -c nbtrar : 1 si arete ajoutee, 2 si arete double supprimee +c nbtrar : 1 si arete ajoutee, 2 si arete double supprimee, 0 si erreur c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c auteur : alain perronnet analyse numerique paris upmc mars 1997 c2345x7..............................................................012 parameter (lchain=6) + common / unites / lecteu, imprim, nunite(30) integer nosoar(mosoar,mxsoar) c c si l'arete n'appartient pas aux aretes de l'etoile naetoi @@ -7274,11 +7356,22 @@ c c arete double de l'etoile. elle est supprimee du chainage na0 = 0 na = n1aeoc + nbpass = 0 c parcours des aretes chainees jusqu'a trouver l'arete noar 10 if( na .ne. noar ) then c passage a la suivante na0 = na na = nosoar( lchain, na ) + if( na .le. 0 ) then + nbtrar = 0 + return + endif + nbpass = nbpass + 1 + if( nbpass .gt. 128 ) then + write(imprim,*)'Pb dans caetoi: boucle infinie evitee' + nbtrar = 0 + return + endif goto 10 endif c @@ -7353,6 +7446,7 @@ c 14 si les lignes fermees se coupent => donnees a revoir c 15 si une seule arete simple frontaliere c 16 si boucle infinie car toutes les aretes simples c de la boule sont frontalieres! +c 17 si boucle infinie dans caetoi c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c auteur : alain perronnet analyse numerique paris upmc mars 1997 c....................................................................012 @@ -7384,7 +7478,12 @@ c l'arete de nosoar a traiter noar = abs( noartr(j,nt) ) call caetoi( noar, mosoar, mxsoar, n1soar, nosoar, % n1aeoc, nbtrar ) -c si arete simple alors suppression du numero de triangle pour cette a + if( nbtrar .le. 0 ) then + ierr = 17 + return + endif +c si arete simple alors suppression du numero de triangle +c pour cette arete if( nbtrar .eq. 1 ) then if( nosoar(4,noar) .eq. nt ) then nosoar(4,noar) = nosoar(5,noar) -- 2.30.2