Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
FldAddWSS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <sstream>
#include <iostream>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/ContField2D.h>
Include dependency graph for FldAddWSS.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 15 of file FldAddWSS.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::LibUtilities::Import(), Nektar::SpatialDomains::MeshGraph::Read(), sign, Vmath::Smul(), Vmath::Svtsvtp(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vsqrt(), Vmath::Vvtvp(), Vmath::Vvtvvtp(), and Nektar::LibUtilities::Write().

{
int i,j;
int surfID, nfiles, nStart;
if(argc != 5)
{
fprintf(stderr,"Usage: FldAddWSS meshfile nfiles infld BoundaryID\n");
exit(1);
}
surfID = boost::lexical_cast<int>(argv[argc - 1]);
nfiles = boost::lexical_cast<int>(argv[2]);
vector<string> infiles(nfiles), outfiles(nfiles);
/* format of files if 1 file: name.fld. Output: name_wss.fld
format of files if n checkpoint files: name_tn.fld. Output: name_tn_wss.fld, n=0,1,2...
*/
if (nfiles == 1)
{
infiles[0] = argv[3];
string basename = argv[3];
basename = basename.substr(0, basename.find_last_of("."));
stringstream filename2;
filename2 << basename << "_wss.fld";
filename2 >> outfiles[0];
}
else {
string basename = argv[3];
basename = basename.substr(basename.find_last_of("t")+1, basename.find_last_of(".")-basename.find_last_of("t"));
stringstream filename3;
filename3 << basename;
filename3 >> nStart;
for (int i = 0; i< nfiles; ++i)
{
basename = argv[3];
string extension = ".fld";
basename = basename.substr(0, basename.find_first_of("_"));
stringstream filename, filename2;
filename << basename << "_t" << i + nStart << extension;
filename >> infiles[i];
filename2 << basename << "_t" << i +nStart << "_wss.fld";
filename2 >> outfiles[i];
}
}
argv[argc - 1] = argv[argc - 2];
argv[argc - 3] = argv[argc - 2];
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(argc, argv);
//----------------------------------------------
// Read in mesh from input file
string meshfile(argv[argc-4]);
SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);
//----------------------------------------------
//----------------------------------------------
//Get viscosity from file
NekDouble m_kinvis;
m_kinvis = vSession->GetParameter("Kinvis");
//----------------------------------------------
//----------------------------------------------
// Import first field file.
string fieldfile(argv[3]);
vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
vector<vector<NekDouble> > fielddata;
LibUtilities::Import(fieldfile,fielddef,fielddata);
//----------------------------------------------
//----------------------------------------------
// Define Expansion
int expdim = graphShPt->GetMeshDimension();
int nfields = fielddef[0]->m_fields.size()-1;
int addfields = (nfields == 3)? 4:3;
int nstress = (nfields == 3)? 6:3;
Array<OneD, MultiRegions::ExpListSharedPtr> vel(nfields), shear(addfields);
MultiRegions::AssemblyMapCGSharedPtr m_locToGlobalMap;
switch(expdim)
{
case 1:
{
ASSERTL0(false,"Expansion dimension not recognised");
}
break;
case 2:
{
ASSERTL0(false,"Not implemented in 2D");
int i=0;
MultiRegions::ContField2DSharedPtr cont2D;
cont2D = MemoryManager<MultiRegions::ContField2D>
::AllocateSharedPtr(vSession,graphShPt,vSession->GetVariable(i));
vel[0] = cont2D;
for(int i = 1; i < nfields; ++i)
{
vel[i] = MemoryManager<MultiRegions::ContField2D>
::AllocateSharedPtr(*cont2D);
}
}
break;
case 3:
{
MultiRegions::ContField3DSharedPtr firstfield =
MemoryManager<MultiRegions::ContField3D>
::AllocateSharedPtr(vSession, graphShPt,
vSession->GetVariable(0));
m_locToGlobalMap = firstfield->GetLocalToGlobalMap();
vel[0] = firstfield;
for(int i = 1; i < nfields; ++i)
{
vel[i] = MemoryManager<MultiRegions::ContField3D>
::AllocateSharedPtr(*firstfield, graphShPt,
vSession->GetVariable(i));
}
for(int i = 0; i < addfields; ++i)
{
shear[i] = MemoryManager<MultiRegions::ContField3D>
::AllocateSharedPtr(*firstfield, graphShPt,
vSession->GetVariable(0));
}
}
break;
default:
ASSERTL0(false,"Expansion dimension not recognised");
break;
}
//----------------------------------------------
// Define arrays for WSS (global array), stress components and gradients (local arrays)
int n, cnt, elmtid, nq, offset, nt, boundary, nfq;
Array<OneD, Array<OneD, NekDouble> > grad(nfields*nfields);
Array<OneD, Array<OneD, NekDouble> > stress(nstress), fstress(nstress);
Array<OneD, Array<OneD, NekDouble> > values(addfields);
// Set up mapping from Boundary condition to element details.
StdRegions::StdExpansionSharedPtr elmt;
StdRegions::StdExpansion2DSharedPtr bc;
Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> > BndExp(addfields);
for (int fileNo = 0; fileNo < nfiles ; ++fileNo)
{
//----------------------------------------------
// Import field file.
fielddef.clear();
fielddata.clear();
LibUtilities::Import(infiles[fileNo],fielddef,fielddata);
//----------------------------------------------
//----------------------------------------------
// Copy data from field file
for(int j = 0; j < nfields; ++j)
{
for(int i = 0; i < fielddata.size(); ++i)
{
vel[j]->ExtractDataToCoeffs(fielddef[i],
fielddata[i],
fielddef[i]->m_fields[j],
vel[j]->UpdateCoeffs());
}
vel[j]->BwdTrans(vel[j]->GetCoeffs(),vel[j]->UpdatePhys());
}
for(int j = 0; j < addfields; ++j)
{
for(int i = 0; i < fielddata.size(); ++i)
{
shear[j]->ExtractDataToCoeffs(fielddef[i],
fielddata[i],
fielddef[i]->m_fields[0],
shear[j]->UpdateCoeffs());
}
shear[j]->BwdTrans(shear[j]->GetCoeffs(),shear[j]->UpdatePhys());
}
//----------------------------------------------
//----------------------------------------------
// Compute WSS for each element on boundary
// Define total quadrature points, and velocity fields.
nt = shear[0]->GetNpoints();
Array<OneD, const NekDouble> U(nt),V(nt),W(nt);
for(int j = 0; j < addfields; ++j)
{
values[j] = Array<OneD, NekDouble>(nt);
}
shear[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
//get boundary expansions for each field
for(int j = 0; j < addfields; ++j)
{
BndExp[j] = shear[j]->GetBndCondExpansions();
}
// loop over the types of boundary conditions
for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
{
// identify boundary which the user wanted
if(n == surfID)
{
for(int i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
{
// find element and face of this expansion.
elmtid = BoundarytoElmtID[cnt];
elmt = shear[0]->GetExp(elmtid);
nq = elmt->GetTotPoints();
offset = shear[0]->GetPhys_Offset(elmtid);
// Initialise local arrays for the velocity gradients, and stress components
// size of total number of quadrature points for each element (hence local).
for(int j = 0; j < nfields*nfields; ++j)
{
grad[j] = Array<OneD, NekDouble>(nq);
}
for(int j = 0; j < nstress; ++j)
{
stress[j] = Array<OneD, NekDouble>(nq);
}
if(nfields == 2)
{
//Not implemented in 2D.
}
else
{
// Get face 2D expansion from element expansion
bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
nfq= bc->GetTotPoints();
//identify boundary of element looking at.
boundary = BoundarytoTraceID[cnt];
//Get face normals
const SpatialDomains::GeomFactorsSharedPtr m_metricinfo=bc->GetMetricInfo();
const Array<OneD, const Array<OneD, NekDouble> > normals
= elmt->GetFaceNormal(boundary);
// initialise arrays
for(int j = 0; j < nstress; ++j)
{
fstress[j] = Array<OneD, NekDouble>(nfq);
}
Array<OneD, NekDouble> values2(nfq), S(nfq), Sx(nfq), Sy(nfq), Sz(nfq);
//Extract Velocities
U = vel[0]->GetPhys() + offset;
V = vel[1]->GetPhys() + offset;
W = vel[2]->GetPhys() + offset;
//Compute gradients (velocity correction scheme method)
elmt->PhysDeriv(U,grad[0],grad[1],grad[2]);
elmt->PhysDeriv(V,grad[3],grad[4],grad[5]);
elmt->PhysDeriv(W,grad[6],grad[7],grad[8]);
//Compute stress component terms
// t_xx = 2.mu.Ux
Vmath::Smul (nq,(2*m_kinvis),grad[0],1,stress[0],1);
// tyy = 2.mu.Vy
Vmath::Smul (nq,(2*m_kinvis),grad[4],1,stress[1],1);
// tzz = 2.mu.Wz
Vmath::Smul (nq,(2*m_kinvis),grad[8],1,stress[2],1);
// txy = mu.(Uy+Vx)
Vmath::Vadd (nq,grad[1],1,grad[3],1,stress[3],1);
Vmath::Smul (nq,m_kinvis,stress[3],1,stress[3],1);
// txz = mu.(Uz+Wx)
Vmath::Vadd (nq,grad[2],1,grad[6],1,stress[4],1);
Vmath::Smul (nq,m_kinvis,stress[4],1,stress[4],1);
// tyz = mu.(Vz+Wy)
Vmath::Vadd (nq,grad[5],1,grad[7],1,stress[5],1);
Vmath::Smul (nq,m_kinvis,stress[5],1,stress[5],1);
// Get face stress values.
for(j = 0; j < nstress; ++j)
{
elmt->GetFacePhysVals(boundary,bc,stress[j],fstress[j]);
}
//calcuate wss, and update velocity coefficients in the elemental boundary expansion
for (int j = 0; j< addfields; j++)
{
values[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[j][n]->GetCoeff_Offset(i);
}
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
// Sx
Vmath::Vvtvvtp(nfq,normals[0],1,fstress[0],1,
normals[1],1,fstress[3],1,Sx,1);
Vmath::Vvtvp (nfq,normals[2],1,fstress[4],1,Sx,1,Sx,1);
// Sy
Vmath::Vvtvvtp(nfq,normals[0],1,fstress[3],1,
normals[1],1,fstress[1],1,Sy,1);
Vmath::Vvtvp (nfq,normals[2],1,fstress[5],1,Sy,1,Sy,1);
// Sz
Vmath::Vvtvvtp(nfq,normals[0],1,fstress[4],1,
normals[1],1,fstress[5],1,Sz,1);
Vmath::Vvtvp (nfq,normals[2],1,fstress[2],1,Sz,1,Sz,1);
}
else
{
// Sx
Vmath::Svtsvtp(nfq,normals[0][0],fstress[0],1,
normals[1][0],fstress[3],1,Sx,1);
Vmath::Svtvp(nfq,normals[2][0],fstress[4],1,Sx,1,Sx,1);
// Sy
Vmath::Svtsvtp(nfq,normals[0][0],fstress[3],1,
normals[1][0],fstress[1],1,Sy,1);
Vmath::Svtvp(nfq,normals[2][0],fstress[5],1,Sy,1,Sy,1);
// Sz
Vmath::Svtsvtp(nfq,normals[0][0],fstress[4],1,
normals[1][0],fstress[5],1,Sz,1);
Vmath::Svtvp(nfq,normals[2][0],fstress[2],1,Sz,1,Sz,1);
}
// T = T - (T.n)n
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vvtvvtp(nfq,normals[0],1,Sx,1,
normals[1],1, Sy,1,values2,1);
Vmath::Vvtvp (nfq,normals[2],1, Sz,1,values2,1,values2,1);
// values2 = (T.n)n
Vmath::Smul(nfq, -1.0, values2, 1, values2, 1);
//Sx
Vmath::Vvtvp(nfq,normals[0],1, values2,1,Sx,1,Sx,1);
bc->FwdTrans(Sx, values[0]);
//Sy
Vmath::Vvtvp(nfq,normals[1],1, values2,1,Sy,1,Sy,1);
bc->FwdTrans(Sy, values[1]);
//Sz
Vmath::Vvtvp(nfq,normals[2],1, values2,1,Sz,1,Sz,1);
bc->FwdTrans(Sz, values[2]);
}
else
{
Vmath::Svtsvtp(nfq,normals[0][0],Sx,1,
normals[1][0],Sy,1,values2,1);
Vmath::Svtvp(nfq,normals[2][0],Sz,1,values2,1,values2,1);
Vmath::Smul(nfq, -1.0, values2, 1, values2, 1);
//Sx
Vmath::Svtvp(nfq,normals[0][0],values2,1,Sx,1,Sx,1);
bc->FwdTrans(Sx, values[0]);
//Sy
Vmath::Svtvp(nfq,normals[1][0],values2,1,Sy,1,Sy,1);
bc->FwdTrans(Sy, values[1]);
//Sz
Vmath::Svtvp(nfq,normals[2][0],values2,1,Sz,1,Sz,1);
bc->FwdTrans(Sz, values[2]);
}
//Tw
Vmath::Vvtvvtp(nfq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, S, 1);
Vmath::Vvtvp(nfq, Sz, 1, Sz, 1, S, 1, S, 1);
Vmath::Vsqrt(nfq, S, 1, S, 1);
bc->FwdTrans(S, values[3]);
}
}
}
else
{
cnt += BndExp[0][n]->GetExpSize();
}
}
for(int j = 0; j < addfields; ++j)
{
int ncoeffs = shear[j]->GetNcoeffs();
Array<OneD, NekDouble> output(ncoeffs);
output=shear[j]->UpdateCoeffs();
int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
Array<OneD, NekDouble> outarray(nGlobal,0.0);
int bndcnt=0;
const Array<OneD,const int>& map = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsMap();
NekDouble sign;
for(int i = 0; i < BndExp[j].num_elements(); ++i)
{
if(i==surfID)
{
const Array<OneD,const NekDouble>& coeffs = BndExp[j][i]->GetCoeffs();
for(int k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
{
sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
outarray[map[bndcnt++]] = sign * coeffs[k];
}
}
else
{
bndcnt += BndExp[j][i]->GetNcoeffs();
}
}
m_locToGlobalMap->GlobalToLocal(outarray,output);
}
//-----------------------------------------------
// Write solution to file with additional computed fields
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
= shear[0]->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
vector<string > outname;
if(addfields == 3)
{
// Not implemented in 2D.
}
else
{
outname.push_back("Tx");
outname.push_back("Ty");
outname.push_back("Tz");
outname.push_back("Tw");
}
for(j = 0; j < nfields + addfields; ++j)
{
for(i = 0; i < FieldDef.size(); ++i)
{
if (j < nfields)
{
FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
vel[j]->AppendFieldData(FieldDef[i], FieldData[i]);
}
else
{
FieldDef[i]->m_fields.push_back(outname[j-nfields]);
shear[j-nfields]->AppendFieldData(FieldDef[i], FieldData[i]);
}
}
}
LibUtilities::Write(outfiles[fileNo], FieldDef, FieldData);
FieldDef.clear();
FieldData.clear();
//-----------------------------------------------
}
return 0;
}
Generated on Sun Mar 15 2015 00:11:39 for Nektar++ by   doxygen 1.8.1.2