mshptg.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

// -*- Mode : c++ -*-
//
// -*- Mode : c++ -*-
//
// SUMMARY  : mshptg mesh generator ( form mshptg.f / inria / emc2 software) 
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curi, Paris,  FRANCE 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//
//

/*
 
 This file is part of Freefem++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */
 

// bof bof FH je ne sais pas pourquoi nnnBAMG_LONG_LONG
// nov 2005 j'ai corrige les problemes 
// il y avait un gros bug de le trie
// overflow long  car on calculais   x*x + y*y 
// pour que cela tienne dans un  long long  
//  x*x + y*y  < 2^63-1 =>  x et y  < 2^31-1 
// --- 
#include <iostream> 
using namespace std;

#ifdef BAMG_LONG_LONG
#define LONG8 long long
#define DLONG8LONG8 1.e17
#define MAXCOOR ((double) 1073741823)// 2^30-1=1073741823. )
#else
#define LONG8 long 
#define DLONG8LONG8 1e9
#define MAXCOOR ((double) 32767. )
#endif
#define DET8(a11,a12,a21,a22) ( (LONG8) (a11) * (LONG8) (a22) - (LONG8) (a21)* (LONG8) (a12))
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
namespace Fem2D {
#define i_nint(x) ((long) x >= 0 ? x + 0.5 : x - 0.5)
#define amin(a, b) ((a) < (b) ? (a) : (b))
#define amax(a, b) ((a) < (b) ? (b) : (a))
#define aabs(a) ((a) < 0 ? -(a) : (a))

static long     nothing = -1073741824L;



/*****************************************************************************
I. Fabrique une triangulation (de nbsmax sommets maxi.) a partir de:
--------------------------------------------------------------------
    -cr: tableau de taille 2*nbs (abscisse & ordonnees des points en entree)
    -arete: tableau de taille 2*nba (depart & arrivee des aretes)
      (decrire le contour dans le sens direct: domaine a
            gauche des aretes exterieures)
    -h: tableau de taille nbs (precision relative aux points en entree)


II. Alloue et affecte une structure t de type triangulation:
------------------------------------------------------------
     -t.np: nombre de sommets en sortie
     -t.nt: nombre de triangles en sortie
     -t.rp[i].x:} abscisse
      t.rp[i].y:} & ordonnee du ieme point du maillage
     -t.tr[i][j]: numero du jeme point du ieme triangle du maillage


III. Renvoie le type d'erreur:
------------------------------
     -  0 : pas d'erreur
     -(-1): pas assez de memoire
     - >0 : erreur mshptg8_
  mshptg erreurs
1: mshptg := 'le nb de point est < 3 or > que le nb max de point';
2: mshptg := 'il y des points confondus ';
3: mshptg := 'tout les points sont align籹 ';
7: mshptg := 'bug dans l''algorithme pour retrouver la fronti焤e (Internal Bug mshfrt)';
8: mshptg := 'Une ar峵e a forcer traverse trop de triangles (Internal Bug : mshfr1)';
9: mshptg := 'la fronti焤e est crois籩 ';
10: mshptg := 'il y a un point qui est sur une ar峵e fronti焤e';
11: mshptg := 'Un sous domaine est reference par aucun triangle ';
20: mshptg := 'mshopt: 3 points confondus (Internal Bug) ';
21: mshptg := 'la pile de mshopt est trop petit (Internal Bug)';



    Remarque: penser a liberer les pointeurs t.rp et t.tr!!!
    --------- un tableau de taille N est indexe de 0 a N-1
        les numeros affectes aux points et aux triangles commencent a 0
******************************************************************************/
/*
long 
MakeTriangulation (triangulation * t, long nbs, long nbsmax, long nba,
		   double *crbdy, double *hbdy, long *arete, int *ngbdy, long *sd, long nbsd, int* flag, int fflag)
{
  int             i, j;
  long            err = 0, nbt, nbsold = nbs;
  long           *c = NULL;
  long           *tri = NULL;
  long           *nu = NULL;
  long           *reft = NULL;
  int            *ngg = NULL;
  double          *cr = NULL;
  double          *h = NULL;

  nbt = 2 * nbsmax;
  nu = new long[6*nbt];
  c = new long[2*nbsmax];
  ngg = new int[nbsmax];
  tri = new long[(4 * nbsmax + 2 * nbsd)];
  reft = new long[nbt];
  cr = new double[(2 * nbsmax + 2)];
  h = new double[nbsmax];

  for (i = 0; i < 2 * nba; i++)
    arete[i]++;
  for (i = 0; i < nbs; i++)
     {
       ngg[i] = ngbdy[i];
       cr[2 * i] = crbdy[2 * i];
       cr[2 * i + 1] = crbdy[2 * i + 1];
       h[i] = hbdy[i];
     }
  for (i = nbs; i < nbsmax; i++)
    ngg[i] = 0;

  mshptg8_ (cr, h, c, nu, &nbs, nbsmax, tri, arete, nba, (long *) sd, nbsd, reft, &nbt, .25, .75, &err);
  for (i = 0; i < 2 * nba; i++)
    arete[i]--;
  if (err)
    goto L1;
  if (*flag)
     {
       delete [] t->rp;t->rp = NULL;
       delete [] t->tr;t->tr = NULL;
       delete [] t->ng;t->ng = NULL;
       delete [] t->ngt;t->ngt = NULL;
     }
  t->rp =new rpoint[nbs];
  t->tr = new triangle[nbt];
  t->ng = new int[nbs];
  t->ngt = new int[nbt];

  *flag = 1;
  t->np = nbs;
  t->nt = nbt;
  for (i = 0; i < nbt; i++)
     {
       for (j = 0; j < 3; j++)
	 t->tr[i][j] = nu[3 * i + j] - 1;
       t->ngt[i] = reft[i] - 1;
     }
  for (i = 0; i < nbs; i++)
     {
       t->rp[i].x = cr[2 * i];
       t->rp[i].y = cr[2 * i + 1];
       if (i < nbsold)
	 t->ng[i] = ngg[i];
       else
	 t->ng[i] = 0;
     }
  renum (t);
  if(!fflag) removeBdyT (t);
L1:
  delete [] nu;nu = NULL;
  delete [] cr;cr = NULL;
  delete [] c;c = NULL;
  delete [] tri;tri = NULL;
  delete [] reft;reft = NULL;
  delete [] ngg;ngg  = NULL;
  delete [] h;h = NULL;
  return err;
}



int 
verifietr (double *cr, int ns)
{
  int             i;
  double           xh, diameter = 1.F;

  if (ns == 0)
    return -1;
  if (ns > 1)
     {
       diameter = 0.F;
       for (i = 0; i < ns; i++)
	 diameter = amax (diameter, aabs (cr[2 * i] - cr[0]) + aabs (cr[2 * i + 1] - cr[1]));
     }
  else
    diameter = 0.001F;
  for (i = 0; i < ns; i++)
     {
       xh = aabs (cr[2 * i] - cr[2 * ns]) + aabs (cr[2 * i + 1] - cr[2 * ns + 1]);
       if (xh < 1e-5 * diameter)
	 return i;
     }
  return -1;
}

*/

int             mshrgl_ (double *c, long *nrfs, long *nbs, long *nu, long *w1,
			 long *w, double omega, long itermx, double eps);
int             mshopt_ (long *c, long *nu, long *t, long a, long *err);
void            mshvoi_ (long *nu, long *w1, long *w, long *nbt, long *nbs);
int             msha1p_ (long *t, long *s, long *c, long *nu, long *reft, long *tete, long *nbt,
			 long *err);
int             mshtri_ (double *cr, long *c, long *nbs, long *tri, LONG8 *nu, double *trfri, long *err);
int             mshcxi_ (long *c, long *nu, long *tri, long *nbs, long *tete, long *err);
int             mshfrt_ (long *c, long *nu, long *nbs, long *arete, long nba, long *sd,
			 long nbsd, long *reft, long *w, long *err);
int             mshgpt_ (long *c, double *cr, long *nu, double *h, long *reft, long *nbs,
    long nbsmx, long *nbt, double coef, double puis, double *trfri, long *err);
long            mshlcl_ (long *c, long *nu, long *tete, long *s);
int             mshtr1_ (LONG8 *criter, long *record, long *n);
int             mshcvx_ (long direct, long *c, long *nu, long *pfold, long *err);
int             mshfr1_ (long *c, long *nu, long *it1, long *ita, long *is1, long *s2, long *err);
int             mshfr2_ (long *c, long *nu, long *lst, long *nbac, long *t, long *ta,
			 long *ss1, long *ss2, long *err);

int 
mshptg8_ (double *cr, double *h, long *c, long *nu, long *nbs, long nbsmx, long *tri,
	 long *arete, long nba, long *sd,
	 long nbsd, long *reft, long *nbt, double coef, double puis, long *err)
{
  /* System generated locals */
  long            i_1;

  /* Local variables */
  static long     tete, i, j, k, t;
  static double    trfri[4];
  static long     nbsgrn, nbtgrn;

/* ----------------------------------------------------------------------- */
/*      but:  construire une triangulation a partir d'un ensemble de */
/*             points et d'un maillage frontalier */
/* ----------------------------------------------------------------------- */
/*     entre : */
/*     ------- */
/*           cr(2,nbsmx)  tableau des coordonnees des nbs points donnes */


/*           h (nbsmx)    tableau du h local voulu autour de chaque point */
/*                          donnes */
/*           nbs          nombre de points donnes */
/*           nbsmx        nombre de points maximal a cree */
/*                        si nbs  = nbsmx ont ne cree pas de points */
/*                        si nbs  < nbsmx => erreur */
/*           arete(2,nba) tableau des aretes du maillage a forcer */
/*                          exemple :la frontiere */
/*           nba          le nombre d'aretes du maillage */
/*           sd(2,nbsd)   tableau definisant les nbsd  sous domaine */
/*                          (reference des triangles gerener) */
/*                          abs(sd(1,i)) =  numero d'une l'arete */
/*                          si sd(1,i) est positive alors le sous domaine */
/*                          est a gauche de l'arete sinon il est a droite */
/*                          sd(2,i) donne le numero du sous doimaine */

/*           puis         coefficient de generation des points */
/*                        .1  => on propage plus loin les rafinement */
/*                               donnes par h */
/*                        .25 => valeur conseillee */
/*           coef         coefficient sur le test arret */
/*                          le valeur conseillee est .75 */
/*                          remarque le nombre de point genere est en */
/*                          O(coef**2) */

/*        tableaux de travail: */
/*        -------------------- */
/*           c(2,nbsmx)    tableau d'entiers (copie de coordonnees) */
/*           tri(ltri)     tableau d'entiers */
/*        out : */
/*        ----- */
/*         nbs         nombre de points   donnes + generes */
/*         nbt         nombre de triangles generes */
/*         cr(1:2,nbs) coordonnees des sommets donnes + generes */
/*         nu(1:3,nbt) sommets des triangles (tableau des connections) */
/*                       telle que les sommets tourne dans le sens direct */
/*         reft(1:nbt) numero de sous domaine de chaque triangle */
/*         err    si err = 0 alors pas de probleme */
/*                sinon nbt = 0 et pas de triangulation */
/*     dimension des tableaux */
/*     ---------------------- */
/*     definition des parameters */
/*     nbtmx = 2*(nbs-1) ,  ltri = max(4*nbs+2*nbsd,nba) */
/*     long : nu(6*nbtmx) , reft(nbtmx) , c(2*nbsmx) , tri(ltri) */
/*     long : arete(2,nba), sd(2,nbsd) */
/*     double    : cr(2*nbsmx) , h(nbsmx) */

/* ---------------------------------------------------------------------- */
/*     programmeur F.Hecht, France */
/*           version 1.0  mars 1986 */
/* ----------------------------------------------------------------------- */

  /* Parameter adjustments */
  --reft;
  sd -= 3;
  arete -= 3;
  --tri;
  --nu;
  c -= 3;
  --h;
  cr -= 3;

  /* Function Body */
  *err = 0;
  *nbt = 0;
  if (*nbs < 3 || nbsmx < *nbs)
     {
       *err = 1;
       printf("mshptg bug in number of points %ld > %ld == max nb points \n",*nbs,nbsmx);
       return 0;
     }
/* ------------------------- */
/* preparation des donnees */
/* ------------------------- */
  mshtri_ (&cr[3], &c[3], nbs, &tri[1], (LONG8*) (void *) &tri[*nbs + 1], trfri, err);
  if (*err != 0)
     {
       return 0;
     }
/* -------------------------------- */
/* maillage de l enveloppe convexe */
/* -------------------------------- */
  mshcxi_ (&c[3], &nu[1], &tri[1], nbs, &tete, err);
/* ----------------------------------------------------------------------- */
/*     definition de tableau nu(1:6,2*nbs-2) */
/* ----------------------------------------------------------------------- */
/*     nu(*,ie) definit soit un element ,soit un sommet frontiere */
/*     si nu(5:6,ie) = (0,0) alors ie est un sommet frontiere */
/*     avec nu(1,ie) = numero du sommet */
/*          nu(2,ie) = 8*t + a */
/*                     ou t est le numero du triangle ayant l'arete */
/*                     frontiere (a) dont le premier sommet est nu(1,ie) */
/*          nu(3,ie) = pointeur dans nu sur sommet frontiere precedent */
/*          nu(4,ie) = pointeur dans nu sur sommet frontiere suivant */
/*     sinon ie est un element : */
/*          nu(1:3,ie) numero des 3 sommets du triangle ie tournant dans */
/*                     le sens direct */
/*          nu(4:6,ie) = (d4,d5,d6) donnee des 3 aretes ai */
/*           ai est forme des sommets nu(i-3,ie),nu(mod(i,3)+1,ie) */
/*           si di < 0 alors arete i est frontiere et -di est pointeur */
/*             sur 1er sommet frontiere de i */
/*           sinon arete est interne et di = 8*ta + ata */
/*              ou ta est le numero du triangle adjacent a l'arete */
/*              et ata est le numero de l'arete dans ta */
/* ------------------------------------------------------------------------ */
  if (*err != 0)
     {
       return 0;
     }

  i_1 = *nbs;
  for (i = 1; i <= i_1; ++i)
     {
       tri[i] = 0;
     }
  i = tete;
L20:
  j = nu[(i - 1) * 6 + 4];
  tri[nu[(i - 1) * 6 + 1]] = nu[(j - 1) * 6 + 1];
  i = j;
  if (i != tete)
     {
       goto L20;
     }
/* ----------------------------- */
/* traitement frontiere */
/* ----------------------------- */
  k = 0;
  mshfrt_ (&c[3], &nu[1], nbs, &arete[3], nba, &sd[3], nbsd, &reft[1], &tri[
								   1], err);
  if (*err != 0)
     {
       return 0;
     }
/* ------------------------------------------------------------------- */
/*       on a modifie nu les sommets frontiere n'ont plus de sens */
/*       ainsi que les pointeurs sur ces elements */
/* ------------------------------------------------------------------- */
  nbsgrn = *nbs;
  mshgpt_ (&c[3], &cr[3], &nu[1], &h[1], &reft[1], &nbsgrn, nbsmx, &nbtgrn,
	   coef, puis, trfri, err);
  if (*err != 0)
     {
       return 0;
     }
/*     construction du tableau nu(1:3,1:nbt) */
/* ------------------------------------------ */
  *nbt = 0;
  k = 0;
  j = 1;
  i_1 = nbtgrn;
  for (t = 1; t <= i_1; ++t)
     {
       if (nu[j + 5] != 0)
	  {
	    ++(*nbt);
	    reft[*nbt] = reft[t];
	    for (i = 0; i <= 2; ++i)
	       {
		 ++k;
		 nu[k] = nu[j + i];
	       }
	  }
       j += 6;
     }
/*     dans nu il y a (s1(t),s2(t),s3(t),t=1,nbt) */
/*     ou s1 s2 s3 sont les 3 sommets de t */
/* ------------------------------------------------ */
  i_1 = *nbs;
  for (i = 1; i <= i_1; ++i)
     {
       tri[i] = 1;
     }
  i_1 = nbsgrn;
  for (i = *nbs + 1; i <= i_1; ++i)
     {
       tri[i] = 0;
     }
  mshvoi_ (&nu[1], &tri[nbsgrn + 1], &nu[*nbt * 3 + 1], nbt, &nbsgrn);
  mshrgl_ (&cr[3], &tri[1], &nbsgrn, &nu[1], &tri[nbsgrn + 1], &nu[*nbt * 3
						      + 1], 1.4F, 20L, .005F);
  *nbs = nbsgrn;
  return 1;
}				/* mshptg_ */

/* ********************************************************************** */
void 
mshvoi_ (long *nu, long *w1, long *w, long *nbt, long *nbs)
{
  /* System generated locals */
  long            i_1;

  /* Local variables */
  static long     i, is;


/* RECHERCHE DU VOISINAGE */
/* ----------------------- */
  /* Parameter adjustments */
  --w;
  --nu;

  /* Function Body */
  i_1 = *nbs;
  for (i = 1; i <= i_1; ++i)
     {
       w1[i] = 0;
     }
  i_1 = *nbt * 3;
  for (i = 1; i <= i_1; ++i)
     {
       ++w1[nu[i]];
     }
  w1[0] = 0;
  i_1 = *nbs;
  for (i = 1; i <= i_1; ++i)
     {
       w1[i] = w1[i - 1] + w1[i];
     }
  i_1 = *nbt * 3;
  for (i = 1; i <= i_1; ++i)
     {
       is = nu[i] - 1;
       ++w1[is];
       w[w1[is]] = i;
     }
  for (i = *nbs; i >= 1; --i)
     {
       w1[i] = w1[i - 1];
     }
  w1[0] = 0;
}				/* mshvoi_ */

/* ********************************************************************** */
int 
mshrgl_ (double *c, long *nrfs, long *nbs, long *nu, long *w1,
	 long *w, double omega, long itermx, double eps)
{
  /* System generated locals */
  long            i_1, i_2, i_3;
  double           r_1, r_2;

  /* Local variables */
  static double    depx, depy;
  static long     iter;
  static double    xmin, ymin, xmax, ymax;
  static long     i, k, i1, i2, ic;
  static double    bx, by, dx;
  static long     is;
  static double    err;


/* REGULARISATION PAR MOYENNE BARYCENTRIQUE */
/* ----------------------------------------- */
  /* Parameter adjustments */
  --w;
  --nu;
  --nrfs;
  c -= 3;

  /* Function Body */
  xmin = c[3];
  ymin = c[4];
  xmax = c[3];
  ymax = c[4];
  i_1 = *nbs;
  for (ic = 1; ic <= i_1; ++ic)
     {
/* Computing MIN */
       r_1 = c[(ic << 1) + 1];
       xmin = amin (r_1, xmin);
/* Computing MIN */
       r_1 = c[(ic << 1) + 2];
       ymin = amin (r_1, ymin);
/* Computing MAX */
       r_1 = c[(ic << 1) + 1];
       xmax = amax (r_1, xmax);
/* Computing MAX */
       r_1 = c[(ic << 1) + 2];
       ymax = amax (r_1, ymax);
     }
/* Computing MAX */
  r_1 = xmax - xmin, r_2 = ymax - ymin;
  dx = amax (r_1, r_2);
  i_1 = itermx;
  for (iter = 1; iter <= i_1; ++iter)
     {
       err = (double) 0.;
       i2 = w1[0];
       i_2 = *nbs;
       for (is = 1; is <= i_2; ++is)
	  {
	    i1 = i2 + 1;
	    i2 = w1[is];
	    if (i2 >= i1 && nrfs[is] == 0)
	       {
		 bx = (double) 0.;
		 by = (double) 0.;
		 i_3 = i2;