-// 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
//
// 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_;
}
//Retourne le temps CPU utilise en secondes
{
tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
- // MESSAGEE( "temps cpu=" << tempsec );
+ //MESSAGE( "temps cpu=" << tempsec );
}
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
// 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
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;
Z *mnarcf2=NULL;
Z *mnarcf3=NULL;
Z *mntrsu=NULL;
- Z *mndalf=NULL;
Z *mnslig=NULL;
Z *mnarst=NULL;
Z *mnlftr=NULL;
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;
// 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
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
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
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 );
//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 )
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 );
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
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;
//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++;
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;
//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
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
{
//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;
}
// 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 );
//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
// 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,
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;
for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
mnlftr[n] = n+1;
-#ifdef DFORTRAN
- TESUEX( nblf, mnlftr,
-#else
- tesuex_( nblf, mnlftr,
-#endif
+ tesuex( nblf, mnlftr,
ndtri0, nbsomm, mnpxyd, mnslig,
mosoar, mxsoar, mnsoar,
moartr, mxartr, n1artr, mnartr, mnarst,
deltacpu_( d );
tcpu += d;
- // MESSAGEE( "Temps de la suppression des triangles externes=" << d );
+ MESSAGE( "Temps de la suppression des triangles externes=" << d );
if( ierr != 0 ) goto ERREUR;
//qualites de la triangulation actuelle
// mise en delaunay de la triangulation
// =====================================================
mnarcf3 = new Z[mxarcf];
- if( mnarcf3 == NULL ) goto ERREUR;
-
-#ifdef DFORTRAN
- TEAMQT( nutysu,
-#else
- teamqt_( nutysu,
-#endif
+ if( mnarcf3 == NULL )
+ {
+ cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
+ 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,
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;
- // MESSAGEE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
+ MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
if( ierr != 0 ) goto ERREUR;
//qualites de la triangulation finale
if( mnartr[nt*moartr-moartr] != 0 )
{
//le numero des 3 sommets du triangle nt
-#ifdef DFORTRAN
- NUSOTR( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
-#else
- nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
-#endif
+ nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
//les 3 sommets du triangle sont actifs
mnarst[ nosotr[0] ] = 1;
mnarst[ nosotr[1] ] = 1;
// -----------------------------------------------------
// 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++)
if( mnartr[i*moartr-moartr] != 0 )
{
//le triangle i est interne => nosotr numero de ses 3 sommets
-#ifdef DFORTRAN
- NUSOTR( i, mosoar, mnsoar, moartr, mnartr, nosotr );
-#else
- nusotr_( i, mosoar, mnsoar, moartr, mnartr, nosotr );
-#endif
+ nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
nust[nbt++] = mnarst[ nosotr[0] ];
nust[nbt++] = mnarst[ nosotr[1] ];
nust[nbt++] = mnarst[ nosotr[2] ];
nust[nbt++] = 0;
}
}
- nbt /= 4; //le nombre final de triangles de la surface
- // MESSAGEE("Nombre de sommets=" << nbst
- // << " Nombre de triangles=" << nbt);
-
+ nbt /= nbsttria; //le nombre final de triangles de la surface
+ MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
+ << nbt << " triangles=" << nbt);
deltacpu_( d );
tcpu += d;
- // MESSAGEE( "Temps total de la triangulation=" << tcpu << " secondes" );
+ MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
// destruction des tableaux auxiliaires
// ------------------------------------
if( mnslig != NULL ) delete [] mnslig;
if( mnsoar != NULL ) delete [] mnsoar;
if( mnpxyd != NULL ) delete [] mnpxyd;
- if( mndalf != NULL ) delete [] mndalf;
if( mntree != NULL ) delete [] mntree;
if( mnqueu != NULL ) delete [] mnqueu;
if( mntrsu != NULL ) delete [] mntrsu;
}
else
{
- // MESSAGEE( "Triangulation non realisee " << ierr );
+ MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
if( ierr == 0 ) ierr=1;
goto NETTOYAGE;
}
nbtria++;
//le numero des 3 sommets du triangle nt
-#ifdef DFORTRAN
- NUSOTR( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
-#else
- nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
-#endif
+ nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
//la qualite du triangle ns1 ns2 ns3
-#ifdef DFORTRAN
- QUTR2D( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
-#else
- qutr2d_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
-#endif
+ qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
qualite );
//la qualite moyenne
quamin = Min( quamin, qualite );
//aire signee du triangle nt
-#ifdef DFORTRAN
- d = SURTD2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
-#else
- d = surtd2_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
-#endif
-
+ d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
if( d<0 )
{
//un triangle d'aire negative de plus
nbtrianeg++;
- // MESSAGEE("ATTENTION: le triangle " << nt << " de sommets:"
- // << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
- // << " a une aire " << d <<"<=0");
+ MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
+ << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
+ << " a une aire " << d <<"<=0");
}
//aire des triangles actuels
//les affichages
quamoy /= nbtria;
- // MESSAGEE("Qualite moyenne=" << quamoy
- // << " Qualite minimale=" << quamin
- // << " des " << nbtria << " triangles de surface totale="
- // << aire);
+ MESSAGE("Qualite moyenne=" << quamoy
+ << " Qualite minimale=" << quamin
+ << " des " << nbtria << " triangles de surface plane totale="
+ << aire);
- //if( nbtrianeg>0 )
- // MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );
+ if( nbtrianeg>0 )
+ MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );
return;
}
// auteur : Alain Perronnet Analyse Numerique Paris UPMC decembre 2001
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#if WIN32 & DFORTRAN
+ #define tempscpu TEMPSCPU
+ #define deltacpu DELTACPU
+ #define insoar INSOAR
+ #define azeroi AZEROI
+ #define fasoar FASOAR
+ #define teajte TEAJTE
+ #define tehote TEHOTE
+ #define tetrte TETRTE
+ #define aisoar AISOAR
+ #define tedela TEDELA
+ #define terefr TEREFR
+ #define tesuex TESUEX
+ #define teamqt TEAMQT
+ #define nusotr NUSOTR
+ #define qutr2d QUTR2D
+ #define surtd2 SURTD2
+
+ #define areteideale ARETEIDEALE
+
+#else
+ #define tempscpu tempscpu_
+ #define deltacpu deltacpu_
+ #define insoar insoar_
+ #define azeroi azeroi_
+ #define fasoar fasoar_
+ #define teajte teajte_
+ #define tehote tehote_
+ #define tetrte tetrte_
+ #define aisoar aisoar_
+ #define tedela tedela_
+ #define terefr terefr_
+ #define tesuex tesuex_
+ #define teamqt teamqt_
+ #define nusotr nusotr_
+ #define qutr2d qutr2d_
+ #define surtd2 surtd2_
+
+ #define areteideale areteideale_
+
+#endif
+
+
extern "C" { void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TEMPSCPU( double & tempsec ); }
-#else
- tempscpu_( double & tempsec ); }
-#endif
+ tempscpu( double & tempsec );
+}
//Retourne le temps CPU utilise en secondes
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- DELTACPU( R & dtcpu ); }
-#else
- deltacpu_( R & dtcpu ); }
-#endif
+ deltacpu( R & dtcpu );
+}
//Retourne le temps CPU utilise en secondes depuis le precedent appel
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- INSOAR( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );}
-#else
- insoar_( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );}
-#endif
+ insoar( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );
+}
//mettre a zero les nb entiers de tab
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- AZEROI( Z & nb, Z * tab );}
-#else
- azeroi_( Z & nb, Z * tab );}
-#endif
+ azeroi( Z & nb, Z * tab );
+}
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- FASOAR( Z & ns1, Z & ns2, Z & nt1, Z & nt2, Z & nolign,
-#else
- fasoar_( Z & ns1, Z & ns2, Z & nt1, Z & nt2, Z & nolign,
-#endif
- Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z * mnarst,
- Z & noar, Z & ierr );}
+ fasoar( Z & ns1, Z & ns2, Z & nt1, Z & nt2, Z & nolign,
+ Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z * mnarst,
+ Z & noar, Z & ierr );
+}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// but : former l'arete de sommet ns1-ns2 dans le hachage du tableau
// ----- nosoar des aretes de la triangulation
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TEAJTE
-#else
- teajte_
-#endif
- ( Z & mxsomm, Z & nbsomm, R3 * mnpxyd, R3 * comxmi,
- R & aretmx, Z & mxtree, Z * letree,
- Z & ierr );}
-
+ teajte( Z & mxsomm, Z & nbsomm, R3 * mnpxyd, R3 * comxmi,
+ R & aretmx, Z & mxtree, Z * letree,
+ Z & ierr );
+}
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TEHOTE
-#else
- tehote_
-#endif
- ( Z & nutysu, Z & nbarpi, Z & mxsomm, Z & nbsomm, R3 * mnpxyd,
- R3 * comxmi, R & aretmx,
- Z * letree, Z & mxqueu, Z * mnqueu,
- Z & ierr );}
+ tehote( Z & nutysu, Z & nbarpi, Z & mxsomm, Z & nbsomm, R3 * mnpxyd,
+ R3 * comxmi, R & aretmx,
+ Z * letree, Z & mxqueu, Z * mnqueu,
+ Z & ierr );
+}
// homogeneisation de l'arbre des te a un saut de taille au plus
// prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TETRTE
-#else
- tetrte_
-#endif
- ( R3 * comxmi, R & aretmx, Z & nbarpi, Z & mxsomm, R3 * mnpxyd,
- Z & mxqueu, Z * mnqueu, Z * mntree,
- Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
- Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
- Z & ierr );}
+ tetrte( R3 * comxmi, R & aretmx, Z & nbarpi, Z & mxsomm, R3 * mnpxyd,
+ Z & mxqueu, Z * mnqueu, Z * mntree,
+ Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
+ Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
+ Z & ierr );
+}
// trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
// et des points de la frontiere, des points internes imposes interieurs
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- AISOAR
-#else
- aisoar_
-#endif
- ( Z & mosoar, Z & mxsoar, Z * mnsoar, Z & na );}
- // formation du chainage 6 des aretes internes a echanger eventuellement
+ aisoar( Z & mosoar, Z & mxsoar, Z * mnsoar, Z & na );
+}
+// formation du chainage 6 des aretes internes a echanger eventuellement
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TEDELA
-#else
- tedela_
-#endif
- ( R3 * mnpxyd, Z * mnarst,
- Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & na,
- Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z & n );}
- // boucle sur les aretes internes (non sur une ligne de la frontiere)
- // avec echange des 2 diagonales afin de rendre la triangulation delaunay
+ tedela( R3 * mnpxyd, Z * mnarst,
+ Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & na,
+ Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z & n );
+}
+// boucle sur les aretes internes (non sur une ligne de la frontiere)
+// avec echange des 2 diagonales afin de rendre la triangulation delaunay
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TEREFR
-#else
- terefr_
-#endif
- ( Z & nbarpi, R3 * mnpxyd,
- Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
- Z & moartr, Z & n1artr, Z * mnartr, Z * mnarst,
- Z & mxarcf, Z * mnarc1, Z * mnarc2,
- Z * mnarc3, Z * mnarc4,
- Z & n, Z & ierr );}
+ terefr( Z & nbarpi, R3 * mnpxyd,
+ Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
+ Z & moartr, Z & n1artr, Z * mnartr, Z * mnarst,
+ Z & mxarcf, Z * mnarc1, Z * mnarc2,
+ Z * mnarc3, Z * mnarc4,
+ Z & n, Z & ierr );
+}
// detection des aretes frontalieres initiales perdues
// triangulation frontale pour les restaurer
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TESUEX
-#else
- tesuex_
-#endif
- ( Z & nblf, Z * nulftr,
- Z & ndtri0, Z & nbsomm, R3 * mnpxyd, Z * mnslig,
- Z & mosoar, Z & mxsoar, Z * mnsoar,
- Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
- Z & nbtria, Z * mntrsu, Z & ierr );}
+ tesuex( Z & nblf, Z * nulftr,
+ Z & ndtri0, Z & nbsomm, R3 * mnpxyd, Z * mnslig,
+ Z & mosoar, Z & mxsoar, Z * mnsoar,
+ Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
+ Z & nbtria, Z * mntrsu, Z & ierr );
+}
// suppression des triangles externes a la surface
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- TEAMQT
-#else
- teamqt_
-#endif
- ( Z & nutysu,
- Z * mnarst, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
- Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr,
- Z & mxarcf, Z * mntrcf, Z * mnstbo,
- Z * n1arcf, Z * mnarcf, Z * mnarc1,
- R3 * comxmi, Z & nbarpi, Z & nbsomm, Z & mxsomm,
- R3 * mnpxyd, Z * mnslig,
- Z & ierr );}
+ teamqt( Z & nutysu,
+ Z * mnarst, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
+ Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr,
+ Z & mxarcf, Z * mntrcf, Z * mnstbo,
+ Z * n1arcf, Z * mnarcf, Z * mnarc1,
+ R3 * comxmi, Z & nbarpi, Z & nbsomm, Z & mxsomm,
+ R3 * mnpxyd, Z * mnslig,
+ Z & ierr );
+}
// amelioration de la qualite de la triangulation par
// barycentrage des sommets internes a la triangulation
// suppression des aretes trop longues ou trop courtes
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- NUSOTR
-#else
- nusotr_
-#endif
- ( Z & nt, Z & mosoar, Z * mnsoar, Z & moartr, Z * mnartr,Z * nosotr );}
+ nusotr( Z & nt, Z & mosoar, Z * mnsoar, Z & moartr, Z * mnartr,Z * nosotr );
+}
//retrouver les numero des 3 sommets du triangle nt
extern "C" {void
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- QUTR2D
-#else
- qutr2d_
-#endif
- ( R3 & p1, R3 & p2, R3 & p3, R & qualite );}
+ qutr2d( R3 & p1, R3 & p2, R3 & p3, R & qualite );
+}
//calculer la qualite d'un triangle de R2 de sommets p1, p2, p3
extern "C" { R
#ifdef WIN32
__stdcall
#endif
-#ifdef DFORTRAN
- SURTD2
-#else
- surtd2_
-#endif
- ( R3 & p1, R3 & p2, R3 & p3 ); }
+ surtd2( R3 & p1, R3 & p2, R3 & p3 );
+}
//calcul de la surface d'un triangle defini par 3 points de r**2
#endif
-c MEFISTO : library to compute 2D triangulation from segmented boundaries
+c MEFISTO2: a library to compute 2D triangulation from segmented boundaries
c
-c Copyright (C) 2003 Laboratoire J.-L. Lions UPMC Paris
+c Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris
c
c This library is free software; you can redistribute it and/or
c modify it under the terms of the GNU Lesser General Public
c License along with this library; if not, write to the Free Software
c Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
c
-c See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
+c See http://www.ann.jussieu.fr/~perronnet or email perronnet@ann.jussieu.fr
c
-c
-c File : trte.f
+c File : trte.f le Fortran du trianguleur plan
c Module : SMESH
-c Author: Alain PERRONNET
+c Author : Alain PERRONNET
+c Date : 16 mars 2006
subroutine qutr2d( p1, p2, p3, qualite )
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xyzd(2) = 0d0
xyzd(3) = 0d0
- longai = areteideale(xyz,xyzd)
+ longai = areteideale()
+c (xyz,xyzd)
if( longai .lt. 0d0 ) then
write(imprim,10000) xyz
10000 format('attention: longueur de areteideale(',
% noarst(*)
integer lapile(1:mxpile)
integer nosotr(3)
+c
+ lhpile = 0
c
c la premiere arete de sommet ns
nar = noarst( ns )
c
c saturation de la pile des triangles
c -----------------------------------
- 9990 write(imprim,*) 'trp1st:saturation pile des triangles autour ',
+ 9990 write(imprim,*)'trp1st: saturation pile des triangles autour ',
%'sommet',ns
goto 9999
c
c erreur triangle ne contenant pas le sommet ns
c ----------------------------------------------
- 9995 write(imprim,*) 'trp1st:triangle ',nta,' st=',
+ 9995 write(imprim,*) 'trp1st: triangle ',nta,' st=',
% (nosotr(nar),nar=1,3),' sans le sommet' ,ns
c
9999 lhpile = 0
% mosoar, mxsoar, n1soar, nosoar,
% moartr, mxartr, n1artr, noartr,
% mxarcf, n1arcf, noarcf, larmin, notrcf, liarcf,
- % nbstsu, ierr )
+ % ierr )
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c but : supprimer de la triangulation les sommets de te trop proches
c ----- soit d'un sommet frontalier ou point interne impose
c
c sortie :
c --------
-c nbstsu : nombre de sommets de te supprimes
c ierr : =0 si pas d'erreur
c >0 si une erreur est survenue
c 11 algorithme defaillant
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
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
% 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
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
% 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
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,*),
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
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
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
% 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
nbstsu = nbstsu + 1
else
c erreur irrecuperable
+ write(imprim,*)
+ % 'teamqs: erreur1 irrecuperable en sortie te1stm'
goto 9999
endif
goto 1000
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
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
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
% 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
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,
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
% 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:
% 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
% 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
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)
% 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
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
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
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
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
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)
