libMesh
equation_systems_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 #include "libmesh/libmesh_common.h"
20 #include "libmesh/libmesh_logging.h"
21 
22 
23 // Local Includes
24 #include "libmesh/libmesh_version.h"
25 #include "libmesh/equation_systems.h"
26 #include "libmesh/mesh_base.h"
27 #include "libmesh/mesh_refinement.h"
28 #include "libmesh/mesh_tools.h"
29 #include "libmesh/parallel.h"
30 #include "libmesh/utility.h"
31 #include "libmesh/xdr_cxx.h"
32 
33 // C++ Includes
34 #include <iomanip> // setfill
35 #include <sstream>
36 #include <string>
37 
38 namespace libMesh
39 {
40 
41 // Forward Declarations
42 
43 // Anonymous namespace for implementation details.
44 namespace {
45 std::string local_file_name (const unsigned int processor_id,
46  std::string_view basename)
47 {
48  std::ostringstream returnval;
49 
50  std::string_view suffix;
51  if (Utility::ends_with(basename, ".bz2"))
52  {
53  basename.remove_suffix(4);
54  suffix = ".bz2";
55  }
56  else if (Utility::ends_with(basename, ".gz"))
57  {
58  basename.remove_suffix(3);
59  suffix = ".gz";
60  }
61 
62  returnval << basename << '.';
63  returnval << std::setfill('0') << std::setw(4);
64  returnval << processor_id;
65  returnval << suffix;
66 
67  return returnval.str();
68 }
69 }
70 
71 
72 
73 
74 // ------------------------------------------------------------
75 // EquationSystem class implementation
76 template <typename InValType>
77 void EquationSystems::read (std::string_view name,
78  const unsigned int read_flags,
79  bool partition_agnostic)
80 {
81  XdrMODE mode = READ;
82  if (name.find(".xdr") != std::string::npos)
83  mode = DECODE;
84  this->read(name, mode, read_flags, partition_agnostic);
85 }
86 
87 
88 
89 template <typename InValType>
90 void EquationSystems::read (std::string_view name,
91  const XdrMODE mode,
92  const unsigned int read_flags,
93  bool partition_agnostic)
94 {
95  // This will unzip a file with .bz2 as the extension, otherwise it
96  // simply returns the name if the file need not be unzipped.
97  Xdr io ((this->processor_id() == 0) ? std::string(name) : "", mode);
98 
99  std::function<std::unique_ptr<Xdr>()> local_io_functor;
100  local_io_functor = [this,&name,&mode]() {
101  return std::make_unique<Xdr>(local_file_name(this->processor_id(), name), mode); };
102 
103  this->read(io, local_io_functor, read_flags, partition_agnostic);
104 }
105 
106 
107 
108 template <typename InValType>
110  std::function<std::unique_ptr<Xdr>()> & local_io_functor,
111  const unsigned int read_flags,
112  bool partition_agnostic)
113 {
177  // Set booleans from the read_flags argument
178  const bool read_header = read_flags & EquationSystems::READ_HEADER;
179  const bool read_data = read_flags & EquationSystems::READ_DATA;
180  const bool read_additional_data = read_flags & EquationSystems::READ_ADDITIONAL_DATA;
181  const bool read_legacy_format = read_flags & EquationSystems::READ_LEGACY_FORMAT;
182  const bool try_read_ifems = read_flags & EquationSystems::TRY_READ_IFEMS;
183  const bool read_basic_only = read_flags & EquationSystems::READ_BASIC_ONLY;
184  bool read_parallel_files = false;
185 
186  std::vector<std::pair<std::string, System *>> xda_systems;
187 
188  libmesh_assert (io.reading());
189 
190  {
191  // 1.)
192  // Read the version header.
193  std::string version = "legacy";
194  if (!read_legacy_format)
195  {
196  if (this->processor_id() == 0) io.data(version);
197  this->comm().broadcast(version);
198 
199  // All processors have the version header, if it does not contain
200  // the libMesh_label string then it is a legacy file.
201  const std::string libMesh_label = "libMesh-";
202  std::string::size_type lm_pos = version.find(libMesh_label);
203  if (lm_pos==std::string::npos)
204  {
205  io.close();
206 
207  // Recursively call this read() function but with the
208  // EquationSystems::READ_LEGACY_FORMAT bit set.
209  this->read (io, local_io_functor, (read_flags | EquationSystems::READ_LEGACY_FORMAT), partition_agnostic);
210  return;
211  }
212 
213  // Figure out the libMesh version that created this file
214  std::istringstream iss(version.substr(lm_pos + libMesh_label.size()));
215  int ver_major = 0, ver_minor = 0, ver_patch = 0;
216  char dot;
217  iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
218  io.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));
219 
220 
221  read_parallel_files = Utility::contains(version, " parallel");
222 
223  // If requested that we try to read infinite element information,
224  // and the string " with infinite elements" is not in the version,
225  // then tack it on. This is for compatibility reading ifem
226  // files written prior to 11/10/2008 - BSK
227  if (try_read_ifems)
228  if (!Utility::contains(version, " with infinite elements"))
229  version += " with infinite elements";
230 
231  }
232  else
233  libmesh_deprecated();
234 
235  LOG_SCOPE("read()", "EquationSystems");
236 
237  // 2.)
238  // Read the number of equation systems
239  unsigned int n_sys=0;
240  if (this->processor_id() == 0) io.data (n_sys);
241  this->comm().broadcast(n_sys);
242 
243  for (unsigned int sys=0; sys<n_sys; sys++)
244  {
245  // 3.)
246  // Read the name of the sys-th equation system
247  std::string sys_name;
248  if (this->processor_id() == 0) io.data (sys_name);
249  this->comm().broadcast(sys_name);
250 
251  // 4.)
252  // Read the type of the sys-th equation system
253  std::string sys_type;
254  if (this->processor_id() == 0) io.data (sys_type);
255  this->comm().broadcast(sys_type);
256 
257  if (read_header)
258  this->add_system (sys_type, sys_name);
259 
260  // 5.) - 9.)
261  // Let System::read_header() do the job
262  System & new_system = this->get_system(sys_name);
263  new_system.read_header (io,
264  version,
265  read_header,
266  read_additional_data,
267  read_legacy_format);
268 
269  xda_systems.emplace_back(sys_name, &new_system);
270 
271  // If we're only creating "basic" systems, we need to tell
272  // each system that before we call init() later.
273  if (read_basic_only)
274  new_system.set_basic_system_only();
275  }
276  }
277 
278 
279 
280  // Now we are ready to initialize the underlying data
281  // structures. This will initialize the vectors for
282  // storage, the dof_map, etc...
283  if (read_header)
284  this->init();
285 
286  // 10.) & 11.)
287  // Read and set the numeric vector values
288  if (read_data)
289  {
290  std::unique_ptr<Xdr> local_io;
291 
292  // the EquationSystems::read() method should look constant from the mesh
293  // perspective, but we need to assign a temporary numbering to the nodes
294  // and elements in the mesh, which requires that we abuse const_cast
295  if (!read_legacy_format && partition_agnostic)
296  {
297  MeshBase & mesh = const_cast<MeshBase &>(this->get_mesh());
299  }
300 
301  for (auto & pr : xda_systems)
302  if (read_legacy_format)
303  {
304  libmesh_deprecated();
305 #ifdef LIBMESH_ENABLE_DEPRECATED
306  pr.second->read_legacy_data (io, read_additional_data);
307 #endif
308  }
309  else
310  if (read_parallel_files)
311  {
312  if (!local_io)
313  {
314  local_io = local_io_functor();
315  libmesh_assert(local_io->reading());
316  }
317  pr.second->read_parallel_data<InValType> (*local_io, read_additional_data);
318  }
319  else
320  pr.second->read_serialized_data<InValType> (io, read_additional_data);
321 
322 
323  // Undo the temporary numbering.
324  if (!read_legacy_format && partition_agnostic)
326  }
327 
328  // Localize each system's data
329  this->update();
330 
331  #ifdef LIBMESH_ENABLE_AMR
332  MeshRefinement mesh_refine(_mesh);
333  mesh_refine.clean_refinement_flags();
334  #endif
335 }
336 
337 
338 
339 void EquationSystems::write(std::string_view name,
340  const unsigned int write_flags,
341  bool partition_agnostic) const
342 {
343  XdrMODE mode = WRITE;
344  if (name.find(".xdr") != std::string::npos)
345  mode = ENCODE;
346  this->write(name, mode, write_flags, partition_agnostic);
347 }
348 
349 
350 
351 void EquationSystems::write(std::string_view name,
352  const XdrMODE mode,
353  const unsigned int write_flags,
354  bool partition_agnostic) const
355 {
356  Xdr io((this->processor_id()==0) ? std::string(name) : "", mode);
357 
358  std::unique_ptr<Xdr> local_io;
359  // open a parallel buffer if warranted
360  if (write_flags & EquationSystems::WRITE_PARALLEL_FILES && write_flags & EquationSystems::WRITE_DATA)
361  local_io = std::make_unique<Xdr>(local_file_name(this->processor_id(),name), mode);
362 
363  this->write(io, write_flags, partition_agnostic, local_io.get());
364 }
365 
366 
367 
369  const unsigned int write_flags,
370  bool partition_agnostic,
371  Xdr * const local_io) const
372 {
436  // the EquationSystems::write() method should look constant,
437  // but we need to assign a temporary numbering to the nodes
438  // and elements in the mesh, which requires that we abuse const_cast
439  if (partition_agnostic)
440  {
441  MeshBase & mesh = const_cast<MeshBase &>(this->get_mesh());
443  }
444 
445  // set booleans from write_flags argument
446  const bool write_data = write_flags & EquationSystems::WRITE_DATA;
447  const bool write_additional_data = write_flags & EquationSystems::WRITE_ADDITIONAL_DATA;
448 
449  // always write parallel files if we're instructed to write in
450  // parallel
451  const bool write_parallel_files =
453  // Even if we're on a distributed mesh, we may or may not have a
454  // consistent way of reconstructing the same mesh partitioning
455  // later, but we need the same mesh partitioning if we want to
456  // reread the parallel solution safely, so let's write a serial file
457  // unless specifically requested not to.
458  // ||
459  // // but also write parallel files if we haven't been instructed to
460  // // write in serial and we're on a distributed mesh
461  // (!(write_flags & EquationSystems::WRITE_SERIAL_FILES) &&
462  // !this->get_mesh().is_serial())
463  ;
464 
465  if (write_parallel_files && write_data)
466  libmesh_assert(local_io);
467 
468  {
469  libmesh_assert (io.writing());
470 
471  LOG_SCOPE("write()", "EquationSystems");
472 
473  const unsigned int proc_id = this->processor_id();
474 
475  unsigned int n_sys = 0;
476  for (auto & pr : _systems)
477  if (!pr.second->hide_output())
478  n_sys++;
479 
480  // set the version number in the Xdr object
481  io.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
482  LIBMESH_MINOR_VERSION,
483  LIBMESH_MICRO_VERSION));
484 
485  // Only write the header information
486  // if we are processor 0.
487  if (proc_id == 0)
488  {
489  std::string comment;
490 
491  // 1.)
492  // Write the version header
493  std::string version("libMesh-" + libMesh::get_io_compatibility_version());
494  if (write_parallel_files) version += " parallel";
495 
496 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
497  version += " with infinite elements";
498 #endif
499  io.data (version, "# File Format Identifier");
500 
501  // 2.)
502  // Write the number of equation systems
503  io.data (n_sys, "# No. of Equation Systems");
504 
505  for (auto & [sys_name, sys] : _systems)
506  {
507  // Ignore this system if it has been marked as hidden
508  if (sys->hide_output()) continue;
509 
510  // 3.)
511  // Write the name of the sys_num-th system
512  {
513  const unsigned int sys_num = sys->number();
514 
515  comment = "# Name, System No. ";
516  comment += std::to_string(sys_num);
517 
518  // Note: There is no Xdr::data overload taking a "const
519  // std::string &" so we need to make a copy.
520  std::string copy = sys_name;
521  io.data (copy, comment);
522  }
523 
524  // 4.)
525  // Write the type of system handled
526  {
527  const unsigned int sys_num = sys->number();
528  std::string sys_type = sys->system_type();
529 
530  comment = "# Type, System No. ";
531  comment += std::to_string(sys_num);
532 
533  io.data (sys_type, comment);
534  }
535 
536  // 5.) - 9.)
537  // Let System::write_header() do the job
538  sys->write_header (io, version, write_additional_data);
539  }
540  }
541 
542  // Start from the first system, again,
543  // to write vectors to disk, if wanted
544  if (write_data)
545  {
546  for (auto & pr : _systems)
547  {
548  // Ignore this system if it has been marked as hidden
549  if (pr.second->hide_output()) continue;
550 
551  // 10.) + 11.)
552  if (write_parallel_files)
553  pr.second->write_parallel_data (*local_io,write_additional_data);
554  else
555  pr.second->write_serialized_data (io,write_additional_data);
556  }
557 
558  if (local_io)
559  local_io->close();
560  }
561 
562  io.close();
563  }
564 
565  // the EquationSystems::write() method should look constant,
566  // but we need to undo the temporary numbering of the nodes
567  // and elements in the mesh, which requires that we abuse const_cast
568  if (partition_agnostic)
569  const_cast<MeshBase &>(_mesh).fix_broken_node_and_element_numbering();
570 }
571 
572 
573 
574 // template specialization
575 
576 template LIBMESH_EXPORT void EquationSystems::read<Number> (Xdr & io, std::function<std::unique_ptr<Xdr>()> & local_io_functor, const unsigned int read_flags, bool partition_agnostic);
577 template LIBMESH_EXPORT void EquationSystems::read<Number> (std::string_view name, const unsigned int read_flags, bool partition_agnostic);
578 template LIBMESH_EXPORT void EquationSystems::read<Number> (std::string_view name, const XdrMODE mode, const unsigned int read_flags, bool partition_agnostic);
579 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
580 template LIBMESH_EXPORT void EquationSystems::read<Real> (Xdr & io, std::function<std::unique_ptr<Xdr>()> & local_io_functor, const unsigned int read_flags, bool partition_agnostic);
581 template LIBMESH_EXPORT void EquationSystems::read<Real> (std::string_view name, const unsigned int read_flags, bool partition_agnostic);
582 template LIBMESH_EXPORT void EquationSystems::read<Real> (std::string_view name, const XdrMODE mode, const unsigned int read_flags, bool partition_agnostic);
583 #endif
584 
585 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
bool writing() const
Definition: xdr_cxx.h:129
void write(std::string_view name, const XdrMODE, const unsigned int write_flags=(WRITE_DATA), bool partition_agnostic=true) const
Write the systems to disk using the XDR data format.
bool ends_with(std::string_view superstring, std::string_view suffix)
Look for a substring at the very end of a string.
Definition: utility.C:213
void close()
Closes the file if it is open.
Definition: xdr_cxx.C:277
MeshBase & mesh
const Parallel::Communicator & comm() const
void clean_refinement_flags()
Sets the refinement flag to Elem::DO_NOTHING for each element in the mesh.
virtual void fix_broken_node_and_element_numbering()=0
There is no reason for a user to ever call this function.
The libMesh namespace provides an interface to certain functionality in the library.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:75
MeshBase & _mesh
The mesh data structure.
XdrMODE
Defines an enum for read/write mode in Xdr format.
Definition: enum_xdr_mode.h:35
void data(T &a, std::string_view comment="")
Inputs or outputs a single value.
Definition: xdr_cxx.C:778
void set_basic_system_only()
Sets the system to be "basic only": i.e.
Definition: system.h:2422
Implements (adaptive) mesh refinement algorithms for a MeshBase.
std::string get_io_compatibility_version()
Specifier for I/O file compatibility features.
void read_header(Xdr &io, std::string_view version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
Reads the basic data header for this System.
Definition: system_io.C:97
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
libmesh_assert(ctx)
void read(std::string_view name, const XdrMODE, const unsigned int read_flags=(READ_HEADER|READ_DATA), bool partition_agnostic=true)
Read & initialize the systems from disk using the XDR data format.
bool contains(std::string_view superstring, std::string_view substring)
Look for a substring within a string.
Definition: utility.C:205
void globally_renumber_nodes_and_elements(MeshBase &)
There is no reason for a user to ever call this function.
Definition: mesh_tools.C:2542
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:67
void set_version(int ver)
Sets the version of the file that is being read.
Definition: xdr_cxx.h:171
const MeshBase & get_mesh() const
bool reading() const
Definition: xdr_cxx.h:123
virtual void init()
Initialize all the systems.
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
processor_id_type processor_id() const
void update()
Updates local values for all the systems.
std::map< std::string, std::unique_ptr< System >, std::less<> > _systems
Data structure holding the systems.