libMesh::Nemesis_IO Class Reference

#include <nemesis_io.h>

Inheritance diagram for libMesh::Nemesis_IO:

List of all members.

Public Member Functions

 Nemesis_IO (MeshBase &mesh)
virtual ~Nemesis_IO ()
virtual void read (const std::string &base_filename)
virtual void write (const std::string &base_filename)
void write_timestep (const std::string &fname, const EquationSystems &es, const int timestep, const Real time)
void write_nodal_data (const std::string &fname, const std::vector< Number > &soln, const std::vector< std::string > &names)
void verbose (bool set_verbosity)
void write_global_data (const std::vector< Number > &, const std::vector< std::string > &)
void write_information_records (const std::vector< std::string > &)
virtual void write_equation_systems (const std::string &, const EquationSystems &, const std::set< std::string > *system_names=NULL)
unsigned int & ascii_precision ()

Protected Member Functions

MeshBasemesh ()
void skip_comment_lines (std::istream &in, const char comment_start)
const MeshBasemesh () const

Protected Attributes

std::vector< bool > elems_of_dimension
const bool _is_parallel_format

Private Attributes

Nemesis_IO_Helpernemhelper
int _timestep
bool _verbose

Detailed Description

The Nemesis_IO class implements reading parallel meshes in the Nemesis file format from Sandia National Labs. Nemesis files are essentially in the Exodus format plus some additional information. All the Nemesis files for a single mesh have the same basename, e.g. cylinder.e, followed by ".size.rank", where size is the total number of files the Mesh is split into and rank is the ID of the processor's elements that were written to the file.

Author:
John Peterson, 2008.

Definition at line 52 of file nemesis_io.h.


Constructor & Destructor Documentation

libMesh::Nemesis_IO::Nemesis_IO ( MeshBase mesh  )  [explicit]

Constructor. Takes a writeable reference to a mesh object. This is the constructor required to read a mesh.

Definition at line 80 of file nemesis_io.C.

00080                                       :
00081   MeshInput<MeshBase> (mesh, /*is_parallel_format=*/true),
00082   MeshOutput<MeshBase> (mesh, /*is_parallel_format=*/true),
00083 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
00084   nemhelper(new Nemesis_IO_Helper),
00085 #endif
00086   _timestep(1),
00087   _verbose (false)
00088 {
00089 }

libMesh::Nemesis_IO::~Nemesis_IO (  )  [virtual]

Destructor.

Definition at line 92 of file nemesis_io.C.

References nemhelper.

00093 {
00094 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
00095 
00096   delete nemhelper;
00097 #endif
00098 }


Member Function Documentation

unsigned int& libMesh::MeshOutput< MeshBase >::ascii_precision (  )  [inherited]

Return/set the precision to use when writing ASCII files.

By default we use numeric_limits<Real>::digits10 + 2, which should be enough to write out to ASCII and get the exact same Real back when reading in.

Referenced by libMesh::TecplotIO::write_ascii(), libMesh::GMVIO::write_ascii_new_impl(), and libMesh::GMVIO::write_ascii_old_impl().

MeshBase & libMesh::MeshInput< MeshBase >::mesh (  )  [protected, inherited]

Returns the object as a writeable reference.

