libMesh
system_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 #include "libmesh/libmesh_common.h"
20 #include "libmesh/parallel.h"
21 
22 // C++ Includes
23 #include <cstdio> // for std::sprintf
24 #include <set>
25 #include <numeric> // for std::partial_sum
26 
27 // Local Include
28 #include "libmesh/libmesh_version.h"
29 #include "libmesh/system.h"
30 #include "libmesh/mesh_base.h"
31 //#include "libmesh/mesh_tools.h"
32 #include "libmesh/elem.h"
33 #include "libmesh/xdr_cxx.h"
34 #include "libmesh/numeric_vector.h"
35 #include "libmesh/dof_map.h"
36 
37 
38 
39 // Anonymous namespace for implementation details.
40 namespace {
41 
42 using libMesh::DofObject;
43 using libMesh::Number;
44 using libMesh::cast_int;
45 
46 // Comments:
47 // ---------
48 // - The max_io_blksize governs how many nodes or elements will be
49 // treated as a single block when performing parallel IO on large
50 // systems.
51 // - This parameter only loosely affects the size of the actual IO
52 // buffer as this depends on the number of components a given
53 // variable has for the nodes/elements in the block.
54 // - When reading/writing each processor uses an ID map which is
55 // 3*io_blksize*sizeof(dof_id_type) bytes long, so with unsigned int
56 // and // io_blksize=256000 we would expect that buffer alone to be
57 // ~3Mb.
58 // - In general, an increase in max_io_blksize should increase the
59 // efficiency of large parallel read/writes by reducing the number
60 // of MPI messages at the expense of memory.
61 // - If the library exhausts memory during IO you might reduce this
62 // parameter.
63 
64 const std::size_t max_io_blksize = 256000;
65 
66 
72 struct CompareDofObjectsByID
73 {
74  bool operator()(const DofObject * a,
75  const DofObject * b) const
76  {
77  libmesh_assert (a);
78  libmesh_assert (b);
79 
80  return a->id() < b->id();
81  }
82 };
83 
87 template <typename InValType>
88 class ThreadedIO
89 {
90 private:
91  libMesh::Xdr & _io;
92  std::vector<InValType> & _data;
93 
94 public:
95  ThreadedIO (libMesh::Xdr & io, std::vector<InValType> & data) :
96  _io(io),
97  _data(data)
98  {}
99 
100  void operator()()
101  {
102  if (_data.empty()) return;
103  _io.data_stream (&_data[0], cast_int<unsigned int>(_data.size()));
104  }
105 };
106 }
107 
108 
109 namespace libMesh
110 {
111 
112 
113 // ------------------------------------------------------------
114 // System class implementation
115 void System::read_header (Xdr & io,
116  const std::string & version,
117  const bool read_header_in,
118  const bool read_additional_data,
119  const bool read_legacy_format)
120 {
121  // This method implements the input of a
122  // System object, embedded in the output of
123  // an EquationSystems<T_sys>. This warrants some
124  // documentation. The output file essentially
125  // consists of 5 sections:
126  //
127  // for this system
128  //
129  // 5.) The number of variables in the system (unsigned int)
130  //
131  // for each variable in the system
132  //
133  // 6.) The name of the variable (string)
134  //
135  // 6.1.) Variable subdomains
136  //
137  // 7.) Combined in an FEType:
138  // - The approximation order(s) of the variable
139  // (Order Enum, cast to int/s)
140  // - The finite element family/ies of the variable
141  // (FEFamily Enum, cast to int/s)
142  //
143  // end variable loop
144  //
145  // 8.) The number of additional vectors (unsigned int),
146  //
147  // for each additional vector in the system object
148  //
149  // 9.) the name of the additional vector (string)
150  //
151  // end system
152  libmesh_assert (io.reading());
153 
154  // Possibly clear data structures and start from scratch.
155  if (read_header_in)
156  this->clear ();
157 
158  // Figure out if we need to read infinite element information.
159  // This will be true if the version string contains " with infinite elements"
160  const bool read_ifem_info =
161  (version.rfind(" with infinite elements") < version.size()) ||
162  libMesh::on_command_line ("--read_ifem_systems");
163 
164 
165  {
166  // 5.)
167  // Read the number of variables in the system
168  unsigned int nv=0;
169  if (this->processor_id() == 0)
170  io.data (nv);
171  this->comm().broadcast(nv);
172 
173  _written_var_indices.clear();
174  _written_var_indices.resize(nv, 0);
175 
176  for (unsigned int var=0; var<nv; var++)
177  {
178  // 6.)
179  // Read the name of the var-th variable
180  std::string var_name;
181  if (this->processor_id() == 0)
182  io.data (var_name);
183  this->comm().broadcast(var_name);
184 
185  // 6.1.)
186  std::set<subdomain_id_type> domains;
187  if (io.version() >= LIBMESH_VERSION_ID(0,7,2))
188  {
189  std::vector<subdomain_id_type> domain_array;
190  if (this->processor_id() == 0)
191  io.data (domain_array);
192  for (std::vector<subdomain_id_type>::iterator it = domain_array.begin(); it != domain_array.end(); ++it)
193  domains.insert(*it);
194  }
195  this->comm().broadcast(domains);
196 
197  // 7.)
198  // Read the approximation order(s) of the var-th variable
199  int order=0;
200  if (this->processor_id() == 0)
201  io.data (order);
202  this->comm().broadcast(order);
203 
204 
205  // do the same for infinite element radial_order
206  int rad_order=0;
207  if (read_ifem_info)
208  {
209  if (this->processor_id() == 0)
210  io.data(rad_order);
211  this->comm().broadcast(rad_order);
212  }
213 
214  // Read the finite element type of the var-th variable
215  int fam=0;
216  if (this->processor_id() == 0)
217  io.data (fam);
218  this->comm().broadcast(fam);
219  FEType type;
220  type.order = static_cast<Order>(order);
221  type.family = static_cast<FEFamily>(fam);
222 
223  // Check for incompatibilities. The shape function indexing was
224  // changed for the monomial and xyz finite element families to
225  // simplify extension to arbitrary p. The consequence is that
226  // old restart files will not be read correctly. This is expected
227  // to be an unlikely occurence, but catch it anyway.
228  if (read_legacy_format)
229  if ((type.family == MONOMIAL || type.family == XYZ) &&
230  ((type.order.get_order() > 2 && this->get_mesh().mesh_dimension() == 2) ||
231  (type.order.get_order() > 1 && this->get_mesh().mesh_dimension() == 3)))
232  {
233  libmesh_here();
234  libMesh::out << "*****************************************************************\n"
235  << "* WARNING: reading a potentially incompatible restart file!!! *\n"
236  << "* contact libmesh-users@lists.sourceforge.net for more details *\n"
237  << "*****************************************************************"
238  << std::endl;
239  }
240 
241  // Read additional information for infinite elements
242  int radial_fam=0;
243  int i_map=0;
244  if (read_ifem_info)
245  {
246  if (this->processor_id() == 0)
247  io.data (radial_fam);
248  this->comm().broadcast(radial_fam);
249  if (this->processor_id() == 0)
250  io.data (i_map);
251  this->comm().broadcast(i_map);
252  }
253 
254 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
255 
256  type.radial_order = static_cast<Order>(rad_order);
257  type.radial_family = static_cast<FEFamily>(radial_fam);
258  type.inf_map = static_cast<InfMapType>(i_map);
259 
260 #endif
261 
262  if (read_header_in)
263  {
264  if (domains.empty())
265  _written_var_indices[var] = this->add_variable (var_name, type);
266  else
267  _written_var_indices[var] = this->add_variable (var_name, type, &domains);
268  }
269  else
270  _written_var_indices[var] = this->variable_number(var_name);
271  }
272  }
273 
274  // 8.)
275  // Read the number of additional vectors.
276  unsigned int nvecs=0;
277  if (this->processor_id() == 0)
278  io.data (nvecs);
279  this->comm().broadcast(nvecs);
280 
281  // If nvecs > 0, this means that write_additional_data
282  // was true when this file was written. We will need to
283  // make use of this fact later.
284  this->_additional_data_written = nvecs;
285 
286  for (unsigned int vec=0; vec<nvecs; vec++)
287  {
288  // 9.)
289  // Read the name of the vec-th additional vector
290  std::string vec_name;
291  if (this->processor_id() == 0)
292  io.data (vec_name);
293  this->comm().broadcast(vec_name);
294 
295  if (read_additional_data)
296  {
297  // Systems now can handle adding post-initialization vectors
298  // libmesh_assert(this->_can_add_vectors);
299  // Some systems may have added their own vectors already
300  // libmesh_assert_equal_to (this->_vectors.count(vec_name), 0);
301 
302  this->add_vector(vec_name);
303  }
304  }
305 }
306 
307 
308 
309 #ifdef LIBMESH_ENABLE_DEPRECATED
310 void System::read_legacy_data (Xdr & io,
311  const bool read_additional_data)
312 {
313  libmesh_deprecated();
314 
315  // This method implements the output of the vectors
316  // contained in this System object, embedded in the
317  // output of an EquationSystems<T_sys>.
318  //
319  // 10.) The global solution vector, re-ordered to be node-major
320  // (More on this later.)
321  //
322  // for each additional vector in the object
323  //
324  // 11.) The global additional vector, re-ordered to be
325  // node-major (More on this later.)
326  libmesh_assert (io.reading());
327 
328  // read and reordering buffers
329  std::vector<Number> global_vector;
330  std::vector<Number> reordered_vector;
331 
332  // 10.)
333  // Read and set the solution vector
334  {
335  if (this->processor_id() == 0)
336  io.data (global_vector);
337  this->comm().broadcast(global_vector);
338 
339  // Remember that the stored vector is node-major.
340  // We need to put it into whatever application-specific
341  // ordering we may have using the dof_map.
342  reordered_vector.resize(global_vector.size());
343 
344  //libMesh::out << "global_vector.size()=" << global_vector.size() << std::endl;
345  //libMesh::out << "this->n_dofs()=" << this->n_dofs() << std::endl;
346 
347  libmesh_assert_equal_to (global_vector.size(), this->n_dofs());
348 
349  dof_id_type cnt=0;
350 
351  const unsigned int sys = this->number();
352  const unsigned int nv = cast_int<unsigned int>
353  (this->_written_var_indices.size());
354  libmesh_assert_less_equal (nv, this->n_vars());
355 
356  for (unsigned int data_var=0; data_var<nv; data_var++)
357  {
358  const unsigned int var = _written_var_indices[data_var];
359 
360  // First reorder the nodal DOF values
361  for (auto & node : this->get_mesh().node_ptr_range())
362  for (unsigned int index=0; index<node->n_comp(sys,var); index++)
363  {
364  libmesh_assert_not_equal_to (node->dof_number(sys, var, index),
365  DofObject::invalid_id);
366 
367  libmesh_assert_less (cnt, global_vector.size());
368 
369  reordered_vector[node->dof_number(sys, var, index)] =
370  global_vector[cnt++];
371  }
372 
373  // Then reorder the element DOF values
374  for (auto & elem : this->get_mesh().active_element_ptr_range())
375  for (unsigned int index=0; index<elem->n_comp(sys,var); index++)
376  {
377  libmesh_assert_not_equal_to (elem->dof_number(sys, var, index),
378  DofObject::invalid_id);
379 
380  libmesh_assert_less (cnt, global_vector.size());
381 
382  reordered_vector[elem->dof_number(sys, var, index)] =
383  global_vector[cnt++];
384  }
385  }
386 
387  *(this->solution) = reordered_vector;
388  }
389 
390  // For each additional vector, simply go through the list.
391  // ONLY attempt to do this IF additional data was actually
392  // written to the file for this system (controlled by the
393  // _additional_data_written flag).
394  if (this->_additional_data_written)
395  {
396  const std::size_t nvecs = this->_vectors.size();
397 
398  // If the number of additional vectors written is non-zero, and
399  // the number of additional vectors we have is non-zero, and
400  // they don't match, then something is wrong and we can't be
401  // sure we're reading data into the correct places.
402  if (read_additional_data && nvecs &&
403  nvecs != this->_additional_data_written)
404  libmesh_error_msg
405  ("Additional vectors in file do not match system");
406 
407  std::map<std::string, NumericVector<Number> *>::iterator
408  pos = this->_vectors.begin();
409 
410  for (std::size_t i = 0; i != this->_additional_data_written; ++i)
411  {
412  // 11.)
413  // Read the values of the vec-th additional vector.
414  // Prior do _not_ clear, but fill with zero, since the
415  // additional vectors _have_ to have the same size
416  // as the solution vector
417  std::fill (global_vector.begin(), global_vector.end(), libMesh::zero);
418 
419  if (this->processor_id() == 0)
420  io.data (global_vector);
421 
422  // If read_additional_data==true and we have additional vectors,
423  // then we will keep this vector data; otherwise we are going to
424  // throw it away.
425  if (read_additional_data && nvecs)
426  {
427  this->comm().broadcast(global_vector);
428 
429  // Remember that the stored vector is node-major.
430  // We need to put it into whatever application-specific
431  // ordering we may have using the dof_map.
432  std::fill (reordered_vector.begin(),
433  reordered_vector.end(),
434  libMesh::zero);
435 
436  reordered_vector.resize(global_vector.size());
437 
438  libmesh_assert_equal_to (global_vector.size(), this->n_dofs());
439 
440  dof_id_type cnt=0;
441 
442  const unsigned int sys = this->number();
443  const unsigned int nv = cast_int<unsigned int>
444  (this->_written_var_indices.size());
445  libmesh_assert_less_equal (nv, this->n_vars());
446 
447  for (unsigned int data_var=0; data_var<nv; data_var++)
448  {
449  const unsigned int var = _written_var_indices[data_var];
450  // First reorder the nodal DOF values
451  for (auto & node : this->get_mesh().node_ptr_range())
452  for (unsigned int index=0; index<node->n_comp(sys,var); index++)
453  {
454  libmesh_assert_not_equal_to (node->dof_number(sys, var, index),
455  DofObject::invalid_id);
456 
457  libmesh_assert_less (cnt, global_vector.size());
458 
459  reordered_vector[node->dof_number(sys, var, index)] =
460  global_vector[cnt++];
461  }
462 
463  // Then reorder the element DOF values
464  for (auto & elem : this->get_mesh().active_element_ptr_range())
465  for (unsigned int index=0; index<elem->n_comp(sys,var); index++)
466  {
467  libmesh_assert_not_equal_to (elem->dof_number(sys, var, index),
468  DofObject::invalid_id);
469 
470  libmesh_assert_less (cnt, global_vector.size());
471 
472  reordered_vector[elem->dof_number(sys, var, index)] =
473  global_vector[cnt++];
474  }
475  }
476 
477  // use the overloaded operator=(std::vector) to assign the values
478  *(pos->second) = reordered_vector;
479  }
480 
481  // If we've got vectors then we need to be iterating through
482  // those too
483  if (pos != this->_vectors.end())
484  ++pos;
485  }
486  } // end if (_additional_data_written)
487 }
488 #endif
489 
490 
491 
492 template <typename InValType>
493 void System::read_parallel_data (Xdr & io,
494  const bool read_additional_data)
495 {
515  // PerfLog pl("IO Performance",false);
516  // pl.push("read_parallel_data");
517  dof_id_type total_read_size = 0;
518 
519  libmesh_assert (io.reading());
520  libmesh_assert (io.is_open());
521 
522  // build the ordered nodes and element maps.
523  // when writing/reading parallel files we need to iterate
524  // over our nodes/elements in order of increasing global id().
525  // however, this is not guaranteed to be ordering we obtain
526  // by using the node_iterators/element_iterators directly.
527  // so build a set, sorted by id(), that provides the ordering.
528  // further, for memory economy build the set but then transfer
529  // its contents to vectors, which will be sorted.
530  std::vector<const DofObject *> ordered_nodes, ordered_elements;
531  {
532  std::set<const DofObject *, CompareDofObjectsByID>
533  ordered_nodes_set (this->get_mesh().local_nodes_begin(),
534  this->get_mesh().local_nodes_end());
535 
536  ordered_nodes.insert(ordered_nodes.end(),
537  ordered_nodes_set.begin(),
538  ordered_nodes_set.end());
539  }
540  {
541  std::set<const DofObject *, CompareDofObjectsByID>
542  ordered_elements_set (this->get_mesh().local_elements_begin(),
543  this->get_mesh().local_elements_end());
544 
545  ordered_elements.insert(ordered_elements.end(),
546  ordered_elements_set.begin(),
547  ordered_elements_set.end());
548  }
549 
550  // std::vector<Number> io_buffer;
551  std::vector<InValType> io_buffer;
552 
553  // 9.)
554  //
555  // Actually read the solution components
556  // for the ith system to disk
557  io.data(io_buffer);
558 
559  total_read_size += cast_int<dof_id_type>(io_buffer.size());
560 
561  const unsigned int sys_num = this->number();
562  const unsigned int nv = cast_int<unsigned int>
563  (this->_written_var_indices.size());
564  libmesh_assert_less_equal (nv, this->n_vars());
565 
566  dof_id_type cnt=0;
567 
568  // Loop over each non-SCALAR variable and each node, and read out the value.
569  for (unsigned int data_var=0; data_var<nv; data_var++)
570  {
571  const unsigned int var = _written_var_indices[data_var];
572  if (this->variable(var).type().family != SCALAR)
573  {
574  // First read the node DOF values
575  for (std::vector<const DofObject *>::const_iterator
576  it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
577  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
578  {
579  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
580  DofObject::invalid_id);
581  libmesh_assert_less (cnt, io_buffer.size());
582  this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
583  }
584 
585  // Then read the element DOF values
586  for (std::vector<const DofObject *>::const_iterator
587  it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
588  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
589  {
590  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
591  DofObject::invalid_id);
592  libmesh_assert_less (cnt, io_buffer.size());
593  this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
594  }
595  }
596  }
597 
598  // Finally, read the SCALAR variables on the last processor
599  for (unsigned int data_var=0; data_var<nv; data_var++)
600  {
601  const unsigned int var = _written_var_indices[data_var];
602  if (this->variable(var).type().family == SCALAR)
603  {
604  if (this->processor_id() == (this->n_processors()-1))
605  {
606  const DofMap & dof_map = this->get_dof_map();
607  std::vector<dof_id_type> SCALAR_dofs;
608  dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
609 
610  for (std::size_t i=0; i<SCALAR_dofs.size(); i++)
611  this->solution->set(SCALAR_dofs[i], io_buffer[cnt++]);
612  }
613  }
614  }
615 
616  // And we're done setting solution entries
617  this->solution->close();
618 
619  // For each additional vector, simply go through the list.
620  // ONLY attempt to do this IF additional data was actually
621  // written to the file for this system (controlled by the
622  // _additional_data_written flag).
623  if (this->_additional_data_written)
624  {
625  const std::size_t nvecs = this->_vectors.size();
626 
627  // If the number of additional vectors written is non-zero, and
628  // the number of additional vectors we have is non-zero, and
629  // they don't match, then something is wrong and we can't be
630  // sure we're reading data into the correct places.
631  if (read_additional_data && nvecs &&
632  nvecs != this->_additional_data_written)
633  libmesh_error_msg
634  ("Additional vectors in file do not match system");
635 
636  std::map<std::string, NumericVector<Number> *>::const_iterator
637  pos = _vectors.begin();
638 
639  for (std::size_t i = 0; i != this->_additional_data_written; ++i)
640  {
641  cnt=0;
642  io_buffer.clear();
643 
644  // 10.)
645  //
646  // Actually read the additional vector components
647  // for the ith system from disk
648  io.data(io_buffer);
649 
650  total_read_size += cast_int<dof_id_type>(io_buffer.size());
651 
652  // If read_additional_data==true and we have additional vectors,
653  // then we will keep this vector data; otherwise we are going to
654  // throw it away.
655  if (read_additional_data && nvecs)
656  {
657  // Loop over each non-SCALAR variable and each node, and read out the value.
658  for (unsigned int data_var=0; data_var<nv; data_var++)
659  {
660  const unsigned int var = _written_var_indices[data_var];
661  if (this->variable(var).type().family != SCALAR)
662  {
663  // First read the node DOF values
664  for (std::vector<const DofObject *>::const_iterator
665  it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
666  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
667  {
668  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
669  DofObject::invalid_id);
670  libmesh_assert_less (cnt, io_buffer.size());
671  pos->second->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
672  }
673 
674  // Then read the element DOF values
675  for (std::vector<const DofObject *>::const_iterator
676  it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
677  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
678  {
679  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
680  DofObject::invalid_id);
681  libmesh_assert_less (cnt, io_buffer.size());
682  pos->second->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
683  }
684  }
685  }
686 
687  // Finally, read the SCALAR variables on the last processor
688  for (unsigned int data_var=0; data_var<nv; data_var++)
689  {
690  const unsigned int var = _written_var_indices[data_var];
691  if (this->variable(var).type().family == SCALAR)
692  {
693  if (this->processor_id() == (this->n_processors()-1))
694  {
695  const DofMap & dof_map = this->get_dof_map();
696  std::vector<dof_id_type> SCALAR_dofs;
697  dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
698 
699  for (std::size_t i=0; i<SCALAR_dofs.size(); i++)
700  pos->second->set(SCALAR_dofs[i], io_buffer[cnt++]);
701  }
702  }
703  }
704 
705  // And we're done setting entries for this variable
706  pos->second->close();
707  }
708 
709  // If we've got vectors then we need to be iterating through
710  // those too
711  if (pos != this->_vectors.end())
712  ++pos;
713  }
714  }
715 
716  // const Real
717  // dt = pl.get_elapsed_time(),
718  // rate = total_read_size*sizeof(Number)/dt;
719 
720  // libMesh::err << "Read " << total_read_size << " \"Number\" values\n"
721  // << " Elapsed time = " << dt << '\n'
722  // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
723 
724  // pl.pop("read_parallel_data");
725 }
726 
727 
728 template <typename InValType>
729 void System::read_serialized_data (Xdr & io,
730  const bool read_additional_data)
731 {
732  // This method implements the input of the vectors
733  // contained in this System object, embedded in the
734  // output of an EquationSystems<T_sys>.
735  //
736  // 10.) The global solution vector, re-ordered to be node-major
737  // (More on this later.)
738  //
739  // for each additional vector in the object
740  //
741  // 11.) The global additional vector, re-ordered to be
742  // node-major (More on this later.)
743  parallel_object_only();
744  std::string comment;
745 
746  // PerfLog pl("IO Performance",false);
747  // pl.push("read_serialized_data");
748  // std::size_t total_read_size = 0;
749 
750  // 10.)
751  // Read the global solution vector
752  {
753  // total_read_size +=
754  this->read_serialized_vector<InValType>(io, this->solution.get());
755 
756  // get the comment
757  if (this->processor_id() == 0)
758  io.comment (comment);
759  }
760 
761  // 11.)
762  // Only read additional vectors if data is available, and only use
763  // that data to fill our vectors if the user requested it.
764  if (this->_additional_data_written)
765  {
766  const std::size_t nvecs = this->_vectors.size();
767 
768  // If the number of additional vectors written is non-zero, and
769  // the number of additional vectors we have is non-zero, and
770  // they don't match, then we can't read additional vectors
771  // and be sure we're reading data into the correct places.
772  if (read_additional_data && nvecs &&
773  nvecs != this->_additional_data_written)
774  libmesh_error_msg
775  ("Additional vectors in file do not match system");
776 
777  std::map<std::string, NumericVector<Number> *>::const_iterator
778  pos = _vectors.begin();
779 
780  for (std::size_t i = 0; i != this->_additional_data_written; ++i)
781  {
782  // Read data, but only put it into a vector if we've been
783  // asked to and if we have a corresponding vector to read.
784 
785  // total_read_size +=
786  this->read_serialized_vector<InValType>
787  (io, (read_additional_data && nvecs) ? pos->second : libmesh_nullptr);
788 
789  // get the comment
790  if (this->processor_id() == 0)
791  io.comment (comment);
792 
793 
794  // If we've got vectors then we need to be iterating through
795  // those too
796  if (pos != this->_vectors.end())
797  ++pos;
798  }
799  }
800 
801  // const Real
802  // dt = pl.get_elapsed_time(),
803  // rate = total_read_size*sizeof(Number)/dt;
804 
805  // libMesh::out << "Read " << total_read_size << " \"Number\" values\n"
806  // << " Elapsed time = " << dt << '\n'
807  // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
808 
809  // pl.pop("read_serialized_data");
810 }
811 
812 
813 
814 template <typename iterator_type, typename InValType>
815 std::size_t System::read_serialized_blocked_dof_objects (const dof_id_type n_objs,
816  const iterator_type begin,
817  const iterator_type end,
818  const InValType ,
819  Xdr & io,
820  const std::vector<NumericVector<Number> *> & vecs,
821  const unsigned int var_to_read) const
822 {
823  //-------------------------------------------------------
824  // General order: (IO format 0.7.4 & greater)
825  //
826  // for (objects ...)
827  // for (vecs ....)
828  // for (vars ....)
829  // for (comps ...)
830  //
831  // where objects are nodes or elements, sorted to be
832  // partition independent,
833  // vecs are one or more *identically distributed* solution
834  // coefficient vectors, vars are one or more variables
835  // to write, and comps are all the components for said
836  // vars on the object.
837 
838  typedef std::vector<NumericVector<Number> *>::const_iterator vec_iterator_type;
839 
840  // variables to read. Unless specified otherwise, defaults to _written_var_indices.
841  std::vector<unsigned int> vars_to_read (_written_var_indices);
842 
843  if (var_to_read != libMesh::invalid_uint)
844  {
845  vars_to_read.clear();
846  vars_to_read.push_back(var_to_read);
847  }
848 
849  const unsigned int
850  sys_num = this->number(),
851  num_vecs = cast_int<unsigned int>(vecs.size());
852  const dof_id_type
853  io_blksize = cast_int<dof_id_type>(std::min(max_io_blksize, static_cast<std::size_t>(n_objs))),
854  num_blks = cast_int<unsigned int>(std::ceil(static_cast<double>(n_objs)/
855  static_cast<double>(io_blksize)));
856 
857  libmesh_assert_less_equal (_written_var_indices.size(), this->n_vars());
858 
859  std::size_t n_read_values=0;
860 
861  std::vector<std::vector<dof_id_type>> xfer_ids(num_blks); // The global IDs and # of components for the local objects in all blocks
862  std::vector<std::vector<Number>> recv_vals(num_blks); // The raw values for the local objects in all blocks
863  std::vector<Parallel::Request>
864  id_requests(num_blks), val_requests(num_blks);
865 
866  // ------------------------------------------------------
867  // First pass - count the number of objects in each block
868  // traverse all the objects and figure out which block they
869  // will ultimately live in.
870  std::vector<std::size_t>
871  xfer_ids_size (num_blks,0),
872  recv_vals_size (num_blks,0);
873 
874 
875  for (iterator_type it=begin; it!=end; ++it)
876  {
877  const dof_id_type
878  id = (*it)->id(),
879  block = id/io_blksize;
880 
881  libmesh_assert_less (block, num_blks);
882 
883  xfer_ids_size[block] += 2; // for each object, we send its id, as well as the total number of components for all variables
884 
885  dof_id_type n_comp_tot=0;
886  for (std::vector<unsigned int>::const_iterator var_it=vars_to_read.begin();
887  var_it!=vars_to_read.end(); ++var_it)
888  n_comp_tot += (*it)->n_comp(sys_num, *var_it); // for each variable, we will receive the nonzero components
889 
890  recv_vals_size[block] += n_comp_tot*num_vecs;
891  }
892 
893  // knowing the recv_vals_size[block] for each processor allows
894  // us to sum them and find the global size for each block.
895  std::vector<std::size_t> tot_vals_size(recv_vals_size);
896  this->comm().sum (tot_vals_size);
897 
898 
899  //------------------------------------------
900  // Collect the ids & number of values needed
901  // for all local objects, binning them into
902  // 'blocks' that will be sent to processor 0
903  for (dof_id_type blk=0; blk<num_blks; blk++)
904  {
905  // Each processor should build up its transfer buffers for its
906  // local objects in [first_object,last_object).
907  const dof_id_type
908  first_object = blk*io_blksize,
909  last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs);
910 
911  // convenience
912  std::vector<dof_id_type> & ids (xfer_ids[blk]);
913  std::vector<Number> & vals (recv_vals[blk]);
914 
915  // we now know the number of values we will store for each block,
916  // so we can do efficient preallocation
917  ids.clear(); ids.reserve (xfer_ids_size[blk]);
918  vals.resize(recv_vals_size[blk]);
919 
920  if (recv_vals_size[blk] != 0) // only if there are nonzero values to receive
921  for (iterator_type it=begin; it!=end; ++it)
922  if (((*it)->id() >= first_object) && // object in [first_object,last_object)
923  ((*it)->id() < last_object))
924  {
925  ids.push_back((*it)->id());
926 
927  unsigned int n_comp_tot=0;
928 
929  for (std::vector<unsigned int>::const_iterator var_it=vars_to_read.begin();
930  var_it!=vars_to_read.end(); ++var_it)
931  n_comp_tot += (*it)->n_comp(sys_num,*var_it);
932 
933  ids.push_back (n_comp_tot*num_vecs);
934  }
935 
936 #ifdef LIBMESH_HAVE_MPI
937  Parallel::MessageTag id_tag = this->comm().get_unique_tag(100*num_blks + blk);
938  Parallel::MessageTag val_tag = this->comm().get_unique_tag(200*num_blks + blk);
939 
940  // nonblocking send the data for this block
941  this->comm().send (0, ids, id_requests[blk], id_tag);
942 
943  // Go ahead and post the receive too
944  this->comm().receive (0, vals, val_requests[blk], val_tag);
945 #endif
946  }
947 
948  //---------------------------------------------------
949  // Here processor 0 will read and distribute the data.
950  // We have to do this block-wise to ensure that we
951  // do not exhaust memory on processor 0.
952 
953  // give these variables scope outside the block to avoid reallocation
954  std::vector<std::vector<dof_id_type>> recv_ids (this->n_processors());
955  std::vector<std::vector<Number>> send_vals (this->n_processors());
956  std::vector<Parallel::Request> reply_requests (this->n_processors());
957  std::vector<unsigned int> obj_val_offsets; // map to traverse entry-wise rather than processor-wise
958  std::vector<Number> input_vals; // The input buffer for the current block
959  std::vector<InValType> input_vals_tmp; // The input buffer for the current block
960 
961  for (dof_id_type blk=0; blk<num_blks; blk++)
962  {
963  // Each processor should build up its transfer buffers for its
964  // local objects in [first_object,last_object).
965  const dof_id_type
966  first_object = blk*io_blksize,
967  last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs),
968  n_objects_blk = last_object - first_object;
969 
970  // Processor 0 has a special job. It needs to gather the requested indices
971  // in [first_object,last_object) from all processors, read the data from
972  // disk, and reply
973  if (this->processor_id() == 0)
974  {
975  // we know the input buffer size for this block and can begin reading it now
976  input_vals.resize(tot_vals_size[blk]);
977  input_vals_tmp.resize(tot_vals_size[blk]);
978 
979  // a ThreadedIO object to perform asynchronous file IO
980  ThreadedIO<InValType> threaded_io(io, input_vals_tmp);
981  Threads::Thread async_io(threaded_io);
982 
983  Parallel::MessageTag id_tag = this->comm().get_unique_tag(100*num_blks + blk);
984  Parallel::MessageTag val_tag = this->comm().get_unique_tag(200*num_blks + blk);
985 
986  // offset array. this will define where each object's values
987  // map into the actual input_vals buffer. this must get
988  // 0-initialized because 0-component objects are not actually sent
989  obj_val_offsets.resize (n_objects_blk); std::fill (obj_val_offsets.begin(), obj_val_offsets.end(), 0);
990  recv_vals_size.resize(this->n_processors()); // reuse this to count how many values are going to each processor
991 
992 #ifndef NDEBUG
993  std::size_t n_vals_blk = 0;
994 #endif
995 
996  // loop over all processors and process their index request
997  for (unsigned int comm_step=0; comm_step<this->n_processors(); comm_step++)
998  {
999 #ifdef LIBMESH_HAVE_MPI
1000  // blocking receive indices for this block, imposing no particular order on processor
1001  Parallel::Status id_status (this->comm().probe (Parallel::any_source, id_tag));
1002  std::vector<dof_id_type> & ids (recv_ids[id_status.source()]);
1003  std::size_t & n_vals_proc (recv_vals_size[id_status.source()]);
1004  this->comm().receive (id_status.source(), ids, id_tag);
1005 #else
1006  // straight copy without MPI
1007  std::vector<dof_id_type> & ids (recv_ids[0]);
1008  std::size_t & n_vals_proc (recv_vals_size[0]);
1009  ids = xfer_ids[blk];
1010 #endif
1011 
1012  n_vals_proc = 0;
1013 
1014  // note its possible we didn't receive values for objects in
1015  // this block if they have no components allocated.
1016  for (std::size_t idx=0; idx<ids.size(); idx+=2)
1017  {
1018  const dof_id_type
1019  local_idx = ids[idx+0]-first_object,
1020  n_vals_tot_allvecs = ids[idx+1];
1021 
1022  libmesh_assert_less (local_idx, n_objects_blk);
1023 
1024  obj_val_offsets[local_idx] = n_vals_tot_allvecs;
1025  n_vals_proc += n_vals_tot_allvecs;
1026  }
1027 
1028 #ifndef NDEBUG
1029  n_vals_blk += n_vals_proc;
1030 #endif
1031  }
1032 
1033  // We need the offsets into the input_vals vector for each object.
1034  // fortunately, this is simply the partial sum of the total number
1035  // of components for each object
1036  std::partial_sum(obj_val_offsets.begin(), obj_val_offsets.end(),
1037  obj_val_offsets.begin());
1038 
1039  libmesh_assert_equal_to (n_vals_blk, obj_val_offsets.back());
1040  libmesh_assert_equal_to (n_vals_blk, tot_vals_size[blk]);
1041 
1042  // Wait for read completion
1043  async_io.join();
1044  // now copy the values back to the main vector for transfer
1045  for (std::size_t i_val=0; i_val<input_vals.size(); i_val++)
1046  input_vals[i_val] = input_vals_tmp[i_val];
1047 
1048  n_read_values += input_vals.size();
1049 
1050  // pack data replies for each processor
1051  for (processor_id_type proc=0; proc<this->n_processors(); proc++)
1052  {
1053  const std::vector<dof_id_type> & ids (recv_ids[proc]);
1054  std::vector<Number> & vals (send_vals[proc]);
1055  const std::size_t & n_vals_proc (recv_vals_size[proc]);
1056 
1057  vals.clear(); vals.reserve(n_vals_proc);
1058 
1059  for (std::size_t idx=0; idx<ids.size(); idx+=2)
1060  {
1061  const dof_id_type
1062  local_idx = ids[idx+0]-first_object,
1063  n_vals_tot_allvecs = ids[idx+1];
1064 
1065  std::vector<Number>::const_iterator in_vals(input_vals.begin());
1066  if (local_idx != 0)
1067  std::advance (in_vals, obj_val_offsets[local_idx-1]);
1068 
1069  for (unsigned int val=0; val<n_vals_tot_allvecs; val++, ++in_vals)
1070  {
1071  libmesh_assert (in_vals != input_vals.end());
1072  //libMesh::out << "*in_vals=" << *in_vals << '\n';
1073  vals.push_back(*in_vals);
1074  }
1075  }
1076 
1077 #ifdef LIBMESH_HAVE_MPI
1078  // send the relevant values to this processor
1079  this->comm().send (proc, vals, reply_requests[proc], val_tag);
1080 #else
1081  recv_vals[blk] = vals;
1082 #endif
1083  }
1084  } // end processor 0 read/reply
1085 
1086  // all processors complete the (already posted) read for this block
1087  {
1088  Parallel::wait (val_requests[blk]);
1089 
1090  const std::vector<Number> & vals (recv_vals[blk]);
1091  std::vector<Number>::const_iterator val_it(vals.begin());
1092 
1093  if (!recv_vals[blk].empty()) // nonzero values to receive
1094  for (iterator_type it=begin; it!=end; ++it)
1095  if (((*it)->id() >= first_object) && // object in [first_object,last_object)
1096  ((*it)->id() < last_object))
1097  // unpack & set the values
1098  for (vec_iterator_type vec_it=vecs.begin(); vec_it!=vecs.end(); ++vec_it)
1099  {
1100  NumericVector<Number> * vec(*vec_it);
1101 
1102  for (std::vector<unsigned int>::const_iterator var_it=vars_to_read.begin();
1103  var_it!=vars_to_read.end(); ++var_it)
1104  {
1105  const unsigned int n_comp = (*it)->n_comp(sys_num,*var_it);
1106 
1107  for (unsigned int comp=0; comp<n_comp; comp++, ++val_it)
1108  {
1109  const dof_id_type dof_index = (*it)->dof_number (sys_num, *var_it, comp);
1110  libmesh_assert (val_it != vals.end());
1111  if (vec)
1112  {
1113  libmesh_assert_greater_equal (dof_index, vec->first_local_index());
1114  libmesh_assert_less (dof_index, vec->last_local_index());
1115  //libMesh::out << "dof_index, *val_it = \t" << dof_index << ", " << *val_it << '\n';
1116  vec->set (dof_index, *val_it);
1117  }
1118  }
1119  }
1120  }
1121  }
1122 
1123  // processor 0 needs to make sure all replies have been handed off
1124  if (this->processor_id () == 0)
1125  Parallel::wait(reply_requests);
1126  }
1127 
1128  return n_read_values;
1129 }
1130 
1131 
1132 
1133 unsigned int System::read_SCALAR_dofs (const unsigned int var,
1134  Xdr & io,
1135  NumericVector<Number> * vec) const
1136 {
1137  unsigned int n_assigned_vals = 0; // the number of values assigned, this will be returned.
1138 
1139  // Processor 0 will read the block from the buffer stream and send it to the last processor
1140  const unsigned int n_SCALAR_dofs = this->variable(var).type().order.get_order();
1141  std::vector<Number> input_buffer(n_SCALAR_dofs);
1142  if (this->processor_id() == 0)
1143  {
1144  io.data_stream(&input_buffer[0], n_SCALAR_dofs);
1145  }
1146 
1147 #ifdef LIBMESH_HAVE_MPI
1148  if (this->n_processors() > 1)
1149  {
1150  const Parallel::MessageTag val_tag = this->comm().get_unique_tag(321);
1151 
1152  // Post the receive on the last processor
1153  if (this->processor_id() == (this->n_processors()-1))
1154  this->comm().receive(0, input_buffer, val_tag);
1155 
1156  // Send the data to processor 0
1157  if (this->processor_id() == 0)
1158  this->comm().send(this->n_processors()-1, input_buffer, val_tag);
1159  }
1160 #endif
1161 
1162  // Finally, set the SCALAR values
1163  if (this->processor_id() == (this->n_processors()-1))
1164  {
1165  const DofMap & dof_map = this->get_dof_map();
1166  std::vector<dof_id_type> SCALAR_dofs;
1167  dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
1168 
1169  for (std::size_t i=0; i<SCALAR_dofs.size(); i++)
1170  {
1171  if (vec)
1172  vec->set (SCALAR_dofs[i], input_buffer[i]);
1173  ++n_assigned_vals;
1174  }
1175  }
1176 
1177  return n_assigned_vals;
1178 }
1179 
1180 
1181 template <typename InValType>
1182 numeric_index_type System::read_serialized_vector (Xdr & io,
1183  NumericVector<Number> * vec)
1184 {
1185  parallel_object_only();
1186 
1187 #ifndef NDEBUG
1188  // In parallel we better be reading a parallel vector -- if not
1189  // we will not set all of its components below!!
1190  if (this->n_processors() > 1 && vec)
1191  {
1192  libmesh_assert (vec->type() == PARALLEL ||
1193  vec->type() == GHOSTED);
1194  }
1195 #endif
1196 
1197  libmesh_assert (io.reading());
1198 
1199  // vector length
1200  unsigned int vector_length=0; // FIXME? size_t would break binary compatibility...
1201 #ifndef NDEBUG
1202  std::size_t n_assigned_vals=0;
1203 #endif
1204 
1205  // Get the buffer size
1206  if (this->processor_id() == 0)
1207  io.data(vector_length, "# vector length");
1208  this->comm().broadcast(vector_length);
1209 
1210  const unsigned int nv = cast_int<unsigned int>
1211  (this->_written_var_indices.size());
1212  const dof_id_type
1213  n_nodes = this->get_mesh().n_nodes(),
1214  n_elem = this->get_mesh().n_elem();
1215 
1216  libmesh_assert_less_equal (nv, this->n_vars());
1217 
1218  // for newer versions, read variables node/elem major
1219  if (io.version() >= LIBMESH_VERSION_ID(0,7,4))
1220  {
1221  //---------------------------------
1222  // Collect the values for all nodes
1223 #ifndef NDEBUG
1224  n_assigned_vals +=
1225 #endif
1226  this->read_serialized_blocked_dof_objects (n_nodes,
1227  this->get_mesh().local_nodes_begin(),
1228  this->get_mesh().local_nodes_end(),
1229  InValType(),
1230  io,
1231  std::vector<NumericVector<Number> *> (1,vec));
1232 
1233 
1234  //------------------------------------
1235  // Collect the values for all elements
1236 #ifndef NDEBUG
1237  n_assigned_vals +=
1238 #endif
1239  this->read_serialized_blocked_dof_objects (n_elem,
1240  this->get_mesh().local_elements_begin(),
1241  this->get_mesh().local_elements_end(),
1242  InValType(),
1243  io,
1244  std::vector<NumericVector<Number> *> (1,vec));
1245  }
1246 
1247  // for older versions, read variables var-major
1248  else
1249  {
1250  // Loop over each variable in the system, and then each node/element in the mesh.
1251  for (unsigned int data_var=0; data_var<nv; data_var++)
1252  {
1253  const unsigned int var = _written_var_indices[data_var];
1254  if (this->variable(var).type().family != SCALAR)
1255  {
1256  //---------------------------------
1257  // Collect the values for all nodes
1258 #ifndef NDEBUG
1259  n_assigned_vals +=
1260 #endif
1261  this->read_serialized_blocked_dof_objects (n_nodes,
1262  this->get_mesh().local_nodes_begin(),
1263  this->get_mesh().local_nodes_end(),
1264  InValType(),
1265  io,
1266  std::vector<NumericVector<Number> *> (1,vec),
1267  var);
1268 
1269 
1270  //------------------------------------
1271  // Collect the values for all elements
1272 #ifndef NDEBUG
1273  n_assigned_vals +=
1274 #endif
1275  this->read_serialized_blocked_dof_objects (n_elem,
1276  this->get_mesh().local_elements_begin(),
1277  this->get_mesh().local_elements_end(),
1278  InValType(),
1279  io,
1280  std::vector<NumericVector<Number> *> (1,vec),
1281  var);
1282  } // end variable loop
1283  }
1284  }
1285 
1286  //-------------------------------------------
1287  // Finally loop over all the SCALAR variables
1288  for (unsigned int data_var=0; data_var<nv; data_var++)
1289  {
1290  const unsigned int var = _written_var_indices[data_var];
1291  if (this->variable(var).type().family == SCALAR)
1292  {
1293 #ifndef NDEBUG
1294  n_assigned_vals +=
1295 #endif
1296  this->read_SCALAR_dofs (var, io, vec);
1297  }
1298  }
1299 
1300  if (vec)
1301  vec->close();
1302 
1303 #ifndef NDEBUG
1304  this->comm().sum (n_assigned_vals);
1305  libmesh_assert_equal_to (n_assigned_vals, vector_length);
1306 #endif
1307 
1308  return vector_length;
1309 }
1310 
1311 
1312 
1313 void System::write_header (Xdr & io,
1314  const std::string & /* version is currently unused */,
1315  const bool write_additional_data) const
1316 {
1350  libmesh_assert (io.writing());
1351 
1352 
1353  // Only write the header information
1354  // if we are processor 0.
1355  if (this->get_mesh().processor_id() != 0)
1356  return;
1357 
1358  std::string comment;
1359  char buf[80];
1360 
1361  // 5.)
1362  // Write the number of variables in the system
1363 
1364  {
1365  // set up the comment
1366  comment = "# No. of Variables in System \"";
1367  comment += this->name();
1368  comment += "\"";
1369 
1370  unsigned int nv = this->n_vars();
1371  io.data (nv, comment.c_str());
1372  }
1373 
1374 
1375  for (unsigned int var=0; var<this->n_vars(); var++)
1376  {
1377  // 6.)
1378  // Write the name of the var-th variable
1379  {
1380  // set up the comment
1381  comment = "# Name, Variable No. ";
1382  std::sprintf(buf, "%u", var);
1383  comment += buf;
1384  comment += ", System \"";
1385  comment += this->name();
1386  comment += "\"";
1387 
1388  std::string var_name = this->variable_name(var);
1389  io.data (var_name, comment.c_str());
1390  }
1391 
1392  // 6.1.) Variable subdomains
1393  {
1394  // set up the comment
1395  comment = "# Subdomains, Variable \"";
1396  std::sprintf(buf, "%s", this->variable_name(var).c_str());
1397  comment += buf;
1398  comment += "\", System \"";
1399  comment += this->name();
1400  comment += "\"";
1401 
1402  const std::set<subdomain_id_type> & domains = this->variable(var).active_subdomains();
1403  std::vector<subdomain_id_type> domain_array;
1404  domain_array.assign(domains.begin(), domains.end());
1405  io.data (domain_array, comment.c_str());
1406  }
1407 
1408  // 7.)
1409  // Write the approximation order of the var-th variable
1410  // in this system
1411  {
1412  // set up the comment
1413  comment = "# Approximation Order, Variable \"";
1414  std::sprintf(buf, "%s", this->variable_name(var).c_str());
1415  comment += buf;
1416  comment += "\", System \"";
1417  comment += this->name();
1418  comment += "\"";
1419 
1420  int order = static_cast<int>(this->variable_type(var).order);
1421  io.data (order, comment.c_str());
1422  }
1423 
1424 
1425 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1426 
1427  // do the same for radial_order
1428  {
1429  comment = "# Radial Approximation Order, Variable \"";
1430  std::sprintf(buf, "%s", this->variable_name(var).c_str());
1431  comment += buf;
1432  comment += "\", System \"";
1433  comment += this->name();
1434  comment += "\"";
1435 
1436  int rad_order = static_cast<int>(this->variable_type(var).radial_order);
1437  io.data (rad_order, comment.c_str());
1438  }
1439 
1440 #endif
1441 
1442  // Write the Finite Element type of the var-th variable
1443  // in this System
1444  {
1445  // set up the comment
1446  comment = "# FE Family, Variable \"";
1447  std::sprintf(buf, "%s", this->variable_name(var).c_str());
1448  comment += buf;
1449  comment += "\", System \"";
1450  comment += this->name();
1451  comment += "\"";
1452 
1453  const FEType & type = this->variable_type(var);
1454  int fam = static_cast<int>(type.family);
1455  io.data (fam, comment.c_str());
1456 
1457 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1458 
1459  comment = "# Radial FE Family, Variable \"";
1460  std::sprintf(buf, "%s", this->variable_name(var).c_str());
1461  comment += buf;
1462  comment += "\", System \"";
1463  comment += this->name();
1464  comment += "\"";
1465 
1466  int radial_fam = static_cast<int>(type.radial_family);
1467  io.data (radial_fam, comment.c_str());
1468 
1469  comment = "# Infinite Mapping Type, Variable \"";
1470  std::sprintf(buf, "%s", this->variable_name(var).c_str());
1471  comment += buf;
1472  comment += "\", System \"";
1473  comment += this->name();
1474  comment += "\"";
1475 
1476  int i_map = static_cast<int>(type.inf_map);
1477  io.data (i_map, comment.c_str());
1478 #endif
1479  }
1480  } // end of the variable loop
1481 
1482  // 8.)
1483  // Write the number of additional vectors in the System.
1484  // If write_additional_data==false, then write zero for
1485  // the number of additional vectors.
1486  {
1487  {
1488  // set up the comment
1489  comment = "# No. of Additional Vectors, System \"";
1490  comment += this->name();
1491  comment += "\"";
1492 
1493  unsigned int nvecs = write_additional_data ? this->n_vectors () : 0;
1494  io.data (nvecs, comment.c_str());
1495  }
1496 
1497  if (write_additional_data)
1498  {
1499  std::map<std::string, NumericVector<Number> *>::const_iterator
1500  vec_pos = this->_vectors.begin();
1501  unsigned int cnt=0;
1502 
1503  for (; vec_pos != this->_vectors.end(); ++vec_pos)
1504  {
1505  // 9.)
1506  // write the name of the cnt-th additional vector
1507  comment = "# Name of ";
1508  std::sprintf(buf, "%d", cnt++);
1509  comment += buf;
1510  comment += "th vector";
1511  std::string vec_name = vec_pos->first;
1512 
1513  io.data (vec_name, comment.c_str());
1514  }
1515  }
1516  }
1517 }
1518 
1519 
1520 
1521 void System::write_parallel_data (Xdr & io,
1522  const bool write_additional_data) const
1523 {
1543  // PerfLog pl("IO Performance",false);
1544  // pl.push("write_parallel_data");
1545  // std::size_t total_written_size = 0;
1546 
1547  std::string comment;
1548 
1549  libmesh_assert (io.writing());
1550 
1551  std::vector<Number> io_buffer; io_buffer.reserve(this->solution->local_size());
1552 
1553  // build the ordered nodes and element maps.
1554  // when writing/reading parallel files we need to iterate
1555  // over our nodes/elements in order of increasing global id().
1556  // however, this is not guaranteed to be ordering we obtain
1557  // by using the node_iterators/element_iterators directly.
1558  // so build a set, sorted by id(), that provides the ordering.
1559  // further, for memory economy build the set but then transfer
1560  // its contents to vectors, which will be sorted.
1561  std::vector<const DofObject *> ordered_nodes, ordered_elements;
1562  {
1563  std::set<const DofObject *, CompareDofObjectsByID>
1564  ordered_nodes_set (this->get_mesh().local_nodes_begin(),
1565  this->get_mesh().local_nodes_end());
1566 
1567  ordered_nodes.insert(ordered_nodes.end(),
1568  ordered_nodes_set.begin(),
1569  ordered_nodes_set.end());
1570  }
1571  {
1572  std::set<const DofObject *, CompareDofObjectsByID>
1573  ordered_elements_set (this->get_mesh().local_elements_begin(),
1574  this->get_mesh().local_elements_end());
1575 
1576  ordered_elements.insert(ordered_elements.end(),
1577  ordered_elements_set.begin(),
1578  ordered_elements_set.end());
1579  }
1580 
1581  const unsigned int sys_num = this->number();
1582  const unsigned int nv = this->n_vars();
1583 
1584  // Loop over each non-SCALAR variable and each node, and write out the value.
1585  for (unsigned int var=0; var<nv; var++)
1586  if (this->variable(var).type().family != SCALAR)
1587  {
1588  // First write the node DOF values
1589  for (std::vector<const DofObject *>::const_iterator
1590  it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
1591  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
1592  {
1593  //libMesh::out << "(*it)->id()=" << (*it)->id() << std::endl;
1594  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
1595  DofObject::invalid_id);
1596 
1597  io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
1598  }
1599 
1600  // Then write the element DOF values
1601  for (std::vector<const DofObject *>::const_iterator
1602  it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
1603  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
1604  {
1605  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
1606  DofObject::invalid_id);
1607 
1608  io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
1609  }
1610  }
1611 
1612  // Finally, write the SCALAR data on the last processor
1613  for (unsigned int var=0; var<this->n_vars(); var++)
1614  if (this->variable(var).type().family == SCALAR)
1615  {
1616  if (this->processor_id() == (this->n_processors()-1))
1617  {
1618  const DofMap & dof_map = this->get_dof_map();
1619  std::vector<dof_id_type> SCALAR_dofs;
1620  dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
1621 
1622  for (std::size_t i=0; i<SCALAR_dofs.size(); i++)
1623  io_buffer.push_back((*this->solution)(SCALAR_dofs[i]));
1624  }
1625  }
1626 
1627  // 9.)
1628  //
1629  // Actually write the reordered solution vector
1630  // for the ith system to disk
1631 
1632  // set up the comment
1633  {
1634  comment = "# System \"";
1635  comment += this->name();
1636  comment += "\" Solution Vector";
1637  }
1638 
1639  io.data (io_buffer, comment.c_str());
1640 
1641  // total_written_size += io_buffer.size();
1642 
1643  // Only write additional vectors if wanted
1644  if (write_additional_data)
1645  {
1646  std::map<std::string, NumericVector<Number> *>::const_iterator
1647  pos = _vectors.begin();
1648 
1649  for (; pos != this->_vectors.end(); ++pos)
1650  {
1651  io_buffer.clear(); io_buffer.reserve( pos->second->local_size());
1652 
1653  // Loop over each non-SCALAR variable and each node, and write out the value.
1654  for (unsigned int var=0; var<nv; var++)
1655  if (this->variable(var).type().family != SCALAR)
1656  {
1657  // First write the node DOF values
1658  for (std::vector<const DofObject *>::const_iterator
1659  it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
1660  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
1661  {
1662  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
1663  DofObject::invalid_id);
1664 
1665  io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
1666  }
1667 
1668  // Then write the element DOF values
1669  for (std::vector<const DofObject *>::const_iterator
1670  it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
1671  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
1672  {
1673  libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
1674  DofObject::invalid_id);
1675 
1676  io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
1677  }
1678  }
1679 
1680  // Finally, write the SCALAR data on the last processor
1681  for (unsigned int var=0; var<this->n_vars(); var++)
1682  if (this->variable(var).type().family == SCALAR)
1683  {
1684  if (this->processor_id() == (this->n_processors()-1))
1685  {
1686  const DofMap & dof_map = this->get_dof_map();
1687  std::vector<dof_id_type> SCALAR_dofs;
1688  dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
1689 
1690  for (std::size_t i=0; i<SCALAR_dofs.size(); i++)
1691  io_buffer.push_back((*pos->second)(SCALAR_dofs[i]));
1692  }
1693  }
1694 
1695  // 10.)
1696  //
1697  // Actually write the reordered additional vector
1698  // for this system to disk
1699 
1700  // set up the comment
1701  {
1702  comment = "# System \"";
1703  comment += this->name();
1704  comment += "\" Additional Vector \"";
1705  comment += pos->first;
1706  comment += "\"";
1707  }
1708 
1709  io.data (io_buffer, comment.c_str());
1710 
1711  // total_written_size += io_buffer.size();
1712  }
1713  }
1714 
1715  // const Real
1716  // dt = pl.get_elapsed_time(),
1717  // rate = total_written_size*sizeof(Number)/dt;
1718 
1719  // libMesh::err << "Write " << total_written_size << " \"Number\" values\n"
1720  // << " Elapsed time = " << dt << '\n'
1721  // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
1722 
1723  // pl.pop("write_parallel_data");
1724 }
1725 
1726 
1727 
1728 void System::write_serialized_data (Xdr & io,
1729  const bool write_additional_data) const
1730 {
1744  parallel_object_only();
1745  std::string comment;
1746 
1747  // PerfLog pl("IO Performance",false);
1748  // pl.push("write_serialized_data");
1749  // std::size_t total_written_size = 0;
1750 
1751  // total_written_size +=
1752  this->write_serialized_vector(io, *this->solution);
1753 
1754  // set up the comment
1755  if (this->processor_id() == 0)
1756  {
1757  comment = "# System \"";
1758  comment += this->name();
1759  comment += "\" Solution Vector";
1760 
1761  io.comment (comment);
1762  }
1763 
1764  // Only write additional vectors if wanted
1765  if (write_additional_data)
1766  {
1767  std::map<std::string, NumericVector<Number> *>::const_iterator
1768  pos = _vectors.begin();
1769 
1770  for (; pos != this->_vectors.end(); ++pos)
1771  {
1772  // total_written_size +=
1773  this->write_serialized_vector(io, *pos->second);
1774 
1775  // set up the comment
1776  if (this->processor_id() == 0)
1777  {
1778  comment = "# System \"";
1779  comment += this->name();
1780  comment += "\" Additional Vector \"";
1781  comment += pos->first;
1782  comment += "\"";
1783  io.comment (comment);
1784  }
1785  }
1786  }
1787 
1788  // const Real
1789  // dt = pl.get_elapsed_time(),
1790  // rate = total_written_size*sizeof(Number)/dt;
1791 
1792  // libMesh::out << "Write " << total_written_size << " \"Number\" values\n"
1793  // << " Elapsed time = " << dt << '\n'
1794  // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
1795 
1796  // pl.pop("write_serialized_data");
1797 
1798 
1799 
1800 
1801  // // test the new method
1802  // {
1803  // std::vector<std::string> names;
1804  // std::vector<NumericVector<Number> *> vectors_to_write;
1805 
1806  // names.push_back("Solution Vector");
1807  // vectors_to_write.push_back(this->solution.get());
1808 
1809  // // Only write additional vectors if wanted
1810  // if (write_additional_data)
1811  // {
1812  // std::map<std::string, NumericVector<Number> *>::const_iterator
1813  // pos = _vectors.begin();
1814 
1815  // for (; pos != this->_vectors.end(); ++pos)
1816  // {
1817  // names.push_back("Additional Vector " + pos->first);
1818  // vectors_to_write.push_back(pos->second);
1819  // }
1820  // }
1821 
1822  // total_written_size =
1823  // this->write_serialized_vectors (io, names, vectors_to_write);
1824 
1825  // const Real
1826  // dt2 = pl.get_elapsed_time(),
1827  // rate2 = total_written_size*sizeof(Number)/(dt2-dt);
1828 
1829  // libMesh::out << "Write (new) " << total_written_size << " \"Number\" values\n"
1830  // << " Elapsed time = " << (dt2-dt) << '\n'
1831  // << " Rate = " << rate2/1.e6 << "(MB/sec)\n\n";
1832 
1833  // }
1834 }
1835 
1836 
1837 
1838 template <typename iterator_type>
1839 std::size_t System::write_serialized_blocked_dof_objects (const std::vector<const NumericVector<Number> *> & vecs,
1840  const dof_id_type n_objs,
1841  const iterator_type begin,
1842  const iterator_type end,
1843  Xdr & io,
1844  const unsigned int var_to_write) const
1845 {
1846  //-------------------------------------------------------
1847  // General order: (IO format 0.7.4 & greater)
1848  //
1849  // for (objects ...)
1850  // for (vecs ....)
1851  // for (vars ....)
1852  // for (comps ...)
1853  //
1854  // where objects are nodes or elements, sorted to be
1855  // partition independent,
1856  // vecs are one or more *identically distributed* solution
1857  // coefficient vectors, vars are one or more variables
1858  // to write, and comps are all the components for said
1859  // vars on the object.
1860 
1861  typedef std::vector<const NumericVector<Number> *>::const_iterator vec_iterator_type;
1862 
1863  // We will write all variables unless requested otherwise.
1864  std::vector<unsigned int> vars_to_write(1, var_to_write);
1865 
1866  if (var_to_write == libMesh::invalid_uint)
1867  {
1868  vars_to_write.clear(); vars_to_write.reserve(this->n_vars());
1869  for (unsigned int var=0; var<this->n_vars(); var++)
1870  vars_to_write.push_back(var);
1871  }
1872 
1873  const dof_id_type io_blksize = cast_int<dof_id_type>
1874  (std::min(max_io_blksize, static_cast<std::size_t>(n_objs)));
1875 
1876  const unsigned int
1877  sys_num = this->number(),
1878  num_vecs = cast_int<unsigned int>(vecs.size()),
1879  num_blks = cast_int<unsigned int>(std::ceil(static_cast<double>(n_objs)/
1880  static_cast<double>(io_blksize)));
1881 
1882  // libMesh::out << "io_blksize = " << io_blksize
1883  // << ", num_objects = " << n_objs
1884  // << ", num_blks = " << num_blks
1885  // << std::endl;
1886 
1887  dof_id_type written_length=0; // The numer of values written. This will be returned
1888  std::vector<std::vector<dof_id_type>> xfer_ids(num_blks); // The global IDs and # of components for the local objects in all blocks
1889  std::vector<std::vector<Number>> send_vals(num_blks); // The raw values for the local objects in all blocks
1890  std::vector<Parallel::Request>
1891  id_requests(num_blks), val_requests(num_blks); // send request handle for each block
1892 
1893  // ------------------------------------------------------
1894  // First pass - count the number of objects in each block
1895  // traverse all the objects and figure out which block they
1896  // will ultimately live in.
1897  std::vector<unsigned int>
1898  xfer_ids_size (num_blks,0),
1899  send_vals_size (num_blks,0);
1900 
1901  for (iterator_type it=begin; it!=end; ++it)
1902  {
1903  const dof_id_type
1904  id = (*it)->id(),
1905  block = id/io_blksize;
1906 
1907  libmesh_assert_less (block, num_blks);
1908 
1909  xfer_ids_size[block] += 2; // for each object, we store its id, as well as the total number of components for all variables
1910 
1911  unsigned int n_comp_tot=0;
1912 
1913  for (std::vector<unsigned int>::const_iterator var_it=vars_to_write.begin();
1914  var_it!=vars_to_write.end(); ++var_it)
1915  n_comp_tot += (*it)->n_comp(sys_num, *var_it); // for each variable, we will store the nonzero components
1916 
1917  send_vals_size[block] += n_comp_tot*num_vecs;
1918  }
1919 
1920  //-----------------------------------------
1921  // Collect the values for all local objects,
1922  // binning them into 'blocks' that will be
1923  // sent to processor 0
1924  for (unsigned int blk=0; blk<num_blks; blk++)
1925  {
1926  // libMesh::out << "Writing object block " << blk << std::endl;
1927 
1928  // Each processor should build up its transfer buffers for its
1929  // local objects in [first_object,last_object).
1930  const dof_id_type
1931  first_object = blk*io_blksize,
1932  last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs);
1933 
1934  // convenience
1935  std::vector<dof_id_type> & ids (xfer_ids[blk]);
1936  std::vector<Number> & vals (send_vals[blk]);
1937 
1938  // we now know the number of values we will store for each block,
1939  // so we can do efficient preallocation
1940  ids.clear(); ids.reserve (xfer_ids_size[blk]);
1941  vals.clear(); vals.reserve (send_vals_size[blk]);
1942 
1943  if (send_vals_size[blk] != 0) // only send if we have nonzero components to write
1944  for (iterator_type it=begin; it!=end; ++it)
1945  if (((*it)->id() >= first_object) && // object in [first_object,last_object)
1946  ((*it)->id() < last_object))
1947  {
1948  ids.push_back((*it)->id());
1949 
1950  // count the total number of nonzeros transferred for this object
1951  {
1952  unsigned int n_comp_tot=0;
1953 
1954  for (std::vector<unsigned int>::const_iterator var_it=vars_to_write.begin();
1955  var_it!=vars_to_write.end(); ++var_it)
1956  n_comp_tot += (*it)->n_comp(sys_num,*var_it);
1957 
1958  ids.push_back (n_comp_tot*num_vecs); // even if 0 - processor 0 has no way of knowing otherwise...
1959  }
1960 
1961  // pack the values to send
1962  for (vec_iterator_type vec_it=vecs.begin(); vec_it!=vecs.end(); ++vec_it)
1963  {
1964  const NumericVector<Number> & vec(**vec_it);
1965 
1966  for (std::vector<unsigned int>::const_iterator var_it=vars_to_write.begin();
1967  var_it!=vars_to_write.end(); ++var_it)
1968  {
1969  const unsigned int n_comp = (*it)->n_comp(sys_num,*var_it);
1970 
1971  for (unsigned int comp=0; comp<n_comp; comp++)
1972  {
1973  libmesh_assert_greater_equal ((*it)->dof_number(sys_num, *var_it, comp), vec.first_local_index());
1974  libmesh_assert_less ((*it)->dof_number(sys_num, *var_it, comp), vec.last_local_index());
1975  vals.push_back(vec((*it)->dof_number(sys_num, *var_it, comp)));
1976  }
1977  }
1978  }
1979  }
1980 
1981 #ifdef LIBMESH_HAVE_MPI
1982  Parallel::MessageTag id_tag = this->comm().get_unique_tag(100*num_blks + blk);
1983  Parallel::MessageTag val_tag = this->comm().get_unique_tag(200*num_blks + blk);
1984 
1985  // nonblocking send the data for this block
1986  this->comm().send (0, ids, id_requests[blk], id_tag);
1987  this->comm().send (0, vals, val_requests[blk], val_tag);
1988 #endif
1989  }
1990 
1991 
1992  if (this->processor_id() == 0)
1993  {
1994  std::vector<std::vector<dof_id_type>> recv_ids (this->n_processors());
1995  std::vector<std::vector<Number>> recv_vals (this->n_processors());
1996  std::vector<unsigned int> obj_val_offsets; // map to traverse entry-wise rather than processor-wise
1997  std::vector<Number> output_vals; // The output buffer for the current block
1998 
1999  // a ThreadedIO object to perform asynchronous file IO
2000  ThreadedIO<Number> threaded_io(io, output_vals);
2001  UniquePtr<Threads::Thread> async_io;
2002 
2003  for (unsigned int blk=0; blk<num_blks; blk++)
2004  {
2005  // Each processor should build up its transfer buffers for its
2006  // local objects in [first_object,last_object).
2007  const dof_id_type
2008  first_object = cast_int<dof_id_type>(blk*io_blksize),
2009  last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs),
2010  n_objects_blk = last_object - first_object;
2011 
2012  // offset array. this will define where each object's values
2013  // map into the actual output_vals buffer. this must get
2014  // 0-initialized because 0-component objects are not actually sent
2015  obj_val_offsets.resize (n_objects_blk); std::fill (obj_val_offsets.begin(), obj_val_offsets.end(), 0);
2016 
2017  std::size_t n_val_recvd_blk=0;
2018 
2019  // tags to select data received
2020  Parallel::MessageTag id_tag (this->comm().get_unique_tag(100*num_blks + blk));
2021  Parallel::MessageTag val_tag (this->comm().get_unique_tag(200*num_blks + blk));
2022 
2023  // receive this block of data from all processors.
2024  for (unsigned int comm_step=0; comm_step<this->n_processors(); comm_step++)
2025  {
2026 #ifdef LIBMESH_HAVE_MPI
2027  // blocking receive indices for this block, imposing no particular order on processor
2028  Parallel::Status id_status (this->comm().probe (Parallel::any_source, id_tag));
2029  std::vector<dof_id_type> & ids (recv_ids[id_status.source()]);
2030  this->comm().receive (id_status.source(), ids, id_tag);
2031 #else
2032  std::vector<dof_id_type> & ids (recv_ids[0]);
2033  ids = xfer_ids[blk];
2034 #endif
2035 
2036  // note its possible we didn't receive values for objects in
2037  // this block if they have no components allocated.
2038  for (std::size_t idx=0; idx<ids.size(); idx+=2)
2039  {
2040  const dof_id_type
2041  local_idx = ids[idx+0]-first_object,
2042  n_vals_tot_allvecs = ids[idx+1];
2043 
2044  libmesh_assert_less (local_idx, n_objects_blk);
2045  libmesh_assert_less (local_idx, obj_val_offsets.size());
2046 
2047  obj_val_offsets[local_idx] = n_vals_tot_allvecs;
2048  }
2049 
2050 #ifdef LIBMESH_HAVE_MPI
2051  // blocking receive values for this block, imposing no particular order on processor
2052  Parallel::Status val_status (this->comm().probe (Parallel::any_source, val_tag));
2053  std::vector<Number> & vals (recv_vals[val_status.source()]);
2054  this->comm().receive (val_status.source(), vals, val_tag);
2055 #else
2056  // straight copy without MPI
2057  std::vector<Number> & vals (recv_vals[0]);
2058  vals = send_vals[blk];
2059 #endif
2060 
2061  n_val_recvd_blk += vals.size();
2062  }
2063 
2064  // We need the offsets into the output_vals vector for each object.
2065  // fortunately, this is simply the partial sum of the total number
2066  // of components for each object
2067  std::partial_sum(obj_val_offsets.begin(), obj_val_offsets.end(),
2068  obj_val_offsets.begin());
2069 
2070  // wait on any previous asynchronous IO - this *must* complete before
2071  // we start messing with the output_vals buffer!
2072  if (async_io.get()) async_io->join();
2073 
2074  // this is the actual output buffer that will be written to disk.
2075  // at ths point we finally know wha size it will be.
2076  output_vals.resize(n_val_recvd_blk);
2077 
2078  // pack data from all processors into output values
2079  for (unsigned int proc=0; proc<this->n_processors(); proc++)
2080  {
2081  const std::vector<dof_id_type> & ids (recv_ids [proc]);
2082  const std::vector<Number> & vals(recv_vals[proc]);
2083  std::vector<Number>::const_iterator proc_vals(vals.begin());
2084 
2085  for (std::size_t idx=0; idx<ids.size(); idx+=2)
2086  {
2087  const dof_id_type
2088  local_idx = ids[idx+0]-first_object,
2089  n_vals_tot_allvecs = ids[idx+1];
2090 
2091  // put this object's data into the proper location
2092  // in the output buffer
2093  std::vector<Number>::iterator out_vals(output_vals.begin());
2094  if (local_idx != 0)
2095  std::advance (out_vals, obj_val_offsets[local_idx-1]);
2096 
2097  for (unsigned int val=0; val<n_vals_tot_allvecs; val++, ++out_vals, ++proc_vals)
2098  {
2099  libmesh_assert (out_vals != output_vals.end());
2100  libmesh_assert (proc_vals != vals.end());
2101  *out_vals = *proc_vals;
2102  }
2103  }
2104  }
2105 
2106  // output_vals buffer is now filled for this block.
2107  // write it to disk
2108  async_io.reset(new Threads::Thread(threaded_io));
2109  written_length += output_vals.size();
2110  }
2111 
2112  // wait on any previous asynchronous IO - this *must* complete before
2113  // our stuff goes out of scope
2114  async_io->join();
2115  }
2116 
2117  Parallel::wait(id_requests);
2118  Parallel::wait(val_requests);
2119 
2120  // we need some synchronization here. Because this method
2121  // can be called for a range of nodes, then a range of elements,
2122  // we need some mechanism to prevent processors from racing past
2123  // to the next range and overtaking ongoing communication. one
2124  // approach would be to figure out unique tags for each range,
2125  // but for now we just impose a barrier here. And might as
2126  // well have it do some useful work.
2127  this->comm().broadcast(written_length);
2128 
2129  return written_length;
2130 }
2131 
2132 
2133 
2134 unsigned int System::write_SCALAR_dofs (const NumericVector<Number> & vec,
2135  const unsigned int var,
2136  Xdr & io) const
2137 {
2138  unsigned int written_length=0;
2139  std::vector<Number> vals; // The raw values for the local objects in the current block
2140  // Collect the SCALARs for the current variable
2141  if (this->processor_id() == (this->n_processors()-1))
2142  {
2143  const DofMap & dof_map = this->get_dof_map();
2144  std::vector<dof_id_type> SCALAR_dofs;
2145  dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
2146  const unsigned int n_scalar_dofs = cast_int<unsigned int>
2147  (SCALAR_dofs.size());
2148 
2149  for (unsigned int i=0; i<n_scalar_dofs; i++)
2150  {
2151  vals.push_back( vec(SCALAR_dofs[i]) );
2152  }
2153  }
2154 
2155 #ifdef LIBMESH_HAVE_MPI
2156  if (this->n_processors() > 1)
2157  {
2158  const Parallel::MessageTag val_tag(1);
2159 
2160  // Post the receive on processor 0
2161  if (this->processor_id() == 0)
2162  {
2163  this->comm().receive(this->n_processors()-1, vals, val_tag);
2164  }
2165 
2166  // Send the data to processor 0
2167  if (this->processor_id() == (this->n_processors()-1))
2168  {
2169  this->comm().send(0, vals, val_tag);
2170  }
2171  }
2172 #endif
2173 
2174  // -------------------------------------------------------
2175  // Write the output on processor 0.
2176  if (this->processor_id() == 0)
2177  {
2178  const unsigned int vals_size =
2179  cast_int<unsigned int>(vals.size());
2180  io.data_stream (&vals[0], vals_size);
2181  written_length += vals_size;
2182  }
2183 
2184  return written_length;
2185 }
2186 
2187 
2188 
2189 dof_id_type System::write_serialized_vector (Xdr & io,
2190  const NumericVector<Number> & vec) const
2191 {
2192  parallel_object_only();
2193 
2194  libmesh_assert (io.writing());
2195 
2196  dof_id_type vec_length = vec.size();
2197  if (this->processor_id() == 0) io.data (vec_length, "# vector length");
2198 
2199  dof_id_type written_length = 0;
2200 
2201  //---------------------------------
2202  // Collect the values for all nodes
2203  written_length += cast_int<dof_id_type>
2204  (this->write_serialized_blocked_dof_objects (std::vector<const NumericVector<Number> *>(1,&vec),
2205  this->get_mesh().n_nodes(),
2206  this->get_mesh().local_nodes_begin(),
2207  this->get_mesh().local_nodes_end(),
2208  io));
2209 
2210  //------------------------------------
2211  // Collect the values for all elements
2212  written_length += cast_int<dof_id_type>
2213  (this->write_serialized_blocked_dof_objects (std::vector<const NumericVector<Number> *>(1,&vec),
2214  this->get_mesh().n_elem(),
2215  this->get_mesh().local_elements_begin(),
2216  this->get_mesh().local_elements_end(),
2217  io));
2218 
2219  //-------------------------------------------
2220  // Finally loop over all the SCALAR variables
2221  for (unsigned int var=0; var<this->n_vars(); var++)
2222  if (this->variable(var).type().family == SCALAR)
2223  {
2224  written_length +=
2225  this->write_SCALAR_dofs (vec, var, io);
2226  }
2227 
2228  if (this->processor_id() == 0)
2229  libmesh_assert_equal_to (written_length, vec_length);
2230 
2231  return written_length;
2232 }
2233 
2234 
2235 template <typename InValType>
2236 std::size_t System::read_serialized_vectors (Xdr & io,
2237  const std::vector<NumericVector<Number> *> & vectors) const
2238 {
2239  parallel_object_only();
2240 
2241  // Error checking
2242  // #ifndef NDEBUG
2243  // // In parallel we better be reading a parallel vector -- if not
2244  // // we will not set all of its components below!!
2245  // if (this->n_processors() > 1)
2246  // {
2247  // libmesh_assert (vec.type() == PARALLEL ||
2248  // vec.type() == GHOSTED);
2249  // }
2250  // #endif
2251 
2252  libmesh_assert (io.reading());
2253 
2254  if (this->processor_id() == 0)
2255  {
2256  // sizes
2257  unsigned int num_vecs=0;
2258  dof_id_type vector_length=0;
2259 
2260  // Get the number of vectors
2261  io.data(num_vecs);
2262  // Get the buffer size
2263  io.data(vector_length);
2264 
2265  libmesh_assert_equal_to (num_vecs, vectors.size());
2266 
2267  if (num_vecs != 0)
2268  {
2269  libmesh_assert_not_equal_to (vectors[0], 0);
2270  libmesh_assert_equal_to (vectors[0]->size(), vector_length);
2271  }
2272  }
2273 
2274  // no need to actually communicate these.
2275  // this->comm().broadcast(num_vecs);
2276  // this->comm().broadcast(vector_length);
2277 
2278  // Cache these - they are not free!
2279  const dof_id_type
2280  n_nodes = this->get_mesh().n_nodes(),
2281  n_elem = this->get_mesh().n_elem();
2282 
2283  std::size_t read_length = 0;
2284 
2285  //---------------------------------
2286  // Collect the values for all nodes
2287  read_length +=
2288  this->read_serialized_blocked_dof_objects (n_nodes,
2289  this->get_mesh().local_nodes_begin(),
2290  this->get_mesh().local_nodes_end(),
2291  InValType(),
2292  io,
2293  vectors);
2294 
2295  //------------------------------------
2296  // Collect the values for all elements
2297  read_length +=
2298  this->read_serialized_blocked_dof_objects (n_elem,
2299  this->get_mesh().local_elements_begin(),
2300  this->get_mesh().local_elements_end(),
2301  InValType(),
2302  io,
2303  vectors);
2304 
2305  //-------------------------------------------
2306  // Finally loop over all the SCALAR variables
2307  for (std::size_t vec=0; vec<vectors.size(); vec++)
2308  for (unsigned int var=0; var<this->n_vars(); var++)
2309  if (this->variable(var).type().family == SCALAR)
2310  {
2311  libmesh_assert_not_equal_to (vectors[vec], 0);
2312 
2313  read_length +=
2314  this->read_SCALAR_dofs (var, io, vectors[vec]);
2315  }
2316 
2317  //---------------------------------------
2318  // last step - must close all the vectors
2319  for (std::size_t vec=0; vec<vectors.size(); vec++)
2320  {
2321  libmesh_assert_not_equal_to (vectors[vec], 0);
2322  vectors[vec]->close();
2323  }
2324 
2325  return read_length;
2326 }
2327 
2328 
2329 
2330 std::size_t System::write_serialized_vectors (Xdr & io,
2331  const std::vector<const NumericVector<Number> *> & vectors) const
2332 {
2333  parallel_object_only();
2334 
2335  libmesh_assert (io.writing());
2336 
2337  // Cache these - they are not free!
2338  const dof_id_type
2339  n_nodes = this->get_mesh().n_nodes(),
2340  n_elem = this->get_mesh().n_elem();
2341 
2342  std::size_t written_length = 0;
2343 
2344  if (this->processor_id() == 0)
2345  {
2346  unsigned int
2347  n_vec = cast_int<unsigned int>(vectors.size());
2348  dof_id_type
2349  vec_size = vectors.empty() ? 0 : vectors[0]->size();
2350  // Set the number of vectors
2351  io.data(n_vec, "# number of vectors");
2352  // Set the buffer size
2353  io.data(vec_size, "# vector length");
2354  }
2355 
2356  //---------------------------------
2357  // Collect the values for all nodes
2358  written_length +=
2359  this->write_serialized_blocked_dof_objects (vectors,
2360  n_nodes,
2361  this->get_mesh().local_nodes_begin(),
2362  this->get_mesh().local_nodes_end(),
2363  io);
2364 
2365  //------------------------------------
2366  // Collect the values for all elements
2367  written_length +=
2368  this->write_serialized_blocked_dof_objects (vectors,
2369  n_elem,
2370  this->get_mesh().local_elements_begin(),
2371  this->get_mesh().local_elements_end(),
2372  io);
2373 
2374  //-------------------------------------------
2375  // Finally loop over all the SCALAR variables
2376  for (std::size_t vec=0; vec<vectors.size(); vec++)
2377  for (unsigned int var=0; var<this->n_vars(); var++)
2378  if (this->variable(var).type().family == SCALAR)
2379  {
2380  libmesh_assert_not_equal_to (vectors[vec], 0);
2381 
2382  written_length +=
2383  this->write_SCALAR_dofs (*vectors[vec], var, io);
2384  }
2385 
2386  return written_length;
2387 }
2388 
2389 
2390 
2391 
2392 template void System::read_parallel_data<Number> (Xdr & io, const bool read_additional_data);
2393 template void System::read_serialized_data<Number> (Xdr & io, const bool read_additional_data);
2394 template numeric_index_type System::read_serialized_vector<Number> (Xdr & io, NumericVector<Number> * vec);
2395 template std::size_t System::read_serialized_vectors<Number> (Xdr & io, const std::vector<NumericVector<Number> *> & vectors) const;
2396 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2397 template void System::read_parallel_data<Real> (Xdr & io, const bool read_additional_data);
2398 template void System::read_serialized_data<Real> (Xdr & io, const bool read_additional_data);
2399 template numeric_index_type System::read_serialized_vector<Real> (Xdr & io, NumericVector<Number> * vec);
2400 template std::size_t System::read_serialized_vectors<Real> (Xdr & io, const std::vector<NumericVector<Number> *> & vectors) const;
2401 #endif
2402 
2403 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
Simple compatibility class for std::thread &#39;concurrent&#39; execution.
Definition: threads.h:69
void data(T &a, const char *comment="")
Inputs or outputs a single value.
Definition: xdr_cxx.C:761
FEFamily family
The type of finite element.
Definition: fe_type.h:203
int get_order() const
Explicitly request the order as an int.
Definition: fe_type.h:77
Status wait(Request &r)
Wait for a non-blocking send or receive to finish.
Definition: parallel.h:565
virtual numeric_index_type size() const =0
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
void join()
Join is a no-op, since the constructor blocked until completion.
Definition: threads.h:83
virtual numeric_index_type last_local_index() const =0
const unsigned int any_source
Accept from any source.
Definition: parallel.h:204
Encapsulates the MPI_Status struct.
Definition: parallel.h:461
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
ImplicitSystem & sys
void comment(std::string &)
Writes or reads (ignores) a comment line.
Definition: xdr_cxx.C:1519
OrderWrapper radial_order
The approximation order in the base of the infinite element.
Definition: fe_type.h:236
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
InfMapType
defines an enum for the types of coordinate mappings available in infinite elements.
libmesh_assert(j)
Tnew cast_int(Told oldvar)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
FEFamily
defines an enum for finite element families.
const dof_id_type n_nodes
Definition: tecplot_io.C:67
dof_id_type numeric_index_type
Definition: id_types.h:92
Encapsulates the MPI tag integers.
Definition: parallel.h:227
bool reading() const
Definition: xdr_cxx.h:113
int version() const
Gets the version of the file that is being read.
Definition: xdr_cxx.h:166
bool writing() const
Definition: xdr_cxx.h:119
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition: fe_type.h:257
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:68
FEFamily radial_family
For InfFE, family contains the radial shape family, while base_family contains the approximation type...
Definition: fe_type.h:249
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Fills the vector di with the global degree of freedom indices corresponding to the SCALAR variable vn...
Definition: dof_map.C:2288
bool is_open() const
Definition: xdr_cxx.C:339
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
OStreamProxy out
ParallelType type() const
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
Definition: dof_object.h:51
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
virtual numeric_index_type first_local_index() const =0
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
IterBase * data
Ideally this private member data should have protected access.
void data_stream(T *val, const unsigned int len, const unsigned int line_break=libMesh::invalid_uint)
Inputs or outputs a raw data stream.
Definition: xdr_cxx.C:825
processor_id_type processor_id()
Definition: libmesh_base.h:96
long double min(long double a, double b)
sys get_dof_map()
status probe(const unsigned int src_processor_id, const MessageTag &tag=any_tag, const Communicator &comm=Communicator_World)
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type n_processors()
Definition: libmesh_base.h:88