file: Nektar2d/src/bwoptim.c
/*---------------------------------------------------------------------------*
* RCS Information *
* *
* $Source: /users/tcew/Hybrid/Nektar2d/src/RCS/bwoptim.C,v $
* $Revision: 1.1 $
* $Date: 1996/06/22 19:49:44 $
* $Author: tcew $
* $State: Exp $
*---------------------------------------------------------------------------*/
#include <math.h>
#include <veclib.h>
#include "nektar.h"
static void RenumberElmts (Element *E, Bsystem *Bsys, int *newmap);
static Fctcon *FacetConnect (Bsystem *B, Element *E); //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
/* this routine optimises the ordering for the unknown degrees of freedom
it assume that the knowns have already been stored.
Note: asumes that the unknown vertices are have gid's 0 < gid < nvs
Note: asumes that the unknown edges are have gid's 0 < gid < nes
Note: asumes that the unknown facess are have gid's 0 < gid < nfs */
void bandwidthopt(Element *E, Bsystem *Bsys){
int nsols, *newmap;
Fctcon *ptcon; //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
#if DIM == 2
nsols = Bsys->nv_solve + Bsys->ne_solve;
#else
nsols = Bsys->nv_solve + Bsys->ne_solve + Bsys->nf_solve;
#endif
if(Bsys->nv_solve > 1 && Bsys->ne_solve > 1){
newmap = ivector(0,nsols-1);
/* Initially set up "point" to "point" connections (a point/facet is //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?); point: Coords.c(?), Curvi.c(?); point: Coords.c(?), Curvi.c(?); facet: nekcomp.h(?), nektar.h(?)
a vertex, edge or face) */
ptcon = FacetConnect(Bsys,E);
/* Given a list of "nsols" points whose interconnectivity is given in
ptcon this routine returns a new mapping in newmap which contains
the order in which the unknowns should be stored */
MinOrdering(nsols,ptcon,newmap);
/* Sort out the numbering according to new ordering in newmap */
RenumberElmts(E, Bsys, newmap);
free_Fctcon(ptcon,nsols);
free(newmap);
}
}
static int Cuthill(Fctcon *ptcon, int *newmap, //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
int *initmap, int nsols, int initpt);
void MinOrdering(int nsols, Fctcon *ptcon, int *newmap){ //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
register int i,j,k,n;
int firstv, *initmap;
/********** Cuthill-McKee optimisation *********/
initmap = ivector(0, nsols-1);
for(i = 0; i < nsols; ++i) initmap[i] = i;
/* to optimise search try algorithm from all points that have
a minimum number of connecting points */
/* order list so that the points with least degree are first */
for(i = 0; i < nsols; ++i){
k = ptcon[initmap[0]].ncon;
n = 0;
for(j = 1; j < nsols-i; ++j)
if(ptcon[initmap[j]].ncon > k){
n = j;
k = ptcon[initmap[j]].ncon;
}
k = initmap[n];
initmap[n] = initmap[nsols-i-1];
initmap[nsols-i-1] = k;
}
/* perform reorganisation for each peripheral vertex */
/* find number of points with minimum connectivity */
for(n=0;n<nsols;++n)
if(ptcon[initmap[n]].ncon != ptcon[initmap[0]].ncon)
break;
/* find Cuthill with mimimum bandwidth using points with min. connect */
for(firstv = 0, k = nsols, i = 0; i < n; ++i){
j = Cuthill(ptcon,newmap,initmap,nsols,i);
if(j < k){
k = j; firstv = i;
}
}
/* finally perform optimal reorganisation */
Cuthill(ptcon,newmap,initmap,nsols,firstv);
free(initmap);
}
static int Cuthill(Fctcon *ptcon, int *newmap, int *initmap, int nsols, //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
int initpt){
register int i,j,k;
int bw=0,start=0,end=1;
int cnt,len,jtmp,beg;
int next,cand,*mark;
Facet *f; //OVERLOAD CALL: Facet: nektar.h(?), nekcomp.h(?)
if(nsols <= 1) return bw;
mark = ivector(0,nsols);
ifill(nsols,1,mark,1);
izero(nsols,newmap,1);
cnt = 1;
newmap[0] = initmap[initpt];
mark[initmap[initpt]]=0;
while(cnt < nsols){
for(i = start; i < end; ++i){
next = newmap[i];
beg = cnt;
for(f = ptcon[next].f; f; f = f->next){
cand = f->id;
if(mark[cand]){
mark[cand] = 0;
newmap[cnt] = cand;
++cnt;
}
}
for(j = beg; j < cnt; ++j){
len = ptcon[newmap[beg]].ncon;
jtmp = beg;
for(k = beg; k < cnt+beg-j; ++k){
if(ptcon[newmap[k]].ncon > len){
len = ptcon[newmap[k]].ncon;
jtmp = k;
}
}
k = newmap[jtmp];
newmap[jtmp] = newmap[cnt+beg-j-1];
newmap[cnt+beg-j-1] = k;
}
if((cnt-beg) > bw) bw = cnt-beg;
}
start= end;
end = cnt;
/* check to see if the tree leaves have been reached
and if so start a new tree */
if(start == end){
i = 0;
while(!mark[initmap[i]]) ++i;
newmap[cnt] = initmap[i];
mark[initmap[i]]=0;
start = cnt;
end = start+1;
++cnt;
}
}
free(mark);
return bw;
}
static void RenumberElmts(Element *U, Bsystem *Bsys, int *newmap){
register int i,n;
int *gbmap, *length, *sort;
int nsols;
const int nvs = Bsys->nv_solve, nes = Bsys->ne_solve;
#if DIM == 3
const int nfs = Bsys->nf_solve;
#endif
Element *E;
nsols = nvs + nes;
#if DIM==3
nsols += nfs;
#endif
length = ivector(0,nsols-1);
gbmap = ivector(0,nsols);
sort = ivector(0,nsols-1);
/* work out global positions and copy into permanent storage */
ifill(nvs,1,length,1);
for(E=U; E; E = E->next){
for(i = 0; i < E->Nedges; ++i)
if((n = E->edge[i].gid) < nes)
length[nvs+n] = E->edge[i].l;
#if DIM==3
for(i = 0; i < Nface; ++i)
if((n = E[k].face[i].gid)< nfs)
length[nvs+nes+n] = Ek].face[i].l*(1+E[k].face[i].l)/2;
#endif
}
for(i = 0, n = 0; i < nsols; ++i){
sort[i] = n;
n += length[newmap[i]];
}
/* need to unshuffles data to store */
for(i = 0; i < nsols; ++i)
gbmap[newmap[i]] = sort[i];
/* replace new vertex gid's */
for(E=U; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
if((n = E->vert[i].gid) < nvs)
E->vert[i].gid = gbmap[n];
/* replace new edge locations */
icopy(nes, gbmap + nvs, 1, Bsys->edge, 1);
#if DIM==3
/* replave new face locations */
icopy(nfs, gbmap + nvs+nes, 1, Bsys->face, 1);
#endif
free(length); free(sort); free(gbmap);
}
/* setup list of connectivity of all connecting points */
static Fctcon *FacetConnect(Bsystem *B, Element *U){ //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
register int i;
const int nvs = B->nv_solve, nes = B->ne_solve;
int nsols = nvs + nes, *pts,n ;
Fctcon *connect; //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
#if DIM == 3
const int nfs = B->nf_solve;
#endif
Element *E;
#if DIM == 2
pts = ivector(0,7); // max set to nedge + nvert for quad
#else
pts = ivector(0,Nvert + Nedge + Nface-1);
nsols += nfs;
#endif
connect = (Fctcon *)calloc(nsols,sizeof(Fctcon)); //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?); Fctcon: nektar.h(?), nekcomp.h(?)
for(E = U; E; E = E->next){
/* find all facets that are unknowns */
n = 0;
for(i = 0; i < E->Nverts; ++i)
if(E->vert[i].gid < nvs)
pts[n++] = E->vert[i].gid;
for(i = 0; i < E->Nedges; ++i)
if(E->edge[i].gid < nes)
pts[n++] = E->edge[i].gid + nvs;
#if DIM == 3
for(i = 0; i < Nface; ++i)
if(E->face[i].gid < nfs)
pts[n++] = E->face[i].gid + nvs + nes;
#endif
addfct(connect,pts,n);
}
free(pts);
return connect;
}
/* assemble connectivity for all points in list pts and put in connect */
void addfct(Fctcon *con, int *pts, int n){ //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
register int i,j;
int trip;
Facet *f,*f1; //OVERLOAD CALL: Facet: nektar.h(?), nekcomp.h(?)
for(i = 0; i < n; ++i){
f = con[pts[i]].f;
for(j = 0; j < n; ++j){
if(i == j) continue;
/* check to see if this points is already included */
for(trip = 1,f1 = f; f1; f1 = f1->next)
if(f1->id == pts[j]) trip = 0;
/* if this point isn't already in list then add it */ //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
if(trip){
f1 = (Facet *)malloc(sizeof(Facet)); //OVERLOAD CALL: Facet: nektar.h(?), nekcomp.h(?); Facet: nektar.h(?), nekcomp.h(?)
f1->id = pts[j];
f1->next = con[pts[i]].f;
con[pts[i]].f = f1;
con[pts[i]].ncon++;
}
}
}
}
void free_Fctcon(Fctcon *con, int n){ //OVERLOAD CALL: Fctcon: nektar.h(?), nekcomp.h(?)
register int i;
Facet *f,*f1; //OVERLOAD CALL: Facet: nektar.h(?), nekcomp.h(?)
for(i = 0; i < n; ++i)
f = con[i].f;
while(f){
f1 = f->next;
free(f);
f = f1;
}
free(con);
}
C++ to HTML Conversion by ctoohtml