Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 <iomanip>
37 #include <set>
38 #include <string>
39 using namespace std;
40 
43 
44 #include "OutputTecplot.h"
45 
46 namespace Nektar
47 {
48 namespace FieldUtils
49 {
50 
51 std::string TecplotZoneTypeMap[] = {
52  "ORDERED",
53  "LINESEG",
54  "TRIANGLE",
55  "QUADRILATERAL",
56  "TETRAHEDRON",
57  "BRICK",
58  "POLYGON",
59  "POLYHEDRON"
60 };
61 
62 ModuleKey OutputTecplot::m_className =
64  ModuleKey(eOutputModule, "dat"),
65  OutputTecplot::create,
66  "Writes a Tecplot file.");
67 ModuleKey OutputTecplotBinary::m_className =
69  ModuleKey(eOutputModule, "plt"),
70  OutputTecplotBinary::create,
71  "Writes a Tecplot file in binary plt format.");
72 
73 OutputTecplot::OutputTecplot(FieldSharedPtr f) : OutputModule(f),
74  m_binary(false),
75  m_oneOutputFile(false)
76 {
77  if (!f->m_setUpEquiSpacedFields)
78  {
79  m_requireEquiSpaced = true;
80  }
81  m_config["writemultiplefiles"] =
82  ConfigOption(true,"0","Write multiple files in parallel");
83 }
84 
86 {
87 }
88 
89 /**
90  * @brief Helper function to write binary data to stream.
91  */
92 template<typename T> void WriteStream(std::ostream &outfile, T data)
93 {
94  T tmp = data;
95  outfile.write(reinterpret_cast<char *>(&tmp), sizeof(T));
96 }
97 
98 /**
99  * @brief Specialisation of WriteStream to support writing std::string.
100  *
101  * Tecplot binary formats represent all strings by writing out their characters
102  * as 32-bit integers, followed by a 32-bit integer null (0) character to denote
103  * the end of the string.
104  */
105 template<> void WriteStream(std::ostream &outfile, std::string data)
106 {
107  // Convert string to array of int32_t
108  for (std::string::size_type i = 0; i < data.size(); ++i)
109  {
110  char strChar = data[i];
111  NekInt32 strCharInt = strChar;
112  WriteStream(outfile, strCharInt);
113  }
114 
115  // Now dump out zero character to terminate
116  WriteStream(outfile, 0);
117 }
118 
119 /**
120  * @brief Specialisation of WriteStream to support writing Nektar::Array
121  * datatype.
122  */
123 template<typename T> void WriteStream(std::ostream &outfile,
124  Array<OneD, T> data)
125 {
126  outfile.write(reinterpret_cast<char *>(&data[0]),
127  data.num_elements() * sizeof(T));
128 }
129 
130 /**
131  * @brief Specialisation of WriteStream to support writing std::vector datatype.
132  */
133 template<typename T> void WriteStream(std::ostream &outfile,
134  std::vector<T> data)
135 {
136  outfile.write(reinterpret_cast<char *>(&data[0]),
137  data.size() * sizeof(T));
138 }
139 
140 /**
141  * @brief Set up member variables to dump Tecplot format output.
142  */
143 void OutputTecplot::Process(po::variables_map &vm)
144 {
145  LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
146 
147  m_numBlocks = 0;
148 
149  // Do nothing if no expansion defined
150  if (fPts == LibUtilities::NullPtsField && !m_f->m_exp.size())
151  {
152  return;
153  }
154 
155 
156  if (m_f->m_verbose)
157  {
158  if (m_f->m_comm->TreatAsRankZero())
159  {
160  cout << "OutputTecplot: Writing file..." << endl;
161  }
162  }
163 
164  // Extract the output filename and extension
165  string filename = m_config["outfile"].as<string>();
166 
167  int nprocs = m_f->m_comm->GetSize();
168  int rank = m_f->m_comm->GetRank();
169 
170  if(m_config["writemultiplefiles"].as<bool>())
171  {
172  m_oneOutputFile = false;
173  }
174  else
175  {
176  m_oneOutputFile = nprocs > 1;
177  }
178 
179  // Amend for parallel output if required
180  if (nprocs != 1 && !m_oneOutputFile)
181  {
182  int dot = filename.find_last_of('.');
183  string ext = filename.substr(dot, filename.length() - dot);
184  string procId = "_P" + boost::lexical_cast<std::string>(rank);
185  string start = filename.substr(0, dot);
186  filename = start + procId + ext;
187  }
188 
189  std::string coordVars[] = { "x", "y", "z" };
190  bool doError = (vm.count("error") == 1) ? true : false;
191 
192  // Open output file
193  ofstream outfile;
194 
195  if ((m_oneOutputFile && rank == 0) || !m_oneOutputFile)
196  {
197  outfile.open(filename.c_str(), m_binary ? ios::binary : ios::out);
198  }
199 
200  std::vector<std::string> var;
201  bool writeHeader = true;
202 
203  if (fPts == LibUtilities::NullPtsField)
204  {
205  // Standard tensor-product element setup.
206  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fDef =
207  m_f->m_fielddef;
208 
209  if (fDef.size())
210  {
211  var = fDef[0]->m_fields;
212  }
213 
214  // Calculate number of FE blocks
216 
217  // Calculate coordinate dimension
218  int nBases = m_f->m_exp[0]->GetExp(0)->GetNumBases();
219  MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
220 
221  m_coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
222 
223  if (HomoExpType == MultiRegions::e3DH1D)
224  {
225  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
226  if (nPlanes == 1) // halfMode case
227  {
228  // do nothing
229  }
230  else
231  {
232  nBases += 1;
233  m_coordim += 1;
234  NekDouble tmp = m_numBlocks * (nPlanes - 1);
235  m_numBlocks = (int)tmp;
236  }
237  }
238  else if (HomoExpType == MultiRegions::e3DH2D)
239  {
240  nBases += 2;
241  m_coordim += 2;
242  }
243 
244  var.insert(var.begin(), coordVars, coordVars + m_coordim);
245 
246  m_zoneType = (TecplotZoneType)(2*(nBases-1) + 1);
247 
248  // Calculate connectivity
250 
251  // Set up storage for output fields
253 
254  // Get coordinates
255  int totpoints = m_f->m_exp[0]->GetTotPoints();
256 
257  for (int i = 0; i < m_coordim; ++i)
258  {
259  m_fields[i] = Array<OneD, NekDouble>(totpoints);
260  }
261 
262  if (m_coordim == 1)
263  {
264  m_f->m_exp[0]->GetCoords(m_fields[0]);
265  }
266  else if (m_coordim == 2)
267  {
268  m_f->m_exp[0]->GetCoords(m_fields[0], m_fields[1]);
269  }
270  else
271  {
272  m_f->m_exp[0]->GetCoords(m_fields[0], m_fields[1], m_fields[2]);
273  }
274 
275  if (var.size() > m_coordim)
276  {
277  // Backward transform all data
278  for (int i = 0; i < m_f->m_exp.size(); ++i)
279  {
280  if (m_f->m_exp[i]->GetPhysState() == false)
281  {
282  m_f->m_exp[i]->BwdTrans(m_f->m_exp[i]->GetCoeffs(),
283  m_f->m_exp[i]->UpdatePhys());
284  }
285  }
286 
287  // Add references to m_fields
288  for (int i = 0; i < m_f->m_exp.size(); ++i)
289  {
290  m_fields[i + m_coordim] = m_f->m_exp[i]->UpdatePhys();
291  }
292  }
293 
294  // Dump L2 errors of fields.
295  if (doError)
296  {
297  for (int i = 0; i < m_fields.num_elements(); ++i)
298  {
299  NekDouble l2err = m_f->m_exp[0]->L2(m_fields[i]);
300  if (rank == 0)
301  {
302  cout << "L 2 error (variable " << var[i] << ") : "
303  << l2err << endl;
304  }
305  }
306  }
307  }
308  else
309  {
310  m_coordim = fPts->GetDim();
311 
312  if (fPts->GetNpoints() == 0)
313  {
314  return;
315  }
316 
317  // Grab connectivity information.
318  fPts->GetConnectivity(m_conn);
319 
320  // Get field names
321  var = fPts->GetFieldNames();
322  var.insert(var.begin(), coordVars, coordVars + m_coordim);
323 
324  switch (fPts->GetPtsType())
325  {
328  m_numPoints.resize(1);
329  m_numPoints[0] = fPts->GetNpoints();
331  break;
333  m_numPoints.resize(2);
334  m_numPoints[0] = fPts->GetPointsPerEdge(0);
335  m_numPoints[1] = fPts->GetPointsPerEdge(1);
337  break;
339  m_numPoints.resize(3);
340  m_numPoints[0] = fPts->GetPointsPerEdge(0);
341  m_numPoints[1] = fPts->GetPointsPerEdge(1);
342  m_numPoints[2] = fPts->GetPointsPerEdge(2);
344  break;
346  {
348  for (int i = 0; i < m_conn.size(); ++i)
349  {
350  m_numBlocks += m_conn[i].num_elements() / 3;
351  }
352  break;
353  }
355  {
357  for (int i = 0; i < m_conn.size(); ++i)
358  {
359  m_numBlocks += m_conn[i].num_elements() / 4;
360  }
361  break;
362  }
363  default:
364  ASSERTL0(false, "This points type is not supported yet.");
365  }
366 
367  // Get fields and coordinates
369 
370  // We can just grab everything from points. This should be a
371  // reference, not a copy.
372  fPts->GetPts(m_fields);
373 
374  // Only write header if we're root or FE block; binary files always
375  // write header
376  writeHeader = (m_zoneType != eOrdered || rank == 0) || m_binary;
377 
378  if (doError)
379  {
380  NekDouble l2err;
381  for (int i = 0; i < m_fields.num_elements(); ++i)
382  {
383  // calculate rms value
384  int npts = m_fields[i].num_elements();
385 
386  l2err = 0.0;
387  for (int j = 0; j < npts; ++j)
388  {
389  l2err += m_fields[i][j] * m_fields[i][j];
390  }
391 
392  m_f->m_comm->AllReduce(l2err, LibUtilities::ReduceSum);
393  m_f->m_comm->AllReduce(npts, LibUtilities::ReduceSum);
394 
395  l2err /= npts;
396  l2err = sqrt(l2err);
397 
398  if (rank == 0)
399  {
400  cout << "L 2 error (variable " << var[i] << ") : "
401  << l2err << endl;
402  }
403  }
404  }
405  }
406 
407  if (m_oneOutputFile)
408  {
409  // Reduce on number of blocks and number of points.
410  m_f->m_comm->AllReduce(m_numBlocks, LibUtilities::ReduceSum);
411  for (int i = 0; i < m_numPoints.size(); ++i)
412  {
413  m_f->m_comm->AllReduce(m_numPoints[i], LibUtilities::ReduceSum);
414  }
415 
416  // Root process needs to know how much data everyone else has for
417  // writing in parallel.
418  m_rankFieldSizes = Array<OneD, int>(nprocs, 0);
419  m_rankConnSizes = Array<OneD, int>(nprocs, 0);
420  m_rankFieldSizes[rank] = m_fields[0].num_elements();
421 
422  m_totConn = 0;
423  for (int i = 0; i < m_conn.size(); ++i)
424  {
425  m_totConn += m_conn[i].num_elements();
426  }
427 
428  m_rankConnSizes[rank] = m_totConn;
429 
430  m_f->m_comm->AllReduce(m_rankFieldSizes, LibUtilities::ReduceSum);
431  m_f->m_comm->AllReduce(m_rankConnSizes, LibUtilities::ReduceSum);
432  }
433 
434  if (writeHeader)
435  {
436  WriteTecplotHeader(outfile, var);
437  }
438 
439  // Write zone data.
440  WriteTecplotZone(outfile);
441 
442  // If we're a FE block format, write connectivity (m_conn will be empty for
443  // point data).
444  WriteTecplotConnectivity(outfile);
445 
446  if ((m_oneOutputFile && rank == 0) || !m_oneOutputFile)
447  {
448  cout << "Written file: " << filename << endl;
449  }
450 }
451 
452 /**
453  * @brief Write Tecplot files header
454  *
455  * @param outfile Output file name
456  * @param var Variables names
457  */
458 void OutputTecplot::WriteTecplotHeader(std::ofstream &outfile,
459  std::vector<std::string> &var)
460 {
461  outfile << "Variables = " << var[0];
462 
463  for (int i = 1; i < var.size(); ++i)
464  {
465  outfile << ", " << var[i];
466  }
467 
468  outfile << std::endl << std::endl;
469 }
470 
471 /**
472  * @brief Write Tecplot files header in binary format
473  *
474  * @param outfile Output file name
475  * @param var Variables names
476  */
477 void OutputTecplotBinary::WriteTecplotHeader(std::ofstream &outfile,
478  std::vector<std::string> &var)
479 {
480  if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
481  {
482  return;
483  }
484 
485  // Version number
486  outfile << "#!TDV112";
487 
488  // Int value of 1 for endian check
489  WriteStream(outfile, 1);
490 
491  // We'll probably write a full solution field
492  WriteStream(outfile, 0);
493 
494  // Title
495  std::string title = "";
496  WriteStream(outfile, title);
497 
498  // Number of variables
499  WriteStream(outfile, (int)var.size());
500 
501  for (int i = 0; i < var.size(); ++i)
502  {
503  WriteStream(outfile, var[i]);
504  }
505 }
506 
507 
508 /**
509  * Write Tecplot zone output in ASCII
510  *
511  * @param outfile Output file name.
512  * @param expansion Expansion that is considered
513  */
514 void OutputTecplot::WriteTecplotZone(std::ofstream &outfile)
515 {
516  // Write either points or finite element block
517  if (m_zoneType != eOrdered)
518  {
519  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) || !m_oneOutputFile)
520  {
521  // Number of points in zone
522  int nPoints = m_oneOutputFile ?
523  Vmath::Vsum(m_f->m_comm->GetSize(), m_rankFieldSizes, 1) :
524  m_fields[0].num_elements();
525 
526  outfile << "Zone, N=" << nPoints << ", E="
527  << m_numBlocks << ", F=FEBlock, ET="
528  << TecplotZoneTypeMap[m_zoneType] << std::endl;
529  }
530 
531 
532  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
533  {
534  for (int j = 0; j < m_fields.num_elements(); ++j)
535  {
536  for (int i = 0; i < m_fields[j].num_elements(); ++i)
537  {
538  if ((!(i % 1000)) && i)
539  {
540  outfile << std::endl;
541  }
542  outfile << m_fields[j][i] << " ";
543  }
544 
545  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
546  {
548  m_f->m_comm->Recv(n, tmp);
549 
550  for (int i = 0; i < m_rankFieldSizes[n]; ++i)
551  {
552  if ((!(i % 1000)) && i)
553  {
554  outfile << std::endl;
555  }
556  outfile << tmp[i] << " ";
557  }
558  }
559  outfile << std::endl;
560  }
561  }
562  else if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
563  {
564  for (int i = 0; i < m_fields.num_elements(); ++i)
565  {
566  m_f->m_comm->Send(0, m_fields[i]);
567  }
568  }
569  else
570  {
571  // Write out coordinates and field data: ordered by field
572  // and then its data.
573  for (int j = 0; j < m_fields.num_elements(); ++j)
574  {
575  for (int i = 0; i < m_fields[j].num_elements(); ++i)
576  {
577  if ((!(i % 1000)) && i)
578  {
579  outfile << std::endl;
580  }
581  outfile << m_fields[j][i] << " ";
582  }
583  outfile << std::endl;
584  }
585  }
586  }
587  else
588  {
589  std::string dirs[] = { "I", "J", "K" };
590  outfile << "Zone";
591  for (int i = 0; i < m_numPoints.size(); ++i)
592  {
593  outfile << ", " << dirs[i] << "=" << m_numPoints[i];
594  }
595  outfile << ", F=POINT" << std::endl;
596 
597  // Write out coordinates and field data: ordered by each point then each
598  // field.
599  for (int i = 0; i < m_fields[0].num_elements(); ++i)
600  {
601  for (int j = 0; j < m_fields.num_elements(); ++j)
602  {
603  outfile << setw(12) << m_fields[j][i] << " ";
604  }
605  outfile << std::endl;
606  }
607  }
608 }
609 
610 /**
611  * @brief Write either double-precision or single-precision output of field
612  * data.
613  *
614  * @param outfile Output file name.
615  * @param expansion Expansion that is considered
616  */
617 void OutputTecplotBinary::WriteDoubleOrFloat(std::ofstream &outfile,
619 {
620  // Data format: either double or single depending on user options
621  bool useDoubles = m_config["double"].m_beenSet;
622 
623  if (useDoubles)
624  {
625  // For doubles, we can just write data.
626  WriteStream(outfile, data);
627  }
628  else
629  {
630  // For single precision, needs typecast first.
631  int nPts = data.num_elements();
632  vector<float> tmp(data.num_elements());
633  std::copy(&data[0], &data[0] + nPts, &tmp[0]);
634  WriteStream(outfile, tmp);
635  }
636 }
637 
638 /**
639  * Write Tecplot zone output in binary
640  *
641  * @param outfile Output file name.
642  * @param expansion Expansion that is considered
643  */
644 void OutputTecplotBinary::WriteTecplotZone(std::ofstream &outfile)
645 {
646  Array<OneD, NekDouble> fieldMin(m_fields.num_elements());
647  Array<OneD, NekDouble> fieldMax(m_fields.num_elements());
648 
649  // Data format: either double or single depending on user options
650  bool useDoubles = m_config["double"].m_beenSet;
651 
652  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) || !m_oneOutputFile)
653  {
654  // Don't bother naming zone
655  WriteStream(outfile, 299.0f); // Zone marker
656 
657  // Write same name as preplot
658  int rank = m_f->m_comm->GetRank();
659  string zonename = "ZONE " + boost::lexical_cast<string>(rank);
660  WriteStream(outfile, zonename);
661 
662  WriteStream(outfile, -1); // No parent zone
663  WriteStream(outfile, -1); // No strand ID
664  WriteStream(outfile, 0.0); // Solution time
665  WriteStream(outfile, -1); // Unused, set to -1
666 
667  // Zone type: 1 = lineseg, 3 = quad, 5 = brick
668  WriteStream(outfile, (int)m_zoneType);
669 
670  WriteStream(outfile, 0); // Data at nodes
671  WriteStream(outfile, 0); // No 1-1 face neighbours
672  WriteStream(outfile, 0); // No user-defined connections
673 
674  if (m_zoneType == eOrdered)
675  {
676  for (int i = 0; i < m_numPoints.size(); ++i)
677  {
678  WriteStream(outfile, m_numPoints[i]);
679  }
680 
681  for (int i = m_numPoints.size(); i < 3; ++i)
682  {
683  WriteStream(outfile, 0);
684  }
685  }
686  else
687  {
688  // Number of points in zone
689  int nPoints = m_oneOutputFile ?
690  Vmath::Vsum(m_f->m_comm->GetSize(), m_rankFieldSizes, 1) :
691  m_fields[0].num_elements();
692 
693  WriteStream(outfile, nPoints); // Total number of points
694  WriteStream(outfile, m_numBlocks); // Number of blocks
695  WriteStream(outfile, 0); // Unused
696  WriteStream(outfile, 0); // Unused
697  WriteStream(outfile, 0); // Unused
698  }
699 
700  WriteStream(outfile, 0); // No auxiliary data names
701 
702  // Finalise header
703  WriteStream(outfile, 357.0f);
704 
705  // Now start to write data section so that we can dump geometry
706  // information
707 
708  // Data marker
709  WriteStream(outfile, 299.0f);
710 
711  for (int j = 0; j < m_fields.num_elements(); ++j)
712  {
713  WriteStream(outfile, useDoubles ? 2 : 1);
714  }
715 
716  // No passive variables or variable sharing, no zone connectivity
717  // sharing (we only dump one zone)
718  WriteStream(outfile, 0);
719  WriteStream(outfile, 0);
720  WriteStream(outfile, -1);
721  }
722 
723  for (int i = 0; i < m_fields.num_elements(); ++i)
724  {
725  fieldMin[i] = Vmath::Vmin(m_fields[i].num_elements(), m_fields[i], 1);
726  fieldMax[i] = Vmath::Vmax(m_fields[i].num_elements(), m_fields[i], 1);
727  }
728 
729  m_f->m_comm->AllReduce(fieldMin, LibUtilities::ReduceMin);
730  m_f->m_comm->AllReduce(fieldMax, LibUtilities::ReduceMax);
731 
732  // Write out min/max of field data
733  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) || !m_oneOutputFile)
734  {
735  for (int i = 0; i < m_fields.num_elements(); ++i)
736  {
737  WriteStream(outfile, fieldMin[i]);
738  WriteStream(outfile, fieldMax[i]);
739  }
740  }
741 
742  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
743  {
744  for (int i = 0; i < m_fields.num_elements(); ++i)
745  {
746  WriteDoubleOrFloat(outfile, m_fields[i]);
747 
748  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
749  {
751  m_f->m_comm->Recv(n, tmp);
752  WriteDoubleOrFloat(outfile, tmp);
753  }
754  }
755  }
756  else if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
757  {
758  for (int i = 0; i < m_fields.num_elements(); ++i)
759  {
760  m_f->m_comm->Send(0, m_fields[i]);
761  }
762  }
763  else
764  {
765  for (int i = 0; i < m_fields.num_elements(); ++i)
766  {
767  WriteDoubleOrFloat(outfile, m_fields[i]);
768  }
769  }
770 }
771 
772 /**
773  * @brief Write Tecplot connectivity information (ASCII)
774  *
775  * @param outfile Output file
776  */
777 void OutputTecplot::WriteTecplotConnectivity(std::ofstream &outfile)
778 {
779  // Ordered data have no connectivity information.
780  if (m_zoneType == eOrdered)
781  {
782  return;
783  }
784 
785  if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
786  {
787  // Need to amalgamate connectivity information
789  for (int i = 0, cnt = 0; i < m_conn.size(); ++i)
790  {
791  Vmath::Vcopy(m_conn[i].num_elements(), &m_conn[i][0], 1,
792  &conn[cnt], 1);
793  cnt += m_conn[i].num_elements();
794  }
795  m_f->m_comm->Send(0, conn);
796  }
797  else
798  {
799  int cnt = 1;
800  for (int i = 0; i < m_conn.size(); ++i)
801  {
802  const int nConn = m_conn[i].num_elements();
803  for (int j = 0; j < nConn; ++j,++cnt)
804  {
805  outfile << m_conn[i][j] + 1 << " ";
806  if (!(cnt % 1000))
807  {
808  outfile << std::endl;
809  }
810  }
811  }
812  outfile << endl;
813 
814  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
815  {
816  int offset = m_rankFieldSizes[0];
817 
818  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
819  {
821  m_f->m_comm->Recv(n, conn);
822  for (int j = 0; j < conn.num_elements(); ++j)
823  {
824  outfile << conn[j] + offset + 1 << " ";
825  if ((!(j % 1000)) && j)
826  {
827  outfile << std::endl;
828  }
829  }
830  offset += m_rankFieldSizes[n];
831  }
832  }
833  }
834 }
835 
837 {
838  if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
839  {
840  // Need to amalgamate connectivity information
842  for (int i = 0, cnt = 0; i < m_conn.size(); ++i)
843  {
844  Vmath::Vcopy(m_conn[i].num_elements(), &m_conn[i][0], 1,
845  &conn[cnt], 1);
846  cnt += m_conn[i].num_elements();
847  }
848  m_f->m_comm->Send(0, conn);
849  }
850  else
851  {
852  for (int i = 0; i < m_conn.size(); ++i)
853  {
854  WriteStream(outfile, m_conn[i]);
855  }
856 
857  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
858  {
859  int offset = m_rankFieldSizes[0];
860 
861  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
862  {
864  m_f->m_comm->Recv(n, conn);
865 
866  for (int j = 0; j < conn.num_elements(); ++j)
867  {
868  conn[j] += offset;
869  }
870 
871  WriteStream(outfile, conn);
872  offset += m_rankFieldSizes[n];
873  }
874  }
875  }
876 }
877 
878 /**
879  * @brief Calculate number of Tecplot blocks.
880  *
881  * @param outfile Output file
882  */
884 {
885  int returnval = 0;
886 
887  if (m_f->m_exp[0]->GetExp(0)->GetNumBases() == 1)
888  {
889  for (int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
890  {
891  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0) - 1);
892  }
893  }
894  else if (m_f->m_exp[0]->GetExp(0)->GetNumBases() == 2)
895  {
896  for (int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
897  {
898  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0) - 1) *
899  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1) - 1);
900  }
901  }
902  else
903  {
904  for (int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
905  {
906  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0) - 1) *
907  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1) - 1) *
908  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(2) - 1);
909  }
910  }
911 
912  return returnval;
913 }
914 
915 /**
916  * @brief Calculate connectivity information for each expansion dimension.
917  *
918  * @param outfile Output file
919  */
921 {
922  int i, j, k, l;
923  int nbase = m_f->m_exp[0]->GetExp(0)->GetNumBases();
924  int cnt = 0;
925 
926  m_conn.resize(m_f->m_exp[0]->GetNumElmts());
927 
928  for (i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
929  {
930  cnt = m_f->m_exp[0]->GetPhys_Offset(i);
931 
932  if (nbase == 1)
933  {
934  int cnt2 = 0;
935  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
936 
937  Array<OneD, int> conn(2 * (np0 - 1));
938 
939  for (k = 1; k < np0; ++k)
940  {
941  conn[cnt2++] = cnt + k;
942  conn[cnt2++] = cnt + k - 1;
943  }
944 
945  m_conn[i] = conn;
946  }
947  else if (nbase == 2)
948  {
949  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
950  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
951  int totPoints = m_f->m_exp[0]->GetTotPoints();
952  int nPlanes = 1;
953  int cnt2 = 0;
954 
955  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
956  {
957  nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
958 
959  // default to 2D case for HalfMode when nPlanes = 1
960  if (nPlanes > 1)
961  {
962  totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
963 
964  Array<OneD, int> conn(8 * (np1 - 1) * (np0 - 1) * (nPlanes - 1));
965 
966  for (int n = 1; n < nPlanes; ++n)
967  {
968  for (j = 1; j < np1; ++j)
969  {
970  for (k = 1; k < np0; ++k)
971  {
972  conn[cnt2++] = cnt + (n - 1) * totPoints +
973  (j - 1) * np0 + k - 1;
974  conn[cnt2++] = cnt + (n - 1) * totPoints +
975  (j - 1) * np0 + k;
976  conn[cnt2++] = cnt + (n - 1) * totPoints +
977  j * np0 + k;
978  conn[cnt2++] = cnt + (n - 1) * totPoints +
979  j * np0 + k - 1;
980  conn[cnt2++] = cnt + n * totPoints +
981  (j - 1) * np0 + k - 1;
982  conn[cnt2++] = cnt + n * totPoints +
983  (j - 1) * np0 + k;
984  conn[cnt2++] = cnt + n * totPoints +
985  j * np0 + k;
986  conn[cnt2++] = cnt + n * totPoints +
987  j * np0 + k - 1;
988  }
989  }
990  }
991  m_conn[i] = conn;
992  }
993  }
994 
995  if (nPlanes == 1)
996  {
997  Array<OneD, int> conn(4 * (np0 - 1) * (np1 - 1));
998  for (j = 1; j < np1; ++j)
999  {
1000  for (k = 1; k < np0; ++k)
1001  {
1002  conn[cnt2++] = cnt + (j - 1) * np0 + k - 1;
1003  conn[cnt2++] = cnt + (j - 1) * np0 + k;
1004  conn[cnt2++] = cnt + j * np0 + k;
1005  conn[cnt2++] = cnt + j * np0 + k - 1;
1006  }
1007  }
1008  m_conn[i] = conn;
1009  }
1010  }
1011  else if (nbase == 3)
1012  {
1013  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
1014  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
1015  int np2 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(2);
1016  int cnt2 = 0;
1017 
1018  Array<OneD, int> conn(8 * (np0 - 1) * (np1 - 1) * (np2 - 1));
1019 
1020  for (j = 1; j < np2; ++j)
1021  {
1022  for (k = 1; k < np1; ++k)
1023  {
1024  for (l = 1; l < np0; ++l)
1025  {
1026  conn[cnt2++] =
1027  cnt + (j - 1) * np0 * np1 + (k - 1) * np0 + l - 1;
1028  conn[cnt2++] =
1029  cnt + (j - 1) * np0 * np1 + (k - 1) * np0 + l;
1030  conn[cnt2++] =
1031  cnt + (j - 1) * np0 * np1 + k * np0 + l;
1032  conn[cnt2++] =
1033  cnt + (j - 1) * np0 * np1 + k * np0 + l - 1;
1034  conn[cnt2++] =
1035  cnt + j * np0 * np1 + (k - 1) * np0 + l - 1;
1036  conn[cnt2++] =
1037  cnt + j * np0 * np1 + (k - 1) * np0 + l;
1038  conn[cnt2++] =
1039  cnt + j * np0 * np1 + k * np0 + l;
1040  conn[cnt2++] =
1041  cnt + j * np0 * np1 + k * np0 + l - 1;
1042  }
1043  }
1044  }
1045 
1046  m_conn[i] = conn;
1047  }
1048  else
1049  {
1050  ASSERTL0(false, "Not set up for this dimension");
1051  }
1052 
1053  }
1054 }
1055 }
1056 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Represents a command-line configuration option.
TecplotZoneType m_zoneType
Tecplot zone type of output.
Definition: OutputTecplot.h:82
boost::int32_t NekInt32
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:871
Array< OneD, int > m_rankConnSizes
Each rank's connectivity sizes.
Definition: OutputTecplot.h:96
virtual void WriteTecplotHeader(std::ofstream &outfile, std::vector< std::string > &var)
Write Tecplot files header in binary format.
STL namespace.
pair< ModuleType, string > ModuleKey
bool m_oneOutputFile
True if writing a single output file.
Definition: OutputTecplot.h:80
void WriteDoubleOrFloat(std::ofstream &outfile, Array< OneD, NekDouble > &data)
Write either double-precision or single-precision output of field data.
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:178
virtual void WriteTecplotConnectivity(std::ofstream &outfile)
Write Tecplot connectivity information (ASCII)
void CalculateConnectivity()
Calculate connectivity information for each expansion dimension.
virtual void WriteTecplotHeader(std::ofstream &outfile, std::vector< std::string > &var)
Write Tecplot files header.
Array< OneD, int > m_rankFieldSizes
Each rank's field sizes.
Definition: OutputTecplot.h:94
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:767
virtual void WriteTecplotZone(std::ofstream &outfile)
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
vector< Array< OneD, int > > m_conn
Connectivty for each block: one per element.
Definition: OutputTecplot.h:92
int m_numBlocks
Number of blocks in Tecplot file.
Definition: OutputTecplot.h:86
int m_totConn
Total number of connectivity entries.
Definition: OutputTecplot.h:90
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:179
std::string TecplotZoneTypeMap[]
int GetNumTecplotBlocks()
Calculate number of Tecplot blocks.
virtual void Process(po::variables_map &vm)
Set up member variables to dump Tecplot format output.
int m_coordim
Coordinate dimension of output.
Definition: OutputTecplot.h:88
Array< OneD, Array< OneD, NekDouble > > m_fields
Field data to output.
Definition: OutputTecplot.h:98
bool m_binary
True if writing binary field output.
Definition: OutputTecplot.h:78
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:737
void WriteStream(std::ostream &outfile, T data)
Helper function to write binary data to stream.
virtual void WriteTecplotConnectivity(std::ofstream &outfile)
Write Tecplot connectivity information (ASCII)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
virtual void WriteTecplotZone(std::ofstream &outfile)
vector< int > m_numPoints
Number of points per block in Tecplot file.
Definition: OutputTecplot.h:84
Abstract base class for output modules.
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215