MPV: new version of MEFISTO integration

This commit is contained in:
mpv 2007-01-10 13:38:28 +00:00
parent fe284ca490
commit 8595d10d2f
6 changed files with 1297 additions and 1251 deletions

View File

@ -1,6 +1,6 @@
// MEFISTO : library to compute 2D triangulation from segmented boundaries // MEFISTO : 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
@ -23,7 +23,7 @@
// File : Rn.h // File : Rn.h
// Module : SMESH // Module : SMESH
// Authors: Frederic HECHT & Alain PERRONNET // Authors: Frederic HECHT & Alain PERRONNET
// // Date : 13 novembre 2006
#ifndef Rn__h #ifndef Rn__h
#define Rn__h #define Rn__h
@ -66,15 +66,7 @@ typedef unsigned long int N;
//le type Z des nombres entiers relatifs //le type Z des nombres entiers relatifs
//========= //=========
// 64-bit porting: "long" replaced with "int". typedef long int Z;
// On 64-bit, C++ long type is 8 byte long. MEFISTO2D C code calls several Fortran subroutines passing
// arguments of this type, however Fortran knows nothing about changed size of arguments,
// therefore stack gets corrupted. With "int" used instead of "long", Fortran calls from C do no harm to the stack
// After this modification, behavior on 32-bit platforms does not change: on all platforms supported by
// SALOME 3, "int" and "long" have the same size of 4 bytes.
//========
//typedef long int Z;
typedef int Z;
//le type R des nombres "reels" //le type R des nombres "reels"
//========= //=========

View File

