Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local includes
19 : #include "libmesh/checkpoint_io.h"
20 : #include "libmesh/boundary_info.h"
21 : #include "libmesh/cell_c0polyhedron.h"
22 : #include "libmesh/distributed_mesh.h"
23 : #include "libmesh/elem.h"
24 : #include "libmesh/enum_to_string.h"
25 : #include "libmesh/enum_xdr_mode.h"
26 : #include "libmesh/face_c0polygon.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_communication.h"
30 : #include "libmesh/mesh_tools.h"
31 : #include "libmesh/node.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/partitioner.h"
34 : #include "libmesh/metis_partitioner.h"
35 : #include "libmesh/remote_elem.h"
36 : #include "libmesh/xdr_io.h"
37 : #include "libmesh/xdr_cxx.h"
38 : #include "libmesh/utility.h"
39 : #include "libmesh/int_range.h"
40 :
41 : // C++ includes
42 : #include <algorithm>
43 : #include <iostream>
44 : #include <iomanip>
45 : #include <cstdio>
46 : #include <vector>
47 : #include <string>
48 : #include <cstring>
49 : #include <filesystem>
50 : #include <fstream>
51 : #include <sstream> // for ostringstream
52 : #include <unordered_map>
53 : #include <unordered_set>
54 : #ifdef LIBMESH_HAVE_DIRECT_H
55 : #include <direct.h> // rmdir() on Windows
56 : #endif
57 : #ifdef LIBMESH_HAVE_UNISTD_H
58 : #include <unistd.h> // rmdir() on Unix
59 : #endif
60 :
61 : namespace
62 : {
63 : // chunking computes the number of chunks and first-chunk-offset when splitting a mesh
64 : // into nsplits pieces using size procs for the given MPI rank. The number of chunks and offset
65 : // are stored in nchunks and first_chunk respectively.
66 0 : void chunking(libMesh::processor_id_type size, libMesh::processor_id_type rank, libMesh::processor_id_type nsplits,
67 : libMesh::processor_id_type & nchunks, libMesh::processor_id_type & first_chunk)
68 : {
69 0 : if (nsplits % size == 0) // the chunks divide evenly over the processors
70 : {
71 0 : nchunks = nsplits / size;
72 0 : first_chunk = libMesh::cast_int<libMesh::processor_id_type>(nchunks * rank);
73 0 : return;
74 : }
75 :
76 0 : libMesh::processor_id_type nextra = nsplits % size;
77 0 : if (rank < nextra) // leftover chunks cause an extra chunk to be added to this processor
78 : {
79 0 : nchunks = libMesh::cast_int<libMesh::processor_id_type>(nsplits / size + 1);
80 0 : first_chunk = libMesh::cast_int<libMesh::processor_id_type>(nchunks * rank);
81 : }
82 : else // no extra chunks, but first chunk is offset by extras on earlier ranks
83 : {
84 0 : nchunks = nsplits / size;
85 : // account for the case where nchunks is zero where we want max int
86 0 : first_chunk = libMesh::cast_int<libMesh::processor_id_type>
87 0 : (std::max((int)((nchunks + 1) * (nsplits % size) + nchunks * (rank - nsplits % size)),
88 0 : (1 - (int)nchunks) * std::numeric_limits<int>::max()));
89 : }
90 : }
91 :
92 2604 : std::string_view extension(std::string_view s)
93 : {
94 2482 : auto pos = s.rfind(".");
95 2604 : if (pos == std::string::npos)
96 0 : return "";
97 2604 : return s.substr(pos, s.size() - pos);
98 : }
99 :
100 3434 : std::string split_dir(const std::string & input_name, libMesh::processor_id_type n_procs)
101 : {
102 6716 : return input_name + "/" + std::to_string(n_procs);
103 : }
104 :
105 :
106 1780 : std::string header_file(const std::string & input_name, libMesh::processor_id_type n_procs)
107 : {
108 6840 : return (split_dir(input_name, n_procs) + "/header").append(extension(input_name));
109 : }
110 :
111 : std::string
112 824 : split_file(const std::string & input_name,
113 : libMesh::processor_id_type n_procs,
114 : libMesh::processor_id_type proc_id)
115 : {
116 3088 : return (split_dir(input_name, n_procs) + "/split-" + std::to_string(n_procs) + "-" +
117 3192 : std::to_string(proc_id)).append(extension(input_name));
118 : }
119 :
120 : // Report, as a comma-separated list, the processor counts for which a split
121 : // already exists as a numbered subdirectory of input_name, to help a user
122 : // pick out why the split they asked for could not be found.
123 12 : std::string existing_splits(const std::string & input_name)
124 : {
125 1 : std::error_code ec;
126 23 : if (!std::filesystem::is_directory(input_name, ec))
127 0 : return "";
128 :
129 3 : std::vector<std::string> splits;
130 59 : for (const auto & entry : std::filesystem::directory_iterator(input_name, ec))
131 : {
132 24 : const auto name = entry.path().filename().string();
133 24 : if (!name.empty() && name.find_first_not_of("0123456789") == std::string::npos &&
134 1 : entry.is_directory(ec))
135 12 : splits.push_back(name);
136 : }
137 12 : std::sort(splits.begin(), splits.end());
138 :
139 2 : std::string list;
140 24 : for (const auto & split : splits)
141 23 : list += (list.empty() ? "" : ", ") + split;
142 11 : return list;
143 10 : }
144 :
145 770 : void make_dir(const std::string & input_name, libMesh::processor_id_type n_procs)
146 : {
147 770 : auto ret = libMesh::Utility::mkdir(input_name.c_str());
148 : // error only if we failed to create dir - don't care if it was already there
149 770 : libmesh_error_msg_if
150 : (ret != 0 && ret != -1,
151 : "Failed to create mesh split directory '" << input_name << "': " << std::strerror(ret));
152 :
153 795 : auto dir_name = split_dir(input_name, n_procs);
154 770 : ret = libMesh::Utility::mkdir(dir_name.c_str());
155 770 : if (ret == -1)
156 : libmesh_warning("In CheckpointIO::write, directory '"
157 : << dir_name << "' already exists, overwriting contents.");
158 : else
159 94 : libmesh_error_msg_if
160 : (ret != 0, "Failed to create mesh split directory '" << dir_name << "': " << std::strerror(ret));
161 770 : }
162 :
163 : } // namespace
164 :
165 : namespace libMesh
166 : {
167 :
168 0 : std::unique_ptr<CheckpointIO> split_mesh(MeshBase & mesh, processor_id_type nsplits)
169 : {
170 : // There is currently an issue with DofObjects not being properly
171 : // reset if the mesh is not first repartitioned onto 1 processor
172 : // *before* being repartitioned onto the desired number of
173 : // processors. So, this is a workaround, but not a particularly
174 : // onerous one.
175 0 : mesh.partition(1);
176 0 : mesh.partition(nsplits);
177 :
178 0 : processor_id_type my_num_chunks = 0;
179 0 : processor_id_type my_first_chunk = 0;
180 0 : chunking(mesh.comm().size(), mesh.comm().rank(), nsplits, my_num_chunks, my_first_chunk);
181 :
182 0 : auto cpr = std::make_unique<CheckpointIO>(mesh);
183 0 : cpr->current_processor_ids().clear();
184 0 : for (processor_id_type i = my_first_chunk; i < my_first_chunk + my_num_chunks; i++)
185 0 : cpr->current_processor_ids().push_back(i);
186 0 : cpr->current_n_processors() = nsplits;
187 0 : cpr->parallel() = true;
188 0 : return cpr;
189 0 : }
190 :
191 :
192 : // ------------------------------------------------------------
193 : // CheckpointIO members
194 1398 : CheckpointIO::CheckpointIO (MeshBase & mesh, const bool binary_in) :
195 : MeshInput<MeshBase> (mesh,/* is_parallel_format = */ true),
196 : MeshOutput<MeshBase>(mesh,/* is_parallel_format = */ true),
197 : ParallelObject (mesh),
198 1306 : _binary (binary_in),
199 1306 : _parallel (false),
200 0 : _version ("checkpoint-1.6"),
201 1490 : _my_processor_ids (1, processor_id()),
202 2750 : _my_n_processors (mesh.is_replicated() ? 1 : n_processors())
203 : {
204 1398 : }
205 :
206 142 : CheckpointIO::CheckpointIO (const MeshBase & mesh, const bool binary_in) :
207 : MeshInput<MeshBase> (), // write-only
208 : MeshOutput<MeshBase>(mesh,/* is_parallel_format = */ true),
209 : ParallelObject (mesh),
210 134 : _binary (binary_in),
211 134 : _parallel (false),
212 0 : _version ("checkpoint-1.6"),
213 150 : _my_processor_ids (1, processor_id()),
214 280 : _my_n_processors (mesh.is_replicated() ? 1 : n_processors())
215 : {
216 142 : }
217 :
218 1440 : CheckpointIO::~CheckpointIO () = default;
219 :
220 770 : processor_id_type CheckpointIO::select_split_config(const std::string & input_name, header_id_type & data_size)
221 : {
222 50 : std::string header_name;
223 :
224 : // We'll read a header file from processor 0 and broadcast.
225 795 : if (this->processor_id() == 0)
226 : {
227 345 : header_name = header_file(input_name, _my_n_processors);
228 :
229 : {
230 : // look for header+splits with nprocs equal to _my_n_processors
231 210 : std::ifstream in (header_name.c_str());
232 180 : if (!in.good())
233 : {
234 : // otherwise fall back to a serial/single-split mesh
235 12 : auto orig_header_name = header_name;
236 33 : header_name = header_file(input_name, 1);
237 14 : std::ifstream in2 (header_name.c_str());
238 13 : const auto splits = existing_splits(input_name);
239 110 : libmesh_error_msg_if(!in2.good(),
240 : "ERROR: Neither one of the following files can be located:\n\t'"
241 : << orig_header_name << "' nor\n\t'" << input_name << "'\n"
242 : << "If you are running a parallel job, double check that you've "
243 : << "created a split for " << _my_n_processors << " ranks.\n"
244 : << "Note: One of paths above may refer to a valid directory on your "
245 : << "system, however we are attempting to read a valid header file.\n"
246 : << (splits.empty()
247 : ? "No processor-count splits were found in '" + input_name + "'."
248 : : "Splits that do exist in '" + input_name +
249 : "': " + splits + "."));
250 10 : }
251 150 : }
252 :
253 445 : Xdr io (header_name, this->binary() ? DECODE : READ);
254 :
255 : // read the version, but don't care about it
256 14 : std::string input_version;
257 168 : io.data(input_version);
258 :
259 : // read the data type
260 168 : io.data (data_size);
261 140 : }
262 :
263 758 : this->comm().broadcast(data_size);
264 758 : this->comm().broadcast(header_name);
265 :
266 : // How many per-processor files are here?
267 : largest_id_type input_n_procs;
268 :
269 758 : switch (data_size) {
270 0 : case 2:
271 0 : input_n_procs = this->read_header<uint16_t>(header_name);
272 0 : break;
273 0 : case 4:
274 0 : input_n_procs = this->read_header<uint32_t>(header_name);
275 0 : break;
276 758 : case 8:
277 758 : input_n_procs = this->read_header<uint64_t>(header_name);
278 24 : break;
279 0 : default:
280 0 : libmesh_error();
281 : }
282 :
283 48 : if (!input_n_procs)
284 8 : input_n_procs = 1;
285 782 : return cast_int<processor_id_type>(input_n_procs);
286 : }
287 :
288 60 : void CheckpointIO::cleanup(const std::string & input_name, processor_id_type n_procs)
289 : {
290 65 : auto header = header_file(input_name, n_procs);
291 60 : auto ret = std::remove(header.c_str());
292 : if (ret != 0)
293 : libmesh_warning("Failed to clean up checkpoint header '" << header << "': " << std::strerror(ret));
294 :
295 132 : for (processor_id_type i = 0; i < n_procs; i++)
296 : {
297 78 : auto split = split_file(input_name, n_procs, i);
298 72 : ret = std::remove(split.c_str());
299 : if (ret != 0)
300 : libmesh_warning("Failed to clean up checkpoint split file '" << split << "': " << std::strerror(ret));
301 : }
302 :
303 65 : auto dir = split_dir(input_name, n_procs);
304 60 : ret = rmdir(dir.c_str());
305 : if (ret != 0)
306 : libmesh_warning("Failed to clean up checkpoint split dir '" << dir << "': " << std::strerror(ret));
307 :
308 : // We expect that this may fail if there are other split configurations still present in this
309 : // directory - so don't bother to check/warn for failure.
310 60 : rmdir(input_name.c_str());
311 60 : }
312 :
313 :
314 1852 : bool CheckpointIO::version_at_least_1_5() const
315 : {
316 : return
317 3704 : (this->version().find("1.5") != std::string::npos) ||
318 1973 : (this->version().find("1.6") != std::string::npos);
319 : }
320 :
321 :
322 752 : bool CheckpointIO::version_at_least_1_6() const
323 : {
324 752 : return (this->version().find("1.6") != std::string::npos);
325 : }
326 :
327 :
328 770 : void CheckpointIO::write (const std::string & name)
329 : {
330 75 : LOG_SCOPE("write()", "CheckpointIO");
331 :
332 : // convenient reference to our mesh
333 50 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
334 :
335 : // FIXME: For backwards compatibility, we'll assume for now that we
336 : // only want to write distributed meshes in parallel. Later we can
337 : // do a gather_to_zero() and support that case too.
338 770 : _parallel = _parallel || !mesh.is_serial();
339 :
340 25 : processor_id_type use_n_procs = 1;
341 770 : if (_parallel)
342 740 : use_n_procs = _my_n_processors;
343 :
344 795 : std::string header_file_name = header_file(name, use_n_procs);
345 770 : make_dir(name, use_n_procs);
346 :
347 : // We'll write a header file from processor 0 to make it easier to do unambiguous
348 : // restarts later:
349 795 : if (this->processor_id() == 0)
350 : {
351 470 : Xdr io (header_file_name, this->binary() ? ENCODE : WRITE);
352 :
353 : // write the version
354 180 : io.data(_version, "# version");
355 :
356 : // write what kind of data type we're using
357 180 : header_id_type data_size = sizeof(largest_id_type);
358 180 : io.data(data_size, "# integer size");
359 :
360 : // Write out the max mesh dimension for backwards compatibility
361 : // with code that sets it independently of element dimensions
362 : {
363 180 : uint16_t mesh_dimension = cast_int<uint16_t>(mesh.mesh_dimension());
364 180 : io.data(mesh_dimension, "# dimensions");
365 : }
366 :
367 : // Write out whether or not this is serial output
368 : {
369 180 : uint16_t parallel = _parallel;
370 180 : io.data(parallel, "# parallel");
371 : }
372 :
373 : // If we're writing out a parallel mesh then we need to write the number of processors
374 : // so we can check it upon reading the file
375 180 : if (_parallel)
376 : {
377 156 : largest_id_type n_procs = _my_n_processors;
378 156 : io.data(n_procs, "# n_procs");
379 : }
380 :
381 : // write subdomain names
382 180 : this->write_subdomain_names(io);
383 :
384 : // write boundary id names
385 15 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
386 180 : write_bc_names(io, boundary_info, true); // sideset names
387 180 : write_bc_names(io, boundary_info, false); // nodeset names
388 :
389 : // write extra integer names
390 180 : const bool write_extra_integers = this->version_at_least_1_5();
391 :
392 180 : if (write_extra_integers)
393 : {
394 180 : largest_id_type n_node_integers = mesh.n_node_integers();
395 180 : io.data(n_node_integers, "# n_extra_integers per node");
396 :
397 45 : std::vector<std::string> node_integer_names;
398 288 : for (unsigned int i=0; i != n_node_integers; ++i)
399 108 : node_integer_names.push_back(mesh.get_node_integer_name(i));
400 :
401 180 : io.data(node_integer_names);
402 :
403 180 : largest_id_type n_elem_integers = mesh.n_elem_integers();
404 180 : io.data(n_elem_integers, "# n_extra_integers per elem");
405 :
406 30 : std::vector<std::string> elem_integer_names;
407 234 : for (unsigned int i=0; i != n_elem_integers; ++i)
408 54 : elem_integer_names.push_back(mesh.get_elem_integer_name(i));
409 :
410 180 : io.data(elem_integer_names);
411 150 : }
412 :
413 :
414 150 : }
415 :
416 : // If this is a serial mesh written to a serial file then we're only
417 : // going to write local data from processor 0. If this is a mesh being
418 : // written in parallel then we're going to write from every
419 : // processor.
420 50 : std::vector<processor_id_type> ids_to_write;
421 :
422 : // We're going to sort elements by pid in one pass, to avoid sending
423 : // predicated iterators through the whole mesh N_p times.
424 : //
425 : // The data type here needs to be a non-const-pointer to whatever
426 : // our element_iterator is a const-pointer to, for compatibility
427 : // later.
428 : typedef std::remove_const<MeshBase::const_element_iterator::value_type>::type nc_v_t;
429 50 : std::unordered_map<processor_id_type, std::vector<nc_v_t>> elements_on_pid;
430 :
431 770 : if (_parallel)
432 : {
433 740 : ids_to_write = _my_processor_ids;
434 1104 : for (processor_id_type p : ids_to_write)
435 18 : elements_on_pid[p].clear();
436 17 : auto eop_end = elements_on_pid.end();
437 13487 : for (auto & elem : mesh.element_ptr_range())
438 : {
439 6240 : const processor_id_type p = elem->processor_id();
440 6468 : if (auto eop_it = elements_on_pid.find(p);
441 228 : eop_it != eop_end)
442 1720 : eop_it->second.push_back(elem);
443 706 : }
444 : }
445 30 : else if (mesh.is_serial())
446 : {
447 38 : if (mesh.processor_id() == 0)
448 : {
449 : // placeholder
450 24 : ids_to_write.push_back(0);
451 : }
452 : }
453 : else
454 : {
455 0 : libmesh_error_msg("Cannot write serial checkpoint from distributed mesh");
456 : }
457 :
458 : // Call build_side_list() and build_node_list() just *once* to avoid
459 : // redundant expensive sorts during mesh splitting.
460 25 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
461 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
462 795 : bc_triples = boundary_info.build_side_list();
463 : std::vector<std::tuple<dof_id_type, boundary_id_type>>
464 795 : bc_tuples = boundary_info.build_node_list();
465 :
466 1158 : for (const auto & my_pid : ids_to_write)
467 : {
468 412 : auto file_name = split_file(name, use_n_procs, my_pid);
469 1004 : Xdr io (file_name, this->binary() ? ENCODE : WRITE);
470 :
471 48 : std::set<const Elem *, CompareElemIdsByLevel> elements;
472 :
473 : // For serial files or for already-distributed meshs, we write
474 : // everything we can see.
475 388 : if (!_parallel || !mesh.is_serial())
476 530 : elements.insert(mesh.elements_begin(), mesh.elements_end());
477 : // For parallel files written from serial meshes we write what
478 : // we'd be required to keep if we were to be deleting remote
479 : // elements. This allows us to write proper parallel files even
480 : // from a ReplicateMesh.
481 : //
482 : // WARNING: If we have a DistributedMesh which used
483 : // "add_extra_ghost_elem" rather than ghosting functors to
484 : // preserve elements and which is *also* currently serialized
485 : // then we're not preserving those elements here. As a quick
486 : // workaround user code should delete_remote_elements() before
487 : // writing the checkpoint; as a long term workaround user code
488 : // should use ghosting functors instead of extra_ghost_elem
489 : // lists.
490 : else
491 : {
492 348 : for (processor_id_type p : {my_pid, DofObject::invalid_processor_id})
493 : {
494 232 : if (const auto elements_vec_it = elements_on_pid.find(p);
495 20 : elements_vec_it != elements_on_pid.end())
496 : {
497 10 : auto & p_elements = elements_vec_it->second;
498 :
499 : // Be compatible with both deprecated and
500 : // corrected MeshBase iterator types
501 : typedef MeshBase::const_element_iterator::value_type v_t;
502 :
503 116 : v_t * elempp = p_elements.data();
504 116 : v_t * elemend = elempp + p_elements.size();
505 :
506 : const MeshBase::const_element_iterator
507 : pid_elements_begin = MeshBase::const_element_iterator
508 126 : (elempp, elemend, Predicates::NotNull<v_t *>()),
509 : pid_elements_end = MeshBase::const_element_iterator
510 126 : (elemend, elemend, Predicates::NotNull<v_t *>()),
511 : active_pid_elements_begin = MeshBase::const_element_iterator
512 126 : (elempp, elemend, Predicates::Active<v_t *>()),
513 : active_pid_elements_end = MeshBase::const_element_iterator
514 222 : (elemend, elemend, Predicates::Active<v_t *>());
515 :
516 : query_ghosting_functors
517 222 : (mesh, p, active_pid_elements_begin,
518 : active_pid_elements_end, elements);
519 222 : connect_children(mesh, pid_elements_begin,
520 : pid_elements_end, elements);
521 : }
522 : }
523 : }
524 :
525 48 : connected_node_set_type connected_nodes;
526 388 : connect_element_dependencies(mesh, elements, connected_nodes);
527 :
528 : // write the nodal locations
529 388 : this->write_nodes (io, connected_nodes);
530 :
531 : // write connectivity
532 388 : this->write_connectivity (io, elements);
533 :
534 : // write remote_elem connectivity
535 388 : this->write_remote_elem (io, elements);
536 :
537 : // write the boundary condition information
538 388 : this->write_bcs (io, elements, bc_triples);
539 :
540 : // write the nodeset information
541 388 : this->write_nodesets (io, connected_nodes, bc_tuples);
542 :
543 : // close it up
544 388 : io.close();
545 340 : }
546 :
547 : // this->comm().barrier();
548 770 : }
549 :
550 180 : void CheckpointIO::write_subdomain_names(Xdr & io) const
551 : {
552 : {
553 30 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
554 :
555 15 : const std::map<subdomain_id_type, std::string> & subdomain_map = mesh.get_subdomain_name_map();
556 :
557 195 : std::vector<largest_id_type> subdomain_ids; subdomain_ids.reserve(subdomain_map.size());
558 210 : std::vector<std::string> subdomain_names; subdomain_names.reserve(subdomain_map.size());
559 :
560 : // We need to loop over the map and make sure that there aren't any invalid entries. Since we
561 : // return writable references in mesh_base, it's possible for the user to leave some entity names
562 : // blank. We can't write those to the XDA file.
563 180 : largest_id_type n_subdomain_names = 0;
564 180 : for (const auto & [id, name] : subdomain_map)
565 0 : if (!name.empty())
566 : {
567 0 : n_subdomain_names++;
568 0 : subdomain_ids.push_back(id);
569 0 : subdomain_names.push_back(name);
570 : }
571 :
572 180 : io.data(n_subdomain_names, "# subdomain id to name map");
573 : // Write out the ids and names in two vectors
574 180 : if (n_subdomain_names)
575 : {
576 0 : io.data(subdomain_ids);
577 0 : io.data(subdomain_names);
578 : }
579 150 : }
580 180 : }
581 :
582 :
583 :
584 388 : void CheckpointIO::write_nodes (Xdr & io,
585 : const connected_node_set_type & nodeset) const
586 : {
587 388 : largest_id_type n_nodes_here = nodeset.size();
588 :
589 388 : io.data(n_nodes_here, "# n_nodes on proc");
590 :
591 388 : const bool write_extra_integers = this->version_at_least_1_5();
592 : const unsigned int n_extra_integers =
593 388 : write_extra_integers ? MeshOutput<MeshBase>::mesh().n_node_integers() : 0;
594 :
595 : // Will hold the node id and pid and extra integers
596 412 : std::vector<largest_id_type> id_pid(2 + n_extra_integers);
597 :
598 : // For the coordinates
599 412 : std::vector<Real> coords(LIBMESH_DIM);
600 :
601 6374 : for (const auto & node : nodeset)
602 : {
603 5986 : id_pid[0] = node->id();
604 5986 : id_pid[1] = node->processor_id();
605 :
606 412 : libmesh_assert_equal_to(n_extra_integers, node->n_extra_integers());
607 12878 : for (unsigned int i=0; i != n_extra_integers; ++i)
608 7348 : id_pid[2+i] = node->get_extra_integer(i);
609 :
610 5986 : io.data_stream(id_pid.data(), 2 + n_extra_integers, 2 + n_extra_integers);
611 :
612 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
613 5986 : largest_id_type unique_id = node->unique_id();
614 :
615 5986 : io.data(unique_id, "# unique id");
616 : #endif
617 :
618 5986 : coords[0] = (*node)(0);
619 :
620 : #if LIBMESH_DIM > 1
621 5986 : coords[1] = (*node)(1);
622 : #endif
623 :
624 : #if LIBMESH_DIM > 2
625 5986 : coords[2] = (*node)(2);
626 : #endif
627 :
628 5986 : io.data_stream(coords.data(), LIBMESH_DIM, 3);
629 : }
630 388 : }
631 :
632 :
633 :
634 388 : void CheckpointIO::write_connectivity (Xdr & io,
635 : const std::set<const Elem *, CompareElemIdsByLevel> & elements) const
636 : {
637 24 : libmesh_assert (io.writing());
638 :
639 388 : const bool write_extra_integers = this->version_at_least_1_5();
640 388 : const bool write_runtime_topology = this->version_at_least_1_6();
641 : const unsigned int n_extra_integers =
642 388 : write_extra_integers ? MeshOutput<MeshBase>::mesh().n_elem_integers() : 0;
643 :
644 : // Put these out here to reduce memory churn
645 : // id type pid subdomain_id parent_id extra_integer_0 ...
646 412 : std::vector<largest_id_type> elem_data(6 + n_extra_integers);
647 48 : std::vector<largest_id_type> conn_data;
648 48 : std::vector<largest_id_type> runtime_topology;
649 :
650 388 : largest_id_type n_elems_here = elements.size();
651 :
652 388 : io.data(n_elems_here, "# number of elements");
653 :
654 3276 : for (const auto & elem : elements)
655 : {
656 2888 : unsigned int n_nodes = elem->n_nodes();
657 :
658 2888 : elem_data[0] = elem->id();
659 2888 : elem_data[1] = elem->type();
660 2888 : elem_data[2] = elem->processor_id();
661 2888 : elem_data[3] = elem->subdomain_id();
662 :
663 : #ifdef LIBMESH_ENABLE_AMR
664 3105 : if (elem->parent() != nullptr)
665 : {
666 0 : elem_data[4] = elem->parent()->id();
667 0 : elem_data[5] = elem->parent()->which_child_am_i(elem);
668 : }
669 : else
670 : #endif
671 : {
672 2888 : elem_data[4] = static_cast<largest_id_type>(-1);
673 2888 : elem_data[5] = static_cast<largest_id_type>(-1);
674 : }
675 :
676 3983 : for (unsigned int i=0; i != n_extra_integers; ++i)
677 1173 : elem_data[6+i] = elem->get_extra_integer(i);
678 :
679 2888 : conn_data.resize(n_nodes);
680 :
681 16232 : for (unsigned int i=0; i<n_nodes; i++)
682 15174 : conn_data[i] = elem->node_id(i);
683 :
684 3105 : io.data_stream(elem_data.data(),
685 : cast_int<unsigned int>(elem_data.size()),
686 : cast_int<unsigned int>(elem_data.size()));
687 :
688 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
689 2888 : largest_id_type unique_id = elem->unique_id();
690 :
691 2888 : io.data(unique_id, "# unique id");
692 : #endif
693 :
694 : #ifdef LIBMESH_ENABLE_AMR
695 3105 : uint16_t p_level = cast_int<uint16_t>(elem->p_level());
696 2888 : io.data(p_level, "# p_level");
697 :
698 2888 : uint16_t rflag = elem->refinement_flag();
699 2888 : io.data(rflag, "# rflag");
700 :
701 2888 : uint16_t pflag = elem->p_refinement_flag();
702 2888 : io.data(pflag, "# pflag");
703 : #endif
704 :
705 2888 : if (elem->runtime_topology())
706 : {
707 48 : libmesh_error_msg_if
708 : (!write_runtime_topology,
709 : "Checkpoint format 1.6 or newer is required to write " <<
710 : Utility::enum_to_string(elem->type()) << " elements.");
711 :
712 4 : runtime_topology.clear();
713 48 : runtime_topology.push_back(n_nodes);
714 48 : runtime_topology.push_back(elem->n_sides());
715 364 : for (auto s : elem->side_index_range())
716 : {
717 338 : const auto side_nodes = elem->nodes_on_side(s);
718 338 : runtime_topology.push_back(side_nodes.size());
719 1416 : for (const auto n : side_nodes)
720 1104 : runtime_topology.push_back(n);
721 : }
722 :
723 48 : io.data(runtime_topology, "# runtime topology");
724 : }
725 :
726 3105 : io.data_stream(conn_data.data(),
727 : cast_int<unsigned int>(conn_data.size()),
728 : cast_int<unsigned int>(conn_data.size()));
729 : }
730 388 : }
731 :
732 :
733 388 : void CheckpointIO::write_remote_elem (Xdr & io,
734 : const std::set<const Elem *, CompareElemIdsByLevel> & elements) const
735 : {
736 24 : libmesh_assert (io.writing());
737 :
738 : // Find the remote_elem neighbor and child links
739 48 : std::vector<largest_id_type> elem_ids, parent_ids;
740 48 : std::vector<uint16_t> elem_sides, child_numbers;
741 :
742 3276 : for (const auto & elem : elements)
743 : {
744 15385 : for (auto n : elem->side_index_range())
745 : {
746 12280 : const Elem * neigh = elem->neighbor_ptr(n);
747 13132 : if (neigh == remote_elem ||
748 852 : (neigh && !elements.count(neigh)))
749 : {
750 708 : elem_ids.push_back(elem->id());
751 708 : elem_sides.push_back(n);
752 : }
753 : }
754 :
755 : #ifdef LIBMESH_ENABLE_AMR
756 2888 : if (elem->has_children())
757 : {
758 0 : for (unsigned short c = 0,
759 0 : nc = cast_int<unsigned short>(elem->n_children());
760 0 : c != nc; ++c)
761 : {
762 0 : const Elem * child = elem->child_ptr(c);
763 0 : if (child == remote_elem ||
764 0 : (child && !elements.count(child)))
765 : {
766 0 : parent_ids.push_back(elem->id());
767 0 : child_numbers.push_back(c);
768 : }
769 : }
770 : }
771 : #endif
772 : }
773 :
774 388 : io.data(elem_ids, "# remote neighbor elem_ids");
775 388 : io.data(elem_sides, "# remote neighbor elem_sides");
776 388 : io.data(parent_ids, "# remote child parent_ids");
777 388 : io.data(child_numbers, "# remote child_numbers");
778 388 : }
779 :
780 :
781 :
782 388 : void CheckpointIO::write_bcs (Xdr & io,
783 : const std::set<const Elem *, CompareElemIdsByLevel> & elements,
784 : const std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> & bc_triples) const
785 : {
786 24 : libmesh_assert (io.writing());
787 :
788 : // Build a list of (elem, side, bc) tuples.
789 48 : std::size_t bc_size = bc_triples.size();
790 :
791 48 : std::vector<largest_id_type> element_id_list;
792 48 : std::vector<uint16_t> side_list;
793 48 : std::vector<largest_id_type> bc_id_list;
794 :
795 388 : element_id_list.reserve(bc_size);
796 388 : side_list.reserve(bc_size);
797 388 : bc_id_list.reserve(bc_size);
798 :
799 24 : std::unordered_set<dof_id_type> elems;
800 3276 : for (auto & e : elements)
801 2888 : elems.insert(e->id());
802 :
803 4276 : for (const auto & t : bc_triples)
804 500 : if (elems.count(std::get<0>(t)))
805 : {
806 3360 : element_id_list.push_back(std::get<0>(t));
807 3360 : side_list.push_back(std::get<1>(t));
808 3360 : bc_id_list.push_back(std::get<2>(t));
809 : }
810 :
811 :
812 388 : io.data(element_id_list, "# element ids for bcs");
813 388 : io.data(side_list, "# sides of elements for bcs");
814 388 : io.data(bc_id_list, "# bc ids");
815 388 : }
816 :
817 :
818 :
819 388 : void CheckpointIO::write_nodesets (Xdr & io,
820 : const connected_node_set_type & nodeset,
821 : const std::vector<std::tuple<dof_id_type, boundary_id_type>> & bc_tuples) const
822 : {
823 24 : libmesh_assert (io.writing());
824 :
825 : // convenient reference to our mesh
826 48 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
827 :
828 : // Build a list of (node, bc) tuples
829 48 : std::size_t nodeset_size = bc_tuples.size();
830 :
831 48 : std::vector<largest_id_type> node_id_list;
832 24 : std::vector<largest_id_type> bc_id_list;
833 :
834 388 : node_id_list.reserve(nodeset_size);
835 388 : bc_id_list.reserve(nodeset_size);
836 :
837 6576 : for (const auto & t : bc_tuples)
838 6244 : if (nodeset.count(mesh.node_ptr(std::get<0>(t))))
839 : {
840 5572 : node_id_list.push_back(std::get<0>(t));
841 5572 : bc_id_list.push_back(std::get<1>(t));
842 : }
843 :
844 388 : io.data(node_id_list, "# node id list");
845 388 : io.data(bc_id_list, "# nodeset bc id list");
846 388 : }
847 :
848 :
849 :
850 360 : void CheckpointIO::write_bc_names (Xdr & io, const BoundaryInfo & info, bool is_sideset) const
851 : {
852 360 : const std::map<boundary_id_type, std::string> & boundary_map = is_sideset ?
853 30 : info.get_sideset_name_map() : info.get_nodeset_name_map();
854 :
855 390 : std::vector<largest_id_type> boundary_ids; boundary_ids.reserve(boundary_map.size());
856 420 : std::vector<std::string> boundary_names; boundary_names.reserve(boundary_map.size());
857 :
858 : // We need to loop over the map and make sure that there aren't any invalid entries. Since we
859 : // return writable references in boundary_info, it's possible for the user to leave some entity names
860 : // blank. We can't write those to the XDA file.
861 360 : largest_id_type n_boundary_names = 0;
862 1416 : for (const auto & [id, name] : boundary_map)
863 1056 : if (!name.empty())
864 : {
865 1056 : n_boundary_names++;
866 1056 : boundary_ids.push_back(id);
867 1056 : boundary_names.push_back(name);
868 : }
869 :
870 360 : if (is_sideset)
871 180 : io.data(n_boundary_names, "# sideset id to name map");
872 : else
873 180 : io.data(n_boundary_names, "# nodeset id to name map");
874 : // Write out the ids and names in two vectors
875 360 : if (n_boundary_names)
876 : {
877 264 : io.data(boundary_ids);
878 264 : io.data(boundary_names);
879 : }
880 660 : }
881 :
882 770 : void CheckpointIO::read (const std::string & input_name)
883 : {
884 75 : LOG_SCOPE("read()","CheckpointIO");
885 :
886 770 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
887 :
888 25 : libmesh_assert(!mesh.n_elem());
889 :
890 : header_id_type data_size;
891 770 : processor_id_type input_n_procs = select_split_config(input_name, data_size);
892 782 : auto header_name = header_file(input_name, input_n_procs);
893 24 : bool input_parallel = input_n_procs > 0;
894 :
895 : // If this is a serial read then we're going to only read the mesh
896 : // on processor 0, then broadcast it
897 758 : if ((input_parallel && !mesh.is_replicated()) || mesh.processor_id() == 0)
898 : {
899 : // If we're trying to read a parallel checkpoint file on a
900 : // replicated mesh, we'll read every file on processor 0 so we
901 : // can broadcast it later. If we're on a distributed mesh then
902 : // we'll read every id to it's own processor and we'll "wrap
903 : // around" with any ids that exceed our processor count.
904 : const processor_id_type begin_proc_id =
905 516 : (input_parallel && !mesh.is_replicated()) ?
906 16 : mesh.processor_id() : 0;
907 : const processor_id_type stride =
908 516 : (input_parallel && !mesh.is_replicated()) ?
909 16 : mesh.n_processors() : 1;
910 :
911 880 : for (processor_id_type proc_id = begin_proc_id; proc_id < input_n_procs;
912 364 : proc_id = cast_int<processor_id_type>(proc_id + stride))
913 : {
914 386 : auto file_name = split_file(input_name, input_n_procs, proc_id);
915 :
916 : {
917 408 : std::ifstream in (file_name.c_str());
918 :
919 364 : libmesh_error_msg_if(!in.good(), "ERROR: cannot locate specified file:\n\t" << file_name);
920 320 : }
921 :
922 : // Do we expect all our files' remote_elem entries to really
923 : // be remote? Only if we're not reading multiple input
924 : // files on the same processor.
925 : const bool expect_all_remote =
926 742 : (input_n_procs <= mesh.n_processors() &&
927 356 : !mesh.is_replicated());
928 :
929 590 : Xdr io (file_name, this->binary() ? DECODE : READ);
930 :
931 364 : switch (data_size) {
932 0 : case 2:
933 0 : this->read_subfile<uint16_t>(io, expect_all_remote);
934 0 : break;
935 0 : case 4:
936 0 : this->read_subfile<uint32_t>(io, expect_all_remote);
937 0 : break;
938 364 : case 8:
939 364 : this->read_subfile<uint64_t>(io, expect_all_remote);
940 22 : break;
941 0 : default:
942 0 : libmesh_error();
943 : }
944 :
945 364 : io.close();
946 320 : }
947 : }
948 :
949 : // If the mesh was only read on processor 0 then we need to broadcast it
950 758 : if (mesh.is_replicated())
951 314 : MeshCommunication().broadcast(mesh);
952 : // If the mesh is really distributed then we need to make sure it
953 : // knows that
954 452 : else if (mesh.n_processors() > 1)
955 402 : mesh.set_distributed();
956 :
957 : // If the mesh isn't getting even critical partitioning then we
958 : // should update cached data from the partitioning we just read in
959 758 : if (mesh.skip_partitioning())
960 : {
961 0 : mesh.recalculate_n_partitions();
962 0 : mesh.update_post_partitioning();
963 : }
964 758 : }
965 :
966 :
967 :
968 : template <typename file_id_type>
969 758 : file_id_type CheckpointIO::read_header (const std::string & name)
970 : {
971 758 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
972 :
973 : // Hack for codes which don't look at all elem dimensions
974 : uint16_t mesh_dimension;
975 :
976 : // Will this be a parallel input file? With how many processors? Stay tuned!
977 : uint16_t input_parallel;
978 : file_id_type input_n_procs;
979 :
980 48 : std::string input_version;
981 72 : std::vector<std::string> node_integer_names, elem_integer_names;
982 :
983 : // We'll write a header file from processor 0 and broadcast.
984 782 : if (this->processor_id() == 0)
985 : {
986 434 : Xdr io (name, this->binary() ? DECODE : READ);
987 :
988 : // read the version
989 168 : io.data(input_version);
990 :
991 : // read the data type, don't care about it this time
992 : header_id_type data_size;
993 168 : io.data (data_size);
994 :
995 : // read the dimension
996 168 : io.data (mesh_dimension);
997 :
998 : // Read whether or not this is a parallel file
999 168 : io.data(input_parallel);
1000 :
1001 : // With how many processors?
1002 168 : if (input_parallel)
1003 144 : io.data(input_n_procs);
1004 :
1005 : // read subdomain names
1006 168 : this->read_subdomain_names<file_id_type>(io);
1007 :
1008 : // read boundary names
1009 14 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1010 :
1011 168 : this->read_bc_names<file_id_type>(io, boundary_info, true); // sideset names
1012 168 : this->read_bc_names<file_id_type>(io, boundary_info, false); // nodeset names
1013 :
1014 : // read extra integer names?
1015 14 : std::swap(input_version, this->version());
1016 168 : const bool read_extra_integers = this->version_at_least_1_5();
1017 14 : std::swap(input_version, this->version());
1018 :
1019 168 : if (read_extra_integers)
1020 28 : this->read_integers_names<file_id_type>
1021 140 : (io, node_integer_names, elem_integer_names);
1022 140 : }
1023 :
1024 : // broadcast data from processor 0, set values everywhere
1025 758 : this->comm().broadcast(input_version);
1026 24 : this->version() = input_version;
1027 :
1028 758 : this->comm().broadcast(mesh_dimension);
1029 1492 : mesh.set_mesh_dimension(cast_int<unsigned char>(mesh_dimension));
1030 :
1031 758 : this->comm().broadcast(input_parallel);
1032 :
1033 758 : if (input_parallel)
1034 728 : this->comm().broadcast(input_n_procs);
1035 : else
1036 30 : input_n_procs = 1;
1037 :
1038 : std::map<subdomain_id_type, std::string> & subdomain_map =
1039 24 : mesh.set_subdomain_name_map();
1040 758 : this->comm().broadcast(subdomain_map);
1041 :
1042 24 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1043 758 : this->comm().broadcast(boundary_info.set_sideset_name_map());
1044 758 : this->comm().broadcast(boundary_info.set_nodeset_name_map());
1045 :
1046 758 : this->comm().broadcast(node_integer_names);
1047 758 : this->comm().broadcast(elem_integer_names);
1048 :
1049 1346 : for (auto & int_name : node_integer_names)
1050 1152 : mesh.add_node_integer(int_name);
1051 :
1052 1052 : for (auto & int_name : elem_integer_names)
1053 576 : mesh.add_elem_integer(int_name);
1054 :
1055 1516 : return input_parallel ? input_n_procs : 0;
1056 710 : }
1057 :
1058 :
1059 :
1060 : template <typename file_id_type>
1061 364 : void CheckpointIO::read_subfile (Xdr & io, bool expect_all_remote)
1062 : {
1063 : // read the nodal locations
1064 364 : this->read_nodes<file_id_type> (io);
1065 :
1066 : // read connectivity
1067 364 : this->read_connectivity<file_id_type> (io);
1068 :
1069 : // read remote_elem connectivity
1070 364 : this->read_remote_elem<file_id_type> (io, expect_all_remote);
1071 :
1072 : // read the boundary conditions
1073 364 : this->read_bcs<file_id_type> (io);
1074 :
1075 : // read the nodesets
1076 364 : this->read_nodesets<file_id_type> (io);
1077 364 : }
1078 :
1079 :
1080 :
1081 : template <typename file_id_type>
1082 168 : void CheckpointIO::read_subdomain_names(Xdr & io)
1083 : {
1084 168 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1085 :
1086 : std::map<subdomain_id_type, std::string> & subdomain_map =
1087 14 : mesh.set_subdomain_name_map();
1088 :
1089 28 : std::vector<file_id_type> subdomain_ids;
1090 168 : subdomain_ids.reserve(subdomain_map.size());
1091 :
1092 42 : std::vector<std::string> subdomain_names;
1093 168 : subdomain_names.reserve(subdomain_map.size());
1094 :
1095 168 : file_id_type n_subdomain_names = 0;
1096 168 : io.data(n_subdomain_names, "# subdomain id to name map");
1097 :
1098 168 : if (n_subdomain_names)
1099 : {
1100 0 : io.data(subdomain_ids);
1101 0 : io.data(subdomain_names);
1102 :
1103 0 : for (auto i : index_range(subdomain_ids))
1104 0 : subdomain_map[cast_int<subdomain_id_type>(subdomain_ids[i])] =
1105 0 : subdomain_names[i];
1106 : }
1107 308 : }
1108 :
1109 :
1110 :
1111 : template <typename file_id_type>
1112 364 : void CheckpointIO::read_nodes (Xdr & io)
1113 : {
1114 : // convenient reference to our mesh
1115 364 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1116 :
1117 : file_id_type n_nodes_here;
1118 364 : io.data(n_nodes_here, "# n_nodes on proc");
1119 :
1120 364 : const bool read_extra_integers = this->version_at_least_1_5();
1121 :
1122 44 : const unsigned int n_extra_integers =
1123 342 : read_extra_integers ? mesh.n_node_integers() : 0;
1124 :
1125 : // Will hold the node id and pid and extra integers
1126 386 : std::vector<file_id_type> id_pid(2 + n_extra_integers);
1127 :
1128 : // For the coordinates
1129 386 : std::vector<Real> coords(LIBMESH_DIM);
1130 :
1131 6134 : for (unsigned int i=0; i<n_nodes_here; i++)
1132 : {
1133 5770 : io.data_stream(id_pid.data(), 2 + n_extra_integers, 2 + n_extra_integers);
1134 :
1135 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1136 5770 : file_id_type unique_id = 0;
1137 5770 : io.data(unique_id, "# unique id");
1138 : #endif
1139 :
1140 5770 : io.data_stream(coords.data(), LIBMESH_DIM, LIBMESH_DIM);
1141 :
1142 394 : Point p;
1143 5770 : p(0) = coords[0];
1144 :
1145 : #if LIBMESH_DIM > 1
1146 5770 : p(1) = coords[1];
1147 : #endif
1148 :
1149 : #if LIBMESH_DIM > 2
1150 5770 : p(2) = coords[2];
1151 : #endif
1152 :
1153 5770 : const dof_id_type id = cast_int<dof_id_type>(id_pid[0]);
1154 :
1155 : // "Wrap around" if we see more processors than we're using.
1156 394 : processor_id_type pid =
1157 6164 : cast_int<processor_id_type>(id_pid[1] % mesh.n_processors());
1158 :
1159 : // If we already have this node (e.g. from another file, when
1160 : // reading multiple distributed CheckpointIO files into a
1161 : // ReplicatedMesh) then we don't want to add it again (because
1162 : // ReplicatedMesh can't handle that) but we do want to assert
1163 : // consistency between what we're reading and what we have.
1164 5770 : const Node * old_node = mesh.query_node_ptr(id);
1165 :
1166 5770 : if (old_node)
1167 : {
1168 60 : libmesh_assert_equal_to(pid, old_node->processor_id());
1169 :
1170 60 : libmesh_assert_equal_to(n_extra_integers, old_node->n_extra_integers());
1171 : #ifndef NDEBUG
1172 60 : for (unsigned int ei=0; ei != n_extra_integers; ++ei)
1173 : {
1174 0 : const dof_id_type extra_int = cast_int<dof_id_type>(id_pid[2+ei]);
1175 0 : libmesh_assert_equal_to(extra_int, old_node->get_extra_integer(ei));
1176 : }
1177 : #endif
1178 :
1179 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1180 60 : libmesh_assert_equal_to(unique_id, old_node->unique_id());
1181 : #endif
1182 : }
1183 : else
1184 : {
1185 : Node * node =
1186 5104 : mesh.add_point(p, id, pid);
1187 :
1188 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1189 5104 : node->set_unique_id(unique_id);
1190 : #endif
1191 :
1192 334 : libmesh_assert_equal_to(n_extra_integers, node->n_extra_integers());
1193 :
1194 11996 : for (unsigned int ei=0; ei != n_extra_integers; ++ei)
1195 : {
1196 7120 : const dof_id_type extra_int = cast_int<dof_id_type>(id_pid[2+ei]);
1197 6892 : node->set_extra_integer(ei, extra_int);
1198 : }
1199 : }
1200 : }
1201 364 : }
1202 :
1203 :
1204 :
1205 : template <typename file_id_type>
1206 364 : void CheckpointIO::read_connectivity (Xdr & io)
1207 : {
1208 : // convenient reference to our mesh
1209 364 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1210 :
1211 364 : const bool read_extra_integers = this->version_at_least_1_5();
1212 364 : const bool read_runtime_topology = this->version_at_least_1_6();
1213 :
1214 44 : const unsigned int n_extra_integers =
1215 342 : read_extra_integers ? mesh.n_elem_integers() : 0;
1216 :
1217 : file_id_type n_elems_here;
1218 364 : io.data(n_elems_here);
1219 :
1220 : // Keep track of the highest dimensional element we've added to the mesh
1221 364 : unsigned int highest_elem_dim = mesh.mesh_dimension();
1222 :
1223 : // RHS: Originally we used invalid_processor_id as a "no parent" tag
1224 : // number, because I'm an idiot. Let's try to support broken files
1225 : // as much as possible.
1226 22 : bool file_is_broken = false;
1227 :
1228 5530 : for (unsigned int i=0; i<n_elems_here; i++)
1229 : {
1230 : // id type pid subdomain_id parent_id
1231 3001 : std::vector<file_id_type> elem_data(6 + n_extra_integers);
1232 627 : io.data_stream
1233 2374 : (elem_data.data(), cast_int<unsigned int>(elem_data.size()),
1234 : cast_int<unsigned int>(elem_data.size()));
1235 :
1236 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1237 2792 : file_id_type unique_id = 0;
1238 2792 : io.data(unique_id, "# unique id");
1239 : #endif
1240 :
1241 : #ifdef LIBMESH_ENABLE_AMR
1242 2792 : uint16_t p_level = 0;
1243 2792 : io.data(p_level, "# p_level");
1244 :
1245 : uint16_t rflag, pflag;
1246 2792 : io.data(rflag, "# rflag");
1247 2792 : io.data(pflag, "# pflag");
1248 : #endif
1249 :
1250 2792 : const ElemType elem_type =
1251 2792 : static_cast<ElemType>(elem_data[1]);
1252 209 : const bool is_c0polygon = (elem_type == C0POLYGON);
1253 209 : const bool is_c0polyhedron = (elem_type == C0POLYHEDRON);
1254 :
1255 2792 : unsigned int n_nodes = Elem::type_to_n_nodes_map[elem_data[1]];
1256 : // Runtime-topology types have no fixed node count in this map.
1257 209 : const bool has_runtime_topology = (n_nodes == invalid_uint);
1258 627 : std::vector<std::vector<unsigned int>> nodes_on_sides;
1259 :
1260 2792 : if (has_runtime_topology)
1261 : {
1262 48 : libmesh_error_msg_if
1263 : (!read_runtime_topology,
1264 : "Checkpoint format 1.6 or newer is required to read " <<
1265 : Utility::enum_to_string(elem_type) << " elements.");
1266 :
1267 8 : std::vector<file_id_type> runtime_topology;
1268 48 : io.data(runtime_topology, "# runtime topology");
1269 52 : libmesh_error_msg_if(runtime_topology.size() < 2,
1270 : "Invalid runtime element topology.");
1271 :
1272 4 : std::size_t topology_index = 0;
1273 44 : n_nodes =
1274 48 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1275 : const unsigned int n_sides =
1276 48 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1277 48 : nodes_on_sides.resize(n_sides);
1278 :
1279 360 : for (auto s : index_range(nodes_on_sides))
1280 : {
1281 338 : libmesh_error_msg_if
1282 : (topology_index == runtime_topology.size(),
1283 : "Incomplete runtime element checkpoint topology.");
1284 :
1285 : const unsigned int n_side_nodes =
1286 312 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1287 312 : libmesh_error_msg_if
1288 : (n_side_nodes > runtime_topology.size() - topology_index,
1289 : "Invalid runtime element side checkpoint topology.");
1290 :
1291 52 : auto & side_nodes = nodes_on_sides[s];
1292 312 : side_nodes.resize(n_side_nodes);
1293 1416 : for (auto n : index_range(side_nodes))
1294 : {
1295 : const unsigned int local_node =
1296 1196 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1297 1104 : libmesh_error_msg_if
1298 : (local_node >= n_nodes,
1299 : "Runtime element side checkpoint topology references "
1300 : "an invalid local node.");
1301 1104 : side_nodes[n] = local_node;
1302 : }
1303 : }
1304 :
1305 52 : libmesh_error_msg_if
1306 : (topology_index != runtime_topology.size(),
1307 : "Extra data in runtime element checkpoint topology.");
1308 :
1309 48 : if (is_c0polygon)
1310 : {
1311 24 : libmesh_error_msg_if
1312 : (n_nodes < 3 || n_sides != n_nodes,
1313 : "Invalid C0POLYGON checkpoint topology.");
1314 144 : for (const auto & side_nodes : nodes_on_sides)
1315 130 : libmesh_error_msg_if
1316 : (side_nodes.size() != 2,
1317 : "Invalid C0POLYGON side checkpoint topology.");
1318 : }
1319 24 : else if (is_c0polyhedron)
1320 : {
1321 24 : libmesh_error_msg_if(n_sides < 4,
1322 : "Invalid C0POLYHEDRON checkpoint topology.");
1323 216 : for (const auto & side_nodes : nodes_on_sides)
1324 208 : libmesh_error_msg_if
1325 : (side_nodes.size() < 3,
1326 : "Invalid C0POLYHEDRON side checkpoint topology.");
1327 : }
1328 : }
1329 :
1330 : // Snag the node ids this element was connected to
1331 3001 : std::vector<file_id_type> conn_data(n_nodes);
1332 627 : io.data_stream
1333 2374 : (conn_data.data(), cast_int<unsigned int>(conn_data.size()),
1334 : cast_int<unsigned int>(conn_data.size()));
1335 :
1336 : const dof_id_type id =
1337 2792 : cast_int<dof_id_type> (elem_data[0]);
1338 209 : const processor_id_type proc_id =
1339 : cast_int<processor_id_type>
1340 3001 : (elem_data[2] % mesh.n_processors());
1341 209 : const subdomain_id_type subdomain_id =
1342 2792 : restrict_int<subdomain_id_type>(elem_data[3]);
1343 :
1344 : // Old broken files used processsor_id_type(-1)...
1345 : // But we *know* our first element will be level 0
1346 2792 : if (i == 0 && elem_data[4] == 65535)
1347 0 : file_is_broken = true;
1348 :
1349 : // On a broken file we can't tell whether a parent of 65535 is a
1350 : // null parent or an actual parent of 65535. Assuming the
1351 : // former will cause less breakage.
1352 2583 : Elem * parent =
1353 209 : (elem_data[4] == static_cast<largest_id_type>(-1) ||
1354 2792 : (file_is_broken && elem_data[4] == 65535)) ?
1355 0 : nullptr : mesh.elem_ptr(cast_int<dof_id_type>(elem_data[4]));
1356 :
1357 209 : const unsigned short int child_num =
1358 209 : (elem_data[5] == static_cast<largest_id_type>(-1) ||
1359 2792 : (file_is_broken && elem_data[5] == 65535)) ?
1360 : static_cast<unsigned short>(-1) :
1361 0 : cast_int<unsigned short>(elem_data[5]);
1362 :
1363 209 : if (!parent)
1364 209 : libmesh_assert_equal_to
1365 : (child_num, static_cast<unsigned short>(-1));
1366 :
1367 2792 : Elem * old_elem = mesh.query_elem_ptr(id);
1368 :
1369 : // If we already have this element (e.g. from another file,
1370 : // when reading multiple distributed CheckpointIO files into
1371 : // a ReplicatedMesh) then we don't want to add it again
1372 : // (because ReplicatedMesh can't handle that) but we do want
1373 : // to assert consistency between what we're reading and what
1374 : // we have.
1375 2792 : if (old_elem)
1376 : {
1377 32 : libmesh_assert_equal_to(elem_type, old_elem->type());
1378 32 : libmesh_assert_equal_to(proc_id, old_elem->processor_id());
1379 32 : libmesh_assert_equal_to(subdomain_id, old_elem->subdomain_id());
1380 32 : if (parent)
1381 0 : libmesh_assert_equal_to(parent, old_elem->parent());
1382 : else
1383 32 : libmesh_assert(!old_elem->parent());
1384 :
1385 32 : libmesh_assert_equal_to(n_extra_integers, old_elem->n_extra_integers());
1386 : #ifndef NDEBUG
1387 32 : for (unsigned int ei=0; ei != n_extra_integers; ++ei)
1388 : {
1389 0 : const dof_id_type extra_int = cast_int<dof_id_type>(elem_data[6+ei]);
1390 0 : libmesh_assert_equal_to(extra_int, old_elem->get_extra_integer(ei));
1391 : }
1392 : #endif
1393 :
1394 32 : libmesh_assert_equal_to(old_elem->n_nodes(), conn_data.size());
1395 :
1396 128 : for (unsigned int n=0,
1397 32 : n_conn = cast_int<unsigned int>(conn_data.size());
1398 160 : n != n_conn; n++)
1399 128 : libmesh_assert_equal_to
1400 : (old_elem->node_id(n),
1401 : cast_int<dof_id_type>(conn_data[n]));
1402 :
1403 32 : if (has_runtime_topology)
1404 : {
1405 0 : libmesh_assert_equal_to(old_elem->n_sides(),
1406 : nodes_on_sides.size());
1407 : #ifndef NDEBUG
1408 0 : for (auto s : index_range(nodes_on_sides))
1409 0 : libmesh_assert(old_elem->nodes_on_side(s) ==
1410 : nodes_on_sides[s]);
1411 : #endif
1412 : }
1413 : }
1414 : else
1415 : {
1416 : // Create the element
1417 2259 : std::unique_ptr<Elem> elem;
1418 2436 : std::unique_ptr<Node> generated_mid_elem_node;
1419 :
1420 2436 : if (is_c0polygon)
1421 26 : elem = std::make_unique<C0Polygon>(n_nodes, parent);
1422 2412 : else if (is_c0polyhedron)
1423 : {
1424 30 : std::vector<std::shared_ptr<Polygon>> sides(nodes_on_sides.size());
1425 216 : for (auto s : index_range(nodes_on_sides))
1426 : {
1427 32 : const auto & side_node_indices = nodes_on_sides[s];
1428 16 : auto side =
1429 : std::make_shared<C0Polygon>
1430 208 : (cast_int<unsigned int>(side_node_indices.size()));
1431 1056 : for (auto n : index_range(side_node_indices))
1432 288 : side->set_node
1433 720 : (n, mesh.node_ptr(cast_int<dof_id_type>
1434 936 : (conn_data[side_node_indices[n]])));
1435 32 : sides[s] = std::move(side);
1436 : }
1437 :
1438 24 : elem = std::make_unique<C0Polyhedron>
1439 : (sides, generated_mid_elem_node, parent);
1440 :
1441 24 : libmesh_error_msg_if
1442 : (elem->n_nodes() != conn_data.size(),
1443 : "C0POLYHEDRON checkpoint topology is incompatible with "
1444 : "this libMesh configuration.");
1445 :
1446 312 : for (auto n : make_range(elem->n_vertices()))
1447 360 : libmesh_error_msg_if
1448 : (elem->node_id(n) !=
1449 : cast_int<dof_id_type>(conn_data[n]),
1450 : "C0POLYHEDRON checkpoint topology has inconsistent "
1451 : "local node ordering.");
1452 20 : }
1453 : else
1454 4603 : elem = Elem::build(elem_type, parent);
1455 :
1456 2436 : if (has_runtime_topology)
1457 : {
1458 48 : libmesh_error_msg_if
1459 : (!elem->runtime_topology() ||
1460 : elem->n_nodes() != conn_data.size() ||
1461 : elem->n_sides() != nodes_on_sides.size(),
1462 : Utility::enum_to_string(elem_type) <<
1463 : " checkpoint topology is incompatible with this "
1464 : "libMesh configuration.");
1465 :
1466 360 : for (auto s : index_range(nodes_on_sides))
1467 624 : libmesh_error_msg_if
1468 : (elem->nodes_on_side(s) != nodes_on_sides[s],
1469 : Utility::enum_to_string(elem_type) <<
1470 : " checkpoint topology has inconsistent side ordering.");
1471 : }
1472 :
1473 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1474 2436 : elem->set_unique_id(unique_id);
1475 : #endif
1476 :
1477 2436 : if (elem->dim() > highest_elem_dim)
1478 0 : highest_elem_dim = elem->dim();
1479 :
1480 2436 : elem->set_id() = id;
1481 2436 : elem->processor_id() = proc_id;
1482 2436 : elem->subdomain_id() = subdomain_id;
1483 :
1484 : #ifdef LIBMESH_ENABLE_AMR
1485 2436 : elem->hack_p_level(p_level);
1486 :
1487 2436 : elem->set_refinement_flag (cast_int<Elem::RefinementState>(rflag));
1488 2436 : elem->set_p_refinement_flag(cast_int<Elem::RefinementState>(pflag));
1489 :
1490 : // Set parent connections
1491 2436 : if (parent)
1492 : {
1493 : // We must specify a child_num, because we will have
1494 : // skipped adding any preceding remote_elem children
1495 0 : parent->add_child(elem.get(), child_num);
1496 : }
1497 : #else
1498 : libmesh_ignore(child_num);
1499 : #endif
1500 :
1501 177 : libmesh_assert(elem->n_nodes() == conn_data.size());
1502 :
1503 : // Connect all the nodes to this element
1504 11536 : for (unsigned int n=0,
1505 354 : n_conn = cast_int<unsigned int>(conn_data.size());
1506 13972 : n != n_conn; n++)
1507 13046 : elem->set_node(n,
1508 11536 : mesh.node_ptr(cast_int<dof_id_type>(conn_data[n])));
1509 :
1510 2790 : Elem * added_elem = mesh.add_elem(std::move(elem));
1511 :
1512 177 : libmesh_assert_equal_to(n_extra_integers, added_elem->n_extra_integers());
1513 3531 : for (unsigned int ei=0; ei != n_extra_integers; ++ei)
1514 : {
1515 1134 : const dof_id_type extra_int = cast_int<dof_id_type>(elem_data[6+ei]);
1516 1095 : added_elem->set_extra_integer(ei, extra_int);
1517 : }
1518 2082 : }
1519 : }
1520 :
1521 364 : mesh.set_mesh_dimension(cast_int<unsigned char>(highest_elem_dim));
1522 364 : }
1523 :
1524 :
1525 : template <typename file_id_type>
1526 364 : void CheckpointIO::read_remote_elem (Xdr & io, bool libmesh_dbg_var(expect_all_remote))
1527 : {
1528 : // convenient reference to our mesh
1529 364 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1530 :
1531 : // Find the remote_elem neighbor links
1532 44 : std::vector<file_id_type> elem_ids;
1533 44 : std::vector<uint16_t> elem_sides;
1534 :
1535 364 : io.data(elem_ids, "# remote neighbor elem_ids");
1536 364 : io.data(elem_sides, "# remote neighbor elem_sides");
1537 :
1538 22 : libmesh_assert_equal_to(elem_ids.size(), elem_sides.size());
1539 :
1540 1072 : for (auto i : index_range(elem_ids))
1541 : {
1542 836 : Elem & elem = mesh.elem_ref(cast_int<dof_id_type>(elem_ids[i]));
1543 836 : if (!elem.neighbor_ptr(elem_sides[i]))
1544 708 : elem.set_neighbor(elem_sides[i],
1545 : const_cast<RemoteElem *>(remote_elem));
1546 : else
1547 0 : libmesh_assert(!expect_all_remote);
1548 : }
1549 :
1550 : // Find the remote_elem children links
1551 44 : std::vector<file_id_type> parent_ids;
1552 44 : std::vector<uint16_t> child_numbers;
1553 :
1554 364 : io.data(parent_ids, "# remote child parent_ids");
1555 364 : io.data(child_numbers, "# remote child_numbers");
1556 :
1557 : #ifdef LIBMESH_ENABLE_AMR
1558 364 : for (auto i : index_range(parent_ids))
1559 : {
1560 0 : Elem & elem = mesh.elem_ref(cast_int<dof_id_type>(parent_ids[i]));
1561 :
1562 : // We'd like to assert that no child pointer already exists to
1563 : // be overwritten by remote_elem, but Elem doesn't actually have
1564 : // an API that will return a child pointer without asserting
1565 : // that it isn't nullptr.
1566 0 : const Elem * child = elem.raw_child_ptr(child_numbers[i]);
1567 :
1568 0 : if (!child)
1569 0 : elem.add_child(const_cast<RemoteElem *>(remote_elem),
1570 0 : child_numbers[i]);
1571 : else
1572 0 : libmesh_assert(!expect_all_remote);
1573 : }
1574 : #endif
1575 364 : }
1576 :
1577 :
1578 :
1579 : template <typename file_id_type>
1580 364 : void CheckpointIO::read_bcs (Xdr & io)
1581 : {
1582 : // convenient reference to our mesh
1583 364 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1584 :
1585 : // and our boundary info object
1586 22 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1587 :
1588 44 : std::vector<file_id_type> element_id_list;
1589 44 : std::vector<uint16_t> side_list;
1590 44 : std::vector<file_id_type> bc_id_list;
1591 :
1592 364 : io.data(element_id_list, "# element ids for bcs");
1593 364 : io.data(side_list, "# sides of elements for bcs");
1594 364 : io.data(bc_id_list, "# bc ids");
1595 :
1596 3532 : for (auto i : index_range(element_id_list))
1597 558 : boundary_info.add_side
1598 3354 : (cast_int<dof_id_type>(element_id_list[i]), side_list[i],
1599 3354 : cast_int<boundary_id_type>(bc_id_list[i]));
1600 364 : }
1601 :
1602 :
1603 :
1604 : template <typename file_id_type>
1605 364 : void CheckpointIO::read_nodesets (Xdr & io)
1606 : {
1607 : // convenient reference to our mesh
1608 364 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1609 :
1610 : // and our boundary info object
1611 22 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1612 :
1613 44 : std::vector<file_id_type> node_id_list;
1614 44 : std::vector<file_id_type> bc_id_list;
1615 :
1616 364 : io.data(node_id_list, "# node id list");
1617 364 : io.data(bc_id_list, "# nodeset bc id list");
1618 :
1619 5648 : for (auto i : index_range(node_id_list))
1620 792 : boundary_info.add_node
1621 5284 : (cast_int<dof_id_type>(node_id_list[i]),
1622 5548 : cast_int<boundary_id_type>(bc_id_list[i]));
1623 364 : }
1624 :
1625 :
1626 :
1627 : template <typename file_id_type>
1628 336 : void CheckpointIO::read_bc_names(Xdr & io, BoundaryInfo & info, bool is_sideset)
1629 : {
1630 336 : std::map<boundary_id_type, std::string> & boundary_map = is_sideset ?
1631 28 : info.set_sideset_name_map() : info.set_nodeset_name_map();
1632 :
1633 56 : std::vector<file_id_type> boundary_ids;
1634 84 : std::vector<std::string> boundary_names;
1635 :
1636 336 : file_id_type n_boundary_names = 0;
1637 :
1638 336 : if (is_sideset)
1639 168 : io.data(n_boundary_names, "# sideset id to name map");
1640 : else
1641 168 : io.data(n_boundary_names, "# nodeset id to name map");
1642 :
1643 336 : if (n_boundary_names)
1644 : {
1645 240 : io.data(boundary_ids);
1646 240 : io.data(boundary_names);
1647 : }
1648 :
1649 : // Add them back into the map
1650 1296 : for (auto i : index_range(boundary_ids))
1651 2000 : boundary_map[cast_int<boundary_id_type>(boundary_ids[i])] =
1652 160 : boundary_names[i];
1653 616 : }
1654 :
1655 :
1656 : template <typename file_id_type>
1657 168 : void CheckpointIO::read_integers_names
1658 : (Xdr & io,
1659 : std::vector<std::string> & node_integer_names,
1660 : std::vector<std::string> & elem_integer_names)
1661 : {
1662 : file_id_type n_node_integers, n_elem_integers;
1663 :
1664 168 : io.data(n_node_integers, "# n_extra_integers per node");
1665 168 : io.data(node_integer_names);
1666 168 : io.data(n_elem_integers, "# n_extra_integers per elem");
1667 168 : io.data(elem_integer_names);
1668 168 : }
1669 :
1670 :
1671 0 : unsigned int CheckpointIO::n_active_levels_in(MeshBase::const_element_iterator begin,
1672 : MeshBase::const_element_iterator end) const
1673 : {
1674 0 : unsigned int max_level = 0;
1675 :
1676 0 : for (const auto & elem : as_range(begin, end))
1677 0 : max_level = std::max(elem->level(), max_level);
1678 :
1679 0 : return max_level + 1;
1680 : }
1681 :
1682 : } // namespace libMesh
|
