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