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