{
int ndim = pV.num_elements();
Array<OneD, Array<OneD, NekDouble> > AdvVel (pV.num_elements());
Array<OneD, NekDouble> Outarray;
int nPointsTot = pFields[0]->GetNpoints();
Array<OneD, NekDouble> grad0,grad1,grad2,wkSp;
{
nPointsTot = pFields[0]->Get1DScaledTotPoints(OneDptscale);
}
grad0 = Array<OneD, NekDouble> (nPointsTot);
int nadv = pV.num_elements();
{
AdvVel[0] = Array<OneD, NekDouble> (nPointsTot*(nadv+1));
for(int i = 0; i < nadv; ++i)
{
if(i)
{
AdvVel[i] = AdvVel[i-1]+nPointsTot;
}
pFields[0]->PhysInterp1DScaled(OneDptscale,pV[i],AdvVel[i]);
}
Outarray = AdvVel[nadv-1] + nPointsTot;
}
else
{
for(int i = 0; i < nadv; ++i)
{
AdvVel[i] = pV[i];
}
Outarray = pOutarray;
}
wkSp = Array<OneD, NekDouble> (nPointsTot);
switch(ndim)
{
case 1:
pFields[0]->PhysDeriv(pU,grad0);
break;
case 2:
{
grad1 = Array<OneD, NekDouble> (nPointsTot);
pFields[0]->PhysDeriv(pU,grad0,grad1);
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
}
Vmath::Vmul (nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,Outarray,1,Outarray,1);
{
pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,pOutarray);
}
}
break;
case 3:
grad1 = Array<OneD, NekDouble> (pFields[0]->GetNpoints());
grad2 = Array<OneD, NekDouble> (pFields[0]->GetNpoints());
{
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
Vmath::Vadd(nPointsTot,grad2,1,pOutarray,1,pOutarray,1);
}
{
pFields[0]->PhysDeriv(pV[pVelocityComponent],grad0,grad1);
pOutarray);
pFields[0]->HomogeneousBwdTrans(pOutarray,grad2);
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
}
else
{
Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
}
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
Outarray,1);
}
else
{
Outarray,1);
}
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad2,wkSp);
Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[2],1,Outarray,1,Outarray,1);
pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,grad2);
pFields[0]->HomogeneousFwdTrans(grad2,pOutarray);
}
else
{
Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,Outarray,1,grad0,1);
pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
}
}
{
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
}
else
{
Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
}
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
Outarray,1);
}
else
{
Outarray,1);
}
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad2,wkSp);
Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[2],1,Outarray,1,Outarray,1);
pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,pOutarray);
}
else
{
Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,Outarray,1,pOutarray,1);
}
}
{
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
pFields[0]->HomogeneousBwdTrans(grad0, pOutarray);
pFields[0]->DealiasedProd(pV[0], pOutarray, grad0,
pFields[0]->HomogeneousBwdTrans(grad1,pOutarray);
pFields[0]->DealiasedProd(pV[1], pOutarray, grad1,
pFields[0]->HomogeneousBwdTrans(grad2,pOutarray);
pFields[0]->DealiasedProd(pV[2], pOutarray, grad2,
pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
}
else
{
ASSERTL0(
false,
"Advection term calculation not implented or "
"possible with the current problem set up");
}
break;
default:
}
}