Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 using namespace std;
29 using namespace Nektar;
30 
31 int main(int argc, char *argv[])
32 {
33  int i,j;
34 
35  if(argc != 3)
36  {
37  fprintf(stderr,"Usage: ./StreamFunction2D file.xml file.fld\n");
38  exit(1);
39  }
40 
42  = LibUtilities::SessionReader::CreateInstance(argc, argv);
43 
44 
45  //----------------------------------------------
46  // Read in mesh from input file
47  string meshfile(argv[argc-2]);
48  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
49  //----------------------------------------------
50 
51  //----------------------------------------------
52  // Import field file.
53  string fieldfile(argv[argc-1]);
54  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
55  vector<vector<NekDouble> > fielddata;
56  LibUtilities::Import(fieldfile,fielddef,fielddata);
57 
58  //----------------------------------------------
59  // Define Expansion
60  int expdim = graphShPt->GetMeshDimension();
61  int nfields = fielddef[0]->m_fields.size();
62  int vorticitydim;
63 
64  vorticitydim = 2;
65 
66 
67 
68  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + vorticitydim);
69 
71 
72  switch(expdim)
73  {
74 
75  case 2:
76  {
77 
78  {
81  ::AllocateSharedPtr(vSession,graphShPt);
82  Exp[0] = Exp2D;
83 
84  FieldVar[0] = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(vSession, graphShPt, "v");
85 
86  for(i = 1; i < nfields + vorticitydim; ++i)
87  {
89  ::AllocateSharedPtr(*Exp2D);
90  }
91  }
92  }
93  break;
94 
95  default:
96  ASSERTL0(false,"The input file must be two-dimensional");
97  break;
98  }
99 
100  //----------------------------------------------
101  // Copy data from field file
102  for(j = 0; j < nfields; ++j)
103  {
104  for(unsigned int i = 0; i < fielddata.size(); ++i)
105  {
106  Exp[j]->ExtractDataToCoeffs(fielddef [i],
107  fielddata[i],
108  fielddef [i]->m_fields[j],
109  Exp[j]->UpdateCoeffs());
110  }
111  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
112  }
113  //----------------------------------------------
114 
115  int nq = Exp[0]->GetNpoints();
116 
117  Array<OneD, NekDouble> Uy(nq);
118  Array<OneD, NekDouble> Vx(nq);
119 
120  Array<OneD, NekDouble> StreamFunc(nq);
121  Array<OneD, NekDouble> Qz(nq);
122 
123  switch(expdim)
124  {
125 
126  case 2:
127  {
128 
129  {
130  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
131  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
132 
133  Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
134 
135  //The Vorticity is stored
136  Exp[3]->FwdTrans(Qz, Exp[3]->UpdateCoeffs());
137 
138  //We now calculate the Stream Function as the solution of the
139  //Poisson equation: Vorticity = - \nabla^2 StreamFunction
141  factor[StdRegions::eFactorLambda] = 0.0;
142 
143  Vmath::Smul(nq,-1.0,Qz,1,Qz,1);
144  Exp[4]->SetPhys(Qz);
145 
146  FieldVar[0]->HelmSolve(Qz, Exp[4]->UpdateCoeffs(), NullFlagList, factor);
147  }
148  }
149  break;
150 
151  default:
152  {
153  ASSERTL0(false,"The input file must be two-dimensional");
154  }
155  break;
156  }
157 
158  //-----------------------------------------------
159  // Write solution to file with additional computed fields
160  string fldfilename(argv[2]);
161  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
162  string endfile("_with_2DstremFunction.fld");
163  out += endfile;
164  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
165  = Exp[0]->GetFieldDefinitions();
166  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
167 
168  for(j = 0; j < nfields + vorticitydim; ++j)
169  {
170  for(i = 0; i < FieldDef.size(); ++i)
171  {
172  if (j >= nfields)
173  {
174  if(j == 4)
175  {
176  FieldDef[i]->m_fields.push_back("StreamFunc");
177  }
178  else
179  {
180  FieldDef[i]->m_fields.push_back("Qz");
181  }
182  }
183  else
184  {
185  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
186  }
187  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
188  }
189  }
190  LibUtilities::Write(out, FieldDef, FieldData);
191  //-----------------------------------------------
192 
193  return 0;
194 }
195 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:279
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
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:213
int main(int argc, char *argv[])
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
Definition: FieldIO.cpp:235
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:343
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
static FlagList NullFlagList
An empty flag list.