MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
datacollection.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 = false;
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);
433  visit_levels_of_detail = std::max(visit_levels_of_detail, 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  {
444  mesh->GetElementBaseGeometry(e),
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 #ifdef MFEM_USE_ZLIB
771  compression = -1; // default zlib compression level, equivalent to 6
772 #else
773  compression = 0;
774 #endif
775 }
776 
777 void ParaViewDataCollection::SetLevelsOfDetail(int levels_of_detail_)
778 {
779  levels_of_detail = levels_of_detail_;
780 }
781 
783 {
784  MFEM_WARNING("ParaViewDataCollection::Load() is not implemented!");
785 }
786 
788 {
790 }
791 
793 {
794  return "Cycle" + to_padded_string(cycle,pad_digits_cycle);
795 }
796 
798 {
799  return GeneratePVTUPath();
800 }
801 
803 {
804  return GetCollectionName() + ".pvd";
805 }
806 
808  const std::string &prefix)
809 {
810  return prefix + ".pvtu";
811 }
812 
814  const std::string &prefix, int rank)
815 {
816  return prefix + to_padded_string(rank, pad_digits_rank) + ".vtu";
817 }
818 
820 {
821  // add a new collection to the PDV file
822 
823  std::string col_path = GenerateCollectionPath();
824  // check if the directories are created
825  {
826  std::string path = col_path + "/" + GenerateVTUPath();
827  int error_code = create_directory(path, mesh, myid);
828  if (error_code)
829  {
830  error = WRITE_ERROR;
831  MFEM_WARNING("Error creating directory: " << path);
832  return; // do not even try to write the mesh
833  }
834  }
835  // the directory is created
836 
837  // create pvd file if needed. If we are not in restart mode, a new pvd file
838  // is always created. In restart mode, we keep any previously defined
839  // timestep values as long as they are less than the currently defined time.
840 
841  if (myid == 0 && !pvd_stream.is_open())
842  {
843  std::string pvdname = col_path + "/" + GeneratePVDFileName();
844 
845  bool write_header = true;
846  std::ifstream pvd_in;
847  if (restart_mode && (pvd_in.open(pvdname,std::ios::binary),pvd_in.good()))
848  {
849  // PVD file exists and restart mode enabled: preserve existing time
850  // steps less than the current time.
851  std::fstream::pos_type pos_begin = pvd_in.tellg();
852  std::fstream::pos_type pos_end = pos_begin;
853 
854  std::regex regexp("timestep=\"([^[:space:]]+)\".*file=\"Cycle(\\d+)");
855  std::smatch match;
856 
857  std::string line;
858  while (getline(pvd_in,line))
859  {
860  if (regex_search(line,match,regexp))
861  {
862  MFEM_ASSERT(match.size() == 3, "Unable to parse DataSet");
863  double tvalue = std::stod(match[1]);
864  if (tvalue >= GetTime()) { break; }
865  int cvalue = std::stoi(match[2]);
866  MFEM_VERIFY(cvalue < GetCycle(), "Cycle " << GetCycle() <<
867  " is too small for restart mode: trying to overwrite"
868  " existing data.");
869  pos_end = pvd_in.tellg();
870  }
871  }
872  // Since pvd_in is opened in binary mode, count will store the number
873  // of bytes from the beginning of the file until the desired insertion
874  // point (in text mode on Windows this is not the case).
875  size_t count = pos_end - pos_begin;
876  if (count != 0)
877  {
878  write_header = false;
879  std::vector<char> buf(count);
880  // Read the contents of the PVD file, from the beginning to the
881  // insertion point.
882  pvd_in.clear();
883  pvd_in.seekg(pos_begin);
884  pvd_in.read(buf.data(), count);
885  pvd_in.close();
886  // Open the PVD file in truncate mode to delete the previous
887  // contents. Open in binary mode to write the data buffer without
888  // converting \r\n to \r\r\n on Windows.
889  pvd_stream.open(pvdname,std::ios::out|std::ios::trunc|std::ios::binary);
890  pvd_stream.write(buf.data(), count);
891  // Close and reopen the file in text mode, appending to the end.
892  pvd_stream.close();
893  pvd_stream.open(pvdname,std::ios::in|std::ios::out|std::ios::ate);
894  }
895  }
896  if (write_header)
897  {
898  // Initialize new pvd file.
899  pvd_stream.open(pvdname,std::ios::out|std::ios::trunc);
900  pvd_stream << "<?xml version=\"1.0\"?>\n";
901  pvd_stream << "<VTKFile type=\"Collection\" version=\"0.1\"";
902  pvd_stream << " byte_order=\"" << VTKByteOrder() << "\">\n";
903  pvd_stream << "<Collection>" << std::endl;
904  }
905  }
906 
907  std::string vtu_prefix = col_path + "/" + GenerateVTUPath() + "/";
908 
909  // Save the local part of the mesh and grid functions fields to the local
910  // VTU file
911  {
912  std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
913  os.precision(precision);
914  SaveDataVTU(os, levels_of_detail);
915  }
916 
917  // Save the local part of the quadrature function fields
918  for (const auto &qfield : q_field_map)
919  {
920  const std::string &field_name = qfield.first;
921  std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
922  qfield.second->SaveVTU(os, pv_data_format, compression);
923  }
924 
925  // MPI rank 0 also creates a "PVTU" file that points to all of the separately
926  // written VTU files.
927  // This file path is then appended to the PVD file.
928  if (myid == 0)
929  {
930  // Create the main PVTU file
931  {
932  std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
933  WritePVTUHeader(pvtu_out);
934 
935  // Grid function fields
936  pvtu_out << "<PPointData>\n";
937  for (auto &field_it : field_map)
938  {
939  int vec_dim = field_it.second->VectorDim();
940  pvtu_out << "<PDataArray type=\"" << GetDataTypeString()
941  << "\" Name=\"" << field_it.first
942  << "\" NumberOfComponents=\"" << vec_dim << "\" "
943  << "format=\"" << GetDataFormatString() << "\" />\n";
944  }
945  pvtu_out << "</PPointData>\n";
946  // Element attributes
947  pvtu_out << "<PCellData>\n";
948  pvtu_out << "\t<PDataArray type=\"Int32\" Name=\"" << "attribute"
949  << "\" NumberOfComponents=\"1\""
950  << " format=\"" << GetDataFormatString() << "\"/>\n";
951  pvtu_out << "</PCellData>\n";
952 
953  WritePVTUFooter(pvtu_out, "proc");
954  }
955 
956  // Add the latest PVTU to the PVD
957  pvd_stream << "<DataSet timestep=\"" << GetTime()
958  << "\" group=\"\" part=\"" << 0 << "\" file=\""
959  << GeneratePVTUPath() + "/" + GeneratePVTUFileName("data")
960  << "\" name=\"mesh\"/>\n";
961 
962  // Create PVTU files for each quadrature field and add them to the PVD
963  // file
964  for (auto &q_field : q_field_map)
965  {
966  const std::string &q_field_name = q_field.first;
967  std::string q_fname = GeneratePVTUPath() + "/"
968  + GeneratePVTUFileName(q_field_name);
969 
970  std::ofstream pvtu_out(col_path + "/" + q_fname);
971  WritePVTUHeader(pvtu_out);
972  int vec_dim = q_field.second->GetVDim();
973  pvtu_out << "<PPointData>\n";
974  pvtu_out << "<PDataArray type=\"" << GetDataTypeString()
975  << "\" Name=\"" << q_field_name
976  << "\" NumberOfComponents=\"" << vec_dim << "\" "
977  << "format=\"" << GetDataFormatString() << "\" />\n";
978  pvtu_out << "</PPointData>\n";
979  WritePVTUFooter(pvtu_out, q_field_name);
980 
981  pvd_stream << "<DataSet timestep=\"" << GetTime()
982  << "\" group=\"\" part=\"" << 0 << "\" file=\""
983  << q_fname << "\" name=\"" << q_field_name << "\"/>\n";
984  }
985  pvd_stream.flush();
986  // Move the insertion point before the closing collection tag, so that
987  // the PVD file is valid even when writing incrementally.
988  std::fstream::pos_type pos = pvd_stream.tellp();
989  pvd_stream << "</Collection>\n";
990  pvd_stream << "</VTKFile>" << std::endl;
991  pvd_stream.seekp(pos);
992  }
993 }
994 
996 {
997  os << "<?xml version=\"1.0\"?>\n";
998  os << "<VTKFile type=\"PUnstructuredGrid\"";
999  os << " version =\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
1000  os << "<PUnstructuredGrid GhostLevel=\"0\">\n";
1001 
1002  os << "<PPoints>\n";
1003  os << "\t<PDataArray type=\"" << GetDataTypeString() << "\" ";
1004  os << " Name=\"Points\" NumberOfComponents=\"3\""
1005  << " format=\"" << GetDataFormatString() << "\"/>\n";
1006  os << "</PPoints>\n";
1007 
1008  os << "<PCells>\n";
1009  os << "\t<PDataArray type=\"Int32\" ";
1010  os << " Name=\"connectivity\" NumberOfComponents=\"1\""
1011  << " format=\"" << GetDataFormatString() << "\"/>\n";
1012  os << "\t<PDataArray type=\"Int32\" ";
1013  os << " Name=\"offsets\" NumberOfComponents=\"1\""
1014  << " format=\"" << GetDataFormatString() << "\"/>\n";
1015  os << "\t<PDataArray type=\"UInt8\" ";
1016  os << " Name=\"types\" NumberOfComponents=\"1\""
1017  << " format=\"" << GetDataFormatString() << "\"/>\n";
1018  os << "</PCells>\n";
1019 }
1020 
1022  const std::string &vtu_prefix)
1023 {
1024  for (int ii=0; ii<num_procs; ii++)
1025  {
1026  std::string vtu_filename = GenerateVTUFileName(vtu_prefix, ii);
1027  os << "<Piece Source=\"" << vtu_filename << "\"/>\n";
1028  }
1029  os << "</PUnstructuredGrid>\n";
1030  os << "</VTKFile>\n";
1031 }
1032 
1033 void ParaViewDataCollection::SaveDataVTU(std::ostream &os, int ref)
1034 {
1035  os << "<VTKFile type=\"UnstructuredGrid\"";
1036  if (compression != 0)
1037  {
1038  os << " compressor=\"vtkZLibDataCompressor\"";
1039  }
1040  os << " version=\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
1041  os << "<UnstructuredGrid>\n";
1042  mesh->PrintVTU(os,ref,pv_data_format,high_order_output,compression);
1043 
1044  // dump out the grid functions as point data
1045  os << "<PointData >\n";
1046  // save the grid functions
1047  // iterate over all grid functions
1048  for (FieldMapIterator it=field_map.begin(); it!=field_map.end(); ++it)
1049  {
1050  SaveGFieldVTU(os,ref,it);
1051  }
1052  os << "</PointData>\n";
1053  // close the mesh
1054  os << "</Piece>\n"; // close the piece open in the PrintVTU method
1055  os << "</UnstructuredGrid>\n";
1056  os << "</VTKFile>" << std::endl;
1057 }
1058 
1059 void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
1060  const FieldMapIterator &it)
1061 {
1062  RefinedGeometry *RefG;
1063  Vector val;
1064  DenseMatrix vval, pmat;
1065  std::vector<char> buf;
1066  int vec_dim = it->second->VectorDim();
1067  os << "<DataArray type=\"" << GetDataTypeString()
1068  << "\" Name=\"" << it->first
1069  << "\" NumberOfComponents=\"" << vec_dim << "\""
1070  << " format=\"" << GetDataFormatString() << "\" >" << '\n';
1071  if (vec_dim == 1)
1072  {
1073  // scalar data
1074  for (int i = 0; i < mesh->GetNE(); i++)
1075  {
1076  RefG = GlobGeometryRefiner.Refine(
1077  mesh->GetElementBaseGeometry(i), ref_, 1);
1078  it->second->GetValues(i, RefG->RefPts, val, pmat);
1079  for (int j = 0; j < val.Size(); j++)
1080  {
1081  WriteBinaryOrASCII(os, buf, val(j), "\n", pv_data_format);
1082  }
1083  }
1084  }
1085  else
1086  {
1087  // vector data
1088  for (int i = 0; i < mesh->GetNE(); i++)
1089  {
1090  RefG = GlobGeometryRefiner.Refine(
1091  mesh->GetElementBaseGeometry(i), ref_, 1);
1092  it->second->GetVectorValues(i, RefG->RefPts, vval, pmat);
1093  for (int jj = 0; jj < vval.Width(); jj++)
1094  {
1095  for (int ii = 0; ii < vval.Height(); ii++)
1096  {
1097  WriteBinaryOrASCII(os, buf, vval(ii,jj), " ", pv_data_format);
1098  }
1099  if (pv_data_format == VTKFormat::ASCII) { os << '\n'; }
1100  }
1101  }
1102  }
1103 
1104  if (IsBinaryFormat())
1105  {
1106  WriteVTKEncodedCompressed(os,buf.data(),buf.size(),compression);
1107  os << '\n';
1108  }
1109  os << "</DataArray>" << std::endl;
1110 }
1111 
1113 {
1114  pv_data_format = fmt;
1115 }
1116 
1118 {
1119  return pv_data_format != VTKFormat::ASCII;
1120 }
1121 
1122 void ParaViewDataCollection::SetHighOrderOutput(bool high_order_output_)
1123 {
1124  high_order_output = high_order_output_;
1125 }
1126 
1128 {
1129  MFEM_ASSERT(compression_level_ >= -1 && compression_level_ <= 9,
1130  "Compression level must be between -1 and 9 (inclusive).");
1131  compression = compression_level_;
1132 }
1133 
1135 {
1136  // If we are enabling compression, and it was disabled previously, use the
1137  // default compression level. Otherwise, leave the compression level
1138  // unchanged.
1139  if (compression_ && compression == 0)
1140  {
1141  SetCompressionLevel(-1);
1142  }
1143 }
1144 
1146 {
1147  restart_mode = restart_mode_;
1148 }
1149 
1151 {
1152  if (pv_data_format == VTKFormat::ASCII)
1153  {
1154  return "ascii";
1155  }
1156  else
1157  {
1158  return "binary";
1159  }
1160 }
1161 
1163 {
1164  if (pv_data_format==VTKFormat::ASCII || pv_data_format==VTKFormat::BINARY)
1165  {
1166  return "Float64";
1167  }
1168  else
1169  {
1170  return "Float32";
1171  }
1172 }
1173 
1174 } // end namespace MFEM
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
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.
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)
double GetTime() const
Get physical time (for time-dependent simulations)
bool appendRankToFileName
Append rank to any output file names.
iterator begin()
Returns a begin iterator to the registered fields.
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.
const char * GetDataTypeString() const
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:73
const std::string & GetCollectionName() const
Get the name of the collection.
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 VectorDim() const
Definition: gridfunc.cpp:324
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
Helper class for VisIt visualization data.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
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.
void ParseVisItRootString(const std::string &json)
Read in a VisIt root file in JSON format.
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1069
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:145
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
int GetNRanks() const
Definition: pmesh.hpp:352
iterator find(const std::string &fname)
Returns an iterator to the field fname.
void SetCompressionLevel(int compression_level_)
const IntegrationRule & GetIntRule(int idx) const
Get the IntegrationRule associated with entity (element or face) idx.
Definition: qfunction.hpp:164
std::string GetMeshShortFileName() const
Mesh * mesh
The (common) mesh for the collected fields.
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
void SaveOneQField(const QFieldMapIterator &it)
Save one q-field to disk, assuming the collection directory exists.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
bool serial
Serial or parallel run? False iff mesh is a ParMesh.
int cycle
Time cycle; for time-dependent simulations cycle &gt;= 0, otherwise = -1.
VisItDataCollection(const std::string &collection_name, Mesh *mesh_=NULL)
Constructor. The collection name is used when saving the data.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
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
int GetCycle() const
Get time cycle (for time-dependent simulations)
void SaveDataVTU(std::ostream &out, int ref)
VTKFormat
Data array format for VTK and VTU files.
Definition: vtk.hpp:96
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)
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
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:282
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:668
int Dimension() const
Definition: mesh.hpp:1006
const char * GetDataFormatString() const
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
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.
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:353
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:10503
int SpaceDimension() const
Definition: mesh.hpp:1007
std::map< std::string, VisItFieldInfo > field_info_map
void UseRestartMode(bool restart_mode_)
int GetMyRank() const
Definition: pmesh.hpp:353
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:443
QFieldMap::iterator QFieldMapIterator
int precision
Precision (number of digits) used for the text output of doubles.
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
bool IsBinaryFormat() const
Returns true if the output format is BINARY or BINARY32, false if ASCII.
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.
void SetMaxLevelsOfDetail(int max_levels_of_detail)
Set VisIt parameter: maximum levels of detail for the MultiresControl.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
virtual ~DataCollection()
Delete the mesh and fields if owned by the collection.
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:70
std::string GetFieldFileName(const std::string &field_name) const
void clear()
Clears the map of registered fields without reclaiming memory.
void SaveRootFile()
Save a VisIt root file for the collection.
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1671
const char * VTKByteOrder()
Determine the byte order and return either &quot;BigEndian&quot; or &quot;LittleEndian&quot;.
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...
int GetNE() const
Return the number of entities.
Definition: qspace.hpp:60
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:2781
virtual void Load(int cycle_=0)
Load the collection. Not implemented in the base class DataCollection.
void SetLevelsOfDetail(int levels_of_detail_)
std::string GetMeshFileName() const
Vector data type.
Definition: vector.hpp:60
Mesh * GetMesh() const
Returns the mesh.
Definition: qspace.hpp:63
void LoadVisItRootFile(const std::string &root_name)
virtual void SetMesh(Mesh *new_mesh)
Set/change the mesh associated with the collection.
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)
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
int num_procs
Number of MPI ranks (in parallel)
void ParPrint(std::ostream &out) const
Save the mesh in a parallel mesh format.
Definition: pmesh.cpp:6238
Class for parallel meshes.
Definition: pmesh.hpp:32
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:23
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.