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 424 : std::string_view extension(std::string_view s)
91 : {
92 310 : auto pos = s.rfind(".");
93 424 : if (pos == std::string::npos)
94 0 : return "";
95 424 : return s.substr(pos, s.size() - pos);
96 : }
97 :
98 526 : std::string split_dir(const std::string & input_name, libMesh::processor_id_type n_procs)
99 : {
100 910 : return input_name + "/" + std::to_string(n_procs);
101 : }
102 :
103 :
104 244 : std::string header_file(const std::string & input_name, libMesh::processor_id_type n_procs)
105 : {
106 712 : return (split_dir(input_name, n_procs) + "/header").append(extension(input_name));
107 : }
108 :
109 : std::string
110 180 : split_file(const std::string & input_name,
111 : libMesh::processor_id_type n_procs,
112 : libMesh::processor_id_type proc_id)
113 : {
114 528 : return (split_dir(input_name, n_procs) + "/split-" + std::to_string(n_procs) + "-" +
115 624 : std::to_string(proc_id)).append(extension(input_name));
116 : }
117 :
118 86 : void make_dir(const std::string & input_name, libMesh::processor_id_type n_procs)
119 : {
120 86 : 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 86 : libmesh_error_msg_if
123 : (ret != 0 && ret != -1,
124 : "Failed to create mesh split directory '" << input_name << "': " << std::strerror(ret));
125 :
126 110 : auto dir_name = split_dir(input_name, n_procs);
127 86 : ret = libMesh::Utility::mkdir(dir_name.c_str());
128 86 : if (ret == -1)
129 : libmesh_warning("In CheckpointIO::write, directory '"
130 : << dir_name << "' already exists, overwriting contents.");
131 : else
132 30 : libmesh_error_msg_if
133 : (ret != 0, "Failed to create mesh split directory '" << dir_name << "': " << std::strerror(ret));
134 86 : }
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 158 : 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 70 : _binary (binary_in),
172 70 : _parallel (false),
173 0 : _version ("checkpoint-1.6"),
174 246 : _my_processor_ids (1, processor_id()),
175 272 : _my_n_processors (mesh.is_replicated() ? 1 : n_processors())
176 : {
177 158 : }
178 :
179 14 : 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 6 : _binary (binary_in),
184 6 : _parallel (false),
185 0 : _version ("checkpoint-1.6"),
186 22 : _my_processor_ids (1, processor_id()),
187 24 : _my_n_processors (mesh.is_replicated() ? 1 : n_processors())
188 : {
189 14 : }
190 :
191 76 : CheckpointIO::~CheckpointIO () = default;
192 :
193 86 : 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 110 : if (this->processor_id() == 0)
199 : {
200 98 : header_name = header_file(input_name, _my_n_processors);
201 :
202 : {
203 : // look for header+splits with nprocs equal to _my_n_processors
204 84 : std::ifstream in (header_name.c_str());
205 56 : 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 28 : }
220 :
221 154 : 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 56 : io.data(input_version);
226 :
227 : // read the data type
228 56 : io.data (data_size);
229 28 : }
230 :
231 86 : this->comm().broadcast(data_size);
232 86 : this->comm().broadcast(header_name);
233 :
234 : // How many per-processor files are here?
235 : largest_id_type input_n_procs;
236 :
237 86 : 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 86 : case 8:
245 86 : 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 110 : return cast_int<processor_id_type>(input_n_procs);
254 : }
255 :
256 16 : void CheckpointIO::cleanup(const std::string & input_name, processor_id_type n_procs)
257 : {
258 20 : auto header = header_file(input_name, n_procs);
259 16 : 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 32 : for (processor_id_type i = 0; i < n_procs; i++)
264 : {
265 20 : auto split = split_file(input_name, n_procs, i);
266 16 : 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 20 : auto dir = split_dir(input_name, n_procs);
272 16 : 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 16 : rmdir(input_name.c_str());
279 16 : }
280 :
281 :
282 440 : bool CheckpointIO::version_at_least_1_5() const
283 : {
284 : return
285 880 : (this->version().find("1.5") != std::string::npos) ||
286 556 : (this->version().find("1.6") != std::string::npos);
287 : }
288 :
289 :
290 164 : bool CheckpointIO::version_at_least_1_6() const
291 : {
292 164 : return (this->version().find("1.6") != std::string::npos);
293 : }
294 :
295 :
296 86 : 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 86 : _parallel = _parallel || !mesh.is_serial();
307 :
308 24 : processor_id_type use_n_procs = 1;
309 86 : if (_parallel)
310 56 : use_n_procs = _my_n_processors;
311 :
312 110 : std::string header_file_name = header_file(name, use_n_procs);
313 86 : 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 110 : if (this->processor_id() == 0)
318 : {
319 154 : Xdr io (header_file_name, this->binary() ? ENCODE : WRITE);
320 :
321 : // write the version
322 56 : io.data(_version, "# version");
323 :
324 : // write what kind of data type we're using
325 56 : header_id_type data_size = sizeof(largest_id_type);
326 56 : 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 56 : uint16_t mesh_dimension = cast_int<uint16_t>(mesh.mesh_dimension());
332 56 : io.data(mesh_dimension, "# dimensions");
333 : }
334 :
335 : // Write out whether or not this is serial output
336 : {
337 56 : uint16_t parallel = _parallel;
338 56 : 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 56 : if (_parallel)
344 : {
345 32 : largest_id_type n_procs = _my_n_processors;
346 32 : io.data(n_procs, "# n_procs");
347 : }
348 :
349 : // write subdomain names
350 56 : this->write_subdomain_names(io);
351 :
352 : // write boundary id names
353 14 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
354 56 : write_bc_names(io, boundary_info, true); // sideset names
355 56 : write_bc_names(io, boundary_info, false); // nodeset names
356 :
357 : // write extra integer names
358 56 : const bool write_extra_integers = this->version_at_least_1_5();
359 :
360 56 : if (write_extra_integers)
361 : {
362 56 : largest_id_type n_node_integers = mesh.n_node_integers();
363 56 : io.data(n_node_integers, "# n_extra_integers per node");
364 :
365 42 : std::vector<std::string> node_integer_names;
366 100 : for (unsigned int i=0; i != n_node_integers; ++i)
367 44 : node_integer_names.push_back(mesh.get_node_integer_name(i));
368 :
369 56 : io.data(node_integer_names);
370 :
371 56 : largest_id_type n_elem_integers = mesh.n_elem_integers();
372 56 : io.data(n_elem_integers, "# n_extra_integers per elem");
373 :
374 28 : std::vector<std::string> elem_integer_names;
375 78 : for (unsigned int i=0; i != n_elem_integers; ++i)
376 22 : elem_integer_names.push_back(mesh.get_elem_integer_name(i));
377 :
378 56 : io.data(elem_integer_names);
379 28 : }
380 :
381 :
382 28 : }
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 86 : if (_parallel)
400 : {
401 56 : ids_to_write = _my_processor_ids;
402 114 : 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 1248 : for (auto & elem : mesh.element_ptr_range())
406 : {
407 800 : const processor_id_type p = elem->processor_id();
408 1024 : if (auto eop_it = elements_on_pid.find(p);
409 224 : eop_it != eop_end)
410 512 : eop_it->second.push_back(elem);
411 24 : }
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 110 : bc_triples = boundary_info.build_side_list();
431 : std::vector<std::tuple<dof_id_type, boundary_id_type>>
432 110 : bc_tuples = boundary_info.build_node_list();
433 :
434 168 : for (const auto & my_pid : ids_to_write)
435 : {
436 104 : auto file_name = split_file(name, use_n_procs, my_pid);
437 227 : 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 82 : if (!_parallel || !mesh.is_serial())
444 90 : 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 90 : for (processor_id_type p : {my_pid, DofObject::invalid_processor_id})
461 : {
462 60 : 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 30 : v_t * elempp = p_elements.data();
472 30 : v_t * elemend = elempp + p_elements.size();
473 :
474 : const MeshBase::const_element_iterator
475 : pid_elements_begin = MeshBase::const_element_iterator
476 38 : (elempp, elemend, Predicates::NotNull<v_t *>()),
477 : pid_elements_end = MeshBase::const_element_iterator
478 38 : (elemend, elemend, Predicates::NotNull<v_t *>()),
479 : active_pid_elements_begin = MeshBase::const_element_iterator
480 38 : (elempp, elemend, Predicates::Active<v_t *>()),
481 : active_pid_elements_end = MeshBase::const_element_iterator
482 52 : (elemend, elemend, Predicates::Active<v_t *>());
483 :
484 : query_ghosting_functors
485 52 : (mesh, p, active_pid_elements_begin,
486 : active_pid_elements_end, elements);
487 52 : connect_children(mesh, pid_elements_begin,
488 : pid_elements_end, elements);
489 : }
490 : }
491 : }
492 :
493 44 : connected_node_set_type connected_nodes;
494 82 : connect_element_dependencies(mesh, elements, connected_nodes);
495 :
496 : // write the nodal locations
497 82 : this->write_nodes (io, connected_nodes);
498 :
499 : // write connectivity
500 82 : this->write_connectivity (io, elements);
501 :
502 : // write remote_elem connectivity
503 82 : this->write_remote_elem (io, elements);
504 :
505 : // write the boundary condition information
506 82 : this->write_bcs (io, elements, bc_triples);
507 :
508 : // write the nodeset information
509 82 : this->write_nodesets (io, connected_nodes, bc_tuples);
510 :
511 : // close it up
512 82 : io.close();
513 38 : }
514 :
515 : // this->comm().barrier();
516 86 : }
517 :
518 56 : 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 70 : std::vector<largest_id_type> subdomain_ids; subdomain_ids.reserve(subdomain_map.size());
526 84 : 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 56 : largest_id_type n_subdomain_names = 0;
532 56 : 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 56 : io.data(n_subdomain_names, "# subdomain id to name map");
541 : // Write out the ids and names in two vectors
542 56 : if (n_subdomain_names)
543 : {
544 0 : io.data(subdomain_ids);
545 0 : io.data(subdomain_names);
546 : }
547 28 : }
548 56 : }
549 :
550 :
551 :
552 82 : void CheckpointIO::write_nodes (Xdr & io,
553 : const connected_node_set_type & nodeset) const
554 : {
555 82 : largest_id_type n_nodes_here = nodeset.size();
556 :
557 82 : io.data(n_nodes_here, "# n_nodes on proc");
558 :
559 82 : const bool write_extra_integers = this->version_at_least_1_5();
560 : const unsigned int n_extra_integers =
561 82 : write_extra_integers ? MeshOutput<MeshBase>::mesh().n_node_integers() : 0;
562 :
563 : // Will hold the node id and pid and extra integers
564 104 : std::vector<largest_id_type> id_pid(2 + n_extra_integers);
565 :
566 : // For the coordinates
567 104 : std::vector<Real> coords(LIBMESH_DIM);
568 :
569 1568 : for (const auto & node : nodeset)
570 : {
571 1486 : id_pid[0] = node->id();
572 1486 : id_pid[1] = node->processor_id();
573 :
574 394 : libmesh_assert_equal_to(n_extra_integers, node->n_extra_integers());
575 2322 : for (unsigned int i=0; i != n_extra_integers; ++i)
576 1292 : id_pid[2+i] = node->get_extra_integer(i);
577 :
578 1486 : io.data_stream(id_pid.data(), 2 + n_extra_integers, 2 + n_extra_integers);
579 :
580 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
581 1486 : largest_id_type unique_id = node->unique_id();
582 :
583 1486 : io.data(unique_id, "# unique id");
584 : #endif
585 :
586 1486 : coords[0] = (*node)(0);
587 :
588 : #if LIBMESH_DIM > 1
589 1486 : coords[1] = (*node)(1);
590 : #endif
591 :
592 : #if LIBMESH_DIM > 2
593 1486 : coords[2] = (*node)(2);
594 : #endif
595 :
596 1486 : io.data_stream(coords.data(), LIBMESH_DIM, 3);
597 : }
598 82 : }
599 :
600 :
601 :
602 82 : void CheckpointIO::write_connectivity (Xdr & io,
603 : const std::set<const Elem *, CompareElemIdsByLevel> & elements) const
604 : {
605 22 : libmesh_assert (io.writing());
606 :
607 82 : const bool write_extra_integers = this->version_at_least_1_5();
608 82 : const bool write_runtime_topology = this->version_at_least_1_6();
609 : const unsigned int n_extra_integers =
610 82 : 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 104 : 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 82 : largest_id_type n_elems_here = elements.size();
619 :
620 82 : io.data(n_elems_here, "# number of elements");
621 :
622 870 : for (const auto & elem : elements)
623 : {
624 788 : unsigned int n_nodes = elem->n_nodes();
625 :
626 788 : elem_data[0] = elem->id();
627 788 : elem_data[1] = elem->type();
628 788 : elem_data[2] = elem->processor_id();
629 788 : elem_data[3] = elem->subdomain_id();
630 :
631 : #ifdef LIBMESH_ENABLE_AMR
632 997 : 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 788 : elem_data[4] = static_cast<largest_id_type>(-1);
641 788 : elem_data[5] = static_cast<largest_id_type>(-1);
642 : }
643 :
644 931 : for (unsigned int i=0; i != n_extra_integers; ++i)
645 221 : elem_data[6+i] = elem->get_extra_integer(i);
646 :
647 788 : conn_data.resize(n_nodes);
648 :
649 4128 : for (unsigned int i=0; i<n_nodes; i++)
650 5106 : conn_data[i] = elem->node_id(i);
651 :
652 997 : 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 788 : largest_id_type unique_id = elem->unique_id();
658 :
659 788 : io.data(unique_id, "# unique id");
660 : #endif
661 :
662 : #ifdef LIBMESH_ENABLE_AMR
663 997 : uint16_t p_level = cast_int<uint16_t>(elem->p_level());
664 788 : io.data(p_level, "# p_level");
665 :
666 788 : uint16_t rflag = elem->refinement_flag();
667 788 : io.data(rflag, "# rflag");
668 :
669 788 : uint16_t pflag = elem->p_refinement_flag();
670 788 : io.data(pflag, "# pflag");
671 : #endif
672 :
673 788 : if (elem->runtime_topology())
674 : {
675 16 : 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 16 : runtime_topology.push_back(n_nodes);
682 16 : runtime_topology.push_back(elem->n_sides());
683 124 : for (auto s : elem->side_index_range())
684 : {
685 130 : const auto side_nodes = elem->nodes_on_side(s);
686 130 : runtime_topology.push_back(side_nodes.size());
687 472 : for (const auto n : side_nodes)
688 368 : runtime_topology.push_back(n);
689 : }
690 :
691 16 : io.data(runtime_topology, "# runtime topology");
692 : }
693 :
694 997 : io.data_stream(conn_data.data(),
695 : cast_int<unsigned int>(conn_data.size()),
696 : cast_int<unsigned int>(conn_data.size()));
697 : }
698 82 : }
699 :
700 :
701 82 : 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 870 : for (const auto & elem : elements)
711 : {
712 4213 : for (auto n : elem->side_index_range())
713 : {
714 3216 : const Elem * neigh = elem->neighbor_ptr(n);
715 4036 : if (neigh == remote_elem ||
716 820 : (neigh && !elements.count(neigh)))
717 : {
718 208 : elem_ids.push_back(elem->id());
719 208 : elem_sides.push_back(n);
720 : }
721 : }
722 :
723 : #ifdef LIBMESH_ENABLE_AMR
724 788 : 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 82 : io.data(elem_ids, "# remote neighbor elem_ids");
743 82 : io.data(elem_sides, "# remote neighbor elem_sides");
744 82 : io.data(parent_ids, "# remote child parent_ids");
745 82 : io.data(child_numbers, "# remote child_numbers");
746 82 : }
747 :
748 :
749 :
750 82 : 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 82 : element_id_list.reserve(bc_size);
764 82 : side_list.reserve(bc_size);
765 82 : bc_id_list.reserve(bc_size);
766 :
767 22 : std::unordered_set<dof_id_type> elems;
768 870 : for (auto & e : elements)
769 788 : elems.insert(e->id());
770 :
771 970 : for (const auto & t : bc_triples)
772 468 : if (elems.count(std::get<0>(t)))
773 : {
774 720 : element_id_list.push_back(std::get<0>(t));
775 720 : side_list.push_back(std::get<1>(t));
776 720 : bc_id_list.push_back(std::get<2>(t));
777 : }
778 :
779 :
780 82 : io.data(element_id_list, "# element ids for bcs");
781 82 : io.data(side_list, "# sides of elements for bcs");
782 82 : io.data(bc_id_list, "# bc ids");
783 82 : }
784 :
785 :
786 :
787 82 : 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 82 : node_id_list.reserve(nodeset_size);
803 82 : bc_id_list.reserve(nodeset_size);
804 :
805 1298 : for (const auto & t : bc_tuples)
806 1272 : if (nodeset.count(mesh.node_ptr(std::get<0>(t))))
807 : {
808 1020 : node_id_list.push_back(std::get<0>(t));
809 1020 : bc_id_list.push_back(std::get<1>(t));
810 : }
811 :
812 82 : io.data(node_id_list, "# node id list");
813 82 : io.data(bc_id_list, "# nodeset bc id list");
814 82 : }
815 :
816 :
817 :
818 112 : void CheckpointIO::write_bc_names (Xdr & io, const BoundaryInfo & info, bool is_sideset) const
819 : {
820 112 : 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 140 : std::vector<largest_id_type> boundary_ids; boundary_ids.reserve(boundary_map.size());
824 168 : 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 112 : largest_id_type n_boundary_names = 0;
830 432 : for (const auto & [id, name] : boundary_map)
831 320 : if (!name.empty())
832 : {
833 320 : n_boundary_names++;
834 320 : boundary_ids.push_back(id);
835 320 : boundary_names.push_back(name);
836 : }
837 :
838 112 : if (is_sideset)
839 56 : io.data(n_boundary_names, "# sideset id to name map");
840 : else
841 56 : io.data(n_boundary_names, "# nodeset id to name map");
842 : // Write out the ids and names in two vectors
843 112 : if (n_boundary_names)
844 : {
845 80 : io.data(boundary_ids);
846 80 : io.data(boundary_names);
847 : }
848 168 : }
849 :
850 86 : void CheckpointIO::read (const std::string & input_name)
851 : {
852 48 : LOG_SCOPE("read()","CheckpointIO");
853 :
854 86 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
855 :
856 24 : libmesh_assert(!mesh.n_elem());
857 :
858 : header_id_type data_size;
859 86 : processor_id_type input_n_procs = select_split_config(input_name, data_size);
860 110 : 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 86 : 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 68 : (input_parallel && !mesh.is_replicated()) ?
874 16 : mesh.processor_id() : 0;
875 : const processor_id_type stride =
876 68 : (input_parallel && !mesh.is_replicated()) ?
877 16 : mesh.n_processors() : 1;
878 :
879 150 : for (processor_id_type proc_id = begin_proc_id; proc_id < input_n_procs;
880 82 : proc_id = cast_int<processor_id_type>(proc_id + stride))
881 : {
882 104 : auto file_name = split_file(input_name, input_n_procs, proc_id);
883 :
884 : {
885 126 : std::ifstream in (file_name.c_str());
886 :
887 82 : libmesh_error_msg_if(!in.good(), "ERROR: cannot locate specified file:\n\t" << file_name);
888 38 : }
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 182 : (input_n_procs <= mesh.n_processors() &&
895 78 : !mesh.is_replicated());
896 :
897 167 : Xdr io (file_name, this->binary() ? DECODE : READ);
898 :
899 82 : 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 82 : case 8:
907 82 : this->read_subfile<uint64_t>(io, expect_all_remote);
908 22 : break;
909 0 : default:
910 0 : libmesh_error();
911 : }
912 :
913 82 : io.close();
914 38 : }
915 : }
916 :
917 : // If the mesh was only read on processor 0 then we need to broadcast it
918 86 : if (mesh.is_replicated())
919 58 : MeshCommunication().broadcast(mesh);
920 : // If the mesh is really distributed then we need to make sure it
921 : // knows that
922 36 : else if (mesh.n_processors() > 1)
923 24 : 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 86 : if (mesh.skip_partitioning())
928 : {
929 0 : mesh.recalculate_n_partitions();
930 0 : mesh.update_post_partitioning();
931 : }
932 86 : }
933 :
934 :
935 :
936 : template <typename file_id_type>
937 86 : file_id_type CheckpointIO::read_header (const std::string & name)
938 : {
939 86 : 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 110 : if (this->processor_id() == 0)
953 : {
954 154 : Xdr io (name, this->binary() ? DECODE : READ);
955 :
956 : // read the version
957 56 : io.data(input_version);
958 :
959 : // read the data type, don't care about it this time
960 : header_id_type data_size;
961 56 : io.data (data_size);
962 :
963 : // read the dimension
964 56 : io.data (mesh_dimension);
965 :
966 : // Read whether or not this is a parallel file
967 56 : io.data(input_parallel);
968 :
969 : // With how many processors?
970 56 : if (input_parallel)
971 32 : io.data(input_n_procs);
972 :
973 : // read subdomain names
974 56 : this->read_subdomain_names<file_id_type>(io);
975 :
976 : // read boundary names
977 14 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
978 :
979 56 : this->read_bc_names<file_id_type>(io, boundary_info, true); // sideset names
980 56 : 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 56 : const bool read_extra_integers = this->version_at_least_1_5();
985 14 : std::swap(input_version, this->version());
986 :
987 56 : if (read_extra_integers)
988 28 : this->read_integers_names<file_id_type>
989 28 : (io, node_integer_names, elem_integer_names);
990 28 : }
991 :
992 : // broadcast data from processor 0, set values everywhere
993 86 : this->comm().broadcast(input_version);
994 24 : this->version() = input_version;
995 :
996 86 : this->comm().broadcast(mesh_dimension);
997 148 : mesh.set_mesh_dimension(cast_int<unsigned char>(mesh_dimension));
998 :
999 86 : this->comm().broadcast(input_parallel);
1000 :
1001 86 : if (input_parallel)
1002 56 : 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 86 : this->comm().broadcast(subdomain_map);
1009 :
1010 24 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1011 86 : this->comm().broadcast(boundary_info.set_sideset_name_map());
1012 86 : this->comm().broadcast(boundary_info.set_nodeset_name_map());
1013 :
1014 86 : this->comm().broadcast(node_integer_names);
1015 86 : this->comm().broadcast(elem_integer_names);
1016 :
1017 162 : for (auto & int_name : node_integer_names)
1018 128 : mesh.add_node_integer(int_name);
1019 :
1020 124 : for (auto & int_name : elem_integer_names)
1021 64 : mesh.add_elem_integer(int_name);
1022 :
1023 172 : return input_parallel ? input_n_procs : 0;
1024 38 : }
1025 :
1026 :
1027 :
1028 : template <typename file_id_type>
1029 82 : void CheckpointIO::read_subfile (Xdr & io, bool expect_all_remote)
1030 : {
1031 : // read the nodal locations
1032 82 : this->read_nodes<file_id_type> (io);
1033 :
1034 : // read connectivity
1035 82 : this->read_connectivity<file_id_type> (io);
1036 :
1037 : // read remote_elem connectivity
1038 82 : this->read_remote_elem<file_id_type> (io, expect_all_remote);
1039 :
1040 : // read the boundary conditions
1041 82 : this->read_bcs<file_id_type> (io);
1042 :
1043 : // read the nodesets
1044 82 : this->read_nodesets<file_id_type> (io);
1045 82 : }
1046 :
1047 :
1048 :
1049 : template <typename file_id_type>
1050 56 : void CheckpointIO::read_subdomain_names(Xdr & io)
1051 : {
1052 56 : 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 56 : subdomain_ids.reserve(subdomain_map.size());
1059 :
1060 42 : std::vector<std::string> subdomain_names;
1061 56 : subdomain_names.reserve(subdomain_map.size());
1062 :
1063 56 : file_id_type n_subdomain_names = 0;
1064 56 : io.data(n_subdomain_names, "# subdomain id to name map");
1065 :
1066 56 : 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 84 : }
1076 :
1077 :
1078 :
1079 : template <typename file_id_type>
1080 82 : void CheckpointIO::read_nodes (Xdr & io)
1081 : {
1082 : // convenient reference to our mesh
1083 82 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1084 :
1085 : file_id_type n_nodes_here;
1086 82 : io.data(n_nodes_here, "# n_nodes on proc");
1087 :
1088 82 : const bool read_extra_integers = this->version_at_least_1_5();
1089 :
1090 44 : const unsigned int n_extra_integers =
1091 60 : read_extra_integers ? mesh.n_node_integers() : 0;
1092 :
1093 : // Will hold the node id and pid and extra integers
1094 104 : std::vector<file_id_type> id_pid(2 + n_extra_integers);
1095 :
1096 : // For the coordinates
1097 104 : std::vector<Real> coords(LIBMESH_DIM);
1098 :
1099 1568 : for (unsigned int i=0; i<n_nodes_here; i++)
1100 : {
1101 1486 : io.data_stream(id_pid.data(), 2 + n_extra_integers, 2 + n_extra_integers);
1102 :
1103 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1104 1486 : file_id_type unique_id = 0;
1105 1486 : io.data(unique_id, "# unique id");
1106 : #endif
1107 :
1108 1486 : io.data_stream(coords.data(), LIBMESH_DIM, LIBMESH_DIM);
1109 :
1110 394 : Point p;
1111 1486 : p(0) = coords[0];
1112 :
1113 : #if LIBMESH_DIM > 1
1114 1486 : p(1) = coords[1];
1115 : #endif
1116 :
1117 : #if LIBMESH_DIM > 2
1118 1486 : p(2) = coords[2];
1119 : #endif
1120 :
1121 1486 : 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 1880 : 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 1486 : const Node * old_node = mesh.query_node_ptr(id);
1133 :
1134 1486 : 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 1276 : mesh.add_point(p, id, pid);
1155 :
1156 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1157 1276 : node->set_unique_id(unique_id);
1158 : #endif
1159 :
1160 334 : libmesh_assert_equal_to(n_extra_integers, node->n_extra_integers());
1161 :
1162 2112 : for (unsigned int ei=0; ei != n_extra_integers; ++ei)
1163 : {
1164 1064 : const dof_id_type extra_int = cast_int<dof_id_type>(id_pid[2+ei]);
1165 836 : node->set_extra_integer(ei, extra_int);
1166 : }
1167 : }
1168 : }
1169 82 : }
1170 :
1171 :
1172 :
1173 : template <typename file_id_type>
1174 82 : void CheckpointIO::read_connectivity (Xdr & io)
1175 : {
1176 : // convenient reference to our mesh
1177 82 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1178 :
1179 82 : const bool read_extra_integers = this->version_at_least_1_5();
1180 82 : const bool read_runtime_topology = this->version_at_least_1_6();
1181 :
1182 44 : const unsigned int n_extra_integers =
1183 60 : read_extra_integers ? mesh.n_elem_integers() : 0;
1184 :
1185 : file_id_type n_elems_here;
1186 82 : io.data(n_elems_here);
1187 :
1188 : // Keep track of the highest dimensional element we've added to the mesh
1189 82 : 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 1240 : for (unsigned int i=0; i<n_elems_here; i++)
1197 : {
1198 : // id type pid subdomain_id parent_id
1199 997 : std::vector<file_id_type> elem_data(6 + n_extra_integers);
1200 627 : io.data_stream
1201 370 : (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 788 : file_id_type unique_id = 0;
1206 788 : io.data(unique_id, "# unique id");
1207 : #endif
1208 :
1209 : #ifdef LIBMESH_ENABLE_AMR
1210 788 : uint16_t p_level = 0;
1211 788 : io.data(p_level, "# p_level");
1212 :
1213 : uint16_t rflag, pflag;
1214 788 : io.data(rflag, "# rflag");
1215 788 : io.data(pflag, "# pflag");
1216 : #endif
1217 :
1218 788 : const ElemType elem_type =
1219 788 : 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 788 : 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 788 : if (has_runtime_topology)
1229 : {
1230 16 : 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 16 : io.data(runtime_topology, "# runtime topology");
1237 20 : libmesh_error_msg_if(runtime_topology.size() < 2,
1238 : "Invalid runtime element topology.");
1239 :
1240 4 : std::size_t topology_index = 0;
1241 12 : n_nodes =
1242 16 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1243 : const unsigned int n_sides =
1244 16 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1245 16 : nodes_on_sides.resize(n_sides);
1246 :
1247 120 : for (auto s : index_range(nodes_on_sides))
1248 : {
1249 130 : 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 104 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1255 104 : 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 104 : side_nodes.resize(n_side_nodes);
1261 472 : for (auto n : index_range(side_nodes))
1262 : {
1263 : const unsigned int local_node =
1264 460 : cast_int<unsigned int>(runtime_topology[topology_index++]);
1265 368 : libmesh_error_msg_if
1266 : (local_node >= n_nodes,
1267 : "Runtime element side checkpoint topology references "
1268 : "an invalid local node.");
1269 368 : side_nodes[n] = local_node;
1270 : }
1271 : }
1272 :
1273 20 : libmesh_error_msg_if
1274 : (topology_index != runtime_topology.size(),
1275 : "Extra data in runtime element checkpoint topology.");
1276 :
1277 16 : if (is_c0polygon)
1278 : {
1279 8 : libmesh_error_msg_if
1280 : (n_nodes < 3 || n_sides != n_nodes,
1281 : "Invalid C0POLYGON checkpoint topology.");
1282 48 : for (const auto & side_nodes : nodes_on_sides)
1283 50 : libmesh_error_msg_if
1284 : (side_nodes.size() != 2,
1285 : "Invalid C0POLYGON side checkpoint topology.");
1286 : }
1287 8 : else if (is_c0polyhedron)
1288 : {
1289 8 : libmesh_error_msg_if(n_sides < 4,
1290 : "Invalid C0POLYHEDRON checkpoint topology.");
1291 72 : for (const auto & side_nodes : nodes_on_sides)
1292 80 : 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 997 : std::vector<file_id_type> conn_data(n_nodes);
1300 627 : io.data_stream
1301 370 : (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 788 : cast_int<dof_id_type> (elem_data[0]);
1306 209 : const processor_id_type proc_id =
1307 : cast_int<processor_id_type>
1308 997 : (elem_data[2] % mesh.n_processors());
1309 209 : const subdomain_id_type subdomain_id =
1310 788 : 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 788 : 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 579 : Elem * parent =
1321 209 : (elem_data[4] == static_cast<largest_id_type>(-1) ||
1322 788 : (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 788 : (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 788 : 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 788 : 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 499 : std::unique_ptr<Elem> elem;
1386 676 : std::unique_ptr<Node> generated_mid_elem_node;
1387 :
1388 676 : if (is_c0polygon)
1389 10 : elem = std::make_unique<C0Polygon>(n_nodes, parent);
1390 668 : else if (is_c0polyhedron)
1391 : {
1392 14 : std::vector<std::shared_ptr<Polygon>> sides(nodes_on_sides.size());
1393 72 : 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 80 : (cast_int<unsigned int>(side_node_indices.size()));
1399 352 : for (auto n : index_range(side_node_indices))
1400 288 : side->set_node
1401 144 : (n, mesh.node_ptr(cast_int<dof_id_type>
1402 360 : (conn_data[side_node_indices[n]])));
1403 32 : sides[s] = std::move(side);
1404 : }
1405 :
1406 8 : elem = std::make_unique<C0Polyhedron>
1407 : (sides, generated_mid_elem_node, parent);
1408 :
1409 8 : libmesh_error_msg_if
1410 : (elem->n_nodes() != conn_data.size(),
1411 : "C0POLYHEDRON checkpoint topology is incompatible with "
1412 : "this libMesh configuration.");
1413 :
1414 104 : for (auto n : make_range(elem->n_vertices()))
1415 168 : 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 4 : }
1421 : else
1422 1147 : elem = Elem::build(elem_type, parent);
1423 :
1424 676 : if (has_runtime_topology)
1425 : {
1426 16 : 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 120 : for (auto s : index_range(nodes_on_sides))
1435 208 : 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 676 : elem->set_unique_id(unique_id);
1443 : #endif
1444 :
1445 676 : if (elem->dim() > highest_elem_dim)
1446 0 : highest_elem_dim = elem->dim();
1447 :
1448 676 : elem->set_id() = id;
1449 676 : elem->processor_id() = proc_id;
1450 676 : elem->subdomain_id() = subdomain_id;
1451 :
1452 : #ifdef LIBMESH_ENABLE_AMR
1453 676 : elem->hack_p_level(p_level);
1454 :
1455 676 : elem->set_refinement_flag (cast_int<Elem::RefinementState>(rflag));
1456 676 : elem->set_p_refinement_flag(cast_int<Elem::RefinementState>(pflag));
1457 :
1458 : // Set parent connections
1459 676 : 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 2892 : for (unsigned int n=0,
1473 354 : n_conn = cast_int<unsigned int>(conn_data.size());
1474 3568 : n != n_conn; n++)
1475 4402 : elem->set_node(n,
1476 2892 : mesh.node_ptr(cast_int<dof_id_type>(conn_data[n])));
1477 :
1478 1030 : 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 819 : for (unsigned int ei=0; ei != n_extra_integers; ++ei)
1482 : {
1483 182 : const dof_id_type extra_int = cast_int<dof_id_type>(elem_data[6+ei]);
1484 143 : added_elem->set_extra_integer(ei, extra_int);
1485 : }
1486 322 : }
1487 : }
1488 :
1489 82 : mesh.set_mesh_dimension(cast_int<unsigned char>(highest_elem_dim));
1490 82 : }
1491 :
1492 :
1493 : template <typename file_id_type>
1494 82 : void CheckpointIO::read_remote_elem (Xdr & io, bool libmesh_dbg_var(expect_all_remote))
1495 : {
1496 : // convenient reference to our mesh
1497 82 : 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 82 : io.data(elem_ids, "# remote neighbor elem_ids");
1504 82 : io.data(elem_sides, "# remote neighbor elem_sides");
1505 :
1506 22 : libmesh_assert_equal_to(elem_ids.size(), elem_sides.size());
1507 :
1508 290 : for (auto i : index_range(elem_ids))
1509 : {
1510 336 : Elem & elem = mesh.elem_ref(cast_int<dof_id_type>(elem_ids[i]));
1511 336 : if (!elem.neighbor_ptr(elem_sides[i]))
1512 208 : 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 82 : io.data(parent_ids, "# remote child parent_ids");
1523 82 : io.data(child_numbers, "# remote child_numbers");
1524 :
1525 : #ifdef LIBMESH_ENABLE_AMR
1526 82 : 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 82 : }
1544 :
1545 :
1546 :
1547 : template <typename file_id_type>
1548 82 : void CheckpointIO::read_bcs (Xdr & io)
1549 : {
1550 : // convenient reference to our mesh
1551 82 : 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 82 : io.data(element_id_list, "# element ids for bcs");
1561 82 : io.data(side_list, "# sides of elements for bcs");
1562 82 : io.data(bc_id_list, "# bc ids");
1563 :
1564 802 : for (auto i : index_range(element_id_list))
1565 558 : boundary_info.add_side
1566 906 : (cast_int<dof_id_type>(element_id_list[i]), side_list[i],
1567 906 : cast_int<boundary_id_type>(bc_id_list[i]));
1568 82 : }
1569 :
1570 :
1571 :
1572 : template <typename file_id_type>
1573 82 : void CheckpointIO::read_nodesets (Xdr & io)
1574 : {
1575 : // convenient reference to our mesh
1576 82 : 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 82 : io.data(node_id_list, "# node id list");
1585 82 : io.data(bc_id_list, "# nodeset bc id list");
1586 :
1587 1102 : for (auto i : index_range(node_id_list))
1588 792 : boundary_info.add_node
1589 1020 : (cast_int<dof_id_type>(node_id_list[i]),
1590 1284 : cast_int<boundary_id_type>(bc_id_list[i]));
1591 82 : }
1592 :
1593 :
1594 :
1595 : template <typename file_id_type>
1596 112 : void CheckpointIO::read_bc_names(Xdr & io, BoundaryInfo & info, bool is_sideset)
1597 : {
1598 112 : 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 112 : file_id_type n_boundary_names = 0;
1605 :
1606 112 : if (is_sideset)
1607 56 : io.data(n_boundary_names, "# sideset id to name map");
1608 : else
1609 56 : io.data(n_boundary_names, "# nodeset id to name map");
1610 :
1611 112 : if (n_boundary_names)
1612 : {
1613 80 : io.data(boundary_ids);
1614 80 : io.data(boundary_names);
1615 : }
1616 :
1617 : // Add them back into the map
1618 432 : for (auto i : index_range(boundary_ids))
1619 720 : boundary_map[cast_int<boundary_id_type>(boundary_ids[i])] =
1620 160 : boundary_names[i];
1621 168 : }
1622 :
1623 :
1624 : template <typename file_id_type>
1625 56 : 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 56 : io.data(n_node_integers, "# n_extra_integers per node");
1633 56 : io.data(node_integer_names);
1634 56 : io.data(n_elem_integers, "# n_extra_integers per elem");
1635 56 : io.data(elem_integer_names);
1636 56 : }
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
|
