MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
datacollection.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
30namespace mfem
31{
32
33// static method
34int 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
75DataCollection::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;
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;
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;
136 }
137#endif
138}
139
140#ifdef MFEM_USE_MPI
141void 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
175void 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
191void 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
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 {
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 {
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
261std::string DataCollection::GetFieldFileName(const std::string &field_name)
262const
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 {
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 {
299 MFEM_WARNING("Error writing q-field to file: " << it->first);
300 }
301}
302
303void 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
312void 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();
336}
337
339{
340 DeleteData();
341}
342
343
344// class VisItDataCollection implementation
345
347{
348 if (mesh)
349 {
352 if (mesh->NURBSext)
353 {
356 }
357 }
358 else
359 {
360 spatial_dim = 0;
361 topo_dim = 0;
362 }
363}
364
365VisItDataCollection::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
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
394}
395#endif
396
398{
399 DataCollection::SetMesh(new_mesh);
402}
403
404#ifdef MFEM_USE_MPI
405void 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
415void 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
436void VisItDataCollection::RegisterQField(const std::string& name,
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
456{
457 visit_levels_of_detail = levels_of_detail;
458}
459
460void 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 {
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.");
508#else
509 if (m_comm == MPI_COMM_NULL)
510 {
511 MFEM_WARNING("Cannot load parallel VisIt root file without MPI"
512 " communicator");
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);
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
548void 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 {
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.
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 {
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
590 MFEM_WARNING("Reading parallel format in serial is not supported");
591 return;
592#endif
593 }
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
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 {
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
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
699void 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 {
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 {
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
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 {
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
1036void 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
1062void 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 {
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 {
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
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
int cycle
Time cycle; for time-dependent simulations cycle >= 0, otherwise = -1.
real_t time
Physical time (for time-dependent simulations)
QFieldMap::iterator QFieldMapIterator
virtual void SetMesh(Mesh *new_mesh)
Set/change the mesh associated with the collection.
void SaveOneQField(const QFieldMapIterator &it)
Save one q-field to disk, assuming the collection directory exists.
virtual void RegisterQField(const std::string &q_field_name, QuadratureFunction *qf)
Add a QuadratureFunction to the collection.
const std::string & GetCollectionName() const
Get the name of the collection.
bool own_data
Should the collection delete its mesh and fields.
DataCollection(const std::string &collection_name, Mesh *mesh_=NULL)
Initialize the collection with its name and Mesh.
static int create_directory(const std::string &dir_name, const Mesh *mesh, int myid)
int GetCycle() const
Get time cycle (for time-dependent simulations)
GFieldMap::iterator FieldMapIterator
void SaveOneField(const FieldMapIterator &it)
Save one field to disk, assuming the collection directory exists.
void DeleteAll()
Delete data owned by the DataCollection including field information.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
static const int precision_default
Default value for precision.
int pad_digits_cycle
Number of digits used for the cycle and MPI rank in filenames.
bool serial
Serial or parallel run? False iff mesh is a ParMesh.
virtual void SetCompression(bool comp)
Set the flag for use of gz compressed files.
std::string prefix_path
A path where the directory with results is saved. If not empty, it has '/' at the end.
std::string GetFieldFileName(const std::string &field_name) const
int myid
MPI rank (in parallel)
static const int pad_digits_default
Default value for pad_digits_*.
real_t time_step
Time step i.e. delta_t (for time-dependent simulations)
int num_procs
Number of MPI ranks (in parallel)
int error
Error state.
std::string GetMeshShortFileName() const
virtual void Load(int cycle_=0)
Load the collection. Not implemented in the base class DataCollection.
virtual void SetFormat(int fmt)
Set the desired output mesh and data format.
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
virtual ~DataCollection()
Delete the mesh and fields if owned by the collection.
std::string GetMeshFileName() const
void DeleteData()
Delete data owned by the DataCollection keeping field information.
bool appendRankToFileName
Append rank to any output file names.
std::string name
Name of the collection, used as a directory name when saving.
virtual void Save()
Save the collection to disk.
virtual void SaveField(const std::string &field_name)
Save one field, assuming the collection directory already exists.
MPI_Comm m_comm
Associated MPI communicator.
Mesh * mesh
The (common) mesh for the collected fields.
virtual void SaveMesh()
Save the mesh, creating the collection directory.
int precision
Precision (number of digits) used for the text output of doubles.
int format
Output mesh format: see the Format enumeration.
virtual void SaveQField(const std::string &q_field_name)
Save one q-field, assuming the collection directory already exists.
real_t GetTime() const
Get physical time (for time-dependent simulations)
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:561
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition: fespace.cpp:3168
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:333
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1136
static int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition: geom.cpp:1849
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:696
int VectorDim() const
Definition: gridfunc.cpp:323
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:256
Mesh data type.
Definition: mesh.hpp:56
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:290
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition: mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1163
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:11695
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1385
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:489
void Register(const std::string &fname, T *field, bool own_data)
Register field field with name fname.
iterator end()
Returns an end iterator to the registered fields.
iterator begin()
Returns a begin iterator to the registered fields.
void DeleteData(bool own_data)
Clear all associations between names and fields.
void clear()
Clears the map of registered fields without reclaiming memory.
iterator find(const std::string &fname)
Returns an iterator to the field fname.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
Class for parallel grid function.
Definition: pgridfunc.hpp:33
Class for parallel meshes.
Definition: pmesh.hpp:34
MPI_Comm GetComm() const
Definition: pmesh.hpp:402
int GetMyRank() const
Definition: pmesh.hpp:404
void ParPrint(std::ostream &out, const std::string &comments="") const
Definition: pmesh.cpp:6314
int GetNRanks() const
Definition: pmesh.hpp:403
void WritePVTUHeader(std::ostream &out)
void SetHighOrderOutput(bool high_order_output_)
void UseRestartMode(bool restart_mode_)
void SetCompression(bool compression_) override
virtual void Load(int cycle_=0) override
Load the collection - not implemented in the ParaView writer.
bool IsBinaryFormat() const
Returns true if the output format is BINARY or BINARY32, false if ASCII.
ParaViewDataCollection(const std::string &collection_name, mfem::Mesh *mesh_=NULL)
Constructor. The collection name is used when saving the data.
void SaveGFieldVTU(std::ostream &out, int ref_, const FieldMapIterator &it)
const char * GetDataFormatString() const
void WritePVTUFooter(std::ostream &out, const std::string &vtu_prefix)
void SaveDataVTU(std::ostream &out, int ref)
void SetCompressionLevel(int compression_level_)
Set the zlib compression level.
void SetLevelsOfDetail(int levels_of_detail_)
int GetCompressionLevel() const
If compression is enabled, return the compression level, otherwise return 0.
std::string GenerateVTUFileName(const std::string &prefix, int rank)
virtual void Save() override
const char * GetDataTypeString() const
std::string GeneratePVTUFileName(const std::string &prefix)
void SetDataFormat(VTKFormat fmt)
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:24
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:82
const IntegrationRule & GetIntRule(int idx) const
Get the IntegrationRule associated with entity (element or face) idx.
Definition: qfunction.hpp:176
Mesh * GetMesh() const
Returns the mesh.
Definition: qspace.hpp:72
IntegrationRule RefPts
Definition: geom.hpp:317
Vector data type.
Definition: vector.hpp:80
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
void SaveRootFile()
Save a VisIt root file for the collection.
std::string GetVisItRootString()
Prepare the VisIt root file in JSON format for the current collection.
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 LoadVisItRootFile(const std::string &root_name)
VisItDataCollection(const std::string &collection_name, Mesh *mesh_=NULL)
Constructor. The collection name is used when saving the data.
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.
virtual void Save() override
Save the collection and a VisIt root file.
void SetMaxLevelsOfDetail(int max_levels_of_detail)
Set VisIt parameter: maximum levels of detail for the MultiresControl.
std::map< std::string, VisItFieldInfo >::iterator FieldInfoMapIterator
virtual void SetMesh(Mesh *new_mesh) override
Set/change the mesh associated with the collection.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
std::map< std::string, VisItFieldInfo > field_info_map
void DeleteAll()
Delete all data owned by VisItDataCollection including field data information.
void ParseVisItRootString(const std::string &json)
Read in a VisIt root file in JSON format.
Helper class for VisIt visualization data.
int main()
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
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1891
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
VTKFormat
Data array format for VTK and VTU files.
Definition: vtk.hpp:99
@ ASCII
Data arrays will be written in ASCII format.
int to_int(const std::string &str)
Convert a string to an int.
Definition: text.hpp:62
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
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
Definition: vtk.cpp:602