MFEM  v4.6.0
Finite element discretization library
datacollection.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "fem.hpp"
13 #include "../mesh/nurbs.hpp"
14 #include "../mesh/vtk.hpp"
15 #include "../general/binaryio.hpp"
16 #include "../general/text.hpp"
17 #include "picojson.h"
18 
19 #include <cerrno> // errno
20 #include <sstream>
21 #include <regex>
22 
23 #ifndef _WIN32
24 #include <sys/stat.h> // mkdir
25 #else
26 #include <direct.h> // _mkdir
27 #define mkdir(dir, mode) _mkdir(dir)
28 #endif
29 
30 namespace mfem
31 {
32 
33 // static method
34 int DataCollection::create_directory(const std::string &dir_name,
35  const Mesh *mesh, int myid)
36 {
37  // create directories recursively
38  const char path_delim = '/';
39  std::string::size_type pos = 0;
40  int err_flag;
41 #ifdef MFEM_USE_MPI
42  const ParMesh *pmesh = dynamic_cast<const ParMesh*>(mesh);
43 #endif
44 
45  do
46  {
47  pos = dir_name.find(path_delim, pos+1);
48  std::string subdir = dir_name.substr(0, pos);
49 
50 #ifndef MFEM_USE_MPI
51  err_flag = mkdir(subdir.c_str(), 0777);
52  err_flag = (err_flag && (errno != EEXIST)) ? 1 : 0;
53 #else
54  if (myid == 0 || pmesh == NULL)
55  {
56  err_flag = mkdir(subdir.c_str(), 0777);
57  err_flag = (err_flag && (errno != EEXIST)) ? 1 : 0;
58  }
59 #endif
60  }
61  while ( pos != std::string::npos );
62 
63 #ifdef MFEM_USE_MPI
64  if (pmesh)
65  {
66  MPI_Bcast(&err_flag, 1, MPI_INT, 0, pmesh->GetComm());
67  }
68 #endif
69 
70  return err_flag;
71 }
72 
73 // class DataCollection implementation
74 
75 DataCollection::DataCollection(const std::string& collection_name, Mesh *mesh_)
76 {
77  std::string::size_type pos = collection_name.find_last_of('/');
78  if (pos == std::string::npos)
79  {
80  name = collection_name;
81  // leave prefix_path empty
82  }
83  else
84  {
85  prefix_path = collection_name.substr(0, pos+1);
86  name = collection_name.substr(pos+1);
87  }
88  mesh = mesh_;
89  myid = 0;
90  num_procs = 1;
91  serial = true;
92  appendRankToFileName = false;
93 
94 #ifdef MFEM_USE_MPI
95  m_comm = MPI_COMM_NULL;
96  ParMesh *par_mesh = dynamic_cast<ParMesh*>(mesh);
97  if (par_mesh)
98  {
99  myid = par_mesh->GetMyRank();
100  num_procs = par_mesh->GetNRanks();
101  m_comm = par_mesh->GetComm();
102  serial = false;
103  appendRankToFileName = true;
104  }
105 #endif
106  own_data = false;
107  cycle = -1;
108  time = 0.0;
109  time_step = 0.0;
112  format = SERIAL_FORMAT; // use serial mesh format
113  compression = 0;
114  error = No_Error;
115 }
116 
118 {
119  if (own_data && new_mesh != mesh) { delete mesh; }
120  mesh = new_mesh;
121  myid = 0;
122  num_procs = 1;
123  serial = true;
124  appendRankToFileName = false;
125 
126 #ifdef MFEM_USE_MPI
127  m_comm = MPI_COMM_NULL;
128  ParMesh *par_mesh = dynamic_cast<ParMesh*>(mesh);
129  if (par_mesh)
130  {
131  myid = par_mesh->GetMyRank();
132  num_procs = par_mesh->GetNRanks();
133  m_comm = par_mesh->GetComm();
134  serial = false;
135  appendRankToFileName = true;
136  }
137 #endif
138 }
139 
140 #ifdef MFEM_USE_MPI
141 void DataCollection::SetMesh(MPI_Comm comm, Mesh *new_mesh)
142 {
143  // This seems to be the cleanest way to accomplish this
144  // and avoid duplicating fine grained details:
145 
146  SetMesh(new_mesh);
147 
148  m_comm = comm;
149  MPI_Comm_rank(comm, &myid);
150  MPI_Comm_size(comm, &num_procs);
151 }
152 #endif
153 
155 {
156  switch (fmt)
157  {
158  case SERIAL_FORMAT: break;
159 #ifdef MFEM_USE_MPI
160  case PARALLEL_FORMAT: break;
161 #endif
162  default: MFEM_ABORT("unknown format: " << fmt);
163  }
164  format = fmt;
165 }
166 
168 {
169  compression = comp;
170 #ifndef MFEM_USE_ZLIB
171  MFEM_VERIFY(!compression, "ZLib not enabled in MFEM build.");
172 #endif
173 }
174 
175 void DataCollection::SetPrefixPath(const std::string& prefix)
176 {
177  if (!prefix.empty())
178  {
179  prefix_path = prefix;
180  if (!prefix_path.empty() && prefix_path[prefix_path.size()-1] != '/')
181  {
182  prefix_path += '/';
183  }
184  }
185  else
186  {
187  prefix_path.clear();
188  }
189 }
190 
191 void DataCollection::Load(int cycle_)
192 {
193  MFEM_ABORT("this method is not implemented");
194 }
195 
197 {
198  SaveMesh();
199 
200  if (error) { return; }
201 
202  for (FieldMapIterator it = field_map.begin(); it != field_map.end(); ++it)
203  {
204  SaveOneField(it);
205  // Even if there is an error, try saving the other fields
206  }
207 
208  for (QFieldMapIterator it = q_field_map.begin(); it != q_field_map.end();
209  ++it)
210  {
211  SaveOneQField(it);
212  }
213 }
214 
216 {
217  std::string dir_name = prefix_path + name;
218  if (cycle != -1)
219  {
220  dir_name += "_" + to_padded_string(cycle, pad_digits_cycle);
221  }
222  int error_code = create_directory(dir_name, mesh, myid);
223  if (error_code)
224  {
225  error = WRITE_ERROR;
226  MFEM_WARNING("Error creating directory: " << dir_name);
227  return; // do not even try to write the mesh
228  }
229 
230  std::string mesh_name = GetMeshFileName();
231  mfem::ofgzstream mesh_file(mesh_name, compression);
232  mesh_file.precision(precision);
233 #ifdef MFEM_USE_MPI
234  const ParMesh *pmesh = dynamic_cast<const ParMesh*>(mesh);
235  if (pmesh && format == PARALLEL_FORMAT)
236  {
237  pmesh->ParPrint(mesh_file);
238  }
239  else
240 #endif
241  {
242  mesh->Print(mesh_file);
243  }
244  if (!mesh_file)
245  {
246  error = WRITE_ERROR;
247  MFEM_WARNING("Error writing mesh to file: " << mesh_name);
248  }
249 }
250 
252 {
253  return (serial || format == SERIAL_FORMAT) ? "mesh" : "pmesh";
254 }
255 
257 {
259 }
260 
261 std::string DataCollection::GetFieldFileName(const std::string &field_name)
262 const
263 {
264  std::string dir_name = prefix_path + name;
265  if (cycle != -1)
266  {
267  dir_name += "_" + to_padded_string(cycle, pad_digits_cycle);
268  }
269  std::string file_name = dir_name + "/" + field_name;
271  {
272  file_name += "." + to_padded_string(myid, pad_digits_rank);
273  }
274  return file_name;
275 }
276 
278 {
279  mfem::ofgzstream field_file(GetFieldFileName(it->first), compression);
280 
281  field_file.precision(precision);
282  (it->second)->Save(field_file);
283  if (!field_file)
284  {
285  error = WRITE_ERROR;
286  MFEM_WARNING("Error writing field to file: " << it->first);
287  }
288 }
289 
291 {
292  mfem::ofgzstream q_field_file(GetFieldFileName(it->first), compression);
293 
294  q_field_file.precision(precision);
295  (it->second)->Save(q_field_file);
296  if (!q_field_file)
297  {
298  error = WRITE_ERROR;
299  MFEM_WARNING("Error writing q-field to file: " << it->first);
300  }
301 }
302 
303 void DataCollection::SaveField(const std::string &field_name)
304 {
305  FieldMapIterator it = field_map.find(field_name);
306  if (it != field_map.end())
307  {
308  SaveOneField(it);
309  }
310 }
311 
312 void DataCollection::SaveQField(const std::string &q_field_name)
313 {
314  QFieldMapIterator it = q_field_map.find(q_field_name);
315  if (it != q_field_map.end())
316  {
317  SaveOneQField(it);
318  }
319 }
320 
322 {
323  if (own_data) { delete mesh; }
324  mesh = NULL;
325 
328  own_data = false;
329 }
330 
332 {
333  DeleteData();
334  field_map.clear();
335  q_field_map.clear();
336 }
337 
339 {
340  DeleteData();
341 }
342 
343 
344 // class VisItDataCollection implementation
345 
347 {
348  if (mesh)
349  {
351  topo_dim = mesh->Dimension();
352  if (mesh->NURBSext)
353  {
356  }
357  }
358  else
359  {
360  spatial_dim = 0;
361  topo_dim = 0;
362  }
363 }
364 
365 VisItDataCollection::VisItDataCollection(const std::string& collection_name,
366  Mesh *mesh)
367  : DataCollection(collection_name, mesh)
368 {
369  appendRankToFileName = true; // always include rank in file names
370  cycle = 0; // always include cycle in directory names
371 
374 
375  UpdateMeshInfo();
376 }
377 
378 #ifdef MFEM_USE_MPI
380  const std::string& collection_name,
381  Mesh *mesh)
382  : DataCollection(collection_name, mesh)
383 {
384  m_comm = comm;
385  MPI_Comm_rank(comm, &myid);
386  MPI_Comm_size(comm, &num_procs);
387  appendRankToFileName = true; // always include rank in file names
388  cycle = 0; // always include cycle in directory names
389 
392 
393  UpdateMeshInfo();
394 }
395 #endif
396 
398 {
399  DataCollection::SetMesh(new_mesh);
400  appendRankToFileName = true;
401  UpdateMeshInfo();
402 }
403 
404 #ifdef MFEM_USE_MPI
405 void VisItDataCollection::SetMesh(MPI_Comm comm, Mesh *new_mesh)
406 {
407  // use VisItDataCollection's custom SetMesh, then set MPI info
408  SetMesh(new_mesh);
409  m_comm = comm;
410  MPI_Comm_rank(comm, &myid);
411  MPI_Comm_size(comm, &num_procs);
412 }
413 #endif
414 
415 void VisItDataCollection::RegisterField(const std::string& name,
416  GridFunction *gf)
417 {
418  int LOD = 1;
419  if (gf->FESpace()->GetNURBSext())
420  {
421  LOD = gf->FESpace()->GetNURBSext()->GetOrder();
422  }
423  else
424  {
425  for (int e=0; e<gf->FESpace()->GetNE(); e++)
426  {
427  LOD = std::max(LOD,gf->FESpace()->GetFE(e)->GetOrder());
428  }
429  }
430 
432  field_info_map[name] = VisItFieldInfo("nodes", gf->VectorDim(), LOD);
434 }
435 
436 void VisItDataCollection::RegisterQField(const std::string& name,
437  QuadratureFunction *qf)
438 {
439  int LOD = -1;
440  Mesh *mesh = qf->GetSpace()->GetMesh();
441  for (int e=0; e<qf->GetSpace()->GetNE(); e++)
442  {
445  qf->GetIntRule(e).GetNPoints());
446 
447  LOD = std::max(LOD,locLOD);
448  }
449 
451  field_info_map[name] = VisItFieldInfo("elements", 1, LOD);
453 }
454 
455 void VisItDataCollection::SetLevelsOfDetail(int levels_of_detail)
456 {
457  visit_levels_of_detail = levels_of_detail;
458 }
459 
460 void VisItDataCollection::SetMaxLevelsOfDetail(int max_levels_of_detail)
461 {
462  visit_max_levels_of_detail = max_levels_of_detail;
463 }
464 
466 {
467  field_info_map.clear();
469 }
470 
472 {
474  SaveRootFile();
475 }
476 
478 {
479  if (myid != 0) { return; }
480 
481  std::string root_name = prefix_path + name + "_" +
483  ".mfem_root";
484  std::ofstream root_file(root_name);
485  root_file << GetVisItRootString();
486  if (!root_file)
487  {
488  error = WRITE_ERROR;
489  MFEM_WARNING("Error writing VisIt root file: " << root_name);
490  }
491 }
492 
494 {
495  DeleteAll();
496  time_step = 0.0;
497  error = No_Error;
498  cycle = cycle_;
499  std::string root_name = prefix_path + name + "_" +
501  ".mfem_root";
502  LoadVisItRootFile(root_name);
503  if (format != SERIAL_FORMAT || num_procs > 1)
504  {
505 #ifndef MFEM_USE_MPI
506  MFEM_WARNING("Cannot load parallel VisIt root file in serial.");
507  error = READ_ERROR;
508 #else
509  if (m_comm == MPI_COMM_NULL)
510  {
511  MFEM_WARNING("Cannot load parallel VisIt root file without MPI"
512  " communicator");
513  error = READ_ERROR;
514  }
515  else
516  {
517  // num_procs was read from the root file, check for consistency with
518  // the associated MPI_Comm, m_comm:
519  int comm_size;
520  MPI_Comm_size(m_comm, &comm_size);
521  if (comm_size != num_procs)
522  {
523  MFEM_WARNING("Processor number mismatch: VisIt root file: "
524  << num_procs << ", MPI_comm: " << comm_size);
525  error = READ_ERROR;
526  }
527  else
528  {
529  // myid was set when setting m_comm
530  }
531  }
532 #endif
533  }
534  if (!error)
535  {
536  LoadMesh(); // sets own_data to true, when there is no error
537  }
538  if (!error)
539  {
540  LoadFields();
541  }
542  if (error)
543  {
544  DeleteAll();
545  }
546 }
547 
548 void VisItDataCollection::LoadVisItRootFile(const std::string& root_name)
549 {
550  std::ifstream root_file(root_name);
551  std::stringstream buffer;
552  buffer << root_file.rdbuf();
553  if (!buffer)
554  {
555  error = READ_ERROR;
556  MFEM_WARNING("Error reading the VisIt root file: " << root_name);
557  }
558  else
559  {
560  ParseVisItRootString(buffer.str());
561  }
562 }
563 
565 {
566  // GetMeshFileName() uses 'serial', so we need to set it in advance.
567  serial = (format == SERIAL_FORMAT);
568  std::string mesh_fname = GetMeshFileName();
569  named_ifgzstream file(mesh_fname);
570  // TODO: in parallel, check for errors on all processors
571  if (!file)
572  {
573  error = READ_ERROR;
574  MFEM_WARNING("Unable to open mesh file: " << mesh_fname);
575  return;
576  }
577  // TODO: 1) load parallel mesh on one processor
578  if (format == SERIAL_FORMAT)
579  {
580  mesh = new Mesh(file, 1, 0, false);
581  serial = true;
582  }
583  else
584  {
585 #ifdef MFEM_USE_MPI
586  mesh = new ParMesh(m_comm, file);
587  serial = false;
588 #else
589  error = READ_ERROR;
590  MFEM_WARNING("Reading parallel format in serial is not supported");
591  return;
592 #endif
593  }
595  topo_dim = mesh->Dimension();
596  own_data = true;
597 }
598 
600 {
601  std::string path_left = prefix_path + name + "_" +
603  std::string path_right = "." + to_padded_string(myid, pad_digits_rank);
604 
605  field_map.clear();
606  for (FieldInfoMapIterator it = field_info_map.begin();
607  it != field_info_map.end(); ++it)
608  {
609  std::string fname = path_left + it->first + path_right;
610  mfem::ifgzstream file(fname);
611  // TODO: in parallel, check for errors on all processors
612  if (!file)
613  {
614  error = READ_ERROR;
615  MFEM_WARNING("Unable to open field file: " << fname);
616  return;
617  }
618  // TODO: 1) load parallel GridFunction on one processor
619  if (serial)
620  {
621  if ((it->second).association == "nodes")
622  {
623  field_map.Register(it->first, new GridFunction(mesh, file), own_data);
624  }
625  else if ((it->second).association == "elements")
626  {
627  q_field_map.Register(it->first, new QuadratureFunction(mesh, file), own_data);
628  }
629  }
630  else
631  {
632 #ifdef MFEM_USE_MPI
633  if ((it->second).association == "nodes")
634  {
636  it->first,
637  new ParGridFunction(dynamic_cast<ParMesh*>(mesh), file), own_data);
638  }
639  else if ((it->second).association == "elements")
640  {
641  q_field_map.Register(it->first, new QuadratureFunction(mesh, file), own_data);
642  }
643 #else
644  error = READ_ERROR;
645  MFEM_WARNING("Reading parallel format in serial is not supported");
646  return;
647 #endif
648  }
649  }
650 }
651 
653 {
654  // Get the path string (relative to where the root file is, i.e. no prefix).
655  std::string path_str =
657 
658  // We have to build the json tree inside out to get all the values in there
659  picojson::object top, dsets, main, mesh, fields, field, mtags, ftags;
660 
661  // Build the mesh data
662  std::string file_ext_format = ".%0" + to_string(pad_digits_rank) + "d";
663  mtags["spatial_dim"] = picojson::value(to_string(spatial_dim));
664  mtags["topo_dim"] = picojson::value(to_string(topo_dim));
665  mtags["max_lods"] = picojson::value(to_string(visit_max_levels_of_detail));
666  mesh["path"] = picojson::value(path_str + GetMeshShortFileName() +
667  file_ext_format);
668  mesh["tags"] = picojson::value(mtags);
669  mesh["format"] = picojson::value(to_string(format));
670 
671  // Build the fields data entries
672  for (FieldInfoMapIterator it = field_info_map.begin();
673  it != field_info_map.end(); ++it)
674  {
675  ftags["assoc"] = picojson::value((it->second).association);
676  ftags["comps"] = picojson::value(to_string((it->second).num_components));
677  ftags["lod"] = picojson::value(to_string((it->second).lod));
678  field["path"] = picojson::value(path_str + it->first + file_ext_format);
679  field["tags"] = picojson::value(ftags);
680  fields[it->first] = picojson::value(field);
681  }
682 
683  main["cycle"] = picojson::value(double(cycle));
684  main["time"] = picojson::value(time);
685  main["time_step"] = picojson::value(time_step);
686  main["domains"] = picojson::value(double(num_procs));
687  main["mesh"] = picojson::value(mesh);
688  if (!field_info_map.empty())
689  {
690  main["fields"] = picojson::value(fields);
691  }
692 
693  dsets["main"] = picojson::value(main);
694  top["dsets"] = picojson::value(dsets);
695 
696  return picojson::value(top).serialize(true);
697 }
698 
699 void VisItDataCollection::ParseVisItRootString(const std::string& json)
700 {
701  picojson::value top, dsets, main, mesh, fields;
702  std::string parse_err = picojson::parse(top, json);
703  if (!parse_err.empty())
704  {
705  error = READ_ERROR;
706  MFEM_WARNING("Unable to parse VisIt root data.");
707  return;
708  }
709 
710  // Process "main"
711  dsets = top.get("dsets");
712  main = dsets.get("main");
713  cycle = int(main.get("cycle").get<double>());
714  time = main.get("time").get<double>();
715  if (main.contains("time_step"))
716  {
717  time_step = main.get("time_step").get<double>();
718  }
719  num_procs = int(main.get("domains").get<double>());
720  mesh = main.get("mesh");
721  fields = main.get("fields");
722 
723  // ... Process "mesh"
724 
725  // Set the DataCollection::name using the mesh path
726  std::string path = mesh.get("path").get<std::string>();
727  size_t right_sep = path.rfind('_');
728  if (right_sep == std::string::npos)
729  {
730  error = READ_ERROR;
731  MFEM_WARNING("Unable to parse VisIt root data.");
732  return;
733  }
734  name = path.substr(0, right_sep);
735 
736  if (mesh.contains("format"))
737  {
738  format = to_int(mesh.get("format").get<std::string>());
739  }
740  spatial_dim = to_int(mesh.get("tags").get("spatial_dim").get<std::string>());
741  topo_dim = to_int(mesh.get("tags").get("topo_dim").get<std::string>());
743  to_int(mesh.get("tags").get("max_lods").get<std::string>());
744 
745  // ... Process "fields"
746  field_info_map.clear();
747  if (fields.is<picojson::object>())
748  {
749  picojson::object fields_obj = fields.get<picojson::object>();
750  for (picojson::object::iterator it = fields_obj.begin();
751  it != fields_obj.end(); ++it)
752  {
753  picojson::value tags = it->second.get("tags");
754  field_info_map[it->first] =
755  VisItFieldInfo(tags.get("assoc").get<std::string>(),
756  to_int(tags.get("comps").get<std::string>()));
757  }
758  }
759 }
760 
762  collection_name,
763  Mesh *mesh_)
764  : DataCollection(collection_name, mesh_),
765  levels_of_detail(1),
766  pv_data_format(VTKFormat::BINARY),
767  high_order_output(false),
768  restart_mode(false)
769 {
770  cycle = 0; // always include a valid cycle index in file names
771 
772  compression_level = -1; // default zlib compression level, equivalent to 6
773 #ifdef MFEM_USE_ZLIB
774  compression = true; // if we have zlib, enable compression
775 #else
776  compression = false; // otherwise, disable compression
777 #endif
778 }
779 
780 void ParaViewDataCollection::SetLevelsOfDetail(int levels_of_detail_)
781 {
782  levels_of_detail = levels_of_detail_;
783 }
784 
786 {
787  MFEM_WARNING("ParaViewDataCollection::Load() is not implemented!");
788 }
789 
791 {
793 }
794 
796 {
797  return "Cycle" + to_padded_string(cycle,pad_digits_cycle);
798 }
799 
801 {
802  return GeneratePVTUPath();
803 }
804 
806 {
807  return GetCollectionName() + ".pvd";
808 }
809 
811  const std::string &prefix)
812 {
813  return prefix + ".pvtu";
814 }
815 
817  const std::string &prefix, int rank)
818 {
819  return prefix + to_padded_string(rank, pad_digits_rank) + ".vtu";
820 }
821 
823 {
824  // add a new collection to the PDV file
825 
826  std::string col_path = GenerateCollectionPath();
827  // check if the directories are created
828  {
829  std::string path = col_path + "/" + GenerateVTUPath();
830  int error_code = create_directory(path, mesh, myid);
831  if (error_code)
832  {
833  error = WRITE_ERROR;
834  MFEM_WARNING("Error creating directory: " << path);
835  return; // do not even try to write the mesh
836  }
837  }
838  // the directory is created
839 
840  // create pvd file if needed. If we are not in restart mode, a new pvd file
841  // is always created. In restart mode, we keep any previously defined
842  // timestep values as long as they are less than the currently defined time.
843 
844  if (myid == 0 && !pvd_stream.is_open())
845  {
846  std::string pvdname = col_path + "/" + GeneratePVDFileName();
847 
848  bool write_header = true;
849  std::ifstream pvd_in;
850  if (restart_mode && (pvd_in.open(pvdname,std::ios::binary),pvd_in.good()))
851  {
852  // PVD file exists and restart mode enabled: preserve existing time
853  // steps less than the current time.
854  std::fstream::pos_type pos_begin = pvd_in.tellg();
855  std::fstream::pos_type pos_end = pos_begin;
856 
857  std::regex regexp("timestep=\"([^[:space:]]+)\".*file=\"Cycle(\\d+)");
858  std::smatch match;
859 
860  std::string line;
861  while (getline(pvd_in,line))
862  {
863  if (regex_search(line,match,regexp))
864  {
865  MFEM_ASSERT(match.size() == 3, "Unable to parse DataSet");
866  double tvalue = std::stod(match[1]);
867  if (tvalue >= GetTime()) { break; }
868  int cvalue = std::stoi(match[2]);
869  MFEM_VERIFY(cvalue < GetCycle(), "Cycle " << GetCycle() <<
870  " is too small for restart mode: trying to overwrite"
871  " existing data.");
872  pos_end = pvd_in.tellg();
873  }
874  }
875  // Since pvd_in is opened in binary mode, count will store the number
876  // of bytes from the beginning of the file until the desired insertion
877  // point (in text mode on Windows this is not the case).
878  size_t count = pos_end - pos_begin;
879  if (count != 0)
880  {
881  write_header = false;
882  std::vector<char> buf(count);
883  // Read the contents of the PVD file, from the beginning to the
884  // insertion point.
885  pvd_in.clear();
886  pvd_in.seekg(pos_begin);
887  pvd_in.read(buf.data(), count);
888  pvd_in.close();
889  // Open the PVD file in truncate mode to delete the previous
890  // contents. Open in binary mode to write the data buffer without
891  // converting \r\n to \r\r\n on Windows.
892  pvd_stream.open(pvdname,std::ios::out|std::ios::trunc|std::ios::binary);
893  pvd_stream.write(buf.data(), count);
894  // Close and reopen the file in text mode, appending to the end.
895  pvd_stream.close();
896  pvd_stream.open(pvdname,std::ios::in|std::ios::out|std::ios::ate);
897  }
898  }
899  if (write_header)
900  {
901  // Initialize new pvd file.
902  pvd_stream.open(pvdname,std::ios::out|std::ios::trunc);
903  pvd_stream << "<?xml version=\"1.0\"?>\n";
904  pvd_stream << "<VTKFile type=\"Collection\" version=\"0.1\"";
905  pvd_stream << " byte_order=\"" << VTKByteOrder() << "\">\n";
906  pvd_stream << "<Collection>" << std::endl;
907  }
908  }
909 
910  std::string vtu_prefix = col_path + "/" + GenerateVTUPath() + "/";
911 
912  // Save the local part of the mesh and grid functions fields to the local
913  // VTU file
914  {
915  std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
916  os.precision(precision);
917  SaveDataVTU(os, levels_of_detail);
918  }
919 
920  // Save the local part of the quadrature function fields
921  for (const auto &qfield : q_field_map)
922  {
923  const std::string &field_name = qfield.first;
924  std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
925  qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
926  }
927 
928  // MPI rank 0 also creates a "PVTU" file that points to all of the separately
929  // written VTU files.
930  // This file path is then appended to the PVD file.
931  if (myid == 0)
932  {
933  // Create the main PVTU file
934  {
935  std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
936  WritePVTUHeader(pvtu_out);
937 
938  // Grid function fields
939  pvtu_out << "<PPointData>\n";
940  for (auto &field_it : field_map)
941  {
942  int vec_dim = field_it.second->VectorDim();
943  pvtu_out << "<PDataArray type=\"" << GetDataTypeString()
944  << "\" Name=\"" << field_it.first
945  << "\" NumberOfComponents=\"" << vec_dim << "\" "
946  << "format=\"" << GetDataFormatString() << "\" />\n";
947  }
948  pvtu_out << "</PPointData>\n";
949  // Element attributes
950  pvtu_out << "<PCellData>\n";
951  pvtu_out << "\t<PDataArray type=\"Int32\" Name=\"" << "attribute"
952  << "\" NumberOfComponents=\"1\""
953  << " format=\"" << GetDataFormatString() << "\"/>\n";
954  pvtu_out << "</PCellData>\n";
955 
956  WritePVTUFooter(pvtu_out, "proc");
957  }
958 
959  // Add the latest PVTU to the PVD
960  pvd_stream << "<DataSet timestep=\"" << GetTime()
961  << "\" group=\"\" part=\"" << 0 << "\" file=\""
962  << GeneratePVTUPath() + "/" + GeneratePVTUFileName("data")
963  << "\" name=\"mesh\"/>\n";
964 
965  // Create PVTU files for each quadrature field and add them to the PVD
966  // file
967  for (auto &q_field : q_field_map)
968  {
969  const std::string &q_field_name = q_field.first;
970  std::string q_fname = GeneratePVTUPath() + "/"
971  + GeneratePVTUFileName(q_field_name);
972 
973  std::ofstream pvtu_out(col_path + "/" + q_fname);
974  WritePVTUHeader(pvtu_out);
975  int vec_dim = q_field.second->GetVDim();
976  pvtu_out << "<PPointData>\n";
977  pvtu_out << "<PDataArray type=\"" << GetDataTypeString()
978  << "\" Name=\"" << q_field_name
979  << "\" NumberOfComponents=\"" << vec_dim << "\" "
980  << "format=\"" << GetDataFormatString() << "\" />\n";
981  pvtu_out << "</PPointData>\n";
982  WritePVTUFooter(pvtu_out, q_field_name);
983 
984  pvd_stream << "<DataSet timestep=\"" << GetTime()
985  << "\" group=\"\" part=\"" << 0 << "\" file=\""
986  << q_fname << "\" name=\"" << q_field_name << "\"/>\n";
987  }
988  pvd_stream.flush();
989  // Move the insertion point before the closing collection tag, so that
990  // the PVD file is valid even when writing incrementally.
991  std::fstream::pos_type pos = pvd_stream.tellp();
992  pvd_stream << "</Collection>\n";
993  pvd_stream << "</VTKFile>" << std::endl;
994  pvd_stream.seekp(pos);
995  }
996 }
997 
999 {
1000  os << "<?xml version=\"1.0\"?>\n";
1001  os << "<VTKFile type=\"PUnstructuredGrid\"";
1002  os << " version =\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
1003  os << "<PUnstructuredGrid GhostLevel=\"0\">\n";
1004 
1005  os << "<PPoints>\n";
1006  os << "\t<PDataArray type=\"" << GetDataTypeString() << "\" ";
1007  os << " Name=\"Points\" NumberOfComponents=\"3\""
1008  << " format=\"" << GetDataFormatString() << "\"/>\n";
1009  os << "</PPoints>\n";
1010 
1011  os << "<PCells>\n";
1012  os << "\t<PDataArray type=\"Int32\" ";
1013  os << " Name=\"connectivity\" NumberOfComponents=\"1\""
1014  << " format=\"" << GetDataFormatString() << "\"/>\n";
1015  os << "\t<PDataArray type=\"Int32\" ";
1016  os << " Name=\"offsets\" NumberOfComponents=\"1\""
1017  << " format=\"" << GetDataFormatString() << "\"/>\n";
1018  os << "\t<PDataArray type=\"UInt8\" ";
1019  os << " Name=\"types\" NumberOfComponents=\"1\""
1020  << " format=\"" << GetDataFormatString() << "\"/>\n";
1021  os << "</PCells>\n";
1022 }
1023 
1025  const std::string &vtu_prefix)
1026 {
1027  for (int ii=0; ii<num_procs; ii++)
1028  {
1029  std::string vtu_filename = GenerateVTUFileName(vtu_prefix, ii);
1030  os << "<Piece Source=\"" << vtu_filename << "\"/>\n";
1031  }
1032  os << "</PUnstructuredGrid>\n";
1033  os << "</VTKFile>\n";
1034 }
1035 
1036 void ParaViewDataCollection::SaveDataVTU(std::ostream &os, int ref)
1037 {
1038  os << "<VTKFile type=\"UnstructuredGrid\"";
1039  if (GetCompressionLevel() != 0)
1040  {
1041  os << " compressor=\"vtkZLibDataCompressor\"";
1042  }
1043  os << " version=\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
1044  os << "<UnstructuredGrid>\n";
1045  mesh->PrintVTU(os,ref,pv_data_format,high_order_output,GetCompressionLevel());
1046 
1047  // dump out the grid functions as point data
1048  os << "<PointData >\n";
1049  // save the grid functions
1050  // iterate over all grid functions
1051  for (FieldMapIterator it=field_map.begin(); it!=field_map.end(); ++it)
1052  {
1053  SaveGFieldVTU(os,ref,it);
1054  }
1055  os << "</PointData>\n";
1056  // close the mesh
1057  os << "</Piece>\n"; // close the piece open in the PrintVTU method
1058  os << "</UnstructuredGrid>\n";
1059  os << "</VTKFile>" << std::endl;
1060 }
1061 
1062 void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
1063  const FieldMapIterator &it)
1064 {
1065  RefinedGeometry *RefG;
1066  Vector val;
1067  DenseMatrix vval, pmat;
1068  std::vector<char> buf;
1069  int vec_dim = it->second->VectorDim();
1070  os << "<DataArray type=\"" << GetDataTypeString()
1071  << "\" Name=\"" << it->first
1072  << "\" NumberOfComponents=\"" << vec_dim << "\""
1073  << " format=\"" << GetDataFormatString() << "\" >" << '\n';
1074  if (vec_dim == 1)
1075  {
1076  // scalar data
1077  for (int i = 0; i < mesh->GetNE(); i++)
1078  {
1079  RefG = GlobGeometryRefiner.Refine(
1080  mesh->GetElementBaseGeometry(i), ref_, 1);
1081  it->second->GetValues(i, RefG->RefPts, val, pmat);
1082  for (int j = 0; j < val.Size(); j++)
1083  {
1084  WriteBinaryOrASCII(os, buf, val(j), "\n", pv_data_format);
1085  }
1086  }
1087  }
1088  else
1089  {
1090  // vector data
1091  for (int i = 0; i < mesh->GetNE(); i++)
1092  {
1093  RefG = GlobGeometryRefiner.Refine(
1094  mesh->GetElementBaseGeometry(i), ref_, 1);
1095  it->second->GetVectorValues(i, RefG->RefPts, vval, pmat);
1096  for (int jj = 0; jj < vval.Width(); jj++)
1097  {
1098  for (int ii = 0; ii < vval.Height(); ii++)
1099  {
1100  WriteBinaryOrASCII(os, buf, vval(ii,jj), " ", pv_data_format);
1101  }
1102  if (pv_data_format == VTKFormat::ASCII) { os << '\n'; }
1103  }
1104  }
1105  }
1106 
1107  if (IsBinaryFormat())
1108  {
1109  WriteVTKEncodedCompressed(os,buf.data(),buf.size(),GetCompressionLevel());
1110  os << '\n';
1111  }
1112  os << "</DataArray>" << std::endl;
1113 }
1114 
1116 {
1117  pv_data_format = fmt;
1118 }
1119 
1121 {
1122  return pv_data_format != VTKFormat::ASCII;
1123 }
1124 
1125 void ParaViewDataCollection::SetHighOrderOutput(bool high_order_output_)
1126 {
1127  high_order_output = high_order_output_;
1128 }
1129 
1131 {
1132  MFEM_ASSERT(compression_level_ >= -1 && compression_level_ <= 9,
1133  "Compression level must be between -1 and 9 (inclusive).");
1134  compression_level = compression_level_;
1135  compression = compression_level_ != 0;
1136 }
1137 
1139 {
1140  compression = compression_;
1141 }
1142 
1144 {
1145  restart_mode = restart_mode_;
1146 }
1147 
1149 {
1150  if (pv_data_format == VTKFormat::ASCII)
1151  {
1152  return "ascii";
1153  }
1154  else
1155  {
1156  return "binary";
1157  }
1158 }
1159 
1161 {
1162  if (pv_data_format==VTKFormat::ASCII || pv_data_format==VTKFormat::BINARY)
1163  {
1164  return "Float64";
1165  }
1166  else
1167  {
1168  return "Float32";
1169  }
1170 }
1171 
1173 {
1174  return compression ? compression_level : 0;
1175 }
1176 
1177 } // end namespace MFEM
virtual void SaveMesh()
Save the mesh, creating the collection directory.
void SaveOneField(const FieldMapIterator &it)
Save one field to disk, assuming the collection directory exists.
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
int GetNE() const
Return the number of entities.
Definition: qspace.hpp:61
void WritePVTUHeader(std::ostream &out)
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int format
Output mesh format: see the Format enumeration.
int error
Error state.
void SetDataFormat(VTKFormat fmt)
bool appendRankToFileName
Append rank to any output file names.
iterator begin()
Returns a begin iterator to the registered fields.
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
int GetCompressionLevel() const
If compression is enabled, return the compression level, otherwise return 0.
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
virtual void RegisterQField(const std::string &q_field_name, QuadratureFunction *qf)
Add a QuadratureFunction to the collection.
void DeleteAll()
Delete data owned by the DataCollection including field information.
iterator end()
Returns an end iterator to the registered fields.
int pad_digits_cycle
Number of digits used for the cycle and MPI rank in filenames.
GFieldMap::iterator FieldMapIterator
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:557
Helper class for VisIt visualization data.
const char * GetDataTypeString() const
std::string GetVisItRootString()
Prepare the VisIt root file in JSON format for the current collection.
bool own_data
Should the collection delete its mesh and fields.
Data arrays will be written in ASCII format.
virtual void SetCompression(bool comp)
Set the flag for use of gz compressed files.
int GetCycle() const
Get time cycle (for time-dependent simulations)
void ParseVisItRootString(const std::string &json)
Read in a VisIt root file in JSON format.
std::string GetMeshShortFileName() const
virtual void SaveQField(const std::string &q_field_name)
Save one q-field, assuming the collection directory already exists.
void WriteBinaryOrASCII(std::ostream &os, std::vector< char > &buf, const T &val, const char *suffix, VTKFormat format)
Write either ASCII data to the stream or binary data to the buffer depending on the given format...
Definition: vtk.hpp:147
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
MPI_Comm m_comm
Associated MPI communicator.
virtual void Load(int cycle_=0) override
Load the collection - not implemented in the ParaView writer.
double time
Physical time (for time-dependent simulations)
std::map< std::string, VisItFieldInfo >::iterator FieldInfoMapIterator
iterator find(const std::string &fname)
Returns an iterator to the field fname.
void SetCompressionLevel(int compression_level_)
Set the zlib compression level.
Mesh * mesh
The (common) mesh for the collected fields.
void SaveOneQField(const QFieldMapIterator &it)
Save one q-field to disk, assuming the collection directory exists.
void ParPrint(std::ostream &out) const
Save the mesh in a parallel mesh format.
Definition: pmesh.cpp:6295
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
bool serial
Serial or parallel run? False iff mesh is a ParMesh.
int cycle
Time cycle; for time-dependent simulations cycle >= 0, otherwise = -1.
VisItDataCollection(const std::string &collection_name, Mesh *mesh_=NULL)
Constructor. The collection name is used when saving the data.
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:413
virtual void Save()
Save the collection to disk.
std::string GenerateVTUFileName(const std::string &prefix, int rank)
std::string to_padded_string(int i, int digits)
Convert an integer to a 0-padded string with the given number of digits.
Definition: text.hpp:54
int to_int(const std::string &str)
Convert a string to an int.
Definition: text.hpp:62
void SaveDataVTU(std::ostream &out, int ref)
const IntegrationRule & GetIntRule(int idx) const
Get the IntegrationRule associated with entity (element or face) idx.
Definition: qfunction.hpp:176
std::string GetFieldFileName(const std::string &field_name) const
VTKFormat
Data array format for VTK and VTU files.
Definition: vtk.hpp:98
void SaveGFieldVTU(std::ostream &out, int ref_, const FieldMapIterator &it)
virtual void SetFormat(int fmt)
Set the desired output mesh and data format.
static int create_directory(const std::string &dir_name, const Mesh *mesh, int myid)
Mesh * GetMesh() const
Returns the mesh.
Definition: qspace.hpp:64
void WritePVTUFooter(std::ostream &out, const std::string &vtu_prefix)
void Register(const std::string &fname, T *field, bool own_data)
Register field field with name fname.
DataCollection(const std::string &collection_name, Mesh *mesh_=NULL)
Initialize the collection with its name and Mesh.
virtual void SetMesh(Mesh *new_mesh) override
Set/change the mesh associated with the collection.
void DeleteData()
Delete data owned by the DataCollection keeping field information.
virtual void Load(int cycle_=0) override
Load the collection based on its VisIt data (described in its root file)
void SetLevelsOfDetail(int levels_of_detail)
Set VisIt parameter: default levels of detail for the MultiresControl.
IntegrationRule RefPts
Definition: geom.hpp:314
void SetHighOrderOutput(bool high_order_output_)
static const int pad_digits_default
Default value for pad_digits_*.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1099
void SetCompression(bool compression_) override
virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition: geom.cpp:1732
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
void WriteVTKEncodedCompressed(std::ostream &os, const void *bytes, uint32_t nbytes, int compression_level)
Outputs encoded binary data in the base 64 format needed by VTK.
Definition: vtk.cpp:559
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
std::string GeneratePVTUFileName(const std::string &prefix)
virtual void RegisterQField(const std::string &q_field_name, QuadratureFunction *qf) override
Add a quadrature function to the collection and update the root file.
void PrintVTU(std::ostream &os, int ref=1, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0, bool bdr_elements=false)
Definition: mesh.cpp:10919
int GetMyRank() const
Definition: pmesh.hpp:353
bool IsBinaryFormat() const
Returns true if the output format is BINARY or BINARY32, false if ASCII.
std::map< std::string, VisItFieldInfo > field_info_map
void UseRestartMode(bool restart_mode_)
QFieldMap::iterator QFieldMapIterator
int precision
Precision (number of digits) used for the text output of doubles.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
virtual void Save() override
Save the collection and a VisIt root file.
int myid
MPI rank (in parallel)
virtual void SaveField(const std::string &field_name)
Save one field, assuming the collection directory already exists.
static const int precision_default
Default value for precision.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
void SetMaxLevelsOfDetail(int max_levels_of_detail)
Set VisIt parameter: maximum levels of detail for the MultiresControl.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
virtual ~DataCollection()
Delete the mesh and fields if owned by the collection.
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:82
int GetNRanks() const
Definition: pmesh.hpp:352
int VectorDim() const
Definition: gridfunc.cpp:323
void clear()
Clears the map of registered fields without reclaiming memory.
void SaveRootFile()
Save a VisIt root file for the collection.
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
Definition: vtk.cpp:602
ParaViewDataCollection(const std::string &collection_name, mfem::Mesh *mesh_=NULL)
Constructor. The collection name is used when saving the data.
std::string prefix_path
A path where the directory with results is saved. If not empty, it has &#39;/&#39; at the end...
virtual void Load(int cycle_=0)
Load the collection. Not implemented in the base class DataCollection.
void SetLevelsOfDetail(int levels_of_detail_)
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
void LoadVisItRootFile(const std::string &root_name)
virtual void SetMesh(Mesh *new_mesh)
Set/change the mesh associated with the collection.
std::string GetMeshFileName() const
virtual void Save() override
std::string name
Name of the collection, used as a directory name when saving.
void DeleteAll()
Delete all data owned by VisItDataCollection including field data information.
void DeleteData(bool own_data)
Clear all associations between names and fields.
Class for parallel grid function.
Definition: pgridfunc.hpp:32
double time_step
Time step i.e. delta_t (for time-dependent simulations)
int num_procs
Number of MPI ranks (in parallel)
Class for parallel meshes.
Definition: pmesh.hpp:32
const std::string & GetCollectionName() const
Get the name of the collection.
const char * GetDataFormatString() const
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:23
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
int main()
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
double GetTime() const
Get physical time (for time-dependent simulations)