ForkMaster/smesh
1
0
Fork 0
mirror of https://git.salome-platform.org/gitpub/modules/smesh.git synced 2025-01-14 21:20:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

Compile new version of MEFISTO with DF6.0

Browse Source
This commit is contained in:
abd 2006-04-26 10:41:24 +00:00
parent cf91f3e22a
commit 1db8a57bb1
3 changed files with 362 additions and 350 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

246
src/MEFISTO2/aptrte.cxx
View File

@ -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 // This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public // 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 // See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
// //
// // File : aptrte.cxx le C++ de l'appel du trianguleur plan
// File : aptrte.cxx
// Module : SMESH // Module : SMESH
// Author: Alain PERRONNET // Author : Alain PERRONNET
// Date : 16 mars 2006
#include "Rn.h" #include "Rn.h"
#include "aptrte.h" #include "aptrte.h"
#ifndef WIN32
#include "utilities.h" #include "utilities.h"
using namespace std; using namespace std;
#endif
extern "C" extern "C"
{ {
R aretemaxface_; R aretemaxface_;
R areteideale_( R3 xyz, R3 direction ) MEFISTO2D_EXPORT
R
#ifdef WIN32
__stdcall
#endif
areteideale()//( R3 xyz, R3 direction )
{ {
return aretemaxface_; return aretemaxface_;
} }
@ -50,7 +53,7 @@ void tempscpu_( double & tempsec )
//Retourne le temps CPU utilise en secondes //Retourne le temps CPU utilise en secondes
{ {
tempsec = ( (double) clock() ) / CLOCKS_PER_SEC; tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
// MESSAGEE( "temps cpu=" << tempsec ); //MESSAGE( "temps cpu=" << tempsec );
} }
@ -60,15 +63,15 @@ void deltacpu_( R & dtcpu )
tempscpu_( cpunew ); tempscpu_( cpunew );
dtcpu = R( cpunew - cpuold ); dtcpu = R( cpunew - cpuold );
cpuold = cpunew; cpuold = cpunew;
// MESSAGEE( "delta temps cpu=" << dtcpu ); //MESSAGE( "delta temps cpu=" << dtcpu );
return; return;
} }
void aptrte( Z nutysu, R aretmx, void aptrte( Z nutysu, R aretmx,
Z nblf, Z * nudslf, R2 * uvslf, Z nblf, Z * nudslf, R2 * uvslf,
Z nbpti, R2 *uvpti, Z nbpti, R2 * uvpti,
Z & nbst, R2 * & uvst, Z & nbt, Z * & nust, Z & nbst, R2 * & uvst,
Z & nbt, Z * & nust,
Z & ierr ) Z & ierr )
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// but : appel de la triangulation par un arbre-4 recouvrant // 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 // ierr : 0 si pas d'erreur
// > 0 sinon // > 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; R d, tcpu=0;
R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici! R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
@ -131,7 +137,6 @@ void aptrte( Z nutysu, R aretmx,
Z *mnarcf2=NULL; Z *mnarcf2=NULL;
Z *mnarcf3=NULL; Z *mnarcf3=NULL;
Z *mntrsu=NULL; Z *mntrsu=NULL;
Z *mndalf=NULL;
Z *mnslig=NULL; Z *mnslig=NULL;
Z *mnarst=NULL; Z *mnarst=NULL;
Z *mnlftr=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 i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn; Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
Z moins1=-1; Z moins1=-1;
R dist;
aretemaxface_ = aretmx; aretemaxface_ = aretmx;
@ -153,19 +159,14 @@ void aptrte( Z nutysu, R aretmx,
// quelques reservations de tableaux pour faire les calculs // 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 // declaration du tableau des coordonnees des sommets de la frontiere
// puis des sommets internes ajoutes // puis des sommets internes ajoutes
// majoration empirique du nombre de sommets de la triangulation // majoration empirique du nombre de sommets de la triangulation
i = 4*nbarfr/10; i = 4*nbarfr/10;
mxsomm = Max( 20000, 64*nbpti+i*i ); mxsomm = Max( 20000, 64*nbpti+i*i );
// MESSAGEE( "APTRTE: Depart de la triangulation avec " ); MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
// MESSAGEE( "nutysu=" << nutysu << " aretmx=" << aretmx << " mxsomm=" << mxsomm ); MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx << " mxsomm=" << mxsomm );
MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
NEWDEPART: NEWDEPART:
//mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets //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]; mnsoar = new Z[mosoar*mxsoar];
if( mnsoar==NULL ) goto ERREUR; if( mnsoar==NULL ) goto ERREUR;
//initialiser le tableau mnsoar pour le hachage des aretes //initialiser le tableau mnsoar pour le hachage des aretes
#ifdef DFORTRAN insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
INSOAR( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
#else
insoar_( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
#endif
// mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
if( mnarst!=NULL ) delete [] mnarst; if( mnarst!=NULL ) delete [] mnarst;
mnarst = new Z[1+mxsomm]; mnarst = new Z[1+mxsomm];
if( mnarst==NULL ) goto ERREUR; if( mnarst==NULL ) goto ERREUR;
n = 1+mxsomm; n = 1+mxsomm;
azeroi( n, mnarst );
#ifdef DFORTRAN
AZEROI( n, mnarst );
#else
azeroi_( n, mnarst );
#endif
// mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
// ou no du point si interne forc'e par l'utilisateur // 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; if( mnslig!=NULL ) delete [] mnslig;
mnslig = new Z[mxsomm]; mnslig = new Z[mxsomm];
if( mnslig==NULL ) goto ERREUR; if( mnslig==NULL ) goto ERREUR;
#ifdef DFORTRAN azeroi( mxsomm, mnslig );
AZEROI( mxsomm, mnslig );
#else
azeroi_( mxsomm, mnslig );
#endif
// initialisation des aretes frontalieres de la triangulation future // initialisation des aretes frontalieres de la triangulation future
// renumerotation des sommets des aretes des lignes pour la triangulation // renumerotation des sommets des aretes des lignes pour la triangulation
@ -238,13 +226,13 @@ void aptrte( Z nutysu, R aretmx,
ns0 = nudslf[n-1]; ns0 = nudslf[n-1];
mnpxyd[ns0].x = uvslf[ns0].x; mnpxyd[ns0].x = uvslf[ns0].x;
mnpxyd[ns0].y = uvslf[ns0].y; 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 // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
// << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z); // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
//carre de la longueur de l'arete 1 de la ligne fermee n //carre de la longueur de l'arete 1 de la ligne fermee n
d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ); d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
d = d + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ; + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
aremin = Min( aremin, d ); aremin = Min( aremin, d );
aremax = Max( aremax, 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 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1]; mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
#ifdef DFORTRAN fasoar( ns1, ns2, moins1, moins1, n,
FASOAR( ns1, ns2, moins1, moins1, n,
#else
fasoar_( ns1, ns2, moins1, moins1, n,
#endif
mosoar, mxsoar, n1soar, mnsoar, mnarst, mosoar, mxsoar, n1soar, mnsoar, mnarst,
noar0, ierr ); noar0, ierr );
//pas de test sur ierr car pas de saturation possible a ce niveau //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 //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 //la nouvelle arete est la suivante de l'arete definie juste avant
if( noar > 0 ) if( noar > 0 )
@ -294,13 +278,13 @@ void aptrte( Z nutysu, R aretmx,
ns = ns1 - 1; ns = ns1 - 1;
mnpxyd[ns].x = uvslf[ns].x; mnpxyd[ns].x = uvslf[ns].x;
mnpxyd[ns].y = uvslf[ns].y; 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 // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
// << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z); // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
//carre de la longueur de l'arete //carre de la longueur de l'arete
d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2); d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
d = d + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2); + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
aremin = Min( aremin, d ); aremin = Min( aremin, d );
aremax = Max( aremax, d ); aremax = Max( aremax, d );
@ -308,11 +292,7 @@ void aptrte( Z nutysu, R aretmx,
mnslig[ns] = 1000000 * n + ns1-nudslf[n-1]; mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
//ajout de l'arete dans la liste //ajout de l'arete dans la liste
#ifdef DFORTRAN fasoar( ns1, ns2, moins1, moins1, n,
FASOAR( ns1, ns2, moins1, moins1, n,
#else
fasoar_( ns1, ns2, moins1, moins1, n,
#endif
mosoar, mxsoar, n1soar, mnsoar, mosoar, mxsoar, n1soar, mnsoar,
mnarst, noar, ierr ); mnarst, noar, ierr );
//pas de test sur ierr car pas de saturation possible a ce niveau //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 aremax = sqrt( aremax ); //longueur maximale d'une arete
aretmx = Min( aretmx, aremax ); //pour homogeneiser aretmx = Min( aretmx, aremax ); //pour homogeneiser
// MESSAGEE("nutysu=" << nutysu << " aretmx=" << aretmx MESSAGE("nutysu=" << nutysu << " aretmx=" << aretmx
// << " arete min=" << aremin << " arete max=" << aremax); << " arete min=" << aremin << " arete max=" << aremax);
//chainage des aretes frontalieres : la derniere arete frontaliere //chainage des aretes frontalieres : la derniere arete frontaliere
mnsoar[ mosoar * noar - mosoar + 5 ] = 0; 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 //les 2 coordonnees du point i de sommet nbs
mnpxyd[ns1].x = uvpti[i].x; mnpxyd[ns1].x = uvpti[i].x;
mnpxyd[ns1].y = uvpti[i].y; 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 //le numero i du point interne
mnslig[ns1] = i+1; mnslig[ns1] = i+1;
ns1++; ns1++;
@ -373,18 +353,14 @@ void aptrte( Z nutysu, R aretmx,
mxtree = 2 * mxsomm; mxtree = 2 * mxsomm;
NEWTREE: //en cas de saturation de l'un des tableaux, on boucle 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; if( mntree != NULL ) delete [] mntree;
nbsomm = nbarpi; nbsomm = nbarpi;
mntree = new Z[motree*(1+mxtree)]; mntree = new Z[motree*(1+mxtree)];
if( mntree==NULL ) goto ERREUR; if( mntree==NULL ) goto ERREUR;
//initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
#ifdef DFORTRAN teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
TEAJTE( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
#else
teajte_( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
#endif
comxmi[0].z=0; comxmi[0].z=0;
comxmi[1].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 //saturation de letree => sa taille est augmentee et relance
mxtree = mxtree * 2; mxtree = mxtree * 2;
ierr = 0; ierr = 0;
// MESSAGEE( "Nouvelle valeur de mxtree=" << mxtree ); MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
goto NEWTREE; goto NEWTREE;
} }
deltacpu_( d ); deltacpu_( d );
tcpu += 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; if( ierr != 0 ) goto ERREUR;
//ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes //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]; mnqueu = new Z[mxqueu];
if( mnqueu==NULL) goto ERREUR; if( mnqueu==NULL) goto ERREUR;
#ifdef DFORTRAN tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
TEHOTE( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
#else
tehote_( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
#endif
comxmi, aretmx, comxmi, aretmx,
mntree, mxqueu, mnqueu, mntree, mxqueu, mnqueu,
ierr ); ierr );
deltacpu_( d ); deltacpu_( d );
tcpu += d; tcpu += d;
// MESSAGEE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE=" MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
// << d << " secondes"); << d << " secondes");
if( ierr != 0 ) if( ierr != 0 )
{ {
//destruction du tableau auxiliaire et de l'arbre //destruction du tableau auxiliaire et de l'arbre
@ -432,7 +404,7 @@ void aptrte( Z nutysu, R aretmx,
{ {
//letree sature //letree sature
mxtree = mxtree * 2; mxtree = mxtree * 2;
// MESSAGEE( "Redemarrage avec la valeur de mxtree=" << mxtree ); MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
ierr = 0; ierr = 0;
goto NEWTREE; goto NEWTREE;
} }
@ -443,11 +415,7 @@ void aptrte( Z nutysu, R aretmx,
// trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
// et des points de la frontiere, des points internes imposes interieurs // et des points de la frontiere, des points internes imposes interieurs
// ========================================================================== // ==========================================================================
#ifdef DFORTRAN tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
TETRTE( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
#else
tetrte_( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
#endif
mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar, mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
moartr, mxartr, n1artr, mnartr, mnarst, moartr, mxartr, n1artr, mnartr, mnarst,
ierr ); ierr );
@ -459,7 +427,7 @@ void aptrte( Z nutysu, R aretmx,
//Temps calcul //Temps calcul
deltacpu_( d ); deltacpu_( d );
tcpu += 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 // ierr =0 si pas d'erreur
// =1 si le tableau mnsoar est sature // =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 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
// ====================================================================== // ======================================================================
// formation du chainage 6 des aretes internes a echanger eventuellement // formation du chainage 6 des aretes internes a echanger eventuellement
#ifdef DFORTRAN aisoar( mosoar, mxsoar, mnsoar, na );
AISOAR( mosoar, mxsoar, mnsoar, na ); tedela( mnpxyd, mnarst,
TEDELA( mnpxyd, mnarst,
#else
aisoar_( mosoar, mxsoar, mnsoar, na );
tedela_( mnpxyd, mnarst,
#endif
mosoar, mxsoar, n1soar, mnsoar, na, mosoar, mxsoar, n1soar, mnsoar, na,
moartr, mxartr, n1artr, mnartr, n ); 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 ); deltacpu_( d );
tcpu += d; tcpu += d;
// MESSAGEE("Temps de la triangulation Delaunay par echange des diagonales=" MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
// << d << " secondes"); << d << " secondes");
//qualites de la triangulation actuelle //qualites de la triangulation actuelle
qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
@ -514,21 +476,17 @@ void aptrte( Z nutysu, R aretmx,
mnarcf2 = new Z[mxarcf]; mnarcf2 = new Z[mxarcf];
if( mnarcf2 == NULL ) goto ERREUR; if( mnarcf2 == NULL ) goto ERREUR;
#ifdef DFORTRAN terefr( nbarpi, mnpxyd,
TEREFR( nbarpi, mnpxyd,
#else
terefr_( nbarpi, mnpxyd,
#endif
mosoar, mxsoar, n1soar, mnsoar, mosoar, mxsoar, n1soar, mnsoar,
moartr, n1artr, mnartr, mnarst, moartr, n1artr, mnartr, mnarst,
mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2, mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
n, ierr ); n, ierr );
// MESSAGEE( "Restauration de " << n << " aretes perdues de la frontiere" ); MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" );
deltacpu_( d ); deltacpu_( d );
tcpu += d; tcpu += d;
// MESSAGEE("Temps de la recuperation des aretes perdues de la frontiere=" MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
// << d << " secondes"); << d << " secondes");
if( ierr != 0 ) goto ERREUR; if( ierr != 0 ) goto ERREUR;
@ -559,11 +517,7 @@ void aptrte( Z nutysu, R aretmx,
for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
mnlftr[n] = n+1; mnlftr[n] = n+1;
#ifdef DFORTRAN tesuex( nblf, mnlftr,
TESUEX( nblf, mnlftr,
#else
tesuex_( nblf, mnlftr,
#endif
ndtri0, nbsomm, mnpxyd, mnslig, ndtri0, nbsomm, mnpxyd, mnslig,
mosoar, mxsoar, mnsoar, mosoar, mxsoar, mnsoar,
moartr, mxartr, n1artr, mnartr, mnarst, moartr, mxartr, n1artr, mnartr, mnarst,
@ -574,7 +528,7 @@ void aptrte( Z nutysu, R aretmx,
deltacpu_( d ); deltacpu_( d );
tcpu += 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; if( ierr != 0 ) goto ERREUR;
//qualites de la triangulation actuelle //qualites de la triangulation actuelle
@ -588,28 +542,27 @@ void aptrte( Z nutysu, R aretmx,
// mise en delaunay de la triangulation // mise en delaunay de la triangulation
// ===================================================== // =====================================================
mnarcf3 = new Z[mxarcf]; mnarcf3 = new Z[mxarcf];
if( mnarcf3 == NULL ) goto ERREUR; if( mnarcf3 == NULL )
{
#ifdef DFORTRAN cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
TEAMQT( nutysu, goto ERREUR;
#else }
teamqt_( nutysu, teamqt( nutysu,
#endif
mnarst, mosoar, mxsoar, n1soar, mnsoar, mnarst, mosoar, mxsoar, n1soar, mnsoar,
moartr, mxartr, n1artr, mnartr, moartr, mxartr, n1artr, mnartr,
mxarcf, mnarcf2, mnarcf3, mxarcf, mnarcf2, mnarcf3,
mn1arcf, mnarcf, mnarcf1, mn1arcf, mnarcf, mnarcf1,
comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig, comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
ierr ); ierr );
if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;} if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;} if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;} if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;} if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
deltacpu_( d ); deltacpu_( d );
tcpu += 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; if( ierr != 0 ) goto ERREUR;
//qualites de la triangulation finale //qualites de la triangulation finale
@ -626,11 +579,7 @@ void aptrte( Z nutysu, R aretmx,
if( mnartr[nt*moartr-moartr] != 0 ) if( mnartr[nt*moartr-moartr] != 0 )
{ {
//le numero des 3 sommets du triangle nt //le numero des 3 sommets du triangle nt
#ifdef DFORTRAN nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
NUSOTR( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
#else
nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
#endif
//les 3 sommets du triangle sont actifs //les 3 sommets du triangle sont actifs
mnarst[ nosotr[0] ] = 1; mnarst[ nosotr[0] ] = 1;
mnarst[ nosotr[1] ] = 1; mnarst[ nosotr[1] ] = 1;
@ -689,7 +638,7 @@ void aptrte( Z nutysu, R aretmx,
// ----------------------------------------------------- // -----------------------------------------------------
// boucle sur les triangles occupes (internes et externes) // boucle sur les triangles occupes (internes et externes)
if( nust != NULL ) delete [] nust; if( nust != NULL ) delete [] nust;
nust = new Z[4*nbt]; nust = new Z[nbsttria*nbt];
if( nust == NULL ) goto ERREUR; if( nust == NULL ) goto ERREUR;
nbt = 0; nbt = 0;
for (i=1; i<=mxartr; i++) for (i=1; i<=mxartr; i++)
@ -698,24 +647,19 @@ void aptrte( Z nutysu, R aretmx,
if( mnartr[i*moartr-moartr] != 0 ) if( mnartr[i*moartr-moartr] != 0 )
{ {
//le triangle i est interne => nosotr numero de ses 3 sommets //le triangle i est interne => nosotr numero de ses 3 sommets
#ifdef DFORTRAN nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
NUSOTR( i, mosoar, mnsoar, moartr, mnartr, nosotr );
#else
nusotr_( i, mosoar, mnsoar, moartr, mnartr, nosotr );
#endif
nust[nbt++] = mnarst[ nosotr[0] ]; nust[nbt++] = mnarst[ nosotr[0] ];
nust[nbt++] = mnarst[ nosotr[1] ]; nust[nbt++] = mnarst[ nosotr[1] ];
nust[nbt++] = mnarst[ nosotr[2] ]; nust[nbt++] = mnarst[ nosotr[2] ];
nust[nbt++] = 0; nust[nbt++] = 0;
} }
} }
nbt /= 4; //le nombre final de triangles de la surface nbt /= nbsttria; //le nombre final de triangles de la surface
// MESSAGEE("Nombre de sommets=" << nbst MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
// << " Nombre de triangles=" << nbt); << nbt << " triangles=" << nbt);
deltacpu_( d ); deltacpu_( d );
tcpu += 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 // destruction des tableaux auxiliaires
// ------------------------------------ // ------------------------------------
@ -725,7 +669,6 @@ void aptrte( Z nutysu, R aretmx,
if( mnslig != NULL ) delete [] mnslig; if( mnslig != NULL ) delete [] mnslig;
if( mnsoar != NULL ) delete [] mnsoar; if( mnsoar != NULL ) delete [] mnsoar;
if( mnpxyd != NULL ) delete [] mnpxyd; if( mnpxyd != NULL ) delete [] mnpxyd;
if( mndalf != NULL ) delete [] mndalf;
if( mntree != NULL ) delete [] mntree; if( mntree != NULL ) delete [] mntree;
if( mnqueu != NULL ) delete [] mnqueu; if( mnqueu != NULL ) delete [] mnqueu;
if( mntrsu != NULL ) delete [] mntrsu; if( mntrsu != NULL ) delete [] mntrsu;
@ -747,7 +690,7 @@ void aptrte( Z nutysu, R aretmx,
} }
else else
{ {
// MESSAGEE( "Triangulation non realisee " << ierr ); MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
if( ierr == 0 ) ierr=1; if( ierr == 0 ) ierr=1;
goto NETTOYAGE; goto NETTOYAGE;
} }
@ -806,18 +749,10 @@ void qualitetrte( R3 *mnpxyd,
nbtria++; nbtria++;
//le numero des 3 sommets du triangle nt //le numero des 3 sommets du triangle nt
#ifdef DFORTRAN nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
NUSOTR( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
#else
nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
#endif
//la qualite du triangle ns1 ns2 ns3 //la qualite du triangle ns1 ns2 ns3
#ifdef DFORTRAN qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
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
qualite ); qualite );
//la qualite moyenne //la qualite moyenne
@ -827,19 +762,14 @@ void qualitetrte( R3 *mnpxyd,
quamin = Min( quamin, qualite ); quamin = Min( quamin, qualite );
//aire signee du triangle nt //aire signee du triangle nt
#ifdef DFORTRAN d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
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
if( d<0 ) if( d<0 )
{ {
//un triangle d'aire negative de plus //un triangle d'aire negative de plus
nbtrianeg++; nbtrianeg++;
// MESSAGEE("ATTENTION: le triangle " << nt << " de sommets:" MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
// << nosotr[0] << " " << nosotr[1] << " " << nosotr[2] << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
// << " a une aire " << d <<"<=0"); << " a une aire " << d <<"<=0");
} }
//aire des triangles actuels //aire des triangles actuels
@ -849,12 +779,12 @@ void qualitetrte( R3 *mnpxyd,
//les affichages //les affichages
quamoy /= nbtria; quamoy /= nbtria;
// MESSAGEE("Qualite moyenne=" << quamoy MESSAGE("Qualite moyenne=" << quamoy
// << " Qualite minimale=" << quamin << " Qualite minimale=" << quamin
// << " des " << nbtria << " triangles de surface totale=" << " des " << nbtria << " triangles de surface plane totale="
// << aire); << aire);
//if( nbtrianeg>0 ) if( nbtrianeg>0 )
// MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg ); MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );
return; return;
} }

189
src/MEFISTO2/aptrte.h
View File

@ -141,15 +141,55 @@ MEFISTO2D_EXPORT
// auteur : Alain Perronnet Analyse Numerique Paris UPMC decembre 2001 // 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 extern "C" { void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN tempscpu( double & tempsec );
TEMPSCPU( double & tempsec ); } }
#else
tempscpu_( double & tempsec ); }
#endif
//Retourne le temps CPU utilise en secondes //Retourne le temps CPU utilise en secondes
@ -157,11 +197,8 @@ extern "C" { void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN deltacpu( R & dtcpu );
DELTACPU( R & dtcpu ); } }
#else
deltacpu_( R & dtcpu ); }
#endif
//Retourne le temps CPU utilise en secondes depuis le precedent appel //Retourne le temps CPU utilise en secondes depuis le precedent appel
@ -170,34 +207,25 @@ extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN insoar( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );
INSOAR( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );} }
#else
insoar_( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );}
#endif
//mettre a zero les nb entiers de tab //mettre a zero les nb entiers de tab
extern "C" {void extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN azeroi( Z & nb, Z * tab );
AZEROI( Z & nb, Z * tab );} }
#else
azeroi_( Z & nb, Z * tab );}
#endif
extern "C" {void extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN fasoar( Z & ns1, Z & ns2, Z & nt1, Z & nt2, Z & nolign,
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 & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z * mnarst,
Z & noar, Z & ierr );} Z & noar, Z & ierr );
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// but : former l'arete de sommet ns1-ns2 dans le hachage du tableau // but : former l'arete de sommet ns1-ns2 dans le hachage du tableau
// ----- nosoar des aretes de la triangulation // ----- nosoar des aretes de la triangulation
@ -246,29 +274,20 @@ extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN teajte( Z & mxsomm, Z & nbsomm, R3 * mnpxyd, R3 * comxmi,
TEAJTE
#else
teajte_
#endif
( Z & mxsomm, Z & nbsomm, R3 * mnpxyd, R3 * comxmi,
R & aretmx, Z & mxtree, Z * letree, R & aretmx, Z & mxtree, Z * letree,
Z & ierr );} Z & ierr );
}
extern "C" {void extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN tehote( Z & nutysu, Z & nbarpi, Z & mxsomm, Z & nbsomm, R3 * mnpxyd,
TEHOTE
#else
tehote_
#endif
( Z & nutysu, Z & nbarpi, Z & mxsomm, Z & nbsomm, R3 * mnpxyd,
R3 * comxmi, R & aretmx, R3 * comxmi, R & aretmx,
Z * letree, Z & mxqueu, Z * mnqueu, Z * letree, Z & mxqueu, Z * mnqueu,
Z & ierr );} Z & ierr );
}
// homogeneisation de l'arbre des te a un saut de taille au plus // 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 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
@ -276,16 +295,12 @@ extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN tetrte( R3 * comxmi, R & aretmx, Z & nbarpi, Z & mxsomm, R3 * mnpxyd,
TETRTE
#else
tetrte_
#endif
( R3 * comxmi, R & aretmx, Z & nbarpi, Z & mxsomm, R3 * mnpxyd,
Z & mxqueu, Z * mnqueu, Z * mntree, Z & mxqueu, Z * mnqueu, Z * mntree,
Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst, Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
Z & ierr );} Z & ierr );
}
// trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
// et des points de la frontiere, des points internes imposes interieurs // et des points de la frontiere, des points internes imposes interieurs
@ -293,44 +308,32 @@ extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN aisoar( Z & mosoar, Z & mxsoar, Z * mnsoar, Z & na );
AISOAR }
#else // formation du chainage 6 des aretes internes a echanger eventuellement
aisoar_
#endif
( Z & mosoar, Z & mxsoar, Z * mnsoar, Z & na );}
// formation du chainage 6 des aretes internes a echanger eventuellement
extern "C" {void extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN tedela( R3 * mnpxyd, Z * mnarst,
TEDELA
#else
tedela_
#endif
( R3 * mnpxyd, Z * mnarst,
Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & na, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & na,
Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z & n );} 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 // 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 extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN terefr( Z & nbarpi, R3 * mnpxyd,
TEREFR
#else
terefr_
#endif
( Z & nbarpi, R3 * mnpxyd,
Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
Z & moartr, Z & n1artr, Z * mnartr, Z * mnarst, Z & moartr, Z & n1artr, Z * mnartr, Z * mnarst,
Z & mxarcf, Z * mnarc1, Z * mnarc2, Z & mxarcf, Z * mnarc1, Z * mnarc2,
Z * mnarc3, Z * mnarc4, Z * mnarc3, Z * mnarc4,
Z & n, Z & ierr );} Z & n, Z & ierr );
}
// detection des aretes frontalieres initiales perdues // detection des aretes frontalieres initiales perdues
// triangulation frontale pour les restaurer // triangulation frontale pour les restaurer
@ -338,35 +341,27 @@ extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN tesuex( Z & nblf, Z * nulftr,
TESUEX
#else
tesuex_
#endif
( Z & nblf, Z * nulftr,
Z & ndtri0, Z & nbsomm, R3 * mnpxyd, Z * mnslig, Z & ndtri0, Z & nbsomm, R3 * mnpxyd, Z * mnslig,
Z & mosoar, Z & mxsoar, Z * mnsoar, Z & mosoar, Z & mxsoar, Z * mnsoar,
Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst, Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
Z & nbtria, Z * mntrsu, Z & ierr );} Z & nbtria, Z * mntrsu, Z & ierr );
}
// suppression des triangles externes a la surface // suppression des triangles externes a la surface
extern "C" {void extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN teamqt( Z & nutysu,
TEAMQT
#else
teamqt_
#endif
( Z & nutysu,
Z * mnarst, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z * mnarst, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr,
Z & mxarcf, Z * mntrcf, Z * mnstbo, Z & mxarcf, Z * mntrcf, Z * mnstbo,
Z * n1arcf, Z * mnarcf, Z * mnarc1, Z * n1arcf, Z * mnarcf, Z * mnarc1,
R3 * comxmi, Z & nbarpi, Z & nbsomm, Z & mxsomm, R3 * comxmi, Z & nbarpi, Z & nbsomm, Z & mxsomm,
R3 * mnpxyd, Z * mnslig, R3 * mnpxyd, Z * mnslig,
Z & ierr );} Z & ierr );
}
// amelioration de la qualite de la triangulation par // amelioration de la qualite de la triangulation par
// barycentrage des sommets internes a la triangulation // barycentrage des sommets internes a la triangulation
// suppression des aretes trop longues ou trop courtes // suppression des aretes trop longues ou trop courtes
@ -377,36 +372,24 @@ extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN nusotr( Z & nt, Z & mosoar, Z * mnsoar, Z & moartr, Z * mnartr,Z * nosotr );
NUSOTR }
#else
nusotr_
#endif
( Z & nt, Z & mosoar, Z * mnsoar, Z & moartr, Z * mnartr,Z * nosotr );}
//retrouver les numero des 3 sommets du triangle nt //retrouver les numero des 3 sommets du triangle nt
extern "C" {void extern "C" {void
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN qutr2d( R3 & p1, R3 & p2, R3 & p3, R & qualite );
QUTR2D }
#else
qutr2d_
#endif
( R3 & p1, R3 & p2, R3 & p3, R & qualite );}
//calculer la qualite d'un triangle de R2 de sommets p1, p2, p3 //calculer la qualite d'un triangle de R2 de sommets p1, p2, p3
extern "C" { R extern "C" { R
#ifdef WIN32 #ifdef WIN32
__stdcall __stdcall
#endif #endif
#ifdef DFORTRAN surtd2( R3 & p1, R3 & p2, R3 & p3 );
SURTD2 }
#else
surtd2_
#endif
( R3 & p1, R3 & p2, R3 & p3 ); }
//calcul de la surface d'un triangle defini par 3 points de r**2 //calcul de la surface d'un triangle defini par 3 points de r**2
#endif #endif

