Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CalcGrad.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>
13 #include <MultiRegions/ExpList2DHomogeneous1D.h>
14 #include <MultiRegions/ExpList3DHomogeneous1D.h>
15 #include <MultiRegions/ExpList3DHomogeneous2D.h>
16 
17 using namespace Nektar;
18 
19 int main(int argc, char *argv[])
20 {
21  int i,j;
22 
23  if(argc != 4)
24  {
25  fprintf(stderr,"Usage: ./CalcGrad ifield file.xml file.fld\n");
26  exit(1);
27  }
28 
29  LibUtilities::SessionReaderSharedPtr vSession
30  = LibUtilities::SessionReader::CreateInstance(argc-1, argv+1);
31 
32 
33  int ifield = atoi(argv[argc-3]);
34 
35  //----------------------------------------------
36  // Read in mesh from input file
37  string meshfile(argv[argc-2]);
38  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
39  //----------------------------------------------
40 
41  //----------------------------------------------
42  // Import field file.
43  string fieldfile(argv[argc-1]);
44  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
45  vector<vector<NekDouble> > fielddata;
46  LibUtilities::Import(fieldfile,fielddef,fielddata);
47  bool useFFT = false;
48  bool dealiasing = false;
49  //----------------------------------------------
50 
51  //----------------------------------------------
52  // Define Expansion
53  int expdim = graphShPt->GetMeshDimension();
54  int nfields = fielddef[0]->m_fields.size();
55  int graddim;
56  if(expdim == 1)
57  {
58  if(fielddef[0]->m_numHomogeneousDir == 2)//3D Homogeneous 2D
59  {
60  graddim = 3;
61  }
62  else // 1D
63  {
64  graddim = 1;
65  }
66  }
67  else if(expdim ==2)
68  {
69  if(fielddef[0]->m_numHomogeneousDir == 1)// 3D Homogeneous 1D
70  {
71  graddim = 3;
72  }
73  else //2D
74  {
75  graddim = 2;
76  }
77 
78  }
79  else // Full 3D
80  {
81  graddim = 3;
82  }
83 
84 
85  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + graddim);
86 
87  switch(expdim)
88  {
89  case 1:
90  {
91  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
92 
93  if(fielddef[0]->m_numHomogeneousDir == 1)
94  {
95  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
96 
97  // Define Homogeneous expansion
98  int nplanes = fielddef[0]->m_numModes[1];
99 
100  // choose points to be at evenly spaced points at
101  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
102  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
103  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
104 
105  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
106  Exp[0] = Exp2DH1;
107 
108  for(i = 1; i < nfields + graddim; ++i)
109  {
110  Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
111  }
112  }
113  else if(fielddef[0]->m_numHomogeneousDir == 2)
114  {
115  MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
116 
117  // Define Homogeneous expansion
118  int nylines = fielddef[0]->m_numModes[1];
119  int nzlines = fielddef[0]->m_numModes[2];
120 
121  // choose points to be at evenly spaced points at
122  const LibUtilities::PointsKey PkeyY(nylines+1,LibUtilities::ePolyEvenlySpaced);
123  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
124 
125  const LibUtilities::PointsKey PkeyZ(nzlines+1,LibUtilities::ePolyEvenlySpaced);
126  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
127 
128  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
129  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
130 
131  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
132  Exp[0] = Exp3DH2;
133 
134  for(i = 1; i < nfields + graddim; ++i)
135  {
136  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
137  }
138  }
139  else
140  {
141  MultiRegions::ExpList1DSharedPtr Exp1D;
142  Exp1D = MemoryManager<MultiRegions::ExpList1D>
143  ::AllocateSharedPtr(vSession,graphShPt);
144  Exp[0] = Exp1D;
145  for(i = 1; i < nfields + graddim; ++i)
146  {
147  Exp[i] = MemoryManager<MultiRegions::ExpList1D>
148  ::AllocateSharedPtr(*Exp1D);
149  }
150  }
151  }
152  break;
153  case 2:
154  {
155  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
156 
157  if(fielddef[0]->m_numHomogeneousDir == 1)
158  {
159  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
160 
161  // Define Homogeneous expansion
162  int nplanes = fielddef[0]->m_numModes[2];
163 
164  // choose points to be at evenly spaced points at
165  // nplanes + 1 points
166  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
167  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
168  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
169 
170  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
171  Exp[0] = Exp3DH1;
172 
173  for(i = 1; i < nfields + graddim; ++i)
174  {
175  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
176  }
177  }
178  else
179  {
180  MultiRegions::ExpList2DSharedPtr Exp2D;
181  Exp2D = MemoryManager<MultiRegions::ExpList2D>
182  ::AllocateSharedPtr(vSession,graphShPt);
183  Exp[0] = Exp2D;
184 
185  for(i = 1; i < nfields + graddim; ++i)
186  {
187  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
188  ::AllocateSharedPtr(*Exp2D);
189  }
190  }
191  }
192  break;
193  case 3:
194  {
195  MultiRegions::ExpList3DSharedPtr Exp3D;
196  Exp3D = MemoryManager<MultiRegions::ExpList3D>
197  ::AllocateSharedPtr(vSession,graphShPt);
198  Exp[0] = Exp3D;
199 
200  for(i = 1; i < nfields + graddim; ++i)
201  {
202  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
203  ::AllocateSharedPtr(*Exp3D);
204  }
205  }
206  break;
207  default:
208  ASSERTL0(false,"Expansion dimension not recognised");
209  break;
210  }
211 
212  //----------------------------------------------
213  // Copy data from field file
214  for(j = 0; j < nfields; ++j)
215  {
216  for(unsigned int i = 0; i < fielddata.size(); ++i)
217  {
218  Exp[j]->ExtractDataToCoeffs(fielddef [i],
219  fielddata[i],
220  fielddef [i]->m_fields[j],
221  Exp[j]->UpdateCoeffs());
222  }
223  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
224  }
225  //----------------------------------------------
226 
227  int nq = Exp[0]->GetNpoints();
228 
229  Array<OneD, NekDouble> Dy(nq);
230  Array<OneD, NekDouble> Dz(nq);
231  Array<OneD, NekDouble> Dx(nq);
232 
233  switch(expdim)
234  {
235  case 1:
236  {
237  if(fielddef[0]->m_numHomogeneousDir == 1)
238  {
239  ASSERTL0(false,"Not implemented yet");
240  }
241  else if(fielddef[0]->m_numHomogeneousDir == 2)
242  {
243  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
244  Exp[ifield]->GetPhys(),Dx);
245  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
246  Exp[ifield]->GetPhys(),Dy);
247  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
248  Exp[ifield]->GetPhys(),Dz);
249  Exp[nfields ]->FwdTrans(Dx,Exp[nfields ]->UpdateCoeffs());
250  Exp[nfields+1]->FwdTrans(Dy,Exp[nfields+1]->UpdateCoeffs());
251  Exp[nfields+2]->FwdTrans(Dz,Exp[nfields+2]->UpdateCoeffs());
252 
253  }
254  else
255  {
256  ASSERTL0(false,"Not implemented yet");
257  }
258  }
259  break;
260  case 2:
261  {
262  if(fielddef[0]->m_numHomogeneousDir == 1)
263  {
264  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[ifield]->GetPhys(),Dx);
265  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[ifield]->GetPhys(),Dy);
266  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[ifield]->GetPhys(),Dz);
267  Exp[nfields ]->FwdTrans(Dx,Exp[nfields ]->UpdateCoeffs());
268  Exp[nfields+1]->FwdTrans(Dy,Exp[nfields+1]->UpdateCoeffs());
269  Exp[nfields+2]->FwdTrans(Dz,Exp[nfields+2]->UpdateCoeffs());
270  }
271  else
272  {
273  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[ifield]->GetPhys(),Dx);
274  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[ifield]->GetPhys(),Dy);
275  Exp[nfields ]->FwdTrans(Dx,Exp[nfields ]->UpdateCoeffs());
276  Exp[nfields+1]->FwdTrans(Dy,Exp[nfields+1]->UpdateCoeffs());
277  }
278  }
279  break;
280  case 3:
281  {
282  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[ifield]->GetPhys(),Dx);
283  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[ifield]->GetPhys(),Dy);
284  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[ifield]->GetPhys(),Dz);
285  Exp[nfields ]->FwdTrans(Dx,Exp[nfields ]->UpdateCoeffs());
286  Exp[nfields+1]->FwdTrans(Dy,Exp[nfields+1]->UpdateCoeffs());
287  Exp[nfields+2]->FwdTrans(Dz,Exp[nfields+2]->UpdateCoeffs());
288  }
289  break;
290  default:
291  {
292  ASSERTL0(false,"Expansion dimension not recognised");
293  }
294  break;
295  }
296 
297  //-----------------------------------------------
298  // Write solution to file with additional computed fields
299  string fldfilename(argv[2]);
300  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
301  string endfile("_with_grad.fld");
302  out += endfile;
303  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
304  = Exp[0]->GetFieldDefinitions();
305  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
306 
307  for(j = 0; j < nfields + graddim; ++j)
308  {
309  for(i = 0; i < FieldDef.size(); ++i)
310  {
311  if (j >= nfields)
312  {
313  if(j == nfields)
314  {
315  std::string val = FieldDef[i]->m_fields[ifield] + "_x";
316  FieldDef[i]->m_fields.push_back(val.c_str());
317  }
318  else if(j == nfields+1)
319  {
320  std::string val = FieldDef[i]->m_fields[ifield] + "_y";
321  FieldDef[i]->m_fields.push_back(val.c_str());
322  }
323  else
324  {
325  std::string val = FieldDef[i]->m_fields[ifield] + "_z";
326  FieldDef[i]->m_fields.push_back(val.c_str());
327  }
328  }
329  else
330  {
331  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
332  }
333  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
334  }
335  }
336  LibUtilities::Write(out, FieldDef, FieldData);
337  //-----------------------------------------------
338 
339  return 0;
340 }
341 
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2