file: Nektar2d/src/forces.c


/*---------------------------------------------------------------------------*
 *                        RCS Information                                    *
 *                                                                           *
 * $RCSfile:
 * $Revision:
 * $Author: 
 * $Date: 
 * $State:
 * ------------------------------------------------------------------------- */
#ifdef FORCES
#include "nektar.h"

static int init=1; /* externals */

static void set_SurGeofac(Bndry *Pbc, Bndry *Ubc);

/*-------------------------------------------------------------------*
 * This is a function to calculate the forces of the fluid acting on *
 * a body pressumed to be represented by the flag 'W'.               *
 *                                                                   *
 *   F_i = P n_i - RHO*KINVIS*T_ij n_j  (RHO=1)                      *
 *                                                                   *
 *-------------------------------------------------------------------*/
#if DIM == 2
void forces(Domain *omega, double time){     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
  register int i;
  int      eid,qa,face,lmax=0;
  Bndry   *B = omega->Ubc;
  double   fp[2],fv[2],*w;
  double   kinvis = dparam("KINVIS");
  double  *wk = dvector(0,max(LGmax+1,QGmax)*QGmax-1);
  double  **D = dmatrix(0,4,0,QGmax*QGmax-1);
  double  *ux = D[0], *uy = D[1];
  double  *vx = D[2], *vy = D[3], *p = D[4];
  Element_List *U  = omega->U, *V = omega->V, *P = omega->P;
  Element *eU, *eV, *eP;
  FILE    *fout = omega->fce_file;
  Basis   *b;
  
  fp[0] = fp[1] = 0.0;
  fv[0] = fv[1] = 0.0;

  /* write out header */
  if(init){
    fprintf(fout,"# Force acting on body\n");
    fprintf(fout,"# \n");
    fprintf(fout,"# Time \t (Fx-press, Fx-visc)Fx \t (Fy-press, Fy-visc) Fy\n");

    /* need to set up surface geometric factors in U from P */
    set_SurGeofac(omega->Pbc,omega->Ubc);

    init = 0;
  }

  for(;B; B = B->next)
#ifdef OSSALE
    if(B->type == 'V'){
#else
    if(B->type == 'W'){
#endif
      face = B->face;
      eid  = B->elmt->id;
      eU   = U->flist[eid];
      eV   = V->flist[eid];
      eP   = P->flist[eid];

      qa   = eU->qa;
      
      if(eU->lmax != lmax){
	b    = eU->getbasis();        //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
	lmax = eU->lmax;
	getzw(qa,&w,&w,'a');
      }

      if((SLVTYPE) iparam("SLVTYPE") == StokesSlv){
#ifndef PCONTBASE
	eP->Obwd(eP->vert->hj,*eP->h,eP->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)
	eP->state = 't';	
#else
	eP->Jbwd(eP,b);         //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
	eP->state = 't';
#endif
      }
      else{
	eP->Jbwd(eP,b);         //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
	eP->state = 't';
      }

      eU->Jbwd(eU,b); eU->state = 't';        //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      eV->Jbwd(eV,b); eV->state = 't';        //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      
      eU->Grad_d( ux, uy, 0, 'a' );        //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)
      eV->Grad_d( vx, vy, 0, 'a' );        //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)

      /* get appropriate face values from grad(V) */
      for(i = 0; i < 4; ++i){
	eU->GetFace(D[i] ,face,wk);        //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
	eU->InterpToFace1(face,wk,D[i]);        //OVERLOAD CALL: InterpToFace1: Interp.c(Tri), Interp.c(Quad), Interp.c(Tet), Interp.c(Pyr), Interp.c(Prism), Interp.c(Hex), Interp.c(Element)
      }

      eP->GetFace(eP->h[0],face,wk);        //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
      eP->InterpToFace1(face,wk,p);        //OVERLOAD CALL: InterpToFace1: Interp.c(Tri), Interp.c(Quad), Interp.c(Tet), Interp.c(Pyr), Interp.c(Prism), Interp.c(Hex), Interp.c(Element)
      
      /* calculate forces */
      if(eU->curvX)
	for(i = 0; i < qa; ++i){
	  fp[0] += p[i]*B->nx.p[i]*B->sjac.p[i]*w[i];
	  fp[1] += p[i]*B->ny.p[i]*B->sjac.p[i]*w[i];
	  fv[0] -= kinvis*(2.0*ux[i]*B->nx.p[i] + (uy[i]+vx[i])*B->ny.p[i]
	             )*B->sjac.p[i]*w[i];
	  fv[1] -= kinvis*(2.0*vy[i]*B->ny.p[i] + (uy[i]+vx[i])*B->nx.p[i]
		     )*B->sjac.p[i]*w[i];
	}
      else
	for(i = 0; i < qa; ++i){
	  fp[0] += p[i]*B->nx.d*B->sjac.d*w[i];
	  fp[1] += p[i]*B->ny.d*B->sjac.d*w[i];
	  fv[0] -= kinvis*(2.0*ux[i]*B->nx.d + (uy[i]+vx[i])*B->ny.d
	             )*B->sjac.d*w[i];
	  fv[1] -= kinvis*(2.0*vy[i]*B->ny.d + (uy[i]+vx[i])*B->nx.d
		     )*B->sjac.d*w[i];
	}
    }

  fprintf(fout,"%lf ( %lf %lf ) %lf ( %lf %lf ) %lf\n",time,fp[0],fv[0],
	  fp[0]+fv[0],fp[1],fv[1],fp[1]+fv[1]);

  free(wk); free_dmatrix(D,0,0);
}
#
#else /* 3D version */
#
void forces(Domain *omega, double time){     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
  register int i,j;
  int      eid, qa, qb, face, lmax=0;
  Bndry   *B = omega->Ubc;
  double  fp[3],fv[3],*wa,*wb;
  double  kinvis = dparam("KINVIS");
  double  *w  = dvector(0,QGmax*QGmax-1);
  double  *wk = dvector(0,3*QGmax*QGmax*QGmax-1);
  double  **D = dmatrix(0,9,0,QGmax*QGmax*QGmax-1);
  double  *ux = D[0], *uy = D[1], *uz = D[2];
  double  *vx = D[3], *vy = D[4], *vz = D[5];
  double  *wx = D[6], *wy = D[7], *wz = D[8], *p = D[9];
  Element *U  = omega->Us, *V = omega->Vs, *W = omega->Ws, *P = omega->P;
  FILE    *fout = omega->fce_file;
  Basis   *b;
  
  fp[0] = fp[1] = fp[2] = 0.0;
  fv[0] = fv[1] = fv[2] = 0.0;
  
  /* print header */
  if(init){
    fprintf(fout,"# Force acting on body\n");
    fprintf(fout,"# \n");
    fprintf(fout,"# Time  (Fx-press, Fx-visc) Fx  (Fx-press, Fx-visc)"
	    "Fy  (Fx-press, Fx-visc) Fz \n");
    
    /* need to set up surface geometric factors in U from P */
    set_SurGeofac(omega->Pbc,omega->Ubc);
    
    init = 0;
  }

  for(;B; B = B->next)
    if(B->type == 'W'){
      face = B->face;
      eid  = B->elmt->id;
      
      if(U[eid].lmax != lmax){
	qa   = U[eid].qa;
	qb   = U[eid].qb;
	b    = getbasis(U+eid);        //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
	lmax = U[eid].lmax;
	getzw(qa,&wa,&wa,'a');
	getzw(qb,&wb,&wb,'b');

	for(i = 0; i < qb; ++i)
	  dcopy(qa,wa,1,w+i*qa,1);
	for(i = 0; i < qb; ++i)
	  dsmul(qa,wb[i],w+i*qa,1,w+i*qa,1);
      }
      
      EJbwd(P+eid,P+eid,b,wk); P[eid].state = 't';        //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      EJbwd(U+eid,U+eid,b,wk); U[eid].state = 't';        //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      EJbwd(V+eid,V+eid,b,wk); V[eid].state = 't';        //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      EJbwd(W+eid,W+eid,b,wk); W[eid].state = 't';        //OVERLOAD CALL: Jbwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      
      Egradv (U + eid, **U[eid].h, ux, uy, uz, 'a', wk);
      Egradv (V + eid, **V[eid].h, vx, vy, vz, 'a', wk);
      Egradv (W + eid, **W[eid].h, wx, wy, wz, 'a', wk);
 
      /* get appropriate face values from D */
      for(i = 0; i < 9; ++i){
	GetFace(U+eid,D[i],face,wk);        //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
	InterpToFace1(U+eid,face,wk,D[i]);        //OVERLOAD CALL: InterpToFace1: Interp.c(Tri), Interp.c(Quad), Interp.c(Tet), Interp.c(Pyr), Interp.c(Prism), Interp.c(Hex), Interp.c(Element)
      }

      GetFace(P+eid,**P[eid].h,face,wk);        //OVERLOAD CALL: GetFace: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
      InterpToFace1(P+eid,face,wk,p);        //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(U[eid].curvX)
	for(i = 0; i < qa*qb; ++i){
	  fp[0] += p[i]*B->nx.p[i]*B->sjac.p[i]*w[i];
	  fp[1] += p[i]*B->ny.p[i]*B->sjac.p[i]*w[i];
	  fp[2] += p[i]*B->nz.p[i]*B->sjac.p[i]*w[i];

	  fv[0] -= kinvis*(2.0*ux[i]*B->nx.p[i] + (uy[i] + vx[i])*B->ny.p[i]
		     + (uz[i] + wx[i])*B->nz.p[i])*B->sjac.p[i]*w[i];
	  fv[1] -= kinvis*(2.0*vy[i]*B->ny.p[i] + (uy[i] + vx[i])*B->nx.p[i]
		     + (vz[i] + wy[i])*B->nz.p[i])*B->sjac.p[i]*w[i];
	  fv[2] -= kinvis*(2.0*wz[i]*B->nz.p[i] + (uz[i] + wx[i])*B->nx.p[i]
		     + (vz[i] + wy[i])*B->ny.p[i])*B->sjac.p[i]*w[i];
	}
      else
	for(i = 0; i < qa*qb; ++i){
	  fp[0] += p[i]*B->nx.d*B->sjac.d*w[i];
	  fp[1] += p[i]*B->ny.d*B->sjac.d*w[i];
	  fp[2] += p[i]*B->nz.d*B->sjac.d*w[i];

	  fv[0] -= kinvis*(2.0*ux[i]*B->nx.d + (uy[i] + vx[i])*B->ny.d
		     + (uz[i] + wx[i])*B->nz.d)*B->sjac.d*w[i];
	  fv[1] -= kinvis*(2.0*vy[i]*B->ny.d + (uy[i] + vx[i])*B->nx.d
		     + (vz[i] + wy[i])*B->nz.d)*B->sjac.d*w[i];
	  fv[2] -= kinvis*(2.0*wz[i]*B->nz.d + (uz[i] + wx[i])*B->nx.d
		     + (vz[i] + wy[i])*B->ny.d)*B->sjac.d*w[i];
	}
    }

  fprintf(fout,"%lf ( %lf %lf ) %lf ( %lf %lf ) %lf ( %lf %lf ) %lf \n",
	  time,fp[0],fv[0],fp[0]+fv[0],fp[1],fv[1],fp[1]+fv[1],fp[2],
	  fv[2],fp[2]+fv[2]);

  free(w); free(wk); free_dmatrix(D,0,0);
}  
#endif

static void set_SurGeofac(Bndry *Pbc, Bndry *Ubc){
  Bndry *Ebc;
  
  if(iparam("SLVTYPE") == StokesSlv){ /* need to setup Ubc directly */
    for(Ebc = Ubc; Ebc; Ebc = Ebc->next)
      if(Ebc->type == 'W')
	Ebc->elmt->Surface_geofac(Ebc);        //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)
  }
  else
    for(;Pbc; Pbc = Pbc->next)
      if(Pbc->type == 'F')
	for(Ebc = Ubc; Ebc; Ebc = Ebc->next)
	  if((Ebc->elmt->id == Pbc->elmt->id)&&(Ebc->face == Pbc->face)){
	    Ebc->sjac = Pbc->sjac;
	    Ebc->nx   = Pbc->nx;
	    Ebc->ny   = Pbc->ny;
#if DIM == 3
	    Ebc->nz   = Pbc->nz;
#endif
	  } 
}
#endif



Back to Source File Index


C++ to HTML Conversion by ctoohtml