225
src/MEFISTO2/trte.f
View File

@ -1,6 +1,6 @@
c MEFISTO : library to compute 2D triangulation from segmented boundaries c MEFISTO2: a library to compute 2D triangulation from segmented boundaries
c c
c Copyright (C) 2003 Laboratoire J.-L. Lions UPMC Paris c Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris
c c
c This library is free software; you can redistribute it and/or 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 modify it under the terms of the GNU Lesser General Public
@ -16,12 +16,12 @@ c You should have received a copy of the GNU Lesser General Public
c License along with this library; if not, write to the Free Software 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 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
c 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 c File : trte.f le Fortran du trianguleur plan
c File : trte.f
c Module : SMESH c Module : SMESH
c Author: Alain PERRONNET c Author : Alain PERRONNET
c Date : 16 mars 2006
subroutine qutr2d( p1, p2, p3, qualite ) subroutine qutr2d( p1, p2, p3, qualite )
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -1257,7 +1257,8 @@ c la direction pour le calcul de la longueur (inactif ici!)
xyzd(2) = 0d0 xyzd(2) = 0d0
xyzd(3) = 0d0 xyzd(3) = 0d0
longai = areteideale(xyz,xyzd) longai = areteideale()
c (xyz,xyzd)
if( longai .lt. 0d0 ) then if( longai .lt. 0d0 ) then
write(imprim,10000) xyz write(imprim,10000) xyz
10000 format('attention: longueur de areteideale(', 10000 format('attention: longueur de areteideale(',
@ -2640,6 +2641,8 @@ c....................................................................012
% noarst(*) % noarst(*)
integer lapile(1:mxpile) integer lapile(1:mxpile)
integer nosotr(3) integer nosotr(3)
c
lhpile = 0
c c
c la premiere arete de sommet ns c la premiere arete de sommet ns
nar = noarst( ns ) nar = noarst( ns )
@ -2817,13 +2820,13 @@ c les triangles de sommet ns ne forment pas une boule de centre ns
c c
c saturation de la pile des triangles c saturation de la pile des triangles
c ----------------------------------- c -----------------------------------
9990 write(imprim,*) 'trp1st:saturation pile des triangles autour ', 9990 write(imprim,*)'trp1st: saturation pile des triangles autour ',
%'sommet',ns %'sommet',ns
goto 9999 goto 9999
c c
c erreur triangle ne contenant pas le sommet ns c erreur triangle ne contenant pas le sommet ns
c ---------------------------------------------- c ----------------------------------------------
9995 write(imprim,*) 'trp1st:triangle ',nta,' st=', 9995 write(imprim,*) 'trp1st: triangle ',nta,' st=',
% (nosotr(nar),nar=1,3),' sans le sommet' ,ns % (nosotr(nar),nar=1,3),' sans le sommet' ,ns
c c
9999 lhpile = 0 9999 lhpile = 0
@ -2883,7 +2886,7 @@ c le sommet nosotr(3 du triangle 123
% mosoar, mxsoar, n1soar, nosoar, % mosoar, mxsoar, n1soar, nosoar,
% moartr, mxartr, n1artr, noartr, % moartr, mxartr, n1artr, noartr,
% mxarcf, n1arcf, noarcf, larmin, notrcf, liarcf, % mxarcf, n1arcf, noarcf, larmin, notrcf, liarcf,
% nbstsu, ierr ) % ierr )
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c but : supprimer de la triangulation les sommets de te trop proches c but : supprimer de la triangulation les sommets de te trop proches
c ----- soit d'un sommet frontalier ou point interne impose c ----- soit d'un sommet frontalier ou point interne impose
@ -2931,7 +2934,6 @@ c liarcf : tableau ( mxarcf ) auxiliaire d'entiers
c c
c sortie : c sortie :
c -------- c --------
c nbstsu : nombre de sommets de te supprimes
c ierr : =0 si pas d'erreur c ierr : =0 si pas d'erreur
c >0 si une erreur est survenue c >0 si une erreur est survenue
c 11 algorithme defaillant c 11 algorithme defaillant
@ -2958,8 +2960,8 @@ c
equivalence (nosotr(1),ns1), (nosotr(2),ns2), equivalence (nosotr(1),ns1), (nosotr(2),ns2),
% (nosotr(3),ns3) % (nosotr(3),ns3)
c c
c le nombre de sommets de te supprimes cccc le nombre de sommets de te supprimes
nbstsu = 0 ccc nbstsu = 0
c c
c initialisation du chainage des aretes des cf => 0 arete de cf c initialisation du chainage des aretes des cf => 0 arete de cf
do 10 narete=1,mxsoar 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 c ils doivent former un contour ferme de type etoile
call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr, call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr,
% mxarcf, nbtrcf, notrcf ) % 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 erreur: impossible de trouver tous les triangles de sommet ns
c seule une partie est a priori retrouvee c seule une partie est a priori retrouvee
nbtrcf = -nbtrcf nbtrcf = -nbtrcf
@ -3050,8 +3053,9 @@ c ==========================================
% mxarcf, n1arcf, noarcf, % mxarcf, n1arcf, noarcf,
% larmin, notrcf, liarcf, ierr ) % larmin, notrcf, liarcf, ierr )
if( ierr .eq. 0 ) then if( ierr .eq. 0 ) then
c un sommet de te supprime de plus cccc un sommet de te supprime de plus
nbstsu = nbstsu + 1 ccc nbstsu = nbstsu + 1
goto 100
else if( ierr .lt. 0 ) then else if( ierr .lt. 0 ) then
c le sommet nste est externe donc non supprime c le sommet nste est externe donc non supprime
c ou bien le sommet nste est le centre d'un cf dont toutes 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 return
endif endif
c c
c boucle jusqu'a obtenir une qualite suffisante cccc boucle jusqu'a obtenir une qualite suffisante
c si triangulation tres irreguliere => cccc si triangulation tres irreguliere =>
c destruction de beaucoup de points internes cccc destruction de beaucoup de points internes
c les 2 variables suivantes brident ces destructions massives cccc les 2 variables suivantes brident ces destructions massives
ccc nbsuns = nbsuns + 1 ccc nbsuns = nbsuns + 1
quaopt = quaopt * 0.8 ccc quaopt = quaopt * 0.8
ccc if( nbsuns .le. 5 ) goto 15 ccc if( nbsuns .lt. 5 ) goto 15
goto 15
endif endif
endif endif
c c
100 continue 100 continue
c
c write(imprim,*)'retrait de',nbstsu,
c % ' sommets de te trop proches de la frontiere'
return
end end
@ -3370,7 +3377,10 @@ c l'arete suivante
% mxtrcf, n1arcf, noarcf, % mxtrcf, n1arcf, noarcf,
% larmin, notrcf, nostbo, % larmin, notrcf, nostbo,
% ierr ) % 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 le sommet ns est externe donc non supprime
c ou bien le sommet ns est le centre d'un cf dont toutes c ou bien le sommet ns est le centre d'un cf dont toutes
c les aretes simples sont frontalieres 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 noarcf : tableau (3,mxtrcf) auxiliaire d'entiers
c larmin : tableau ( mxtrcf ) auxiliaire d'entiers c larmin : tableau ( mxtrcf ) auxiliaire d'entiers
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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 c....................................................................012
parameter (lchain=6) parameter (lchain=6)
common / unites / lecteu, imprim, nunite(30) common / unites / lecteu, imprim, nunite(30)
double precision pxyd(3,*) double precision pxyd(3,*)
double precision ponder, ponde1, xbar, ybar, x, y, d, dmin, dmax double precision ponder, ponde1, xbar, ybar, x, y, d, dmin, dmax
double precision surtd2
double precision d2d3(3,3) double precision d2d3(3,3)
real origin(3), xyz(3) real origin(3), xyz(3)
integer noartr(moartr,*), integer noartr(moartr,*),
@ -3745,7 +3756,8 @@ c les compteurs de passage sur les differents cas
nbst8 = 0 nbst8 = 0
c c
c coefficient de ponderation croissant avec les iterations 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 ponde1 = 1d0 - ponder
c c
c l'ordre du parcours dans le sens croissant ou decroissant 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 c ils doivent former un contour ferme de type etoile
call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr, call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr,
% mxtrcf, nbtrcf, notrcf ) % mxtrcf, nbtrcf, notrcf )
if( nbtrcf .le. 0 ) goto 1000 if( nbtrcf .le. 2 ) goto 1000
c c
c boucle sur les triangles qui forment une boule autour du sommet ns c boucle sur les triangles qui forment une boule autour du sommet ns
nbstbo = 0 nbstbo = 0
@ -3840,9 +3852,14 @@ c pas de modification de la topologie lors de la derniere iteration
c ================================================================= c =================================================================
if( iter .ge. nbitaq ) goto 200 if( iter .ge. nbitaq ) goto 200
c c
c si la boule de ns contient 3 ou 4 triangles le sommet ns est detruit c si la boule de ns contient au plus 3 triangles
c ==================================================================== c => pas de changement de topologie
if( nbtrcf .le. 4 ) then 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
c remise a -1 du chainage des aretes peripheriques de la boule ns c remise a -1 du chainage des aretes peripheriques de la boule ns
noar = noar0 noar = noar0
@ -3861,7 +3878,10 @@ c l'arete suivante
% mxtrcf, n1arcf, noarcf, % mxtrcf, n1arcf, noarcf,
% larmin, notrcf, nostbo, % larmin, notrcf, nostbo,
% ierr ) % 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 le sommet ns est externe donc non supprime
c ou bien le sommet ns est le centre d'un cf dont toutes c ou bien le sommet ns est le centre d'un cf dont toutes
c les aretes simples sont frontalieres 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 nbstsu = nbstsu + 1
else else
c erreur irrecuperable c erreur irrecuperable
write(imprim,*)
% 'teamqs: erreur1 irrecuperable en sortie te1stm'
goto 9999 goto 9999
endif endif
goto 1000 goto 1000
@ -3910,6 +3932,9 @@ c l'arete suivante
endif endif
c c
c le point ns1 devient le milieu de l'arete ns-ns1 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 do 75 j=1,3
pxyd(j,ns1) = (pxyd(j,ns) + pxyd(j,ns1)) * 0.5d0 pxyd(j,ns1) = (pxyd(j,ns) + pxyd(j,ns1)) * 0.5d0
75 continue 75 continue
@ -3931,16 +3956,22 @@ c suppression du point ns et mise en delaunay
if( ierr .lt. 0 ) then if( ierr .lt. 0 ) then
c le sommet ns est externe donc non supprime c le sommet ns est externe donc non supprime
c ou bien le sommet ns est le centre d'un cf dont toutes c ou bien le sommet ns est le centre d'un cf dont toutes
c les aretes simples sont frontalieres c les aretes simples sont frontalieres ou erreur
c dans les 2 cas le sommet ns n'est pas supprime 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 ierr = 0
goto 200 goto 1000
else if( ierr .eq. 0 ) then else if( ierr .eq. 0 ) then
nbstsu = nbstsu + 1 nbstsu = nbstsu + 1
nbst5 = nbst5 + 1 nbst5 = nbst5 + 1
goto 1000 goto 1000
else else
c erreur irrecuperable c erreur irrecuperable
write(imprim,*)
% 'teamqs: erreur2 irrecuperable en sortie te1stm'
goto 9999 goto 9999
endif endif
endif endif
@ -3983,9 +4014,9 @@ c le numero pxyd de ses 3 sommets
c c
c ajout du nouveau barycentre c ajout du nouveau barycentre
if( nbsomm .ge. mxsomm ) then if( nbsomm .ge. mxsomm ) then
write(imprim,*) 'saturation du tableau pxyd' write(imprim,*) 'teamqs: saturation du tableau pxyd'
c abandon de l'amelioration c abandon de l'amelioration
goto 1100 goto 9999
endif endif
nbsomm = nbsomm + 1 nbsomm = nbsomm + 1
do 120 i=1,3 do 120 i=1,3
@ -4024,7 +4055,11 @@ c protection a ne pas modifier sinon erreur!
% moartr, mxartr, n1artr, noartr, % moartr, mxartr, n1artr, noartr,
% noarst, % noarst,
% nosotr, ierr ) % 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 140 continue
c c
nbst8 = nbst8 + 1 nbst8 = nbst8 + 1
@ -4041,10 +4076,56 @@ c simples de la boule du sommet ns
200 xbar = xbar / nbstbo 200 xbar = xbar / nbstbo
ybar = ybar / nbstbo ybar = ybar / nbstbo
c 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(1,ns) = ponde1 * pxyd(1,ns) + ponder * xbar
pxyd(2,ns) = ponde1 * pxyd(2,ns) + ponder * ybar pxyd(2,ns) = ponde1 * pxyd(2,ns) + ponder * ybar
c 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 c les aretes chainees de la boule sont rendues delaunay
300 call tedela( pxyd, noarst, 300 call tedela( pxyd, noarst,
% mosoar, mxsoar, n1soar, nosoar, noar0, % mosoar, mxsoar, n1soar, nosoar, noar0,
@ -4052,11 +4133,9 @@ c les aretes chainees de la boule sont rendues delaunay
c c
1000 continue 1000 continue
c c
c trace de la triangulation actuelle et calcul de la qualite ccc write(imprim,11000) iter, nbitaq, nbst4, nbst5, nbst8
1100 continue ccc11000 format( 'teamqs iter=',i2,' max iter=',i2,':',
c ccc % i7,' sommets de 4t',
ccc write(imprim,11000) nbst4, nbst5, nbst8
ccc11000 format( i7,' sommets de 4t',
ccc % i7,' sommets 5t+5t', ccc % i7,' sommets 5t+5t',
ccc % i7,' sommets >7t' ) ccc % i7,' sommets >7t' )
c c
@ -4083,7 +4162,7 @@ c
% comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign, % comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign,
% ierr ) % ierr )
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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 c
c entrees: c entrees:
@ -4153,10 +4232,8 @@ c ==============================================================
% mosoar, mxsoar, n1soar, nosoar, % mosoar, mxsoar, n1soar, nosoar,
% moartr, mxartr, n1artr, noartr, % moartr, mxartr, n1artr, noartr,
% mxarcf, n1arcf, noarcf, larmin, notrcf, nostbo, % mxarcf, n1arcf, noarcf, larmin, notrcf, nostbo,
% nbstsu, ierr ) % ierr )
if( ierr .ne. 0 ) goto 9999 if( ierr .ne. 0 ) goto 9999
c write(imprim,*) 'retrait de',nbstsu,
c % ' sommets de te trop proches de la frontiere'
c c
c ajustage des tailles moyennes des aretes avec ampli=1.34d0 entre c ajustage des tailles moyennes des aretes avec ampli=1.34d0 entre
c ampli/2 x taille_souhaitee et ampli x taille_souhaitee c ampli/2 x taille_souhaitee et ampli x taille_souhaitee
@ -4171,18 +4248,18 @@ c ================================================================
% ierr ) % ierr )
if( ierr .ne. 0 ) goto 9999 if( ierr .ne. 0 ) goto 9999
c c
c modification de la topologie autour des sommets frontaliers cccc modification de la topologie autour des sommets frontaliers
c pour avoir un nombre de triangles egal a l'angle/60 degres cccc pour avoir un nombre de triangles egal a l'angle/60 degres
c et mise en triangulation delaunay locale cccc et mise en triangulation delaunay locale
c =========================================================== cccc ===========================================================
call teamsf( nutysu, ccc call teamsf( nutysu,
% noarst, mosoar, mxsoar, n1soar, nosoar, ccc % noarst, mosoar, mxsoar, n1soar, nosoar,
% moartr, mxartr, n1artr, noartr, ccc % moartr, mxartr, n1artr, noartr,
% mxarcf, notrcf, nostbo, ccc % mxarcf, notrcf, nostbo,
% n1arcf, noarcf, larmin, ccc % n1arcf, noarcf, larmin,
% comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign, ccc % comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign,
% ierr ) ccc % ierr )
if( ierr .ne. 0 ) goto 9999 ccc if( ierr .ne. 0 ) goto 9999
c c
c quelques iterations de barycentrage des points internes c quelques iterations de barycentrage des points internes
c modification de la topologie pour avoir 4 ou 5 ou 6 triangles 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 >0 si une erreur est survenue
c =11 algorithme defaillant c =11 algorithme defaillant
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c auteur : alain perronnet analyse numerique paris upmc mars 1997 c auteur : alain perronnet analyse numerique paris upmc mars 2006
c....................................................................012 c....................................................................012
parameter ( lchain=6, quamal=0.3) parameter ( lchain=6, quamal=0.3)
common / unites / lecteu,imprim,intera,nunite(29) 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, % moartr, n1artr, noartr,
% nbarcf, n1arcf, noarcf, % nbarcf, n1arcf, noarcf,
% ierr ) % 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
c ici le sommet nsasup appartient a aucune arete c ici le sommet nsasup appartient a aucune arete
noarst( nsasup ) = 0 noarst( nsasup ) = 0
@ -7247,11 +7328,12 @@ c n1aeoc : numero dans nosoar de la premiere arete simple de l'etoile
c c
c sortie : c sortie :
c -------- 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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c auteur : alain perronnet analyse numerique paris upmc mars 1997 c auteur : alain perronnet analyse numerique paris upmc mars 1997
c2345x7..............................................................012 c2345x7..............................................................012
parameter (lchain=6) parameter (lchain=6)
common / unites / lecteu, imprim, nunite(30)
integer nosoar(mosoar,mxsoar) integer nosoar(mosoar,mxsoar)
c c
c si l'arete n'appartient pas aux aretes de l'etoile naetoi 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 c arete double de l'etoile. elle est supprimee du chainage
na0 = 0 na0 = 0
na = n1aeoc na = n1aeoc
nbpass = 0
c parcours des aretes chainees jusqu'a trouver l'arete noar c parcours des aretes chainees jusqu'a trouver l'arete noar
10 if( na .ne. noar ) then 10 if( na .ne. noar ) then
c passage a la suivante c passage a la suivante
na0 = na na0 = na
na = nosoar( lchain, 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 goto 10
endif endif
c 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 15 si une seule arete simple frontaliere
c 16 si boucle infinie car toutes les aretes simples c 16 si boucle infinie car toutes les aretes simples
c de la boule sont frontalieres! c de la boule sont frontalieres!
c 17 si boucle infinie dans caetoi
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c auteur : alain perronnet analyse numerique paris upmc mars 1997 c auteur : alain perronnet analyse numerique paris upmc mars 1997
c....................................................................012 c....................................................................012
@ -7384,7 +7478,12 @@ c l'arete de nosoar a traiter
noar = abs( noartr(j,nt) ) noar = abs( noartr(j,nt) )
call caetoi( noar, mosoar, mxsoar, n1soar, nosoar, call caetoi( noar, mosoar, mxsoar, n1soar, nosoar,
% n1aeoc, nbtrar ) % 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( nbtrar .eq. 1 ) then
if( nosoar(4,noar) .eq. nt ) then if( nosoar(4,noar) .eq. nt ) then
nosoar(4,noar) = nosoar(5,noar) nosoar(4,noar) = nosoar(5,noar)