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