legacy_xdr_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 
20 // C++ includes
21 #include <iostream>
22 #include <iomanip>
23 
24 #include <vector>
25 #include <string>
26 
27 // Local includes
28 #include "libmesh/legacy_xdr_io.h"
29 #include "libmesh/mesh_base.h"
30 #include "libmesh/mesh_data.h"
31 #include "libmesh/mesh_tools.h" // MeshTools::n_levels(mesh)
32 #include "libmesh/parallel.h" // - which makes write parallel-only
33 #include "libmesh/cell_hex27.h" // Needed for MGF-style Hex27 meshes
34 #include "libmesh/boundary_info.h"
36 #include "libmesh/xdr_mgf.h"
37 #include "libmesh/xdr_mesh.h"
38 #include "libmesh/xdr_mhead.h"
39 #include "libmesh/xdr_soln.h"
40 #include "libmesh/xdr_shead.h"
41 
42 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
43 #include "libmesh/utility.h"
44 #endif
45 
46 
47 namespace libMesh
48 {
49 
50 
51 
52 
53 
54 
55 // ------------------------------------------------------------
56 // LegacyXdrIO members
57 LegacyXdrIO::LegacyXdrIO (MeshBase& mesh, const bool binary_in) :
58  MeshInput<MeshBase> (mesh),
59  MeshOutput<MeshBase>(mesh),
60  _binary (binary_in)
61 {
62 }
63 
64 
65 
66 
67 LegacyXdrIO::LegacyXdrIO (const MeshBase& mesh, const bool binary_in) :
68  MeshOutput<MeshBase>(mesh),
69  _binary (binary_in)
70 {
71 }
72 
73 
74 
75 
77 {
78 }
79 
80 
81 
82 
84 {
85  return _binary;
86 }
87 
88 
89 
90 
91 bool LegacyXdrIO::binary () const
92 {
93  return _binary;
94 }
95 
96 
97 
98 void LegacyXdrIO::read (const std::string& name)
99 {
100  // This is a serial-only process for now;
101  // the Mesh should be read on processor 0 and
102  // broadcast later
104  return;
105 
106  if (this->binary())
107  this->read_binary (name);
108  else
109  this->read_ascii (name);
110 }
111 
112 
113 
114 void LegacyXdrIO::read_mgf (const std::string& name)
115 {
116  if (this->binary())
117  this->read_binary (name, LegacyXdrIO::MGF);
118  else
119  this->read_ascii (name, LegacyXdrIO::MGF);
120 }
121 
122 
123 
124 void LegacyXdrIO::write (const std::string& name)
125 {
126  if (this->binary())
127  this->write_binary (name);
128  else
129  this->write_ascii (name);
130 }
131 
132 
133 
134 void LegacyXdrIO::write_mgf (const std::string& name)
135 {
136  if (this->binary())
137  this->write_binary (name, LegacyXdrIO::MGF);
138  else
139  this->write_ascii (name, LegacyXdrIO::MGF);
140 }
141 
142 
143 
144 void LegacyXdrIO::read_mgf_soln (const std::string& name,
145  std::vector<Number>& soln,
146  std::vector<std::string>& var_names) const
147 {
148  libmesh_deprecated();
149 
150 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
151 
152  // buffer for writing separately
153  std::vector<Real> real_soln;
154  std::vector<Real> imag_soln;
155 
156  Utility::prepare_complex_data (soln, real_soln, imag_soln);
157 
158  this->read_soln (Utility::complex_filename(name, 0),
159  real_soln,
160  var_names);
161 
162  this->read_soln (Utility::complex_filename(name, 1),
163  imag_soln,
164  var_names);
165 
166 #else
167 
168  this->read_soln (name, soln, var_names);
169 
170 #endif
171 }
172 
173 
174 
175 void LegacyXdrIO::write_mgf_soln (const std::string& name,
176  std::vector<Number>& soln,
177  std::vector<std::string>& var_names) const
178 {
179  libmesh_deprecated();
180 
181 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
182 
183  // buffer for writing separately
184  std::vector<Real> real_soln;
185  std::vector<Real> imag_soln;
186 
187  Utility::prepare_complex_data (soln, real_soln, imag_soln);
188 
189  this->write_soln (Utility::complex_filename(name, 0),
190  real_soln,
191  var_names);
192 
193  this->write_soln (Utility::complex_filename(name, 1),
194  imag_soln,
195  var_names);
196 
197 #else
198 
199  this->write_soln (name, soln, var_names);
200 
201 #endif
202 }
203 
204 
205 
206 void LegacyXdrIO::read_ascii (const std::string& name, const LegacyXdrIO::FileFormat originator)
207 {
208  // get a writeable reference to the underlying mesh
210 
211  // clear any existing mesh data
212  mesh.clear();
213 
214  // read the mesh
215  this->read_mesh (name, originator);
216 }
217 
218 
219 
220 #ifndef LIBMESH_HAVE_XDR
221 void LegacyXdrIO::read_binary (const std::string& name, const LegacyXdrIO::FileFormat)
222 {
223 
224  libMesh::err << "WARNING: Compiled without XDR binary support.\n"
225  << "Will try ASCII instead" << std::endl << std::endl;
226 
227  this->read_ascii (name);
228 }
229 #else
230 void LegacyXdrIO::read_binary (const std::string& name, const LegacyXdrIO::FileFormat originator)
231 {
232  // get a writeable reference to the underlying mesh
234 
235  // clear any existing mesh data
236  mesh.clear();
237 
238  // read the mesh
239  this->read_mesh (name, originator);
240 }
241 #endif
242 
243 
244 
245 void LegacyXdrIO::write_ascii (const std::string& name, const LegacyXdrIO::FileFormat originator)
246 {
247  this->write_mesh (name, originator);
248 }
249 
250 
251 
252 #ifndef LIBMESH_HAVE_XDR
253 void LegacyXdrIO::write_binary (const std::string& name, const LegacyXdrIO::FileFormat)
254 {
255  libMesh::err << "WARNING: Compiled without XDR binary support.\n"
256  << "Will try ASCII instead" << std::endl << std::endl;
257 
258  this->write_ascii (name);
259 }
260 #else
261 void LegacyXdrIO::write_binary (const std::string& name, const LegacyXdrIO::FileFormat originator)
262 {
263  this->write_mesh (name, originator);
264 }
265 #endif
266 
267 
268 
269 void LegacyXdrIO::read_mesh (const std::string& name,
270  const LegacyXdrIO::FileFormat originator,
271  MeshData* mesh_data)
272 {
273  // This is a serial-only process for now;
274  // the Mesh should be read on processor 0 and
275  // broadcast later
276  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
277 
278  // get a writeable reference to the mesh
280 
281  // clear any data in the mesh
282  mesh.clear();
283 
284  // Keep track of what kinds of elements this file contains
285  elems_of_dimension.clear();
286  elems_of_dimension.resize(4, false);
287 
288  // Create an XdrMESH object.
289  XdrMESH m;
290 
291  // Create a pointer
292  // to an XdrMESH file
293  // header.
294  XdrMHEAD mh;
295 
296  // Open the XDR file for reading.
297  // Note 1: Provide an additional argument
298  // to specify the dimension.
299  //
300  // Note 2: Has to do the right thing for
301  // both binary and ASCII files.
302  m.set_orig_flag(originator);
303  m.init((this->binary() ? XdrMGF::DECODE : XdrMGF::R_ASCII), name.c_str(), 0); // mesh files are always number 0 ...
304 
305  // From here on, things depend
306  // on whether we are reading or
307  // writing! First, we define
308  // header variables that may
309  // be read OR written.
310  unsigned int n_blocks = 0;
311  unsigned int n_levels = 0;
312 
313  if (m.get_orig_flag() == LegacyXdrIO::LIBM)
314  n_levels = m.get_num_levels();
315 
316 
317  std::vector<ElemType> etypes;
318  std::vector<unsigned int> neeb;
319 
320  // Get the information from
321  // the header, and place it
322  // in the header pointer.
323  m.header(&mh);
324 
325  // Read information from the
326  // file header. This depends on
327  // whether its a libMesh or MGF mesh.
328  const int numElem = mh.getNumEl();
329  const int numNodes = mh.getNumNodes();
330  const int totalWeight = mh.getSumWghts();
331  const int numBCs = mh.getNumBCs();
332 
333  // If a libMesh-type mesh, read the augmented mesh information
335  {
336  // Read augmented header
337  n_blocks = mh.get_n_blocks();
338 
339  etypes.resize(n_blocks);
340  mh.get_block_elt_types(etypes);
341 
342  mh.get_num_elem_each_block(neeb);
343  }
344 
345 
346 
347  // Read the connectivity
348  std::vector<int> conn;
349 
350  // Now that we know the
351  // number of nodes and elements,
352  // we can resize the
353  // appropriate vectors if we are
354  // reading information in.
355  mesh.reserve_nodes (numNodes);
356  mesh.reserve_elem (numElem);
357 
358  // Each element stores two extra
359  // locations: one which tells
360  // what type of element it is,
361  // and one which tells how many
362  // nodes it has. Therefore,
363  // the total number of nodes
364  // (totalWeight) must be augmented
365  // by 2 times the number of elements
366  // in order to read in the entire
367  // connectivity array.
368 
369  // Note: This section now depends on
370  // whether we are reading an old-style libMesh,
371  // MGF, or a new-style libMesh mesh.
372  if (m.get_orig_flag() == LegacyXdrIO::DEAL)
373  {
374  conn.resize(totalWeight);
375  m.Icon(&conn[0], 1, totalWeight);
376  }
377 
378  else if (m.get_orig_flag() == LegacyXdrIO::MGF)
379  {
380  conn.resize(totalWeight+(2*numElem));
381  m.Icon(&conn[0], 1, totalWeight+(2*numElem));
382  }
383 
384  else if (m.get_orig_flag() == LegacyXdrIO::LIBM)
385  {
386  conn.resize(totalWeight);
387  m.Icon(&conn[0], 1, totalWeight);
388  }
389 
390  else
391  {
392  // I don't know what type of mesh it is.
393  libmesh_error();
394  }
395 
396  // read in the nodal coordinates and form points.
397  {
398  std::vector<Real> coords(numNodes*3); // Always use three coords per node
399  m.coord(&coords[0], 3, numNodes);
400 
401 
402 
403  // Form Nodes out of
404  // the coordinates. If the
405  // MeshData object is active,
406  // add the nodes and ids also
407  // to its map.
408  for (int innd=0; innd<numNodes; ++innd)
409  {
410  Node* node = mesh.add_point (Point(coords[0+innd*3],
411  coords[1+innd*3],
412  coords[2+innd*3]), innd);
413 
414  if (mesh_data != NULL)
415  if (mesh_data->active())
416  {
417  // add the id to the MeshData, so that
418  // it knows the foreign id, even when
419  // the underlying mesh got re-numbered,
420  // refined, elements/nodes added...
421  mesh_data->add_foreign_node_id(node, innd);
422  }
423  }
424  }
425 
426 
427 
428  // Build the elements.
429  // Note: If the originator was MGF, we don't
430  // have to do much checking ...
431  // all the elements are Hex27.
432  // If the originator was
433  // this code, we have to loop over
434  // et and neeb to read in all the
435  // elements correctly.
436  //
437  // (This used to be before the coords block, but it
438  // had to change now that elements store pointers to
439  // nodes. The nodes must exist before we assign them to
440  // the elements. BSK, 1/13/2003)
442  {
443  unsigned int lastConnIndex = 0;
444  unsigned int lastFaceIndex = 0;
445 
446  // This map keeps track of elements we've previously
447  // constructed, to avoid O(n) lookup times for parent pointers
448  // and to enable elements to be added in ascending ID order
449  std::map<unsigned int, Elem*> parents;
450 
451  {
452  // Keep track of Element ids in MGF-style meshes;
453  unsigned int next_elem_id = 0;
454 
455  for (unsigned int level=0; level<=n_levels; level++)
456  {
457  for (unsigned int idx=0; idx<n_blocks; idx++)
458  {
459  for (unsigned int e=lastFaceIndex; e<lastFaceIndex+neeb[level*n_blocks+idx]; e++)
460  {
461  // Build a temporary element of the right type, so we know how
462  // connectivity entries will be on the line for this element.
463  AutoPtr<Elem> temp_elem = Elem::build(etypes[idx]);
464 
465  // A pointer to the element which will eventually be added to the mesh.
466  Elem* elem;
467 
468  // New-style libMesh mesh
469  if (m.get_orig_flag() == LegacyXdrIO::LIBM)
470  {
471  unsigned int self_ID = conn[lastConnIndex + temp_elem->n_nodes()];
472 
473 #ifdef LIBMESH_ENABLE_AMR
474  unsigned int parent_ID = conn[lastConnIndex + temp_elem->n_nodes()+1];
475 
476  if (level > 0)
477  {
478  // Do a linear search for the parent
479  Elem* my_parent;
480 
481  // Search for parent in the parents map (log(n))
482  START_LOG("log(n) search for parent", "LegacyXdrIO::read_mesh");
483  std::map<unsigned int, Elem*>::iterator it = parents.find(parent_ID);
484  STOP_LOG("log(n) search for parent", "LegacyXdrIO::read_mesh");
485 
486  // If the parent was not previously added, we cannot continue.
487  if (it == parents.end())
488  {
489  libMesh::err << "Parent element with ID " << parent_ID
490  << " not found." << std::endl;
491  libmesh_error();
492  }
493 
494  // Set the my_parent pointer
495  my_parent = (*it).second;
496 
497  // my_parent is now INACTIVE, since he has children
499 
500  // Now that we know the parent, build the child
501  elem = Elem::build(etypes[idx],my_parent).release();
502 
503  // The new child is marked as JUST_REFINED
505 
506  // Tell the parent about his new child
507  my_parent->add_child(elem);
508 
509  // sanity check
510  libmesh_assert_equal_to (my_parent->type(), elem->type());
511  }
512 
513  // Add level-0 elements to the mesh
514  else
515 #endif // #ifdef LIBMESH_ENABLE_AMR
516  {
517  elem = Elem::build(etypes[idx]).release();
518  }
519 
520  // Assign the newly-added element's ID so that future
521  // children which may be added can find it correctly.
522  elem->set_id() = self_ID;
523 
524  // Add this element to the map, it may be a parent for a future element
525  START_LOG("insert elem into map", "LegacyXdrIO::read_mesh");
526  parents[self_ID] = elem;
527  STOP_LOG("insert elem into map", "LegacyXdrIO::read_mesh");
528  }
529 
530  // MGF-style meshes
531  else
532  {
533  elem = Elem::build(etypes[idx]).release();
534  elem->set_id(next_elem_id++);
535 
536  elems_of_dimension[elem->dim()] = true;
537 
538  mesh.add_elem(elem);
539  }
540 
541  // Add elements with the same id as in libMesh.
542  // Provided the data files that MeshData reads
543  // were only written with MeshData, then this
544  // should work properly. This is an inline
545  // function, so that for disabled MeshData, this
546  // should not induce too much cost
547  if (mesh_data != NULL)
548  mesh_data->add_foreign_elem_id (elem, e);
549 
550  // Set the node pointers of the newly-created element
551  for (unsigned int innd=0; innd < elem->n_nodes(); innd++)
552  {
553  elem->set_node(innd) = mesh.node_ptr(conn[innd+lastConnIndex]);
554  }
555 
556  lastConnIndex += (m.get_orig_flag() == LegacyXdrIO::LIBM) ? (elem->n_nodes()+2) : elem->n_nodes();
557  }
558  lastFaceIndex += neeb[idx];
559  }
560  }
561  }
562 
563  if (m.get_orig_flag() == LegacyXdrIO::LIBM)
564  {
565  {
566  // Iterate in ascending elem ID order
567  unsigned int next_elem_id = 0;
568  for (std::map<unsigned int, Elem *>::iterator i =
569  parents.begin();
570  i != parents.end(); ++i)
571  {
572  Elem *elem = i->second;
573  if (elem)
574  {
575  elem->set_id(next_elem_id++);
576 
577  elems_of_dimension[elem->dim()] = true;
578 
579  mesh.add_elem(elem);
580  }
581  else
582  // We can probably handle this, but we don't expect it
583  libmesh_error();
584  }
585  }
586  }
587  }
588 
589  // MGF-style (1) Hex27 mesh
590  else if (m.get_orig_flag() == LegacyXdrIO::MGF)
591  {
592 
593 #ifdef DEBUG
594  if (mesh_data != NULL)
595  if (mesh_data->active())
596  {
597  libMesh::err << "ERROR: MeshData not implemented for MGF-style mesh."
598  << std::endl;
599  libmesh_error();
600  }
601 #endif
602 
603  for (int ielm=0; ielm < numElem; ++ielm)
604  {
605  Elem* elem = new Hex27;
606  elem->set_id(ielm);
607 
608  elems_of_dimension[elem->dim()] = true;
609 
610  mesh.add_elem(elem);
611 
612  for (int innd=0; innd < 27; ++innd)
613  elem->set_node(innd) = mesh.node_ptr(conn[innd+2+(27+2)*ielm]);
614  }
615  }
616 
617  // Set the mesh dimension to the largest encountered for an element
618  for (unsigned int i=0; i!=4; ++i)
619  if (elems_of_dimension[i])
620  mesh.set_mesh_dimension(i);
621 
622 #if LIBMESH_DIM < 3
623  if (mesh.mesh_dimension() > LIBMESH_DIM)
624  {
625  libMesh::err << "Cannot open dimension " <<
626  mesh.mesh_dimension() <<
627  " mesh file when configured without " <<
628  mesh.mesh_dimension() << "D support." <<
629  std::endl;
630  libmesh_error();
631  }
632 #endif
633 
634  // tell the MeshData object that we are finished
635  // reading data
636  if (mesh_data != NULL)
637  mesh_data->close_foreign_id_maps ();
638 
639  // Free memory used in
640  // the connectivity
641  // vector.
642  conn.clear();
643 
644 
645  // If we are reading,
646  // read in the BCs
647  // from the mesh file,
648  // otherwise write the
649  // boundary conditions
650  // if the BoundaryInfo
651  // object exists.
652  if (numBCs > 0)
653  {
654  std::vector<int> bcs(numBCs*3);
655 
656  // Read the BCs from the XDR file
657  m.BC(&bcs[0], numBCs);
658 
659  // Add to the boundary_info
660  for (int ibc=0; ibc < numBCs; ibc++)
661  mesh.boundary_info->add_side(bcs[0+ibc*3], bcs[1+ibc*3], bcs[2+ibc*3]);
662  }
663 }
664 
665 
666 
667 void LegacyXdrIO::write_mesh (const std::string& name,
668  const LegacyXdrIO::FileFormat originator)
669 {
670  // get a read-only reference to the mesh
671  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();
672 
673  // n_levels is a parallel-only method
674  libmesh_parallel_only(mesh.comm());
675  const unsigned int n_levels = MeshTools::n_levels(mesh);
676 
677  // The Legacy Xdr IO code only works if we have a serialized mesh
678  libmesh_assert (mesh.is_serial());
679 
680  // In which case only processor 0 needs to do any writing
681  if (MeshOutput<MeshBase>::mesh().processor_id() != 0)
682  return;
683 
684  // Create an XdrMESH object.
685  XdrMESH m;
686 
687  // Create a pointer
688  // to an XdrMESH file
689  // header.
690  XdrMHEAD mh;
691 
692  // Open the XDR file for writing.
693  // Note 1: Provide an additional argument
694  // to specify the dimension.
695  //
696  // Note 2: Has to do the right thing for
697  // both binary and ASCII files.
698  m.set_orig_flag(originator);
699 
700  // From here on, things depend
701  // on whether we are reading or
702  // writing! First, we define
703  // header variables that may
704  // be read OR written.
705  std::vector<unsigned int> neeb;
706  std::vector<ElemType> etypes;
707 
708 
709  int n_non_subactive = 0;
710  int non_subactive_weight = 0;
711 
712  // This map will associate
713  // the distance from the beginning of the set
714  // to each node ID with the node ID itself.
715  std::map<dof_id_type, dof_id_type> node_map;
716 
717  {
718  // For each non-subactive element:
719  // 1.) Increment the number of non subactive elements
720  // 2.) Accumulate the total weight
721  // 3.) Add the node ids to a set of non subactive node ids
722  std::set<dof_id_type> not_subactive_node_ids;
724  const MeshBase::const_element_iterator end_el = mesh.elements_end();
725  for( ; el != end_el; ++el)
726  {
727  const Elem* elem = (*el);
728  if(!elem->subactive())
729  {
730  n_non_subactive++;
731  non_subactive_weight += elem->n_nodes();
732 
733  for (unsigned int n=0; n<elem->n_nodes(); ++n)
734  not_subactive_node_ids.insert(elem->node(n));
735  }
736  }
737 
738  // Now that the set is built, most of the hard work is done. We build
739  // the map next and let the set go out of scope.
740  std::set<dof_id_type>::iterator it = not_subactive_node_ids.begin();
741  const std::set<dof_id_type>::iterator end = not_subactive_node_ids.end();
742  dof_id_type cnt=0;
743  for (; it!=end; ++it)
744  node_map[*it] = cnt++;
745  }
746 
747 
748  const int numElem = n_non_subactive;
749  const int numBCs = mesh.boundary_info->n_boundary_conds();
750 
751  // Fill the etypes vector with all of the element types found in the mesh
752  MeshTools::elem_types(mesh, etypes);
753 
754  // store number of elements in each block at each refinement level
755  neeb.resize((n_levels+1)*etypes.size());
756 
757  // Store a variable for the number of element types
758  const unsigned int n_el_types =
759  libmesh_cast_int<unsigned int>(etypes.size());
760 
761  m.set_num_levels(n_levels);
762 
763  // The last argument is zero because mesh files are always number 0 ...
764  m.init((this->binary() ? XdrMGF::ENCODE : XdrMGF::W_ASCII), name.c_str(), 0);
765 
766  // Loop over all levels and all element types to set the entries of neeb
767  for(unsigned int level=0; level<=n_levels; level++)
768  for (unsigned int el_type=0; el_type<n_el_types; el_type++)
769  neeb[level*n_el_types + el_type] =
770  MeshTools::n_non_subactive_elem_of_type_at_level(mesh, etypes[el_type], level);
771  // gotta change this function name!!!
772 
773 
774  // Now we check to see if we're doing
775  // MGF-style headers or libMesh-style
776  // "augmented" headers. An
777  // augmented header contains
778  // information about mesh blocks,
779  // allowing us to optimize storage
780  // and minimize IO requirements
781  // for these meshes.
783  {
784  mh.set_n_blocks(etypes.size());
785  mh.set_block_elt_types(etypes);
786  mh.set_num_elem_each_block(neeb);
787  }
788  else
789  libmesh_assert_equal_to (etypes.size(), 1);
790 
791  mh.setNumEl(numElem);
792  mh.setNumNodes(node_map.size());
793  mh.setStrSize(65536);
794 
795  // set a local variable for the total weight of the mesh
796  int totalWeight =0;
797 
798  if (m.get_orig_flag() == LegacyXdrIO::DEAL) // old-style LibMesh
799  totalWeight=MeshTools::total_weight(mesh);
800 
801  else if (m.get_orig_flag() == LegacyXdrIO::MGF) // MGF-style
802  totalWeight = MeshTools::total_weight(mesh)+2*numElem;
803 
804  else if (m.get_orig_flag() == LegacyXdrIO::LIBM) // new-style LibMesh
805  totalWeight = non_subactive_weight+2*numElem;
806 
807  else
808  libmesh_error();
809 
810  // Set the total weight in the header
811  mh.setSumWghts(totalWeight);
812 
813  mh.setNumBCs(numBCs);
814  mh.setId("Id String"); // You can put whatever you want, it will be ignored
815  mh.setTitle("Title String"); // You can put whatever you want, it will be ignored
816 
817  // Put the information
818  // in the XDR file.
819  m.header(&mh);
820 
821 
822  // Write the connectivity
823  {
824  std::vector<int> conn;
825  LegacyXdrIO::FileFormat orig_type = m.get_orig_flag();
826 
827  // Resize the connectivity vector to hold all the connectivity for the mesh
828  conn.resize(totalWeight);
829 
830  unsigned int lastConnIndex = 0;
831  unsigned int nn = 0;
832 
833  // Loop over levels and types again, write connectivity information to conn.
834  for (unsigned int level=0; level<=n_levels; level++)
835  for (unsigned int idx=0; idx<etypes.size(); idx++)
836  {
837  nn = lastConnIndex = 0;
838 
839  for (unsigned int e=0; e<mesh.n_elem(); e++)
840  if ((mesh.elem(e)->type() == etypes[idx]) &&
841  (mesh.elem(e)->level() == level) &&
842  !mesh.elem(e)->subactive())
843  {
844  int nstart=0;
845 
846  if (orig_type == LegacyXdrIO::DEAL)
847  nn = mesh.elem(e)->n_nodes();
848 
849  else if (orig_type == LegacyXdrIO::MGF)
850  {
851  nstart=2; // ignore the 27 and 0 entries
852  nn = mesh.elem(e)->n_nodes()+2;
853  conn[lastConnIndex + 0] = 27;
854  conn[lastConnIndex + 1] = 0;
855  }
856 
857  else if (orig_type == LegacyXdrIO::LIBM) // LIBMESH format
858  nn = mesh.elem(e)->n_nodes() + 2;
859 
860  else
861  libmesh_error();
862 
863  // Loop over the connectivity entries for this element and write to conn.
864  START_LOG("set connectivity", "LegacyXdrIO::write_mesh");
865  const unsigned int loopmax = (orig_type==LegacyXdrIO::LIBM) ? nn-2 : nn;
866  for (unsigned int n=nstart; n<loopmax; n++)
867  {
868  unsigned int connectivity_value=0;
869 
870  // old-style Libmesh and MGF meshes
871  if (orig_type != LegacyXdrIO::LIBM)
872  connectivity_value = mesh.elem(e)->node(n-nstart);
873 
874  // new-style libMesh meshes: compress the connectivity entries to account for
875  // subactive nodes that will not be in the mesh we write out.
876  else
877  {
878  std::map<dof_id_type, dof_id_type>::iterator pos =
879  node_map.find(mesh.elem(e)->node(n-nstart));
880 
881  libmesh_assert (pos != node_map.end());
882 
883  connectivity_value = (*pos).second;
884  }
885  conn[lastConnIndex + n] = connectivity_value;
886  }
887  STOP_LOG("set connectivity", "LegacyXdrIO::write_mesh");
888 
889  // In the case of an adaptive mesh, set last 2 entries to this ID and parent ID
890  if (orig_type == LegacyXdrIO::LIBM)
891  {
892  int self_ID = mesh.elem(e)->id();
893  int parent_ID = -1;
894  if(level != 0)
895  parent_ID = mesh.elem(e)->parent()->id();
896 
897  // Self ID is the second-to-last entry, Parent ID is the last
898  // entry on each connectivity line
899  conn[lastConnIndex+nn-2] = self_ID;
900  conn[lastConnIndex+nn-1] = parent_ID;
901  }
902 
903  lastConnIndex += nn;
904  }
905 
906  // Send conn to the XDR file. If there are no elements of this level and type,
907  // then nn will be zero, and we there is no connectivity to write.
908  if (nn != 0)
909  m.Icon(&conn[0], nn, lastConnIndex/nn);
910  }
911  }
912 
913  // create the vector of coords and send
914  // it to the XDR file.
915  {
916  std::vector<Real> coords;
917 
918  coords.resize(3*node_map.size());
919  int lastIndex=0;
920 
921  std::map<dof_id_type,dof_id_type>::iterator it = node_map.begin();
922  const std::map<dof_id_type,dof_id_type>::iterator end = node_map.end();
923  for (; it != end; ++it)
924  {
925  const Point& p = mesh.node((*it).first);
926 
927  coords[lastIndex+0] = p(0);
928  coords[lastIndex+1] = p(1);
929  coords[lastIndex+2] = p(2);
930  lastIndex += 3;
931  }
932 
933  // Put the nodes in the XDR file
934  m.coord(&coords[0], 3, node_map.size());
935  }
936 
937 
938  // write the
939  // boundary conditions
940  // if the BoundaryInfo
941  // object exists.
942  if (numBCs > 0)
943  {
944  std::vector<int> bcs(numBCs*3);
945 
946  //libMesh::out << "numBCs=" << numBCs << std::endl;
947 
948  //libMesh::out << "Preparing to write boundary conditions." << std::endl;
949  std::vector<dof_id_type> elem_list;
950  std::vector<unsigned short int> side_list;
951  std::vector<boundary_id_type> elem_id_list;
952 
953  mesh.boundary_info->build_side_list (elem_list, side_list, elem_id_list);
954 
955  for (int ibc=0; ibc<numBCs; ibc++)
956  {
957  bcs[0+ibc*3] = elem_list[ibc];
958  bcs[1+ibc*3] = side_list[ibc];
959  bcs[2+ibc*3] = elem_id_list[ibc];
960  }
961 
962  // Put the BCs in the XDR file
963  m.BC(&bcs[0], numBCs);
964  }
965 }
966 
967 
968 
969 void LegacyXdrIO::read_soln (const std::string& name,
970  std::vector<Real>& soln,
971  std::vector<std::string>& var_names) const
972 {
973  // Create an XdrSOLN object.
974  XdrSOLN s;
975 
976  // Create an XdrSHEAD object.
977  XdrSHEAD sh;
978 
979  // Open the XDR file for
980  // reading or writing.
981  // Note 1: Provide an additional argument
982  // to specify the dimension.
983  //
984  // Note 2: Has to do the right thing for
985  // both binary and ASCII files.
986  s.init((this->binary() ? XdrMGF::DECODE : XdrMGF::R_ASCII), name.c_str(), 0); // mesh files are always number 0 ...
987 
988  // From here on, things depend
989  // on whether we are reading or
990  // writing! First, we define
991  // header variables that may
992  // be read OR written.
993  int numVar = 0;
994  int numNodes = 0;
995  const char* varNames;
996 
997  // Get the information from
998  // the header, and place it
999  // in the header pointer.
1000  s.header(&sh);
1001 
1002  // Read information from the
1003  // file header. This depends on
1004  // whether its a libMesh or MGF mesh.
1005  numVar = sh.getWrtVar();
1006  numNodes = sh.getNumNodes();
1007  varNames = sh.getVarTitle();
1008 
1009  // Get the variable names
1010  {
1011  var_names.resize(numVar);
1012 
1013  const char* p = varNames;
1014 
1015  for (int i=0; i<numVar; i++)
1016  {
1017  var_names[i] = p;
1018  p += std::strlen(p) + 1;
1019  }
1020  }
1021 
1022  // Read the soln vector
1023  soln.resize(numVar*numNodes);
1024 
1025  s.values(&soln[0], numNodes);
1026 }
1027 
1028 
1029 
1030 void LegacyXdrIO::write_soln (const std::string& name,
1031  std::vector<Real>& soln,
1032  std::vector<std::string>& var_names) const
1033 {
1034  // get a read-only reference to the mesh
1035  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();
1036 
1037  // Create an XdrSOLN object.
1038  XdrSOLN s;
1039 
1040  // Create an XdrSHEAD object.
1041  XdrSHEAD sh;
1042 
1043  // Open the XDR file for
1044  // reading or writing.
1045  // Note 1: Provide an additional argument
1046  // to specify the dimension.
1047  //
1048  // Note 2: Has to do the right thing for
1049  // both binary and ASCII files.
1050  s.init((this->binary() ? XdrMGF::ENCODE : XdrMGF::W_ASCII), name.c_str(), 0); // mesh files are always number 0 ...
1051 
1052  // Build the header
1053  sh.setWrtVar(var_names.size());
1054  sh.setNumVar(var_names.size());
1055  sh.setNumNodes(mesh.n_nodes());
1056  sh.setNumBCs(mesh.boundary_info->n_boundary_conds());
1057  sh.setMeshCnt(0);
1058  sh.setKstep(0);
1059  sh.setTime(0.);
1060  sh.setStrSize(65536);
1061  sh.setId("Id String"); // Ignored
1062  sh.setTitle("Title String"); // Ignored
1063  sh.setUserTitle("User Title String"); // Ignored
1064 
1065  // create the variable array
1066  {
1067  std::string var_title;
1068 
1069  for (unsigned int var=0; var<var_names.size(); var++)
1070  {
1071  for (unsigned int c=0; c<var_names[var].size(); c++)
1072  var_title += var_names[var][c];
1073 
1074  var_title += '\0';
1075  }
1076 
1077  sh.setVarTitle(var_title.c_str(), var_title.size());
1078  }
1079 
1080  // Put the informationin the XDR file.
1081  s.header(&sh); // Needs to work for both types of file
1082 
1083  // Write the solution vector
1084  libmesh_assert_equal_to (soln.size(), var_names.size()*mesh.n_nodes());
1085 
1086  s.values(&soln[0], mesh.n_nodes());
1087 }
1088 
1089 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo