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

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

Hosted By:
SourceForge.net Logo