Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
CalcDivergence.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 CalcDivergence.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 21 of file CalcDivergence.cpp.

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

22 {
23  int i,j;
24 
25  if(argc != 3)
26  {
27  fprintf(stderr,"Usage: ./CalcDivergence file.xml file.fld\n");
28  exit(1);
29  }
30 
31  LibUtilities::SessionReaderSharedPtr vSession
32  = LibUtilities::SessionReader::CreateInstance(argc, argv);
33 
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 
56  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + 1);
57 
58  switch(expdim)
59  {
60  case 1:
61  {
62  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
63 
64  if(fielddef[0]->m_numHomogeneousDir == 1)
65  {
66  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
67 
68  // Define Homogeneous expansion
69  //int nplanes = fielddef[0]->m_numModes[1];
70 
71  int nplanes;
72  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
73 
74  // choose points to be at evenly spaced points at
75  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
76  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
77  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
78 
79  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
80  Exp[0] = Exp2DH1;
81 
82  for(i = 1; i < nfields + 1; ++i)
83  {
84  Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
85  }
86  }
87  else if(fielddef[0]->m_numHomogeneousDir == 2)
88  {
89  MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
90 
91  // Define Homogeneous expansion
92  //int nylines = fielddef[0]->m_numModes[1];
93  //int nzlines = fielddef[0]->m_numModes[2];
94 
95  int nylines;
96  int nzlines;
97  vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
98  vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
99 
100  // choose points to be at evenly spaced points at
101  const LibUtilities::PointsKey PkeyY(nylines+1,LibUtilities::ePolyEvenlySpaced);
102  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
103 
104  const LibUtilities::PointsKey PkeyZ(nzlines+1,LibUtilities::ePolyEvenlySpaced);
105  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
106 
107  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
108  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
109 
110  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
111  Exp[0] = Exp3DH2;
112 
113  for(i = 1; i < nfields + 1; ++i)
114  {
115  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
116  }
117  }
118  else
119  {
120  MultiRegions::ExpList1DSharedPtr Exp1D;
121  Exp1D = MemoryManager<MultiRegions::ExpList1D>
122  ::AllocateSharedPtr(vSession,graphShPt);
123  Exp[0] = Exp1D;
124  for(i = 1; i < nfields + 1; ++i)
125  {
126  Exp[i] = MemoryManager<MultiRegions::ExpList1D>
127  ::AllocateSharedPtr(*Exp1D);
128  }
129  }
130  }
131  break;
132  case 2:
133  {
134  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
135 
136  if(fielddef[0]->m_numHomogeneousDir == 1)
137  {
138  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
139 
140  // Define Homogeneous expansion
141  //int nplanes = fielddef[0]->m_numModes[2];
142 
143  int nplanes;
144  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
145 
146  // choose points to be at evenly spaced points at
147  // nplanes + 1 points
148  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
149  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
150  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
151 
152  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
153  Exp[0] = Exp3DH1;
154 
155  for(i = 1; i < nfields + 1; ++i)
156  {
157  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
158  }
159  }
160  else
161  {
162  MultiRegions::ExpList2DSharedPtr Exp2D;
163  Exp2D = MemoryManager<MultiRegions::ExpList2D>
164  ::AllocateSharedPtr(vSession,graphShPt);
165  Exp[0] = Exp2D;
166 
167  for(i = 1; i < nfields + 1; ++i)
168  {
169  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
170  ::AllocateSharedPtr(*Exp2D);
171  }
172  }
173  }
174  break;
175  case 3:
176  {
177  MultiRegions::ExpList3DSharedPtr Exp3D;
178  Exp3D = MemoryManager<MultiRegions::ExpList3D>
179  ::AllocateSharedPtr(vSession,graphShPt);
180  Exp[0] = Exp3D;
181 
182  for(i = 1; i < nfields + 1; ++i)
183  {
184  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
185  ::AllocateSharedPtr(*Exp3D);
186  }
187  }
188  break;
189  default:
190  ASSERTL0(false,"Expansion dimension not recognised");
191  break;
192  }
193 
194  //----------------------------------------------
195  // Copy data from field file
196  for(j = 0; j < nfields; ++j)
197  {
198  for(unsigned int i = 0; i < fielddata.size(); ++i)
199  {
200  Exp[j]->ExtractDataToCoeffs(fielddef [i],
201  fielddata[i],
202  fielddef [i]->m_fields[j],
203  Exp[j]->UpdateCoeffs());
204  }
205  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
206  }
207  //----------------------------------------------
208 
209  int nq = Exp[0]->GetNpoints();
210 
211  Array<OneD, NekDouble> Ux(nq);
212  Array<OneD, NekDouble> Vy(nq);
213  Array<OneD, NekDouble> Div(nq);
214 
215  switch(expdim)
216  {
217  case 1:
218  {
219  if(fielddef[0]->m_numHomogeneousDir == 1)
220  {
221  ASSERTL0(false,"Not implemented yet");
222  }
223  else if(fielddef[0]->m_numHomogeneousDir == 2)
224  {
225  ASSERTL0(false,"Not implemented yet");
226  }
227  else
228  {
229  ASSERTL0(false,"Not implemented yet");
230  }
231  }
232  break;
233  case 2:
234  {
235  if(fielddef[0]->m_numHomogeneousDir == 1)
236  {
237  ASSERTL0(false,"Not implemented yet");
238  }
239  else
240  {
241  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[0]->GetPhys(),Ux);
242  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[1]->GetPhys(),Vy);
243 
244  Vmath::Vsub(nq,Ux,1,Vy,1,Div,1);
245 
246  Exp[3]->FwdTrans(Div,Exp[3]->UpdateCoeffs());
247  }
248  }
249  break;
250  case 3:
251  {
252  ASSERTL0(false,"Not implemented yet");
253  }
254  break;
255  default:
256  {
257  ASSERTL0(false,"Expansion dimension not recognised");
258  }
259  break;
260  }
261 
262  //-----------------------------------------------
263  // Write solution to file with additional computed fields
264  string fldfilename(argv[2]);
265  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
266  string endfile("_with_divergence.fld");
267  out += endfile;
268  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
269  = Exp[0]->GetFieldDefinitions();
270  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
271 
272  for(j = 0; j < nfields + 1; ++j)
273  {
274  for(i = 0; i < FieldDef.size(); ++i)
275  {
276  if (j >= nfields)
277  {
278  FieldDef[i]->m_fields.push_back("Div");
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 }
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
Vmath::Vsub
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
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 Jul 16 2016 10:00:54 for Nektar++ by   doxygen 1.8.8