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 
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 
59  if(f->m_setUpEquiSpacedFields)
60  {
62  }
63  else
64  {
65  m_requireEquiSpaced = true;
67  }
68 
69 }
70 
72 {
73 }
74 
75 void OutputTecplot::Process(po::variables_map &vm)
76 {
77  LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
78 
79  m_doError = (vm.count("error") == 1)? true: false;
80 
81  // Do nothing if no expansion defined
82  if (fPts == LibUtilities::NullPtsField && !m_f->m_exp.size())
83  {
84  return;
85  }
86 
87  if (m_f->m_verbose)
88  {
89  cout << "OutputTecplot: Writing file..." << endl;
90  }
91 
92 
93  // Extract the output filename and extension
94  string filename = m_config["outfile"].as<string>();
95 
96  int nprocs = m_f->m_comm->GetSize();
97  int rank = m_f->m_comm->GetRank();
98  // Amend for parallel output if required
99  if(nprocs != 1)
100  {
101  int dot = filename.find_last_of('.');
102  string ext = filename.substr(dot,filename.length()-dot);
103  string procId = "_P" + boost::lexical_cast<std::string>(rank);
104  string start = filename.substr(0,dot);
105  filename = start + procId + ext;
106  }
107 
108  if(fPts != LibUtilities::NullPtsField)
109  {
110  int i = 0;
111  int j = 0;
112  int dim = fPts->GetDim();
113 
114  if(fPts->GetNpoints() == 0)
115  {
116  return;
117  }
118 
119  // Write solution.
120  ofstream outfile(filename.c_str());
121 
122  // points type
123  LibUtilities::PtsType pType = fPts->GetPtsType();
124 
125  vector<Array<OneD, int> > ptsConn;
126  fPts->GetConnectivity(ptsConn);
127 
128  // only dump header info for all proces if ptsType is for
129  // TriBlock or TetBlock
130  if((pType > LibUtilities::ePtsBox)||(rank == 0))
131  {
132  switch(dim)
133  {
134  case 1:
135  outfile << "VARIABLES = x";
136  break;
137  case 2:
138  outfile << "VARIABLES = x,y";
139  break;
140  case 3:
141  outfile << "VARIABLES = x,y,z";
142  break;
143  }
144 
145  for(i = 0; i < fPts->GetNFields(); ++i)
146  {
147  outfile << "," << fPts->GetFieldName(i);
148  }
149  outfile << endl;
150  }
151 
152  bool DumpAsFEPoint = true;
153  switch(fPts->GetPtsType())
154  {
157  {
158  if(rank == 0)
159  {
160  outfile << " ZONE I="
161  << fPts->GetNpoints()
162  << " F=POINT" << endl;
163  }
164  break;
165  }
167  {
168  if(rank == 0)
169  {
170  outfile << " ZONE I=" << fPts->GetPointsPerEdge(0)
171  << " J=" << fPts->GetPointsPerEdge(1)
172  << " F=POINT" << endl;
173  }
174  break;
175  }
177  {
178  if(rank == 0)
179  {
180  outfile << " ZONE I=" << fPts->GetPointsPerEdge(0)
181  << " J=" << fPts->GetPointsPerEdge(1)
182  << " K=" << fPts->GetPointsPerEdge(2)
183  << " F=POINT" << endl;
184  }
185  break;
186  }
187  break;
189  {
190  int numBlocks = 0;
191  for(i = 0; i < ptsConn.size(); ++i)
192  {
193  numBlocks +=
194  ptsConn[i].num_elements()/3;
195  }
196  outfile << "Zone, N="
197  << fPts->GetNpoints()
198  << ", E=" << numBlocks
199  << ", F=FEBlock" << ", ET=TRIANGLE"
200  << std::endl;
201  DumpAsFEPoint = false;
202  break;
203  }
205  {
206  int numBlocks = 0;
207  for(i = 0; i < ptsConn.size(); ++i)
208  {
209  numBlocks +=
210  ptsConn[i].num_elements()/4;
211  }
212  outfile << "Zone, N="
213  << fPts->GetNpoints()
214  << ", E=" << numBlocks
215  << ", F=FEBlock" << ", ET=TETRAHEDRON"
216  << std::endl;
217  DumpAsFEPoint = false;
218  break;
219  }
220  default:
221  ASSERTL0(false, "ptsType not supported yet.");
222  }
223 
224  if(DumpAsFEPoint) // dump in point format
225  {
226  for(i = 0; i < fPts->GetNpoints(); ++i)
227  {
228  for(j = 0; j < dim; ++j)
229  {
230  outfile << std::setw(12)
231  << fPts->GetPointVal(j, i) << " ";
232  }
233 
234  for(j = 0; j < fPts->GetNFields(); ++j)
235  {
236  outfile << std::setw(12)
237  << m_f->m_data[j][i] << " ";
238  }
239  outfile << endl;
240  }
241 
242 
243 
244  if (m_doError)
245  {
246  NekDouble l2err;
247  std::string coordval[] = {"x","y","z"};
248  int rank = m_f->m_comm->GetRank();
249 
250  for(i = 0; i < dim; ++i)
251  {
252  // calculate rms value
253  l2err = 0.0;
254  for(j = 0; j < fPts->GetNpoints(); ++j)
255  {
256  l2err += fPts->GetPointVal(i,j)*fPts->GetPointVal(i,j);
257  }
258  m_f->m_comm->AllReduce(l2err, LibUtilities::ReduceSum);
259 
260  int npts = fPts->GetNpoints();
261  m_f->m_comm->AllReduce(npts, LibUtilities::ReduceSum);
262 
263  l2err /= npts;
264  l2err = sqrt(l2err);
265 
266  if(rank == 0)
267  {
268  cout << "L 2 error (variable "
269  << coordval[i] << ") : "
270  << l2err << endl;
271  }
272  }
273 
274  for(i = 0; i < fPts->GetNFields(); ++i)
275  {
276  // calculate rms value
277  l2err = 0.0;
278  for(j = 0; j < fPts->GetNpoints(); ++j)
279  {
280  l2err += m_f->m_data[i][j]*m_f->m_data[i][j];
281  }
282  m_f->m_comm->AllReduce(l2err,
284 
285  int npts = fPts->GetNpoints();
286  m_f->m_comm->AllReduce(npts,
288 
289  l2err /= npts;
290  l2err = sqrt(l2err);
291 
292  if(rank == 0)
293  {
294  cout << "L 2 error (variable "
295  << fPts->GetFieldName(i) << ") : "
296  << l2err << endl;
297  }
298  }
299  }
300 
301  m_f->m_comm->Block();
302  }
303  else // dump in block format
304  {
305  for(j = 0; j < dim + fPts->GetNFields(); ++j)
306  {
307  for(i = 0; i < fPts->GetNpoints(); ++i)
308  {
309  outfile << fPts->GetPointVal(j, i) << " ";
310  if((!(i % 1000))&&i)
311  {
312  outfile << std::endl;
313  }
314  }
315  outfile << endl;
316  }
317 
318  // dump connectivity data if it exists
319  for(i = 0; i < ptsConn.size();++i)
320  {
321  for(j = 0; j < ptsConn[i].num_elements(); ++j)
322  {
323  outfile << ptsConn[i][j] +1 << " ";
324  if( ( !(j % 10 * dim) ) && j )
325  {
326  outfile << std::endl;
327  }
328  }
329  }
330 
331  if (m_doError)
332  {
333  NekDouble l2err;
334  std::string coordval[] = {"x","y","z"};
335  int rank = m_f->m_comm->GetRank();
336 
337  for(int i = 0; i < dim + fPts->GetNFields(); ++i)
338  {
339  // calculate rms value
340  l2err = 0.0;
341  for(j = 0; j < fPts->GetNpoints(); ++j)
342  {
343  l2err += fPts->GetPointVal(i,j)*fPts->GetPointVal(i,j);
344  }
345  m_f->m_comm->AllReduce(l2err,
347 
348  int npts = fPts->GetNpoints();
349  m_f->m_comm->AllReduce(npts,
351 
352  l2err /= npts;
353  l2err = sqrt(l2err);
354 
355  if(rank == 0)
356  {
357  if(i < dim)
358  {
359  cout << "L 2 error (variable "
360  << coordval[i] << ") : "
361  << l2err << endl;
362  }
363  else
364  {
365  cout << "L 2 error (variable "
366  << fPts->GetFieldName(i-dim) << ") : "
367  << l2err << endl;
368  }
369  }
370  }
371  }
372  }
373  }
374  else
375  {
376 
377  // Write solution.
378  ofstream outfile(filename.c_str());
379  std::string var;
380  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fDef =
381  m_f->m_fielddef;
382  if (fDef.size())
383  {
384  var = fDef[0]->m_fields[0];
385 
386  for (int j = 1; j < fDef[0]->m_fields.size(); ++j)
387  {
388  var = var + ", " + fDef[0]->m_fields[j];
389  }
390  }
391 
392  WriteTecplotHeader(outfile,var);
393  WriteTecplotZone(outfile);
394  if(var.length()) // see if any variables are defined
395  {
396  for(int j = 0; j < m_f->m_exp.size(); ++j)
397  {
398  WriteTecplotField(j,outfile);
399  }
400  }
401 
402  WriteTecplotConnectivity(outfile);
403  }
404 
405  cout << "Written file: " << filename << endl;
406 }
407 
408 /**
409  * Write Tecplot Files Header
410  * @param outfile Output file name.
411  * @param var variables names
412  */
413 void OutputTecplot::WriteTecplotHeader(std::ofstream &outfile,
414  std::string var)
415 {
416  int coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
417  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
418 
419  if(HomoExpType == MultiRegions::e3DH1D)
420  {
421  if(m_f->m_session->DefinesSolverInfo("ModeType")&&
422  boost::iequals(m_f->m_session->GetSolverInfo("ModeType"),"HalfMode"))
423  { // turn off for half mode case
424  }
425  else
426  {
427  coordim +=1;
428  }
429  }
430  else if (HomoExpType == MultiRegions::e3DH2D)
431  {
432  coordim += 2;
433  }
434 
435  outfile << "Variables = x";
436 
437  if(coordim == 2)
438  {
439  outfile << ", y";
440  }
441  else if (coordim == 3)
442  {
443  outfile << ", y, z";
444  }
445 
446  if(var.length())
447  {
448  outfile << ", "<< var << std::endl << std::endl;
449  }
450  else
451  {
452  outfile << std::endl << std::endl;
453  }
454 }
455 
456 
457 /**
458  * Write Tecplot Files Zone
459  * @param outfile Output file name.
460  * @param expansion Expansion that is considered
461  */
462 void OutputTecplot::WriteTecplotZone(std::ofstream &outfile)
463 {
464  switch(m_outputType)
465  {
466  case eFullBlockZone: //write as full block zone
467  {
468  int i,j;
469  int coordim = m_f->m_exp[0]->GetCoordim(0);
470  int totpoints = m_f->m_exp[0]->GetTotPoints();
471  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
472 
473  Array<OneD,NekDouble> coords[3];
474 
475  coords[0] = Array<OneD,NekDouble>(totpoints);
476  coords[1] = Array<OneD,NekDouble>(totpoints);
477  coords[2] = Array<OneD,NekDouble>(totpoints);
478 
479  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
480 
481  if (m_doError)
482  {
483  NekDouble l2err;
484  std::string coordval[] = {"x","y","z"};
485  int rank = m_f->m_comm->GetRank();
486 
487  for(int i = 0; i < coordim; ++i)
488  {
489  l2err = m_f->m_exp[0]->L2(coords[i]);
490  if(rank == 0)
491  {
492  cout << "L 2 error (variable "
493  << coordval[i] << ") : " << l2err << endl;
494  }
495  }
496  }
497 
498  int numBlocks = GetNumTecplotBlocks();
499  int nBases = m_f->m_exp[0]->GetExp(0)->GetNumBases();
500 
501  if (HomoExpType == MultiRegions::e3DH1D)
502  {
503  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
504  if(nPlanes == 1) // halfMode case
505  {
506  // do nothing
507  }
508  else
509  {
510  nBases += 1;
511  coordim += 1;
512  NekDouble tmp = numBlocks * (nPlanes-1);
513  numBlocks = (int)tmp;
514  }
515  }
516  else if (HomoExpType == MultiRegions::e3DH2D)
517  {
518  nBases += 2;
519  coordim += 1;
520  }
521 
522  outfile << "Zone, N=" << totpoints << ", E="<<
523  numBlocks << ", F=FEBlock" ;
524 
525  switch(nBases)
526  {
527  case 1:
528  outfile << ", ET=LINESEG" << std::endl;
529  break;
530  case 2:
531  outfile << ", ET=QUADRILATERAL" << std::endl;
532  break;
533  case 3:
534  outfile << ", ET=BRICK" << std::endl;
535  break;
536  }
537 
538  // write out coordinates in block format
539  for(j = 0; j < coordim; ++j)
540  {
541  for(i = 0; i < totpoints; ++i)
542  {
543  outfile << coords[j][i] << " ";
544  if((!(i % 1000))&&i)
545  {
546  outfile << std::endl;
547  }
548  }
549  outfile << std::endl;
550  }
551  break;
552  }
553  case eSeperateZones:
554  {
555  for(int i = 0; i < m_f->m_exp[0]->GetExpSize(); ++i)
556  {
557  m_f->m_exp[0]->WriteTecplotZone(outfile,i);
558  }
559  break;
560  }
562  ASSERTL0(false,
563  "Should not have this option in this method");
564  break;
565  }
566 }
567 
569 {
570  int returnval = 0;
571 
572  if(m_f->m_exp[0]->GetExp(0)->GetNumBases() == 1)
573  {
574  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
575  {
576  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1);
577  }
578  }
579  else if(m_f->m_exp[0]->GetExp(0)->GetNumBases() == 2)
580  {
581  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
582  {
583  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1)*
584  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1)-1);
585  }
586  }
587  else
588  {
589  for(int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
590  {
591  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0)-1)*
592  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1)-1)*
593  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(2)-1);
594  }
595  }
596 
597  return returnval;
598 }
599 
600 
601 /**
602  * Write Tecplot Files Field
603  * @param outfile Output file name.
604  * @param expansion Expansion that is considered
605  */
606 void OutputTecplot::WriteTecplotField(const int field,
607  std::ofstream &outfile)
608 {
609 
610  if(m_outputType == eFullBlockZone) //write as full block zone
611  {
612  int totpoints = m_f->m_exp[0]->GetTotPoints();
613 
614  if(m_f->m_exp[field]->GetPhysState() == false)
615  {
616  m_f->m_exp[field]->BwdTrans(m_f->m_exp[field]->GetCoeffs(),
617  m_f->m_exp[field]->UpdatePhys());
618  }
619 
620  if (m_doError)
621  {
622  NekDouble l2err = m_f->m_exp[0]->L2(m_f->m_exp[field]->UpdatePhys());
623 
624  if(m_f->m_comm->GetRank() == 0)
625  {
626  cout << "L 2 error (variable "
627  << m_f->m_fielddef[0]->m_fields[field] << ") : "
628  << l2err << endl;
629  }
630  }
631  else
632  {
633  for(int i = 0; i < totpoints; ++i)
634  {
635  outfile << m_f->m_exp[field]->GetPhys()[i] << " ";
636  if((!(i % 1000))&&i)
637  {
638  outfile << std::endl;
639  }
640  }
641  }
642  outfile << std::endl;
643  }
644  else
645  {
646  for(int e = 0; e < m_f->m_exp[field]->GetExpSize(); ++e)
647  {
648  m_f->m_exp[field]->WriteTecplotField(outfile,e);
649  }
650  }
651 }
652 
653 void OutputTecplot::WriteTecplotConnectivity(std::ofstream &outfile)
654 {
655  int i,j,k,l;
656  int nbase = m_f->m_exp[0]->GetExp(0)->GetNumBases();
657  int cnt = 0;
658 
659  for(i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
660  {
661  cnt = m_f->m_exp[0]->GetPhys_Offset(i);
662 
663  if(nbase == 1)
664  {
665  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
666 
667  for(k = 1; k < np0; ++k)
668  {
669  outfile << cnt + k +1 << " ";
670  outfile << cnt + k << endl;
671  }
672 
673  cnt += np0;
674  }
675  else if(nbase == 2)
676  {
677  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
678  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
679  int totPoints = m_f->m_exp[0]->GetTotPoints();
680  int nPlanes = 1;
681 
682  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
683  {
684  nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
685 
686  if(nPlanes > 1) // default to 2D case for HalfMode when nPlanes = 1
687  {
688  totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
689 
690 
691  for(int n = 1; n < nPlanes; ++n)
692  {
693  for(j = 1; j < np1; ++j)
694  {
695  for(k = 1; k < np0; ++k)
696  {
697  outfile << cnt + (n-1)*totPoints + (j-1)*np0 + k
698  << " ";
699  outfile << cnt + (n-1)*totPoints + (j-1)*np0 + k + 1
700  << " ";
701  outfile << cnt + (n-1)*totPoints + j*np0 + k + 1
702  << " ";
703  outfile << cnt + (n-1)*totPoints + j*np0 + k
704  << " ";
705 
706  outfile << cnt + n*totPoints + (j-1)*np0 + k
707  << " ";
708  outfile << cnt + n*totPoints + (j-1)*np0 + k + 1
709  << " ";
710  outfile << cnt + n*totPoints + j*np0 + k + 1
711  << " ";
712  outfile << cnt + n*totPoints + j*np0 + k << endl;
713  }
714  }
715  }
716  cnt += np0*np1;
717  }
718  }
719 
720  if(nPlanes == 1)
721  {
722  for(j = 1; j < np1; ++j)
723  {
724  for(k = 1; k < np0; ++k)
725  {
726  outfile << cnt + (j-1)*np0 + k << " ";
727  outfile << cnt + (j-1)*np0 + k +1 << " ";
728  outfile << cnt + j*np0 + k +1 << " ";
729  outfile << cnt + j*np0 + k << endl;
730  }
731  }
732  }
733  }
734  else if(nbase == 3)
735  {
736  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
737  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
738  int np2 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(2);
739 
740  for(j = 1; j < np2; ++j)
741  {
742  for(k = 1; k < np1; ++k)
743  {
744  for(l = 1; l < np0; ++l)
745  {
746  outfile << cnt + (j-1)*np0*np1 + (k-1)*np0 + l << " ";
747  outfile << cnt + (j-1)*np0*np1 + (k-1)*np0 + l +1 << " ";
748  outfile << cnt + (j-1)*np0*np1 + k*np0 + l +1 << " ";
749  outfile << cnt + (j-1)*np0*np1 + k*np0 + l << " ";
750 
751  outfile << cnt + j*np0*np1 + (k-1)*np0 + l << " ";
752  outfile << cnt + j*np0*np1 + (k-1)*np0 + l +1 << " ";
753  outfile << cnt + j*np0*np1 + k*np0 + l +1 << " ";
754  outfile << cnt + j*np0*np1 + k*np0 + l << endl;
755  }
756  }
757  }
758  }
759  else
760  {
761  ASSERTL0(false,"Not set up for this dimension");
762  }
763 
764  }
765 }
766 
767 }
768 }
769 
770 
771 
772 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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:262
void WriteTecplotHeader(std::ofstream &outfile, std::string var)
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:263
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