Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FldToVtk.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 using namespace Nektar;
12 
13 #include <sys/stat.h>
14 
15 int fexist( const char *filename ) {
16  struct stat buffer ;
17  if ( stat( filename, &buffer ) ) return 0 ;
18  return 1 ;
19 }
20 
21 int main(int argc, char *argv[])
22 {
23  unsigned int i,j;
24 
25  if(argc < 3)
26  {
27  fprintf(stderr,"Usage: FldToVtk meshfile fieldfile(s)\n");
28  exit(1);
29  }
30 
31 
32  int nExtraPoints;
33  LibUtilities::SessionReaderSharedPtr vSession
34  = LibUtilities::SessionReader::CreateInstance(argc, argv);
35 
36  vSession->LoadParameter("OutputExtraPoints", nExtraPoints,0);
37 
38  //----------------------------------------------
39  // Read in mesh from input file
40  SpatialDomains::MeshGraphSharedPtr graphShPt
41  = SpatialDomains::MeshGraph::Read(vSession);
42  //----------------------------------------------
43 
44  for (int n = 1; n < argc; ++n)
45  {
46  string fname = std::string(argv[n]);
47  if (fname.substr(fname.length() - 1, 1) == "/")
48  {
49  fname = fname.substr(0, fname.length() - 1);
50  }
51  int fdot = fname.find_last_of('.');
52  if (fdot != std::string::npos)
53  {
54  string ending = fname.substr(fdot);
55 
56  // If .chk or .fld we exchange the extension in the output file.
57  // For all other files (e.g. .bse) we append the extension to avoid
58  // conflicts.
59  if (ending == ".chk" || ending == ".fld")
60  {
61  fname = fname.substr(0,fdot);
62  }
63  else if (ending == ".gz")
64  {
65  fname = fname.substr(0,fdot);
66  fdot = fname.find_last_of('.');
67  ASSERTL0(fdot != std::string::npos,
68  "Error: expected file extension before .gz.");
69  ending = fname.substr(fdot);
70  ASSERTL0(ending == ".xml",
71  "Compressed non-xml files are not supported.");
72  continue;
73  }
74  else if (ending == ".xml")
75  {
76  continue;
77  }
78  }
79  fname = fname + ".vtu";
80 
81  if (argc > 3)
82  {
83  if (fexist(fname.c_str()))
84  {
85  cout << "Skipping converted file: " << argv[n] << endl;
86  continue;
87  }
88  cout << "Processing " << argv[n] << endl;
89  }
90 
91  //----------------------------------------------
92  // Import field file.
93  string fieldfile(argv[n]);
94  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
95  vector<vector<NekDouble> > fielddata;
96  LibUtilities::Import(fieldfile,fielddef,fielddata);
97  bool useFFT = false;
98  bool dealiasing = false;
99  //----------------------------------------------
100 
101  //----------------------------------------------
102  // Set up Expansion information
103  for(i = 0; i < fielddef.size(); ++i)
104  {
105  vector<LibUtilities::PointsType> ptype;
106  for(j = 0; j < 3; ++j)
107  {
108  ptype.push_back(LibUtilities::ePolyEvenlySpaced);
109  }
110 
111  fielddef[i]->m_pointsDef = true;
112  fielddef[i]->m_points = ptype;
113 
114  vector<unsigned int> porder;
115  if(fielddef[i]->m_numPointsDef == false)
116  {
117  for(j = 0; j < fielddef[i]->m_numModes.size(); ++j)
118  {
119  porder.push_back(fielddef[i]->m_numModes[j]+nExtraPoints);
120  }
121 
122  fielddef[i]->m_numPointsDef = true;
123  }
124  else
125  {
126  for(j = 0; j < fielddef[i]->m_numPoints.size(); ++j)
127  {
128  porder.push_back(fielddef[i]->m_numPoints[j]+nExtraPoints);
129  }
130  }
131  fielddef[i]->m_numPoints = porder;
132 
133  }
134  graphShPt->SetExpansions(fielddef);
135  //----------------------------------------------
136 
137 
138  //----------------------------------------------
139  // Define Expansion
140  int expdim = graphShPt->GetMeshDimension();
141  int nfields = fielddef[0]->m_fields.size();
142  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
143 
144  switch(expdim)
145  {
146  case 1:
147  {
148  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"NumHomogeneousDir is only set up for 1 or 2");
149 
150  if(fielddef[0]->m_numHomogeneousDir == 1)
151  {
152  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
153 
154  // Define Homogeneous expansion
155  //int nplanes = fielddef[0]->m_numModes[1];
156  int nplanes;
157  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
158 
159  // choose points to be at evenly spaced points at
160  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
161  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
162  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
163 
164  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
165 
166  for(i = 1; i < nfields; ++i)
167  {
168  Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
169  }
170  }
171  else if(fielddef[0]->m_numHomogeneousDir == 2)
172  {
173  MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
174 
175  // Define Homogeneous expansion
176  //int nylines = fielddef[0]->m_numModes[1];
177  //int nzlines = fielddef[0]->m_numModes[2];
178 
179  int nylines;
180  int nzlines;
181  vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
182  vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
183 
184 
185  // choose points to be at evenly spaced points at
186  const LibUtilities::PointsKey PkeyY(nylines+1,LibUtilities::ePolyEvenlySpaced);
187  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
188 
189  const LibUtilities::PointsKey PkeyZ(nzlines+1,LibUtilities::ePolyEvenlySpaced);
190  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
191 
192  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
193  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
194 
195  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
196  Exp[0] = Exp3DH2;
197 
198  for(i = 1; i < nfields; ++i)
199  {
200  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
201  }
202  }
203  else
204  {
205  MultiRegions::ExpList1DSharedPtr Exp1D;
206  Exp1D = MemoryManager<MultiRegions::ExpList1D>
207  ::AllocateSharedPtr(vSession,graphShPt);
208  Exp[0] = Exp1D;
209  for(i = 1; i < nfields; ++i)
210  {
211  Exp[i] = MemoryManager<MultiRegions::ExpList1D>
212  ::AllocateSharedPtr(*Exp1D);
213  }
214  }
215  }
216  break;
217  case 2:
218  {
219  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
220 
221  if(fielddef[0]->m_numHomogeneousDir == 1)
222  {
223  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
224 
225  // Define Homogeneous expansion
226  //int nplanes = fielddef[0]->m_numModes[2];
227 
228  int nplanes;
229  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
230 
231  // choose points to be at evenly spaced points at
232  // nplanes + 1 points
233  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
234  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
235  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
236 
237  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
238  Exp[0] = Exp3DH1;
239 
240  for(i = 1; i < nfields; ++i)
241  {
242  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
243  }
244  }
245  else
246  {
247  MultiRegions::ExpList2DSharedPtr Exp2D;
248  Exp2D = MemoryManager<MultiRegions::ExpList2D>
249  ::AllocateSharedPtr(vSession,graphShPt);
250  Exp[0] = Exp2D;
251 
252  for(i = 1; i < nfields; ++i)
253  {
254  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
255  ::AllocateSharedPtr(*Exp2D);
256  }
257  }
258  }
259  break;
260  case 3:
261  {
262  MultiRegions::ExpList3DSharedPtr Exp3D;
263  Exp3D = MemoryManager<MultiRegions::ExpList3D>
264  ::AllocateSharedPtr(vSession,graphShPt);
265  Exp[0] = Exp3D;
266 
267  for(i = 1; i < nfields; ++i)
268  {
269  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
270  ::AllocateSharedPtr(*Exp3D);
271  }
272  }
273  break;
274  default:
275  ASSERTL0(false,"Expansion dimension not recognised");
276  break;
277  }
278  //----------------------------------------------
279 
280  //----------------------------------------------
281  // Copy data from field file
282  for(j = 0; j < nfields; ++j)
283  {
284  for(int i = 0; i < fielddata.size(); ++i)
285  {
286  Exp[j]->ExtractDataToCoeffs(fielddef [i],
287  fielddata[i],
288  fielddef [i]->m_fields[j],
289  Exp[j]->UpdateCoeffs());
290  }
291  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
292  }
293  //----------------------------------------------
294 
295  //----------------------------------------------
296  // Write solution
297  //string outname(strtok(argv[n],"."));
298  //outname += ".vtu";
299  ofstream outfile(fname.c_str());
300  Exp[0]->WriteVtkHeader(outfile);
301 
302  // For each field write out field data for each expansion.
303  for(i = 0; i < Exp[0]->GetNumElmts(); ++i)
304  {
305  Exp[0]->WriteVtkPieceHeader(outfile,i);
306  // For this expansion, write out each field.
307  for(j = 0; j < Exp.num_elements(); ++j)
308  {
309  Exp[j]->WriteVtkPieceData(outfile,i, fielddef[0]->m_fields[j]);
310  }
311  Exp[0]->WriteVtkPieceFooter(outfile,i);
312  }
313  Exp[0]->WriteVtkFooter(outfile);
314  cout << "Written file: " << fname << endl;
315  //----------------------------------------------
316  }
317  return 0;
318 }
319 
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2