file: Nektar2d/src/prepost.c
/*---------------------------------------------------------------------------*
* RCS Information *
* *
* $Source:
* $Revision:
* $Date:
* $Author:
* $State:
*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include "nektar.h"
/* local functions */
static Bsystem *gen_bsystem(Element_List *E, Bndry *Ebc);
static void Summary (void);
int bandwidth(Element *E, Bsystem *Bsys);
static void Stokes_Pbsys(Bsystem **Pbsys, Bsystem *Ubsys,
Element_List *U, Bndry *Ubc);
struct ssn_tag {
char name[FILENAME_MAX]; /* name of session without postfix */
char fld [FILENAME_MAX]; /* Output (field) file */
char his [FILENAME_MAX]; /* History point file */ //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
char rea [FILENAME_MAX]; /* Input file (from pre-processor) */
#ifdef FORCES
char fce [FILENAME_MAX]; /* Forces file */
#endif
} session;
static struct { /* Default options for Nekscal */
char *name;
int val;
char *descrip;
} nektar_ops[] = {
"order" ,0, "-n # ... run with # of modes ",
"binary" ,1, "-a ... write output in ascii format",
"verbose" ,0, "-v ... plot divergence error",
"variable" ,0, "-var ... run with variable order (uses [].ord file)",
"checkpt" ,0, "-chk ... dump checkpoint fields at \'iostep\' steps",
"iterative" ,0, "-i ... iterative solve",
"Initcond" ,0, "-I ... set up initial condition",
"recursive" ,0, "-r # ... recursive static condesation with # recursions",
"timeavg" ,0, "-T ... evaluate time average field",
"scalar" ,0, "",
0, 0, 0
};
void parse_args(int argc, char *argv[])
{
int i;
char c;
/* initialize the parser and install the defaults */
manager_init();
for (i = 0; nektar_ops[i].name; i++)
option_set(nektar_ops[i].name, nektar_ops[i].val);
if(argc == 1) goto showUsage;
while (--argc && (*++argv)[0] == '-') {
if (strcmp (*argv+1, "chk") == 0) {
option_set("checkpt",1);
continue;
}
else if (strcmp (*argv+1, "var") == 0) {
option_set("variable",1);
continue;
}
while (c = *++argv[0])
switch (c) {
case 'a':
option_set("binary",0);
break;
case 'b':
fprintf(stderr,"Option b is now automatic\n");
break;
case 'B':
option_set("Normalise",1);
break;
case 'c':
{ int n;
if (*++argv[0])
n = atoi (*argv);
else {
n = atoi (*++argv);
argc--;
}
option_set("TORDER",n);
(*argv)[1] = '\0';
break;
}
case 'd':
option_set("FullCG",1);
break;
case 'i':
option_set("iterative",1);
break;
case 'I':
option_set("Initcond",1);
break;
case 'm':
option_set("mixediter",1);
break;
case 'n':
{ int n;
if (*++argv[0])
n = atoi (*argv);
else {
n = atoi (*++argv);
argc--;
}
option_set("NORDER.REQ",n);
(*argv)[1] = '\0';
}
break;
case 'r':
{ int n;
if (*++argv[0])
n = atoi (*argv);
else {
n = atoi (*++argv);
argc--;
}
option_set("recursive",n);
(*argv)[1] = '\0';
}
break;
case 'S':
{
option_set("SLICES",1);
break;
}
case 't':
{ double n;
if (*++argv[0])
n = atof (*argv);
else {
n = atof (*++argv);
argc--;
}
dparam_set("THETA",n);
(*argv)[1] = '\0';
}
break;
case 'T':
option_set("timeavg",1);
break;
case 'v':
option_set("verbose",2);
break;
case 'V':
option_set("tvarying",1);
break;
default:
goto showUsage;
}
}
return;
showUsage:
fputs("usage: nektar [options] file[.rea]\n\n"
"options:\n", stderr);
for (i = 0; nektar_ops[i].name; i++)
fprintf(stderr, "%s\n", nektar_ops[i].descrip);
exit(-1);
}
/* ------------------------------------------------------------------------ *
* PreProcess() - Create a new domain *
* *
* This function creates a fully independent problem domain. *
* *
* ------------------------------------------------------------------------ */
void CorrectRobinFluxes(Bndry *B);
static void Set_Oseen(FILE *rea_file,Bsystem *PB, Bsystem *B, Element_List *U);
Domain *PreProcess(int argc, char **argv) //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
{
FILE *rea_file;
Element_List *U,*P;
Bndry *Ubc,*Pbc;
Bsystem *Ubsys,*Vbsys,*Pbsys;
Domain *omega; //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
ACTION Eqtype; //OVERLOAD CALL: ACTION: nekcomp.h(?), nektar.h(?)
SLVTYPE slvtype;
int i,k,Je;
double Re,dt;
Element *E;
double kinvis;
/* Create the new domain and open the input files */
parse_args(argc, argv);
omega = (Domain*) malloc(sizeof(Domain)); //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?); Domain: nekcomp.h(?), nektar.h(?)
sprintf(session.name,"%s" , strtok(argv[argc-1],"."));
sprintf(session.rea, "%s.rea" , argv[argc-1]);
sprintf(session.fld, "%s.fld" , argv[argc-1]);
sprintf(session.his, "%s.his" , argv[argc-1]);
#ifdef FORCES
sprintf(session.fce, "%s.fce" , argv[argc-1]);
#endif
#ifdef DEBUG
debug_out = stdout;
#endif
omega->name = argv[argc-1];
rea_file = fopen(session.rea,"r");
if (rea_file == (FILE*) NULL)
error_msg(PreProcess(): Could not open input file(s)) //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
else
printf ("Reading from %s\n\n", session.rea);
omega->fld_file = fopen(session.fld,"w");
omega->his_file = fopen(session.his,"w");
#ifdef FORCES
omega->fce_file = fopen(session.fce,"w");
#endif
/* Read the input parameters */
ReadParams (rea_file);
ReadPscals (rea_file);
ReadLogics (rea_file);
Je = iparam("INTYPE");
Eqtype = (ACTION) iparam("EQTYPE"); //OVERLOAD CALL: ACTION: nekcomp.h(?), nektar.h(?)
slvtype = (SLVTYPE) iparam("SLVTYPE");
kinvis = dparam("KINVIS");
if(Eqtype == StokesS)
iparam_set("NSTEPS",1);
else if(Eqtype == Oseen){
iparam_set("NSTEPS",1);
option_set("Oseen",1);
option_set("FAMOFF",1);
}
if(slvtype == StokesSlv){
int i;
int *size = ivector(0, MAXFACETS-1);
Element *E;
/* Build the mesh */
/* set velocity to L+1 */
iparam_set("MODES",iparam("MODES")+1);
U = ReadMesh(rea_file, session.name);
ReadSetLink(rea_file,U);
Summary();
init_ortho_basis();
/* set up pressure system */
P = U->gen_aux_field('p');
option_set("variable",1);
/* set pressure field to L-1 */
free(P->base_h); free(P->base_hj);
for(E=P->fhead;E;E=E->next){
for(i=0;i<E->Nedges;++i)
size[i] = max(0,E->edge[i].l-2);
size[i] = max(0,E->face[0].l-2);
E->Mem_J(size, 'n'); //OVERLOAD CALL: Mem_J: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
E->dgL = E->lmax;
}
P->Cat_mem();
#if 0 /* if not turned on then have Q_P+1 Q_P-1 bases */
/* set edges of velocity to L */
free(U->base_h); free(U->base_hj);
for(E=U->fhead;E;E=E->next){
for(i=0;i<E->Nedges;++i)
size[i] = max(0,E->edge[i].l-1);
size[i] = max(0,E->face[0].l);
E->Mem_J(size, 'n'); //OVERLOAD CALL: Mem_J: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
}
U->Cat_mem();
#endif
free(size);
/* set up velocity system using U */
Ubc = ReadBCs(rea_file,U->fhead); /* also sets up edge links system */
Ubsys = gen_bsystem(U,Ubc); /* set up velocity matrix structure */
Stokes_Pbsys(&Pbsys,Ubsys,U,Ubc);
Pbsys->lambda = (Metric *) calloc(P->nel, sizeof(Metric));
Ubsys->lambda = (Metric *) calloc(U->nel, sizeof(Metric));
if(Eqtype == Oseen)
Set_Oseen(rea_file,Pbsys,Ubsys,U);
}
else{
/* Build the mesh */
U = ReadMesh(rea_file, session.name);
ReadSetLink(rea_file,U);
Summary();
/* set up velocity system using U */
Ubc = ReadBCs (rea_file,U->fhead); /* also sets up edge links system */
Ubsys = gen_bsystem(U,Ubc); /* set up velocity matrix structure */
/* set up pressure system */
P = U->gen_aux_field('p');
Pbc = BuildPBCs (P,Ubc);
Pbsys = gen_bsystem (P,Pbc);
}
omega->Pf = P->gen_aux_field('p');
// Set up multistep storage
omega->u = (double**) malloc(Je*sizeof(double*));
for(k = 0; k < Je; ++k){
omega->u[k] = dvector(0,U->htot*U->nz-1);
dzero(U->htot*U->nz, omega->u[k], 1);
}
omega->us = (double**) calloc(1,sizeof(double*));
omega->us[0] = dvector(0, U->hjtot-1);
dzero(U->hjtot, *omega->us, 1);
omega->Uf = U->gen_aux_field('p');
omega->uf = (double**) malloc(Je*sizeof(double*));
for(k = 0; k < Je; ++k){
omega->uf[k] = dvector(0,omega->Uf->htot*omega->Uf->nz-1);
dzero(U->htot*U->nz, omega->uf[k], 1);
}
omega->V = U->gen_aux_field ('v');
if(option("REFLECT"))
for(E=omega->V->fhead;E;E=E->next)
for(i=0;i<E->Nverts;++i)
E->vert[i].solve = 1;
omega->Vbc = ReadBCs (rea_file,omega->V->fhead);
if(option("REFLECT"))
Vbsys = gen_bsystem(omega->V,omega->Vbc);
else
Vbsys = Ubsys;
omega->Vf = omega->V->gen_aux_field('v');
omega->v = (double**) malloc(Je*sizeof(double*));
for(k = 0; k < Je; ++k){
omega->v[k] = dvector(0,omega->V->htot*omega->V->nz-1);
dzero(U->htot*U->nz, omega->v[k], 1);
}
omega->vs = (double**) calloc(1,sizeof(double*));
omega->vs[0] = dvector(0, U->hjtot-1);
dzero(U->hjtot, *omega->vs, 1);
omega->vf = (double**) malloc(Je*sizeof(double*));
for(k = 0; k < Je; ++k){
omega->vf[k] = dvector(0,omega->Vf->htot*omega->Vf->nz-1);
dzero(U->htot*U->nz, omega->vf[k], 1);
}
omega->U = U; omega->Ubc = Ubc;
omega->P = P; omega->Pbc = Pbc;
dt = dparam("DELT");
omega->soln = NULL;
ReadSoln(rea_file, omega);
int skip = 0;
Re = 1.0 / kinvis;
ReadKinvis (omega);
if(slvtype == StokesSlv){
omega->Pbc = (Bndry *)NULL;
for(k = 0; k < U->nel; ++k)
Ubsys->lambda[k].d = kinvis;
if((Eqtype == StokesS)||(Eqtype == Oseen)){
for(k = 0; k < U->nel; ++k)
Pbsys->lambda[k].d = 0.0;
}
else{
double theta = dparam("THETA");
for(k = 0; k < U->nel; ++k)
Pbsys->lambda[k].d = 1.0/dt/kinvis/(1.0-theta);
}
fprintf(stdout,"Generating stokes system [."); fflush(stdout);
GenMat_Stokes (U,P,Ubsys,Pbsys,Pbsys->lambda);
fprintf(stdout,"]\n"); fflush(stdout);
}
else{
Pbsys->lambda = (Metric *) calloc(P->nel, sizeof(Metric));
fprintf(stdout,"Generating pressure system [."); fflush(stdout);
GenMat (P,Pbc,Pbsys,Pbsys->lambda,Helm);
fprintf(stdout,"]\n"); fflush(stdout);
Ubsys->lambda = (Metric *) calloc(U->nel, sizeof(Metric));
if(omega->kinvis->p){
Ubsys->lambda->p = dvector(0, U->htot-1);
dcopy( U->htot, omega->kinvis->p,1,Ubsys->lambda->p, 1);
}
skip = 0;
for(k=0;k<U->nel;++k){
Ubsys->lambda[k].d = Re*getgamma(1)/dt; //OVERLOAD CALL: getgamma: Integrate.c(?), Integrate.c(?)
if(omega->kinvis[k].p)
Ubsys->lambda[k].p = Ubsys->lambda->p + skip;
skip += U->flist[k]->qtot;
}
fprintf(stdout,"Generating velocity system [."); fflush(stdout);
GenMat (U,Ubc,Ubsys,Ubsys->lambda,Helm);
fprintf(stdout,"]\n"); fflush(stdout);
if(option("REFLECT")){
Vbsys->lambda = Ubsys->lambda;
fprintf(stdout,"Generating velocity system (y)[."); fflush(stdout);
GenMat(omega->V,omega->Vbc,Vbsys,Vbsys->lambda,Helm);
}
}
omega->Usys = Ubsys;
omega->Vsys = Vbsys;
omega->Pressure_sys = Pbsys;
if(!option("tvarying")){
if(slvtype == StokesSlv){
Element_List *V[DIM+1];
Bndry *Vbc[DIM+1];
Bsystem *Vsys[DIM+1];
V[0] = U; V[1] = omega->V; V[2] = omega->P;
Vbc[0] = Ubc; Vbc[1] = omega->Vbc; Vbc[2] = omega->Pbc;
Vsys[0] = Ubsys; Vsys[1] = Vbsys; Vsys[2] = Pbsys;
DirBCs_Stokes(V,Vbc,Vsys);
}
else{
DirBCs(U,Ubc,Ubsys,Helm);
DirBCs(omega->V,omega->Vbc,omega->Vsys,Helm);
U->zerofield();
omega->V->zerofield();
}
}
ReadICs (rea_file, omega);
ReadDF (rea_file, DIM);
ReadHisData (rea_file, omega);
ReadDFunc (rea_file, omega);
fclose(rea_file);
return omega;
}
void PostProcess(Domain *omega, int step, double time){ //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
register int i;
Element_List **V;
Element *E, *F;
FILE *fld_file = omega->fld_file;
int nfields;
int *size = ivector(0, MAXFACETS-1);
nfields = DIM+1;
V = (Element_List**) malloc(nfields*sizeof(Element_List*));
V[0] = omega->U;
V[1] = omega->V;
V[2] = omega->Uf;
#ifndef PCONTBASE
if((SLVTYPE) iparam("SLVTYPE") == StokesSlv)
for(E=omega->P->fhead;E;E=E->next){
E->Obwd(E->vert->hj,*E->h,E->dgL); //OVERLOAD CALL: Obwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
E->Jfwd(E); //OVERLOAD CALL: Jfwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
}
#endif
for(E=omega->P->fhead;E;E=E->next){
for(i=0;i<E->Nedges;++i)
size[i] = E->edge[i].l;
size[i] = E->face[0].l;
F = V[2]->flist[E->id];
F->Copy_field(E->vert[0].hj, size); //OVERLOAD CALL: Copy_field: Fieldfiles.c(Tri), Fieldfiles.c(Quad), Fieldfiles.c(Tet), Fieldfiles.c(Pyr), Fieldfiles.c(Prism), Fieldfiles.c(Hex), Fieldfiles.c(Element)
}
V[2]->Set_state('t');
Writefield(fld_file, omega->name, step, time, nfields, V);
free(size);
fclose(fld_file);
fclose(omega->his_file);
#ifdef FORCES
fclose(omega->fce_file);
#endif
}
/* ----------------------------------------------------------------------- *
* Summary() - Print a summary of the input data *
* *
* This function collects the echo of the input data which used to be *
* scattered throughout this file. *
* *
* ----------------------------------------------------------------------- */
static void Summary (void)
{
printf ("Input File : %s\n",session.name);
printf ("Reynolds number : %g\n",1./dparam("KINVIS"));
#ifdef OSSCYL
printf ("KC number : %g\n",dparam("KC"));
printf ("beta number : %g\n",1./(dparam("KC")*dparam("KINVIS")));
#endif
printf ("Time step : %g\n",dparam("DELT"));
printf ("Integration order : %d\n",iparam("INTYPE"));
if(!option("variable"))
printf("Number of modes : %d\n",iparam("MODES"));
else
printf("Number of modes : variable\n");
printf ("Number of elements : %d\n",iparam("ELEMENTS"));
printf ("Number of Families : %d\n",iparam("FAMILIES"));
if(option("iterative")){
printf("Solver Type : Iterative ");
switch((int) dparam("PRECON")){
case Pre_Diag:
printf("(Precon = diagonal)\n");
break;
case Pre_Block:
printf("(Precon = block)\n");
break;
case Pre_None:
printf("(Precon = none)\n");
break;
case Pre_Diag_Stokes:
printf("(Precon = Stokes System Diagonal)\n");
break;
case Pre_Block_Stokes:
printf("(Precon = Stokes System Block)\n");
break;
}
}
else{
printf("Solver Type : Direct ");
if(iparam("SLVTYPE") == Splitting)
printf("Splitting scheme ");
else
printf("Full Stokes Solver ");
if(option("recursive"))
printf("(Mutlilevel) \n");
else
printf("(RCM) \n");
}
fputs ("Equation type : ", stdout);
switch (iparam("EQTYPE")) {
case Rotational:
puts("Navier-Stokes (rotational)");
break;
case Convective:
puts("Navier-Stokes (convective)");
break;
case Stokes:
puts("Stokes flow");
break;
case StokesS:
puts("Steady Stokes flow");
break;
case Oseen:
puts("Steady Oseen flow");
break;
default:
puts("undefined");
break;
}
printf("Integration time : %g, or %d steps\n",
dparam("FINTIME"), iparam("NSTEPS"));
printf("I/O time for saves : %g, or %d steps",
dparam("IOTIME"), iparam("IOSTEP"));
if (option("checkpt"))
fputs (" [checkpoint]\n", stdout);
else {
putchar ('\n');
if(iparam("NSTEPS"))
if (iparam("NSTEPS") / iparam("IOSTEP") > 10)
fputs ("Summary: " "You have more than 10 dumps..."
"did you want to use -chk?\n", stderr);
}
if(dparam("SPONLEN")){
printf("Sponge position : %g\n", dparam("SPONPOS"));
printf("Sponge length : %g\n", dparam("SPONLEN"));
printf("Sponge magnitude : %g\n", dparam("SPONMAG"));
}
return;
}
/* This is a function to sort out the global numbering scheme for *
* vertices, edges and faces. It also sets up the list of cumulative *
* indices for edges and vertices. */
static void set_bmap (Element *, Bsystem *); //OVERLOAD CALL: set_bmap: prepost.c(?), prepost.c(?)
static int suminterior (Element *E);
static void BasicNumScheme (Element_List *E, Bndry *Ebc, Bsystem *Bsys);
static Bsystem *gen_bsystem(Element_List *UL, Bndry *Ebc){
Bsystem *Bsys;
Element *E = UL->fhead;
Bsys = (Bsystem *)calloc(1,sizeof(Bsystem));
Bsys->Gmat = (MatSys *)malloc( sizeof(MatSys) );
if(option("iterative")){
Bsys->smeth = iterative;
int pt = (int) dparam("PRECON");
Bsys->Precon = (PreType) pt; // dparam("PRECON");
}
/* basis numbering scheme for vertices, edges and faces */
BasicNumScheme(UL,Ebc,Bsys);
if(option("recursive"))
/* Recursive_SC_decom(E,Bsys);*/
Mlevel_SC_decom(UL,Bsys);
else if(Bsys->smeth == direct){
bandwidthopt(E,Bsys);
fprintf(stdout,"rcm bandwidth (%c) : %d [%d (%d)] \n",E->type,
bandwidth(E,Bsys),Bsys->nsolve,suminterior(E));
fflush(stdout);
}
/* set up boundary maps for element edges */
set_bmap(E,Bsys); //OVERLOAD CALL: set_bmap: prepost.c(?), prepost.c(?)
Bsys->families = iparam("FAMILIES");
return Bsys;
}
static void setGid(Element_List *E);
static void BasicNumScheme(Element_List *UL, Bndry *Ebc, Bsystem *Bsys){
register int i;
int nvg,nvs,neg,nes,scnt,ncnt;
int *gsolve,*gbmap,l,l1;
Bndry *Ubc;
Vert *v;
Edge *e;
Element *E, *U=UL->fhead;
setGid (UL); /* setup a global numbering scheme i.e. without boundaries*/
/* This part of the routine re-orders the vertices, edges and then
the faces so that the knowns are listed first. For the edges and
the faces it also initialises a cummalative list stored in Bsys. */
/*--------------------*/
/* Vertex re-ordering */
/*--------------------*/
/* find maximum number of global vertices; */
for(E=U,nvg=0; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
nvg = (E->vert[i].gid > nvg)? E->vert[i].gid:nvg;
++nvg;
gsolve = ivector(0,nvg-1);
gbmap = ivector(0,nvg-1);
ifill(nvg, 1, gsolve,1);
/* Assemble vertex solve mask to sort out multiplicity */
for(E=U; E; E = E->next){
for(i = 0; i < E->Nverts; ++i)
gsolve[E->vert[i].gid] &= E->vert[i].solve;
}
/* copy back mask */
for(E=U; E; E = E->next){
v = E->vert;
for(i = 0; i < E->Nverts; ++i)
v[i].solve = gsolve[E->vert[i].gid];
}
scnt = 0; ncnt = 0;
for(i = 0; i < nvg; ++i)
gbmap[i] = gsolve[i]? scnt++:ncnt++;
nvs = scnt;
/* place unknowns at the end of the pile */
for(i = 0; i < nvg; i++)
gbmap[i] += gsolve[i]? 0:nvs;
/* replace vertices numbering into vertex structures */
for(E=U; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
E->vert[i].gid = gbmap[E->vert[i].gid];
free(gsolve); free(gbmap);
/*------------------*/
/* Edge re-ordering */
/*------------------*/
/* find maximum number of global edges; */
for(E = U,neg=0; E; E = E->next)
for(i = 0; i < E->Nedges; ++i)
neg = (E->edge[i].gid > neg)? E->edge[i].gid:neg;
++neg;
gsolve = ivector(0,neg-1);
gbmap = ivector(0,neg-1);
/* form edge gsolve and gbmap */
ifill(neg, 1, gsolve,1);
for(Ubc = Ebc; Ubc; Ubc = Ubc->next)
if((Ubc->type == 'V')||(Ubc->type == 'W')||(Ubc->type == 'o')){
l = Ubc->face;
gsolve[Ubc->elmt->edge[l].gid] = 0;
}
scnt = 0; ncnt = 0;
for(i = 0; i < neg; ++i)
gbmap[i] = gsolve[i]? (scnt++):(ncnt++);
nes = scnt;
/* place unknowns at the end of the pile */
for(i = 0; i < neg; ++i)
gbmap[i] += gsolve[i]? 0:nes;
/* replace sort gid's */
for(E = U; E; E = E->next){
e = E->edge;
for(i = 0; i < E->Nedges; ++i)
e[i].gid = gbmap[E->edge[i].gid];
}
/* set up cumulative edge list */
Bsys->edge = ivector(0,neg);
for(E = U; E; E = E->next)
for(i = 0; i < E->Nedges; ++i)
Bsys->edge[E->edge[i].gid] = E->edge[i].l;
l = Bsys->edge[0];
Bsys->edge[0] = 0;
for(i = 1; i < neg; ++i){
l1 = Bsys->edge[i];
Bsys->edge[i] = Bsys->edge[i-1]+l;
l = l1;
}
Bsys->edge[i] = Bsys->edge[i-1]+l;
free(gsolve); free(gbmap);
Bsys->nv_solve = nvs;
Bsys->ne_solve = nes;
Bsys->nel = countelements(U);
Bsys->nsolve = nvs + Bsys->edge[nes];
Bsys->nglobal = nvg + Bsys->edge[neg];
Bsys->nf_solve = 0;
/* make numbering scheme consequative the unknowns are listed by
vertices, edges and then faces */
/* place known degree of freedom at end of list */
for(E = U; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
if(E->vert[i].gid >= nvs)
E->vert[i].gid += Bsys->nsolve-nvs;
for(i = 0; i < nes; ++i)
Bsys->edge[i] += nvs;
ncnt = Bsys->nsolve + nvg-nvs - Bsys->edge[nes];
for(i = nes; i <= neg; ++i)
Bsys->edge[i] += ncnt;
}
typedef struct vertnum {
int id;
struct vertnum *base;
struct vertnum *link;
} Vertnum;
int tri_Vnum [][2] = {{0,1},{1,2},{2,0}};
int tri_Vnum1 [][2] = {{1,0},{2,1},{0,2}};
int quad_Vnum [][2] = {{0,1},{1,2},{2,3},{3,0}};
int quad_Vnum1 [][2] = {{1,0},{2,1},{3,2},{0,3}};
static void setGid(Element_List *UL){
register int i,j;
int nvg, face, eid, edgeid, faceid, edgmax;
const int nel = UL->nel;
Edge *e,*ed;
Face *f;
Vertnum *V,*vb,*v;
Element *E, *U=UL->fhead;
/* set vector of consequative numbers */ //OVERLOAD CALL: vector: Coords.c(?), Curvi.c(?)
/* set up vertex list */
edgmax = 4;
V = (Vertnum *) calloc(edgmax*nel,sizeof(Vertnum));
for(E = U; E; E = E->next)
for(i = 0; i < E->Nedges; ++i){
for(j = 0; j < E->Nfverts(i); ++j){ //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
if(E->Nedges == 3)
v = V + E->id*edgmax + tri_Vnum[i][j];
if(E->Nedges == 4)
v = V + E->id*edgmax + quad_Vnum[i][j];
if(E->edge[i].base){
if(E->edge[i].link){
eid = E->edge[i].link->eid;
face = E->edge[i].link->id;
}
else{
eid = E->edge[i].base->eid;
face = E->edge[i].base->id;
}
if(UL->flist[eid]->Nedges == 3)
vb = V[eid*edgmax + tri_Vnum1[face][j]].base;
if(UL->flist[eid]->Nedges == 4)
vb = V[eid*edgmax + quad_Vnum1[face][j]].base;
if(eid < E->id){ /* connect to lower element */
if(!v->base) v->base = v;
/* search through all points and assign to same base */
for(;vb->link;vb = vb->link);
if(vb->base != v->base) vb->link = v->base;
for(v = v->base;v; v = v->link) v->base = vb->base;
}
else if(!v->base) v->base = v;
}
else if(!v->base) v->base = v;
}
}
/* number vertices consequatively */
for(E = U, nvg = 1; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
if(!V[E->id*edgmax+i].base->id)
V[E->id*edgmax+i].id = V[E->id*edgmax+i].base->id = nvg++;
else
V[E->id*edgmax+i].id = V[E->id*edgmax+i].base->id;
nvg--;
for(E = U; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
E->vert[i].gid = V[E->id*edgmax+i].id-1;
/* set gid's to -1 */
for(E = U; E; E = E->next){
for(i = 0; i < E->Nedges; ++i) E->edge[i].gid = -1;
for(i = 0; i < E->Nfaces; ++i) E->face[i].gid = -1;
}
/* at present just number edge and faces consequatively */
faceid = 0; edgeid = 0;
for(E = U; E; E = E->next){
e = E->edge;
for(i = 0; i < E->Nedges; ++i){
if(e[i].gid==-1)
if(e[i].base){
for(ed = e[i].base; ed; ed = ed->link)
ed->gid = edgeid;
++edgeid;
}
else
e[i].gid = edgeid++;
}
f = E->face;
for(i = 0; i < E->Nfaces; ++i)
if(f[i].gid==-1)
if(f[i].link){
f[i].gid = faceid;
f[i].link->gid = faceid++;
}
else
f[i].gid = faceid++;
}
free(V);
return;
}
static int suminterior(Element *E){
int sum=0;
for(;E;E = E->next)
sum += E->Nmodes - E->Nbmodes;
return sum;
}
#ifdef NODAL
static void set_bmap(Element *U, Bsystem *B){
register int i,j,k,n;
int l;
const int nel = B->nel;
int **bmap;
Element *E;
int con_edge_id, con[4][4] = {{1,1,0,0},{1,1,0,0},{0,0,1,1},{0,0,1,1}};
int con_edge_eid;
/* declare memory */
bmap = (int **) malloc(nel*sizeof(int *));
for(E=U,l=0;E;E=E->next) l += E->Nbmodes;
bmap[0] = ivector(0,l-1);
for(i = 0, E=U; i < nel-1; ++i, E=E->next)
bmap[i+1] = bmap[i] + E->Nbmodes;
/* fill with bmaps */
for(E=U;E;E=E->next){
for(j = 0; j < E->Nverts; ++j)
bmap[E->id][j] = E->vert[j].gid;
for(j = 0,n = E->Nverts; j < E->Nedges; ++j,n+=l){
l = E->edge[j].l;
if(E->edge[j].base){
if(E->edge[j].link){
con_edge_id = E->edge[j].link->id;
con_edge_eid = E->edge[j].link->eid;
}
else{
con_edge_id = E->edge[j].base->id;
con_edge_eid = E->edge[j].base->eid;
}
if((E->id > con_edge_eid)&&(con[j][con_edge_id]))
for(k = 0; k < l; ++k)
bmap[E->id][n+k] = B->edge[E->edge[j].gid] + l-1-k;
else
for(k = 0; k < l; ++k)
bmap[E->id][n+k] = B->edge[E->edge[j].gid] + k;
}
else
for(k = 0; k < l; ++k)
bmap[E->id][n+k] = B->edge[E->edge[j].gid] + k;
}
}
B->bmap = bmap;
}
#else
static void set_bmap(Element *U, Bsystem *B){
register int i,j,k,n;
int l;
const int nel = B->nel;
int **bmap;
Element *E;
/* declare memory */
bmap = (int **) malloc(nel*sizeof(int *));
for(E=U,l=0;E;E=E->next) l += E->Nbmodes;
bmap[0] = ivector(0,l-1);
for(i = 0, E=U; i < nel-1; ++i, E=E->next)
bmap[i+1] = bmap[i] + E->Nbmodes;
/* fill with bmaps */
for(E=U;E;E=E->next){
for(j = 0; j < E->Nverts; ++j)
bmap[E->id][j] = E->vert[j].gid;
for(j = 0,n = E->Nverts; j < E->Nedges; ++j,n+=l){
l = E->edge[j].l;
for(k = 0; k < l; ++k)
bmap[E->id][n+k] = B->edge[E->edge[j].gid] + k;
}
}
B->bmap = bmap;
}
static void Stokes_Pbsys(Bsystem **Sys, Bsystem *Ubsys, Element_List *U,
Bndry *Ubc){
register int i,j,d;
int l,cnt;
const int nel = Ubsys->nel;
int **bmap;
int nsolve = Ubsys->nsolve;
Element *E;
Bsystem *Pbsys;
Sys[0] = Pbsys = (Bsystem *)calloc(1,sizeof(Bsystem));
Pbsys->Gmat = (MatSys *)malloc( sizeof(MatSys) );
Pbsys->nsolve = DIM*Ubsys->nsolve + nel;
Pbsys->nglobal = DIM*Ubsys->nglobal + nel;
Pbsys->nel = Ubsys->nel;
Pbsys->smeth = Ubsys->smeth;
Pbsys->Precon = Ubsys->Precon;
/* declare memory */
bmap = (int **) malloc(nel*sizeof(int *));
for(E=U->fhead,l=0;E;E=E->next) l += E->Nbmodes;
l = DIM*l + nel;
bmap[0] = ivector(0,l-1);
for(i = 0, E=U->fhead; i < nel-1; ++i, E=E->next)
bmap[i+1] = bmap[i] + DIM*E->Nbmodes + 1;
/* fill with bmaps */
for(E=U->fhead;E;E=E->next){
for(i = 0; i < E->Nbmodes; ++i){
if(Ubsys->bmap[E->id][i] < nsolve)
for(j = 0; j < DIM; ++j)
bmap[E->id][i+j*E->Nbmodes] = DIM*Ubsys->bmap[E->id][i] + j;
else
for(j = 0; j < DIM; ++j)
bmap[E->id][i+j*E->Nbmodes] = DIM*(Ubsys->bmap[E->id][i] - nsolve) +
j + Pbsys->nsolve;
}
/* reset edges so that u velocity is ordered first then v then w */
cnt = E->Nverts;
for(i = 0; i < E->Nedges; ++i){
l = E->edge[i].l;
if(E->edge[i].gid < Ubsys->ne_solve){
for(d = 0; d < DIM; ++d)
for(j = 0; j < l; ++j)
bmap[E->id][cnt + j + d*E->Nbmodes] = bmap[E->id][cnt] + j + d*l;
}
cnt += l;
}
/* add pressure dof */
bmap[E->id][DIM*E->Nbmodes] = DIM*Ubsys->nsolve + E->id;
}
/* set pressure singular mode to be value in element iparam("IESING")*/
Pbsys->singular = bmap[iparam("IESING")][DIM*U->fhead->Nbmodes]+1;
Pbsys->bmap = bmap;
/* reset recursive numbering system and set up in Pbsys */
if(Ubsys->rslv){
register int k,n;
Rsolver *Vrslv = Ubsys->rslv;
Rsolver *Prslv;
int nrecur = Vrslv->nrecur;
Recur *Vrdat = Vrslv->rdata;
Recur *Prdat;
int len,cstart,top,alen,maxid,eid,maplen,*nelmt,*elmtid,*nelmt1,*elmtid1;
int *mapping,*nptchold,*n2optch, *swap, *swap1;
elmtid = ivector(0,nel-1);
nelmt = ivector(0,nel);
elmtid1 = ivector(0,nel-1);
nelmt1 = ivector(0,Vrdat[0].npatch);
n2optch = ivector(0,nel-1);
nptchold = ivector(0,Vrdat[0].npatch);
mapping = ivector(0,Pbsys->nsolve-1);
/* set up elmt numbering for original patch i.e. individual elements */
top = 0; eid = 1;
iramp(nel+1,&top,&eid,nelmt,1);
iramp(nel ,&top,&eid,elmtid,1);
Pbsys->rslv = Prslv = (Rsolver *)calloc(1,sizeof(Rsolver));
Prslv->nrecur = nrecur;
Prdat = Prslv->rdata = (Recur *)calloc(nrecur,sizeof(Recur));
/* declare rdata in Pbsys for new system */
for(i = 0; i < nrecur; ++i){
ifill(Pbsys->nsolve,-1,mapping,1);
Prdat[i].npatch = Vrdat[i].npatch;
Prdat[i].patchlen_a = ivector(0,Prdat[i].npatch-1);
Prdat[i].patchlen_c = ivector(0,Prdat[i].npatch-1);
Prdat[i].map = (int **)malloc(Vrdat[i].npatch*sizeof(int *));
Prdat[i].pmap = Vrdat[i].pmap;
if(i){
top = Prdat[i-1].cstart;
maplen = Prdat[i-1].npatch;
}
else{
top = Pbsys->nsolve;
maplen = nel;
}
/* determine new to old mapping of patches */
nptchold[0] = 0;
for(j = 0; j < Vrdat[i].npatch; ++j){
for(k = 0,cnt=0; k < maplen; ++k)
if(Vrdat[i].pmap[k] == j)
n2optch[nptchold[j]+ cnt++] = k;
nptchold[j+1] = nptchold[j] + cnt;
}
nelmt1[0] = 0;
for(j = 0; j < Vrdat[i].npatch; ++j){
nelmt1[j+1] = nelmt1[j];
for(k=nptchold[j]; k < nptchold[j+1]; ++k){
len = nelmt[n2optch[k]+1]-nelmt[n2optch[k]];
icopy(len,elmtid+nelmt[n2optch[k]],1,elmtid1+nelmt1[j+1],1);
nelmt1[j+1] += len;
}
}
swap = nelmt;
swap1 = elmtid;
nelmt = nelmt1;
elmtid = elmtid1;
nelmt1 = swap;
elmtid1 = swap1;
/* reshuffle bmaps so that pressure dof are lumped with patch blocks */
cstart = DIM*Vrdat[i].cstart;
cnt = Vrdat[i].npatch;
for(j = 0; j < Vrdat[i].npatch; ++j){
for(k = nelmt[j],maxid=0; k < nelmt[j+1]; ++k){
eid = elmtid[k];
alen = DIM*U->flist[eid]->Nbmodes+1;
for(n = 0; n < alen-1; ++n)
if((bmap[eid][n] >= cstart)&&(bmap[eid][n] < top)){
mapping[bmap[eid][n]] = bmap[eid][n] + cnt;
bmap[eid][n] += cnt;
maxid = max(bmap[eid][n],maxid);
}
}
/* Redefine pressure dof putting first element at
beginning of cstart so it is carried over to boundary
solve. Otherwise define the dof to be in the interior
block */
eid = elmtid1[nelmt1[n2optch[nptchold[j]]]];
alen = DIM*U->flist[eid]->Nbmodes+1;
mapping[bmap[eid][alen-1]] = cstart+j;
bmap[eid][alen-1] = cstart+j;
for(k = nptchold[j]+1,n=1; k < nptchold[j+1]; ++k,++n){
eid = elmtid1[nelmt1[n2optch[k]]];
alen = DIM*U->flist[eid]->Nbmodes+1;
mapping[bmap[eid][alen-1]] = maxid+n;
bmap[eid][alen-1] = maxid+n;
}
len = nptchold[j+1]-nptchold[j];
cnt += len-1;
Prdat[i].patchlen_a[j] = DIM*Vrdat[i].patchlen_a[j] + 1;
Prslv->max_asize = max(Prslv->max_asize,Prdat[i].patchlen_a[j]);
Prdat[i].map[j] = ivector(0,Prdat[i].patchlen_a[j]-1);
Prdat[i].patchlen_c[j] = DIM*Vrdat[i].patchlen_c[j] + len-1;
}
Prdat[i].cstart = cstart + Prdat[i].npatch;
/* fill out mapping for each patch to new local A system */
cstart = Prdat[i].cstart;
/* go through and redefine previous Prdat[i].map which will have
been altered due to new recursion */
for(j = 0; j < i; ++j){
for(k = 0; k < Prdat[j].npatch; ++k)
for(n = 0; n < Prdat[j].patchlen_a[k]; ++n)
if(mapping[Prdat[j].map[k][n]]+1)
Prdat[j].map[k][n] = mapping[Prdat[j].map[k][n]];
}
for(j = 0; j < Prdat[i].npatch; ++j){
ifill(cstart,-1,mapping,1);
for(k = nelmt[j],maxid=0; k < nelmt[j+1]; ++k){
eid = elmtid[k];
alen = DIM*U->flist[eid]->Nbmodes+1;
for(n = 0; n < alen; ++n)
if(bmap[eid][n] < cstart)
mapping[bmap[eid][n]] = bmap[eid][n];
}
for(k = 0,cnt=0; k < cstart; ++k)
if(mapping[k]+1)
Prdat[i].map[j][cnt++] = mapping[k];
if(cnt != Prdat[i].patchlen_a[j])
error_msg(error in reset_pbmap); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
}
}
/* find the patch id that bmap[IEsing] belongs to */
k = iparam("IESING");
for(i = 0; i < nrecur; ++i)
k = Prdat[i].pmap[k];
/* work out new patch id */
if(k >= Prdat[nrecur-1].npatch)
error_msg(error sorting singular pressure in Pbsys_bmap); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
/* set pressure singular mode */
Pbsys->singular = Prdat[nrecur-1].cstart-k;
free(mapping);
free(nelmt); free(elmtid);
free(nelmt1); free(elmtid1);
free(nptchold); free(n2optch);
}
/* finally check to see if an outflow is defined - if so turn off
singular pressure mode */
for(;Ubc;Ubc = Ubc->next)
if(Ubc->type == 'O'){
Pbsys->singular = 0;
break;
}
}
static void Set_Oseen(FILE *rea_file,Bsystem *PB, Bsystem *B, Element_List *U){
register int i,j;
int eDIM = U->fhead->dim(),qtot; //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
double **wave,*s;
wave = (double **)malloc(eDIM*B->nel*sizeof(double*));
s = dvector(0,eDIM*U->htot);
dzero(eDIM*U->htot,s,1);
for(i = 0; i < B->nel; ++i){
qtot = U->flist[i]->qtot;
for(j = 0; j < eDIM; ++j, s+= qtot)
wave[i*eDIM + j] = s;
B ->lambda[i].wave = wave + i*eDIM;
PB->lambda[i].wave = wave + i*eDIM;
}
ReadWave(rea_file,wave,U);
/* divide through by the viscosity to account for way that
matrix is formed */
dscal(eDIM*U->htot,1.0/dparam("KINVIS"),wave[0],1);
}
#endif
C++ to HTML Conversion by ctoohtml