Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
OutputTecplot.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: OutputTecplot.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Dat file format output.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <set>
37 #include <string>
38 #include <iomanip>
39 using namespace std;
40 
41 #include "OutputTecplot.h"
42 
43 namespace Nektar
44 {
45  namespace Utilities
46  {
47 
48  ModuleKey OutputTecplot::m_className =
49  GetModuleFactory().RegisterCreatorFunction(
50  ModuleKey(eOutputModule, "dat"), OutputTecplot::create,
51  "Writes a Tecplot file.");
52 
53  OutputTecplot::OutputTecplot(FieldSharedPtr f) : OutputModule(f)
54  {
55  m_requireEquiSpaced = true;
56  if(f->m_setUpEquiSpacedFields)
57  {
58  m_outputType = eFullBlockZoneEquiSpaced;
59  }
60  else
61  {
62  m_outputType = eFullBlockZone;
63  }
64  }
65 
66  OutputTecplot::~OutputTecplot()
67  {
68  }
69 
70  void OutputTecplot::Process(po::variables_map &vm)
71  {
72  FieldPtsSharedPtr f = m_f->m_fieldPts;
73 
74  m_doError = (vm.count("error") == 1)? true: false;
75 
76  if (m_f->m_verbose)
77  {
78  cout << "OutputTecplot: Writing file..." << endl;
79  }
80  // Do nothing if no expansion defined
81  if (f == NullFieldPts && !m_f->m_exp.size())
82  {
83  return;
84  }
85 
86  // Extract the output filename and extension
87  string filename = m_config["outfile"].as<string>();
88 
89  if(f != NullFieldPts)
90  {
91  int i = 0;
92  int j = 0;
93  int dim = f->m_ptsDim;
94 
95  if(f->m_pts[0].num_elements() == 0)
96  {
97  return;
98  }
99 
100  // Write solution.
101  ofstream outfile(filename.c_str());
102 
103  switch(dim)
104  {
105  case 1:
106  outfile << "VARIABLES = x";
107  break;
108  case 2:
109  outfile << "VARIABLES = x,y";
110  break;
111  case 3:
112  outfile << "VARIABLES = x,y,z";
113  break;
114  }
115 
116  for(i = 0; i < f->m_fields.size(); ++i)
117  {
118  outfile << "," << f->m_fields[i];
119  }
120  outfile << endl;
121  bool DumpAsFEPoint = true;
122  switch(f->m_ptype)
123  {
124  case Utilities::ePtsFile:
125  case Utilities::ePtsLine:
126  outfile << " ZONE I="
127  << f->m_pts[0].num_elements()
128  << " F=POINT" << endl;
129  break;
130  case Utilities::ePtsPlane:
131  outfile << " ZONE I=" << f->m_npts[0]
132  << " J=" << f->m_npts[1]
133  << " F=POINT" << endl;
134  break;
135  case Utilities::ePtsTriBlock:
136  {
137  int numBlocks = 0;
138  for(i = 0; i < f->m_ptsConn.size(); ++i)
139  {
140  numBlocks +=
141  f->m_ptsConn[i].num_elements()/3;
142  }
143  outfile << "Zone, N="
144  << f->m_pts[0].num_elements()
145  << ", E=" << numBlocks
146  << ", F=FEBlock" << ", ET=TRIANGLE"
147  << std::endl;
148  DumpAsFEPoint = false;
149  break;
150  }
151  case Utilities::ePtsTetBlock:
152  {
153  int numBlocks = 0;
154  for(i = 0; i < f->m_ptsConn.size(); ++i)
155  {
156  numBlocks +=
157  f->m_ptsConn[i].num_elements()/4;
158  }
159  outfile << "Zone, N="
160  << f->m_pts[0].num_elements()
161  << ", E=" << numBlocks
162  << ", F=FEBlock" << ", ET=TETRAHEDRON"
163  << std::endl;
164  DumpAsFEPoint = false;
165  break;
166  }
167  }
168 
169  if(DumpAsFEPoint) // dump in point format
170  {
171  for(i = 0; i < f->m_pts[0].num_elements(); ++i)
172  {
173  for(j = 0; j < dim; ++j)
174  {
175  outfile << std::setw(12)
176  << f->m_pts[j][i] << " ";
177  }
178 
179  for(j = 0; j < f->m_fields.size(); ++j)
180  {
181  outfile << std::setw(12)
182  << m_f->m_data[j][i] << " ";
183  }
184  outfile << endl;
185  }
186  }
187  else // dump in block format
188  {
189  for(j = 0; j < dim + f->m_fields.size(); ++j)
190  {
191  for(i = 0; i < f->m_pts[0].num_elements(); ++i)
192  {
193  outfile << f->m_pts[j][i] << " ";
194  if((!(i % 1000))&&i)
195  {
196  outfile << std::endl;
197  }
198  }
199  outfile << endl;
200  }
201 
202 
203  // dump connectivity data if it exists
204  for(i = 0; i < f->m_ptsConn.size();++i)
205  {
206  for(j = 0; j < f->m_ptsConn[i].num_elements(); ++j)
207  {
208  outfile << f->m_ptsConn[i][j] +1 << " ";
209  if( ( !(j % 10 * dim) ) && j )
210  {
211  outfile << std::endl;
212  }
213  }
214  }
215  }
216  }
217  else
218  {
219  // Amend for parallel output if required
220  if(m_f->m_session->GetComm()->GetSize() != 1)
221  {
222  int dot = filename.find_last_of('.');
223  string ext = filename.substr(dot,filename.length()-dot);
224  string procId = "_P" + boost::lexical_cast<std::string>(
225  m_f->m_session->GetComm()->GetRank());
226  string start = filename.substr(0,dot);
227  filename = start + procId + ext;
228  }
229 
230  // Write solution.
231  ofstream outfile(filename.c_str());
232  std::string var;
233  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fDef =
234  m_f->m_fielddef;
235  if (fDef.size())
236  {
237  var = fDef[0]->m_fields[0];
238 
239  for (int j = 1; j < fDef[0]->m_fields.size(); ++j)
240  {
241  var = var + ", " + fDef[0]->m_fields[j];
242  }
243  }
244 
245  WriteTecplotHeader(outfile,var);
246  WriteTecplotZone(outfile);
247  if(var.length()) // see if any variables are defined
248  {
249  for(int j = 0; j < m_f->m_exp.size(); ++j)
250  {
251  WriteTecplotField(j,outfile);
252  }
253  }
254 
255  WriteTecplotConnectivity(outfile);
256  }
257 
258  cout << "Written file: " << filename << endl;
259  }
260 
261  /**
262  * Write Tecplot Files Header
263  * @param outfile Output file name.
264  * @param var variables names
265  */
266  void OutputTecplot::WriteTecplotHeader(std::ofstream &outfile,
267  std::string var)
268  {
269  int coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
270  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
271 
272  if(HomoExpType == MultiRegions::e3DH1D)
273  {
274  coordim +=1;
275  }
276  else if (HomoExpType == MultiRegions::e3DH2D)
277  {
278  coordim += 2;
279  }
280 
281  outfile << "Variables = x";
282 
283  if(coordim == 2)
284  {
285  outfile << ", y";
286  }
287  else if (coordim == 3)
288  {
289  outfile << ", y, z";
290  }
291 
292  if(var.length())
293  {
294  outfile << ", "<< var << std::endl << std::endl;
295  }
296  else
297  {
298  outfile << std::endl << std::endl;
299  }
300  }
301 
302 
303  /**
304  * Write Tecplot Files Zone
305  * @param outfile Output file name.
306  * @param expansion Expansion that is considered
307  */
308  void OutputTecplot::WriteTecplotZone(std::ofstream &outfile)
309  {
310  switch(m_outputType)
311  {
312  case eFullBlockZone: //write as full block zone
313  {
314  int i,j;
315  int coordim = m_f->m_exp[0]->GetCoordim(0);
316  int totpoints = m_f->m_exp[0]->GetTotPoints();
317  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
318 
319  Array<OneD,NekDouble> coords[3];
320 
321  coords[0] = Array<OneD,NekDouble>(totpoints);
322  coords[1] = Array<OneD,NekDouble>(totpoints);
323  coords[2] = Array<OneD,NekDouble>(totpoints);
324 
325  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
326 
327  if (m_doError)
328  {
329  NekDouble l2err;
330  std::string coordval[] = {"x","y","z"};
331  int rank = m_f->m_session->GetComm()->GetRank();
332 
333  for(int i = 0; i < coordim; ++i)
334  {
335  l2err = m_f->m_exp[0]->L2(coords[i]);
336  if(rank == 0)
337  {
338  cout << "L 2 error (variable "
339  << coordval[i] << ") : " << l2err << endl;
340  }
341  }
342  }
343 
344  int numBlocks = GetNumTecplotBlocks();
345  int nBases = m_f->m_exp[0]->GetExp(0)->GetNumBases();
346 
347  if (HomoExpType == MultiRegions::e3DH1D)
348  {
349  nBases += 1;
350  coordim += 1;
351  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
352  NekDouble tmp = numBlocks * (nPlanes-1);
353  numBlocks = (int)tmp;
354  }
355  else if (HomoExpType == MultiRegions::e3DH2D)
356  {
357  nBases += 2;
358  coordim += 1;
359  }
360 
361  outfile << "Zone, N=" << totpoints << ", E="<<
362  numBlocks << ", F=FEBlock" ;
363 
364  switch(nBases)
365  {
366  case 1:
367  outfile << ", ET=LINESEG" << std::endl;
368  break;
369  case 2:
370  outfile << ", ET=QUADRILATERAL" << std::endl;
371  break;
372  case 3:
373  outfile << ", ET=BRICK" << std::endl;
374  break;
375  }
376 
377  // write out coordinates in block format
378  for(j = 0; j < coordim; ++j)
379  {
380  for(i = 0; i < totpoints; ++i)
381  {
382  outfile << coords[j][i] << " ";
383  if((!(i % 1000))&&i)
384  {
385  outfile << std::endl;
386  }
387  }
388  outfile << std::endl;
389  }
390  break;
391  }
392  case eSeperateZones:
393  {
394  for(int i = 0; i < m_f->m_exp[0]->GetExpSize(); ++i)
395  {
396  m_f->m_exp[0]->WriteTecplotZone(outfile,i);
397  }
398  break;
399  }
400  case eFullBlockZoneEquiSpaced:
401  ASSERTL0(false,
402  "Should not have this option in this method");
403  break;
404  }
405  }
406 
407  int OutputTecplot::GetNumTecplotBlocks(void)
408  {
409  int returnval = 0;
410 
411  if(m_f->m_exp[0]->GetExp(0)->GetNumBases() == 1)
412  {
413  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
414  {
415  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1);
416  }
417  }
418  else if(m_f->m_exp[0]->GetExp(0)->GetNumBases() == 2)
419  {
420  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
421  {
422  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1)*
423  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1)-1);
424  }
425  }
426  else
427  {
428  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
429  {
430  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1)*
431  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1)-1)*
432  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(2)-1);
433  }
434  }
435 
436  return returnval;
437  }
438 
439 
440  /**
441  * Write Tecplot Files Field
442  * @param outfile Output file name.
443  * @param expansion Expansion that is considered
444  */
445  void OutputTecplot::WriteTecplotField(const int field,
446  std::ofstream &outfile)
447  {
448 
449  if(m_outputType == eFullBlockZone) //write as full block zone
450  {
451  int totpoints = m_f->m_exp[0]->GetTotPoints();
452 
453  if(m_f->m_exp[field]->GetPhysState() == false)
454  {
455  m_f->m_exp[field]->BwdTrans(m_f->m_exp[field]->GetCoeffs(),
456  m_f->m_exp[field]->UpdatePhys());
457  }
458 
459  if (m_doError)
460  {
461  NekDouble l2err = m_f->m_exp[0]->L2(m_f->m_exp[field]->UpdatePhys());
462 
463  if(m_f->m_session->GetComm()->GetRank() == 0)
464  {
465  cout << "L 2 error (variable "
466  << m_f->m_fielddef[0]->m_fields[field] << ") : "
467  << l2err << endl;
468  }
469  }
470  else
471  {
472  for(int i = 0; i < totpoints; ++i)
473  {
474  outfile << m_f->m_exp[field]->GetPhys()[i] << " ";
475  if((!(i % 1000))&&i)
476  {
477  outfile << std::endl;
478  }
479  }
480  }
481  outfile << std::endl;
482  }
483  else
484  {
485  for(int e = 0; e < m_f->m_exp[field]->GetExpSize(); ++e)
486  {
487  m_f->m_exp[field]->WriteTecplotField(outfile,e);
488  }
489  }
490  }
491 
492  void OutputTecplot::WriteTecplotConnectivity(std::ofstream &outfile)
493  {
494  int i,j,k,l;
495  int nbase = m_f->m_exp[0]->GetExp(0)->GetNumBases();
496  int cnt = 0;
497 
498  for(i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
499  {
500  if(nbase == 1)
501  {
502  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
503 
504  for(k = 1; k < np0; ++k)
505  {
506  outfile << cnt + k +1 << " ";
507  outfile << cnt + k << endl;
508  }
509 
510  cnt += np0;
511  }
512  else if(nbase == 2)
513  {
514  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
515  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
516  int totPoints = m_f->m_exp[0]->GetTotPoints();
517  int nPlanes = 1;
518 
519  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
520  {
521  nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
522  totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
523 
524 
525  for(int n = 1; n < nPlanes; ++n)
526  {
527  for(j = 1; j < np1; ++j)
528  {
529  for(k = 1; k < np0; ++k)
530  {
531  outfile << cnt + (n-1)*totPoints + (j-1)*np0 + k
532  << " ";
533  outfile << cnt + (n-1)*totPoints + (j-1)*np0 + k + 1
534  << " ";
535  outfile << cnt + (n-1)*totPoints + j*np0 + k + 1
536  << " ";
537  outfile << cnt + (n-1)*totPoints + j*np0 + k
538  << " ";
539 
540  outfile << cnt + n*totPoints + (j-1)*np0 + k
541  << " ";
542  outfile << cnt + n*totPoints + (j-1)*np0 + k + 1
543  << " ";
544  outfile << cnt + n*totPoints + j*np0 + k + 1
545  << " ";
546  outfile << cnt + n*totPoints + j*np0 + k << endl;
547  }
548  }
549  }
550  cnt += np0*np1;
551  }
552  else
553  {
554  for(j = 1; j < np1; ++j)
555  {
556  for(k = 1; k < np0; ++k)
557  {
558  outfile << cnt + (j-1)*np0 + k << " ";
559  outfile << cnt + (j-1)*np0 + k +1 << " ";
560  outfile << cnt + j*np0 + k +1 << " ";
561  outfile << cnt + j*np0 + k << endl;
562  }
563  }
564  cnt += np0*np1;
565  }
566  }
567  else if(nbase == 3)
568  {
569  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
570  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
571  int np2 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(2);
572 
573  for(j = 1; j < np2; ++j)
574  {
575  for(k = 1; k < np1; ++k)
576  {
577  for(l = 1; l < np0; ++l)
578  {
579  outfile << cnt + (j-1)*np0*np1 + (k-1)*np0 + l << " ";
580  outfile << cnt + (j-1)*np0*np1 + (k-1)*np0 + l +1 << " ";
581  outfile << cnt + (j-1)*np0*np1 + k*np0 + l +1 << " ";
582  outfile << cnt + (j-1)*np0*np1 + k*np0 + l << " ";
583 
584  outfile << cnt + j*np0*np1 + (k-1)*np0 + l << " ";
585  outfile << cnt + j*np0*np1 + (k-1)*np0 + l +1 << " ";
586  outfile << cnt + j*np0*np1 + k*np0 + l +1 << " ";
587  outfile << cnt + j*np0*np1 + k*np0 + l << endl;
588  }
589  }
590  }
591  cnt += np0*np1*np2;
592  }
593  else
594  {
595  ASSERTL0(false,"Not set up for this dimension");
596  }
597 
598  }
599  }
600  }
601 }
602 
603 
604 
605 
Generated on Wed Mar 4 2015 16:33:16 for Nektar++ by   doxygen 1.8.1.2