checkpoint_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 #include "libmesh/checkpoint_io.h"
19 
20 // C++ includes
21 #include <iostream>
22 #include <iomanip>
23 #include <cstdio>
24 #include <vector>
25 #include <string>
26 #include <fstream>
27 #include <sstream> // for ostringstream
28 
29 // Local includes
30 #include "libmesh/xdr_io.h"
31 #include "libmesh/legacy_xdr_io.h"
32 #include "libmesh/xdr_cxx.h"
33 #include "libmesh/enum_xdr_mode.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/node.h"
36 #include "libmesh/elem.h"
37 #include "libmesh/boundary_info.h"
38 #include "libmesh/parallel.h"
39 #include "libmesh/mesh_tools.h"
40 #include "libmesh/partitioner.h"
43 #include "libmesh/parallel_mesh.h"
44 
45 namespace libMesh
46 {
47 
48 // ------------------------------------------------------------
49 // CheckpointIO members
50 CheckpointIO::CheckpointIO (MeshBase& mesh, const bool binary_in) :
51  MeshInput<MeshBase> (mesh,/* is_parallel_format = */ true),
52  MeshOutput<MeshBase>(mesh,/* is_parallel_format = */ true),
53  ParallelObject (mesh),
54  _binary (binary_in),
55  _version ("checkpoint-1.0")
56 {
57 }
58 
59 
60 
61 CheckpointIO::CheckpointIO (const MeshBase& mesh, const bool binary_in) :
62  MeshOutput<MeshBase>(mesh,/* is_parallel_format = */ true),
63  ParallelObject (mesh),
64  _binary (binary_in)
65 {
66 }
67 
68 
69 
71 {
72 }
73 
74 
75 
76 
77 void CheckpointIO::write (const std::string& name)
78 {
79  START_LOG("write()","CheckpointIO");
80 
81  // convenient reference to our mesh
83 
84  // Try to dynamic cast the mesh to see if it's a ParallelMesh object
85  // Note: Just using is_serial() is not good enough because the Mesh won't
86  // have been prepared yet when is when that flag gets set to false... sigh.
87  bool parallel_mesh = dynamic_cast<const ParallelMesh*>(&mesh);
88 
89  // If this is a serial mesh then we're only going to write it on processor 0
90  if(parallel_mesh || this->processor_id() == 0)
91  {
92  std::ostringstream file_name_stream;
93 
94  file_name_stream << name;
95 
96  if(parallel_mesh)
97  file_name_stream << "-" << this->processor_id();
98 
99  Xdr io (file_name_stream.str(), this->binary() ? ENCODE : WRITE);
100 
101  // write the version
102  io.data(_version, "# version");
103 
104  // Write out whether or not this is a serial mesh (helps with error checking on read)
105  {
106  unsigned int parallel = parallel_mesh;
107  io.data(parallel, "# parallel");
108  }
109 
110  // If we're writing out a parallel mesh then we need to write the number of processors
111  // so we can check it upon reading the file
112  if(parallel_mesh)
113  {
114  largest_id_type n_procs = this->n_processors();
115  io.data(n_procs, "# n_procs");
116  }
117 
118  // write subdomain names
119  this->write_subdomain_names(io);
120 
121  // write the nodal locations
122  this->write_nodes (io);
123 
124  // write connectivity
125  this->write_connectivity (io);
126 
127  // write the boundary condition information
128  this->write_bcs (io);
129 
130  // write the nodeset information
131  this->write_nodesets (io);
132 
133  // pause all processes until the writing ends -- this will
134  // protect for the pathological case where a write is
135  // followed immediately by a read. The write must be
136  // guaranteed to complete first.
137  io.close();
138  }
139 
140  this->comm().barrier();
141 
142  STOP_LOG("write()","CheckpointIO");
143 }
144 
145 
146 
148 {
149  {
151 
152  const std::map<subdomain_id_type, std::string> & subdomain_map = mesh.get_subdomain_name_map();
153 
154  std::vector<header_id_type> subdomain_ids; subdomain_ids.reserve(subdomain_map.size());
155  std::vector<std::string> subdomain_names; subdomain_names.reserve(subdomain_map.size());
156 
157  // We need to loop over the map and make sure that there aren't any invalid entries. Since we
158  // return writable references in mesh_base, it's possible for the user to leave some entity names
159  // blank. We can't write those to the XDA file.
160  header_id_type n_subdomain_names = 0;
161  std::map<subdomain_id_type, std::string>::const_iterator it_end = subdomain_map.end();
162  for (std::map<subdomain_id_type, std::string>::const_iterator it = subdomain_map.begin(); it != it_end; ++it)
163  {
164  if (!it->second.empty())
165  {
166  n_subdomain_names++;
167  subdomain_ids.push_back(it->first);
168  subdomain_names.push_back(it->second);
169  }
170  }
171 
172  io.data(n_subdomain_names, "# subdomain id to name map");
173  // Write out the ids and names in two vectors
174  if (n_subdomain_names)
175  {
176  io.data(subdomain_ids);
177  io.data(subdomain_names);
178  }
179  }
180 }
181 
182 
183 
185 {
186  // convenient reference to our mesh
188 
190  it = mesh.nodes_begin(),
191  end = mesh.nodes_end();
192 
193  unsigned int n_nodes_here = MeshTools::n_nodes(it, end);
194 
195  io.data(n_nodes_here, "# n_nodes on proc");
196 
197  it = mesh.nodes_begin();
198 
199  // Will hold the node id and pid
200  std::vector<largest_id_type> id_pid(2);
201 
202  // For the coordinates
203  std::vector<Real> coords(LIBMESH_DIM);
204 
205  for(; it != end; ++it)
206  {
207  Node & node = *(*it);
208 
209  id_pid[0] = node.id();
210  id_pid[1] = node.processor_id();
211 
212  io.data_stream(&id_pid[0], 2, 2);
213 
214 #ifdef LIBMESH_ENABLE_UNIQUE_ID
215  largest_id_type unique_id = node.unique_id();
216 
217  io.data(unique_id, "# unique id");
218 #endif
219 
220  coords[0] = node(0);
221 
222 #if LIBMESH_DIM > 1
223  coords[1] = node(1);
224 #endif
225 
226 #if LIBMESH_DIM > 2
227  coords[2] = node(2);
228 #endif
229 
230  io.data_stream(&coords[0], LIBMESH_DIM, 3);
231  }
232 }
233 
234 
235 
237 {
238  libmesh_assert (io.writing());
239 
240  // convenient reference to our mesh
242 
243  // We will only write active elements and their parents.
244  unsigned int n_active_levels = n_active_levels_on_processor(mesh);
245 
246  std::vector<xdr_id_type> n_elem_at_level(n_active_levels);
247 
248  // Find the number of elements at each level
249  for (unsigned int level=0; level<n_active_levels; level++)
250  {
251  MeshBase::const_element_iterator it = mesh.level_elements_begin(level);
252  MeshBase::const_element_iterator end = mesh.level_elements_end(level);
253 
254  n_elem_at_level[level] = MeshTools::n_elem(it, end);
255  }
256 
257  typedef std::map<dof_id_type, std::pair<processor_id_type, dof_id_type> > id_map_type;
258  id_map_type parent_id_map, child_id_map;
259 
260  io.data(n_active_levels, "# n_active_levels");
261 
262  for(unsigned int level=0; level < n_active_levels; level++)
263  {
264  std::ostringstream comment;
265  comment << "# n_elem at level ";
266  comment << level ;
267  io.data (n_elem_at_level[level], comment.str().c_str());
268 
269  MeshBase::const_element_iterator it = mesh.level_elements_begin(level);
270  MeshBase::const_element_iterator end = mesh.level_elements_end(level);
271  for (; it != end; ++it)
272  {
273  Elem & elem = *(*it);
274 
275  unsigned int n_nodes = elem.n_nodes();
276 
277  // id type pid subdomain_id parent_id
278  std::vector<largest_id_type> elem_data(5);
279 
280  elem_data[0] = elem.id();
281  elem_data[1] = elem.type();
282  elem_data[2] = elem.processor_id();
283  elem_data[3] = elem.subdomain_id();
284 
285  if(elem.parent() != NULL)
286  elem_data[4] = elem.parent()->id();
287  else
288  elem_data[4] = DofObject::invalid_processor_id;
289 
290  std::vector<largest_id_type> conn_data(n_nodes);
291 
292  for(unsigned int i=0; i<n_nodes; i++)
293  conn_data[i] = elem.node(i);
294 
295  io.data_stream(&elem_data[0], elem_data.size(), elem_data.size());
296 
297 #ifdef LIBMESH_ENABLE_UNIQUE_ID
298  largest_id_type unique_id = elem.unique_id();
299 
300  io.data(unique_id, "# unique id");
301 #endif
302 
303 #ifdef LIBMESH_ENABLE_AMR
304  unsigned int p_level = elem.p_level();
305 
306  io.data(p_level, "# p_level");
307 #endif
308 
309  io.data_stream(&conn_data[0], conn_data.size(), conn_data.size());
310  }
311  }
312 }
313 
314 
315 
316 void CheckpointIO::write_bcs (Xdr &io) const
317 {
318  libmesh_assert (io.writing());
319 
320  // convenient reference to our mesh
322 
323  // and our boundary info object
324  const BoundaryInfo &boundary_info = *mesh.boundary_info;
325 
326  // Version 0.9.2+ introduces entity names
327  write_bc_names(io, boundary_info, true); // sideset names
328  write_bc_names(io, boundary_info, false); // nodeset names
329 
330  std::vector<dof_id_type> element_id_list;
331  std::vector<unsigned short int> side_list;
332  std::vector<boundary_id_type> bc_id_list;
333 
334  boundary_info.build_side_list(element_id_list, side_list, bc_id_list);
335 
336  io.data(element_id_list, "# element ids for bcs");
337  io.data(side_list, "# sides of elements for bcs");
338  io.data(bc_id_list, "# bc ids");
339 }
340 
341 
342 
344 {
345  libmesh_assert (io.writing());
346 
347  // convenient reference to our mesh
349 
350  // and our boundary info object
351  const BoundaryInfo &boundary_info = *mesh.boundary_info;
352 
353  std::vector<dof_id_type> node_id_list;
354  std::vector<boundary_id_type> bc_id_list;
355 
356  boundary_info.build_node_list(node_id_list, bc_id_list);
357 
358  io.data(node_id_list, "# node id list");
359  io.data(bc_id_list, "# nodeset bc id list");
360 }
361 
362 
363 
364 void CheckpointIO::write_bc_names (Xdr &io, const BoundaryInfo & info, bool is_sideset) const
365 {
366  const std::map<boundary_id_type, std::string> & boundary_map = is_sideset ?
368 
369  std::vector<boundary_id_type> boundary_ids; boundary_ids.reserve(boundary_map.size());
370  std::vector<std::string> boundary_names; boundary_names.reserve(boundary_map.size());
371 
372  // We need to loop over the map and make sure that there aren't any invalid entries. Since we
373  // return writable references in boundary_info, it's possible for the user to leave some entity names
374  // blank. We can't write those to the XDA file.
375  header_id_type n_boundary_names = 0;
376  std::map<boundary_id_type, std::string>::const_iterator it_end = boundary_map.end();
377  for (std::map<boundary_id_type, std::string>::const_iterator it = boundary_map.begin(); it != it_end; ++it)
378  {
379  if (!it->second.empty())
380  {
381  n_boundary_names++;
382  boundary_ids.push_back(it->first);
383  boundary_names.push_back(it->second);
384  }
385  }
386 
387  if (is_sideset)
388  io.data(n_boundary_names, "# sideset id to name map");
389  else
390  io.data(n_boundary_names, "# nodeset id to name map");
391  // Write out the ids and names in two vectors
392  if (n_boundary_names)
393  {
394  io.data(boundary_ids);
395  io.data(boundary_names);
396  }
397 }
398 
399 
400 
401 void CheckpointIO::read (const std::string& name)
402 {
403  START_LOG("read()","CheckpointIO");
404 
406 
407  // Try to dynamic cast the mesh to see if it's a ParallelMesh object
408  // Note: Just using is_serial() is not good enough because the Mesh won't
409  // have been prepared yet when is when that flag gets set to false... sigh.
410  bool parallel_mesh = dynamic_cast<ParallelMesh*>(&mesh);
411 
412  // If this is a serial mesh then we're going to only read it on processor 0 and broadcast it
413  if(parallel_mesh || this->processor_id() == 0)
414  {
415  std::ostringstream file_name_stream;
416 
417  file_name_stream << name;
418 
419  if(parallel_mesh)
420  file_name_stream << "-" << this->processor_id();
421 
422  {
423  std::ifstream in (file_name_stream.str().c_str());
424 
425  if (!in.good())
426  {
427  libMesh::err << "ERROR: cannot locate specified file:\n\t"
428  << file_name_stream.str()
429  << std::endl;
430  libmesh_error();
431  }
432  }
433 
434  Xdr io (file_name_stream.str(), this->binary() ? DECODE : READ);
435 
436  // read the version
437  io.data (_version);
438 
439  // Check if the mesh we're reading is the same as the one that was written
440  {
441  unsigned int parallel;
442  io.data(parallel, "# parallel");
443 
444  if(parallel_mesh != parallel)
445  {
446  libMesh::err << "Attempted to utilize a checkpoint file with an incompatible mesh distribution!" << std::endl;
447  libmesh_error();
448  }
449  }
450 
451  // If this is a parallel mesh then we need to check to ensure we're reading this on the same number of procs
452  if(parallel_mesh)
453  {
454  largest_id_type n_procs;
455  io.data(n_procs, "# n_procs");
456 
457  if(n_procs != this->n_processors())
458  {
459  libMesh::err << "Attempted to utilize a checkpoint file on " << this->n_processors() << " processors but it was written using " << n_procs << "!!";
460  libmesh_error();
461  }
462  }
463 
464  // read subdomain names
465  this->read_subdomain_names(io);
466 
467  // read the nodal locations
468  this->read_nodes (io);
469 
470  // read connectivity
471  this->read_connectivity (io);
472 
473  // read the boundary conditions
474  this->read_bcs (io);
475 
476  // read the nodesets
477  this->read_nodesets (io);
478 
479  io.close();
480  }
481 
482  // If the mesh is serial then we only read it on processor 0 so we need to broadcast it
483  if(!parallel_mesh)
485 
486  STOP_LOG("read()","CheckpointIO");
487 }
488 
489 
490 
492 {
494 
495  std::map<subdomain_id_type, std::string> & subdomain_map = mesh.set_subdomain_name_map();
496 
497  std::vector<header_id_type> subdomain_ids; subdomain_ids.reserve(subdomain_map.size());
498  std::vector<std::string> subdomain_names; subdomain_names.reserve(subdomain_map.size());
499 
500  header_id_type n_subdomain_names = 0;
501  io.data(n_subdomain_names, "# subdomain id to name map");
502 
503  if(n_subdomain_names)
504  {
505  io.data(subdomain_ids);
506  io.data(subdomain_names);
507 
508  for(unsigned int i=0; i<subdomain_ids.size(); i++)
509  subdomain_map[subdomain_ids[i]] = subdomain_names[i];
510  }
511 }
512 
513 
514 
516 {
517  // convenient reference to our mesh
519 
520  unsigned int n_nodes_here;
521  io.data(n_nodes_here, "# n_nodes on proc");
522 
523  // Will hold the node id and pid
524  std::vector<largest_id_type> id_pid(2);
525 
526  // For the coordinates
527  std::vector<Real> coords(LIBMESH_DIM);
528 
529  for(unsigned int i=0; i<n_nodes_here; i++)
530  {
531  io.data_stream(&id_pid[0], 2, 2);
532 
533 #ifdef LIBMESH_ENABLE_UNIQUE_ID
534  largest_id_type unique_id = 0;
535  io.data(unique_id, "# unique id");
536 #endif
537 
538  io.data_stream(&coords[0], LIBMESH_DIM, LIBMESH_DIM);
539 
540  Point p;
541  p(0) = coords[0];
542 
543 #if LIBMESH_DIM > 1
544  p(1) = coords[1];
545 #endif
546 
547 #if LIBMESH_DIM > 2
548  p(2) = coords[2];
549 #endif
550 
551 #ifdef LIBMESH_ENABLE_UNIQUE_ID
552  Node * node =
553 #endif
554  mesh.add_point(p, id_pid[0], id_pid[1]);
555 
556 #ifdef LIBMESH_ENABLE_UNIQUE_ID
557  node->set_unique_id() = unique_id;
558 #endif
559  }
560 }
561 
562 
563 
565 {
566  // convenient reference to our mesh
568 
569  unsigned int n_active_levels;
570  io.data(n_active_levels, "# n_active_levels");
571 
572  // Keep track of the highest dimensional element we've added to the mesh
573  unsigned int highest_elem_dim = 1;
574 
575  for(unsigned int level=0; level < n_active_levels; level++)
576  {
577  xdr_id_type n_elem_at_level = 0;
578  io.data (n_elem_at_level, "");
579 
580  for (unsigned int i=0; i<n_elem_at_level; i++)
581  {
582  // id type pid subdomain_id parent_id
583  std::vector<largest_id_type> elem_data(5);
584  io.data_stream(&elem_data[0], elem_data.size(), elem_data.size());
585 
586 #ifdef LIBMESH_ENABLE_UNIQUE_ID
587  largest_id_type unique_id = 0;
588  io.data(unique_id, "# unique id");
589 #endif
590 
591 #ifdef LIBMESH_ENABLE_AMR
592  unsigned int p_level = 0;
593 
594  io.data(p_level, "# p_level");
595 #endif
596 
597  unsigned int n_nodes = Elem::type_to_n_nodes_map[elem_data[1]];
598 
599  // Snag the node ids this element was connected to
600  std::vector<largest_id_type> conn_data(n_nodes);
601  io.data_stream(&conn_data[0], conn_data.size(), conn_data.size());
602 
603  const dof_id_type id = elem_data[0];
604  const ElemType elem_type = static_cast<ElemType>(elem_data[1]);
605  const processor_id_type processor_id = elem_data[2];
606  const subdomain_id_type subdomain_id = elem_data[3];
607  const dof_id_type parent_id = elem_data[4];
608 
609  Elem * parent = (parent_id == DofObject::invalid_processor_id) ? NULL : mesh.elem(parent_id);
610 
611  // Create the element
612  Elem * elem = Elem::build(elem_type, parent).release();
613 
614 #ifdef LIBMESH_ENABLE_UNIQUE_ID
615  elem->set_unique_id() = unique_id;
616 #endif
617 
618  if(elem->dim() > highest_elem_dim)
619  highest_elem_dim = elem->dim();
620 
621  elem->set_id() = id;
622  elem->processor_id() = processor_id;
623  elem->subdomain_id() = subdomain_id;
624 
625 #ifdef LIBMESH_ENABLE_AMR
626  elem->hack_p_level(p_level);
627 
628  // Set parent connections
629  if(parent)
630  {
631  parent->add_child(elem);
634  }
635 #endif
636 
637  libmesh_assert(elem->n_nodes() == conn_data.size());
638 
639  // Connect all the nodes to this element
640  for (unsigned int n=0; n<conn_data.size(); n++)
641  elem->set_node(n) = mesh.node_ptr(conn_data[n]);
642 
643  mesh.add_elem(elem);
644  }
645  }
646 
647  mesh.set_mesh_dimension(highest_elem_dim);
648 }
649 
650 
651 
653 {
654  // convenient reference to our mesh
656 
657  // and our boundary info object
658  BoundaryInfo &boundary_info = *mesh.boundary_info;
659 
660  // Version 0.9.2+ introduces entity names
661  read_bc_names(io, boundary_info, true); // sideset names
662  read_bc_names(io, boundary_info, false); // nodeset names
663 
664  std::vector<dof_id_type> element_id_list;
665  std::vector<unsigned short int> side_list;
666  std::vector<boundary_id_type> bc_id_list;
667 
668  io.data(element_id_list, "# element ids for bcs");
669  io.data(side_list, "# sides of elements for bcs");
670  io.data(bc_id_list, "# bc ids");
671 
672  for(unsigned int i=0; i<element_id_list.size(); i++)
673  boundary_info.add_side(element_id_list[i], side_list[i], bc_id_list[i]);
674 }
675 
676 
677 
679 {
680  // convenient reference to our mesh
682 
683  // and our boundary info object
684  BoundaryInfo &boundary_info = *mesh.boundary_info;
685 
686  std::vector<dof_id_type> node_id_list;
687  std::vector<boundary_id_type> bc_id_list;
688 
689  io.data(node_id_list, "# node id list");
690  io.data(bc_id_list, "# nodeset bc id list");
691 
692  for(unsigned int i=0; i<node_id_list.size(); i++)
693  boundary_info.add_node(node_id_list[i], bc_id_list[i]);
694 }
695 
696 
697 
698 void CheckpointIO::read_bc_names(Xdr &io, BoundaryInfo & info, bool is_sideset)
699 {
700  std::map<boundary_id_type, std::string> & boundary_map = is_sideset ?
702 
703  std::vector<boundary_id_type> boundary_ids;
704  std::vector<std::string> boundary_names;
705 
706  header_id_type n_boundary_names = 0;
707 
708  if (is_sideset)
709  io.data(n_boundary_names, "# sideset id to name map");
710  else
711  io.data(n_boundary_names, "# nodeset id to name map");
712 
713  if (n_boundary_names)
714  {
715  io.data(boundary_ids);
716  io.data(boundary_names);
717  }
718 
719  // Add them back into the map
720  for(unsigned int i=0; i<boundary_ids.size(); i++)
721  boundary_map[boundary_ids[i]] = boundary_names[i];
722 }
723 
724 
726 {
727  unsigned int max_level = 0;
728 
731 
732  for( ; el != end_el; ++el)
733  max_level = std::max((*el)->level(), max_level);
734 
735  return max_level + 1;
736 }
737 
738 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo