gibbs.cpp

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

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// 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
 */
#include <cmath>
#include "error.hpp"
#include <iostream>
#include <fstream>
#include "RNM.hpp"
#include "rgraph.hpp"
#include "fem.hpp"
using namespace Fem2D;
#include "FESpacen.hpp"
#include "FESpace.hpp"

#define mmax(a,b)(a>b?a:b)
#define mmin(a,b)(a<b?a:b)
#define ffalse 0
#define ttrue 1
extern long verbosity;
/*  -- translated by f2c (version of 23 May 1992  14:18:33).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/


#define integer long
#define logical long

int gibbs1_(integer* n,integer*  record,integer*  ptvois);
int gibbs2_(integer* n,integer*  record,integer*  criter);
int gibbsa_(integer* n,integer*  ptvois,integer*  vois,integer*  r,integer*  m,
             integer*  nv,integer*  nx,integer*  ny,integer*  nn,integer*  w1,integer*  w2, 
	integer* pfold,integer*  pfnew,integer*  impre,integer*  nfout);
int gibbsb_(integer* x,integer*  y,integer*  n,integer*  ptvois,
						integer*  vois,integer*  nx,integer*  ny,integer*  nv,integer*  nn,integer*  m,
						integer*  wh,integer*  wl,integer* r, integer* impre, integer* nfout);
int gibbsc_(integer* nz,integer*  nv,integer*  niveau,integer*  n,integer* );
int gibbsd_(integer* racine,integer*  n,integer*  ptvois,integer*  
							vois,integer*  nv,integer*  r,integer*  niveau);
int gibbst_(integer* n,integer*  p,integer*  nv,integer*  nn,integer*  ptvois,integer*  vois,
						integer*  m,integer*  r,integer*  new_,integer*  option, 
						integer* pfnew,integer*  impre,integer*  nfout);

/* Subroutine */ int gibbs1_(integer* n,integer*  record,integer*  ptvois)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer crit, i, j, l, r, rec;

/* -----------------------------------------------------------------------
 */
/*     routine appele par gibbs0 */
/* -----------------------------------------------------------------------
 */
/*     but: trie record (ensemble de n sommet) telle que l'ordre des somm 
*/
/*     soit croissant (ordre du sommet i est ptvois(i+1)-ptvois(i)) */
/* -----------------------------------------------------------------------
 */

    /* Parameter adjustments */
    --ptvois;
    --record;

    /* Function Body */
    if (*n <= 1) {
	return 0;
    }
    l = *n / 2 + 1;
    r = *n;
L2:
    if (l <= 1) {
	goto L20;
    }
    --l;
    rec = record[l];
    crit = ptvois[record[l] + 1] - ptvois[record[l]];
    goto L3;
L20:
    rec = record[r];
    crit = ptvois[record[r] + 1] - ptvois[record[r]];
    record[r] = record[1];
    --r;
    if (r == 1) {
	goto L999;
    }
L3:
    j = l;
L4:
    i = j;
    j <<= 1;
    if ((i__1 = j - r) < 0) {
	goto L5;
    } else if (i__1 == 0) {
	goto L6;
    } else {
	goto L8;
    }
L5:
    if (ptvois[record[j] + 1] - ptvois[record[j]] < ptvois[record[j + 1] + 1] 
	    - ptvois[record[j + 1]]) {
	++j;
    }
L6:
    if (crit >= ptvois[record[j] + 1] - ptvois[record[j]]) {
	goto L8;
    }
    record[i] = record[j];
    goto L4;
L8:
    record[i] = rec;
    goto L2;
L999:
    record[1] = rec;
    return 0;
} /* gibbs1_ */

/* Subroutine */ int gibbs2_(integer* n,integer*  record,integer*  criter)
{
    static integer crit, i, j, l, r, rec;


/*     trie record selon les valeurs de criter(record(.)) croissantes */


    /* Parameter adjustments */
    --criter;
    --record;

    /* Function Body */
    if (*n <= 1) {
	return 0;
    }
    l = *n / 2 + 1;
    r = *n;
L2:
    if (l <= 1) {
	goto L20;
    }
    --l;
    rec = record[l];
    crit = criter[rec];
    goto L3;
L20:
    rec = record[r];
    crit = criter[rec];
    record[r] = record[1];
    --r;
    if (r == 1) {
	goto L999;
    }
L3:
    j = l;
L4:
    i = j;
    j <<= 1;
    if (j - r < 0) {
	goto L5;
    } else if (j == r) {
	goto L6;
    } else {
	goto L8;
    }
L5:
    if (criter[record[j]] < criter[record[j + 1]]) {
	++j;
    }
L6:
    if (crit >= criter[record[j]]) {
	goto L8;
    }
    record[i] = record[j];
    goto L4;
L8:
    record[i] = rec;
    goto L2;
L999:
    record[1] = rec;
    return 0;
} /* gibbs2_ */

/* Subroutine */ 
int gibbsa_(integer* n,integer*  ptvois,integer*  vois,integer*  r,integer*  m,
             integer*  nv,integer*  nx,integer*  ny,integer*  nn,integer*  w1,integer*  w2, 
	integer* pfold,integer*  pfnew,integer*  impre,integer*  nfout)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;

    /* Builtin functions */
    /* Subroutine */ int s_stop();

    /* Local variables */
    static integer nbcc, degi, bold, bnew, i, j, k, p, degre, x, y, p1, p2;
/*    extern  Subroutine  int gibbs1_();*/
    static integer pf;
/*    extern  Subroutine int gibbsb_(), gibbsd_(), gibbst_();*/
    static integer nbpass, niveau, pf1, option, old, new_, opt, new1;

