Nektar++
OutputVtkBase.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: OutputVtkBase.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: VTK file format output base class.
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 #include <boost/format.hpp>
42 
44 
45 #include "OutputVtkBase.h"
46 
47 namespace Nektar
48 {
49 namespace FieldUtils
50 {
51 
52 // Disable the base VTK factory if using the VTK library, this is so we can
53 // register the same extension for both.
54 #if !NEKTAR_USING_VTK
55 ModuleKey OutputVtkBase::m_className =
57  OutputVtkBase::create,
58  "Writes a VTU file.");
59 #endif
60 
61 OutputVtkBase::OutputVtkBase(FieldSharedPtr f) : OutputFileBase(f)
62 {
63  m_requireEquiSpaced = true;
64 }
65 
67 {
68 }
69 
70 void OutputVtkBase::OutputFromPts(po::variables_map &vm)
71 {
72  int i, j;
73  LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
74 
75  // Extract the output filename and extension
76  string filename = PrepareOutput(vm);
77 
78  // Write solution.
79  ofstream outfile(filename.c_str());
80  WriteVtkHeader(outfile);
81  int nfields = 0;
82  int dim = fPts->GetDim();
83 
84  int nvert = 1;
85  int vtktype = 1;
86  switch (fPts->GetPtsType())
87  {
90  {
92  "VTK output needs setting up for ePtsFile and ePtsLine");
93  break;
94  }
96  {
98  "VTK output needs setting up for PtsPlane");
99  break;
100  }
102  {
104  "VTK output needs setting up for PtsBox");
105  break;
106  }
108  {
109  nvert = 2;
110  vtktype = 3;
111  break;
112  }
114  {
115  nvert = 3;
116  vtktype = 5;
117  break;
118  }
120  {
121  nvert = 4;
122  vtktype = 10;
123  break;
124  }
125  default:
126  NEKERROR(ErrorUtil::efatal, "ptsType not supported yet.");
127  }
128 
129  vector<Array<OneD, int>> ptsConn;
130  fPts->GetConnectivity(ptsConn);
131 
132  nfields = fPts->GetNFields();
133 
134  int nPts = fPts->GetNpoints();
135  int numBlocks = 0;
136  for (i = 0; i < ptsConn.size(); ++i)
137  {
138  numBlocks += ptsConn[i].size() / nvert;
139  }
140 
141  // write out pieces of data.
142  outfile << " <Piece NumberOfPoints=\"" << nPts << "\" NumberOfCells=\""
143  << numBlocks << "\">" << endl;
144  outfile << " <Points>" << endl;
145  outfile << " <DataArray type=\"Float64\" "
146  << "NumberOfComponents=\"" << 3 << "\" format=\"ascii\">" << endl;
147  for (i = 0; i < nPts; ++i)
148  {
149  for (j = 0; j < dim; ++j)
150  {
151  outfile << " " << setprecision(8) << scientific
152  << fPts->GetPointVal(j, i) << " ";
153  }
154  for (j = dim; j < 3; ++j)
155  {
156  // pack to 3D since paraview does not seem to handle 2D
157  outfile << " 0.000000";
158  }
159  outfile << endl;
160  }
161  outfile << " </DataArray>" << endl;
162  outfile << " </Points>" << endl;
163  outfile << " <Cells>" << endl;
164  outfile << " <DataArray type=\"Int32\" "
165  << "Name=\"connectivity\" format=\"ascii\">" << endl;
166 
167  // dump connectivity data if it exists
168  outfile << " ";
169  int cnt = 1;
170  for (i = 0; i < ptsConn.size(); ++i)
171  {
172  for (j = 0; j < ptsConn[i].size(); ++j)
173  {
174  outfile << ptsConn[i][j] << " ";
175  if ((!(cnt % nvert)) && cnt)
176  {
177  outfile << std::endl;
178  outfile << " ";
179  }
180  cnt++;
181  }
182  }
183  outfile << " </DataArray>" << endl;
184  outfile << " <DataArray type=\"Int32\" "
185  << "Name=\"offsets\" format=\"ascii\">" << endl;
186 
187  outfile << " ";
188  for (i = 0; i < numBlocks; ++i)
189  {
190  outfile << i * nvert + nvert << " ";
191  }
192  outfile << endl;
193  outfile << " </DataArray>" << endl;
194  outfile << " <DataArray type=\"UInt8\" "
195  << "Name=\"types\" format=\"ascii\">" << endl;
196  outfile << " ";
197  for (i = 0; i < numBlocks; ++i)
198  {
199  outfile << vtktype << " ";
200  }
201  outfile << endl;
202  outfile << " </DataArray>" << endl;
203  outfile << " </Cells>" << endl;
204  outfile << " <PointData>" << endl;
205 
206  // printing the fields
207  for (j = 0; j < nfields; ++j)
208  {
209  outfile << " <DataArray type=\"Float64\" Name=\""
210  << m_f->m_variables[j] << "\">" << endl;
211  outfile << " ";
212  for (i = 0; i < fPts->GetNpoints(); ++i)
213  {
214  outfile << fPts->GetPointVal(dim + j, i) << " ";
215  }
216  outfile << endl;
217  outfile << " </DataArray>" << endl;
218  }
219 
220  outfile << " </PointData>" << endl;
221  outfile << " </Piece>" << endl;
222 
223  WriteVtkFooter(outfile);
224  cout << "Written file: " << filename << endl;
225 
226  // output parallel outline info if necessary
227  if ((m_f->m_comm->GetRank() == 0) && (m_f->m_comm->GetSize() != 1))
228  {
229  WritePVtu(vm);
230  cout << "Written file: " << filename << endl;
231  }
232 }
233 
234 void OutputVtkBase::OutputFromExp(po::variables_map &vm)
235 {
236  int i, j;
237  // Extract the output filename and extension
238  string filename = PrepareOutput(vm);
239 
240  // Write solution.
241  ofstream outfile(filename.c_str());
242  WriteVtkHeader(outfile);
243  int nfields = m_f->m_variables.size();
244 
245  int nstrips;
246  m_f->m_session->LoadParameter("Strip_Z", nstrips, 1);
247 
248  // Homogeneous strip variant
249  for (int s = 0; s < nstrips; ++s)
250  {
251  // For each field write out field data for each expansion.
252  for (i = 0; i < m_f->m_exp[0]->GetNumElmts(); ++i)
253  {
254  m_f->m_exp[0]->WriteVtkPieceHeader(outfile, i, s);
255 
256  // For this expansion write out each field.
257  for (j = 0; j < nfields; ++j)
258  {
259  m_f->m_exp[s * nfields + j]->WriteVtkPieceData(
260  outfile, i, m_f->m_variables[j]);
261  }
262  m_f->m_exp[0]->WriteVtkPieceFooter(outfile, i);
263  }
264  }
265 
266  if (m_f->m_exp[0]->GetNumElmts() == 0)
267  {
268  WriteEmptyVtkPiece(outfile);
269  }
270 
271  WriteVtkFooter(outfile);
272  cout << "Written file: " << filename << endl;
273 
274  // output parallel outline info if necessary
275  if ((m_f->m_comm->GetRank() == 0) && (m_f->m_comm->GetSize() != 1))
276  {
277  WritePVtu(vm);
278  }
279 }
280 
281 void OutputVtkBase::OutputFromData(po::variables_map &vm)
282 {
283  boost::ignore_unused(vm);
284  NEKERROR(ErrorUtil::efatal, "OutputVtk can't write using only FieldData.");
285 }
286 
287 fs::path OutputVtkBase::GetPath(std::string &filename, po::variables_map &vm)
288 {
289  boost::ignore_unused(vm);
290 
291  int nprocs = m_f->m_comm->GetSize();
292  fs::path specPath;
293  if (nprocs == 1)
294  {
295  specPath = fs::path(filename);
296  }
297  else
298  {
299  // replace .vtu by _vtu
300  int dot = filename.find_last_of('.');
301  string path = filename.substr(0, dot) + "_vtu";
302  specPath = fs::path(path);
303  }
304  return fs::path(specPath);
305 }
306 
307 fs::path OutputVtkBase::GetFullOutName(std::string &filename,
308  po::variables_map &vm)
309 {
310  int nprocs = m_f->m_comm->GetSize();
311 
312  fs::path fulloutname;
313  if (nprocs == 1)
314  {
315  fulloutname = filename;
316  }
317  else
318  {
319  // Guess at filename that might belong to this process.
320  boost::format pad("P%1$07d.%2$s");
321  pad % m_f->m_comm->GetRank() % "vtu";
322 
323  // Generate full path name
324  fs::path specPath = GetPath(filename, vm);
325  fs::path poutfile(pad.str());
326  fulloutname = specPath / poutfile;
327  }
328  return fulloutname;
329 }
330 
331 void OutputVtkBase::WriteVtkHeader(std::ostream &outfile)
332 {
333  outfile << "<?xml version=\"1.0\"?>" << endl;
334  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
335  << "byte_order=\"LittleEndian\">" << endl;
336  outfile << " <UnstructuredGrid>" << endl;
337 }
338 
339 void OutputVtkBase::WriteVtkFooter(std::ostream &outfile)
340 {
341  outfile << " </UnstructuredGrid>" << endl;
342  outfile << "</VTKFile>" << endl;
343 }
344 
345 void OutputVtkBase::WriteEmptyVtkPiece(std::ofstream &outfile)
346 {
347  // write out empty piece of data.
348  outfile << " <Piece NumberOfPoints=\"" << 0 << "\" NumberOfCells=\"" << 0
349  << "\">" << endl;
350  outfile << " <Points>" << endl;
351  outfile << " <DataArray type=\"Float64\" "
352  << "NumberOfComponents=\"" << 3 << "\" format=\"ascii\">" << endl;
353  outfile << " </DataArray>" << endl;
354  outfile << " </Points>" << endl;
355  outfile << " <Cells>" << endl;
356  outfile << " <DataArray type=\"Int32\" "
357  << "Name=\"connectivity\" format=\"ascii\">" << endl;
358  outfile << " </DataArray>" << endl;
359  outfile << " <DataArray type=\"Int32\" "
360  << "Name=\"offsets\" format=\"ascii\">" << endl;
361 
362  outfile << " ";
363  outfile << endl;
364  outfile << " </DataArray>" << endl;
365  outfile << " <DataArray type=\"UInt8\" "
366  << "Name=\"types\" format=\"ascii\">" << endl;
367  outfile << " ";
368  outfile << endl;
369  outfile << " </DataArray>" << endl;
370  outfile << " </Cells>" << endl;
371  outfile << " <PointData>" << endl;
372 
373  outfile << " </PointData>" << endl;
374  outfile << " </Piece>" << endl;
375 }
376 
377 void OutputVtkBase::WritePVtu(po::variables_map &vm)
378 {
379  string filename = m_config["outfile"].as<string>();
380  int dot = filename.find_last_of('.');
381  string body = filename.substr(0, dot);
382  filename = body + ".pvtu";
383 
384  ofstream outfile(filename.c_str());
385 
386  int nprocs = m_f->m_comm->GetSize();
387  string path = LibUtilities::PortablePath(GetPath(filename, vm));
388 
389  outfile << "<?xml version=\"1.0\"?>" << endl;
390  outfile << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" "
391  << "byte_order=\"LittleEndian\">" << endl;
392  outfile << "<PUnstructuredGrid GhostLevel=\"0\">" << endl;
393  outfile << "<PPoints> " << endl;
394  outfile << "<PDataArray type=\"Float64\" NumberOfComponents=\"" << 3
395  << "\"/> " << endl;
396  outfile << "</PPoints>" << endl;
397  outfile << "<PCells>" << endl;
398  outfile << "<PDataArray type=\"Int32\" Name=\"connectivity\" "
399  "NumberOfComponents=\"1\"/>"
400  << endl;
401  outfile << "<PDataArray type=\"Int32\" Name=\"offsets\" "
402  "NumberOfComponents=\"1\"/>"
403  << endl;
404  outfile << "<PDataArray type=\"UInt8\" Name=\"types\" "
405  "NumberOfComponents=\"1\"/>"
406  << endl;
407  outfile << "</PCells>" << endl;
408  outfile << "<PPointData Scalars=\"Material\">" << endl;
409  for (int i = 0; i < m_f->m_variables.size(); ++i)
410  {
411  outfile << "<PDataArray type=\"Float64\" Name=\"" << m_f->m_variables[i]
412  << "\"/>" << endl;
413  }
414  outfile << "</PPointData>" << endl;
415 
416  for (int i = 0; i < nprocs; ++i)
417  {
418  boost::format pad("P%1$07d.vtu");
419  pad % i;
420  outfile << "<Piece Source=\"" << path << "/" << pad.str() << "\"/>"
421  << endl;
422  }
423  outfile << "</PUnstructuredGrid>" << endl;
424  outfile << "</VTKFile>" << endl;
425 
426  cout << "Written file: " << filename << endl;
427 }
428 
429 std::string OutputVtkBase::PrepareOutput(po::variables_map &vm)
430 {
431  // Extract the output filename and extension
432  string filename = m_config["outfile"].as<string>();
433 
434  fs::path specPath = GetPath(filename, vm);
435  fs::path fulloutname = GetFullOutName(filename, vm);
436  filename = LibUtilities::PortablePath(fulloutname);
437 
438  if (m_f->m_comm->GetSize() != 1)
439  {
440  if (m_f->m_comm->TreatAsRankZero())
441  {
442  try
443  {
444  fs::create_directory(specPath);
445  }
446  catch (fs::filesystem_error &e)
447  {
448  ASSERTL0(false, "Filesystem error: " + string(e.what()));
449  }
450  cout << "Writing files to directory: " << specPath << endl;
451  }
452  m_f->m_comm->Block();
453  }
454  else
455  {
456  cout << "Writing: " << specPath << endl;
457  }
458  return filename;
459 }
460 
461 } // namespace FieldUtils
462 } // 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.
std::string PrepareOutput(po::variables_map &vm)
void WriteEmptyVtkPiece(std::ofstream &outfile)
virtual fs::path GetPath(std::string &filename, po::variables_map &vm)
virtual void OutputFromPts(po::variables_map &vm)
Write from pts to output file.
virtual fs::path GetFullOutName(std::string &filename, po::variables_map &vm)
void WriteVtkHeader(std::ostream &outfile)
void WritePVtu(po::variables_map &vm)
virtual void OutputFromData(po::variables_map &vm)
Write from data to output file.
virtual void OutputFromExp(po::variables_map &vm)
Write from m_exp to output file.
void WriteVtkFooter(std::ostream &outfile)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
std::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:989
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:285
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