file: Nektar2d/src/io.c
/*
* Input functions based on the NEKTON spectral element data model
*/
/*---------------------------------------------------------------------------*
* RCS Information *
* *
* $RCSfile: io.C,v $
* $Revision: 1.1 $
* $Author: tcew $
* $Date: 1996/12/11 16:29:57 $
* $State: Exp $
* ------------------------------------------------------------------------- */
// set up for 2d only
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <veclib.h>
#include "nektar.h"
#include "Quad.h"
#include "Tri.h"
/* ..... Section names ..... */
#define SEC_PARAMS sections[0] /* Parameters */
#define SEC_PSCALS sections[1] /* Passive scalars */
#define SEC_LOGICS sections[2] /* Logical switches */
#define SEC_MESH sections[3] /* Mesh definition */
#define SEC_CURVE sections[4] /* Curved sides */
#define SEC_BCS sections[5] /* Boundary conditions */
#define SEC_ICS sections[6] /* Initial conditions */
#define SEC_DF sections[7] /* Drive force */
#define SEC_HIST sections[8] /* History points */
#define SEC_T_BCS sections[9] /* Thermo boundary conditions */
#define SEC_DFUNCS sections[10] /* Drive force */
#define SEC_WAVE sections[11] /* Wave Prop in Oseen eqn */
static int checklist (char *name, char *list[]);
static char *findSection (char *name, char *buf, FILE *fp);
static int bedgcmp (const void *b1, const void *b2);
static char *sections[] = {
"Paramaters", "Passive Scalars", "LOGICAL", "MESH", "Curved sides",
"BOUNDARY", "INITIAL CONDITIONS", "DRIVE FORCE", "HISTORY", "THERMAL",
"DRIVE FUNCT", "WAVE"};
/* ------------------------------------------------------------------------ *
* Spectral Element Data Model *
* *
* This file contains functions to read a spectral element input file in *
* the NEKTON format. The file format is fairly specific, expecting the *
* following sections in order: *
* *
* PARAMETERS -> Passive scalars -> Logical switches -> MESH -> *
* BC's -> INITIAL COND.'s -> DRIVE FORCE -> Variables -> *
* HISTORY POINTS *
* *
* The uppercase sections are the ones which are actually used. In the *
* parameters section, the following should be defined if you are solving *
* Stokes or Navier-Stokes: *
* *
* KINVIS Viscosity, or 1/Re *
* NSTEP | FINTIME Number of steps or simulation time. *
* DELT Time step *
* *
* Return value: none *
* ------------------------------------------------------------------------ */
void ReadParams (FILE *fp)
{
int n;
char buf[BUFSIZ], value[25], name[25], c;
double val, dt;
static char *dspecial[] = { "KINVIS", "LAMBDA", "IOTIME", "KC", "LZ",
"FRNDX","FRNDY","FRNDZ" , "PECLET", 0 };
static char *ispecial[] = { "TEQN", "EQTYPE", "CFLSTEP", "HISSTEP",
"SLVTYPE", 0 };
rewind (fp);
for (n = 0; n < 4; n++) fgets(buf, BUFSIZ, fp);
if (sscanf(buf, "%d", &n) != 1)
{fputs("ReadParams: can't read # of parameters", stderr);exit(-1);}
while (n--) {
fgets (buf, BUFSIZ, fp);
if(sscanf(buf, "%25s%25s", value, name) == 2) {
if (checklist (name, dspecial))
dparam_set (name, scalar(value));
else if (checklist (name, ispecial))
iparam_set (name, scalar(value));
else if (isupper(c = *name) && 'I' <= c && c <= 'N')
iparam_set(name, scalar(value));
else
dparam_set(name, scalar(value));
}
}
/* The following section is only for unsteady problems */
if (!option("scalar")) {
if ((dt=dparam("DELT")) == 0.) {
if ((dt=dparam("DT")) != 0.)
dparam_set("DELT", dt);
else
{fputs("ReadParams: no time step specified", stderr);exit(-1);}
}
/* FINTIME overrides NSTEPS */
if ((val=dparam("FINTIME")) > 0.)
iparam_set("NSTEPS", n = (val+.5*dt)/dt);
else
dparam_set("FINTIME", val = (n=iparam("NSTEPS"))*dt);
/* IOTIME overrides IOSTEP */
n = clamp(iparam("IOSTEP"), 0 , n);
val = clamp(dparam("IOTIME"), 0., val);
if (val == 0. && n == 0) {
n = iparam("NSTEPS");
val = dparam("FINTIME");
} else if (val > 0.)
n = (val + .5 * dt) / dt;
else
val = n * dt;
iparam_set("IOSTEP", n);
dparam_set("IOTIME", val);
}
/* check to see if there is a command line input for norder */
/* set modes to be equal to norder if it is not already set */
if(n = option("NORDER.REQ")) iparam_set("MODES",n);
if(!iparam("MODES") ) iparam_set("MODES",iparam("NORDER"));
return;
}
static int checklist (char *name, char *list[])
{
do
if (strcmp(name,*list) == 0)
return 1;
while
(*++list);
return 0;
}
/* The next two functions are optional */
void ReadPscals (FILE *fp)
{
int n;
char buf[BUFSIZ];
fgets(buf, BUFSIZ, fp);
if (sscanf(buf,"%d", &n) != 1)
{fputs("ReadPscals: no passive scalar data", stderr);exit(-1);}
while (n--) fgets(buf, BUFSIZ, fp);
return;
}
void ReadLogics (FILE *fp)
{
int n;
char buf[BUFSIZ];
fgets(buf, BUFSIZ, fp);
if (sscanf(buf, "%d", &n) != 1 || !strstr(buf, SEC_LOGICS))
{fputs("ReadLogics: no logical switches found\n", stderr); exit(-1);}
while (n--) fgets(buf, BUFSIZ, fp);
return;
}
static void ReadCurve(FILE *fp, Element_List *new_E);
Element_List *ReadMesh(FILE *fp, char *session_name)
{
int nel, L, qa, qb, k;
char buf[BUFSIZ];
Element **new_E;
register int i;
/* set up modes and quadrature points */
if(!( L = iparam("MODES")))
{fputs("ReadMesh: Number of modes not specified\n",stderr);exit(-1);}
/* note quadrature order reset for variable order runs */
if(qa = iparam("LQUAD"));
else qa = L + 1;
if(qb = iparam("MQUAD"));
else qb = L;
if (!findSection (SEC_MESH, buf, fp))
{fputs("ReadMesh: Section not found\n", stderr), exit(-1);}
if (sscanf(fgets(buf,BUFSIZ,fp), "%d%d", &nel,&k) != 2)
{fputs("ReadMesh: unable to get the number of elements\n",stderr);exit(1);}
if (k != DIM)
{fputs("ReadMesh: Read file is wrong dimension\n", stderr);exit(1);}
iparam_set("ELEMENTS", nel);
/* Set up a new element vector */ //OVERLOAD CALL: vector: Coords.c(?), Curvi.c(?)
QGmax = qa;
LGmax = L;
new_E = (Element**) malloc(nel*sizeof(Element*));
Coord X; //OVERLOAD CALL: Coord: nekstruct.h(?), hotel.h(?)
X.x = dvector(0,3);
X.y = dvector(0,3);
int nvs;
/* Read in mesh information */
for(k = 0; k < nel; k++) {
fgets(buf, BUFSIZ, fp); /* element header */
if(strstr(buf,"Qua") || strstr(buf,"qua"))
nvs = 4;
else
nvs = 3;
/* -X- Coordinates */
for (i = 0; i < nvs; i++)
fscanf (fp, "%lf", X.x+i);
fgets(buf, BUFSIZ, fp); /* finish the line */
/* -Y- Coordinates */
for (i = 0; i < nvs; i++)
fscanf (fp, "%lf", X.y+i);
fgets(buf, BUFSIZ, fp); /* finish the line */
#ifndef NODAL
if(nvs == 4)
new_E[k] = new Quad(k,'u', L,qa,qa,0, &X); //OVERLOAD CALL: Quad: Quad.c(Quad), Quad.c(Quad), Quad.c(Quad), Quad.c(Quad)
else
new_E[k] = new Tri(k,'u', L,qa,qb,0, &X); //OVERLOAD CALL: Tri: Tri.c(Tri), Tri.c(Tri), Tri.c(Tri), Tri.c(Tri)
#else
if(nvs == 4)
new_E[k] = new Nodal_Quad(k,'u', L,qa,qa,0, &X); //OVERLOAD CALL: Quad: Quad.c(Quad), Quad.c(Quad), Quad.c(Quad), Quad.c(Quad)
else
new_E[k] = new Nodal_Tri (k,'u', L,qa,qb,0, &X); //OVERLOAD CALL: Tri: Tri.c(Tri), Tri.c(Tri), Tri.c(Tri), Tri.c(Tri)
#endif
}
for(k = 0; k < nel-1; ++k) new_E[k]->next = new_E[k+1];
new_E[k]->next = (Element*) NULL;
Element_List* E_List = new Element_List(new_E, nel);
// Need to fix
if(option("variable"))
ReadOrderFile (session_name,E_List);
option_set("NZ",1);
option_set("NZTOT", 1);
E_List->nz = 1;
E_List->nztot = 1;
E_List->Cat_mem();
/* Read the curved side information */
ReadCurve(fp,E_List);
double sponlen = dparam("SPONLEN");
double sponpos = dparam("SPONPOS");
for(k = 0; k < nel; ++k) {
if(new_E[k]->identify() == Nek_Quad && new_E[k]->curvX == NULL) //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
new_E[k]->set_curved_elmt(E_List); //OVERLOAD CALL: set_curved_elmt: Coords.c(Tri), Coords.c(Quad), Coords.c(Tet), Coords.c(Pyr), Coords.c(Prism), Coords.c(Hex), Coords.c(Element)
#if 1
if(sponlen){
for(i=0;i<new_E[k]->Nverts;++i)
if(new_E[k]->vert[i].x > sponpos-1.){
option_set("FAMOFF",1);
break;
}
if(i==new_E[k]->Nverts)
#endif
option_set("FAMOFF",0);
}
new_E[k]->set_geofac(); //OVERLOAD CALL: set_geofac: Geofac.c(Tri), Geofac.c(Quad), Geofac.c(Tet), Geofac.c(Pyr), Geofac.c(Prism), Geofac.c(Hex), Geofac.c(Element)
}
return E_List;
}
static void ReadCurve(FILE *fp, Element_List *new_E){
register int i,j;
int k;
int eid,face;
char type[2], *p, buf[BUFSIZ];
Curve *curve,*ctmp;
i = 2; while (i--) fgets (buf, BUFSIZ, fp);
if(strstr(buf,"type")){
int ntags;
char *tagid;
sscanf(buf, "%d", &ntags);
if(ntags){
ctmp = (Curve *) malloc((ntags+1)*sizeof(Curve));
tagid = (char *) malloc(BUFSIZ*sizeof(char));
/* read in curved info */
for(i = 0; i < ntags; ++i){
fgets (buf, BUFSIZ, fp);
if(strstr(buf,"Str")||strstr(buf,"str")){ /* straight sided */
fgets (buf, BUFSIZ, fp);
sscanf(buf,"%1s",tagid+i);
ctmp[i].type = T_Straight;
}
else if(strstr(buf,"Cir")||strstr(buf,"cir")){ /* circle */
fgets (buf, BUFSIZ, fp);
sscanf (buf, "%*lf%*lf%lf%1s",&(ctmp[i].info.arc.radius),tagid+i);
ctmp[i].type = T_Arc;
}
else if(strstr(buf,"Fil")||strstr(buf,"fil")){ /* spline fit - file */
char file[BUFSIZ];
fgets (buf, BUFSIZ, fp);
sscanf (buf, "%s",file);
ctmp[i].info.file.name = strdup(file);
sscanf (buf, "%*s%1s",tagid+i);
ctmp[i].type = T_File;
}
else if(strstr(buf,"Sin")||strstr(buf,"sin")){
fgets (buf, BUFSIZ, fp);
sscanf (buf, "%lf%lf%lf%lf%1s",&(ctmp[i].info.sin.xo),
&(ctmp[i].info.sin.yo),&(ctmp[i].info.sin.amp),
&(ctmp[i].info.sin.wavelength),tagid+i);
ctmp[i].type = T_Sin;
}
else if(strstr(buf,"Ell")||strstr(buf,"ell")){
fgets (buf, BUFSIZ, fp);
sscanf (buf, "%lf%lf%lf%lf%1s",
&(ctmp[i].info.ellipse.xo),
&(ctmp[i].info.ellipse.yo),
&(ctmp[i].info.ellipse.rmaj),
&(ctmp[i].info.ellipse.rmin),
tagid+i);
ctmp[i].type = T_Ellipse;
}
}
fgets (buf, BUFSIZ, fp);
sscanf(buf, "%d", &k);
if (k) {
if(!ntags) error_msg(No curved type information specified); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
curve = (Curve *) malloc(k*sizeof(Curve));
for(i = 0; i < k; ++i){
if (fscanf(fp, "%d%d%1s", &face, &eid, type) != 3)
error_msg(ReadMesh -- unable to read curved information); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
new_E->flist[eid-1]->curve = curve + i;
for(j = 0; j < ntags; ++j)
if(tagid[j] == type[0]){
curve[i] = ctmp[j];
// memcpy(curve+i,ctmp+j,sizeof(Curve));
break;
}
if(j == ntags){
fprintf(stderr,"can not find tag type `%c\' called in element"
" %d face %d\n" ,type[0],eid,face);
exit(-1);
}
curve[i].face = --face;
new_E->flist[--eid]->set_curved_elmt(new_E); //OVERLOAD CALL: set_curved_elmt: Coords.c(Tri), Coords.c(Quad), Coords.c(Tet), Coords.c(Pyr), Coords.c(Prism), Coords.c(Hex), Coords.c(Element)
}
}
free((char*)tagid);
}
}
else{
sscanf(buf, "%d", &k);
if (k) {
curve = (Curve *) calloc(k,sizeof(Curve));
for(i = 0; i < k; ++i){
fgets(buf, BUFSIZ, fp);
if (sscanf(buf, "%d%d", &face, &eid) != 2)
error_msg(ReadMesh -- unable to read curved information); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
new_E->flist[--eid]->curve = curve + i;
curve[i].face = --face;
p = buf + strlen(buf); while (isspace(*--p));
type[0] = toupper(*p);
switch (type[0]) {
case 'S': /* Straight side */
curve[i].type = T_Straight;
break;
case 'C': /* Circular arc */
sscanf (buf, "%*d%*d%lf",&(curve[i].info.arc.radius));
curve[i].type = T_Arc;
break;
default:
sprintf (buf, "unknown curve type -- %c", type[0]);
fputs(buf,stderr);
exit(-1);
break;
}
}
}
}
}
static char *findSection (char *name, char *buf, FILE *fp)
{
char *p;
while (p = fgets (buf, BUFSIZ, fp))
if (strstr (p, name))
break;
return p;
}
/* ----------------------------------------------------------------------- *
* ReadICs() - Read Initial Conditions *
* *
* This function reads in the initial conditions for the velocity fields *
* from the .rea file. The form of the velocity initial conditions is *
* determined by one of the following keywords in the line following the *
* beginning of the INITIAL CONDITIONS section: *
* *
* Default Initial velocity field is at rest *
* Restart File name to read IC's from is given *
* Given V(x,y,z) (t = 0) specified *
* ----------------------------------------------------------------------- */
static void Restart (char *name, int nfields, Element_List *U[]);
void ReadICs (FILE *fp, Domain *omega) //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
{
register int i;
int n, type, nfields, tot;
char buf[BUFSIZ];
Element_List *U [MAXFIELDS];
const int nel=omega->U->nel;
static char *ictypes[] = { "PrepOnly", "Restart", "Default", "Given",
"Hardwired", "Minion"};
/* initialize the Element pointers */
#if DIM == 2
#ifndef THERMO
U[0] = omega->U; U[1] = omega->V; U[2] = omega->P;
nfields = DIM;
#else
U[0] = omega->U; U[1] = omega->V; U[2] = omega->T; U[3] = omega->P;
nfields = DIM+1;
#endif
#else
U[0] = omega->U; U[1] = omega->V; U[2] = omega->W; U[3] = omega->P;
nfields = DIM;
#endif
/* zero pressure fields */
for(i = 0; i < nel; ++i){
#if DIM == 3
tot = U[3]->flist[i]->qa*U[3]->flist[i]->qb*U[3]->flist[i]->qc;
dzero(tot, *U[3]->flist[i]->h[0], 1);
#else
tot = U[nfields]->flist[i]->qa*U[nfields]->flist[i]->qb;
dzero(tot, *U[nfields]->flist[i]->h, 1);
#endif
dzero(U[nfields]->flist[i]->Nmodes,U[nfields]->flist[i]->vert->hj,1);
}
if (!findSection (SEC_ICS, buf, fp))
{fprintf(stderr,"ReadICs: no initial conditions specified "
"(last line read below) %s", buf); exit(-1);}
if (sscanf(buf, "%d", &n) != 1)
{fprintf(stderr,"ReadICs: # of initial conditions not specified "
"(last line read below) %s", buf);exit(-1);}
if (option("prep")) /* Preprocessor */
type = 0;
else if (n == 0) /* Default */
type = 2;
else { /* Find it... */
n--; fgets (buf, BUFSIZ, fp);
for (i = 1; i < 6; i++)
if (strstr(buf, ictypes[type = i])) break;
}
switch (type) {
case 0:
break; /* Pre-processor only */
case 1: { /* Restart from a file */
fgets (buf, BUFSIZ, fp); n--;
Restart (buf, DIM+1, U);
break;
}
case 2: { /* Default Initial Conditions */
for (i = 0; i < nfields; i++){
U[i]->zerofield();
set_state(U[i]->fhead,'t');
}
printf("Initial solution : Default\n");
break;
}
case 3: { /* Function(x,y,z) Initial conditions */
char *s;
printf("Initial solution :");
for (i = 0; i < nfields && n--; i++) {
if(!(s = strchr(fgets(buf,BUFSIZ,fp),'=')))
{fprintf(stderr,"ReadICs: the following function is invalid %s",s);
exit(-1);}
while (isspace(*++s));
U[i]->Set_field(s); //OVERLOAD CALL: Set_field: Element_List.c(Element_List), Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(!i) printf (" %c = %s", U[i]->fhead->type, s);
else printf (" %c = %s", U[i]->fhead->type, s);
}
break;
}
case 5: {
int j;
double x, y, u, v,
eps = dparam("SHEAR_EPS"),
delta = dparam("SHEAR_DELTA");
Element* E;
Coord X; //OVERLOAD CALL: Coord: nekstruct.h(?), hotel.h(?)
X.x = dvector(0,QGmax*QGmax-1);
X.y = dvector(0,QGmax*QGmax-1);
for(i = 0; i < nel; ++i){
E=U[0]->flist[i];
E->coord(&X); //OVERLOAD CALL: coord: Coords.c(Tet), Coords.c(Pyr), Coords.c(Prism), Coords.c(Hex), Coords.c(Element), Coords.c(Tri), Coords.c(Quad)
for(j = 0; j < E->qa*E->qb; ++j){
x = X.x[j];
y = X.y[j];
if (y < 0.5) u = tanh (eps*(y-0.25));
else u = tanh (eps*(0.75-y));
v = delta * cos (2*M_PI*x);
(*U[0]->flist[i]->h)[j] = u;
(*U[1]->flist[i]->h)[j] = v;
}
}
set_state(U[0]->fhead,'p');
set_state(U[1]->fhead,'p');
printf("Initial solution : Minion\n");
break;
}
default:
fputs("ReadICs: invalid initial conditions", stderr);
exit(-1);
break;
}
while (n--) fgets(buf, BUFSIZ, fp); /* read to the end of the section */
return;
}
static void Restart (char *name, int nfields, Element_List *U[]){
Field f;
FILE *fp;
char *p, fname[FILENAME_MAX], buf[BUFSIZ];
register int i, pos, dump, nel;
if (sscanf(name, "%s", fname) != 1)
{fputs("ReadICs: couldn't read the restart file name", stderr);exit(-1);}
if ((fp = fopen(fname,"r")) == (FILE *) NULL)
{fprintf(stderr,"ReadICs: restart file %s not found",fname); exit(-1);}
nel = U[0]->nel;
dump = 0;
memset(&f, '\0', sizeof (Field));
while (readField (fp, &f, U[0]->fhead)) dump++;
if (!dump)
{fputs("Restart: no dumps read from restart file", stderr);exit(-1);}
if (f.dim != DIM)
{fputs("Restart: file if wrong dimension", stderr); exit(-1);}
if (f.nel != nel)
{fputs("Restart: file is the wrong size", stderr); exit(-1);}
printf("Initial solution : %s at t = %g", fname, f.time);
putc ('\n', stdout);
for (i = 0; i < nfields; i++) {
if (!(p = strchr (f.type, U[i]->fhead->type))) {
sprintf (buf, "variable %c is not on file", U[i]->fhead->type);
fprintf(stderr,"Restart: %s", buf); exit(-1);
}
pos = (int) (p - f.type);
copyfield (&f,pos,U[i]->fhead);
}
dparam_set("STARTIME", f.time);
dparam_set("FINTIME", dparam("FINTIME") + f.time);
freeField (&f);
fclose (fp);
return;
}
/* ----------------------------------------------------------------------- *
* ReadDF() - Read Drive Force data *
* *
* This is really one function designed to handle two different types of *
* drive force specification. If the "steady" flag is active, the drive *
* force can be an arbitary function. Otherwise, it must be constant but *
* can have any direction. *
* *
* If the flowrate option is active, this section is skipped entirely. *
* *
* Example: *
* *
* ***** DRIVE FORCE DATA ***** PRESSURE GRAD, FLOW, Q *
* 4 Lines of Drive force data follow *
* FFX = 0. *
* FFY = 0. *
* FFZ = 2. * KINVIS *
* C *
* *
* The names given to these lines does not matter, but their order is *
* taken as (x,y,z). The '=' MUST be present, and everything to the right *
* of it determines the forcing functions. *
* ----------------------------------------------------------------------- */
void ReadDFunc(FILE *fp, Domain *Omega){ //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
int nlines;
char buf[BUFSIZ], *p;
if (!findSection (SEC_DFUNCS, buf, fp))
{
fputs("ReadDFunc: section not found\n", stderr);
Omega->ForceFuncs = (double**)0;
rewind(fp);
return;
}
if (sscanf(fgets(buf, BUFSIZ, fp), "%d", &nlines) != 1)
{fputs("ReadDFunc: can't read the number of lines", stderr); exit(-1);}
Omega->ForceFuncs = dmatrix(0, DIM-1, 0, Omega->U->htot-1);
if(p = strchr (fgets(buf, BUFSIZ, fp), '=')) {
while (isspace(*++p));
Omega->Uf->Set_field(p); //OVERLOAD CALL: Set_field: Element_List.c(Element_List), Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dcopy(Omega->Uf->htot, Omega->Uf->base_h, 1, Omega->ForceFuncs[0],1);
}
if(p = strchr (fgets(buf, BUFSIZ, fp), '=')) {
while (isspace(*++p));
Omega->Uf->Set_field(p); //OVERLOAD CALL: Set_field: Element_List.c(Element_List), Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dcopy(Omega->Uf->htot, Omega->Uf->base_h, 1, Omega->ForceFuncs[1],1);
}
}
void ReadDF(FILE *fp, int nforces, ...)
{
register int i;
int nlines;
char buf[BUFSIZ], *p;
static char *ftypes[] = { "FFX", "FFY", "FFZ" };
Element *E;
if (!findSection (SEC_DF, buf, fp))
{fputs("ReadDF: section not found", stderr); exit(-1);}
if (sscanf(fgets(buf, BUFSIZ, fp), "%d", &nlines) != 1)
{fputs("ReadDF: can't read the number of lines", stderr); exit(-1);}
/* Don't process this section if running in pre-processor mode */
if (option("prep")) {
while (nlines--) fgets(buf, BUFSIZ, fp);
return;
}
/* Check to see if the problem is steady scalar or vector */ //OVERLOAD CALL: vector: Coords.c(?), Curvi.c(?)
if (option("scalar")) {
va_list ap;
Element_List *U [MAXFIELDS];
va_start (ap, nforces);
for (i = 0; i < nforces; i++) U[i] = va_arg(ap, Element_List *);
va_end (ap);
fputs ("Forcing function : ", stdout);
switch (nlines ? iparam("EQTYPE") : Laplace) {
case Poisson: case Helmholtz: default:
for (i = 0; i < nforces && nlines--; i++) {
if (p = strchr (fgets(buf, BUFSIZ, fp), '=')) {
while (isspace(*++p));
if (i) fputs (" ", stdout);
printf ("%c = %s", U[i]->fhead->type, p);
U[i]->Set_field(p); //OVERLOAD CALL: Set_field: Element_List.c(Element_List), Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
}
else
fprintf (stderr,"ReadDF -- invalid function (below) %s", buf);
}
break;
case Laplace:
puts ("none");
for (i = 0; i < nforces; i++){
for(E=U[i]->fhead;E;E=E->next)
dzero(E->qa * E->qb , E->h[0], 1);
}
break;
}
} else { /* ..... Force specification for unsteady problems ..... */
if (dparam("FLOWATE") != 0.) /* check for flowrate */
option_set ("flowrate", 1);
else if (nlines) { /* read the list... */
for (i = 0; i < nforces && nlines--; i++) {
for(p = fgets(buf, BUFSIZ, fp); isspace(*p); p++);
if (toupper(*p) == 'C') {
/* do nothing */
} else {
if (p = strchr(buf,'='))
dparam_set(ftypes[i], scalar(++p));
else
fprintf(stderr,"ReadDF -- invalid function (below)");
}
}
}
if (option("flowrate"))
printf("Flowrate : %g\n", dparam("FLOWRATE"));
else {
printf("Drive Force : FFx = %g\n", dparam("FFX"));
printf(" FFy = %g\n", dparam("FFY"));
}
}
while (nlines>=0){
fgets(buf, BUFSIZ, fp); /* finish the section */
nlines--;
}
return;
}
/* ------------------------------------------------------------------------ *
* ReadBCs() - Read Boundary Conditions *
* *
* This function reads boundary information for a spectral element field *
* and sets up the appropriate boundary structures. The following boundary *
* conditions are currently recognized: *
* *
* Flag Type Field #1 Field #2 Field #3 *
* ----- -------- ---------- ---------- ---------- *
* V Velocity U-velocity V-Velocity W-Velocity *
* W Wall = 0 = 0 = 0 *
* v v U(x,y,z) V(x,y,z) W(x,y,z) *
* F Flux U' = f1 V' = f2 W' = f3 *
* f flux U'(x,y,z) V'(x,y,z) W'(x,y,z) *
* O Outflow U' = 0 V' = 0 W' = 0 *
* E Element w/ ID w/ EDGE *
* P Periodic w/ ID w/ EDGE *
* *
* NOTES: *
* *
* - On the first call, this function records the beginning of the BC *
* section. Subsequent calls will re-read the BC's and assign data *
* from the next available parameter. There should be a total of DIM *
* parameters available (DIM <= 3). *
* *
* - For boundaries with a given function of (x,y,z), the first line *
* following is taken to be U(x,y,z), the second is V(x,y,z), etc. *
* Here is an example (element 1, edge 1): *
* *
* 1 1 v 0.0 0.0 0.0 *
* ux = (1 - y^2)*sin(t) *
* uy = (1 - y^2)*cos(t) *
* uz = 0. *
* 1 2 E 2.0 1.0 0.0 *
* *
* The "=" MUST be present, and the function is defined by the string *
* to the right of it. *
* *
* - NO TIME DEPENDENT BOUNDARY CONDITIONS. *
* *
* - This routine will read any type of boundary conditions, FLUID or *
* otherwise, as long as the file is positioned correctly. *
* *
* ------------------------------------------------------------------------ */
#define on 1
#define off 0
static fpos_t bc_start, bc_stop;
static void prescan (FILE *fp),
posscan (FILE *fp);
Bndry *ReadBCs (FILE *fp, Element *U)
{
register int n,i;
char bc, buf[BUFSIZ];
int iel, iside, nbcs=0;
Element *E;
Bndry *bndry_list = (Bndry *) NULL,
*new_bndry = (Bndry *) NULL;
double f[3];
static int fnum = 0; /* this tracks the field numbers read */
if (!fnum) prescan(fp);
fsetpos(fp, &bc_start);
for(E=U;E;E=E->next){ /* ....... Begin Reading Boundary Conditions ...... */
for (n = 0; n < E->Nedges; ++n){
char *p = fgets (buf, BUFSIZ, fp);
/* vertex solve mask is on by default */
while (isspace(*p)) p++; bc = *p++;
sscanf(p, "%d%d%lf%lf%lf", &iel, &iside, f, f+1, f+2);
if (iel - E->id != 1 || iside - n != 1) {
sprintf(buf, "Mismatched element/side -- got %d,%d, expected %d,%d",
iel, iside, E->id+1, n+1);
fprintf(stderr,"ReadBCs : %s\n", buf); exit(-1);
}
switch (bc) {
case 'W':
f[fnum] = 0.;
case 'V':
E->set_solve(n,off); //OVERLOAD CALL: set_solve: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
new_bndry = E->gen_bndry(bc, n, f[fnum]); //OVERLOAD CALL: gen_bndry: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
new_bndry->usrtype = 'V';
break;
case 'Z':
f[fnum] = 0.;
if(E->type == 'u')
new_bndry = E->gen_bndry('F', n, f[fnum]); //OVERLOAD CALL: gen_bndry: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
else{
E->set_solve(n,off); //OVERLOAD CALL: set_solve: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
new_bndry = E->gen_bndry('W', n, f[fnum]); //OVERLOAD CALL: gen_bndry: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
}
new_bndry->usrtype = bc;
option_set("REFLECT", 1);
break;
case 'I': case 'O':
f[fnum] = 0.;
case 'F': case 'R':
new_bndry = E->gen_bndry(bc, n, f[fnum]); //OVERLOAD CALL: gen_bndry: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
new_bndry->usrtype = bc;
break;
/* Functional boundaries */
case 't': case 'v':
E->set_solve(n,off); //OVERLOAD CALL: set_solve: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
case 'f': case 'r': {
fpos_t pos;
for (i = 0; i <= fnum; i++) fgets(buf, BUFSIZ, fp);
new_bndry = E->gen_bndry(bc, n, buf); //OVERLOAD CALL: gen_bndry: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
new_bndry->usrtype = bc;
do {
fgetpos (fp, &pos); /* Search through the following */
fgets (buf, BUFSIZ, fp); /* lines for the first one without */
} while (strchr (buf, '=')); /* an '=' (function specification) */
fsetpos (fp, &pos);
break;
}
/* Element and Periodic boundaries */
case 'E': case 'P':
new_bndry = (Bndry *) NULL;
break;
default:
fputs("ReadBCs: read a strange boundary condition", stderr);
exit(-1);
break;
}
if (new_bndry) {
nbcs++;
new_bndry->next = bndry_list;
bndry_list = new_bndry;
}
}
if(E->identify() == Nek_Tri || E->identify() == Nek_Nodal_Tri) //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?); identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
fgets (buf, BUFSIZ, fp); /* skip extra side */
}
if (fnum++ == 0) posscan(fp); fsetpos(fp, &bc_stop);
return bsort(bndry_list,nbcs); /* Return the sorted list */
}
/* ------------------------------------------------------------------------ *
* bsort() - Sort Boundary Conditions *
* *
* The following function sorts boundary conditions so they are set accord- *
* ing to precedence and not mesh order. The following is the "hierarchy" *
* of boundary conditions : *
* *
* W : Walls (zero velocity) are the most binding *
* V : Velocity boundary conditions (constant) *
* v : Velocity function boundary conditions *
* all others *
* *
* This function returns the sordid boundary list (ha ha). *
* *
* ------------------------------------------------------------------------ */
Bndry *bsort(Bndry *bndry_list, int nbcs){
Bndry **tmp;
Bndry *bedg = bndry_list;
register int i;
if(!nbcs) return bndry_list;
tmp = (Bndry**) malloc(nbcs*sizeof(Bndry*));
/* Copy the linked list into a regular array and sort it */
for(i = 0 ; i < nbcs; ++i, bedg = bedg->next) tmp[i] = bedg;
qsort(tmp, nbcs, sizeof(Bndry*), bedgcmp);
/* Create a new version of the ordered list */
bedg = (Bndry*) malloc (nbcs * sizeof(Bndry));
for(i = 0; i < nbcs; ++i) {
// memcpy (bedg + i, tmp[i], sizeof(Bndry));
bedg[i] = *(tmp[i]);
bedg[i].id = i;
bedg[i].next = bedg + i + 1;
}
bedg[nbcs-1].next = (Bndry*) NULL;
/* Free up the space occupied by the old list */
bndry_list = bedg;
// for(i = 0; i < nbcs; ++i) free (tmp[i]);
// free (tmp);
return bndry_list;
}
/*
* Boundaries are sorted by their ASCII character values except for
* types 'D' and 'N', which are sorted by their element ID numbers.
*/
static int bedgcmp(const void *b1, const void *b2)
{
Bndry *bedg_1 = *((Bndry**) b1),
*bedg_2 = *((Bndry**) b2);
char btype_1 = bedg_1 -> type,
btype_2 = bedg_2 -> type;
/* Convert {D,N} -> {G,H} to get the ASCII precedence right */
if (btype_1 == 'D') btype_1 = 'G';
if (btype_2 == 'D') btype_2 = 'G';
if (btype_1 == 'N') btype_1 = 'H';
if (btype_2 == 'N') btype_2 = 'H';
if (btype_1 > btype_2) /* Check ASCII code */
return 1;
else if (btype_1 < btype_2)
return -1;
else /* Check element ID */
{
int id_1 = bedg_1 -> elmt -> id,
id_2 = bedg_2 -> elmt -> id;
if (id_1 > id_2)
return 1;
else if (id_1 < id_2)
return -1;
else /* Check edge ID */
{
#if DIM == 3
id_1 = bedg_1 -> face;
id_2 = bedg_2 -> face;
#else
id_1 = bedg_1 -> face;
id_2 = bedg_2 -> face;
#endif
if (id_1 > id_2)
return 1;
else if (id_1 < id_2)
return -1;
}
}
/* If we get this far there is an error in the boundary conditions */
fputs("bedgcmp(): " "Duplicate boundary\n",stderr); exit(-1);
return 0;
}
static void prescan(FILE *fp)
{
char buf[BUFSIZ];
rewind(fp);
findSection(SEC_BCS,buf,fp);
if (!strstr(buf, SEC_BCS))
error_msg(ReadBCs: do not see any boundary conditions); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
do
fgets(buf, BUFSIZ, fp);
while
(strstr(buf, "NO"));
fgetpos (fp, &bc_start);
return;
}
static void posscan (FILE *fp)
{
char buf[BUFSIZ];
fgetpos (fp, &bc_stop);
while (strstr(fgets(buf,BUFSIZ,fp), SEC_BCS))
fgetpos (fp, &bc_stop);
return;
}
/* set up inter element links using the coundary connectivity information */
static void set_link (Element_List *U, int eid1, int id1, int eid2, int id2); //OVERLOAD CALL: set_link: Refine.c(?), io.c(?)
void ReadSetLink(FILE *fp, Element_List *U)
{
register int n;
char bc, buf[BUFSIZ];
int iel, iside;
Element *E;
double f[3];
rewind(fp);
findSection(SEC_BCS,buf,fp);
if (!strstr(buf, SEC_BCS))
error_msg(ReadSetLink: do not see any boundary conditions); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
fgets(buf,BUFSIZ,fp);
/* ....... Search through Boundary Conditions ...... */
char *p = fgets(buf,BUFSIZ, fp);
for(E=U->fhead;E;E=E->next){
for (n = 0; n < E->Nedges; ++n){
while (isspace(*p)) p++; bc = *p++;
sscanf(p,"%d%d%lf%lf%lf", &iel, &iside, f, f+1, f+2);
if (iel - E->id != 1 || iside - n != 1) {
sprintf(buf, "Mismatched element/side -- got %d,%d, expected %d,%d",
iel, iside, E->id+1, n+1);
fprintf(stderr,"ReadSetLink : %s\n", buf); exit(-1);
}
if ((bc == 'E')||(bc == 'P')) /* Element and Periodic boundaries */
set_link(U,iel-1,iside-1,(int)f[0]-1,(int)f[1]-1); //OVERLOAD CALL: set_link: Refine.c(?), io.c(?)
if ((bc == 'v')||(bc == 'f')||(bc == 'r'))
while(strchr(p=fgets(buf,BUFSIZ,fp),'='));
else
p=fgets(buf,BUFSIZ,fp);
}
#if DIM == 2
if(E->Nverts==3){
p=fgets(buf,BUFSIZ,fp);
}
#endif
}
return;
}
#ifdef NODAL
int con[4][4] = {{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
#else
int con[4][4] = {{1,1,0,0},{1,1,0,0},{0,0,1,1},{0,0,1,1}};
#endif
static void set_link(Element_List *U, int eid1, int id1, int eid2, int id2){
/* set up connectivity inverses if required */
/* like faces meet */
Element *E=U->flist[eid1],*F=U->flist[eid2];
if(eid1>eid2)
E->edge[id1].con = con[id1][id2];
else
F->edge[id2].con = con[id1][id2];
/* set links */
if(eid1 < eid2){
E->edge[id1].base = E->edge + id1;
E->edge[id1].link = F->edge + id2;
F->edge[id2].base = E->edge + id1;
}
}
void Tri_work();
void Quad_work();
void Tet_work();
void Prism_work();
void Hex_work();
void Pyr_work();
void Quad_red_work();
void ReadOrderFile(char *name, Element_List *E){
register int i,k;
int *size,n;
const int nel = E->nel;
char buf[BUFSIZ],*p;
FILE *fp;
sprintf(buf,"%s.ord",name);
if(!(fp = fopen(buf,"r"))){
fprintf(stderr,"ERROR: file \"%s\" does not exist\n",buf);
exit(-1);
}
size = ivector(0,5);
/* skip comments */
while(p=fgets(buf,BUFSIZ,fp))
if(strstr(p,"Modes")){
for(k=0;k<nel;++k){
fscanf(fp,"%*d");
for(i = 0; i < E->flist[i]->Nedges; ++i)
fscanf(fp,"%d",size+i);
fscanf(fp,"%d",size+i);
E->flist[k]->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->flist[k]->Mem_Q(); //OVERLOAD CALL: Mem_Q: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
}
break;
}
else if (strstr(p,"Function")){
#if DIM == 2
double x[5],y[5],f[5];
#else
double x[11],y[11],z[11],f[11];
extern int ednum[][3];
#endif
p = fgets(buf,BUFSIZ,fp);
if((p = (strchr(p,'='))) == (char *) NULL)
error_msg(CheckOrderFIle: Illegal function definition); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
#if DIM == 2
vector_def("x y",++p);
#else
vector_def("x y z",++p);
#endif
int max_f;
Element *F;
for(k = 0; k < nel; ++k){
F=E->flist[k];
#if DIM == 2
for(i = 0; i < F->Nverts; ++i){
x[i] = 0.5*(F->vert[i].x + F->vert[(i+1)%F->Nverts].x);
y[i] = 0.5*(F->vert[i].y + F->vert[(i+1)%F->Nverts].y);
}
vector_set(F->Nverts,x,y,f);
/* set fourth side to be consistent with highest edge */
max_f=0;
for(i = 0; i < F->Nedges; ++i)
max_f = (int) ((f[i]>max_f) ? f[i]:max_f);
f[F->Nverts] = (F->identify()==Nek_Tri) ? max_f-1:max_f; //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
if(n=option("aposterr"))
for(i = 0; i < F->Nedges+1; ++i) size[i] = (int)max(f[i],0.0) + n;
else
for(i = 0; i < F->Nedges+1; ++i) size[i] = (int)max(f[i],0.0);
#else
for(i = 0; i < 3; ++i){
x[i] = 0.5*(E[k].vert[i].x + E[k].vert[(i+1)%3].x);
y[i] = 0.5*(E[k].vert[i].y + E[k].vert[(i+1)%3].y);
z[i] = 0.5*(E[k].vert[i].z + E[k].vert[(i+1)%3].z);
}
for(i = 0; i < 3; ++i){
x[i+3] = 0.5*(E[k].vert[i].x + E[k].vert[3].x);
y[i+3] = 0.5*(E[k].vert[i].y + E[k].vert[3].y);
z[i+3] = 0.5*(E[k].vert[i].z + E[k].vert[3].z);
}
vector_set(6,x,y,z,f);
/* set faces to be maximun of surrounding edges */
/* and edges to be maximun of faces */
f[10] = 0;
for(i = 0; i < 4; ++i){
f[i+6] = max(max(f[ednum[i][0]],f[ednum[i][1]]),f[ednum[i][2]])-1;
f[10] = max(f[10],f[i+6]);
}
--f[10];
for(i = 0; i < 11; ++i) size[i] = (int)max(f[i],0.0);
#endif
F->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)
F->Mem_Q(); //OVERLOAD CALL: Mem_Q: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
}
break;
}
option_set("variable",1);
Tri_work();
Quad_work();
Tet_work();
Prism_work();
Hex_work();
Pyr_work();
Quad_red_work();
free(size);
}
/* ----------------------------------------------------------------------- *
* ReadHisData() - Read History and Integral data *
* *
* This function reads in and sets up structures for the processing of *
* history points. This is comprised of a linked list of structures which *
* are processed by the Analyzer(). *
* *
* ----------------------------------------------------------------------- */
static HisPoint *appendHisPoint (HisPoint *list, HisPoint *hp); //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?); HisPoint: nektar.h(?), nekcomp.h(?); HisPoint: nektar.h(?), nekcomp.h(?)
void ReadHisData (FILE *fp, Domain *omega) //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
{
char buf[BUFSIZ], *p;
int cnt, npts, nel;
HisPoint *hp, *his_list = (HisPoint*)NULL; //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?); HisPoint: nektar.h(?), nekcomp.h(?)
if (!findSection(SEC_HIST, buf, fp))
{fputs("ReadHisData: Section not found", stderr); exit(-1);}
if (sscanf(fgets(buf, BUFSIZ, fp), "%d", &npts) != 1)
{fputs("ReadHisData: can't read the number of points", stderr); exit(-1);}
nel = iparam("ELEMENTS");
while (npts--) {
hp = (HisPoint*) calloc(1, sizeof(HisPoint)); //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?); HisPoint: nektar.h(?), nekcomp.h(?)
for (p = fgets(buf, BUFSIZ, fp); !isdigit(*p); p++)
if (isalpha(*p=tolower(*p)))
switch(*p){
case 'h': /* treat as vertex point */ //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
hp->mode = TransVert;
break;
case 'e': /* treat as mid point of edge */ //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
hp->mode = TransEdge;
break;
default:
strncat(hp->flags, p, 1);
}
hp->i = (int) strtol (p, &p, 10);
hp->j = (int) strtol (p, &p, 10);
hp->k = (int) strtol (p, &p, 10);
hp->id = atoi (p);
if ((--hp->i) > Max_Nedges || (hp->id--) > nel)
{ free(hp); continue; }
his_list = appendHisPoint (his_list, hp);
}
/* This history step can be overridden by setting HISSTEP */
option_set ("hisstep", iparam ("HISSTEP") > 0 ?
iparam ("HISSTEP") : max(1, iparam("NSTEPS") / 1000));
if (omega->his_list = his_list) {
/* Echo the history points */
puts("\nHistory points:");
for (hp = his_list, cnt = 0; hp; hp = hp->next) {
if(hp->mode == TransVert)
printf ("%2d: vertice = %2d, elmt = %3d, fields = %s\n", ++cnt,
hp->i + 1, hp->id + 1, hp->flags);
else if(hp->mode == TransEdge)
printf ("%2d: edge = %2d, elmt = %3d, fields = %s\n", ++cnt,
hp->i + 1, hp->id + 1, hp->flags);
}
}
putchar ('\n');
return;
}
static HisPoint *appendHisPoint (HisPoint *list, HisPoint *hp) //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?); HisPoint: nektar.h(?), nekcomp.h(?); HisPoint: nektar.h(?), nekcomp.h(?)
{
HisPoint *h = list; //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?)
if (h) {
while (h->next) h = h->next;
h->next = hp;
} else
list = hp;
return list;
}
void ReadSoln(FILE* fp, Domain* omega){ //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
char buf[BUFSIZ], *s;
int i;
#ifdef THERMO
int nfields = DIM+2;
#else
int nfields = DIM+1;
#endif
rewind(fp);
omega->soln = (char**) malloc(nfields*sizeof(char*));
if(findSection("Exact", buf, fp)){
for(i=0;i<nfields;++i){
if(s = strchr(fgets(buf,BUFSIZ,fp),'=')){
while (isspace(*++s));
omega->soln[i] = (char*) malloc(BUFSIZ*sizeof(char));
strcpy(omega->soln[i],s);
}
else
sprintf(omega->soln[i],"0.0");
}
fprintf(stdout,"Exact U : %s",omega->soln[0]);
fprintf(stdout,"Exact V : %s",omega->soln[1]);
#ifdef THERMO
fprintf(stdout,"Exact T : %s",omega->soln[2]);
fprintf(stdout,"Exact P : %s",omega->soln[3]);
#else
fprintf(stdout,"Exact P : %s",omega->soln[2]);
#endif
}
else
for(i=0;i<nfields;++i){
omega->soln[i] = (char*) malloc(BUFSIZ*sizeof(char));
sprintf(omega->soln[i],"0.0");
}
rewind(fp);
}
void ReadKinvis(Domain* omega){ //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
int i,k, skip;
Element_List *U = omega->U;
Element *E;
double sponpos = dparam("SPONPOS");
double sponlen = dparam("SPONLEN");
double sponmag = dparam("SPONMAG");
omega->kinvis = (Metric*) calloc(U->nel, sizeof(Metric));
if(sponlen)
omega->kinvis->p = dvector(0, U->htot-1);
Coord X; //OVERLOAD CALL: Coord: nekstruct.h(?), hotel.h(?)
X.x = dvector(0, QGmax*QGmax-1);
X.y = dvector(0, QGmax*QGmax-1);
for(skip=0, k=0;k<U->nel;++k){
E = U->flist[k];
E->coord(&X); //OVERLOAD CALL: coord: Coords.c(Tet), Coords.c(Pyr), Coords.c(Prism), Coords.c(Hex), Coords.c(Element), Coords.c(Tri), Coords.c(Quad)
omega->kinvis[k].d = dparam("KINVIS");
if(sponlen){
omega->kinvis[k].p = omega->kinvis->p + skip;
for(i=0;i<E->qtot;++i)
omega->kinvis[k].p[i] = (X.x[i] < sponpos) ? 1.:
// (1+cos(M_PI*(X.x[i]-sponpos)/sponlen))*.5;
(1. + sponmag*(X.x[i]-sponpos)*(X.x[i]-sponpos));
}
skip += E->qtot;
}
if(sponlen)
fprintf(stdout, "ReadKinvis: sponge length: %lf sponge position: %lf \n",
sponlen, sponpos);
free(X.x); free(X.y);
}
void ReadWave(FILE *fp, double **wave, Element_List *U){
int i;
int eDIM = U->fhead->dim(); //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
char buf[BUFSIZ], *p;
if (!findSection (SEC_WAVE, buf, fp)){
fputs("ReadWave: section not found\n", stderr);
rewind(fp);
return;
}
if(p = strchr (fgets(buf, BUFSIZ, fp), '=')) {
while (isspace(*++p));
U->Set_field(p); //OVERLOAD CALL: Set_field: Element_List.c(Element_List), Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
for(i = 0; i < U->nel; ++i)
dcopy(U->flist[i]->qtot,U->flist[i]->h[0],1,wave[i*eDIM],1);
}
if(p = strchr (fgets(buf, BUFSIZ, fp), '=')) {
while (isspace(*++p));
U->Set_field(p); //OVERLOAD CALL: Set_field: Element_List.c(Element_List), Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
for(i = 0; i < U->nel; ++i)
dcopy(U->flist[i]->qtot,U->flist[i]->h[0],1,wave[i*eDIM+1],1);
}
}
C++ to HTML Conversion by ctoohtml