@ -1,5 +1,6 @@
// MEFISTO2: a library to compute 2D triangulation from segmented boundaries // MEFISTO2: a library to compute 2D triangulation from segmented boundaries
// //
//
// Copyright (C) 2006 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
@ -21,7 +22,7 @@
// File : aptrte.cxx le C++ de l'appel du trianguleur plan // File : aptrte.cxx le C++ de l'appel du trianguleur plan
// Module : SMESH // Module : SMESH
// Author : Alain PERRONNET // Author : Alain PERRONNET
// Date : 16 mars 2006 // Date : 13 novembre 2006
#include "Rn.h" #include "Rn.h"
#include "aptrte.h" #include "aptrte.h"
@ -139,15 +140,14 @@ void aptrte( Z nutysu, R aretmx,
R3 comxmi[2]; //coordonnees UV Min et Maximales R3 comxmi[2]; //coordonnees UV Min et Maximales
R aremin, aremax; //longueur minimale et maximale des aretes R aremin, aremax; //longueur minimale et maximale des aretes
R airemx; //aire maximale souhaitee d'un triangle
R quamoy, quamin; R quamoy, quamin;
Z noar0, noar, na; Z noar0, noar, na;
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; Z nuds = 0;
aretemaxface_ = aretmx;
// initialisation du temps cpu // initialisation du temps cpu
deltacpu_( d ); deltacpu_( d );
@ -161,7 +161,8 @@ void aptrte( Z nutysu, R aretmx,
i = 4*nbarfr/10; i = 4*nbarfr/10;
mxsomm = Max( 20000, 64*nbpti+i*i ); mxsomm = Max( 20000, 64*nbpti+i*i );
MESSAGE( "APTRTE: Debut de la triangulation plane avec " ); MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx << " mxsomm=" << mxsomm ); MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx
<< " mxsomm=" << mxsomm );
MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes"); MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
NEWDEPART: NEWDEPART:
@ -272,6 +273,9 @@ void aptrte( Z nutysu, R aretmx,
//l'arete precedente est dotee de sa suivante:celle cree ensuite //l'arete precedente est dotee de sa suivante:celle cree ensuite
//les 2 coordonnees du sommet ns2 de la ligne //les 2 coordonnees du sommet ns2 de la ligne
ns = ns1 - 1; ns = ns1 - 1;
//debut ajout 5/10/2006 ................................................
nuds = Max( nuds, ns ); //le numero du dernier sommet traite
//fin ajout 5/10/2006 ................................................
mnpxyd[ns].x = uvslf[ns].x; mnpxyd[ns].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 );
@ -284,6 +288,14 @@ void aptrte( Z nutysu, R aretmx,
aremin = Min( aremin, d ); aremin = Min( aremin, d );
aremax = Max( aremax, d ); aremax = Max( aremax, d );
//debut ajout du 5/10/2006 .............................................
//la longueur de l'arete ns1-ns2
d = sqrt( d );
//longueur arete = Min ( aretmx, aretes incidentes )
mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d );
mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
//fin ajout du 5/10/2006 ...............................................
//le numero n de la ligne du sommet et son numero ns1 dans la ligne //le numero n de la ligne du sommet et son numero ns1 dans la ligne
mnslig[ns] = 1000000 * n + ns1-nudslf[n-1]; mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
@ -307,9 +319,34 @@ void aptrte( Z nutysu, R aretmx,
aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
aremax = sqrt( aremax ); //longueur maximale d'une arete aremax = sqrt( aremax ); //longueur maximale d'une arete
aretmx = Min( aretmx, aremax ); //pour homogeneiser //debut ajout 9/11/2006 ................................................
MESSAGE("nutysu=" << nutysu << " aretmx=" << aretmx // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
<< " arete min=" << aremin << " arete max=" << aremax);
// protection contre une arete max desiree trop grande ou trop petite
if( aretmx > aremax*2.05 ) aretmx = aremax;
// protection contre une arete max desiree trop petite
if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
aretmx =(aremin+aremax*2)/3.0;
if( aretmx < aremin && aremin > 0 )
aretmx = aremin;
//sauvegarde pour la fonction areteideale_
aretemaxface_ = aretmx;
//aire maximale souhaitee des triangles
airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral
for(i=0; i<=nuds; i++ )
mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
//MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
//fin ajout 9/11/2006 .................................................
MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax );
MESSAGE("Triangulation: arete mx=" << aretmx
<< " triangle aire mx=" << airemx );
//chainage des aretes frontalieres : la derniere arete frontaliere //chainage des aretes frontalieres : la derniere arete frontaliere
mnsoar[ mosoar * noar - mosoar + 5 ] = 0; mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
@ -433,7 +470,7 @@ void aptrte( Z nutysu, R aretmx,
if( ierr != 0 ) goto ERREUR; if( ierr != 0 ) goto ERREUR;
//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,
nbt, quamoy, quamin ); nbt, quamoy, quamin );
// boucle sur les aretes internes (non sur une ligne de la frontiere) // boucle sur les aretes internes (non sur une ligne de la frontiere)
@ -452,7 +489,7 @@ void aptrte( Z nutysu, R aretmx,
<< 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,
nbt, quamoy, quamin ); nbt, quamoy, quamin );
// detection des aretes frontalieres initiales perdues // detection des aretes frontalieres initiales perdues
@ -474,11 +511,11 @@ void aptrte( Z nutysu, R aretmx,
terefr_( nbarpi, mnpxyd, terefr_( nbarpi, mnpxyd,
mosoar, mxsoar, n1soar, mnsoar, mosoar, mxsoar, n1soar, mnsoar,
moartr, n1artr, mnartr, mnarst, moartr, mxartr, n1artr, mnartr, mnarst,
mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2, mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
n, ierr ); n, ierr );
MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" ); MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr );
deltacpu_( d ); deltacpu_( d );
tcpu += d; tcpu += d;
MESSAGE("Temps de la recuperation des aretes perdues de la frontiere=" MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
@ -487,7 +524,7 @@ void aptrte( Z nutysu, R aretmx,
if( ierr != 0 ) goto ERREUR; if( ierr != 0 ) goto ERREUR;
//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,
nbt, quamoy, quamin ); nbt, quamoy, quamin );
// fin de la triangulation avec respect des aretes initiales frontalieres // fin de la triangulation avec respect des aretes initiales frontalieres
@ -524,11 +561,11 @@ void aptrte( Z nutysu, R aretmx,
deltacpu_( d ); deltacpu_( d );
tcpu += d; tcpu += d;
MESSAGE( "Temps de la suppression des triangles externes=" << d ); MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
if( ierr != 0 ) goto ERREUR; if( ierr != 0 ) goto ERREUR;
//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,
nbt, quamoy, quamin ); nbt, quamoy, quamin );
// amelioration de la qualite de la triangulation par // amelioration de la qualite de la triangulation par
@ -543,12 +580,12 @@ void aptrte( Z nutysu, R aretmx,
cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl; cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
goto ERREUR; goto ERREUR;
} }
teamqt_( nutysu, teamqt_( nutysu, aretmx, airemx,
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, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
ierr ); ierr );
if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;} if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;} if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
@ -559,10 +596,11 @@ void aptrte( Z nutysu, R aretmx,
deltacpu_( d ); deltacpu_( d );
tcpu += d; tcpu += d;
MESSAGE( "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 == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
if( ierr != 0 ) goto ERREUR; if( ierr != 0 ) goto ERREUR;
//qualites de la triangulation finale //qualites de la triangulation finale
qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr, qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
nbt, quamoy, quamin ); nbt, quamoy, quamin );
// renumerotation des sommets internes: mnarst(i)=numero final du sommet // renumerotation des sommets internes: mnarst(i)=numero final du sommet
@ -652,7 +690,7 @@ void aptrte( Z nutysu, R aretmx,
} }
nbt /= nbsttria; //le nombre final de triangles de la surface nbt /= nbsttria; //le nombre final de triangles de la surface
MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et " MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
<< nbt << " triangles=" << nbt); << nbt << " triangles" );
deltacpu_( d ); deltacpu_( d );
tcpu += d; tcpu += d;
MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" ); MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
@ -693,7 +731,7 @@ void aptrte( Z nutysu, R aretmx,
} }
void qualitetrte( R3 *mnpxyd, void qualitetrte_( R3 *mnpxyd,
Z & mosoar, Z & mxsoar, Z *mnsoar, Z & mosoar, Z & mxsoar, Z *mnsoar,
Z & moartr, Z & mxartr, Z *mnartr, Z & moartr, Z & mxartr, Z *mnartr,
Z & nbtria, R & quamoy, R & quamin ) Z & nbtria, R & quamoy, R & quamin )
@ -727,13 +765,14 @@ void qualitetrte( R3 *mnpxyd,
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
{ {
R d, aire, qualite; R d, aire, qualite;
Z nosotr[3], mn, nbtrianeg, nt; Z nosotr[3], mn, nbtrianeg, nt, ntqmin;
aire = 0; aire = 0;
quamoy = 0; quamoy = 0;
quamin = 2.0; quamin = 2.0;
nbtria = 0; nbtria = 0;
nbtrianeg = 0; nbtrianeg = 0;
ntqmin = 0;
mn = -moartr; mn = -moartr;
for ( nt=1; nt<=mxartr; nt++ ) for ( nt=1; nt<=mxartr; nt++ )
@ -755,7 +794,11 @@ void qualitetrte( R3 *mnpxyd,
quamoy += qualite; quamoy += qualite;
//la qualite minimale //la qualite minimale
quamin = Min( quamin, qualite ); if( qualite < quamin )
{
quamin = qualite;
ntqmin = nt;
}
//aire signee du triangle nt //aire signee du triangle nt
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] );
@ -780,7 +823,20 @@ void qualitetrte( R3 *mnpxyd,
<< " des " << nbtria << " triangles de surface plane totale=" << " des " << nbtria << " triangles de surface plane totale="
<< aire); << aire);
if( quamin<0.3 )
{
//le numero des 3 sommets du triangle ntqmin de qualite minimale
nusotr_( ntqmin, mosoar, mnsoar, moartr, mnartr, nosotr );
MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
<<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
for (int i=0;i<3;i++)
MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
<<" y="<< mnpxyd[nosotr[i]-1].y);
}
if( nbtrianeg>0 ) if( nbtrianeg>0 )
MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg ); MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );
MESSAGE(" ");
return; return;
} }

