Nektar++
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 
43 
44 #include "OutputTecplot.h"
45 
46 namespace Nektar
47 {
48 namespace Utilities
49 {
50 
51 ModuleKey OutputTecplot::m_className =
53  ModuleKey(eOutputModule, "dat"), OutputTecplot::create,
54  "Writes a Tecplot file.");
55 
56 OutputTecplot::OutputTecplot(FieldSharedPtr f) : OutputModule(f)
57 {
58  if(f->m_setUpEquiSpacedFields)
59  {
61  }
62  else
63  {
64  m_requireEquiSpaced = true;
66  }
67 }
68 
70 {
71 }
72 
73 void OutputTecplot::Process(po::variables_map &vm)
74 {
75  LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
76 
77  m_doError = (vm.count("error") == 1)? true: false;
78 
79  if (m_f->m_verbose)
80  {
81  cout << "OutputTecplot: Writing file..." << endl;
82  }
83  // Do nothing if no expansion defined
84  if (fPts == LibUtilities::NullPtsField && !m_f->m_exp.size())
85  {
86  return;
87  }
88 
89 
90  // Extract the output filename and extension
91  string filename = m_config["outfile"].as<string>();
92 
93  // Amend for parallel output if required
94  if(m_f->m_session->GetComm()->GetSize() != 1)
95  {
96  int dot = filename.find_last_of('.');
97  string ext = filename.substr(dot,filename.length()-dot);
98  string procId = "_P" + boost::lexical_cast<std::string>(
99  m_f->m_session->GetComm()->GetRank());
100  string start = filename.substr(0,dot);
101  filename = start + procId + ext;
102  }
103 
104  if(fPts != LibUtilities::NullPtsField)
105  {
106  int i = 0;
107  int j = 0;
108  int dim = fPts->GetDim();
109 
110  if(fPts->GetNpoints() == 0)
111  {
112  return;
113  }
114 
115  // Write solution.
116  ofstream outfile(filename.c_str());
117 
118  switch(dim)
119  {
120  case 1:
121  outfile << "VARIABLES = x";
122  break;
123  case 2:
124  outfile << "VARIABLES = x,y";
125  break;
126  case 3:
127  outfile << "VARIABLES = x,y,z";
128  break;
129  }
130 
131  vector<Array<OneD, int> > ptsConn;
132  fPts->GetConnectivity(ptsConn);
133 
134  for(i = 0; i < fPts->GetNFields(); ++i)
135  {
136  outfile << "," << fPts->GetFieldName(i);
137  }
138  outfile << endl;
139  bool DumpAsFEPoint = true;
140  switch(fPts->GetPtsType())
141  {
144  {
145  outfile << " ZONE I="
146  << fPts->GetNpoints()
147  << " F=POINT" << endl;
148  break;
149  }
151  {
152  outfile << " ZONE I=" << fPts->GetPointsPerEdge(0)
153  << " J=" << fPts->GetPointsPerEdge(1)
154  << " F=POINT" << endl;
155  break;
156  }
158  {
159  int numBlocks = 0;
160  for(i = 0; i < ptsConn.size(); ++i)
161  {
162  numBlocks +=
163  ptsConn[i].num_elements()/3;
164  }
165  outfile << "Zone, N="
166  << fPts->GetNpoints()
167  << ", E=" << numBlocks
168  << ", F=FEBlock" << ", ET=TRIANGLE"
169  << std::endl;
170  DumpAsFEPoint = false;
171  break;
172  }
174  {
175  int numBlocks = 0;
176  for(i = 0; i < ptsConn.size(); ++i)
177  {
178  numBlocks +=
179  ptsConn[i].num_elements()/4;
180  }
181  outfile << "Zone, N="
182  << fPts->GetNpoints()
183  << ", E=" << numBlocks
184  << ", F=FEBlock" << ", ET=TETRAHEDRON"
185  << std::endl;
186  DumpAsFEPoint = false;
187  break;
188  }
189  default:
190  ASSERTL0(false, "ptsType not supported yet.");
191  }
192 
193  if(DumpAsFEPoint) // dump in point format
194  {
195  for(i = 0; i < fPts->GetNpoints(); ++i)
196  {
197  for(j = 0; j < dim; ++j)
198  {
199  outfile << std::setw(12)
200  << fPts->GetPointVal(j, i) << " ";
201  }
202 
203  for(j = 0; j < fPts->GetNFields(); ++j)
204  {
205  outfile << std::setw(12)
206  << m_f->m_data[j][i] << " ";
207  }
208  outfile << endl;
209  }
210  }
211  else // dump in block format
212  {
213  for(j = 0; j < dim + fPts->GetNFields(); ++j)
214  {
215  for(i = 0; i < fPts->GetNpoints(); ++i)
216  {
217  outfile << fPts->GetPointVal(j, i) << " ";
218  if((!(i % 1000))&&i)
219  {
220  outfile << std::endl;
221  }
222  }
223  outfile << endl;
224  }
225 
226  // dump connectivity data if it exists
227  for(i = 0; i < ptsConn.size();++i)
228  {
229  for(j = 0; j < ptsConn[i].num_elements(); ++j)
230  {
231  outfile << ptsConn[i][j] +1 << " ";
232  if( ( !(j % 10 * dim) ) && j )
233  {
234  outfile << std::endl;
235  }
236  }
237  }
238  }
239  }
240  else
241  {
242 
243  // Write solution.
244  ofstream outfile(filename.c_str());
245  std::string var;
246  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fDef =
247  m_f->m_fielddef;
248  if (fDef.size())
249  {
250  var = fDef[0]->m_fields[0];
251 
252  for (int j = 1; j < fDef[0]->m_fields.size(); ++j)
253  {
254  var = var + ", " + fDef[0]->m_fields[j];
255  }
256  }
257 
258  WriteTecplotHeader(outfile,var);
259  WriteTecplotZone(outfile);
260  if(var.length()) // see if any variables are defined
261  {
262  for(int j = 0; j < m_f->m_exp.size(); ++j)
263  {
264  WriteTecplotField(j,outfile);
265  }
266  }
267 
268  WriteTecplotConnectivity(outfile);
269  }
270 
271  cout << "Written file: " << filename << endl;
272 }
273 
274 /**
275  * Write Tecplot Files Header
276  * @param outfile Output file name.
277  * @param var variables names
278  */
279 void OutputTecplot::WriteTecplotHeader(std::ofstream &outfile,
280  std::string var)
281 {
282  int coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
283  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
284 
285  if(HomoExpType == MultiRegions::e3DH1D)
286  {
287  coordim +=1;
288  }
289  else if (HomoExpType == MultiRegions::e3DH2D)
290  {
291  coordim += 2;
292  }
293 
294  outfile << "Variables = x";
295 
296  if(coordim == 2)
297  {
298  outfile << ", y";
299  }
300  else if (coordim == 3)
301  {
302  outfile << ", y, z";
303  }
304 
305  if(var.length())
306  {
307  outfile << ", "<< var << std::endl << std::endl;
308  }
309  else
310  {
311  outfile << std::endl << std::endl;
312  }
313 }
314 
315 
316 /**
317  * Write Tecplot Files Zone
318  * @param outfile Output file name.
319  * @param expansion Expansion that is considered
320  */
321 void OutputTecplot::WriteTecplotZone(std::ofstream &outfile)
322 {
323  switch(m_outputType)
324  {
325  case eFullBlockZone: //write as full block zone
326  {
327  int i,j;
328  int coordim = m_f->m_exp[0]->GetCoordim(0);
329  int totpoints = m_f->m_exp[0]->GetTotPoints();
330  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
331 
332  Array<OneD,NekDouble> coords[3];
333 
334  coords[0] = Array<OneD,NekDouble>(totpoints);
335  coords[1] = Array<OneD,NekDouble>(totpoints);
336  coords[2] = Array<OneD,NekDouble>(totpoints);
337 
338  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
339 
340  if (m_doError)
341  {
342  NekDouble l2err;
343  std::string coordval[] = {"x","y","z"};
344  int rank = m_f->m_session->GetComm()->GetRank();
345 
346  for(int i = 0; i < coordim; ++i)
347  {
348  l2err = m_f->m_exp[0]->L2(coords[i]);
349  if(rank == 0)
350  {
351  cout << "L 2 error (variable "
352  << coordval[i] << ") : " << l2err << endl;
353  }
354  }
355  }
356 
357  int numBlocks = GetNumTecplotBlocks();
358  int nBases = m_f->m_exp[0]->GetExp(0)->GetNumBases();
359 
360  if (HomoExpType == MultiRegions::e3DH1D)
361  {
362  nBases += 1;
363  coordim += 1;
364  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
365  NekDouble tmp = numBlocks * (nPlanes-1);
366  numBlocks = (int)tmp;
367  }
368  else if (HomoExpType == MultiRegions::e3DH2D)
369  {
370  nBases += 2;
371  coordim += 1;
372  }
373 
374  outfile << "Zone, N=" << totpoints << ", E="<<
375  numBlocks << ", F=FEBlock" ;
376 
377  switch(nBases)
378  {
379  case 1:
380  outfile << ", ET=LINESEG" << std::endl;
381  break;
382  case 2:
383  outfile << ", ET=QUADRILATERAL" << std::endl;
384  break;
385  case 3:
386  outfile << ", ET=BRICK" << std::endl;
387  break;
388  }
389 
390  // write out coordinates in block format
391  for(j = 0; j < coordim; ++j)
392  {
393  for(i = 0; i < totpoints; ++i)
394  {
395  outfile << coords[j][i] << " ";
396  if((!(i % 1000))&&i)
397  {
398  outfile << std::endl;
399  }
400  }
401  outfile << std::endl;
402  }
403  break;
404  }
405  case eSeperateZones:
406  {
407  for(int i = 0; i < m_f->m_exp[0]->GetExpSize(); ++i)
408  {
409  m_f->m_exp[0]->WriteTecplotZone(outfile,i);
410  }
411  break;
412  }
414  ASSERTL0(false,
415  "Should not have this option in this method");
416  break;
417  }
418 }
419 
421 {
422  int returnval = 0;
423 
424  if(m_f->m_exp[0]->GetExp(0)->GetNumBases() == 1)
425  {
426  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
427  {
428  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1);
429  }
430  }
431  else if(m_f->m_exp[0]->GetExp(0)->GetNumBases() == 2)
432  {
433  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
434  {
435  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1)*
436  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1)-1);
437  }
438  }
439  else
440  {
441  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
442  {
443  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1)*
444  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1)-1)*
445  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(2)-1);
446  }
447  }
448 
449  return returnval;
450 }
451 
452 
453 /**
454  * Write Tecplot Files Field
455  * @param outfile Output file name.
456  * @param expansion Expansion that is considered
457  */
458 void OutputTecplot::WriteTecplotField(const int field,
459  std::ofstream &outfile)
460 {
461 
462  if(m_outputType == eFullBlockZone) //write as full block zone
463  {
464  int totpoints = m_f->m_exp[0]->GetTotPoints();
465 
466  if(m_f->m_exp[field]->GetPhysState() == false)
467  {
468  m_f->m_exp[field]->BwdTrans(m_f->m_exp[field]->GetCoeffs(),
469  m_f->m_exp[field]->UpdatePhys());
470  }
471 
472  if (m_doError)
473  {
474  NekDouble l2err = m_f->m_exp[0]->L2(m_f->m_exp[field]->UpdatePhys());
475 
476  if(m_f->m_session->GetComm()->GetRank() == 0)
477  {
478  cout << "L 2 error (variable "
479  << m_f->m_fielddef[0]->m_fields[field] << ") : "
480  << l2err << endl;
481  }
482  }
483  else
484  {
485  for(int i = 0; i < totpoints; ++i)
486  {
487  outfile << m_f->m_exp[field]->GetPhys()[i] << " ";
488  if((!(i % 1000))&&i)
489  {
490  outfile << std::endl;
491  }
492  }
493  }
494  outfile << std::endl;
495  }
496  else
497  {
498  for(int e = 0; e < m_f->m_exp[field]->GetExpSize(); ++e)
499  {
500  m_f->m_exp[field]->WriteTecplotField(outfile,e);
501  }
502  }
503 }
504 
505 void OutputTecplot::WriteTecplotConnectivity(std::ofstream &outfile)
506 {
507  int i,j,k,l;
508  int nbase = m_f->m_exp[0]->GetExp(0)->GetNumBases();
509  int cnt = 0;
510 
511  for(i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
512  {
513  if(nbase == 1)
514  {
515  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
516 
517  for(k = 1; k < np0; ++k)
518  {
519  outfile << cnt + k +1 << " ";
520  outfile << cnt + k << endl;
521  }
522 
523  cnt += np0;
524  }
525  else if(nbase == 2)
526  {
527  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
528  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
529  int totPoints = m_f->m_exp[0]->GetTotPoints();
530  int nPlanes = 1;
531 
532  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
533  {
534  nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
535  totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
536 
537 
538  for(int n = 1; n < nPlanes; ++n)
539  {
540  for(j = 1; j < np1; ++j)
541  {
542  for(k = 1; k < np0; ++k)
543  {
544  outfile << cnt + (n-1)*totPoints + (j-1)*np0 + k
545  << " ";
546  outfile << cnt + (n-1)*totPoints + (j-1)*np0 + k + 1
547  << " ";
548  outfile << cnt + (n-1)*totPoints + j*np0 + k + 1
549  << " ";
550  outfile << cnt + (n-1)*totPoints + j*np0 + k
551  << " ";
552 
553  outfile << cnt + n*totPoints + (j-1)*np0 + k
554  << " ";
555  outfile << cnt + n*totPoints + (j-1)*np0 + k + 1
556  << " ";
557  outfile << cnt + n*totPoints + j*np0 + k + 1
558  << " ";
559  outfile << cnt + n*totPoints + j*np0 + k << endl;
560  }
561  }
562  }
563  cnt += np0*np1;
564  }
565  else
566  {
567  for(j = 1; j < np1; ++j)
568  {
569  for(k = 1; k < np0; ++k)
570  {
571  outfile << cnt + (j-1)*np0 + k << " ";
572  outfile << cnt + (j-1)*np0 + k +1 << " ";
573  outfile << cnt + j*np0 + k +1 << " ";
574  outfile << cnt + j*np0 + k << endl;
575  }
576  }
577  cnt += np0*np1;
578  }
579  }
580  else if(nbase == 3)
581  {
582  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
583  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
584  int np2 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(2);
585 
586  for(j = 1; j < np2; ++j)
587  {
588  for(k = 1; k < np1; ++k)
589  {
590  for(l = 1; l < np0; ++l)
591  {
592  outfile << cnt + (j-1)*np0*np1 + (k-1)*np0 + l << " ";
593  outfile << cnt + (j-1)*np0*np1 + (k-1)*np0 + l +1 << " ";
594  outfile << cnt + (j-1)*np0*np1 + k*np0 + l +1 << " ";
595  outfile << cnt + (j-1)*np0*np1 + k*np0 + l << " ";
596 
597  outfile << cnt + j*np0*np1 + (k-1)*np0 + l << " ";
598  outfile << cnt + j*np0*np1 + (k-1)*np0 + l +1 << " ";
599  outfile << cnt + j*np0*np1 + k*np0 + l +1 << " ";
600  outfile << cnt + j*np0*np1 + k*np0 + l << endl;
601  }
602  }
603  }
604  cnt += np0*np1*np2;
605  }
606  else
607  {
608  ASSERTL0(false,"Not set up for this dimension");
609  }
610 
611  }
612 }
613 
614 }
615 }
616 
617 
618 
619 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
Abstract base class for output modules.
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
void WriteTecplotZone(std::ofstream &outfile)
FieldSharedPtr m_f
Field object.
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:248
void WriteTecplotHeader(std::ofstream &outfile, std::string var)
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:677
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:249
void WriteTecplotField(const int field, std::ofstream &outfile)
void WriteTecplotConnectivity(std::ofstream &outfile)
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215