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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Dat file format output.
32 //
33 ////////////////////////////////////////////////////////////////////////////////
34 
35 #include <iomanip>
36 #include <set>
37 #include <string>
38 using namespace std;
39 
40 #include <boost/core/ignore_unused.hpp>
41 
44 
45 #include "OutputTecplot.h"
46 
47 namespace Nektar
48 {
49 namespace FieldUtils
50 {
51 
52 std::string TecplotZoneTypeMap[] = {"ORDERED", "LINESEG", "TRIANGLE",
53  "QUADRILATERAL", "TETRAHEDRON", "BRICK",
54  "POLYGON", "POLYHEDRON"};
55 
56 ModuleKey OutputTecplot::m_className =
58  OutputTecplot::create,
59  "Writes a Tecplot file.");
60 ModuleKey OutputTecplotBinary::m_className =
62  ModuleKey(eOutputModule, "plt"), OutputTecplotBinary::create,
63  "Writes a Tecplot file in binary plt format.");
64 
65 OutputTecplot::OutputTecplot(FieldSharedPtr f)
66  : OutputFileBase(f), m_binary(false), m_oneOutputFile(false)
67 {
68  m_requireEquiSpaced = true;
69  m_config["double"] = ConfigOption(true, "0",
70  "Write double-precision format data:"
71  "more accurate but more disk space"
72  " required");
73 }
74 
76 {
77 }
78 
79 void OutputTecplot::Process(po::variables_map &vm)
80 {
81 
82  if (m_config["writemultiplefiles"].as<bool>())
83  {
84  m_oneOutputFile = false;
85  }
86  else
87  {
88  m_oneOutputFile = (m_f->m_comm->GetSize() > 1);
89  }
90 
92 }
93 
94 /**
95  * @brief Helper function to write binary data to stream.
96  */
97 template <typename T> void WriteStream(std::ostream &outfile, T data)
98 {
99  T tmp = data;
100  outfile.write(reinterpret_cast<char *>(&tmp), sizeof(T));
101 }
102 
103 /**
104  * @brief Specialisation of WriteStream to support writing std::string.
105  *
106  * Tecplot binary formats represent all strings by writing out their characters
107  * as 32-bit integers, followed by a 32-bit integer null (0) character to denote
108  * the end of the string.
109  */
110 template <> void WriteStream(std::ostream &outfile, std::string data)
111 {
112  // Convert string to array of int32_t
113  for (std::string::size_type i = 0; i < data.size(); ++i)
114  {
115  char strChar = data[i];
116  NekInt32 strCharInt = strChar;
117  WriteStream(outfile, strCharInt);
118  }
119 
120  // Now dump out zero character to terminate
121  WriteStream(outfile, 0);
122 }
123 
124 /**
125  * @brief Specialisation of WriteStream to support writing Nektar::Array
126  * datatype.
127  */
128 template <typename T>
129 void WriteStream(std::ostream &outfile, Array<OneD, T> data)
130 {
131  if (data.size())
132  {
133  outfile.write(reinterpret_cast<char *>(&data[0]),
134  data.size() * sizeof(T));
135  }
136 }
137 
138 /**
139  * @brief Specialisation of WriteStream to support writing std::vector datatype.
140  */
141 template <typename T>
142 void WriteStream(std::ostream &outfile, std::vector<T> data)
143 {
144  if (data.size())
145  {
146  outfile.write(reinterpret_cast<char *>(&data[0]),
147  data.size() * sizeof(T));
148  }
149 }
150 
151 void OutputTecplot::OutputFromPts(po::variables_map &vm)
152 {
153  LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
154 
155  // do not output if zone is empty
156  if (fPts->GetNpoints() == 0)
157  {
158  return;
159  }
160 
161  int rank = m_f->m_comm->GetRank();
162  m_numBlocks = 0;
163 
164  m_coordim = fPts->GetDim();
165 
166  // Grab connectivity information.
167  fPts->GetConnectivity(m_conn);
168 
169  switch (fPts->GetPtsType())
170  {
172  m_numPoints.resize(1);
173  m_numPoints[0] = fPts->GetNpoints();
174  m_f->m_comm->AllReduce(m_numPoints[0], LibUtilities::ReduceSum);
176  break;
178  m_numPoints.resize(1);
179  m_numPoints[0] = fPts->GetPointsPerEdge(0);
181  break;
183  m_numPoints.resize(2);
184  m_numPoints[0] = fPts->GetPointsPerEdge(0);
185  m_numPoints[1] = fPts->GetPointsPerEdge(1);
187  break;
189  m_numPoints.resize(3);
190  m_numPoints[0] = fPts->GetPointsPerEdge(0);
191  m_numPoints[1] = fPts->GetPointsPerEdge(1);
192  m_numPoints[2] = fPts->GetPointsPerEdge(2);
194  break;
196  {
198  for (int i = 0; i < m_conn.size(); ++i)
199  {
200  m_numBlocks += m_conn[i].size() / 3;
201  }
202  break;
203  }
205  {
207  for (int i = 0; i < m_conn.size(); ++i)
208  {
209  m_numBlocks += m_conn[i].size() / 4;
210  }
211  break;
212  }
213  default:
214  ASSERTL0(false, "This points type is not supported yet.");
215  }
216 
217  // Get fields and coordinates
218  m_fields = Array<OneD, Array<OneD, NekDouble>>(m_f->m_variables.size() +
219  m_coordim);
220 
221  // We can just grab everything from points. This should be a
222  // reference, not a copy.
223  fPts->GetPts(m_fields);
224 
225  // Only write header if we're root or FE block; binary files always
226  // write header
227  m_writeHeader = (m_zoneType != eOrdered || rank == 0) || m_binary;
228 
229  WriteTecplotFile(vm);
230 }
231 
232 void OutputTecplot::OutputFromExp(po::variables_map &vm)
233 {
234  m_numBlocks = 0;
235  m_writeHeader = true;
236 
237  // Calculate number of FE blocks
239 
240  // Calculate coordinate dimension
241  int nBases = m_f->m_exp[0]->GetExp(0)->GetNumBases();
242 
243  m_coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
244  int totpoints = m_f->m_exp[0]->GetTotPoints();
245 
246  if (m_f->m_numHomogeneousDir > 0)
247  {
248  int nPlanes = m_f->m_exp[0]->GetZIDs().size();
249  if (nPlanes == 1) // halfMode case
250  {
251  // do nothing
252  }
253  else
254  {
255  // If Fourier points, output extra plane to fill domain
256  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
257  {
258  nPlanes += 1;
259  totpoints += m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
260  }
261  nBases += m_f->m_numHomogeneousDir;
262  m_coordim += m_f->m_numHomogeneousDir;
263  NekDouble tmp = m_numBlocks * (nPlanes - 1);
264  m_numBlocks = (int)tmp;
265  }
266  }
267 
268  m_zoneType = (TecplotZoneType)(2 * (nBases - 1) + 1);
269 
270  // Calculate connectivity
272 
273  // Set up storage for output fields
274  m_fields = Array<OneD, Array<OneD, NekDouble>>(m_f->m_variables.size() +
275  m_coordim);
276 
277  // Get coordinates
278  for (int i = 0; i < m_coordim; ++i)
279  {
280  m_fields[i] = Array<OneD, NekDouble>(totpoints);
281  }
282 
283  if (m_coordim == 1)
284  {
285  m_f->m_exp[0]->GetCoords(m_fields[0]);
286  }
287  else if (m_coordim == 2)
288  {
289  m_f->m_exp[0]->GetCoords(m_fields[0], m_fields[1]);
290  }
291  else
292  {
293  m_f->m_exp[0]->GetCoords(m_fields[0], m_fields[1], m_fields[2]);
294  }
295 
296  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
297  {
298  // Copy values
299  for (int i = 0; i < m_f->m_variables.size(); ++i)
300  {
301  m_fields[i + m_coordim] = Array<OneD, NekDouble>(totpoints);
302  Vmath::Vcopy(m_f->m_exp[0]->GetTotPoints(),
303  m_f->m_exp[i]->UpdatePhys(), 1,
304  m_fields[i + m_coordim], 1);
305  }
306  }
307  else
308  {
309  // Add references to m_fields
310  for (int i = 0; i < m_f->m_variables.size(); ++i)
311  {
312  m_fields[i + m_coordim] = m_f->m_exp[i]->UpdatePhys();
313  }
314  }
315 
316  // If Fourier, fill extra plane with data
317  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
318  {
319  int points_on_plane = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
320  const int offset = totpoints - points_on_plane;
321  NekDouble z = m_fields[m_coordim - 1][totpoints - 2 * points_on_plane] +
322  (m_fields[m_coordim - 1][points_on_plane] -
323  m_fields[m_coordim - 1][0]);
324  // x and y
325  Array<OneD, NekDouble> tmp = m_fields[0] + offset;
326  Vmath::Vcopy(points_on_plane, m_fields[0], 1, tmp, 1);
327  tmp = m_fields[1] + offset;
328  Vmath::Vcopy(points_on_plane, m_fields[1], 1, tmp, 1);
329  // z coordinate
330  tmp = m_fields[2] + offset;
331  Vmath::Vcopy(points_on_plane, m_fields[2], 1, tmp, 1);
332  Vmath::Sadd(points_on_plane, z, m_fields[2], 1, tmp, 1);
333 
334  // variables
335  for (int i = 0; i < m_f->m_variables.size(); ++i)
336  {
337  tmp = m_fields[i + m_coordim] + offset;
338  Vmath::Vcopy(points_on_plane, m_fields[i + m_coordim], 1, tmp, 1);
339  }
340  }
341 
342  WriteTecplotFile(vm);
343 }
344 
345 void OutputTecplot::OutputFromData(po::variables_map &vm)
346 {
347  boost::ignore_unused(vm);
348 
350  "OutputTecplot can't write using only FieldData.");
351 }
352 
353 fs::path OutputTecplot::GetPath(std::string &filename, po::variables_map &vm)
354 {
355  boost::ignore_unused(vm);
356 
357  int nprocs = m_f->m_comm->GetSize();
358  string returnstr(filename);
359 
360  // Amend for parallel output if required
361  if (nprocs != 1 && !m_oneOutputFile)
362  {
363  int rank = m_f->m_comm->GetRank();
364  int dot = filename.find_last_of('.');
365  string ext = filename.substr(dot, filename.length() - dot);
366  string procId = "_P" + boost::lexical_cast<std::string>(rank);
367  string start = filename.substr(0, dot);
368  returnstr = start + procId + ext;
369  }
370  return fs::path(returnstr);
371 }
372 
373 fs::path OutputTecplot::GetFullOutName(std::string &filename,
374  po::variables_map &vm)
375 {
376  return GetPath(filename, vm);
377 }
378 
379 void OutputTecplot::WriteTecplotFile(po::variables_map &vm)
380 {
381  // Variable names
382  std::string coordVars[] = {"x", "y", "z"};
383  std::vector<string> variables = m_f->m_variables;
384  variables.insert(variables.begin(), coordVars, coordVars + m_coordim);
385 
386  int nprocs = m_f->m_comm->GetSize();
387  int rank = m_f->m_comm->GetRank();
388 
389  // Extract the output filename and extension
390  string filename = m_config["outfile"].as<string>();
391  string outFile = LibUtilities::PortablePath(GetFullOutName(filename, vm));
392  // Open output file
393  ofstream outfile;
394  if ((m_oneOutputFile && rank == 0) || !m_oneOutputFile)
395  {
396  outfile.open(outFile.c_str(), m_binary ? ios::binary : ios::out);
397  }
398 
399  if (m_oneOutputFile)
400  {
401  // Reduce on number of blocks and number of points.
402  m_f->m_comm->AllReduce(m_numBlocks, LibUtilities::ReduceSum);
403 
404  // Root process needs to know how much data everyone else has for
405  // writing in parallel.
406  m_rankFieldSizes = Array<OneD, int>(nprocs, 0);
407  m_rankConnSizes = Array<OneD, int>(nprocs, 0);
408  m_rankFieldSizes[rank] = m_fields[0].size();
409 
410  m_totConn = 0;
411  for (int i = 0; i < m_conn.size(); ++i)
412  {
413  m_totConn += m_conn[i].size();
414  }
415 
416  m_rankConnSizes[rank] = m_totConn;
417 
418  m_f->m_comm->AllReduce(m_rankFieldSizes, LibUtilities::ReduceSum);
419  m_f->m_comm->AllReduce(m_rankConnSizes, LibUtilities::ReduceSum);
420  }
421 
422  if (m_writeHeader)
423  {
424  WriteTecplotHeader(outfile, variables);
425  }
426 
427  // Write zone data.
428  WriteTecplotZone(outfile);
429 
430  // If we're a FE block format, write connectivity (m_conn will be empty for
431  // point data).
432  WriteTecplotConnectivity(outfile);
433 
434  if ((m_oneOutputFile && rank == 0) || !m_oneOutputFile)
435  {
436  cout << "Written file: " << GetFullOutName(filename, vm) << endl;
437  }
438 }
439 
440 /**
441  * @brief Write Tecplot files header
442  *
443  * @param outfile Output file name
444  * @param var Variables names
445  */
446 void OutputTecplot::WriteTecplotHeader(std::ofstream &outfile,
447  std::vector<std::string> &var)
448 {
449  outfile << "Variables = " << var[0];
450 
451  for (int i = 1; i < var.size(); ++i)
452  {
453  outfile << ", " << var[i];
454  }
455 
456  outfile << std::endl << std::endl;
457 }
458 
459 /**
460  * @brief Write Tecplot files header in binary format
461  *
462  * @param outfile Output file name
463  * @param var Variables names
464  */
465 void OutputTecplotBinary::WriteTecplotHeader(std::ofstream &outfile,
466  std::vector<std::string> &var)
467 {
468  if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
469  {
470  return;
471  }
472 
473  // Version number
474  outfile << "#!TDV112";
475 
476  // Int value of 1 for endian check
477  WriteStream(outfile, 1);
478 
479  // We'll probably write a full solution field
480  WriteStream(outfile, 0);
481 
482  // Title
483  std::string title = "";
484  WriteStream(outfile, title);
485 
486  // Number of variables
487  WriteStream(outfile, (int)var.size());
488 
489  for (int i = 0; i < var.size(); ++i)
490  {
491  WriteStream(outfile, var[i]);
492  }
493 }
494 
495 /**
496  * Write Tecplot zone output in ASCII
497  *
498  * @param outfile Output file name.
499  * @param expansion Expansion that is considered
500  */
501 void OutputTecplot::WriteTecplotZone(std::ofstream &outfile)
502 {
503  bool useDoubles = m_config["double"].as<bool>();
504 
505  if (useDoubles)
506  {
507  int precision = std::numeric_limits<double>::max_digits10;
508  outfile << std::setprecision(precision);
509  }
510 
511  // Write either points or finite element block
512  if (m_zoneType != eOrdered)
513  {
514  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) ||
516  {
517  // Number of points in zone
518  int nPoints = m_oneOutputFile ? Vmath::Vsum(m_f->m_comm->GetSize(),
519  m_rankFieldSizes, 1)
520  : m_fields[0].size();
521 
522  outfile << "Zone, N=" << nPoints << ", E=" << m_numBlocks
523  << ", F=FEBlock, ET=" << TecplotZoneTypeMap[m_zoneType]
524  << std::endl;
525  }
526 
527  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
528  {
529  for (int j = 0; j < m_fields.size(); ++j)
530  {
531  for (int i = 0; i < m_fields[j].size(); ++i)
532  {
533  if ((!(i % 1000)) && i)
534  {
535  outfile << std::endl;
536  }
537  outfile << m_fields[j][i] << " ";
538  }
539 
540  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
541  {
542  if (m_rankFieldSizes[n])
543  {
545  m_f->m_comm->Recv(n, tmp);
546 
547  for (int i = 0; i < m_rankFieldSizes[n]; ++i)
548  {
549  if ((!(i % 1000)) && i)
550  {
551  outfile << std::endl;
552  }
553  outfile << tmp[i] << " ";
554  }
555  }
556  }
557  outfile << std::endl;
558  }
559  }
560  else if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
561  {
562  if (m_fields[0].size())
563  {
564  for (int i = 0; i < m_fields.size(); ++i)
565  {
566  m_f->m_comm->Send(0, m_fields[i]);
567  }
568  }
569  }
570  else
571  {
572  // Write out coordinates and field data: ordered by field
573  // and then its data.
574  for (int j = 0; j < m_fields.size(); ++j)
575  {
576  for (int i = 0; i < m_fields[j].size(); ++i)
577  {
578  if ((!(i % 1000)) && i)
579  {
580  outfile << std::endl;
581  }
582  outfile << m_fields[j][i] << " ";
583  }
584  outfile << std::endl;
585  }
586  }
587  }
588  else
589  {
590  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) ||
592  {
593  std::string dirs[] = {"I", "J", "K"};
594  outfile << "Zone";
595  for (int i = 0; i < m_numPoints.size(); ++i)
596  {
597  outfile << ", " << dirs[i] << "=" << m_numPoints[i];
598  }
599  outfile << ", F=POINT" << std::endl;
600  }
601 
602  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
603  {
604  Array<OneD, NekDouble> tmp(m_fields.size());
605  for (int i = 0; i < m_fields[0].size(); ++i)
606  {
607  for (int j = 0; j < m_fields.size(); ++j)
608  {
609  outfile << setw(12) << m_fields[j][i] << " ";
610  }
611  outfile << std::endl;
612  }
613 
614  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
615  {
616  for (int i = 0; i < m_rankFieldSizes[n]; ++i)
617  {
618  m_f->m_comm->Recv(n, tmp);
619  for (int j = 0; j < m_fields.size(); ++j)
620  {
621  outfile << setw(12) << tmp[j] << " ";
622  }
623  outfile << std::endl;
624  }
625  }
626  }
627  else if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
628  {
629  Array<OneD, NekDouble> tmp(m_fields.size());
630  for (int i = 0; i < m_fields[0].size(); ++i)
631  {
632  for (int j = 0; j < m_fields.size(); ++j)
633  {
634  tmp[j] = m_fields[j][i];
635  }
636  m_f->m_comm->Send(0, tmp);
637  }
638  }
639  else
640  {
641  // Write out coordinates and field data: ordered by each
642  // point then each field.
643  for (int i = 0; i < m_fields[0].size(); ++i)
644  {
645  for (int j = 0; j < m_fields.size(); ++j)
646  {
647  outfile << setw(12) << m_fields[j][i] << " ";
648  }
649  outfile << std::endl;
650  }
651  }
652  }
653 }
654 
655 /**
656  * @brief Write either double-precision or single-precision output of field
657  * data.
658  *
659  * @param outfile Output file name.
660  * @param expansion Expansion that is considered
661  */
662 void OutputTecplotBinary::WriteDoubleOrFloat(std::ofstream &outfile,
664 {
665  // Data format: either double or single depending on user options
666  bool useDoubles = m_config["double"].as<bool>();
667 
668  if (useDoubles)
669  {
670  // For doubles, we can just write data.
671  WriteStream(outfile, data);
672  }
673  else
674  {
675  // For single precision, needs typecast first.
676  int nPts = data.size();
677  std::vector<float> tmp(data.size());
678  std::copy(&data[0], &data[0] + nPts, &tmp[0]);
679  WriteStream(outfile, tmp);
680  }
681 }
682 
683 /**
684  * Write Tecplot zone output in binary
685  *
686  * @param outfile Output file name.
687  * @param expansion Expansion that is considered
688  */
689 void OutputTecplotBinary::WriteTecplotZone(std::ofstream &outfile)
690 {
691  Array<OneD, NekDouble> fieldMin(m_fields.size());
692  Array<OneD, NekDouble> fieldMax(m_fields.size());
693 
694  // Data format: either double or single depending on user options
695  bool useDoubles = m_config["double"].as<bool>();
696 
697  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) || !m_oneOutputFile)
698  {
699  // Don't bother naming zone
700  WriteStream(outfile, 299.0f); // Zone marker
701 
702  // Write same name as preplot
703  int rank = m_f->m_comm->GetRank();
704  string zonename = "ZONE " + boost::lexical_cast<string>(rank);
705  WriteStream(outfile, zonename);
706 
707  WriteStream(outfile, -1); // No parent zone
708  WriteStream(outfile, -1); // No strand ID
709  WriteStream(outfile, 0.0); // Solution time
710  WriteStream(outfile, -1); // Unused, set to -1
711 
712  // Zone type: 1 = lineseg, 3 = quad, 5 = brick
713  WriteStream(outfile, (int)m_zoneType);
714 
715  WriteStream(outfile, 0); // Data at nodes
716  WriteStream(outfile, 0); // No 1-1 face neighbours
717  WriteStream(outfile, 0); // No user-defined connections
718 
719  if (m_zoneType == eOrdered)
720  {
721  for (int i = 0; i < m_numPoints.size(); ++i)
722  {
723  WriteStream(outfile, m_numPoints[i]);
724  }
725 
726  for (int i = m_numPoints.size(); i < 3; ++i)
727  {
728  WriteStream(outfile, 0);
729  }
730  }
731  else
732  {
733  // Number of points in zone
734  int nPoints = m_oneOutputFile ? Vmath::Vsum(m_f->m_comm->GetSize(),
735  m_rankFieldSizes, 1)
736  : m_fields[0].size();
737 
738  WriteStream(outfile, nPoints); // Total number of points
739  WriteStream(outfile, m_numBlocks); // Number of blocks
740  WriteStream(outfile, 0); // Unused
741  WriteStream(outfile, 0); // Unused
742  WriteStream(outfile, 0); // Unused
743  }
744 
745  WriteStream(outfile, 0); // No auxiliary data names
746 
747  // Finalise header
748  WriteStream(outfile, 357.0f);
749 
750  // Now start to write data section so that we can dump geometry
751  // information
752 
753  // Data marker
754  WriteStream(outfile, 299.0f);
755 
756  for (int j = 0; j < m_fields.size(); ++j)
757  {
758  WriteStream(outfile, useDoubles ? 2 : 1);
759  }
760 
761  // No passive variables or variable sharing, no zone connectivity
762  // sharing (we only dump one zone)
763  WriteStream(outfile, 0);
764  WriteStream(outfile, 0);
765  WriteStream(outfile, -1);
766  }
767 
768  for (int i = 0; i < m_fields.size(); ++i)
769  {
770  fieldMin[i] = Vmath::Vmin(m_fields[i].size(), m_fields[i], 1);
771  fieldMax[i] = Vmath::Vmax(m_fields[i].size(), m_fields[i], 1);
772  }
773 
774  m_f->m_comm->AllReduce(fieldMin, LibUtilities::ReduceMin);
775  m_f->m_comm->AllReduce(fieldMax, LibUtilities::ReduceMax);
776 
777  // Write out min/max of field data
778  if ((m_oneOutputFile && m_f->m_comm->GetRank() == 0) || !m_oneOutputFile)
779  {
780  for (int i = 0; i < m_fields.size(); ++i)
781  {
782  WriteStream(outfile, fieldMin[i]);
783  WriteStream(outfile, fieldMax[i]);
784  }
785  }
786 
787  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
788  {
789  for (int i = 0; i < m_fields.size(); ++i)
790  {
791  WriteDoubleOrFloat(outfile, m_fields[i]);
792 
793  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
794  {
796  m_f->m_comm->Recv(n, tmp);
797  WriteDoubleOrFloat(outfile, tmp);
798  }
799  }
800  }
801  else if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
802  {
803  for (int i = 0; i < m_fields.size(); ++i)
804  {
805  m_f->m_comm->Send(0, m_fields[i]);
806  }
807  }
808  else
809  {
810  for (int i = 0; i < m_fields.size(); ++i)
811  {
812  WriteDoubleOrFloat(outfile, m_fields[i]);
813  }
814  }
815 }
816 
817 /**
818  * @brief Write Tecplot connectivity information (ASCII)
819  *
820  * @param outfile Output file
821  */
822 void OutputTecplot::WriteTecplotConnectivity(std::ofstream &outfile)
823 {
824  // Ordered data have no connectivity information.
825  if (m_zoneType == eOrdered)
826  {
827  return;
828  }
829 
830  if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
831  {
832  // Need to amalgamate connectivity information
833  if (m_totConn)
834  {
836  for (int i = 0, cnt = 0; i < m_conn.size(); ++i)
837  {
838  if (m_conn[i].size())
839  {
840  Vmath::Vcopy(m_conn[i].size(), &m_conn[i][0], 1, &conn[cnt],
841  1);
842  cnt += m_conn[i].size();
843  }
844  }
845  m_f->m_comm->Send(0, conn);
846  }
847  }
848  else
849  {
850  int cnt = 1;
851  for (int i = 0; i < m_conn.size(); ++i)
852  {
853  const int nConn = m_conn[i].size();
854  for (int j = 0; j < nConn; ++j, ++cnt)
855  {
856  outfile << m_conn[i][j] + 1 << " ";
857  if (!(cnt % 1000))
858  {
859  outfile << std::endl;
860  }
861  }
862  }
863  outfile << endl;
864 
865  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
866  {
867  int offset = m_rankFieldSizes[0];
868 
869  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
870  {
871  if (m_rankConnSizes[n])
872  {
874  m_f->m_comm->Recv(n, conn);
875  for (int j = 0; j < conn.size(); ++j)
876  {
877  outfile << conn[j] + offset + 1 << " ";
878  if ((!(j % 1000)) && j)
879  {
880  outfile << std::endl;
881  }
882  }
883  }
884  offset += m_rankFieldSizes[n];
885  }
886  }
887  }
888 }
889 
891 {
892  if (m_oneOutputFile && m_f->m_comm->GetRank() > 0)
893  {
894  // Need to amalgamate connectivity information
896  for (int i = 0, cnt = 0; i < m_conn.size(); ++i)
897  {
898  Vmath::Vcopy(m_conn[i].size(), &m_conn[i][0], 1, &conn[cnt], 1);
899  cnt += m_conn[i].size();
900  }
901  m_f->m_comm->Send(0, conn);
902  }
903  else
904  {
905  for (int i = 0; i < m_conn.size(); ++i)
906  {
907  WriteStream(outfile, m_conn[i]);
908  }
909 
910  if (m_oneOutputFile && m_f->m_comm->GetRank() == 0)
911  {
912  int offset = m_rankFieldSizes[0];
913 
914  for (int n = 1; n < m_f->m_comm->GetSize(); ++n)
915  {
917  m_f->m_comm->Recv(n, conn);
918 
919  for (int j = 0; j < conn.size(); ++j)
920  {
921  conn[j] += offset;
922  }
923 
924  WriteStream(outfile, conn);
925  offset += m_rankFieldSizes[n];
926  }
927  }
928  }
929 }
930 
931 /**
932  * @brief Calculate number of Tecplot blocks.
933  *
934  * @param outfile Output file
935  */
937 {
938  int returnval = 0;
939 
940  if (m_f->m_exp[0]->GetExp(0)->GetNumBases() == 1)
941  {
942  for (int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
943  {
944  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0) - 1);
945  }
946  }
947  else if (m_f->m_exp[0]->GetExp(0)->GetNumBases() == 2)
948  {
949  for (int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
950  {
951  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0) - 1) *
952  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1) - 1);
953  }
954  }
955  else
956  {
957  for (int i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
958  {
959  returnval += (m_f->m_exp[0]->GetExp(i)->GetNumPoints(0) - 1) *
960  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(1) - 1) *
961  (m_f->m_exp[0]->GetExp(i)->GetNumPoints(2) - 1);
962  }
963  }
964 
965  return returnval;
966 }
967 
968 /**
969  * @brief Calculate connectivity information for each expansion dimension.
970  *
971  * @param outfile Output file
972  */
974 {
975  int i, j, k, l;
976  int nbase = m_f->m_exp[0]->GetExp(0)->GetNumBases();
977  int cnt = 0;
978 
979  m_conn.resize(m_f->m_exp[0]->GetNumElmts());
980 
981  for (i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
982  {
983  cnt = m_f->m_exp[0]->GetPhys_Offset(i);
984 
985  if (nbase == 1)
986  {
987  int cnt2 = 0;
988  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
989  int nPlanes = 1;
990 
991  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e2DH1D)
992  {
993  nPlanes = m_f->m_exp[0]->GetZIDs().size();
994 
995  if (nPlanes > 1)
996  {
997  int totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
998 
999  Array<OneD, int> conn(4 * (np0 - 1) * (nPlanes - 1));
1000  for (int n = 1; n < nPlanes; ++n)
1001  {
1002  for (k = 1; k < np0; ++k)
1003  {
1004  conn[cnt2++] = cnt + (n - 1) * totPoints + k;
1005  conn[cnt2++] = cnt + (n - 1) * totPoints + k - 1;
1006  conn[cnt2++] = cnt + n * totPoints + k - 1;
1007  conn[cnt2++] = cnt + n * totPoints + k;
1008  }
1009  }
1010  m_conn[i] = conn;
1011  }
1012  }
1013 
1014  if (nPlanes == 1)
1015  {
1016  Array<OneD, int> conn(2 * (np0 - 1));
1017 
1018  for (k = 1; k < np0; ++k)
1019  {
1020  conn[cnt2++] = cnt + k;
1021  conn[cnt2++] = cnt + k - 1;
1022  }
1023 
1024  m_conn[i] = conn;
1025  }
1026  }
1027  else if (nbase == 2)
1028  {
1029  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
1030  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
1031  int totPoints = m_f->m_exp[0]->GetTotPoints();
1032  int nPlanes = 1;
1033  int cnt2 = 0;
1034 
1035  if (m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D)
1036  {
1037  nPlanes = m_f->m_exp[0]->GetZIDs().size();
1038 
1039  // default to 2D case for HalfMode when nPlanes = 1
1040  if (nPlanes > 1)
1041  {
1042  // If Fourier points, output extra plane to fill domain
1043  nPlanes += 1;
1044  totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
1045 
1046  Array<OneD, int> conn(8 * (np1 - 1) * (np0 - 1) *
1047  (nPlanes - 1));
1048 
1049  for (int n = 1; n < nPlanes; ++n)
1050  {
1051  for (j = 1; j < np1; ++j)
1052  {
1053  for (k = 1; k < np0; ++k)
1054  {
1055  conn[cnt2++] = cnt + (n - 1) * totPoints +
1056  (j - 1) * np0 + k - 1;
1057  conn[cnt2++] = cnt + (n - 1) * totPoints +
1058  (j - 1) * np0 + k;
1059  conn[cnt2++] =
1060  cnt + (n - 1) * totPoints + j * np0 + k;
1061  conn[cnt2++] =
1062  cnt + (n - 1) * totPoints + j * np0 + k - 1;
1063  conn[cnt2++] =
1064  cnt + n * totPoints + (j - 1) * np0 + k - 1;
1065  conn[cnt2++] =
1066  cnt + n * totPoints + (j - 1) * np0 + k;
1067  conn[cnt2++] =
1068  cnt + n * totPoints + j * np0 + k;
1069  conn[cnt2++] =
1070  cnt + n * totPoints + j * np0 + k - 1;
1071  }
1072  }
1073  }
1074  m_conn[i] = conn;
1075  }
1076  }
1077 
1078  if (nPlanes == 1)
1079  {
1080  Array<OneD, int> conn(4 * (np0 - 1) * (np1 - 1));
1081  for (j = 1; j < np1; ++j)
1082  {
1083  for (k = 1; k < np0; ++k)
1084  {
1085  conn[cnt2++] = cnt + (j - 1) * np0 + k - 1;
1086  conn[cnt2++] = cnt + (j - 1) * np0 + k;
1087  conn[cnt2++] = cnt + j * np0 + k;
1088  conn[cnt2++] = cnt + j * np0 + k - 1;
1089  }
1090  }
1091  m_conn[i] = conn;
1092  }
1093  }
1094  else if (nbase == 3)
1095  {
1096  int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
1097  int np1 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(1);
1098  int np2 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(2);
1099  int cnt2 = 0;
1100 
1101  Array<OneD, int> conn(8 * (np0 - 1) * (np1 - 1) * (np2 - 1));
1102 
1103  for (j = 1; j < np2; ++j)
1104  {
1105  for (k = 1; k < np1; ++k)
1106  {
1107  for (l = 1; l < np0; ++l)
1108  {
1109  conn[cnt2++] =
1110  cnt + (j - 1) * np0 * np1 + (k - 1) * np0 + l - 1;
1111  conn[cnt2++] =
1112  cnt + (j - 1) * np0 * np1 + (k - 1) * np0 + l;
1113  conn[cnt2++] = cnt + (j - 1) * np0 * np1 + k * np0 + l;
1114  conn[cnt2++] =
1115  cnt + (j - 1) * np0 * np1 + k * np0 + l - 1;
1116  conn[cnt2++] =
1117  cnt + j * np0 * np1 + (k - 1) * np0 + l - 1;
1118  conn[cnt2++] = cnt + j * np0 * np1 + (k - 1) * np0 + l;
1119  conn[cnt2++] = cnt + j * np0 * np1 + k * np0 + l;
1120  conn[cnt2++] = cnt + j * np0 * np1 + k * np0 + l - 1;
1121  }
1122  }
1123  }
1124 
1125  m_conn[i] = conn;
1126  }
1127  else
1128  {
1129  ASSERTL0(false, "Not set up for this dimension");
1130  }
1131  }
1132 }
1133 } // namespace FieldUtils
1134 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
FieldSharedPtr m_f
Field object.
Definition: Module.h:225
std::map< std::string, ConfigOption > m_config
List of configuration values.
Definition: Module.h:228
Converter from fld to vtk.
virtual void Process(po::variables_map &vm)
Write fld to output file.
void WriteDoubleOrFloat(std::ofstream &outfile, Array< OneD, NekDouble > &data)
Write either double-precision or single-precision output of field data.
virtual void WriteTecplotZone(std::ofstream &outfile)
virtual void WriteTecplotHeader(std::ofstream &outfile, std::vector< std::string > &var)
Write Tecplot files header in binary format.
virtual void WriteTecplotConnectivity(std::ofstream &outfile)
Write Tecplot connectivity information (ASCII)
int m_totConn
Total number of connectivity entries.
std::vector< Array< OneD, int > > m_conn
Connectivty for each block: one per element.
Array< OneD, int > m_rankConnSizes
Each rank's connectivity sizes.
bool m_binary
True if writing binary field output.
Definition: OutputTecplot.h:97
bool m_oneOutputFile
True if writing a single output file.
Definition: OutputTecplot.h:99
virtual fs::path GetFullOutName(std::string &filename, po::variables_map &vm)
virtual void WriteTecplotZone(std::ofstream &outfile)
void CalculateConnectivity()
Calculate connectivity information for each expansion dimension.
virtual void WriteTecplotConnectivity(std::ofstream &outfile)
Write Tecplot connectivity information (ASCII)
virtual void Process(po::variables_map &vm)
Write fld to output file.
bool m_writeHeader
True if writing header.
virtual void OutputFromExp(po::variables_map &vm)
Write from m_exp to output file.
TecplotZoneType m_zoneType
Tecplot zone type of output.
Array< OneD, int > m_rankFieldSizes
Each rank's field sizes.
std::vector< int > m_numPoints
Number of points per block in Tecplot file.
virtual fs::path GetPath(std::string &filename, po::variables_map &vm)
virtual void WriteTecplotHeader(std::ofstream &outfile, std::vector< std::string > &var)
Write Tecplot files header.
void WriteTecplotFile(po::variables_map &vm)
virtual void OutputFromPts(po::variables_map &vm)
Write from pts to output file.
int m_numBlocks
Number of blocks in Tecplot file.
int m_coordim
Coordinate dimension of output.
int GetNumTecplotBlocks()
Calculate number of Tecplot blocks.
Array< OneD, Array< OneD, NekDouble > > m_fields
Field data to output.
virtual void OutputFromData(po::variables_map &vm)
Write from data to output file.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
def copy(self)
Definition: pycml.py:2663
void WriteStream(std::ostream &outfile, T data)
Helper function to write binary data to stream.
std::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:989
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:285
std::string TecplotZoneTypeMap[]
ModuleFactory & GetModuleFactory()
Definition: Module.cpp:49
std::string PortablePath(const boost::filesystem::path &path)
create portable path on different platforms for boost::filesystem path
Definition: FileSystem.cpp:41
std::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:190
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::int32_t NekInt32
double NekDouble
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:1050
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384
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:945
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
Represents a command-line configuration option.
Definition: Module.h:131