file: Nektar2d/src/pressure.c
/*---------------------------------------------------------------------------*
* RCS Information *
* *
* $Source: /disk2/tcew/Hybrid/Nektar2d/src/RCS/pressure.C,v $
* $Revision: 1.1 $
* $Date: 1996/06/10 20:43:25 $
* $Author: tcew $
* $State: Exp $
*---------------------------------------------------------------------------*/
#include "nektar.h"
#include "string.h"
/*
* Build boundary conditions for the pressure
*/
Bndry *BuildPBCs(Element_List *P, Bndry *temp)
{
int nbcs = 0, Je = iparam("INTYPE");
Bndry *Pbc;
Element *E;
register int i,j;
/* Count the number of boundaries to create */
if(!temp)
return (Bndry*)0;
while (temp[nbcs++].next);
Pbc = (Bndry*)calloc(nbcs, sizeof(Bndry));
for(i = 0; i < nbcs; ++i) {
Pbc[i].id = i;
Pbc[i].type = temp[i].type;
Pbc[i].usrtype = temp[i].usrtype;
if(Pbc[i].usrtype == 'Z') // override for symmetry
Pbc[i].type = 'V';
Pbc[i].elmt = P->flist[temp[i].elmt->id];
Pbc[i].face = temp[i].face;
Pbc[i].next = Pbc + i + 1;
}
Pbc[nbcs-1].next = (Bndry*) NULL;
/* clear vertex solve mask */
for(E = P->fhead; E; E = E->next)
for(i = 0; i < E->Nverts; ++i)
E->vert[i].solve = 1;
/* Translate the boundary conditions */
for(i = 0; i < nbcs; ++i) {
switch (Pbc[i].type) {
case 'O':
Pbc[i].type = 'o'; /* fall through */
E = Pbc[i].elmt;
for(j = 0; j < DIM; ++j)
E->vert[E->vnum(Pbc[i].face,j)].solve = 0; //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
E->MemBndry(Pbc+i,Pbc[i].face,(iparam("EQTYPE") == Rotational)? Je:1); //OVERLOAD CALL: MemBndry: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
break;
case 'V': case 'v': /* allocate additional multi-step storage */
case 'W':
Pbc[i].type = 'F';
Pbc[i].elmt->Surface_geofac(Pbc+i); //OVERLOAD CALL: Surface_geofac: Coords.c(Tri), Coords.c(Quad), Coords.c(Tet), Coords.c(Pyr), Coords.c(Prism), Coords.c(Hex), Coords.c(Element)
Pbc[i].elmt->MemBndry(Pbc+i,Pbc[i].face,Je); //OVERLOAD CALL: MemBndry: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
break;
default:
Pbc[i].elmt->MemBndry(Pbc+i,Pbc[i].face,1); //OVERLOAD CALL: MemBndry: Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
break;
}
}
return bsort(Pbc,nbcs);
}
/* ------------------------------------------------------------------------- *
* SetPBCs() - Set boundary conditions for the pressure *
* *
* Je 1 *
* dP/dn = SUM beta[q] * [ N(u) - - curl ( curl u ) ] * n *
* q=0 R *
* *
* where n is the unit outward normal along the edge, u is the velocity *
* field, and N(u) are the non-linear terms in the momentum equation. This *
* routine computes the RHS (already integrated) of the above equation at *
* the current time level and saves it in the corresponding Bedge for Press. *
* ------------------------------------------------------------------------- */
static void CalcPbc (Bndry *B, Element_List *V[3], Element_List *N[3], Metric *lambda);
static void IntPbc (Bndry *B, int Je);
static void outflow (Bndry *B, Domain *omega); //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
#ifdef OSSCYL
static void Addforce(Bndry *B);
#endif
void SetPBCs(Domain *omega){ //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
Bndry *Pbc = omega->Pbc;
Element_List *V[3], *N[3];
int Je = iparam("INTYPE");
V[0] = omega->U, V[1] = omega->V;
N[0] = omega->Uf, N[1] = omega->Vf;
/* Get the integration coefficients */
while (Pbc) {
switch (Pbc->type) {
case 'D': case 'N': case 'P':
break;
case 'o':
if (iparam("EQTYPE") == Rotational){
outflow (Pbc, omega);
IntPbc (Pbc, Je);
}
break;
case 'F': case 'R': {
// fprintf(stderr, " %d " , Pbc->elmt->id);
CalcPbc(Pbc,V,N,omega->kinvis+Pbc->elmt->id);
IntPbc(Pbc, Je);
#ifdef OSSCYL
Addforce(Pbc);
#endif
break;
}
default:
error_msg(SetPBCs -- unknown pressure b.c.) //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
break;
}
Pbc = Pbc->next;
}
return;
}
static void IntPbc(Bndry *B, int Je){
register int i,j;
const int face = B->face;
int ll;
double beta[3];
double *tmp = dvector(0,LGmax-1);
getbeta (beta);
/* integrate vertices */
dzero(DIM,tmp,1);
for(i = 0; i < DIM; ++i)
for(j = 0; j < Je; ++j)
tmp[i] += beta[j]*B->bvert[i+j*DIM];
dcopy(DIM * (Je-1), B->bvert, 1, B->bvert+DIM,1);
dcopy(DIM , tmp, 1, B->bvert,1);
ll = B->elmt->edge[face].l;
dzero(ll,tmp,1);
for(j = 0; j < Je; ++j)
daxpy(ll, beta[j], B->bedge[0] + j*ll, 1, tmp, 1);
dcopy(ll * (Je-1), B->bedge[0], 1, B->bedge[0] + ll, 1);
dcopy(ll , tmp , 1, B->bedge[0] , 1);
free(tmp);
}
#ifdef OSSCYL
/* add acceleration foring term onto pressure boundary conditions */
static void Addforce(Bndry *B){
if(!(B->elmt->curvX)&&(B->ny.d == 0.0)){
/* this is a fix to get end bndry's */
const int id = B->elmt->id, qa = B->elmt->qa;
double time = dparam("TIME"), dt = dparam("DT"),KC = dparam("KC");
double *Q = dvector(0,qa-1),
*store = dvector(0,LGmax-1),f;
f = (sin(2*M_PI/KC*(time+dt)) - sin(2*M_PI/KC*time))/dt;
dfill(qa,f*B->nx.d,Q,1);
/* store old values */
store[0] = B->bvert[0];
store[1] = B->bvert[1];
dcopy(B->elmt->edge[B->face].l,B->bedge[0],1,store+2,1);
B->elmt->MakeFlux(B,0,Q); //OVERLOAD CALL: MakeFlux: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
/* add stored values */
B->bvert[0] += store[0];
B->bvert[1] += store[1];
dvadd(B->elmt->edge[B->face].l,store+2,1,B->bedge[0],1,B->bedge[0],1);
free(Q); free(store);
}
}
#endif
/* --------------------------------------------------------------------- *
* This is a function to caclulate the flux integral of *
* dP/dn = Int {g, n * [N(u) - 1/Re Curl (Curl U)]} *
* Given N and U. The value is then stored in B. *
* --------------------------------------------------------------------- */
static void CalcPbc(Bndry *B, Element_List *V[3],
Element_List *N[3], Metric *kinvis){
const int id = B->elmt->id, qa = V[0]->flist[id]->qa,
tot = qa*V[0]->flist[id]->qb;
double *Q = dvector(0,tot-1),
*Uy = dvector(0,tot-1),
*Vx = dvector(0,tot-1),
*GxWx = dvector(0,tot-1),
*GxWy = dvector(0,tot-1);
V[0]->flist[id]->Grad_d ( 0, Uy, 0,'y'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
V[1]->flist[id]->Grad_d ( Vx, 0, 0,'x'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvsub (tot, Vx, 1, Uy, 1, Q, 1);
V[0]->flist[id]->Grad_h(Q, GxWy, GxWx,(double*)0, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dsvtvp (tot, -kinvis->d, GxWx, 1, *N[0]->flist[id]->h, 1, GxWx, 1);
dsvtvp (tot, kinvis->d, GxWy, 1, *N[1]->flist[id]->h, 1, GxWy, 1);
V[0]->flist[id]->GetFace(GxWx, B->face, GxWx); //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
V[0]->flist[id]->GetFace(GxWy, B->face, GxWy); //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
V[0]->flist[id]->InterpToFace1(B->face,GxWx,Q); //OVERLOAD CALL: InterpToFace1: Interp.c(Tri), Interp.c(Quad), Interp.c(Tet), Interp.c(Pyr), Interp.c(Prism), Interp.c(Hex), Interp.c(Element)
V[0]->flist[id]->InterpToFace1(B->face,GxWy,GxWx); //OVERLOAD CALL: InterpToFace1: Interp.c(Tri), Interp.c(Quad), Interp.c(Tet), Interp.c(Pyr), Interp.c(Prism), Interp.c(Hex), Interp.c(Element)
if(V[0]->flist[id]->curvX){
dvmul (qa,B->nx.p,1,Q ,1,Q,1);
dvvtvp (qa,B->ny.p,1,GxWx,1,Q,1,Q,1);
}
else {
dscal (qa,B->nx.d,Q ,1);
dsvtvp (qa,B->ny.d,GxWx,1,Q,1,Q,1);
}
B->elmt->MakeFlux(B,0,Q); //OVERLOAD CALL: MakeFlux: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
free(Uy); free(Vx); free(Q);
free(GxWx); free(GxWy);
}
/*
* PI = p + 1/2 U.U (outflow boundary conditions)
*/
static void outflow (Bndry *B, Domain *omega){ //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
const int id = B->elmt->id, face = B->face;
Element *U,*V;
double *u = dvector(0,QGmax-1), *v = dvector(0,QGmax-1);
int q;
U = omega->U->flist[id];
V = omega->V->flist[id];
if(U->identify() == Nek_Tri) //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
if(face == 0){
q = U->qa;
B->bvert[1] = 0.5*(U->vert[1].hj[0]*U->vert[1].hj[0] +
V->vert[1].hj[0]*V->vert[1].hj[0]);
}
else{
q = U->qb;
B->bvert[1] = 0.5*(U->vert[2].hj[0]*U->vert[2].hj[0] +
V->vert[2].hj[0]*V->vert[2].hj[0]);
}
else
switch(face){
case 0:
q = U->qa;
B->bvert[1] = 0.5*(U->vert[1].hj[0]*U->vert[1].hj[0] +
V->vert[1].hj[0]*V->vert[1].hj[0]);
break;
case 1:
q = U->qb;
B->bvert[1] = 0.5*(U->vert[2].hj[0]*U->vert[2].hj[0] +
V->vert[2].hj[0]*V->vert[2].hj[0]);
break;
case 2:
q = U->qa;
B->bvert[1] = 0.5*(U->vert[2].hj[0]*U->vert[2].hj[0] +
V->vert[2].hj[0]*V->vert[2].hj[0]);
case 3:
q = U->qb;
B->bvert[1] = 0.5*(U->vert[3].hj[0]*U->vert[3].hj[0] +
V->vert[3].hj[0]*V->vert[3].hj[0]);
break;
}
U->GetFace(*U->h, face, u); //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
V->GetFace(*V->h, face, v); //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (q, u, 1, u, 1, u, 1);
dvvtvp (q, v, 1, v, 1, u, 1, u, 1);
dscal (q, 0.5, u, 1);
B->bvert[0] = u[0];
B->elmt->JtransEdge(B,face,0,u); //OVERLOAD CALL: JtransEdge: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
free(u); free(v);
return;
}
C++ to HTML Conversion by ctoohtml