X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEFISTO2%2Faptrte.cxx;h=4c550a880434507ac8356f6c8f17df3e7a8259a5;hb=b2fe6ee8896fca4876f426c77f6e20a7a843e0af;hp=58da98168175e2835bd00c924cc51a64c728ab1f;hpb=bef9beee88cac57394b8dc3bc914381c1a2fff83;p=modules%2Fsmesh.git diff --git a/src/MEFISTO2/aptrte.cxx b/src/MEFISTO2/aptrte.cxx index 58da98168..4c550a880 100755 --- a/src/MEFISTO2/aptrte.cxx +++ b/src/MEFISTO2/aptrte.cxx @@ -1,8 +1,35 @@ -using namespace std; +// MEFISTO2: a library to compute 2D triangulation from segmented boundaries +// +// +// 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 +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr +// +// File : aptrte.cxx le C++ de l'appel du trianguleur plan +// Module : SMESH +// Author : Alain PERRONNET +// Date : 13 novembre 2006 + #include "Rn.h" #include "aptrte.h" #include "utilities.h" +using namespace std; + extern "C" { R aretemaxface_; @@ -37,10 +64,11 @@ void deltacpu_( R & dtcpu ) } -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 @@ -84,9 +112,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 @@ -95,7 +126,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; @@ -103,21 +134,20 @@ 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; R3 comxmi[2]; //coordonnees UV Min et Maximales R aremin, aremax; //longueur minimale et maximale des aretes + R airemx; //aire maximale souhaitee d'un triangle R quamoy, quamin; Z noar0, noar, na; Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt; Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn; Z moins1=-1; - - aretemaxface_ = aretmx; + Z nuds = 0; // initialisation du temps cpu deltacpu_( d ); @@ -125,19 +155,15 @@ 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 ); - MESSAGE( "APTRTE: Depart de la triangulation avec " ); - MESSAGE( "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 @@ -198,8 +224,8 @@ void aptrte( Z nutysu, R aretmx, mnpxyd[ns0].x = uvslf[ns0].x; mnpxyd[ns0].y = uvslf[ns0].y; mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction ); -// cout << "Sommet " << ns0 << ": " << mnpxyd[ns0].x -// << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z << endl; +// 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 ) @@ -225,7 +251,7 @@ void aptrte( Z nutysu, R aretmx, //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 ) @@ -247,11 +273,14 @@ void aptrte( Z nutysu, R aretmx, //l'arete precedente est dotee de sa suivante:celle cree ensuite //les 2 coordonnees du sommet ns2 de la ligne ns = ns1 - 1; +//debut ajout 5/10/2006 ................................................ + nuds = Max( nuds, ns ); //le numero du dernier sommet traite +//fin ajout 5/10/2006 ................................................ mnpxyd[ns].x = uvslf[ns].x; mnpxyd[ns].y = uvslf[ns].y; mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction ); -// cout << "Sommet " << ns << ": " << mnpxyd[ns].x -// << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z << endl; +// 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) @@ -259,6 +288,14 @@ void aptrte( Z nutysu, R aretmx, aremin = Min( aremin, d ); aremax = Max( aremax, d ); +//debut ajout du 5/10/2006 ............................................. + //la longueur de l'arete ns1-ns2 + d = sqrt( d ); + //longueur arete = Min ( aretmx, aretes incidentes ) + mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d ); + mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d ); +//fin ajout du 5/10/2006 ............................................... + //le numero n de la ligne du sommet et son numero ns1 dans la ligne mnslig[ns] = 1000000 * n + ns1-nudslf[n-1]; @@ -282,9 +319,34 @@ void aptrte( Z nutysu, R aretmx, aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees aremax = sqrt( aremax ); //longueur maximale d'une arete - aretmx = Min( aretmx, aremax ); //pour homogeneiser - cout << "nutysu=" << nutysu << " aretmx=" << aretmx - << " arete min=" << aremin << " arete max=" << aremax << endl; +//debut ajout 9/11/2006 ................................................ + // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser + + // protection contre une arete max desiree trop grande ou trop petite + if( aretmx > aremax*2.05 ) aretmx = aremax; + + // protection contre une arete max desiree trop petite + if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 ) + aretmx =(aremin+aremax*2)/3.0; + + if( aretmx < aremin && aremin > 0 ) + aretmx = aremin; + + //sauvegarde pour la fonction areteideale_ + aretemaxface_ = aretmx; + + //aire maximale souhaitee des triangles + airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral + + for(i=0; i<=nuds; i++ ) + mnpxyd[i].z = Min( mnpxyd[i].z, aretmx ); + //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1); +//fin ajout 9/11/2006 ................................................. + + + MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax ); + MESSAGE("Triangulation: arete mx=" << aretmx + << " triangle aire mx=" << airemx ); //chainage des aretes frontalieres : la derniere arete frontaliere mnsoar[ mosoar * noar - mosoar + 5 ] = 0; @@ -366,8 +428,8 @@ void aptrte( Z nutysu, R aretmx, deltacpu_( d ); tcpu += d; - cout << "Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE=" - << d << " secondes" << endl; + 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 @@ -408,8 +470,8 @@ void aptrte( Z nutysu, R aretmx, if( ierr != 0 ) goto ERREUR; //qualites de la triangulation actuelle - qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, - nbt, quamoy, quamin ); + qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, + nbt, quamoy, quamin ); // boucle sur les aretes internes (non sur une ligne de la frontiere) // avec echange des 2 diagonales afin de rendre la triangulation delaunay @@ -423,12 +485,12 @@ void aptrte( Z nutysu, R aretmx, MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n ); deltacpu_( d ); tcpu += d; - cout << "Temps de la triangulation Delaunay par echange des diagonales=" - << d << " secondes" << endl; + 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, - nbt, quamoy, quamin ); + qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, + nbt, quamoy, quamin ); // detection des aretes frontalieres initiales perdues // triangulation frontale pour les restaurer @@ -449,21 +511,21 @@ void aptrte( Z nutysu, R aretmx, terefr_( nbarpi, mnpxyd, mosoar, mxsoar, n1soar, mnsoar, - moartr, n1artr, mnartr, mnarst, + moartr, mxartr, n1artr, mnartr, mnarst, mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2, n, ierr ); - MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" ); + MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr ); deltacpu_( d ); tcpu += d; - cout << "Temps de la recuperation des aretes perdues de la frontiere=" - << d << " secondes" << endl; + MESSAGE("Temps de la recuperation des aretes perdues de la frontiere=" + << d << " secondes"); if( ierr != 0 ) goto ERREUR; //qualites de la triangulation actuelle - qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, - nbt, quamoy, quamin ); + qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, + nbt, quamoy, quamin ); // fin de la triangulation avec respect des aretes initiales frontalieres @@ -499,12 +561,12 @@ void aptrte( Z nutysu, R aretmx, deltacpu_( d ); tcpu += d; - MESSAGE( "Temps de la suppression des triangles externes=" << d ); + MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr ); if( ierr != 0 ) goto ERREUR; //qualites de la triangulation actuelle - qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, - nbt, quamoy, quamin ); + qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, + nbt, quamoy, quamin ); // amelioration de la qualite de la triangulation par // barycentrage des sommets internes a la triangulation @@ -513,29 +575,33 @@ void aptrte( Z nutysu, R aretmx, // mise en delaunay de la triangulation // ===================================================== mnarcf3 = new Z[mxarcf]; - if( mnarcf3 == NULL ) goto ERREUR; - - teamqt_( nutysu, - mnarst, mosoar, mxsoar, n1soar, mnsoar, - moartr, mxartr, n1artr, mnartr, - mxarcf, mnarcf2, mnarcf3, - mn1arcf, mnarcf, mnarcf1, - comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig, + if( mnarcf3 == NULL ) + { + cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl; + goto ERREUR; + } + teamqt_( nutysu, aretmx, airemx, + mnarst, mosoar, mxsoar, n1soar, mnsoar, + moartr, mxartr, n1artr, mnartr, + mxarcf, mnarcf2, mnarcf3, + mn1arcf, mnarcf, mnarcf1, + nbarpi, nbsomm, mxsomm, mnpxyd, mnslig, ierr ); + if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;} if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;} if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;} if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;} if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;} - if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;} deltacpu_( d ); tcpu += d; MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d ); - if( ierr != 0 ) goto ERREUR; + if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi + if( ierr != 0 ) goto ERREUR; //qualites de la triangulation finale - qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, - nbt, quamoy, quamin ); + qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, + nbt, quamoy, quamin ); // renumerotation des sommets internes: mnarst(i)=numero final du sommet // =================================== @@ -606,7 +672,7 @@ void aptrte( Z nutysu, R aretmx, // ----------------------------------------------------- // boucle sur les triangles occupes (internes et externes) if( nust != NULL ) delete [] nust; - nust = new Z[4*nbt]; + nust = new Z[nbsttria*nbt]; if( nust == NULL ) goto ERREUR; nbt = 0; for (i=1; i<=mxartr; i++) @@ -622,13 +688,12 @@ void aptrte( Z nutysu, R aretmx, nust[nbt++] = 0; } } - nbt /= 4; //le nombre final de triangles de la surface - cout << "Nombre de sommets=" << nbst - << " Nombre de triangles=" << nbt << endl; - + nbt /= nbsttria; //le nombre final de triangles de la surface + MESSAGE( "APTRTE: Fin de la triangulation plane avec "<