Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FldAddVort.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 
4 #include <MultiRegions/ExpList.h>
5 #include <MultiRegions/ExpList1D.h>
6 #include <MultiRegions/ExpList2D.h>
7 #include <MultiRegions/ExpList3D.h>
8 #include <MultiRegions/ExpList2DHomogeneous1D.h>
9 #include <MultiRegions/ExpList3DHomogeneous1D.h>
10 #include <MultiRegions/ExpList3DHomogeneous2D.h>
11 
12 using namespace Nektar;
13 
14 int main(int argc, char *argv[])
15 {
16  int i,j;
17 
18  if(argc != 4)
19  {
20  fprintf(stderr,"Usage: FldAddVort meshfile infld outfld\n");
21  exit(1);
22  }
23 
24  LibUtilities::SessionReaderSharedPtr vSession
25  = LibUtilities::SessionReader::CreateInstance(argc, argv);
26 
27 
28  //----------------------------------------------
29  // Read in mesh from input file
30  string meshfile(argv[argc-3]);
31  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
32  //----------------------------------------------
33 
34  //----------------------------------------------
35  // Import field file.
36  string fieldfile(argv[argc-2]);
37  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
38  vector<vector<NekDouble> > fielddata;
39  LibUtilities::Import(fieldfile,fielddef,fielddata);
40  bool useFFT = false;
41  bool dealiasing = false;
42  //----------------------------------------------
43 
44  //----------------------------------------------
45  // Define Expansion
46  int expdim = graphShPt->GetMeshDimension();
47  int nfields = fielddef[0]->m_fields.size();
48  int addfields = (nfields == 4)? 3:1;
49  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + addfields);
50 
51  switch(expdim)
52  {
53  case 1:
54  {
55  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
56 
57  if(fielddef[0]->m_numHomogeneousDir == 1)
58  {
59  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
60 
61  // Define Homogeneous expansion
62  //int nplanes = fielddef[0]->m_numModes[1];
63  int nplanes;
64  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
65 
66  // choose points to be at evenly spaced points at
67  // nplanes points
68  const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
69  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
70  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
71 
72  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
73  Exp[0] = Exp2DH1;
74 
75  for(i = 1; i < nfields; ++i)
76  {
77  Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
78  }
79  }
80  else if(fielddef[0]->m_numHomogeneousDir == 2)
81  {
82  MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
83 
84  // Define Homogeneous expansion
85  //int nylines = fielddef[0]->m_numModes[1];
86  //int nzlines = fielddef[0]->m_numModes[2];
87 
88  int nylines;
89  int nzlines;
90  vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
91  vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
92 
93  // choose points to be at evenly spaced points at
94  // nplanes points
95  const LibUtilities::PointsKey PkeyY(nylines,LibUtilities::ePolyEvenlySpaced);
96  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
97 
98  const LibUtilities::PointsKey PkeyZ(nzlines,LibUtilities::ePolyEvenlySpaced);
99  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
100 
101  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
102  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
103 
104  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
105  Exp[0] = Exp3DH2;
106 
107  for(i = 1; i < nfields; ++i)
108  {
109  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
110  }
111  }
112  else
113  {
114  MultiRegions::ExpList1DSharedPtr Exp1D;
115  Exp1D = MemoryManager<MultiRegions::ExpList1D>
116  ::AllocateSharedPtr(vSession,graphShPt);
117  Exp[0] = Exp1D;
118  for(i = 1; i < nfields + addfields; ++i)
119  {
120  Exp[i] = MemoryManager<MultiRegions::ExpList1D>
121  ::AllocateSharedPtr(*Exp1D);
122  }
123  }
124  }
125  break;
126  case 2:
127  {
128  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
129 
130  if(fielddef[0]->m_numHomogeneousDir == 1)
131  {
132  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
133 
134  // Define Homogeneous expansion
135  //int nplanes = fielddef[0]->m_numModes[2];
136 
137  int nplanes;
138  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
139 
140  // choose points to be at evenly spaced points at
141  // nplanes points
142  const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
143  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
144  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
145 
146  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
147  Exp[0] = Exp3DH1;
148 
149  for(i = 1; i < nfields + addfields; ++i)
150  {
151  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
152 
153  }
154  }
155  else
156  {
157  MultiRegions::ExpList2DSharedPtr Exp2D;
158  Exp2D = MemoryManager<MultiRegions::ExpList2D>
159  ::AllocateSharedPtr(vSession,graphShPt);
160  Exp[0] = Exp2D;
161 
162  for(i = 1; i < nfields + addfields; ++i)
163  {
164  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
165  ::AllocateSharedPtr(*Exp2D);
166  }
167  }
168  }
169  break;
170  case 3:
171  {
172  MultiRegions::ExpList3DSharedPtr Exp3D;
173  Exp3D = MemoryManager<MultiRegions::ExpList3D>
174  ::AllocateSharedPtr(vSession,graphShPt);
175  Exp[0] = Exp3D;
176 
177  for(i = 1; i < nfields + addfields; ++i)
178  {
179  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
180  ::AllocateSharedPtr(*Exp3D);
181  }
182  }
183  break;
184  default:
185  ASSERTL0(false,"Expansion dimension not recognised");
186  break;
187  }
188  //----------------------------------------------
189 
190  //----------------------------------------------
191  // Copy data from field file
192  for(j = 0; j < nfields; ++j)
193  {
194  for(int i = 0; i < fielddata.size(); ++i)
195  {
196  Exp[j]->ExtractDataToCoeffs(fielddef [i],
197  fielddata[i],
198  fielddef [i]->m_fields[j],
199  Exp[j]->UpdateCoeffs());
200  }
201  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
202  }
203  //----------------------------------------------
204 
205  //----------------------------------------------
206  // Compute gradients of fields
207  ASSERTL0(nfields >= 3, "Need two fields (u,v) to add reentricity");
208  int nq = Exp[0]->GetNpoints();
209  Array<OneD, Array<OneD, NekDouble> > grad(nfields*nfields);
210  Array<OneD, Array<OneD, NekDouble> > outfield(addfields);
211 
212  for(i = 0; i < nfields*nfields; ++i)
213  {
214  grad[i] = Array<OneD, NekDouble>(nq);
215  }
216 
217  for(i = 0; i < addfields; ++i)
218  {
219  outfield[i] = Array<OneD, NekDouble>(nq);
220  }
221 
222  // Calculate Gradient & Vorticity
223  if(nfields == 3)
224  {
225  for(i = 0; i < nfields; ++i)
226  {
227  Exp[i]->PhysDeriv(Exp[i]->GetPhys(), grad[i*nfields],grad[i*nfields+1]);
228  }
229  // W_z = Vx - Uy
230  Vmath::Vsub(nq,grad[1*nfields+0],1,grad[0*nfields+1],1,outfield[0],1);
231  }
232  else
233  {
234  for(i = 0; i < nfields; ++i)
235  {
236  Exp[i]->PhysDeriv(Exp[i]->GetPhys(), grad[i*nfields],grad[i*nfields+1],grad[i*nfields+2]);
237  }
238 
239  // W_x = Wy - Vz
240  Vmath::Vsub(nq,grad[2*nfields+1],1,grad[1*nfields+2],1,outfield[0],1);
241  // W_y = Uz - Wx
242  Vmath::Vsub(nq,grad[0*nfields+2],1,grad[2*nfields+0],1,outfield[1],1);
243  // W_z = Vx - Uy
244  Vmath::Vsub(nq,grad[1*nfields+0],1,grad[0*nfields+1],1,outfield[2],1);
245  }
246 
247  for (i = 0; i < addfields; ++i)
248  {
249  Exp[nfields + i]->FwdTrans(outfield[i], Exp[nfields+i]->UpdateCoeffs());
250  }
251 
252  //-----------------------------------------------
253  // Write solution to file with additional computed fields
254  string out(argv[argc-1]);
255  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
256  = Exp[0]->GetFieldDefinitions();
257  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
258 
259  vector<string > outname;
260 
261  if(addfields == 1)
262  {
263  outname.push_back("W_z");
264  }
265  else
266  {
267  outname.push_back("W_x");
268  outname.push_back("W_y");
269  outname.push_back("W_z");
270  }
271 
272  for(j = 0; j < nfields + addfields; ++j)
273  {
274  for(i = 0; i < FieldDef.size(); ++i)
275  {
276  if (j >= nfields)
277  {
278  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
279  }
280  else
281  {
282  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
283  }
284  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
285  }
286  }
287  LibUtilities::Write(out, FieldDef, FieldData);
288  //-----------------------------------------------
289 
290  return 0;
291 }
292 
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2