Nektar++
StreamFunction2D.cpp
Go to the documentation of this file.
1 /**
2  * This function calculate the vorticity vector starting from an .fld file.
3  * It is meant to be used with solutions produced by the incompressible Navier-Stokes solver.
4  * To use it with solutions coming form another solver further generalisations are required.
5  */
6 #include <cstdio>
7 #include <cstdlib>
8 
9 #include <MultiRegions/ExpList.h>
10 #include <MultiRegions/ExpList1D.h>
11 #include <MultiRegions/ExpList2D.h>
12 #include <MultiRegions/ExpList3D.h>
16 
17 
19 
20 #include <LocalRegions/MatrixKey.h>
27 
28 
29 
30 using namespace Nektar;
31 
32 int main(int argc, char *argv[])
33 {
34  int i,j;
35 
36  if(argc != 3)
37  {
38  fprintf(stderr,"Usage: ./StreamFunction2D file.xml file.fld\n");
39  exit(1);
40  }
41 
44 
45 
46  //----------------------------------------------
47  // Read in mesh from input file
48  string meshfile(argv[argc-2]);
50  //----------------------------------------------
51 
52  //----------------------------------------------
53  // Import field file.
54  string fieldfile(argv[argc-1]);
55  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
56  vector<vector<NekDouble> > fielddata;
57  LibUtilities::Import(fieldfile,fielddef,fielddata);
58 
59  //----------------------------------------------
60  // Define Expansion
61  int expdim = graphShPt->GetMeshDimension();
62  int nfields = fielddef[0]->m_fields.size();
63  int vorticitydim;
64 
65  vorticitydim = 2;
66 
67 
68 
69  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + vorticitydim);
70 
72 
73  switch(expdim)
74  {
75 
76  case 2:
77  {
78 
79  {
82  ::AllocateSharedPtr(vSession,graphShPt);
83  Exp[0] = Exp2D;
84 
85  FieldVar[0] = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(vSession, graphShPt, "v");
86 
87  for(i = 1; i < nfields + vorticitydim; ++i)
88  {
90  ::AllocateSharedPtr(*Exp2D);
91  }
92  }
93  }
94  break;
95 
96  default:
97  ASSERTL0(false,"The input file must be two-dimensional");
98  break;
99  }
100 
101  //----------------------------------------------
102  // Copy data from field file
103  for(j = 0; j < nfields; ++j)
104  {
105  for(unsigned int i = 0; i < fielddata.size(); ++i)
106  {
107  Exp[j]->ExtractDataToCoeffs(fielddef [i],
108  fielddata[i],
109  fielddef [i]->m_fields[j],
110  Exp[j]->UpdateCoeffs());
111  }
112  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
113  }
114  //----------------------------------------------
115 
116  int nq = Exp[0]->GetNpoints();
117 
118  Array<OneD, NekDouble> Uy(nq);
119  Array<OneD, NekDouble> Vx(nq);
120 
121  Array<OneD, NekDouble> StreamFunc(nq);
122  Array<OneD, NekDouble> Qz(nq);
123 
124  switch(expdim)
125  {
126 
127  case 2:
128  {
129 
130  {
131  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
132  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
133 
134  Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
135 
136  //The Vorticity is stored
137  Exp[3]->FwdTrans(Qz, Exp[3]->UpdateCoeffs());
138 
139  //We now calculate the Stream Function as the solution of the
140  //Poisson equation: Vorticity = - \nabla^2 StreamFunction
142  factor[StdRegions::eFactorLambda] = 0.0;
143 
144  Vmath::Smul(nq,-1.0,Qz,1,Qz,1);
145  Exp[4]->SetPhys(Qz);
146 
147  FieldVar[0]->HelmSolve(Qz, Exp[4]->UpdateCoeffs(), NullFlagList, factor);
148  }
149  }
150  break;
151 
152  default:
153  {
154  ASSERTL0(false,"The input file must be two-dimensional");
155  }
156  break;
157  }
158 
159  //-----------------------------------------------
160  // Write solution to file with additional computed fields
161  string fldfilename(argv[2]);
162  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
163  string endfile("_with_2DstremFunction.fld");
164  out += endfile;
165  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
166  = Exp[0]->GetFieldDefinitions();
167  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
168 
169  for(j = 0; j < nfields + vorticitydim; ++j)
170  {
171  for(i = 0; i < FieldDef.size(); ++i)
172  {
173  if (j >= nfields)
174  {
175  if(j == 4)
176  {
177  FieldDef[i]->m_fields.push_back("StreamFunc");
178  }
179  else
180  {
181  FieldDef[i]->m_fields.push_back("Qz");
182  }
183  }
184  else
185  {
186  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
187  }
188  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
189  }
190  }
191  LibUtilities::Write(out, FieldDef, FieldData);
192  //-----------------------------------------------
193 
194  return 0;
195 }
196 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:119
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
Definition: FieldIO.cpp:106
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
int main(int argc, char *argv[])
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:72
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
static FlagList NullFlagList
An empty flag list.