Referenced by libMesh::GMVIO::_read_materials(), libMesh::GMVIO::_read_nodes(), libMesh::GMVIO::_read_one_cell(), libMesh::AbaqusIO::assign_boundary_node_ids(), libMesh::AbaqusIO::assign_sideset_ids(), libMesh::AbaqusIO::assign_subdomain_ids(), libMesh::VTKIO::cells_to_vtk(), libMesh::GMVIO::copy_nodal_solution(), libMesh::UNVIO::element_in(), libMesh::TetGenIO::element_in(), libMesh::UNVIO::element_out(), libMesh::UNVIO::node_in(), libMesh::TetGenIO::node_in(), libMesh::UNVIO::node_out(), libMesh::VTKIO::nodes_to_vtk(), libMesh::XdrIO::read(), libMesh::VTKIO::read(), libMesh::TetGenIO::read(), read(), libMesh::GMVIO::read(), libMesh::ExodusII_IO::read(), libMesh::AbaqusIO::read(), libMesh::LegacyXdrIO::read_ascii(), libMesh::AbaqusIO::read_elements(), libMesh::UNVIO::read_implementation(), libMesh::UCDIO::read_implementation(), libMesh::LegacyXdrIO::read_mesh(), libMesh::GmshIO::read_mesh(), libMesh::AbaqusIO::read_nodes(), libMesh::XdrIO::read_serialized_bcs(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::OFFIO::read_stream(), libMesh::MatlabIO::read_stream(), libMesh::XdrIO::write(), libMesh::TetGenIO::write(), write(), libMesh::GMVIO::write_ascii_new_impl(), libMesh::GMVIO::write_ascii_old_impl(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::UNVIO::write_implementation(), libMesh::UCDIO::write_implementation(), libMesh::LegacyXdrIO::write_mesh(), libMesh::GmshIO::write_mesh(), libMesh::UCDIO::write_nodal_data(), write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::XdrIO::write_parallel(), libMesh::GmshIO::write_post(), libMesh::XdrIO::write_serialized_bcs(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::LegacyXdrIO::write_soln().

void libMesh::Nemesis_IO::read ( const std::string &  base_filename  )  [virtual]

Implements reading the mesh from several different files. You provide the basename, then LibMesh appends the ".size.rank" depending on libMesh::n_processors() and libMesh::processor_id().

Implements libMesh::MeshInput< MeshBase >.

Definition at line 116 of file nemesis_io.C.

References _verbose, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::Parallel::Communicator::allgather(), libMesh::Parallel::Communicator::alltoall(), libMesh::Parallel::any_source, libMesh::ExodusII_IO_Helper::ElementMaps::assign_conversion(), libMesh::ExodusII_IO_Helper::block_ids, libMesh::MeshBase::boundary_info, libMesh::Elem::build(), libMesh::Parallel::Communicator_World, libMesh::CommWorld, libMesh::ExodusII_IO_Helper::connect, libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::Utility::deallocate(), libMesh::MeshBase::delete_remote_elements(), libMesh::Elem::dim(), libMesh::MeshBase::elem(), libMesh::ExodusII_IO_Helper::elem_type, libMesh::MeshInput< MeshBase >::elems_of_dimension, libMesh::err, libMesh::ExodusII_IO_Helper::Conversion::get_canonical_type(), libMesh::Nemesis_IO_Helper::get_cmap_params(), libMesh::Nemesis_IO_Helper::get_eb_info_global(), libMesh::Nemesis_IO_Helper::get_elem_cmap(), libMesh::ExodusII_IO_Helper::get_elem_list(), libMesh::ExodusII_IO_Helper::get_id_list(), libMesh::Nemesis_IO_Helper::get_init_global(), libMesh::Nemesis_IO_Helper::get_loadbal_param(), libMesh::Nemesis_IO_Helper::get_node_cmap(), libMesh::ExodusII_IO_Helper::get_node_list(), libMesh::ExodusII_IO_Helper::Conversion::get_node_map(), libMesh::Nemesis_IO_Helper::get_node_map(), libMesh::ExodusII_IO_Helper::get_nodeset_id(), libMesh::Nemesis_IO_Helper::get_ns_param_global(), libMesh::ExodusII_IO_Helper::get_num_node_sets(), libMesh::ExodusII_IO_Helper::get_num_side_sets(), libMesh::ExodusII_IO_Helper::get_num_sides_per_set(), libMesh::ExodusII_IO_Helper::get_side_list(), libMesh::ExodusII_IO_Helper::Conversion::get_side_map(), libMesh::Nemesis_IO_Helper::get_ss_param_global(), libMesh::Nemesis_IO_Helper::global_sideset_ids, libMesh::DofObject::id(), libMesh::MeshTools::Generation::Private::idx(), libMesh::MeshBase::is_serial(), libMesh::Parallel::Communicator::max(), libMesh::MeshInput< MeshBase >::mesh(), libMesh::MeshBase::mesh_dimension(), std::min(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::n_processors(), nemhelper, libMesh::Nemesis_IO_Helper::node_cmap_ids, libMesh::Nemesis_IO_Helper::node_cmap_node_cnts, libMesh::Nemesis_IO_Helper::node_cmap_node_ids, libMesh::Nemesis_IO_Helper::node_cmap_proc_ids, libMesh::Nemesis_IO_Helper::node_mapb, libMesh::Nemesis_IO_Helper::node_mape, libMesh::Nemesis_IO_Helper::node_mapi, libMesh::ExodusII_IO_Helper::node_num_map, libMesh::MeshBase::node_ptr(), libMesh::Nemesis_IO_Helper::num_border_elems, libMesh::Nemesis_IO_Helper::num_border_nodes, libMesh::ExodusII_IO_Helper::num_elem, libMesh::ExodusII_IO_Helper::num_elem_all_sidesets, libMesh::ExodusII_IO_Helper::num_elem_blk, libMesh::Nemesis_IO_Helper::num_elem_cmaps, libMesh::ExodusII_IO_Helper::num_elem_this_blk, libMesh::Nemesis_IO_Helper::num_elems_global, libMesh::Nemesis_IO_Helper::num_external_nodes, libMesh::Nemesis_IO_Helper::num_global_side_counts, libMesh::Nemesis_IO_Helper::num_internal_elems, libMesh::Nemesis_IO_Helper::num_internal_nodes, libMesh::Nemesis_IO_Helper::num_node_cmaps, libMesh::ExodusII_IO_Helper::num_node_sets, libMesh::ExodusII_IO_Helper::num_nodes, libMesh::Nemesis_IO_Helper::num_nodes_global, libMesh::ExodusII_IO_Helper::num_nodes_per_elem, libMesh::ExodusII_IO_Helper::num_side_sets, libMesh::ExodusII_IO_Helper::open(), libMesh::out, libMesh::MeshBase::parallel_n_elem(), libMesh::MeshBase::parallel_n_nodes(), libMesh::ExodusII_IO_Helper::print_header(), libMesh::Parallel::Communicator::probe(), libMesh::DofObject::processor_id(), libMesh::processor_id(), libMesh::ExodusII_IO_Helper::read_block_info(), libMesh::ExodusII_IO_Helper::read_elem_in_block(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::ExodusII_IO_Helper::read_header(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::ExodusII_IO_Helper::read_nodes(), libMesh::ExodusII_IO_Helper::read_nodeset(), libMesh::ExodusII_IO_Helper::read_nodeset_info(), libMesh::ExodusII_IO_Helper::read_sideset(), libMesh::ExodusII_IO_Helper::read_sideset_info(), libMesh::Parallel::Communicator::receive(), libMesh::Parallel::Communicator::send(), libMesh::DofObject::set_id(), libMesh::MeshBase::set_mesh_dimension(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), libMesh::Parallel::Communicator::sum(), libMesh::Elem::type(), libMesh::MeshBase::update_post_partitioning(), libMesh::Parallel::wait(), libMesh::ExodusII_IO_Helper::x, libMesh::ExodusII_IO_Helper::y, and libMesh::ExodusII_IO_Helper::z.

00117 {
00118   // On one processor, Nemesis and ExodusII should be equivalent, so
00119   // let's cowardly defer to that implementation...
00120   if (libMesh::n_processors() == 1)
00121     {
00122       // We can do this in one line but if the verbose flag was set in this
00123       // object, it will no longer be set... thus no extra print-outs for serial runs.
00124       // ExodusII_IO(this->mesh()).read (base_filename); // ambiguous when Nemesis_IO is multiply-inherited
00125 
00126       MeshBase& mesh = MeshInput<MeshBase>::mesh();
00127       ExodusII_IO(mesh).read (base_filename);
00128       return;
00129     }
00130 
00131   START_LOG ("read()","Nemesis_IO");
00132 
00133   // This function must be run on all processors at once
00134   parallel_only();
00135 
00136   if (_verbose)
00137     {
00138       libMesh::out << "[" << libMesh::processor_id() << "] ";
00139       libMesh::out << "Reading Nemesis file on processor: " << libMesh::processor_id() << std::endl;
00140     }
00141 
00142   // Construct the Nemesis filename based on the number of processors and the
00143   // current processor ID.
00144   std::string nemesis_filename = nemhelper->construct_nemesis_filename(base_filename);
00145 
00146   if (_verbose)
00147     libMesh::out << "Opening file: " << nemesis_filename << std::endl;
00148 
00149   // Open the Exodus file
00150   nemhelper->open(nemesis_filename.c_str());
00151 
00152   // Get a reference to the mesh.  We need to be specific
00153   // since Nemesis_IO is multiply-inherited
00154   // MeshBase& mesh = this->mesh();
00155   MeshBase& mesh = MeshInput<MeshBase>::mesh();
00156 
00157   // Local information: Read the following information from the standard Exodus header
00158   //  title[0]
00159   //  num_dim
00160   //  num_nodes
00161   //  num_elem
00162   //  num_elem_blk
00163   //  num_node_sets
00164   //  num_side_sets
00165   nemhelper->read_header();
00166   nemhelper->print_header();
00167 
00168   // Get global information: number of nodes, elems, blocks, nodesets and sidesets
00169   nemhelper->get_init_global();
00170 
00171   // Get "load balance" information.  This includes the number of internal & border
00172   // nodes and elements as well as the number of communication maps.
00173   nemhelper->get_loadbal_param();
00174 
00175   // Do some error checking
00176   if (nemhelper->num_external_nodes)
00177     {
00178       libMesh::err << "ERROR: there should be no external nodes in an element-based partitioning!"
00179                     << std::endl;
00180       libmesh_error();
00181     }
00182 
00183   libmesh_assert_equal_to (nemhelper->num_nodes,
00184                            (nemhelper->num_internal_nodes +
00185                             nemhelper->num_border_nodes));
00186 
00187   libmesh_assert_equal_to (nemhelper->num_elem,
00188                            (nemhelper->num_internal_elems +
00189                             nemhelper->num_border_elems));
00190 
00191   libmesh_assert_less_equal (nemhelper->num_nodes, nemhelper->num_nodes_global);
00192   libmesh_assert_less_equal (nemhelper->num_elem, nemhelper->num_elems_global);
00193 
00194   // Read nodes from the exodus file: this fills the nemhelper->x,y,z arrays.
00195   nemhelper->read_nodes();
00196 
00197   // Reads the nemhelper->node_num_map array, node_num_map[i] is the global node number for
00198   // local node number i.
00199   nemhelper->read_node_num_map();
00200 
00201   // The get_cmap_params() function reads in the:
00202   //  node_cmap_ids[],
00203   //  node_cmap_node_cnts[],
00204   //  elem_cmap_ids[],
00205   //  elem_cmap_elem_cnts[],
00206   nemhelper->get_cmap_params();
00207 
00208   // Read the IDs of the interior, boundary, and external nodes.  This function
00209   // fills the vectors:
00210   //  node_mapi[],
00211   //  node_mapb[],
00212   //  node_mape[]
00213   nemhelper->get_node_map();
00214 
00215   // Read each node communication map for this processor.  This function
00216   // fills the vectors of vectors named:
00217   //  node_cmap_node_ids[][]
00218   //  node_cmap_proc_ids[][]
00219   nemhelper->get_node_cmap();
00220 
00221   libmesh_assert_equal_to (to_uint(nemhelper->num_node_cmaps), nemhelper->node_cmap_node_cnts.size());
00222   libmesh_assert_equal_to (to_uint(nemhelper->num_node_cmaps), nemhelper->node_cmap_node_ids.size());
00223   libmesh_assert_equal_to (to_uint(nemhelper->num_node_cmaps), nemhelper->node_cmap_proc_ids.size());
00224 
00225 #ifndef NDEBUG
00226   // We expect the communication maps to be symmetric - e.g. if processor i thinks it
00227   // communicates with processor j, then processor j should also be expecting to
00228   // communicate with i.  We can assert that here easily enough with an alltoall,
00229   // but let's only do it when not in optimized mode to limit unnecessary communication.
00230   {
00231     std::vector<unsigned char> pid_send_partner (libMesh::n_processors(), 0);
00232 
00233     // strictly speaking, we should expect to communicate with ourself...
00234     pid_send_partner[libMesh::processor_id()] = 1;
00235 
00236     // mark each processor id we reference with a node cmap
00237     for (unsigned int cmap=0; cmap<to_uint(nemhelper->num_node_cmaps); cmap++)
00238       {
00239         libmesh_assert_less (to_uint(nemhelper->node_cmap_ids[cmap]), libMesh::n_processors());
00240 
00241         pid_send_partner[nemhelper->node_cmap_ids[cmap]] = 1;
00242       }
00243 
00244     // Copy the send pairing so we can catch the receive paring and
00245     // test for equality
00246     const std::vector<unsigned char> pid_recv_partner (pid_send_partner);
00247 
00248     CommWorld.alltoall (pid_send_partner);
00249 
00250     libmesh_assert (pid_send_partner == pid_recv_partner);
00251   }
00252 #endif
00253 
00254   // We now have enough information to infer node ownership.  We start by assuming
00255   // we own all the nodes on this processor.  We will then interrogate the
00256   // node cmaps and see if a lower-rank processor is associated with any of
00257   // our nodes.  If so, then that processor owns the node, not us...
00258   std::vector<unsigned short int> node_ownership (nemhelper->num_internal_nodes +
00259                                                   nemhelper->num_border_nodes,
00260                                                   libMesh::processor_id());
00261 
00262   // a map from processor id to cmap number, to be used later
00263   std::map<unsigned int, unsigned int> pid_to_cmap_map;
00264 
00265   // For each node_cmap...
00266   for (unsigned int cmap=0; cmap<to_uint(nemhelper->num_node_cmaps); cmap++)
00267     {
00268       // Good time for error checking...
00269       libmesh_assert_equal_to (to_uint(nemhelper->node_cmap_node_cnts[cmap]),
00270                                nemhelper->node_cmap_node_ids[cmap].size());
00271 
00272       libmesh_assert_equal_to (to_uint(nemhelper->node_cmap_node_cnts[cmap]),
00273                                nemhelper->node_cmap_proc_ids[cmap].size());
00274 
00275       // In all the samples I have seen, node_cmap_ids[cmap] is the processor
00276       // rank of the remote processor...
00277       const unsigned short int adjcnt_pid_idx = nemhelper->node_cmap_ids[cmap];
00278 
00279       libmesh_assert_less (adjcnt_pid_idx, libMesh::n_processors());
00280       libmesh_assert_not_equal_to (adjcnt_pid_idx, libMesh::processor_id());
00281 
00282       // We only expect one cmap per adjacent processor
00283       libmesh_assert (!pid_to_cmap_map.count(adjcnt_pid_idx));
00284 
00285       pid_to_cmap_map[adjcnt_pid_idx] = cmap;
00286 
00287       // ...and each node in that cmap...
00288       for (unsigned int idx=0; idx<to_uint(nemhelper->node_cmap_node_cnts[cmap]); idx++)
00289         {
00290           //  Are the node_cmap_ids and node_cmap_proc_ids really redundant?
00291           libmesh_assert_equal_to (adjcnt_pid_idx, nemhelper->node_cmap_proc_ids[cmap][idx]);
00292 
00293           // we are expecting the exodus node numbering to be 1-based...
00294           const unsigned int local_node_idx = nemhelper->node_cmap_node_ids[cmap][idx]-1;
00295 
00296           libmesh_assert_less (local_node_idx, node_ownership.size());
00297 
00298           // if the adjacent processor is lower rank than the current
00299           // owner for this node, then it will get the node...
00300           node_ownership[local_node_idx] =
00301             std::min(node_ownership[local_node_idx], adjcnt_pid_idx);
00302         }
00303     } // We now should have established proper node ownership.
00304 
00305   // now that ownership is established, we can figure out how many nodes we
00306   // will be responsible for numbering.
00307   unsigned int num_nodes_i_must_number = 0;
00308 
00309   for (unsigned int idx=0; idx<node_ownership.size(); idx++)
00310     if (node_ownership[idx] == libMesh::processor_id())
00311       num_nodes_i_must_number++;
00312 
00313   // more error checking...
00314   libmesh_assert_greater_equal (num_nodes_i_must_number, to_uint(nemhelper->num_internal_nodes));
00315   libmesh_assert (num_nodes_i_must_number <= to_uint(nemhelper->num_internal_nodes +
00316                                                   nemhelper->num_border_nodes));
00317   if (_verbose)
00318     libMesh::out << "[" << libMesh::processor_id() << "] "
00319                   << "num_nodes_i_must_number="
00320                   << num_nodes_i_must_number
00321                   << std::endl;
00322 
00323   // The call to get_loadbal_param() gets 7 pieces of information.  We allgather
00324   // these now across all processors to determine some global numberings. We should
00325   // also gather the number of nodes each processor thinks it will number so that
00326   // we can (i) determine our offset, and (ii) do some error checking.
00327   std::vector<int> all_loadbal_data ( 8 );
00328   all_loadbal_data[0] = nemhelper->num_internal_nodes;
00329   all_loadbal_data[1] = nemhelper->num_border_nodes;
00330   all_loadbal_data[2] = nemhelper->num_external_nodes;
00331   all_loadbal_data[3] = nemhelper->num_internal_elems;
00332   all_loadbal_data[4] = nemhelper->num_border_elems;
00333   all_loadbal_data[5] = nemhelper->num_node_cmaps;
00334   all_loadbal_data[6] = nemhelper->num_elem_cmaps;
00335   all_loadbal_data[7] = num_nodes_i_must_number;
00336 
00337   CommWorld.allgather (all_loadbal_data, /* identical_buffer_sizes = */ true);
00338 
00339   // OK, we are now in a position to request new global indices for all the nodes
00340   // we do not own
00341 
00342   // Let's get a unique message tag to use for send()/receive()
00343   Parallel::MessageTag nodes_tag = Parallel::Communicator_World.get_unique_tag(12345);
00344 
00345   std::vector<std::vector<int> >
00346     needed_node_idxs (nemhelper->num_node_cmaps); // the indices we will ask for
00347 
00348   std::vector<Parallel::Request>
00349     needed_nodes_requests (nemhelper->num_node_cmaps);
00350 
00351   for (unsigned int cmap=0; cmap<to_uint(nemhelper->num_node_cmaps); cmap++)
00352     {
00353       // We know we will need no more indices than there are nodes
00354       // in this cmap, but that number is an upper bound in general
00355       // since the neighboring processor associated with the cmap
00356       //  may not actually own it
00357       needed_node_idxs[cmap].reserve   (nemhelper->node_cmap_node_cnts[cmap]);
00358 
00359       const unsigned int adjcnt_pid_idx = nemhelper->node_cmap_ids[cmap];
00360 
00361       // ...and each node in that cmap...
00362       for (unsigned int idx=0; idx<to_uint(nemhelper->node_cmap_node_cnts[cmap]); idx++)
00363         {
00364           const unsigned int
00365             local_node_idx  = nemhelper->node_cmap_node_ids[cmap][idx]-1,
00366             owning_pid_idx  = node_ownership[local_node_idx];
00367 
00368           // add it to the request list for its owning processor.
00369           if (owning_pid_idx == adjcnt_pid_idx)
00370             {
00371               const unsigned int
00372                 global_node_idx = nemhelper->node_num_map[local_node_idx]-1;
00373               needed_node_idxs[cmap].push_back(global_node_idx);
00374             }
00375         }
00376       // now post the send for this cmap
00377       CommWorld.send (adjcnt_pid_idx,              // destination
00378                       needed_node_idxs[cmap],      // send buffer
00379                       needed_nodes_requests[cmap], // request
00380                       nodes_tag);
00381     } // all communication requests for getting updated global indices for border
00382       // nodes have been initiated
00383 
00384   // Figure out how many nodes each processor thinks it will number and make sure
00385   // that it adds up to the global number of nodes. Also, set up global node
00386   // index offsets for each processor.
00387   std::vector<unsigned int>
00388     all_num_nodes_i_must_number (libMesh::n_processors());
00389 
00390   for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
00391     all_num_nodes_i_must_number[pid] = all_loadbal_data[8*pid + 7];
00392 
00393   // The sum of all the entries in this vector should sum to the number of global nodes
00394   libmesh_assert (std::accumulate(all_num_nodes_i_must_number.begin(),
00395                                   all_num_nodes_i_must_number.end(),
00396                                   0) == nemhelper->num_nodes_global);
00397 
00398   unsigned int my_next_node = 0;
00399   for (unsigned int pid=0; pid<libMesh::processor_id(); pid++)
00400     my_next_node += all_num_nodes_i_must_number[pid];
00401 
00402   const unsigned int my_node_offset = my_next_node;
00403 
00404   if (_verbose)
00405     libMesh::out << "[" << libMesh::processor_id() << "] "
00406                   << "my_node_offset="
00407                   << my_node_offset
00408                   << std::endl;
00409 
00410   // Add internal nodes to the ParallelMesh, using the node ID offset we
00411   // computed and the current processor's ID.
00412   for (unsigned int i=0; i<to_uint(nemhelper->num_internal_nodes); ++i)
00413     {
00414       const unsigned int local_node_idx  = nemhelper->node_mapi[i]-1;
00415 #ifndef NDEBUG
00416       const unsigned int global_node_idx = nemhelper->node_num_map[local_node_idx]-1;
00417       const unsigned int owning_pid_idx  = node_ownership[local_node_idx];
00418 #endif
00419 
00420       // an internal node we do not own? huh??
00421       libmesh_assert_equal_to (owning_pid_idx, libMesh::processor_id());
00422       libmesh_assert_less (global_node_idx, to_uint(nemhelper->num_nodes_global));
00423 
00424       // "Catch" the node pointer after addition, make sure the
00425       // ID matches the requested value.
00426       Node* added_node =
00427         mesh.add_point (Point(nemhelper->x[local_node_idx],
00428                               nemhelper->y[local_node_idx],
00429                               nemhelper->z[local_node_idx]),
00430                         my_next_node,
00431                         libMesh::processor_id());
00432 
00433       // Make sure the node we added has the ID we thought it would
00434       if (added_node->id() != my_next_node)
00435         {
00436           libMesh::err << "Error, node added with ID " << added_node->id()
00437                        << ", but we wanted ID " << my_next_node << std::endl;
00438         }
00439 
00440       // update the local->global index map, when we are done
00441       // it will be 0-based.
00442       nemhelper->node_num_map[local_node_idx] = my_next_node++;
00443     }
00444 
00445   // Now, for the boundary nodes...  We may very well own some of them,
00446   // but there may be others for which we have requested the new global
00447   // id.  We expect to be asked for the ids of the ones we own, so
00448   // we need to create a map from the old global id to the new one
00449   // we are about to create.
00450   typedef std::vector<std::pair<unsigned int, unsigned int> > global_idx_mapping_type;
00451   global_idx_mapping_type old_global_to_new_global_map;
00452   old_global_to_new_global_map.reserve (num_nodes_i_must_number // total # i will have
00453                                         - (my_next_node         // amount i have thus far
00454                                            - my_node_offset));  // this should be exact!
00455   CompareGlobalIdxMappings global_idx_mapping_comp;
00456 
00457   for (unsigned int i=0; i<to_uint(nemhelper->num_border_nodes); ++i)
00458     {
00459       const unsigned int
00460         local_node_idx  = nemhelper->node_mapb[i]-1,
00461         owning_pid_idx  = node_ownership[local_node_idx];
00462 
00463       // if we own it...
00464       if (owning_pid_idx == libMesh::processor_id())
00465         {
00466           const unsigned int
00467             global_node_idx = nemhelper->node_num_map[local_node_idx]-1;
00468 
00469           // we will number it, and create a mapping from its old global index to
00470           // the new global index, for lookup purposes when neighbors come calling
00471           old_global_to_new_global_map.push_back(std::make_pair(global_node_idx,
00472                                                                 my_next_node));
00473 
00474           // "Catch" the node pointer after addition, make sure the
00475           // ID matches the requested value.
00476           Node* added_node =
00477             mesh.add_point (Point(nemhelper->x[local_node_idx],
00478                                   nemhelper->y[local_node_idx],
00479                                   nemhelper->z[local_node_idx]),
00480                             my_next_node,
00481                             libMesh::processor_id());
00482 
00483           // Make sure the node we added has the ID we thought it would
00484           if (added_node->id() != my_next_node)
00485             {
00486               libMesh::err << "Error, node added with ID " << added_node->id()
00487                            << ", but we wanted ID " << my_next_node << std::endl;
00488             }
00489 
00490           // update the local->global index map, when we are done
00491           // it will be 0-based.
00492           nemhelper->node_num_map[local_node_idx] = my_next_node++;
00493         }
00494     }
00495   // That should cover numbering all the nodes which belong to us...
00496   libmesh_assert_equal_to (num_nodes_i_must_number, (my_next_node - my_node_offset));
00497 
00498   // Let's sort the mapping so we can efficiently answer requests
00499   std::sort (old_global_to_new_global_map.begin(),
00500              old_global_to_new_global_map.end(),
00501              global_idx_mapping_comp);
00502 
00503   // and it had better be unique...
00504   libmesh_assert (std::unique (old_global_to_new_global_map.begin(),
00505                                old_global_to_new_global_map.end(),
00506                                global_idx_mapping_equality) ==
00507                   old_global_to_new_global_map.end());
00508 
00509   // We can now catch incoming requests and process them. for efficiency
00510   // let's do whatever is available next
00511   std::map<unsigned int, std::vector<int> > requested_node_idxs; // the indices asked of us
00512 
00513   std::vector<Parallel::Request> requested_nodes_requests(nemhelper->num_node_cmaps);
00514 
00515   // We know we will receive the request from a given processor before
00516   // we receive its reply to our request. However, we may receive
00517   // a request and a response from one processor before getting
00518   // a request from another processor.  So what we are doing here
00519   // is processing whatever message comes next, while recognizing
00520   // we will receive a request from a processor before receiving
00521   // its reply
00522   std::vector<bool> processed_cmap (nemhelper->num_node_cmaps, false);
00523 
00524   for (unsigned int comm_step=0; comm_step<2*to_uint(nemhelper->num_node_cmaps); comm_step++)
00525     {
00526       // query the first message which is available
00527       const Parallel::Status
00528         status (CommWorld.probe (Parallel::any_source,
00529                                  nodes_tag));
00530       const unsigned int
00531         requesting_pid_idx = status.source(),
00532         source_pid_idx     = status.source();
00533 
00534       // this had better be from a processor we are expecting...
00535       libmesh_assert (pid_to_cmap_map.count(requesting_pid_idx));
00536 
00537       // the local cmap which corresponds to the source processor
00538       const unsigned int cmap = pid_to_cmap_map[source_pid_idx];
00539 
00540       if (!processed_cmap[cmap])
00541         {
00542           processed_cmap[cmap] = true;
00543 
00544           // we should only get one request per paired processor
00545           libmesh_assert (!requested_node_idxs.count(requesting_pid_idx));
00546 
00547           // get a reference to the request buffer for this processor to
00548           // avoid repeated map lookups
00549           std::vector<int> &xfer_buf (requested_node_idxs[requesting_pid_idx]);
00550 
00551           // actually receive the message.
00552           CommWorld.receive (requesting_pid_idx, xfer_buf, nodes_tag);
00553 
00554           // Fill the request
00555           for (unsigned int i=0; i<xfer_buf.size(); i++)
00556             {
00557               // the requested old global node index, *now 0-based*
00558               const unsigned int old_global_node_idx = xfer_buf[i];
00559 
00560               // find the new global node index for the requested node -
00561               // note that requesting_pid_idx thinks we own this node,
00562               // so we better!
00563               const global_idx_mapping_type::const_iterator it =
00564                 std::lower_bound (old_global_to_new_global_map.begin(),
00565                                   old_global_to_new_global_map.end(),
00566                                   old_global_node_idx,
00567                                   global_idx_mapping_comp);
00568 
00569               libmesh_assert (it != old_global_to_new_global_map.end());
00570               libmesh_assert_equal_to (it->first, old_global_node_idx);
00571               libmesh_assert_greater_equal (it->second, my_node_offset);
00572               libmesh_assert_less (it->second, my_next_node);
00573 
00574               // overwrite the requested old global node index with the new global index
00575               xfer_buf[i] = it->second;
00576             }
00577 
00578           // and send the new global indices back to the processor which asked for them
00579           CommWorld.send (requesting_pid_idx,
00580                           xfer_buf,
00581                           requested_nodes_requests[cmap],
00582                           nodes_tag);
00583         } // done processing the request
00584 
00585       // this is the second time we have heard from this processor,
00586       // so it must be its reply to our request
00587       else
00588         {
00589           // a long time ago, we sent off our own requests.  now it is time to catch the
00590           // replies and get the new global node numbering.  note that for any reply
00591           // we receive, the corresponding nonblocking send from above *must* have been
00592           // completed, since the reply is in response to that request!!
00593 
00594           // if we have received a reply, our send *must* have completed
00595           // (note we never actually need to wait on the request)
00596           libmesh_assert (needed_nodes_requests[cmap].test());
00597           libmesh_assert_equal_to (to_uint(nemhelper->node_cmap_ids[cmap]), source_pid_idx);
00598 
00599           // now post the receive for this cmap
00600           CommWorld.receive (source_pid_idx,
00601                              needed_node_idxs[cmap],
00602                              nodes_tag);
00603 
00604           libmesh_assert_less_equal (needed_node_idxs[cmap].size(),
00605                                      nemhelper->node_cmap_node_ids[cmap].size());
00606 
00607           for (unsigned int i=0,j=0; i<nemhelper->node_cmap_node_ids[cmap].size(); i++)
00608             {
00609               const unsigned int
00610                 local_node_idx  = nemhelper->node_cmap_node_ids[cmap][i]-1,
00611                 owning_pid_idx  = node_ownership[local_node_idx];
00612 
00613               // if this node is owned by source_pid_idx, its new global id
00614               // is in the buffer we just received
00615               if (owning_pid_idx == source_pid_idx)
00616                 {
00617                   libmesh_assert_less (j, needed_node_idxs[cmap].size());
00618 
00619                   const unsigned int // now 0-based!
00620                     global_node_idx = needed_node_idxs[cmap][j++];
00621 
00622                   // "Catch" the node pointer after addition, make sure the
00623                   // ID matches the requested value.
00624                   Node* added_node =
00625                     mesh.add_point (Point(nemhelper->x[local_node_idx],
00626                                           nemhelper->y[local_node_idx],
00627                                           nemhelper->z[local_node_idx]),
00628                                     global_node_idx,
00629                                     source_pid_idx);
00630 
00631                   // Make sure the node we added has the ID we thought it would
00632                   if (added_node->id() != global_node_idx)
00633                     {
00634                       libMesh::err << "Error, node added with ID " << added_node->id()
00635                                    << ", but we wanted ID " << global_node_idx << std::endl;
00636                     }
00637 
00638                   // update the local->global index map, when we are done
00639                   // it will be 0-based.
00640                   nemhelper->node_num_map[local_node_idx] = global_node_idx;
00641 
00642                   // we are not really going to use my_next_node again, but we can
00643                   // keep incrimenting it to track how many nodes we have added
00644                   // to the mesh
00645                   my_next_node++;
00646                 }
00647             }
00648         }
00649     } // end of node index communication loop
00650 
00651   // we had better have added all the nodes we need to!
00652   libmesh_assert_equal_to ((my_next_node - my_node_offset), to_uint(nemhelper->num_nodes));
00653 
00654   // After all that, we should be done with all node-related arrays *except* the
00655   // node_num_map, which we have transformed to use our new numbering...
00656   // So let's clean up the arrays we are done with.
00657   {
00658     Utility::deallocate (nemhelper->node_mapi);
00659     Utility::deallocate (nemhelper->node_mapb);
00660     Utility::deallocate (nemhelper->node_mape);
00661     Utility::deallocate (nemhelper->node_cmap_ids);
00662     Utility::deallocate (nemhelper->node_cmap_node_cnts);
00663     Utility::deallocate (nemhelper->node_cmap_node_ids);
00664     Utility::deallocate (nemhelper->node_cmap_proc_ids);
00665     Utility::deallocate (nemhelper->x);
00666     Utility::deallocate (nemhelper->y);
00667     Utility::deallocate (nemhelper->z);
00668     Utility::deallocate (needed_node_idxs);
00669     Utility::deallocate (node_ownership);
00670   }
00671 
00672   Parallel::wait (requested_nodes_requests);
00673   requested_node_idxs.clear();
00674 
00675   // See what the node count is up to now.
00676   if (_verbose)
00677     {
00678       // Report the number of nodes which have been added locally
00679       libMesh::out << "[" << libMesh::processor_id() << "] ";
00680       libMesh::out << "mesh.n_nodes()=" << mesh.n_nodes() << std::endl;
00681 
00682       // Reports the number of nodes that have been added in total.
00683       libMesh::out << "[" << libMesh::processor_id() << "] ";
00684       libMesh::out << "mesh.parallel_n_nodes()=" << mesh.parallel_n_nodes() << std::endl;
00685     }
00686 
00687 
00688 
00689   // --------------------------------------------------------------------------------
00690   // --------------------------------------------------------------------------------
00691   // --------------------------------------------------------------------------------
00692 
00693 
00694   // We can now read in the elements...Exodus stores them in blocks in which all
00695   // elements have the same geometric type.  This code is adapted directly from exodusII_io.C
00696 
00697   // Assertion: The sum of the border and internal elements on all processors
00698   // should equal nemhelper->num_elems_global
00699 #ifndef NDEBUG
00700   {
00701     int sum_internal_elems=0, sum_border_elems=0;
00702     for (unsigned int j=3,c=0; c<libMesh::n_processors(); j+=8,++c)
00703       sum_internal_elems += all_loadbal_data[j];
00704 
00705     for (unsigned int j=4,c=0; c<libMesh::n_processors(); j+=8,++c)
00706       sum_border_elems += all_loadbal_data[j];
00707 
00708     if (_verbose)
00709       {
00710         libMesh::out << "[" << libMesh::processor_id() << "] ";
00711         libMesh::out << "sum_internal_elems=" << sum_internal_elems << std::endl;
00712 
00713         libMesh::out << "[" << libMesh::processor_id() << "] ";
00714         libMesh::out << "sum_border_elems=" << sum_border_elems << std::endl;
00715       }
00716 
00717     libmesh_assert_equal_to (sum_internal_elems+sum_border_elems, nemhelper->num_elems_global);
00718   }
00719 #endif
00720 
00721   // We need to set the mesh dimension, but the following...
00722   // mesh.set_mesh_dimension(static_cast<unsigned int>(nemhelper->num_dim));
00723 
00724   // ... is not sufficient since some codes report num_dim==3 for two dimensional
00725   // meshes living in 3D, even though all the elements are of 2D type.  Therefore,
00726   // we instead use the dimension of the highest element found for the Mesh dimension,
00727   // similar to what is done by the Exodus reader, except here it requires a
00728   // parallel communication.
00729   elems_of_dimension.resize(4, false); // will use 1-based
00730 
00731   // Compute my_elem_offset, the amount by which to offset the local elem numbering
00732   // on my processor.
00733   unsigned int my_next_elem = 0;
00734   for (unsigned int pid=0; pid<libMesh::processor_id(); ++pid)
00735     my_next_elem += (all_loadbal_data[8*pid + 3]+  // num_internal_elems, proc pid
00736                      all_loadbal_data[8*pid + 4]); // num_border_elems, proc pid
00737   const unsigned int my_elem_offset = my_next_elem;
00738 
00739   if (_verbose)
00740     libMesh::out << "[" << libMesh::processor_id() << "] "
00741               << "my_elem_offset=" << my_elem_offset << std::endl;
00742 
00743 
00744   // Fills in the:
00745   // global_elem_blk_ids[] and
00746   // global_elem_blk_cnts[] arrays.
00747   nemhelper->get_eb_info_global();
00748 
00749 //   // Fills in the vectors
00750 //   // elem_mapi[num_internal_elems]
00751 //   // elem_mapb[num_border_elems  ]
00752 //   // These tell which of the (locally-numbered) elements are internal and which are border elements.
00753 //   // In our test example these arrays are sorted (but non-contiguous), which makes it possible to
00754 //   // binary search for each element ID... however I don't think we need to distinguish between the
00755 //   // two types, since either can have nodes the boundary!
00756 //   nemhelper->get_elem_map();
00757 
00758   // Fills in the vectors of vectors:
00759   // elem_cmap_elem_ids[][]
00760   // elem_cmap_side_ids[][]
00761   // elem_cmap_proc_ids[][]
00762   // These arrays are of size num_elem_cmaps * elem_cmap_elem_cnts[i], i = 0..num_elem_cmaps
00763   nemhelper->get_elem_cmap();
00764 
00765   // Get information about the element blocks:
00766   // (read in the array nemhelper->block_ids[])
00767   nemhelper->read_block_info();
00768 
00769   // Reads the nemhelper->elem_num_map array, elem_num_map[i] is the global element number for
00770   // local element number i.
00771   nemhelper->read_elem_num_map();
00772 
00773   // Instantiate the ElementMaps interface.  This is what translates LibMesh's
00774   // element numbering scheme to Exodus's.
00775   ExodusII_IO_Helper::ElementMaps em;
00776 
00777   // Read in the element connectivity for each block by
00778   // looping over all the blocks.
00779   for (unsigned int i=0; i<to_uint(nemhelper->num_elem_blk); i++)
00780     {
00781       // Read the information for block i:  For nemhelper->block_ids[i], reads
00782       // elem_type
00783       // num_elem_this_blk
00784       // num_nodes_per_elem
00785       // num_attr
00786       // connect <-- the nodal connectivity array for each element in the block.
00787       nemhelper->read_elem_in_block(i);
00788 
00789       // Note that with parallel files it is possible we have no elements in
00790       // this block!
00791       if (!nemhelper->num_elem_this_blk) continue;
00792 
00793       // Set subdomain ID based on the block ID.
00794       int subdomain_id = nemhelper->block_ids[i];
00795 
00796       // Create a type string (this uses the null-terminated string ctor).
00797       const std::string type_str ( &(nemhelper->elem_type[0]) );
00798 
00799       // Set any relevant node/edge maps for this element
00800       const ExodusII_IO_Helper::Conversion conv = em.assign_conversion(type_str);
00801 
00802       if (_verbose)
00803         libMesh::out << "Reading a block of " << type_str << " elements." << std::endl;
00804 
00805       // Loop over all the elements in this block
00806       for (unsigned int j=0; j<to_uint(nemhelper->num_elem_this_blk); j++)
00807         {
00808           Elem* elem = Elem::build (conv.get_canonical_type()).release();
00809           libmesh_assert (elem);
00810 
00811           // Assign subdomain and processor ID to the newly-created Elem.
00812           // Assigning the processor ID beforehand ensures that the Elem is
00813           // not added as an "unpartitioned" element.  Note that the element
00814           // numbering in Exodus is also 1-based.
00815           elem->subdomain_id() = subdomain_id;
00816           elem->processor_id() = libMesh::processor_id();
00817           elem->set_id()       = my_next_elem++;
00818 
00819           // Mark that we have seen an element of the current element's
00820           // dimension.
00821           elems_of_dimension[elem->dim()] = true;
00822 
00823           // Add the created Elem to the Mesh, catch the Elem
00824           // pointer that the Mesh throws back.
00825           elem = mesh.add_elem (elem);
00826 
00827           // We are expecting the element "thrown back" by libmesh to have the ID we specified for it...
00828           // Check to see that really is the case.  Note that my_next_elem was post-incremented, so
00829           // subtract 1 when performing the check.
00830           if (elem->id() != my_next_elem-1)
00831             {
00832               libMesh::err << "Unexpected ID "
00833                         << elem->id()
00834                         << " set by parallel mesh. (expecting "
00835                         << my_next_elem-1
00836                         << ")." << std::endl;
00837               libmesh_error();
00838             }
00839 
00840           // Set all the nodes for this element
00841           if (_verbose)
00842             libMesh::out << "[" << libMesh::processor_id() << "] "
00843                           << "Setting nodes for Elem " << elem->id() << std::endl;
00844 
00845           for (unsigned int k=0; k<to_uint(nemhelper->num_nodes_per_elem); k++)
00846             {
00847               const unsigned int
00848                 gi              = (j*nemhelper->num_nodes_per_elem +       // index into connectivity array
00849                                    conv.get_node_map(k)),
00850                 local_node_idx  = nemhelper->connect[gi]-1,                // local node index
00851                 global_node_idx = nemhelper->node_num_map[local_node_idx]; // new global node index
00852 
00853               // Set node number
00854               elem->set_node(k) = mesh.node_ptr(global_node_idx);
00855             }
00856         } // for (unsigned int j=0; j<nemhelper->num_elem_this_blk; j++)
00857     } // end for (unsigned int i=0; i<nemhelper->num_elem_blk; i++)
00858 
00859   libmesh_assert_equal_to ((my_next_elem - my_elem_offset), to_uint(nemhelper->num_elem));
00860 
00861   if (_verbose)
00862     {
00863       // Print local elems_of_dimension information
00864       for (unsigned int i=1; i<elems_of_dimension.size(); ++i)
00865         libMesh::out << "[" << libMesh::processor_id() << "] "
00866                      << "elems_of_dimension[" << i << "]=" << elems_of_dimension[i] << std::endl;
00867     }
00868 
00869   // Get the max dimension seen on the current processor
00870   unsigned int max_dim_seen = 0;
00871   for (unsigned int i=1; i<elems_of_dimension.size(); ++i)
00872     if (elems_of_dimension[i])
00873       max_dim_seen = i;
00874 
00875   // Do a global max to determine the max dimension seen by all processors.
00876   // It should match -- I don't think we even support calculations on meshes
00877   // with elements of different dimension...
00878   CommWorld.max(max_dim_seen);
00879 
00880   if (_verbose)
00881     {
00882       // Print the max element dimension from all processors
00883       libMesh::out << "[" << libMesh::processor_id() << "] "
00884                    << "max_dim_seen=" << max_dim_seen << std::endl;
00885     }
00886 
00887   // Set the mesh dimension to the largest encountered for an element
00888   mesh.set_mesh_dimension(max_dim_seen);
00889 
00890 #if LIBMESH_DIM < 3
00891   if (mesh.mesh_dimension() > LIBMESH_DIM)
00892     {
00893       libMesh::err << "Cannot open dimension " <<
00894                       mesh.mesh_dimension() <<
00895                       " mesh file when configured without " <<
00896                       mesh.mesh_dimension() << "D support." <<
00897                       std::endl;
00898       libmesh_error();
00899     }
00900 #endif
00901 
00902 
00903   // Global sideset information, they are distributed as well, not sure if they will require communication...?
00904   nemhelper->get_ss_param_global();
00905 
00906   if (_verbose)
00907     {
00908       libMesh::out << "[" << libMesh::processor_id() << "] "
00909                    << "Read global sideset parameter information." << std::endl;
00910 
00911       // These global values should be the same on all processors...
00912       libMesh::out << "[" << libMesh::processor_id() << "] "
00913                    << "Number of global sideset IDs: " << nemhelper->global_sideset_ids.size() << std::endl;
00914     }
00915 
00916   // Read *local* sideset info the same way it is done in
00917   // exodusII_io_helper.  May be called any time after
00918   // nem_helper->read_header(); This sets num_side_sets and resizes
00919   // elem_list, side_list, and id_list to num_elem_all_sidesets.  Note
00920   // that there appears to be the same number of sidesets in each file
00921   // but they all have different numbers of entries (some are empty).
00922   // Note that the sum of "nemhelper->num_elem_all_sidesets" over all
00923   // processors should equal the sum of the entries in the "num_global_side_counts" array
00924   // filled up by nemhelper->get_ss_param_global()
00925   nemhelper->read_sideset_info();
00926 
00927   if (_verbose)
00928     {
00929       libMesh::out << "[" << libMesh::processor_id() << "] "
00930                    << "nemhelper->num_side_sets = " << nemhelper->num_side_sets << std::endl;
00931 
00932       libMesh::out << "[" << libMesh::processor_id() << "] "
00933                    << "nemhelper->num_elem_all_sidesets = " << nemhelper->num_elem_all_sidesets << std::endl;
00934     }
00935 
00936 #ifdef DEBUG
00937   {
00938     // In DEBUG mode, check that the global number of sidesets reported
00939     // in each nemesis file matches the sum of all local sideset counts
00940     // from each processor.  This requires a small communication, so only
00941     // do it in DEBUG mode.
00942     int sum_num_global_side_counts = std::accumulate(nemhelper->num_global_side_counts.begin(),
00943                                                      nemhelper->num_global_side_counts.end(),
00944                                                      0);
00945 
00946     // MPI sum up the local files contributions
00947     int sum_num_elem_all_sidesets = nemhelper->num_elem_all_sidesets;
00948     CommWorld.sum(sum_num_elem_all_sidesets);
00949 
00950     if (sum_num_global_side_counts != sum_num_elem_all_sidesets)
00951       {
00952         libMesh::err << "Error! global side count reported by Nemesis does not "
00953                      << "match the side count reported by the individual files!" << std::endl;
00954         libmesh_error();
00955       }
00956   }
00957 #endif
00958 
00959   // Note that exodus stores sidesets in separate vectors but we want to pack
00960   // them all into a single vector.  So when we call read_sideset(), we pass an offset
00961   // into the single vector of all IDs
00962   for (int offset=0, i=0; i<nemhelper->get_num_side_sets(); i++)
00963     {
00964       offset += (i > 0 ? nemhelper->get_num_sides_per_set(i-1) : 0); // Compute new offset
00965       nemhelper->read_sideset (i, offset);
00966     }
00967 
00968   // Now that we have the lists of elements, sides, and IDs, we are ready to set them
00969   // in the BoundaryInfo object of our Mesh object.  This is slightly different in parallel...
00970   // For example, I think the IDs in each of the split Exodus files are numbered locally,
00971   // and we have to know the appropriate ID for this processor to be able to set the
00972   // entry in BoundaryInfo.  This offset should be given by my_elem_offset determined in
00973   // this function...
00974 
00975   // Get references to the elem, side, and ID lists for this processor.  We should be
00976   // able to get pointers to these elements on the local processor, and add them to the
00977   // BoundaryInfo object
00978   const std::vector<int>& elem_list = nemhelper->get_elem_list();
00979   const std::vector<int>& side_list = nemhelper->get_side_list();
00980   const std::vector<int>& id_list   = nemhelper->get_id_list();
00981 
00982 
00983   // Debugging:
00984   // Print entries of elem_list
00985   // libMesh::out << "[" << libMesh::processor_id() << "] "
00986   //           << "elem_list = ";
00987   // for (unsigned int e=0; e<elem_list.size(); e++)
00988   //   {
00989   //     libMesh::out << elem_list[e] << ", ";
00990   //   }
00991   // libMesh::out << std::endl;
00992 
00993   // Print entries of side_list
00994   // libMesh::out << "[" << libMesh::processor_id() << "] "
00995   //           << "side_list = ";
00996   // for (unsigned int e=0; e<side_list.size(); e++)
00997   //   {
00998   //     libMesh::out << side_list[e] << ", ";
00999   //   }
01000   // libMesh::out << std::endl;
01001 
01002 
01003   // Loop over the entries of the elem_list, get their pointers from the
01004   // Mesh data structure, and assign the appropriate side to the BoundaryInfo object.
01005   for (unsigned int e=0; e<elem_list.size(); e++)
01006     {
01007       // Calling mesh.elem() feels wrong, for example, in
01008       // ParallelMesh, Mesh::elem() can return NULL so we have to check for this...
01009       //
01010       // Perhaps we should iterate over elements and look them up in
01011       // the elem list instead?  Note that the IDs in elem_list are
01012       // not necessarily in order, so if we did instead loop over the
01013       // mesh, we would have to search the (unsorted) elem_list vector
01014       // for each entry!  We'll settle for doing some error checking instead.
01015       Elem* elem = mesh.elem(my_elem_offset + (elem_list[e]-1)/*Exodus numbering is 1-based!*/);
01016 
01017       if (elem == NULL)
01018         {
01019           libMesh::err << "Mesh returned a NULL pointer when asked for element "
01020                        << my_elem_offset << " + " << elem_list[e] << " = " << my_elem_offset+elem_list[e] << std::endl;
01021           libmesh_error();
01022         }
01023 
01024       // The side numberings in libmesh and exodus are not 1:1, so we need to map
01025       // whatever side number is stored in Exodus into a libmesh side number using
01026       // a conv object...
01027       const ExodusII_IO_Helper::Conversion conv =
01028         em.assign_conversion(elem->type());
01029 
01030       // Finally, we are ready to add the element and its side to the BoundaryInfo object.
01031       // Call the version of add_side which takes a pointer, since we have already gone to
01032       // the trouble of getting said pointer...
01033       mesh.boundary_info->add_side (elem,
01034                                     conv.get_side_map(side_list[e]-1/*Exodus numbering is 1-based*/),
01035                                     id_list[e]);
01036     }
01037 
01038   // Debugging: make sure there are as many boundary conditions in the
01039   // boundary ID object as expected.  Note that, at this point, the
01040   // mesh still thinks it's serial, so n_boundary_conds() returns the
01041   // local number of boundary conditions (and is therefore cheap)
01042   // which should match elem_list.size().
01043   {
01044     std::size_t nbcs = mesh.boundary_info->n_boundary_conds();
01045     if (nbcs != elem_list.size())
01046       {
01047         libMesh::err << "[" << libMesh::processor_id() << "] ";
01048         libMesh::err << "BoundaryInfo contains "
01049                      << nbcs
01050                      << " boundary conditions, while the Exodus file had "
01051                      << elem_list.size()
01052                      << std::endl;
01053         libmesh_error();
01054       }
01055   }
01056 
01057   // Read global nodeset parameters?  We might be able to use this to verify
01058   // something about the local files, but I haven't figured out what yet...
01059   nemhelper->get_ns_param_global();
01060 
01061   // Read local nodeset info
01062   nemhelper->read_nodeset_info();
01063 
01064   if (_verbose)
01065     {
01066       libMesh::out << "[" << libMesh::processor_id() << "] ";
01067       libMesh::out << "nemhelper->num_node_sets=" << nemhelper->num_node_sets << std::endl;
01068     }
01069 
01070 //  // Debugging, what is currently in nemhelper->node_num_map anyway?
01071 //  libMesh::out << "[" << libMesh::processor_id() << "] "
01072 //             << "nemhelper->node_num_map = ";
01073 //
01074 //  for (unsigned int i=0; i<nemhelper->node_num_map.size(); ++i)
01075 //    libMesh::out << nemhelper->node_num_map[i] << ", ";
01076 //  libMesh::out << std::endl;
01077 
01078   // For each nodeset,
01079   for (int nodeset=0; nodeset<nemhelper->get_num_node_sets(); nodeset++)
01080     {
01081       // Get the user-defined ID associcated with the nodeset
01082       int nodeset_id = nemhelper->get_nodeset_id(nodeset);
01083 
01084       if (_verbose)
01085         {
01086           libMesh::out << "[" << libMesh::processor_id() << "] ";
01087           libMesh::out << "nemhelper->get_nodeset_id(" << nodeset << ")=" << nodeset_id << std::endl;
01088         }
01089 
01090       // Read the nodeset from file, store them in a vector
01091       nemhelper->read_nodeset(nodeset);
01092 
01093       // Get access to that vector
01094       const std::vector<int>& node_list = nemhelper->get_node_list();
01095 
01096       // Add nodes from the node_list to the BoundaryInfo object
01097       for(unsigned int node=0; node<node_list.size(); node++)
01098         {
01099           // Don't run past the end of our node map!
01100           if (to_uint(node_list[node]-1) >= nemhelper->node_num_map.size())
01101             {
01102               libMesh::err << "Error, index is past the end of node_num_map array!" << std::endl;
01103               libmesh_error();
01104             }
01105 
01106           // We should be able to use the node_num_map data structure set up previously to determine
01107           // the proper global node index.
01108           unsigned global_node_id = nemhelper->node_num_map[ node_list[node]-1 /*Exodus is 1-based!*/ ];
01109 
01110           if (_verbose)
01111             {
01112               libMesh::out << "[" << libMesh::processor_id() << "] "
01113                            << "nodeset " << nodeset
01114                            << ", local node number: " << node_list[node]-1
01115                            << ", global node id: " << global_node_id
01116                            << std::endl;
01117             }
01118 
01119           // Add the node to the BoundaryInfo object with the proper nodeset_id
01120           mesh.boundary_info->add_node(global_node_id, nodeset_id);
01121         }
01122     }
01123 
01124   // See what the elem count is up to now.
01125   if (_verbose)
01126     {
01127       // Report the number of elements which have been added locally
01128       libMesh::out << "[" << libMesh::processor_id() << "] ";
01129       libMesh::out << "mesh.n_elem()=" << mesh.n_elem() << std::endl;
01130 
01131       // Reports the number of elements that have been added in total.
01132       libMesh::out << "[" << libMesh::processor_id() << "] ";
01133       libMesh::out << "mesh.parallel_n_elem()=" << mesh.parallel_n_elem() << std::endl;
01134     }
01135 
01136   STOP_LOG ("read()","Nemesis_IO");
01137 
01138   // For ParallelMesh, it seems that _is_serial is true by default.  A hack to
01139   // make the Mesh think it's parallel might be to call:
01140   mesh.update_post_partitioning();
01141   mesh.delete_remote_elements();
01142 
01143   // And if that didn't work, then we're actually reading into a
01144   // SerialMesh, so forget about gathering neighboring elements
01145   if (mesh.is_serial())
01146     return;
01147 
01148   // Gather neighboring elements so that the mesh has the proper "ghost" neighbor information.
01149   MeshCommunication().gather_neighboring_elements(libmesh_cast_ref<ParallelMesh&>(mesh));
01150 }

void libMesh::MeshInput< MeshBase >::skip_comment_lines ( std::istream &  in,
const char  comment_start 
) [protected, inherited]

Reads input from in, skipping all the lines that start with the character comment_start.

Referenced by libMesh::TetGenIO::read(), and libMesh::UCDIO::read_implementation().

void libMesh::Nemesis_IO::verbose ( bool  set_verbosity  ) 

Set the flag indicationg if we should be verbose.

Definition at line 102 of file nemesis_io.C.

References _verbose, nemhelper, and libMesh::ExodusII_IO_Helper::verbose().

00103 {
00104   _verbose = set_verbosity;
00105 
00106 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
00107   // Set the verbose flag in the helper object
00108   // as well.
00109   nemhelper->verbose(_verbose);
00110 #endif
00111 }

void libMesh::Nemesis_IO::write ( const std::string &  base_filename  )  [virtual]

This method implements writing a mesh to a specified file.

Implements libMesh::MeshOutput< MeshBase >.

Definition at line 1168 of file nemesis_io.C.

References libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::Nemesis_IO_Helper::create(), libMesh::ExodusII_IO_Helper::ex_err, libMesh::ExodusII_IO_Helper::ex_id, ex_update(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::MeshInput< MeshBase >::mesh(), nemhelper, libMesh::Nemesis_IO_Helper::write_elements(), libMesh::Nemesis_IO_Helper::write_nodal_coordinates(), libMesh::Nemesis_IO_Helper::write_nodesets(), and libMesh::Nemesis_IO_Helper::write_sidesets().

01169 {
01170   // Get a constant reference to the mesh for writing
01171   const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
01172 
01173   // Create the filename for this processor given the base_filename passed in.
01174   std::string nemesis_filename = nemhelper->construct_nemesis_filename(base_filename);
01175 
01176   nemhelper->create(nemesis_filename);
01177 
01178   // Initialize data structures and write some global Nemesis-specifc data, such as
01179   // communication maps, to file.
01180   nemhelper->initialize(nemesis_filename,mesh);
01181 
01182   // Call the Nemesis-specialized version of write_nodal_coordinates() to write
01183   // the nodal coordinates.
01184   nemhelper->write_nodal_coordinates(mesh);
01185 
01186   // Call the Nemesis-specialized version of write_elements() to write
01187   // the elements.  Note: Must write a zero if a given global block ID has no
01188   // elements...
01189   nemhelper->write_elements(mesh);
01190 
01191   // Call our specialized function to write the nodesets
01192   nemhelper->write_nodesets(mesh);
01193 
01194   // Call our specialized write_sidesets() function to write the sidesets to file
01195   nemhelper->write_sidesets(mesh);
01196 
01197   // Not sure if this is really necessary, but go ahead and flush the file
01198   // once we have written all this stuff.
01199   nemhelper->ex_err = exII::ex_update(nemhelper->ex_id);
01200 }

virtual void libMesh::MeshOutput< MeshBase >::write_equation_systems ( const std::string &  ,
const EquationSystems ,
const std::set< std::string > *  system_names = NULL 
) [virtual, inherited]

This method implements writing a mesh with data to a specified file where the data is taken from the EquationSystems object.

Referenced by write_timestep().

void libMesh::Nemesis_IO::write_global_data ( const std::vector< Number > &  soln,
const std::vector< std::string > &  names 
)

Write out global variables.

Definition at line 1292 of file nemesis_io.C.

References _timestep, libMesh::ExodusII_IO_Helper::created(), libMesh::err, libMesh::ExodusII_IO_Helper::initialize_global_variables(), nemhelper, and libMesh::ExodusII_IO_Helper::write_global_values().

01294 {
01295   if (!nemhelper->created())
01296     {
01297       libMesh::err << "ERROR, Nemesis file must be initialized "
01298                    << "before outputting global variables.\n"
01299                    << std::endl;
01300       libmesh_error();
01301     }
01302 
01303   // Call the Exodus writer implementation
01304   nemhelper->initialize_global_variables( names );
01305   nemhelper->write_global_values( soln, _timestep);
01306 }

void libMesh::Nemesis_IO::write_information_records ( const std::vector< std::string > &  records  ) 

Write out information records.

Definition at line 1324 of file nemesis_io.C.

References libMesh::ExodusII_IO_Helper::created(), libMesh::err, nemhelper, and libMesh::ExodusII_IO_Helper::write_information_records().

01325 {
01326   if (!nemhelper->created())
01327     {
01328       libMesh::err << "ERROR, Nemesis file must be initialized "
01329                    << "before outputting information records.\n"
01330                    << std::endl;
01331       libmesh_error();
01332     }
01333 
01334   // Call the Exodus writer implementation
01335   nemhelper->write_information_records( records );
01336 }

void libMesh::Nemesis_IO::write_nodal_data ( const std::string &  fname,
const std::vector< Number > &  soln,
const std::vector< std::string > &  names 
) [virtual]

Output a nodal solution.

Reimplemented from libMesh::MeshOutput< MeshBase >.

Definition at line 1241 of file nemesis_io.C.

References _timestep, libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::Nemesis_IO_Helper::create(), libMesh::ExodusII_IO_Helper::created(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::ExodusII_IO_Helper::initialize_nodal_variables(), libMesh::MeshInput< MeshBase >::mesh(), nemhelper, libMesh::Nemesis_IO_Helper::write_elements(), libMesh::Nemesis_IO_Helper::write_nodal_coordinates(), libMesh::Nemesis_IO_Helper::write_nodal_solution(), libMesh::Nemesis_IO_Helper::write_nodesets(), and libMesh::Nemesis_IO_Helper::write_sidesets().

01244 {
01245   START_LOG("write_nodal_data()", "Nemesis_IO");
01246 
01247   const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
01248 
01249   std::string nemesis_filename = nemhelper->construct_nemesis_filename(base_filename);
01250 
01251   if (!nemhelper->created())
01252   {
01253     nemhelper->create(nemesis_filename);
01254     nemhelper->initialize(nemesis_filename,mesh);
01255     nemhelper->write_nodal_coordinates(mesh);
01256     nemhelper->write_elements(mesh);
01257     nemhelper->write_nodesets(mesh);
01258     nemhelper->write_sidesets(mesh);
01259 
01260     // If we don't have any nodes written out on this processor,
01261     // Exodus seems to like us better if we don't try to write out any
01262     // variable names too...
01263     nemhelper->initialize_nodal_variables(names);
01264   }
01265 
01266   nemhelper->write_nodal_solution(soln, names, _timestep);
01267 
01268   STOP_LOG("write_nodal_data()", "Nemesis_IO");
01269 }

void libMesh::Nemesis_IO::write_timestep ( const std::string &  fname,
const EquationSystems es,
const int  timestep,
const Real  time 
)

Write one timestep's worth of the solution.

Definition at line 1215 of file nemesis_io.C.

References _timestep, nemhelper, libMesh::MeshOutput< MeshBase >::write_equation_systems(), and libMesh::ExodusII_IO_Helper::write_timestep().

01219 {
01220   _timestep=timestep;
01221   write_equation_systems(fname,es);
01222 
01223   nemhelper->write_timestep(timestep, time);
01224 }


Member Data Documentation

const bool libMesh::MeshOutput< MeshBase >::_is_parallel_format [protected, inherited]

Flag specifying whether this format is parallel-capable. If this is false (default) I/O is only permitted when the mesh has been serialized.

Definition at line 126 of file mesh_output.h.

Referenced by libMesh::PostscriptIO::write(), libMesh::FroIO::write(), libMesh::EnsightIO::write(), and libMesh::DivaIO::write().

Definition at line 112 of file nemesis_io.h.

Referenced by write_global_data(), write_nodal_data(), and write_timestep().

Definition at line 114 of file nemesis_io.h.

Referenced by read(), and verbose().


The documentation for this class was generated from the following files:

Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:30 UTC

Hosted By:
SourceForge.net Logo