Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
FldAddVort.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
Include dependency graph for FldAddVort.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 14 of file FldAddVort.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::Import(), Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Vsub(), and Nektar::LibUtilities::Write().

{
int i,j;
if(argc != 4)
{
fprintf(stderr,"Usage: FldAddVort meshfile infld outfld\n");
exit(1);
}
= LibUtilities::SessionReader::CreateInstance(argc, argv);
//----------------------------------------------
// Read in mesh from input file
string meshfile(argv[argc-3]);
SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
//----------------------------------------------
//----------------------------------------------
// Import field file.
string fieldfile(argv[argc-2]);
vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
vector<vector<NekDouble> > fielddata;
LibUtilities::Import(fieldfile,fielddef,fielddata);
bool useFFT = false;
bool dealiasing = false;
//----------------------------------------------
//----------------------------------------------
// Define Expansion
int expdim = graphShPt->GetMeshDimension();
int nfields = fielddef[0]->m_fields.size();
int addfields = (nfields == 4)? 3:1;
Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + addfields);
switch(expdim)
{
case 1:
{
ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
if(fielddef[0]->m_numHomogeneousDir == 1)
{
// Define Homogeneous expansion
//int nplanes = fielddef[0]->m_numModes[1];
int nplanes;
vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
// choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
Exp[0] = Exp2DH1;
for(i = 1; i < nfields; ++i)
{
}
}
else if(fielddef[0]->m_numHomogeneousDir == 2)
{
// Define Homogeneous expansion
//int nylines = fielddef[0]->m_numModes[1];
//int nzlines = fielddef[0]->m_numModes[2];
int nylines;
int nzlines;
vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
// choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
Exp[0] = Exp3DH2;
for(i = 1; i < nfields; ++i)
{
}
}
else
{
::AllocateSharedPtr(vSession,graphShPt);
Exp[0] = Exp1D;
for(i = 1; i < nfields + addfields; ++i)
{
}
}
}
break;
case 2:
{
ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
if(fielddef[0]->m_numHomogeneousDir == 1)
{
// Define Homogeneous expansion
//int nplanes = fielddef[0]->m_numModes[2];
int nplanes;
vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
// choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
Exp[0] = Exp3DH1;
for(i = 1; i < nfields + addfields; ++i)
{
}
}
else
{
::AllocateSharedPtr(vSession,graphShPt);
Exp[0] = Exp2D;
for(i = 1; i < nfields + addfields; ++i)
{
}
}
}
break;
case 3:
{
::AllocateSharedPtr(vSession,graphShPt);
Exp[0] = Exp3D;
for(i = 1; i < nfields + addfields; ++i)
{
}
}
break;
default:
ASSERTL0(false,"Expansion dimension not recognised");
break;
}
//----------------------------------------------
//----------------------------------------------
// Copy data from field file
for(j = 0; j < nfields; ++j)
{
for(int i = 0; i < fielddata.size(); ++i)
{
Exp[j]->ExtractDataToCoeffs(fielddef [i],
fielddata[i],
fielddef [i]->m_fields[j],
Exp[j]->UpdateCoeffs());
}
Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
}
//----------------------------------------------
//----------------------------------------------
// Compute gradients of fields
ASSERTL0(nfields >= 3, "Need two fields (u,v) to add reentricity");
int nq = Exp[0]->GetNpoints();
Array<OneD, Array<OneD, NekDouble> > grad(nfields*nfields);
Array<OneD, Array<OneD, NekDouble> > outfield(addfields);
for(i = 0; i < nfields*nfields; ++i)
{
grad[i] = Array<OneD, NekDouble>(nq);
}
for(i = 0; i < addfields; ++i)
{
outfield[i] = Array<OneD, NekDouble>(nq);
}
// Calculate Gradient & Vorticity
if(nfields == 3)
{
for(i = 0; i < nfields; ++i)
{
Exp[i]->PhysDeriv(Exp[i]->GetPhys(), grad[i*nfields],grad[i*nfields+1]);
}
// W_z = Vx - Uy
Vmath::Vsub(nq,grad[1*nfields+0],1,grad[0*nfields+1],1,outfield[0],1);
}
else
{
for(i = 0; i < nfields; ++i)
{
Exp[i]->PhysDeriv(Exp[i]->GetPhys(), grad[i*nfields],grad[i*nfields+1],grad[i*nfields+2]);
}
// W_x = Wy - Vz
Vmath::Vsub(nq,grad[2*nfields+1],1,grad[1*nfields+2],1,outfield[0],1);
// W_y = Uz - Wx
Vmath::Vsub(nq,grad[0*nfields+2],1,grad[2*nfields+0],1,outfield[1],1);
// W_z = Vx - Uy
Vmath::Vsub(nq,grad[1*nfields+0],1,grad[0*nfields+1],1,outfield[2],1);
}
for (i = 0; i < addfields; ++i)
{
Exp[nfields + i]->FwdTrans(outfield[i], Exp[nfields+i]->UpdateCoeffs());
}
//-----------------------------------------------
// Write solution to file with additional computed fields
string out(argv[argc-1]);
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
= Exp[0]->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
vector<string > outname;
if(addfields == 1)
{
outname.push_back("W_z");
}
else
{
outname.push_back("W_x");
outname.push_back("W_y");
outname.push_back("W_z");
}
for(j = 0; j < nfields + addfields; ++j)
{
for(i = 0; i < FieldDef.size(); ++i)
{
if (j >= nfields)
{
FieldDef[i]->m_fields.push_back(outname[j-nfields]);
}
else
{
FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
}
Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
}
}
LibUtilities::Write(out, FieldDef, FieldData);
//-----------------------------------------------
return 0;
}