/* -----------------------------------------------------------------------
 */
/*  but: calculer une renumerotation des sommets d'un graphe defini par: 
*/
/*     par la methode de gibbs */
/* -----------------------------------------------------------------------
 */
/*  entree */
/* -------- */
/*     n = nb de sommet du graphe */
/*      les voisins d'un sommet i ont pour numero : */
/*     ( vois(j) , j=ptvois(i),ptvois(i+1)-1 ) */

/*     impre   parametre d'impression */
/*     nfout   numero du fichier pour impression */

/*  sortie */
/*  ------ */
/*     r(1:n) tableau donnant la nouvelle numerotation: */
/*       r(i) = nouveau numero du sommet i */
/*     pfolf = ancien  profile */
/*     pfnew = nouveau profile */

/*  tableau de travail : */
/*  -------------------- */
/*     m(n) */
/*     nv(0:n+n) */
/*     nx(n) */
/*     ny(n) */
/*     nn(0:n) */
/*     w1(n) */
/*     w2(n) */

/* -----------------------------------------------------------------------
 */
/*     programmeur f. hecht  le 3/02/1987 */
/* -----------------------------------------------------------------------
 */

/*     tri des voisins d'un sommet du graphe par degre croissant */
/* --------------------------------------------------------------- */
    /* Parameter adjustments */
    --w2;
    --w1;
    --ny;
    --nx;
    --m;
    --r;
    --vois;
    --ptvois;

    /* Function Body */
    p2 = ptvois[1] - 1;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	p1 = p2 + 1;
	p2 = ptvois[i + 1] - 1;
	i__2 = p2 - p1 + 1;
	gibbs1_(&i__2, &vois[p1], &ptvois[1]);
/*       if(impre.le.-9) then */
/*        write (nfout,*) 'les voisin de ',i,'sont: ', (vois(j),j=p1,p
2) */
/*       endif */
/* L10: */
    }
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	r[i] = 0;
/* L20: */
    }
/*     boucle sur les composante connexe du graphe */
    new_ = 0;
    nbcc = 0;
L30:
    if (new_ < *n) {
	++nbcc;
/*       recherche d'une racine y (un sommet non numerote) de degree m
ini */
	y = 0;
	degre = *n + 1;
	i__1 = *n;
	for (i = 1; i <= i__1; ++i) {
	    if (r[i] <= 0) {
		degi = ptvois[i + 1] - ptvois[i];
		if (degi < degre) {
		    degre = degi;
		    y = i;
		}
	    }
/* L40: */
	}
	if (y == 0) {
	   return -3;/*  s_stop("fatal erreur  gibbs 2 : pb racine", 33L); */
	}
	gibbsd_(&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
	nbpass = 0;
L50:
	++nbpass;
	x = y;
	p = niveau;
	k = 0;
	i__1 = nv[p + 1];
	for (i = nv[p] + 1; i <= i__1; ++i) {
	    ++k;
	    m[k] = nv[i];
/* L60: */
	}
	gibbs1_(&k, &m[1], &ptvois[1]);
	i__1 = k;
	for (i = 1; i <= i__1; ++i) {
	    y = m[i];
	    gibbsd_(&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
	    if (niveau > p) {
		goto L50;
	    }
/* L70: */
	}
	y = m[1];
/*        if(impre.lt.0) then */
/*          write(nfout,*) */
/*     +  '    nb de pass pour trouver le pseudo diametre',nbpass */
/*     +         ,' x=',x,',y=',y,' de la composante connexe ',nbcc */

/*          write (nfout,*) ('-',i=1,78) */
/*        endif */
/*       optimisation de la descendance de la numerotation */
/*       ------------------------------------------------- */
	gibbsb_(&x, &y, n, &ptvois[1], &vois[1], &nx[1], &ny[1], nv, nn, &m[1]
		, &w1[1], &w2[1], &r[1], impre, nfout);

/*     renumerotation de cuthill mac kee avec la meilleur des 4 option
s */
/*     --------------------------------------------------------------
--- */
	pf = 1073741824;
	option = -2;
	new1 = new_;
	for (opt = -2; opt <= 2; ++opt) {
	    new_ = new1;
	    if (opt != 0) {
		gibbst_(n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &
			new_, &opt, &pf1, impre, nfout);
		if (pf1 < pf) {
		    pf = pf1;
		    option = opt;
		}
	    }
/* L80: */
	}
/*        if(impre.ne.0) write (nfout,*) '    on a choisi l''option ',
 */
/*     +             option,', new =',new */
	new_ = new1;
	gibbst_(n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &new_, &
		option, &pf1, impre, nfout);
	goto L30;
    }
/*      if(impre.ne.0) write(nfout,*) */
/*     +       '   nb de composante connexe du graphe =',nbcc */
/*     calcul du profile */
    *pfold = 0;
    *pfnew = 0;
    bnew = 0;
    bold = 0;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	old = i;
	new_ = r[i];
	i__2 = ptvois[i + 1] - 1;
	for (j = ptvois[i]; j <= i__2; ++j) {
/* Computing MIN */
	    i__3 = old, i__4 = vois[j];
	    old = mmin(i__3,i__4);
/* Computing MIN */
	    i__3 = new_, i__4 = r[vois[j]];
	    new_ = mmin(i__3,i__4);
/* L100: */
	}
	*pfold = *pfold + i - old + 1;
/* Computing MAX */
	i__2 = bold, i__3 = i - old + 1;
	bold = mmax(i__2,i__3);
	*pfnew = *pfnew + r[i] - new_ + 1;
/* Computing MAX */
	i__2 = bnew, i__3 = r[i] - new_ + 1;