Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2798) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2799) @@ -814,7 +814,4 @@ // Triangles(Triangles &){ std::cerr << " BUG call copy opretor of Triangles" << std::endl;MeshError(111);} Triangles(const Triangles &,const int *flag,const int *bb); // truncature - - int gibbsv (long int* ptvoi, long int* vois,long int* lvois,long int* w,long int* v); - int gibbs(void); void SetIntCoor(const char * from =0); Index: sm/trunk/src/c/Bamgx/Meshgibbs.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Meshgibbs.cpp (revision 2798) +++ (revision ) @@ -1,1064 +1,0 @@ -// -*- Mode : c++ -*- -// -// SUMMARY : -// USAGE : -// ORG : -// AUTHOR : Frederic Hecht -// E-MAIL : hecht@ann.jussieu.fr -// - -/* - - 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 "./Mesh2.h" - -#include -#include -#include -//#include "vect.h" - - -#define mmax(a,b)(a>b?a:b) -#define mmin(a,b)(a= 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; - bnew = mmax(i__2,i__3); -/* L110: */ - } -/* if(impre.ne.0) then */ -/* write(nfout,*)'profile old = ',pfold,', profile new = ',pfnew - */ -/* write(nfout,*)'1/2 bande old = ',bold ,', 1/2 band new = ',bnew -*/ -/* endif */ -return 0; -} /* gibbsa_ */ - -/* Subroutine */ 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) -{ - /* System generated locals */ - integer i__1, i__2; - - /* Local variables */ - static int flag_; - static integer i, j, k, p, s, h0, i1, l0, i2; -/* extern Subroutine int gibbs1_(); */ - static integer lg; -/* extern Subroutine int gibbsd_(), gibbsc_();*/ - static integer niveau, mxcanx, mxcany, nbc; - -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - */ -/* ...................................................................... -*/ -/* attention on met la descente optimiser dans r <0 ou nulle */ -/* ....................................................................... - */ - /* Parameter adjustments */ - --r; - --m; - --ny; - --nx; - --vois; - --ptvois; - - /* Function Body */ - gibbsd_(y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau); - gibbsc_(&ny[1], nv, &niveau, n, &mxcany); - gibbsd_(x, n, &ptvois[1], &vois[1], nv, &r[1], &niveau); - p = niveau; - gibbsc_(&nx[1], nv, &niveau, n, &mxcanx); - flag_ = ffalse; - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - if (nx[i] + ny[i] == p) { - r[i] = -nx[i]; - } else if (nx[i] >= 0) { - flag_ = ttrue; - r[i] = -1073741824; - } else { - if (r[i] <= 0) { - r[i] = -1073741822; - } - } -/* L20: */ - } - if (flag_) { -/* calcul des composantes connexe du graphe sans les sommets de -nn */ -/* ------------------------------------------------------------ ---- */ - j = *n; - k = 0; - nbc = 0; - nv[nbc] = j; -L30: - ++k; - if (k <= *n) { - if (r[k] == -1073741824) { -/* recherche de la fermeture transitive partant de k - */ - ++nbc; - i = -1; - s = k; -L40: - ++i; - wl[i] = ptvois[s]; - wh[i] = ptvois[s + 1]; - ++j; - nv[j] = s; - r[s] = -1073741823; -L50: - if (i >= 0) { - if (wl[i] < wh[i]) { - s = vois[wl[i]]; - ++wl[i]; - if (r[s] == -1073741824) { - goto L40; - } - goto L50; - } - --i; - goto L50; - } - nv[nbc] = j; - m[nbc] = nbc; - } - goto L30; - } -/* if(impre.lt.0) write(nfout,*) */ -/* + ' nb de composante connexe du graphe reduit =',nbc */ - -/* --------------- fin de construction des composantes connexes------ ---- */ -/* nv(0)=n */ -/* if(impre.le.-10) write(nfout,5555)'nv(0:n+n) = ',(nv(i),i=0, -n+n) */ - gibbs1_(&nbc, &m[1], nv); -/* if(impre.le.-10)write(nfout,5555)'trie m =',(m(i),i=1,nbc) -*/ - i__1 = p; - for (i = 0; i <= i__1; ++i) { - nn[i] = 0; -/* L60: */ - } - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - j = -r[i]; - if (j >= 0 && j <= p) { - ++nn[j]; - } -/* L70: */ - } - -/* boucle sur les composante connexes par ordre croissantes */ -/* -------------------------------------------------------- */ - for (k = nbc; k >= 1; --k) { - i = m[k]; - i1 = nv[i - 1] + 1; - i2 = nv[i]; - lg = i2 - i1 + 1; -/* if(impre.le.-7) */ -/* + write(nfout,*) k,' composante ',i,',lg=',lg,',i1,i2 -=',i1,i2 */ -/* if(impre.le.-8) */ -/* + write (nfout,5555)' ',(nv(i),i=i1,i2) */ - h0 = 0; - l0 = 0; - i__1 = p; - for (j = 0; j <= i__1; ++j) { - wh[j] = nn[j]; - wl[j] = nn[j]; -/* L90: */ - } - i__1 = i2; - for (i = i1; i <= i__1; ++i) { - s = nv[i]; - ++wh[nx[s]]; - ++wl[p - ny[s]]; -/* L100: */ - } - i__1 = p; - for (j = 0; j <= i__1; ++j) { - if (wh[j] != nn[j]) { -/* Computing MAX */ - i__2 = wh[j]; - h0 = mmax(i__2,h0); - } - if (wl[j] != nn[j]) { -/* Computing MAX */ - i__2 = wl[j]; - l0 = mmax(i__2,l0); - } -/* L110: */ - } - if (h0 < l0 || h0 == l0 && mxcanx <= mxcany) { -/* if(impre.le.-2) write(nfout,*) */ -/* + ' h0 = ',h0,',l0 = ',l0,' ------- XXXX - --------' */ - i__1 = i2; - for (i = i1; i <= i__1; ++i) { - s = nv[i]; - r[s] = -nx[s]; - ++nn[-r[s]]; -/* L120: */ - } - } else { -/* if (impre.le.-2) write(nfout,*) */ -/* + ' h0 = ',h0,',l0 = ',l0,' ------- YYYY - --------' */ - i__1 = i2; - for (i = i1; i <= i__1; ++i) { - s = nv[i]; - r[s] = -p + ny[s]; - ++nn[-r[s]]; -/* L130: */ - } - } -/* L140: */ - } - } -/* on met les nouveaux niveaux de la descendance optimiser dans nn */ -/* ----------------------------------------------------------------- -*/ - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - if (r[i] > 0) { - nn[i] = -1; - } else if (r[i] == -1073741822) { - nn[i] = -2; - } else { - nn[i] = -r[i]; - } -/* L150: */ - } -/* if(impre.le.-10)write (nfout,5555)' nn(i)=',(nn(i),i=1,n) */ -/* 5555 format(' -------- ',a,/,5(15x,10(i5)/)) */ -return 0;} /* gibbsb_ */ - -/* Subroutine */ int gibbsc_(integer* nz,integer* nv,integer* niveau,integer* n,integer* mxz) -{ - /* System generated locals */ - integer i__1, i__2, i__3; - - /* Local variables */ - static integer i, j; - -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - */ - /* Parameter adjustments */ - --nz; - - /* Function Body */ - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - nz[i] = -1; -/* L10: */ - } - *mxz = 0; - i__1 = *niveau; - for (i = 0; i <= i__1; ++i) { -/* Computing MAX */ - i__2 = *mxz, i__3 = nv[i + 1] - nv[i]; - *mxz = mmax(i__2,i__3); - i__2 = nv[i + 1]; - for (j = nv[i] + 1; j <= i__2; ++j) { - nz[nv[j]] = i; -/* L20: */ - } - } -return 0;} /* gibbsc_ */ - -/* Subroutine */ int gibbsd_(integer* racine,integer* n,integer* ptvois,integer* vois,integer* nv,integer* r,integer* niveau) -{ - /* System generated locals */ - integer i__1, i__2; - - /* Local variables */ - static integer i, k, s, sv, stk, stk1, stk2; - -/* ----------------------------------------------------------------------- - */ -/* but construire la structure des descendant de racine du graphe */ -/* ----------------------------------------------------------------------- - */ -/* sortie : */ -/* -------- */ -/* nv est la structure des niveaux */ -/* les sommets du niveau (i =0,niveau_ sont defini par : */ -/* (nv(j),j=nv(i),nv(i+1)-1) */ - -/* le tableau r(i) n'est modifier que sur les sommets */ -/* de la composante connexe du graphe contenant la racine */ - -/* ----------------------------------------------------------------------- - */ - -/* on demark tout les sommets non remuneroter */ -/* -------------------------------------------------- */ - /* Parameter adjustments */ - --r; - --vois; - --ptvois; - - /* Function Body */ - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - if (r[i] < 0) { - r[i] = 0; - } -/* L10: */ - } - -/* initialisation */ - - stk = *n - 1; - nv[0] = stk; - stk2 = stk; - *niveau = 0; - ++stk; - nv[stk] = *racine; - r[*racine] = -1; -L20: - if (stk2 < stk) { - ++(*niveau); - stk1 = stk2 + 1; - nv[*niveau] = stk; - stk2 = stk; -/* print *,' ------- niveau =',niveau,' stk=',stk1,stk2 */ - i__1 = stk2; - for (k = stk1; k <= i__1; ++k) { - s = nv[k]; -/* print *,'----------------- s=',s */ - i__2 = ptvois[s + 1] - 1; - for (i = ptvois[s]; i <= i__2; ++i) { -/* pour tout les sommets (sv) voisin */ -/* d'un sommet (s) du niveau precedent */ - sv = vois[i]; -/* print *,' voisin =',sv */ -/* si le sommet n'est pas marque on le marque et - on l'ajout */ - if (r[sv] == 0) { - ++stk; - nv[stk] = sv; - r[sv] = -1; - } -/* L30: */ - } -/* L40: */ - } - goto L20; - } - --(*niveau); -/* call pnv(' gibbsd ',n,nv,niveau) */ -return 0;} /* gibbsd_ */ - - -/* Subroutine */ 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) -{ - /* System generated locals */ - integer i__1, i__2, i__3, i__4, i__5; - - /* Local variables */ - static integer nbsc, bnew, knew, step, plus, i, j, k, s, debut, i1, i2; -/* extern Subroutine int gibbs2_();*/ - static integer fin; - - -/* construction de la stucture de niveau dans nv a partir de nn */ -/* ------------------------------------------------------------ */ - /* Parameter adjustments */ - --r; - --m; - --vois; - --ptvois; - - /* Function Body */ - nv[0] = *n; - i__1 = *p + 1; - for (i = 1; i <= i__1; ++i) { - nv[i] = 0; -/* L150: */ - } - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - if (nn[i] >= 0) { - ++nv[nn[i] + 1]; - } -/* L160: */ - } - i__1 = *p; - for (i = 0; i <= i__1; ++i) { - nv[i + 1] += nv[i]; -/* L170: */ - } - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - if (nn[i] >= 0) { - j = nn[i]; - ++nv[j]; - nv[nv[j]] = i; - } -/* L180: */ - } - for (i = *p; i >= 0; --i) { - nv[i + 1] = nv[i]; -/* L190: */ - } - nv[0] = *n; - nbsc = nv[*p + 1] - nv[0]; -/* --- fin de la construction ------------------------------------ */ - if (*option == -2) { - i__1 = *impre - 1; - } - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - m[i] = *n * 3 + ptvois[i + 1] - ptvois[i]; -/* L10: */ - } - if ((((int)*option) == 1)||(((int)*option) == -1)) { - debut = 0; - fin = *p; - step = 1; - } else { - debut = *p; - fin = 0; - step = -1; - } - i__1 = fin; - i__2 = step; - for (i = debut; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) { - i1 = nv[i] + 1; - i2 = nv[i + 1]; - i__3 = i2 - i1 + 1; - gibbs2_(&i__3, &nv[i1], &m[1]); - i__3 = i2; - for (j = i1; j <= i__3; ++j) { - s = nv[j]; - i__4 = ptvois[s + 1] - 1; - for (k = ptvois[s]; k <= i__4; ++k) { -/* Computing MIN */ - i__5 = m[vois[k]]; - m[vois[k]] = mmin(i__5,j); -/* L20: */ - } -/* L30: */ - } -/* L40: */ - } - if (*option > 0) { - knew = *new_; - plus = 1; - } else { - knew = *new_ + nbsc + 1; - plus = -1; - } - *new_ += nbsc; -/* if(option.gt.0) then */ -/* do 60 k = debut , fin , step */ -/* do 60 j = nv(k+1),nv(k)+1,-1 */ -/* knew = knew + plus */ -/* r(nv(j)) = knew */ -/* 60 continue */ -/* else */ - i__2 = fin; - i__1 = step; - for (k = debut; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) { - i__3 = nv[k + 1]; - for (j = nv[k] + 1; j <= i__3; ++j) { - knew += plus; - r[nv[j]] = knew; -/* L70: */ - } - } -/* endif */ - *pfnew = 0; - bnew = 0; - i__3 = *n; - for (i = 1; i <= i__3; ++i) { - k = r[i]; - if (k > 0) { - i__1 = ptvois[i + 1] - 1; - for (j = ptvois[i]; j <= i__1; ++j) { - if (r[vois[j]] > 0) { -/* Computing MIN */ - i__2 = k, i__4 = r[vois[j]]; - k = mmin(i__2,i__4); - } -/* L100: */ - } - *pfnew = *pfnew + r[i] - k + 1; -/* Computing MAX */ - i__1 = bnew, i__2 = r[i] - k + 1; - bnew = mmax(i__1,i__2); - } -/* L110: */ - } -/* if(impre.lt.0.or.impre.gt.2) then */ -/* write(nfout,*) ' option =',option,', profile =',pfnew */ -/* + ,', 1/2 bande =',bnew,', new=',new,', nbss composante=',nbsc - */ -/* endif */ -return 0;} /* gibbst_ */ - -/* function */ -int bamg::Triangles::gibbsv (integer* ptvoi, integer* vois,integer* lvois,integer* w,integer* v){ - /* System generated locals */ - integer i__2; - - /* Local variables */ - integer i, j, k, T, ss, iii, ptv, ptv1; - integer nbss = nbv; - /*--- Prepare les donees pour gibbsa en construisant ptvoi, vois, lvois - - ------------*/ - /* in */ - /* --- nbnt =3 pour des triangles 2D, */ - /* nbt = nb de triangle */ - /* nbss = nb de sommets */ - /* nsea = numeros de 3 sommets de chaque triangle (me) */ - /* out */ - /* --- ptvoi, vois, lvois, err */ - /* tableaux de travail w, v */ - /*----------------------------------------------------------------------- - ----------*/ - /* Parameter adjustments */ - --v; - --w; - --vois; - --ptvoi; - long nt = nbt-NbOutT; - /* Function Body */ - for (i = 1; i <= nbss; ++i) { - w[i] = -1; - ptvoi[i] = 0; } - ptvoi[nbss + 1] = 0; - for (i = 0; i < nt; ++i) - { - assert(triangles[i].link); - for (j = 0; j < 3; ++j) - { - ss = Number(triangles[i][j])+1; - ++ptvoi[ss + 1]; - w[ss] = 0; - } - } - - for (i = 1; i <= nbss; ++i) - ptvoi[i + 1] += ptvoi[i]; - - for (i = 0; i < nt; ++i) - if (triangles[i].link) - for (j = 0; j < 3; ++j) - { - ss = Number(triangles[i][j])+1; - ++ptvoi[ss]; - v[ptvoi[ss]] = i; - } - - ptv1 = 0; - iii = 1; - for (i = 1; i <= nbss; ++i) { - ptv = ptv1 + 1; - ptv1 = ptvoi[i]; - ptvoi[i] = iii; - i__2 = ptv1; - for (j = ptv; j <= i__2; ++j) { - T = v[j]; - for (k = 0; k < 3; ++k) { - ss = Number(triangles[T][k])+1; /* nsea[k + T * nsea_dim1]; */ - if (w[ss] != i) { - w[ss] = i; - if (iii > *lvois) return 2 ; - /* print*,'pas assez de place memoire' */ - - vois[iii] = ss; - ++iii;} - } - } - } - ptvoi[nbss + 1] = iii; - *lvois = iii - 1; - return 0; /* OK */ - return 0;} /* gibbsv_ */ - -int bamg::Triangles::gibbs(void) -/* -------- - renumber vertices by gibbs method; updates triangle and edge array - in: mesh - out: mesh - auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f - all of size nv+1 except vois (10(nv+1)) and nv (2(nv+1)) - err = -1 : memory alloc pb; err = -3: fatal erreur gibbs 2 : pb racine -*/ -{ - long nv = nbv; - long nt = nbt-NbOutT; - long i, j, pfold, pfnew; - long* ptvois=NULL; - long* vois=NULL; - long* nn =NULL; - long* r =NULL; - long* m =NULL; - long* nnv =NULL; - long* nx =NULL; - long* ny =NULL; - long* w1 =NULL; - long* w2=NULL; - long nbvoisin = 10*nv; - long printint=0, iodev=6; - int err=0; - ptvois = new long[nv+1]; //(long*)calloc((long)(nv + 1) , sizeof(long)); - nn = new long[3*nt]; //(long*)calloc(3 * nt ,sizeof(long)); - vois = new long[nbvoisin+10]; //(long*)calloc((long)(nbvoisin + 10) , sizeof(long)); - r = new long[nv+1]; //(long*)calloc((long)(nv + 1) , sizeof(long)); - if((!ptvois)||(!nn)||(!vois)||(!r)) return -1; - err = gibbsv(ptvois,vois,&nbvoisin,r,nn) ; - delete [] nn; // free(nn); - if(err==0) - { - m = new long[nv+1]; - nn = new long[nv+1]; - nnv = new long[(nv+1)<<1]; - nx = new long[nv+1]; - ny = new long[nv+1]; - w1 = new long[nv+1]; - w2 = new long[nv+1]; - long lnv = nv; - err = gibbsa_ (&lnv, ptvois, vois, r, m, nnv, nx, ny, nn, w1, w2, &pfold, &pfnew, - &printint, &iodev); - delete [] m; - delete [] nnv; - delete [] nn; - delete [] nx; - delete [] ny; - delete [] w1; - delete [] w2; - } - - delete [] vois; - delete [] ptvois; - /* - if (err == 0 && (pfnew <= pfold)) - { - A f(nv); - for (i = 0; i < nv; ++i) - { f[i].x = v[i].x; - f[i].y = v[i].y; - f[i].where = v[i].where; - } - for (i = 0; i < nv; ++i) - { v[r[i] - 1].x = f[i].x; - v[r[i] - 1].y = f[i].y; - v[r[i] - 1].where = f[i].where; - } - - for (j = 0; j < nt; ++j) // updates triangle array - for (i = 0; i < 3; i++) - t[j].v[i] = &v[r[no(t[j].v[i])] - 1]; - - for (j = 0; j < ne; ++j) // updates edge array - { - e[j].in = &v[r[no(e[j].in)] - 1]; - e[j].out = &v[r[no(e[j].out)] - 1]; - } - f.destroy(); - if (!NumThinGrid) - { NumThinGrid= new int [nv]; - for (i=0;i