Line data Source code
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 0 : std::string local_file_name (const unsigned int processor_id,
46 : std::string_view basename)
47 : {
48 0 : std::ostringstream returnval;
49 :
50 0 : std::string_view suffix;
51 0 : if (Utility::ends_with(basename, ".bz2"))
52 : {
53 0 : basename.remove_suffix(4);
54 0 : suffix = ".bz2";
55 : }
56 0 : else if (Utility::ends_with(basename, ".gz"))
57 : {
58 0 : basename.remove_suffix(3);
59 0 : suffix = ".gz";
60 : }
61 :
62 0 : returnval << basename << '.';
63 0 : returnval << std::setfill('0') << std::setw(4);
64 0 : returnval << processor_id;
65 0 : returnval << suffix;
66 :
67 0 : return returnval.str();
68 0 : }
69 : }
70 :
71 :
72 :
73 :
74 : // ------------------------------------------------------------
75 : // EquationSystem class implementation
76 : template <typename InValType>
77 71 : void EquationSystems::read (std::string_view name,
78 : const unsigned int read_flags,
79 : bool partition_agnostic)
80 : {
81 2 : XdrMODE mode = READ;
82 71 : if (name.find(".xdr") != std::string::npos)
83 0 : mode = DECODE;
84 71 : this->read(name, mode, read_flags, partition_agnostic);
85 71 : }
86 :
87 :
88 :
89 : template <typename InValType>
90 10578 : 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 23887 : Xdr io ((this->processor_id() == 0) ? std::string(name) : "", mode);
98 :
99 608 : std::function<std::unique_ptr<Xdr>()> local_io_functor;
100 10578 : local_io_functor = [this,&name,&mode]() {
101 0 : return std::make_unique<Xdr>(local_file_name(this->processor_id(), name), mode); };
102 :
103 10578 : this->read(io, local_io_functor, read_flags, partition_agnostic);
104 10578 : }
105 :
106 :
107 :
108 : template <typename InValType>
109 10578 : void EquationSystems::read (Xdr & io,
110 : std::function<std::unique_ptr<Xdr>()> & local_io_functor,
111 : const unsigned int read_flags,
112 : bool partition_agnostic)
113 : {
114 : /**
115 : * This program implements the output of an
116 : * EquationSystems object. This warrants some
117 : * documentation. The output file essentially
118 : * consists of 11 sections:
119 : \verbatim
120 : 1.) A version header (for non-'legacy' formats, libMesh-0.7.0 and greater).
121 : 2.) The number of individual equation systems (unsigned int)
122 :
123 : for each system
124 :
125 : 3.) The name of the system (string)
126 : 4.) The type of the system (string)
127 :
128 : handled through System::read():
129 :
130 : +-------------------------------------------------------------+
131 : | 5.) The number of variables in the system (unsigned int) |
132 : | |
133 : | for each variable in the system |
134 : | |
135 : | 6.) The name of the variable (string) |
136 : | |
137 : | 7.) Combined in an FEType: |
138 : | - The approximation order(s) of the variable (Order |
139 : | Enum, cast to int/s) |
140 : | - The finite element family/ies of the variable |
141 : | (FEFamily Enum, cast to int/s) |
142 : | |
143 : | end variable loop |
144 : | |
145 : | 8.) The number of additional vectors (unsigned int), |
146 : | |
147 : | for each additional vector in the equation system object |
148 : | |
149 : | 9.) the name of the additional vector (string) |
150 : +-------------------------------------------------------------+
151 :
152 : end system loop
153 :
154 :
155 : for each system, handled through System::read_{serialized,parallel}_data():
156 :
157 : +--------------------------------------------------------------+
158 : | 10.) The global solution vector, re-ordered to be node-major |
159 : | (More on this later.) |
160 : | |
161 : | for each additional vector in the equation system object |
162 : | |
163 : | 11.) The global additional vector, re-ordered to be |
164 : | node-major (More on this later.) |
165 : +--------------------------------------------------------------+
166 :
167 : end system loop
168 : \endverbatim
169 : *
170 : * Note that the actual IO is handled through the Xdr class
171 : * (to be renamed later?) which provides a uniform interface to
172 : * both the XDR (eXternal Data Representation) interface and standard
173 : * ASCII output. Thus this one section of code will read XDR or ASCII
174 : * files with no changes.
175 : */
176 :
177 : // Set booleans from the read_flags argument
178 10578 : const bool read_header = read_flags & EquationSystems::READ_HEADER;
179 10578 : const bool read_data = read_flags & EquationSystems::READ_DATA;
180 10578 : const bool read_additional_data = read_flags & EquationSystems::READ_ADDITIONAL_DATA;
181 10578 : const bool read_legacy_format = read_flags & EquationSystems::READ_LEGACY_FORMAT;
182 10578 : const bool try_read_ifems = read_flags & EquationSystems::TRY_READ_IFEMS;
183 10578 : const bool read_basic_only = read_flags & EquationSystems::READ_BASIC_ONLY;
184 304 : bool read_parallel_files = false;
185 :
186 608 : std::vector<std::pair<std::string, System *>> xda_systems;
187 :
188 304 : libmesh_assert (io.reading());
189 :
190 : {
191 : // 1.)
192 : // Read the version header.
193 10578 : std::string version = "legacy";
194 10578 : if (!read_legacy_format)
195 : {
196 11034 : if (this->processor_id() == 0) io.data(version);
197 10578 : 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 10578 : const std::string libMesh_label = "libMesh-";
202 304 : std::string::size_type lm_pos = version.find(libMesh_label);
203 10578 : if (lm_pos==std::string::npos)
204 : {
205 0 : io.close();
206 :
207 : // Recursively call this read() function but with the
208 : // EquationSystems::READ_LEGACY_FORMAT bit set.
209 0 : this->read (io, local_io_functor, (read_flags | EquationSystems::READ_LEGACY_FORMAT), partition_agnostic);
210 0 : return;
211 : }
212 :
213 : // Figure out the libMesh version that created this file
214 11186 : std::istringstream iss(version.substr(lm_pos + libMesh_label.size()));
215 10578 : int ver_major = 0, ver_minor = 0, ver_patch = 0;
216 : char dot;
217 10578 : iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
218 10578 : io.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));
219 :
220 :
221 10578 : 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 10578 : if (try_read_ifems)
228 0 : if (!Utility::contains(version, " with infinite elements"))
229 0 : version += " with infinite elements";
230 :
231 9970 : }
232 : else
233 : libmesh_deprecated();
234 :
235 608 : LOG_SCOPE("read()", "EquationSystems");
236 :
237 : // 2.)
238 : // Read the number of equation systems
239 10578 : unsigned int n_sys=0;
240 11034 : if (this->processor_id() == 0) io.data (n_sys);
241 10578 : this->comm().broadcast(n_sys);
242 :
243 21507 : for (unsigned int sys=0; sys<n_sys; sys++)
244 : {
245 : // 3.)
246 : // Read the name of the sys-th equation system
247 628 : std::string sys_name;
248 11400 : if (this->processor_id() == 0) io.data (sys_name);
249 10929 : this->comm().broadcast(sys_name);
250 :
251 : // 4.)
252 : // Read the type of the sys-th equation system
253 628 : std::string sys_type;
254 11400 : if (this->processor_id() == 0) io.data (sys_type);
255 10929 : this->comm().broadcast(sys_type);
256 :
257 10929 : if (read_header)
258 875 : this->add_system (sys_type, sys_name);
259 :
260 : // 5.) - 9.)
261 : // Let System::read_header() do the job
262 10929 : System & new_system = this->get_system(sys_name);
263 11243 : new_system.read_header (io,
264 : version,
265 : read_header,
266 : read_additional_data,
267 : read_legacy_format);
268 :
269 10929 : 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 10929 : if (read_basic_only)
274 0 : 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 10578 : if (read_header)
284 498 : this->init();
285 :
286 : // 10.) & 11.)
287 : // Read and set the numeric vector values
288 10578 : if (read_data)
289 : {
290 10274 : 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 10578 : if (!read_legacy_format && partition_agnostic)
296 : {
297 608 : MeshBase & mesh = const_cast<MeshBase &>(this->get_mesh());
298 10578 : MeshTools::Private::globally_renumber_nodes_and_elements(mesh);
299 : }
300 :
301 21507 : for (auto & pr : xda_systems)
302 : {
303 10929 : libmesh_error_msg_if(read_legacy_format,
304 : "Reading legacy format XDR files is officially no longer supported.");
305 :
306 10929 : if (read_parallel_files)
307 : {
308 0 : if (!local_io)
309 : {
310 0 : local_io = local_io_functor();
311 0 : libmesh_assert(local_io->reading());
312 : }
313 0 : pr.second->read_parallel_data<InValType> (*local_io, read_additional_data);
314 : }
315 : else
316 10929 : pr.second->read_serialized_data<InValType> (io, read_additional_data);
317 : }
318 :
319 : // Undo the temporary numbering.
320 10578 : if (!read_legacy_format && partition_agnostic)
321 10578 : _mesh.fix_broken_node_and_element_numbering();
322 9970 : }
323 :
324 : // Localize each system's data
325 10578 : this->update();
326 :
327 : #ifdef LIBMESH_ENABLE_AMR
328 11186 : MeshRefinement mesh_refine(_mesh);
329 10578 : mesh_refine.clean_refinement_flags();
330 : #endif
331 9970 : }
332 :
333 :
334 :
335 71 : void EquationSystems::write(std::string_view name,
336 : const unsigned int write_flags,
337 : bool partition_agnostic) const
338 : {
339 2 : XdrMODE mode = WRITE;
340 71 : if (name.find(".xdr") != std::string::npos)
341 0 : mode = ENCODE;
342 71 : this->write(name, mode, write_flags, partition_agnostic);
343 71 : }
344 :
345 :
346 :
347 8344 : void EquationSystems::write(std::string_view name,
348 : const XdrMODE mode,
349 : const unsigned int write_flags,
350 : bool partition_agnostic) const
351 : {
352 10752 : Xdr io((this->processor_id()==0) ? std::string(name) : "", mode);
353 :
354 8102 : std::unique_ptr<Xdr> local_io;
355 : // open a parallel buffer if warranted
356 8344 : if (write_flags & EquationSystems::WRITE_PARALLEL_FILES && write_flags & EquationSystems::WRITE_DATA)
357 0 : local_io = std::make_unique<Xdr>(local_file_name(this->processor_id(),name), mode);
358 :
359 8344 : this->write(io, write_flags, partition_agnostic, local_io.get());
360 8344 : }
361 :
362 :
363 :
364 8344 : void EquationSystems::write(Xdr & io,
365 : const unsigned int write_flags,
366 : bool partition_agnostic,
367 : Xdr * const local_io) const
368 : {
369 : /**
370 : * This program implements the output of an
371 : * EquationSystems object. This warrants some
372 : * documentation. The output file essentially
373 : * consists of 11 sections:
374 : \verbatim
375 : 1.) The version header.
376 : 2.) The number of individual equation systems (unsigned int)
377 :
378 : for each system
379 :
380 : 3.) The name of the system (string)
381 : 4.) The type of the system (string)
382 :
383 : handled through System::read():
384 :
385 : +-------------------------------------------------------------+
386 : | 5.) The number of variables in the system (unsigned int) |
387 : | |
388 : | for each variable in the system |
389 : | |
390 : | 6.) The name of the variable (string) |
391 : | |
392 : | 7.) Combined in an FEType: |
393 : | - The approximation order(s) of the variable (Order |
394 : | Enum, cast to int/s) |
395 : | - The finite element family/ies of the variable |
396 : | (FEFamily Enum, cast to int/s) |
397 : | |
398 : | end variable loop |
399 : | |
400 : | 8.) The number of additional vectors (unsigned int), |
401 : | |
402 : | for each additional vector in the equation system object |
403 : | |
404 : | 9.) the name of the additional vector (string) |
405 : +-------------------------------------------------------------+
406 :
407 : end system loop
408 :
409 :
410 : for each system, handled through System::write_{serialized,parallel}_data():
411 :
412 : +--------------------------------------------------------------+
413 : | 10.) The global solution vector, re-ordered to be node-major |
414 : | (More on this later.) |
415 : | |
416 : | for each additional vector in the equation system object |
417 : | |
418 : | 11.) The global additional vector, re-ordered to be |
419 : | node-major (More on this later.) |
420 : +--------------------------------------------------------------+
421 :
422 : end system loop
423 : \endverbatim
424 : *
425 : * Note that the actual IO is handled through the Xdr class
426 : * (to be renamed later?) which provides a uniform interface to
427 : * both the XDR (eXternal Data Representation) interface and standard
428 : * ASCII output. Thus this one section of code will write XDR or ASCII
429 : * files with no changes.
430 : */
431 :
432 : // the EquationSystems::write() method should look constant,
433 : // but we need to assign a temporary numbering to the nodes
434 : // and elements in the mesh, which requires that we abuse const_cast
435 8344 : if (partition_agnostic)
436 : {
437 484 : MeshBase & mesh = const_cast<MeshBase &>(this->get_mesh());
438 8344 : MeshTools::Private::globally_renumber_nodes_and_elements(mesh);
439 : }
440 :
441 : // set booleans from write_flags argument
442 8344 : const bool write_data = write_flags & EquationSystems::WRITE_DATA;
443 8344 : const bool write_additional_data = write_flags & EquationSystems::WRITE_ADDITIONAL_DATA;
444 :
445 : // always write parallel files if we're instructed to write in
446 : // parallel
447 242 : const bool write_parallel_files =
448 8344 : (write_flags & EquationSystems::WRITE_PARALLEL_FILES)
449 : // Even if we're on a distributed mesh, we may or may not have a
450 : // consistent way of reconstructing the same mesh partitioning
451 : // later, but we need the same mesh partitioning if we want to
452 : // reread the parallel solution safely, so let's write a serial file
453 : // unless specifically requested not to.
454 : // ||
455 : // // but also write parallel files if we haven't been instructed to
456 : // // write in serial and we're on a distributed mesh
457 : // (!(write_flags & EquationSystems::WRITE_SERIAL_FILES) &&
458 : // !this->get_mesh().is_serial())
459 : ;
460 :
461 242 : if (write_parallel_files && write_data)
462 0 : libmesh_assert(local_io);
463 :
464 : {
465 242 : libmesh_assert (io.writing());
466 :
467 484 : LOG_SCOPE("write()", "EquationSystems");
468 :
469 484 : const unsigned int proc_id = this->processor_id();
470 :
471 8344 : unsigned int n_sys = 0;
472 16759 : for (auto & pr : _systems)
473 8415 : if (!pr.second->hide_output())
474 8415 : n_sys++;
475 :
476 : // set the version number in the Xdr object
477 242 : io.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
478 : LIBMESH_MINOR_VERSION,
479 : LIBMESH_MICRO_VERSION));
480 :
481 : // Only write the header information
482 : // if we are processor 0.
483 8344 : if (proc_id == 0)
484 : {
485 242 : std::string comment;
486 :
487 : // 1.)
488 : // Write the version header
489 1440 : std::string version("libMesh-" + libMesh::get_io_compatibility_version());
490 1319 : if (write_parallel_files) version += " parallel";
491 :
492 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
493 121 : version += " with infinite elements";
494 : #endif
495 1319 : io.data (version, "# File Format Identifier");
496 :
497 : // 2.)
498 : // Write the number of equation systems
499 1319 : io.data (n_sys, "# No. of Equation Systems");
500 :
501 2650 : for (auto & [sys_name, sys] : _systems)
502 : {
503 : // Ignore this system if it has been marked as hidden
504 1331 : if (sys->hide_output()) continue;
505 :
506 : // 3.)
507 : // Write the name of the sys_num-th system
508 : {
509 244 : const unsigned int sys_num = sys->number();
510 :
511 122 : comment = "# Name, System No. ";
512 1331 : comment += std::to_string(sys_num);
513 :
514 : // Note: There is no Xdr::data overload taking a "const
515 : // std::string &" so we need to make a copy.
516 1331 : std::string copy = sys_name;
517 1331 : io.data (copy, comment);
518 : }
519 :
520 : // 4.)
521 : // Write the type of system handled
522 : {
523 244 : const unsigned int sys_num = sys->number();
524 1453 : std::string sys_type = sys->system_type();
525 :
526 122 : comment = "# Type, System No. ";
527 2418 : comment += std::to_string(sys_num);
528 :
529 1331 : io.data (sys_type, comment);
530 : }
531 :
532 : // 5.) - 9.)
533 : // Let System::write_header() do the job
534 1453 : sys->write_header (io, version, write_additional_data);
535 : }
536 : }
537 :
538 : // Start from the first system, again,
539 : // to write vectors to disk, if wanted
540 8344 : if (write_data)
541 : {
542 16759 : for (auto & pr : _systems)
543 : {
544 : // Ignore this system if it has been marked as hidden
545 8415 : if (pr.second->hide_output()) continue;
546 :
547 : // 10.) + 11.)
548 8415 : if (write_parallel_files)
549 0 : pr.second->write_parallel_data (*local_io,write_additional_data);
550 : else
551 8415 : pr.second->write_serialized_data (io,write_additional_data);
552 : }
553 :
554 8344 : if (local_io)
555 0 : local_io->close();
556 : }
557 :
558 8344 : io.close();
559 : }
560 :
561 : // the EquationSystems::write() method should look constant,
562 : // but we need to undo the temporary numbering of the nodes
563 : // and elements in the mesh, which requires that we abuse const_cast
564 8344 : if (partition_agnostic)
565 8344 : const_cast<MeshBase &>(_mesh).fix_broken_node_and_element_numbering();
566 8344 : }
567 :
568 :
569 :
570 : // template specialization
571 :
572 : 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);
573 : template LIBMESH_EXPORT void EquationSystems::read<Number> (std::string_view name, const unsigned int read_flags, bool partition_agnostic);
574 : template LIBMESH_EXPORT void EquationSystems::read<Number> (std::string_view name, const XdrMODE mode, const unsigned int read_flags, bool partition_agnostic);
575 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
576 : 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);
577 : template LIBMESH_EXPORT void EquationSystems::read<Real> (std::string_view name, const unsigned int read_flags, bool partition_agnostic);
578 : template LIBMESH_EXPORT void EquationSystems::read<Real> (std::string_view name, const XdrMODE mode, const unsigned int read_flags, bool partition_agnostic);
579 : #endif
580 :
581 : } // namespace libMesh
|