View File

@ -1,6 +1,6 @@
// SMESH MEFISTO2 : algorithm for meshing // SMESH MEFISTO2 : algorithm for meshing
// //
// 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
@ -23,6 +23,7 @@
// File : aptrte.h // File : aptrte.h
// Author : Alain PERRONNET // Author : Alain PERRONNET
// Module : SMESH // Module : SMESH
// Date : 13 novembre 2006
#ifndef aptrte__h #ifndef aptrte__h
#define aptrte__h #define aptrte__h
@ -41,10 +42,11 @@
#include <sys/types.h> #include <sys/types.h>
#include <sys/time.h> #include <sys/time.h>
void qualitetrte( R3 *mnpxyd, extern "C" {
void qualitetrte_( R3 *mnpxyd,
Z & mosoar, Z & mxsoar, Z *mnsoar, Z & mosoar, Z & mxsoar, Z *mnsoar,
Z & moartr, Z & mxartr, Z *mnartr, Z & moartr, Z & mxartr, Z *mnartr,
Z & nbtria, R & quamoy, R & quamin ); Z & nbtria, R & quamoy, R & quamin ); }
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// but : calculer la qualite moyenne et minimale de la triangulation // but : calculer la qualite moyenne et minimale de la triangulation
// ----- actuelle definie par les tableaux nosoar et noartr // ----- actuelle definie par les tableaux nosoar et noartr
@ -214,7 +216,7 @@ extern "C" {void tedela_( R3 * mnpxyd, Z * mnarst,
extern "C" {void terefr_( Z & nbarpi, R3 * mnpxyd, extern "C" {void terefr_( 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 & mxartr, 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 );}
@ -228,12 +230,12 @@ extern "C" {void tesuex_( Z & nblf, Z * nulftr,
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 teamqt_( Z & nutysu, extern "C" {void teamqt_( Z & nutysu, R & aretmx, R & airemx,
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, 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

View File

@ -1,6 +1,6 @@
c MEFISTO : library to compute 2D triangulation from segmented boundaries c MEFISTO : 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
@ -22,6 +22,7 @@ c
c File : areteideale.f c File : areteideale.f
c Module : SMESH c Module : SMESH
c Author : Alain PERRONNET c Author : Alain PERRONNET
c Date : 13 novembre 2006
double precision function areteideale( xyz, direction ) double precision function areteideale( xyz, direction )
double precision xyz(3), direction(3) double precision xyz(3), direction(3)

File diff suppressed because it is too large Load Diff

View File

@ -159,7 +159,7 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis
if (_maxElementArea > 0) if (_maxElementArea > 0)
{ {
// _edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant // _edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant
_edgeLength = 2 * sqrt(_maxElementArea/sqrt(3.0)); _edgeLength = sqrt(2. * _maxElementArea/sqrt(3.0));
isOk = true; isOk = true;
} }
else else