Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
CalcGrad.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
Include dependency graph for CalcGrad.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 20 of file CalcGrad.cpp.

References ASSERTL0, Nektar::MultiRegions::DirCartesianMap, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::Import(), and Nektar::LibUtilities::Write().

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