libMesh
nemesis_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 // LibMesh includes
20 #include "libmesh/distributed_mesh.h"
21 #include "libmesh/dof_map.h" // local_index
22 #include "libmesh/elem.h"
23 #include "libmesh/exodusII_io.h"
24 #include "libmesh/libmesh_logging.h"
25 #include "libmesh/nemesis_io.h"
26 #include "libmesh/nemesis_io_helper.h"
27 #include "libmesh/node.h"
28 #include "libmesh/parallel.h"
29 #include "libmesh/utility.h" // deallocate
30 #include "libmesh/boundary_info.h"
31 #include "libmesh/mesh_communication.h"
32 #include "libmesh/fe_type.h"
33 #include "libmesh/equation_systems.h"
34 #include "libmesh/numeric_vector.h"
35 #include "libmesh/int_range.h"
36 
37 // C++ includes
38 #include <memory>
39 #include <numeric> // std::accumulate
40 
41 namespace libMesh
42 {
43 
44 
45 //-----------------------------------------------
46 // anonymous namespace for implementation details
47 namespace {
48 
49 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
50 struct CompareGlobalIdxMappings
51 {
52  // strict weak ordering for a.first -> a.second mapping. since we can only map to one
53  // value only order the first entry
54  bool operator()(const std::pair<unsigned int, unsigned int> & a,
55  const std::pair<unsigned int, unsigned int> & b) const
56  { return a.first < b.first; }
57 
58  // strict weak ordering for a.first -> a.second mapping. lookups will
59  // be in terms of a single integer, which is why we need this method.
60  bool operator()(const std::pair<unsigned int, unsigned int> & a,
61  const unsigned int b) const
62  { return a.first < b; }
63 };
64 #endif // defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
65 
66 // Nemesis & ExodusII use int for all integer values, even the ones which
67 // should never be negative. we like to use unsigned as a force of habit,
68 // this trivial little method saves some typing & also makes sure something
69 // is not horribly wrong.
70 template <typename T>
71 inline unsigned int to_uint ( const T & t )
72 {
73  libmesh_assert_equal_to (t, static_cast<T>(static_cast<unsigned int>(t)));
74 
75  return static_cast<unsigned int>(t);
76 }
77 
78 // test equality for a.first -> a.second mapping. since we can only map to one
79 // value only test the first entry
80 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API) && !defined(NDEBUG)
81 inline bool global_idx_mapping_equality (const std::pair<unsigned int, unsigned int> & a,
82  const std::pair<unsigned int, unsigned int> & b)
83 {
84  return a.first == b.first;
85 }
86 #endif
87 
88 }
89 
90 
91 
92 // ------------------------------------------------------------
93 // Nemesis_IO class members
95  bool single_precision) :
96  MeshInput<MeshBase> (mesh, /*is_parallel_format=*/true),
97  MeshOutput<MeshBase> (mesh, /*is_parallel_format=*/true),
99 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
100  nemhelper(std::make_unique<Nemesis_IO_Helper>(*this, false, single_precision)),
101  _timestep(1),
102 #endif
103  _verbose (false),
104  _append(false),
105  _allow_empty_variables(false)
106 {
107  // if !LIBMESH_HAVE_EXODUS_API, we didn't use this
108  libmesh_ignore(single_precision);
109 }
110 
111 
112 
114  bool single_precision) :
115  MeshInput<MeshBase> (),
116  MeshOutput<MeshBase> (mesh, /*is_parallel_format=*/true),
118 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
119  nemhelper(std::make_unique<Nemesis_IO_Helper>(*this, false, single_precision)),
120  _timestep(1),
121 #endif
122  _verbose (false),
123  _append(false),
124  _allow_empty_variables(false)
125 {
126  // if !LIBMESH_HAVE_EXODUS_API, we didn't use this
127  libmesh_ignore(single_precision);
128 }
129 
130 
131 
132 // Destructor. Defined in the C file so we can be sure to get away
133 // with a forward declaration of Nemesis_IO_Helper in the header file.
134 Nemesis_IO::~Nemesis_IO () = default;
135 
136 
137 
138 void Nemesis_IO::verbose (bool set_verbosity)
139 {
140  _verbose = set_verbosity;
141 
142 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
143  // Set the verbose flag in the helper object
144  // as well.
145  nemhelper->verbose = _verbose;
146 #endif
147 }
148 
149 
150 
152 {
153 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
154  nemhelper->write_complex_abs = val;
155 #endif
156  libmesh_ignore(val);
157 }
158 
159 
160 
161 void Nemesis_IO::append(bool val)
162 {
163  _append = val;
164 }
165 
166 
167 
168 void Nemesis_IO::set_output_variables(const std::vector<std::string> & output_variables,
169  bool allow_empty)
170 {
171  _output_variables = output_variables;
172  _allow_empty_variables = allow_empty;
173 }
174 
175 
176 
178 {
179 #ifndef NDEBUG
180 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
181  // We expect the communication maps to be symmetric - e.g. if processor i thinks it
182  // communicates with processor j, then processor j should also be expecting to
183  // communicate with i. We can assert that here easily enough with an alltoall,
184  // but let's only do it when not in optimized mode to limit unnecessary communication.
185  {
186  std::vector<unsigned char> pid_send_partner (this->n_processors(), 0);
187 
188  // strictly speaking, we should expect to communicate with ourself...
189  pid_send_partner[this->processor_id()] = 1;
190 
191  // mark each processor id we reference with a node cmap
192  for (unsigned int cmap=0; cmap<to_uint(nemhelper->num_node_cmaps); cmap++)
193  {
194  libmesh_assert_less (nemhelper->node_cmap_ids[cmap], this->n_processors());
195 
196  pid_send_partner[nemhelper->node_cmap_ids[cmap]] = 1;
197  }
198 
199  // Copy the send pairing so we can catch the receive paring and
200  // test for equality
201  const std::vector<unsigned char> pid_recv_partner (pid_send_partner);
202 
203  this->comm().alltoall (pid_send_partner);
204 
205  libmesh_assert (pid_send_partner == pid_recv_partner);
206  }
207 #endif
208 #endif
209 }
210 
211 
212 
213 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
214 void Nemesis_IO::read (const std::string & base_filename)
215 {
216  LOG_SCOPE ("read()","Nemesis_IO");
217 
218  // This function must be run on all processors at once
219  parallel_object_only();
220 
221  if (_verbose)
222  {
223  libMesh::out << "[" << this->processor_id() << "] ";
224  libMesh::out << "Reading Nemesis file on processor: " << this->processor_id() << std::endl;
225  }
226 
227  // Construct the Nemesis filename based on the number of processors and the
228  // current processor ID.
229  std::string nemesis_filename = nemhelper->construct_nemesis_filename(base_filename);
230 
231  if (_verbose)
232  libMesh::out << "Opening file: " << nemesis_filename << std::endl;
233 
234  // Open the Exodus file in EX_READ mode
235  nemhelper->open(nemesis_filename.c_str(), /*read_only=*/true);
236 
237  // Get a reference to the mesh. We need to be specific
238  // since Nemesis_IO is multiply-inherited
239  // MeshBase & mesh = this->mesh();
241 
242  // We're reading a file on each processor, so our mesh is
243  // partitioned into that many parts as it's created
244  this->set_n_partitions(this->n_processors());
245 
246  // Local information: Read the following information from the standard Exodus header
247  // title[0]
248  // num_dim
249  // num_nodes
250  // num_elem
251  // num_elem_blk
252  // num_node_sets
253  // num_side_sets
254  nemhelper->read_and_store_header_info();
255  nemhelper->print_header();
256 
257  // Get global information: number of nodes, elems, blocks, nodesets and sidesets
258  nemhelper->get_init_global();
259 
260  // Get "load balance" information. This includes the number of internal & border
261  // nodes and elements as well as the number of communication maps.
262  nemhelper->get_loadbal_param();
263 
264  // Do some error checking
265  libmesh_error_msg_if(nemhelper->num_external_nodes,
266  "ERROR: there should be no external nodes in an element-based partitioning!");
267 
268  libmesh_assert_equal_to (nemhelper->num_nodes,
269  (nemhelper->num_internal_nodes +
270  nemhelper->num_border_nodes));
271 
272  libmesh_assert_equal_to (nemhelper->num_elem,
273  (nemhelper->num_internal_elems +
274  nemhelper->num_border_elems));
275 
276  libmesh_assert_less_equal (nemhelper->num_nodes, nemhelper->num_nodes_global);
277  libmesh_assert_less_equal (nemhelper->num_elem, nemhelper->num_elems_global);
278 
279  // Read nodes from the exodus file: this fills the nemhelper->x,y,z arrays.
280  nemhelper->read_nodes();
281 
282  // Reads the nemhelper->node_num_map array, node_num_map[i] is the global node number for
283  // local node number i.
284  nemhelper->read_node_num_map();
285 
286  // The get_cmap_params() function reads in the:
287  // node_cmap_ids[],
288  // node_cmap_node_cnts[],
289  // elem_cmap_ids[],
290  // elem_cmap_elem_cnts[],
291  nemhelper->get_cmap_params();
292 
293  // Read the IDs of the interior, boundary, and external nodes. This function
294  // fills the vectors:
295  // node_mapi[],
296  // node_mapb[],
297  // node_mape[]
298  nemhelper->get_node_map();
299 
300  // Read each node communication map for this processor. This function
301  // fills the vectors of vectors named:
302  // node_cmap_node_ids[][]
303  // node_cmap_proc_ids[][]
304  nemhelper->get_node_cmap();
305 
306  libmesh_assert_equal_to (to_uint(nemhelper->num_node_cmaps), nemhelper->node_cmap_node_cnts.size());
307  libmesh_assert_equal_to (to_uint(nemhelper->num_node_cmaps), nemhelper->node_cmap_node_ids.size());
308  libmesh_assert_equal_to (to_uint(nemhelper->num_node_cmaps), nemhelper->node_cmap_proc_ids.size());
309 
310  this->assert_symmetric_cmaps();
311 
312  // We now have enough information to infer node ownership. We start by assuming
313  // we own all the nodes on this processor. We will then interrogate the
314  // node cmaps and see if a lower-rank processor is associated with any of
315  // our nodes. If so, then that processor owns the node, not us...
316  std::vector<processor_id_type> node_ownership (nemhelper->num_internal_nodes +
317  nemhelper->num_border_nodes,
318  this->processor_id());
319 
320  // a map from processor id to cmap number, to be used later
321  std::map<unsigned int, unsigned int> pid_to_cmap_map;
322 
323  // For each node_cmap...
324  for (unsigned int cmap=0; cmap<to_uint(nemhelper->num_node_cmaps); cmap++)
325  {
326  // Good time for error checking...
327  libmesh_assert_equal_to (to_uint(nemhelper->node_cmap_node_cnts[cmap]),
328  nemhelper->node_cmap_node_ids[cmap].size());
329 
330  libmesh_assert_equal_to (to_uint(nemhelper->node_cmap_node_cnts[cmap]),
331  nemhelper->node_cmap_proc_ids[cmap].size());
332 
333  // In all the samples I have seen, node_cmap_ids[cmap] is the processor
334  // rank of the remote processor...
335  const processor_id_type adjcnt_pid_idx =
336  cast_int<processor_id_type>(nemhelper->node_cmap_ids[cmap]);
337 
338  libmesh_assert_less (adjcnt_pid_idx, this->n_processors());
339  libmesh_assert_not_equal_to (adjcnt_pid_idx, this->processor_id());
340 
341  // We only expect one cmap per adjacent processor
342  libmesh_assert (!pid_to_cmap_map.count(adjcnt_pid_idx));
343 
344  pid_to_cmap_map[adjcnt_pid_idx] = cmap;
345 
346  // ...and each node in that cmap...
347  for (unsigned int idx=0; idx<to_uint(nemhelper->node_cmap_node_cnts[cmap]); idx++)
348  {
349  // Are the node_cmap_ids and node_cmap_proc_ids really redundant?
350  libmesh_assert_equal_to
351  (adjcnt_pid_idx,
352  cast_int<processor_id_type>(nemhelper->node_cmap_proc_ids[cmap][idx]));
353 
354  // we are expecting the exodus node numbering to be 1-based...
355  const unsigned int local_node_idx = nemhelper->node_cmap_node_ids[cmap][idx]-1;
356 
357  libmesh_assert_less (local_node_idx, node_ownership.size());
358 
359  // if the adjacent processor is lower rank than the current
360  // owner for this node, then it will get the node...
361  node_ownership[local_node_idx] =
362  std::min(node_ownership[local_node_idx], adjcnt_pid_idx);
363  }
364  } // We now should have established proper node ownership.
365 
366  // now that ownership is established, we can figure out how many nodes we
367  // will be responsible for numbering.
368  unsigned int num_nodes_i_must_number = 0;
369 
370  for (const auto & pid : node_ownership)
371  if (pid == this->processor_id())
372  num_nodes_i_must_number++;
373 
374  // more error checking...
375  libmesh_assert_greater_equal (num_nodes_i_must_number, nemhelper->num_internal_nodes);
376  libmesh_assert (num_nodes_i_must_number <= to_uint(nemhelper->num_internal_nodes +
377  nemhelper->num_border_nodes));
378  if (_verbose)
379  libMesh::out << "[" << this->processor_id() << "] "
380  << "num_nodes_i_must_number="
381  << num_nodes_i_must_number
382  << std::endl;
383 
384  // The call to get_loadbal_param() gets 7 pieces of information. We allgather
385  // these now across all processors to determine some global numberings. We should
386  // also gather the number of nodes each processor thinks it will number so that
387  // we can (i) determine our offset, and (ii) do some error checking.
388  std::vector<int> all_loadbal_data ( 8 );
389  all_loadbal_data[0] = nemhelper->num_internal_nodes;
390  all_loadbal_data[1] = nemhelper->num_border_nodes;
391  all_loadbal_data[2] = nemhelper->num_external_nodes;
392  all_loadbal_data[3] = nemhelper->num_internal_elems;
393  all_loadbal_data[4] = nemhelper->num_border_elems;
394  all_loadbal_data[5] = nemhelper->num_node_cmaps;
395  all_loadbal_data[6] = nemhelper->num_elem_cmaps;
396  all_loadbal_data[7] = num_nodes_i_must_number;
397 
398  this->comm().allgather (all_loadbal_data, /* identical_buffer_sizes = */ true);
399 
400  // OK, we are now in a position to request new global indices for all the nodes
401  // we do not own
402 
403  // Let's get a unique message tag to use for send()/receive()
405 
406  std::vector<std::vector<int>>
407  needed_node_idxs (nemhelper->num_node_cmaps); // the indices we will ask for
408 
409  std::vector<Parallel::Request>
410  needed_nodes_requests (nemhelper->num_node_cmaps);
411 
412  for (unsigned int cmap=0; cmap<to_uint(nemhelper->num_node_cmaps); cmap++)
413  {
414  // We know we will need no more indices than there are nodes
415  // in this cmap, but that number is an upper bound in general
416  // since the neighboring processor associated with the cmap
417  // may not actually own it
418  needed_node_idxs[cmap].reserve (nemhelper->node_cmap_node_cnts[cmap]);
419 
420  const unsigned int adjcnt_pid_idx = nemhelper->node_cmap_ids[cmap];
421 
422  // ...and each node in that cmap...
423  for (unsigned int idx=0; idx<to_uint(nemhelper->node_cmap_node_cnts[cmap]); idx++)
424  {
425  const unsigned int
426  local_node_idx = nemhelper->node_cmap_node_ids[cmap][idx]-1,
427  owning_pid_idx = node_ownership[local_node_idx];
428 
429  // add it to the request list for its owning processor.
430  if (owning_pid_idx == adjcnt_pid_idx)
431  {
432  const unsigned int
433  global_node_idx = nemhelper->node_num_map[local_node_idx]-1;
434  needed_node_idxs[cmap].push_back(global_node_idx);
435  }
436  }
437  // now post the send for this cmap
438  this->comm().send (adjcnt_pid_idx, // destination
439  needed_node_idxs[cmap], // send buffer
440  needed_nodes_requests[cmap], // request
441  nodes_tag);
442  } // all communication requests for getting updated global indices for border
443  // nodes have been initiated
444 
445  // Figure out how many nodes each processor thinks it will number and make sure
446  // that it adds up to the global number of nodes. Also, set up global node
447  // index offsets for each processor.
448  std::vector<unsigned int>
449  all_num_nodes_i_must_number (this->n_processors());
450 
451  for (auto pid : make_range(this->n_processors()))
452  all_num_nodes_i_must_number[pid] = all_loadbal_data[8*pid + 7];
453 
454  // The sum of all the entries in this vector should sum to the number of global nodes
455  libmesh_assert (std::accumulate(all_num_nodes_i_must_number.begin(),
456  all_num_nodes_i_must_number.end(),
457  0) == nemhelper->num_nodes_global);
458 
459  unsigned int my_next_node = 0;
460  for (auto pid : make_range(this->processor_id()))
461  my_next_node += all_num_nodes_i_must_number[pid];
462 
463  const unsigned int my_node_offset = my_next_node;
464 
465  if (_verbose)
466  libMesh::out << "[" << this->processor_id() << "] "
467  << "my_node_offset="
468  << my_node_offset
469  << std::endl;
470 
471  // Add internal nodes to the DistributedMesh, using the node ID offset we
472  // computed and the current processor's ID.
473  for (unsigned int i=0; i<to_uint(nemhelper->num_internal_nodes); ++i)
474  {
475  const unsigned int local_node_idx = nemhelper->node_mapi[i]-1;
476 #ifndef NDEBUG
477  const unsigned int owning_pid_idx = node_ownership[local_node_idx];
478 #endif
479 
480  // an internal node we do not own? huh??
481  libmesh_assert_equal_to (owning_pid_idx, this->processor_id());
482  libmesh_assert_less (my_next_node, nemhelper->num_nodes_global);
483 
484  // "Catch" the node pointer after addition, make sure the
485  // ID matches the requested value.
486  Node * added_node =
487  mesh.add_point (Point(nemhelper->x[local_node_idx],
488  nemhelper->y[local_node_idx],
489  nemhelper->z[local_node_idx]),
490  my_next_node,
491  this->processor_id());
492 
493  // Make sure the node we added has the ID we thought it would
494  if (added_node->id() != my_next_node)
495  {
496  libMesh::err << "Error, node added with ID " << added_node->id()
497  << ", but we wanted ID " << my_next_node << std::endl;
498  }
499 
500  // Set a unique_id ourselves since ReplicatedMesh can't handle
501  // distributed unique_id generation. Make sure it doesn't
502  // overlap element unique_id() values either.
503 #ifdef LIBMESH_ENABLE_UNIQUE_ID
504  added_node->set_unique_id(added_node->id() + nemhelper->num_elems_global);
505 #endif
506 
507  // update the local->global index map, keeping it 1-based
508  nemhelper->node_num_map[local_node_idx] = my_next_node++ + 1;
509  }
510 
511  // Now, for the boundary nodes... We may very well own some of them,
512  // but there may be others for which we have requested the new global
513  // id. We expect to be asked for the ids of the ones we own, so
514  // we need to create a map from the old global id to the new one
515  // we are about to create.
516  typedef std::vector<std::pair<unsigned int, unsigned int>> global_idx_mapping_type;
517  global_idx_mapping_type old_global_to_new_global_map;
518  old_global_to_new_global_map.reserve (num_nodes_i_must_number // total # i will have
519  - (my_next_node // amount i have thus far
520  - my_node_offset)); // this should be exact!
521  CompareGlobalIdxMappings global_idx_mapping_comp;
522 
523  for (unsigned int i=0; i<to_uint(nemhelper->num_border_nodes); ++i)
524  {
525  const unsigned int
526  local_node_idx = nemhelper->node_mapb[i]-1,
527  owning_pid_idx = node_ownership[local_node_idx];
528 
529  // if we own it...
530  if (owning_pid_idx == this->processor_id())
531  {
532  const unsigned int
533  global_node_idx = nemhelper->node_num_map[local_node_idx]-1;
534 
535  // we will number it, and create a mapping from its old global index to
536  // the new global index, for lookup purposes when neighbors come calling
537  old_global_to_new_global_map.emplace_back(global_node_idx, my_next_node);
538 
539  // "Catch" the node pointer after addition, make sure the
540  // ID matches the requested value.
541  Node * added_node =
542  mesh.add_point (Point(nemhelper->x[local_node_idx],
543  nemhelper->y[local_node_idx],
544  nemhelper->z[local_node_idx]),
545  my_next_node,
546  this->processor_id());
547 
548  // Make sure the node we added has the ID we thought it would
549  if (added_node->id() != my_next_node)
550  {
551  libMesh::err << "Error, node added with ID " << added_node->id()
552  << ", but we wanted ID " << my_next_node << std::endl;
553  }
554 
555  // Set a unique_id ourselves since ReplicatedMesh can't handle
556  // distributed unique_id generation. Make sure it doesn't
557  // overlap element unique_id() values either.
558 #ifdef LIBMESH_ENABLE_UNIQUE_ID
559  added_node->set_unique_id(added_node->id() + nemhelper->num_elems_global);
560 #endif
561 
562  // update the local->global index map, keeping it 1-based
563  nemhelper->node_num_map[local_node_idx] = my_next_node++ + 1;
564  }
565  }
566  // That should cover numbering all the nodes which belong to us...
567  libmesh_assert_equal_to (num_nodes_i_must_number, (my_next_node - my_node_offset));
568 
569  // Let's sort the mapping so we can efficiently answer requests
570  std::sort (old_global_to_new_global_map.begin(),
571  old_global_to_new_global_map.end(),
572  global_idx_mapping_comp);
573 
574  // and it had better be unique...
575  libmesh_assert (std::unique (old_global_to_new_global_map.begin(),
576  old_global_to_new_global_map.end(),
577  global_idx_mapping_equality) ==
578  old_global_to_new_global_map.end());
579 
580  // We can now catch incoming requests and process them. for efficiency
581  // let's do whatever is available next
582  std::map<unsigned int, std::vector<int>> requested_node_idxs; // the indices asked of us
583 
584  std::vector<Parallel::Request> requested_nodes_requests(nemhelper->num_node_cmaps);
585 
586  // We know we will receive the request from a given processor before
587  // we receive its reply to our request. However, we may receive
588  // a request and a response from one processor before getting
589  // a request from another processor. So what we are doing here
590  // is processing whatever message comes next, while recognizing
591  // we will receive a request from a processor before receiving
592  // its reply
593  std::vector<bool> processed_cmap (nemhelper->num_node_cmaps, false);
594 
595  for (unsigned int comm_step=0; comm_step<2*to_uint(nemhelper->num_node_cmaps); comm_step++)
596  {
597  // query the first message which is available
598  const Parallel::Status
599  status (this->comm().probe (Parallel::any_source,
600  nodes_tag));
601  const unsigned int
602  requesting_pid_idx = status.source(),
603  source_pid_idx = status.source();
604 
605  // this had better be from a processor we are expecting...
606  libmesh_assert (pid_to_cmap_map.count(requesting_pid_idx));
607 
608  // the local cmap which corresponds to the source processor
609  const unsigned int cmap = pid_to_cmap_map[source_pid_idx];
610 
611  if (!processed_cmap[cmap])
612  {
613  processed_cmap[cmap] = true;
614 
615  // we should only get one request per paired processor
616  libmesh_assert (!requested_node_idxs.count(requesting_pid_idx));
617 
618  // get a reference to the request buffer for this processor to
619  // avoid repeated map lookups
620  std::vector<int> & xfer_buf (requested_node_idxs[requesting_pid_idx]);
621 
622  // actually receive the message.
623  this->comm().receive (requesting_pid_idx, xfer_buf, nodes_tag);
624 
625  // Fill the request
626  for (auto i : index_range(xfer_buf))
627  {
628  // the requested old global node index, *now 0-based*
629  const unsigned int old_global_node_idx = xfer_buf[i];
630 
631  // find the new global node index for the requested node -
632  // note that requesting_pid_idx thinks we own this node,
633  // so we better!
634  const global_idx_mapping_type::const_iterator it =
635  std::lower_bound (old_global_to_new_global_map.begin(),
636  old_global_to_new_global_map.end(),
637  old_global_node_idx,
638  global_idx_mapping_comp);
639 
640  libmesh_assert (it != old_global_to_new_global_map.end());
641  libmesh_assert_equal_to (it->first, old_global_node_idx);
642  libmesh_assert_greater_equal (it->second, my_node_offset);
643  libmesh_assert_less (it->second, my_next_node);
644 
645  // overwrite the requested old global node index with the new global index
646  xfer_buf[i] = it->second;
647  }
648 
649  // and send the new global indices back to the processor which asked for them
650  this->comm().send (requesting_pid_idx,
651  xfer_buf,
652  requested_nodes_requests[cmap],
653  nodes_tag);
654  } // done processing the request
655 
656  // this is the second time we have heard from this processor,
657  // so it must be its reply to our request
658  else
659  {
660  // a long time ago, we sent off our own requests. now it is time to catch the
661  // replies and get the new global node numbering. note that for any reply
662  // we receive, the corresponding nonblocking send from above *must* have been
663  // completed, since the reply is in response to that request!!
664 
665  // if we have received a reply, our send *must* have completed
666  // (note we never actually need to wait on the request)
667  libmesh_assert (needed_nodes_requests[cmap].test());
668  libmesh_assert_equal_to (to_uint(nemhelper->node_cmap_ids[cmap]), source_pid_idx);
669 
670  // now post the receive for this cmap
671  this->comm().receive (source_pid_idx,
672  needed_node_idxs[cmap],
673  nodes_tag);
674 
675  libmesh_assert_less_equal (needed_node_idxs[cmap].size(),
676  nemhelper->node_cmap_node_ids[cmap].size());
677 
678  for (std::size_t i=0, j=0, ncnis=nemhelper->node_cmap_node_ids[cmap].size(); i < ncnis; i++)
679  {
680  const unsigned int
681  local_node_idx = nemhelper->node_cmap_node_ids[cmap][i]-1,
682  owning_pid_idx = node_ownership[local_node_idx];
683 
684  // if this node is owned by source_pid_idx, its new global id
685  // is in the buffer we just received
686  if (owning_pid_idx == source_pid_idx)
687  {
688  libmesh_assert_less (j, needed_node_idxs[cmap].size());
689 
690  const unsigned int // now 0-based!
691  global_node_idx = needed_node_idxs[cmap][j++];
692 
693  // "Catch" the node pointer after addition, make sure the
694  // ID matches the requested value.
695  Node * added_node =
696  mesh.add_point (Point(nemhelper->x[local_node_idx],
697  nemhelper->y[local_node_idx],
698  nemhelper->z[local_node_idx]),
699  cast_int<dof_id_type>(global_node_idx),
700  cast_int<processor_id_type>(source_pid_idx));
701 
702  // Make sure the node we added has the ID we thought it would
703  if (added_node->id() != global_node_idx)
704  {
705  libMesh::err << "Error, node added with ID " << added_node->id()
706  << ", but we wanted ID " << global_node_idx << std::endl;
707  }
708 
709  // Set a unique_id ourselves since ReplicatedMesh can't handle
710  // distributed unique_id generation. Make sure it doesn't
711  // overlap element unique_id() values either.
712 #ifdef LIBMESH_ENABLE_UNIQUE_ID
713  added_node->set_unique_id(added_node->id() + nemhelper->num_elems_global);
714 #endif
715 
716  // update the local->global index map, keeping it 1-based
717  nemhelper->node_num_map[local_node_idx] = global_node_idx + 1;
718 
719  // we are not really going to use my_next_node again, but we can
720  // keep incrementing it to track how many nodes we have added
721  // to the mesh
722  my_next_node++;
723  }
724  }
725  }
726  } // end of node index communication loop
727 
728  // we had better have added all the nodes we need to!
729  libmesh_assert_equal_to ((my_next_node - my_node_offset), to_uint(nemhelper->num_nodes));
730 
731  // After all that, we should be done with all node-related arrays
732  // *except* the node_num_map.
733  // So let's clean up the arrays we are done with.
734  {
735  Utility::deallocate (nemhelper->node_mapi);
736  Utility::deallocate (nemhelper->node_mapb);
737  Utility::deallocate (nemhelper->node_mape);
738  Utility::deallocate (nemhelper->node_cmap_ids);
739  Utility::deallocate (nemhelper->node_cmap_node_cnts);
740  Utility::deallocate (nemhelper->node_cmap_node_ids);
741  Utility::deallocate (nemhelper->node_cmap_proc_ids);
745  Utility::deallocate (needed_node_idxs);
746  Utility::deallocate (node_ownership);
747  }
748 
749  Parallel::wait (needed_nodes_requests);
750  Parallel::wait (requested_nodes_requests);
751  requested_node_idxs.clear();
752 
753  // See what the node count is up to now.
754  if (_verbose)
755  {
756  // Report the number of nodes which have been added locally
757  libMesh::out << "[" << this->processor_id() << "] ";
758  libMesh::out << "mesh.n_nodes()=" << mesh.n_nodes() << std::endl;
759 
760  // Reports the number of nodes that have been added in total.
761  libMesh::out << "[" << this->processor_id() << "] ";
762  libMesh::out << "mesh.parallel_n_nodes()=" << mesh.parallel_n_nodes() << std::endl;
763  }
764 
765 
766 
767  // --------------------------------------------------------------------------------
768  // --------------------------------------------------------------------------------
769  // --------------------------------------------------------------------------------
770 
771 
772  // We can now read in the elements...Exodus stores them in blocks in which all
773  // elements have the same geometric type. This code is adapted directly from exodusII_io.C
774 
775  // Assertion: The sum of the border and internal elements on all processors
776  // should equal nemhelper->num_elems_global
777 #ifndef NDEBUG
778  {
779  int sum_internal_elems=0, sum_border_elems=0;
780  for (unsigned int j=3,c=0; c<this->n_processors(); j+=8,++c)
781  sum_internal_elems += all_loadbal_data[j];
782 
783  for (unsigned int j=4,c=0; c<this->n_processors(); j+=8,++c)
784  sum_border_elems += all_loadbal_data[j];
785 
786  if (_verbose)
787  {
788  libMesh::out << "[" << this->processor_id() << "] ";
789  libMesh::out << "sum_internal_elems=" << sum_internal_elems << std::endl;
790 
791  libMesh::out << "[" << this->processor_id() << "] ";
792  libMesh::out << "sum_border_elems=" << sum_border_elems << std::endl;
793  }
794 
795  libmesh_assert_equal_to (sum_internal_elems+sum_border_elems, nemhelper->num_elems_global);
796  }
797 #endif
798 
799  // We need to set the mesh dimension, but the following...
800  // mesh.set_mesh_dimension(static_cast<unsigned int>(nemhelper->num_dim));
801 
802  // ... is not sufficient since some codes report num_dim==3 for two dimensional
803  // meshes living in 3D, even though all the elements are of 2D type. Therefore,
804  // we instead use the dimension of the highest element found for the Mesh dimension,
805  // similar to what is done by the Exodus reader, except here it requires a
806  // parallel communication.
807  elems_of_dimension.resize(4, false); // will use 1-based
808 
809  // Fills in the:
810  // global_elem_blk_ids[] and
811  // global_elem_blk_cnts[] arrays.
812  nemhelper->get_eb_info_global();
813 
814  // // Fills in the vectors
815  // // elem_mapi[num_internal_elems]
816  // // elem_mapb[num_border_elems ]
817  // // These tell which of the (locally-numbered) elements are internal and which are border elements.
818  // // In our test example these arrays are sorted (but non-contiguous), which makes it possible to
819  // // binary search for each element ID... however I don't think we need to distinguish between the
820  // // two types, since either can have nodes the boundary!
821  // nemhelper->get_elem_map();
822 
823  // Fills in the vectors of vectors:
824  // elem_cmap_elem_ids[][]
825  // elem_cmap_side_ids[][]
826  // elem_cmap_proc_ids[][]
827  // These arrays are of size num_elem_cmaps * elem_cmap_elem_cnts[i], i = 0..num_elem_cmaps
828  nemhelper->get_elem_cmap();
829 
830  // Get information about the element blocks:
831  // (read in the array nemhelper->block_ids[])
832  nemhelper->read_block_info();
833 
834  // Reads the nemhelper->elem_num_map array.
835  // elem_num_map[i] is the exodus element number for local element
836  // number i, which makes elem_num_map[i]-1 the libMesh element
837  // number.
838  nemhelper->read_elem_num_map();
839 
840  std::size_t local_elem_num = 0;
841 
842  // Read in the element connectivity for each block by
843  // looping over all the blocks.
844  for (unsigned int i=0; i<to_uint(nemhelper->num_elem_blk); i++)
845  {
846  // Read the information for block i: For nemhelper->block_ids[i], reads
847  // elem_type
848  // num_elem_this_blk
849  // num_nodes_per_elem
850  // num_attr
851  // connect <-- the nodal connectivity array for each element in the block.
852  nemhelper->read_elem_in_block(i);
853 
854  // Note that with parallel files it is possible we have no elements in
855  // this block!
856  if (!nemhelper->num_elem_this_blk) continue;
857 
858  // Set subdomain ID based on the block ID.
859  subdomain_id_type subdomain_id =
860  cast_int<subdomain_id_type>(nemhelper->block_ids[i]);
861 
862  // Create a type string (this uses the null-terminated string ctor).
863  const std::string type_str ( nemhelper->elem_type.data() );
864 
865  // Set any relevant node/edge maps for this element
866  const auto & conv = nemhelper->get_conversion(type_str);
867 
868  if (_verbose)
869  libMesh::out << "Reading a block of " << type_str << " elements." << std::endl;
870 
871  // Loop over all the elements in this block
872  for (unsigned int j=0; j<to_uint(nemhelper->num_elem_this_blk); j++)
873  {
874  auto uelem = Elem::build (conv.libmesh_elem_type());
875 
876  // Assign subdomain and processor ID to the newly-created Elem.
877  // Assigning the processor ID beforehand ensures that the Elem is
878  // not added as an "unpartitioned" element. Note that the element
879  // numbering in Exodus is also 1-based.
880  uelem->subdomain_id() = subdomain_id;
881  uelem->processor_id() = this->processor_id();
882  uelem->set_id() = nemhelper->elem_num_map[local_elem_num++]-1;
883 
884  // Handle unique_id numbering, just in case we're using a
885  // ReplicatedMesh that doesn't know how to handle it in
886  // parallel.
887 #ifdef LIBMESH_ENABLE_UNIQUE_ID
888  uelem->set_unique_id(uelem->id());
889 #endif
890 
891  // Mark that we have seen an element of the current element's
892  // dimension.
893  elems_of_dimension[uelem->dim()] = true;
894 
895  // Add the created Elem to the Mesh, catch the Elem
896  // pointer that the Mesh throws back.
897  Elem * elem = mesh.add_elem(std::move(uelem));
898 
899  // We are expecting the element "thrown back" by libmesh to have the ID we specified for it...
900  // Check to see that really is the case. Note that local_elem_num was post-incremented, so
901  // subtract 1 when performing the check.
902  libmesh_assert_equal_to(elem->id(),
903  cast_int<dof_id_type>(nemhelper->elem_num_map[local_elem_num-1]-1));
904 
905  // Set all the nodes for this element
906  if (_verbose)
907  libMesh::out << "[" << this->processor_id() << "] "
908  << "Setting nodes for Elem " << elem->id() << std::endl;
909 
910  for (unsigned int k=0; k<to_uint(nemhelper->num_nodes_per_elem); k++)
911  {
912  const unsigned int
913  gi = (j*nemhelper->num_nodes_per_elem + // index into connectivity array
914  conv.get_node_map(k)),
915  local_node_idx = nemhelper->connect[gi]-1, // local node index
916  global_node_idx = nemhelper->node_num_map[local_node_idx]-1; // new global node index
917 
918  // Set node number
919  elem->set_node(k, mesh.node_ptr(global_node_idx));
920  }
921  } // for (unsigned int j=0; j<nemhelper->num_elem_this_blk; j++)
922  } // end for (unsigned int i=0; i<nemhelper->num_elem_blk; i++)
923 
924  for (const auto & [id, name] : nemhelper->id_to_block_names)
925  if (name != "")
926  mesh.subdomain_name(id) = name;
927 
928  if (_verbose)
929  {
930  // Print local elems_of_dimension information
931  for (auto i : IntRange<std::size_t>(1, elems_of_dimension.size()))
932  libMesh::out << "[" << this->processor_id() << "] "
933  << "elems_of_dimension[" << i << "]=" << elems_of_dimension[i] << std::endl;
934  }
935 
936  // Get the max dimension seen on the current processor
937  unsigned char max_dim_seen = 0;
938  for (auto i : IntRange<std::size_t>(1, elems_of_dimension.size()))
939  if (elems_of_dimension[i])
940  max_dim_seen = static_cast<unsigned char>(i);
941 
942  // Do a global max to determine the max dimension seen by all processors.
943  // It should match -- I don't think we even support calculations on meshes
944  // with elements of different dimension...
945  this->comm().max(max_dim_seen);
946 
947  if (_verbose)
948  {
949  // Print the max element dimension from all processors
950  libMesh::out << "[" << this->processor_id() << "] "
951  << "max_dim_seen=" << +max_dim_seen << std::endl;
952  }
953 
954  // Set the mesh dimension to the largest encountered for an element
955  mesh.set_mesh_dimension(max_dim_seen);
956 
957 #if LIBMESH_DIM < 3
958  libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
959  "Cannot open dimension "
960  << mesh.mesh_dimension()
961  << " mesh file when configured without "
962  << mesh.mesh_dimension()
963  << "D support." );
964 #endif
965 
966 
967  // Global sideset information, they are distributed as well, not sure if they will require communication...?
968  nemhelper->get_ss_param_global();
969 
970  if (_verbose)
971  {
972  libMesh::out << "[" << this->processor_id() << "] "
973  << "Read global sideset parameter information." << std::endl;
974 
975  // These global values should be the same on all processors...
976  libMesh::out << "[" << this->processor_id() << "] "
977  << "Number of global sideset IDs: " << nemhelper->global_sideset_ids.size() << std::endl;
978  }
979 
980  // Read *local* sideset info the same way it is done in
981  // exodusII_io_helper. May be called any time after
982  // nemhelper->read_and_store_header_info(); This sets num_side_sets and resizes
983  // elem_list, side_list, and id_list to num_elem_all_sidesets. Note
984  // that there appears to be the same number of sidesets in each file
985  // but they all have different numbers of entries (some are empty).
986  // Note that the sum of "nemhelper->num_elem_all_sidesets" over all
987  // processors should equal the sum of the entries in the "num_global_side_counts" array
988  // filled up by nemhelper->get_ss_param_global()
989  nemhelper->read_sideset_info();
990 
991  if (_verbose)
992  {
993  libMesh::out << "[" << this->processor_id() << "] "
994  << "nemhelper->num_side_sets = " << nemhelper->num_side_sets << std::endl;
995 
996  libMesh::out << "[" << this->processor_id() << "] "
997  << "nemhelper->num_elem_all_sidesets = " << nemhelper->num_elem_all_sidesets << std::endl;
998 
999  if (nemhelper->num_side_sets > 0)
1000  {
1001  libMesh::out << "Sideset names are: ";
1002  for (const auto & [id, name] : nemhelper->id_to_ss_names)
1003  libMesh::out << "(" << id << "," << name << ") ";
1004  libMesh::out << std::endl;
1005  }
1006  }
1007 
1008 #ifdef DEBUG
1009  {
1010  // In DEBUG mode, check that the global number of sidesets reported
1011  // in each nemesis file matches the sum of all local sideset counts
1012  // from each processor. This requires a small communication, so only
1013  // do it in DEBUG mode.
1014  int sum_num_global_side_counts = std::accumulate(nemhelper->num_global_side_counts.begin(),
1015  nemhelper->num_global_side_counts.end(),
1016  0);
1017 
1018  // MPI sum up the local files contributions
1019  int sum_num_elem_all_sidesets = nemhelper->num_elem_all_sidesets;
1020  this->comm().sum(sum_num_elem_all_sidesets);
1021 
1022  libmesh_error_msg_if(sum_num_global_side_counts != sum_num_elem_all_sidesets,
1023  "Error! global side count reported by Nemesis does not "
1024  "match the side count reported by the individual files!");
1025  }
1026 #endif
1027 
1028  // Note that exodus stores sidesets in separate vectors but we want to pack
1029  // them all into a single vector. So when we call read_sideset(), we pass an offset
1030  // into the single vector of all IDs
1031  for (int offset=0, i=0; i<nemhelper->num_side_sets; i++)
1032  {
1033  offset += (i > 0 ? nemhelper->num_sides_per_set[i-1] : 0); // Compute new offset
1034  nemhelper->read_sideset (i, offset);
1035  }
1036 
1037  // Now that we have the lists of elements, sides, and IDs, we are ready to set them
1038  // in the BoundaryInfo object of our Mesh object. This is slightly different in parallel...
1039  // For example, I think the IDs in each of the split Exodus files are numbered locally,
1040  // and we have to know the appropriate ID for this processor to be able to set the
1041  // entry in BoundaryInfo. This id should be given by
1042  // elem_num_map[i]-1 for the local index i
1043 
1044  // Debugging:
1045  // Print entries of elem_list
1046  // libMesh::out << "[" << this->processor_id() << "] "
1047  // << "elem_list = ";
1048  // for (const auto & id : nemhelper->elem_list)
1049  // libMesh::out << id << ", ";
1050  // libMesh::out << std::endl;
1051 
1052  // Print entries of side_list
1053  // libMesh::out << "[" << this->processor_id() << "] "
1054  // << "side_list = ";
1055  // for (const auto & id : nemhelper->side_list)
1056  // libMesh::out << id << ", ";
1057  // libMesh::out << std::endl;
1058 
1059 
1060  // Loop over the entries of the elem_list, get their pointers from the
1061  // Mesh data structure, and assign the appropriate side to the BoundaryInfo object.
1062  for (auto e : index_range(nemhelper->elem_list))
1063  {
1064  // Exodus numbering is 1-based
1065  const std::size_t local_id = nemhelper->elem_list[e]-1;
1066  const dof_id_type elem_id = nemhelper->elem_num_map[local_id]-1;
1067 
1068  Elem * elem = mesh.elem_ptr(elem_id);
1069 
1070  // The side numberings in libmesh and exodus are not 1:1, so we need to map
1071  // whatever side number is stored in Exodus into a libmesh side number using
1072  // a conv object...
1073  const auto & conv = nemhelper->get_conversion(elem->type());
1074 
1075  // Finally, we are ready to add the element and its side to the BoundaryInfo object.
1076  // Call the version of add_side which takes a pointer, since we have already gone to
1077  // the trouble of getting said pointer...
1079  cast_int<unsigned short>(conv.get_side_map(nemhelper->side_list[e]-1)), // Exodus numbering is 1-based
1080  cast_int<boundary_id_type>(nemhelper->id_list[e]));
1081  }
1082 
1083  for (const auto & [id, name] : nemhelper->id_to_ss_names)
1084  if (name != "")
1086 
1087  // Debugging: make sure there are as many boundary conditions in the
1088  // boundary ID object as expected. Note that, at this point, the
1089  // mesh still thinks it's serial, so n_boundary_conds() returns the
1090  // local number of boundary conditions (and is therefore cheap)
1091  // which should match nemhelper->elem_list.size().
1092  {
1093  std::size_t nbcs = mesh.get_boundary_info().n_boundary_conds();
1094  libmesh_error_msg_if(nbcs != nemhelper->elem_list.size(),
1095  "[" << this->processor_id() << "] "
1096  << "BoundaryInfo contains "
1097  << nbcs
1098  << " boundary conditions, while the Exodus file had "
1099  << nemhelper->elem_list.size());
1100  }
1101 
1102  // Read global nodeset parameters? We might be able to use this to verify
1103  // something about the local files, but I haven't figured out what yet...
1104  nemhelper->get_ns_param_global();
1105 
1106  // Read local nodeset info
1107  nemhelper->read_nodeset_info();
1108 
1109  if (_verbose)
1110  {
1111  libMesh::out << "[" << this->processor_id() << "] ";
1112  libMesh::out << "nemhelper->num_node_sets=" << nemhelper->num_node_sets << std::endl;
1113  if (nemhelper->num_node_sets > 0)
1114  {
1115  libMesh::out << "Nodeset names are: ";
1116  for (const auto & [id, name] : nemhelper->id_to_ns_names)
1117  libMesh::out << "(" << id << "," << name << ") ";
1118  libMesh::out << std::endl;
1119  }
1120  }
1121 
1122  // // Debugging, what is currently in nemhelper->node_num_map anyway?
1123  // libMesh::out << "[" << this->processor_id() << "] "
1124  // << "nemhelper->node_num_map = ";
1125  //
1126  // for (const auto & id : nemhelper->node_num_map)
1127  // libMesh::out << id << ", ";
1128  // libMesh::out << std::endl;
1129 
1130  // For each nodeset,
1131  for (int nodeset=0; nodeset<nemhelper->num_node_sets; nodeset++)
1132  {
1133  // Get the user-defined ID associated with the nodeset
1134  int nodeset_id = nemhelper->nodeset_ids[nodeset];
1135 
1136  if (_verbose)
1137  {
1138  libMesh::out << "[" << this->processor_id() << "] ";
1139  libMesh::out << "nemhelper->nodeset_ids[" << nodeset << "]=" << nodeset_id << std::endl;
1140  }
1141 
1142  // Read the nodeset from file, store them in a vector
1143  nemhelper->read_nodeset(nodeset);
1144 
1145  // Add nodes from the node_list to the BoundaryInfo object
1146  for (auto node : index_range(nemhelper->node_list))
1147  {
1148  // Don't run past the end of our node map!
1149  libmesh_error_msg_if(to_uint(nemhelper->node_list[node]-1) >= nemhelper->node_num_map.size(),
1150  "Error, index is past the end of node_num_map array!");
1151 
1152  // We should be able to use the node_num_map data structure set up previously to determine
1153  // the proper global node index.
1154  unsigned global_node_id = nemhelper->node_num_map[ nemhelper->node_list[node]-1 /*Exodus is 1-based!*/ ]-1;
1155 
1156  if (_verbose)
1157  {
1158  libMesh::out << "[" << this->processor_id() << "] "
1159  << "nodeset " << nodeset
1160  << ", local node number: " << nemhelper->node_list[node]-1
1161  << ", global node id: " << global_node_id
1162  << std::endl;
1163  }
1164 
1165  // Add the node to the BoundaryInfo object with the proper nodeset_id
1167  (cast_int<dof_id_type>(global_node_id),
1168  cast_int<boundary_id_type>(nodeset_id));
1169  }
1170  }
1171 
1172  for (const auto & [id, name] : nemhelper->id_to_ns_names)
1173  if (name != "")
1175 
1176  // See what the elem count is up to now.
1177  if (_verbose)
1178  {
1179  // Report the number of elements which have been added locally
1180  libMesh::out << "[" << this->processor_id() << "] ";
1181  libMesh::out << "mesh.n_elem()=" << mesh.n_elem() << std::endl;
1182 
1183  // Reports the number of elements that have been added in total.
1184  libMesh::out << "[" << this->processor_id() << "] ";
1185  libMesh::out << "mesh.parallel_n_elem()=" << mesh.parallel_n_elem() << std::endl;
1186  }
1187 
1188  // For DistributedMesh, it seems that _is_serial is true by default. A hack to
1189  // make the Mesh think it's parallel might be to call:
1193 
1194  // If that didn't work, then we're actually reading into a
1195  // ReplicatedMesh, so we want to gather *all* elements
1196  if (mesh.is_serial())
1197  // Don't just use mesh.allgather(); that's a no-op, since
1198  // ReplicatedMesh didn't expect to be distributed in the first
1199  // place!
1201  else
1202  // Gather neighboring elements so that a distributed mesh has the
1203  // proper "ghost" neighbor information.
1204  MeshCommunication().gather_neighboring_elements(cast_ref<DistributedMesh &>(mesh));
1205 
1206 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1207  // We've been setting unique_ids by hand; let's make sure that later
1208  // ones are consistent with them.
1210 #endif
1211 }
1212 
1213 #else
1214 
1215 void Nemesis_IO::read (const std::string &)
1216 {
1217  libmesh_error_msg("ERROR, Nemesis API is not defined!");
1218 }
1219 
1220 #endif // #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1221 
1222 
1223 
1224 
1225 
1226 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1227 
1228 void Nemesis_IO::write (const std::string & base_filename)
1229 {
1230  // Get a constant reference to the mesh for writing
1232 
1233  // Create the filename for this processor given the base_filename passed in.
1234  std::string nemesis_filename = nemhelper->construct_nemesis_filename(base_filename);
1235 
1236  // If the user has set the append flag here, it doesn't really make
1237  // sense: the intent of this function is to write a Mesh with no
1238  // data, while "appending" is really intended to add data to an
1239  // existing file. If we're verbose, print a message to this effect.
1240  if (_append && _verbose)
1241  libmesh_warning("Warning: Appending in Nemesis_IO::write() does not make sense.\n"
1242  "Creating a new file instead!");
1243 
1244  nemhelper->create(nemesis_filename);
1245 
1246  // Initialize data structures and write some global Nemesis-specific data, such as
1247  // communication maps, to file.
1248  nemhelper->initialize(nemesis_filename,mesh);
1249 
1250  // Make sure we're writing communication maps we can reuse as
1251  // expected when reading
1252  this->assert_symmetric_cmaps();
1253 
1254  // Call the Nemesis-specialized version of write_nodal_coordinates() to write
1255  // the nodal coordinates.
1256  nemhelper->write_nodal_coordinates(mesh);
1257 
1258  // Call the Nemesis-specialized version of write_elements() to write
1259  // the elements. Note: Must write a zero if a given global block ID has no
1260  // elements...
1261  nemhelper->write_elements(mesh);
1262 
1263  // Call our specialized function to write the nodesets
1264  nemhelper->write_nodesets(mesh);
1265 
1266  // Call our specialized write_sidesets() function to write the sidesets to file
1267  nemhelper->write_sidesets(mesh);
1268 
1269  // Not sure if this is really necessary, but go ahead and flush the file
1270  // once we have written all this stuff.
1271  nemhelper->update();
1272 
1273  if ((mesh.get_boundary_info().n_edge_conds() > 0) && _verbose)
1274  libmesh_warning("Warning: Mesh contains edge boundary IDs, but these "
1275  "are not supported by the Nemesis format.");
1276 }
1277 
1278 #else
1279 
1280 void Nemesis_IO::write (const std::string & )
1281 {
1282  libmesh_error_msg("ERROR, Nemesis API is not defined!");
1283 }
1284 
1285 #endif // #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1286 
1287 
1288 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1289 
1290 void Nemesis_IO::write_timestep (const std::string & fname,
1291  const EquationSystems & es,
1292  const int timestep,
1293  const Real time)
1294 {
1295  _timestep=timestep;
1296  write_equation_systems(fname,es);
1297 
1298  nemhelper->write_timestep(timestep, time);
1299 }
1300 
1301 #else
1302 
1303 void Nemesis_IO::write_timestep (const std::string &,
1304  const EquationSystems &,
1305  const int,
1306  const Real)
1307 {
1308  libmesh_error_msg("ERROR, Nemesis API is not defined!");
1309 }
1310 
1311 #endif // #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1312 
1313 
1314 
1315 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1316 
1317 void Nemesis_IO::prepare_to_write_nodal_data (const std::string & fname,
1318  const std::vector<std::string> & names)
1319 {
1321 
1322  std::string nemesis_filename = nemhelper->construct_nemesis_filename(fname);
1323 
1324  if (!nemhelper->opened_for_writing)
1325  {
1326  // If we're appending, open() the file with read_only=false,
1327  // otherwise create() it and write the contents of the mesh to
1328  // it.
1329  if (_append)
1330  {
1331  nemhelper->open(nemesis_filename.c_str(), /*read_only=*/false);
1332  // After opening the file, read the header so that certain
1333  // fields, such as the number of nodes and the number of
1334  // elements, are correctly initialized for the subsequent
1335  // call to write the nodal solution.
1336  nemhelper->read_and_store_header_info();
1337 
1338  // ...and reading the block info
1339  nemhelper->read_block_info();
1340 
1341  // ...and rebuild the "exodus_node_num_to_libmesh" map
1342  nemhelper->compute_num_global_elem_blocks(mesh);
1343  nemhelper->build_element_and_node_maps(mesh);
1344  }
1345  else
1346  {
1347  nemhelper->create(nemesis_filename);
1348  nemhelper->initialize(nemesis_filename,mesh);
1349 
1350  // Make sure we're writing communication maps we can reuse
1351  // as expected when reading
1352  this->assert_symmetric_cmaps();
1353 
1354  nemhelper->write_nodal_coordinates(mesh);
1355  nemhelper->write_elements(mesh);
1356  nemhelper->write_nodesets(mesh);
1357  nemhelper->write_sidesets(mesh);
1358 
1359  if ((mesh.get_boundary_info().n_edge_conds() > 0) && _verbose)
1360  libmesh_warning("Warning: Mesh contains edge boundary IDs, but these "
1361  "are not supported by the ExodusII format.");
1362  }
1363  }
1364 
1365  // Even if we were already open for writing, we might not have
1366  // initialized the nodal variable names yet. Even if we did, it
1367  // should not hurt to call this twice because the routine sets a
1368  // flag the first time it is called.
1369 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1370  std::vector<std::string> complex_names =
1371  nemhelper->get_complex_names(names, nemhelper->write_complex_abs);
1372  nemhelper->initialize_nodal_variables(complex_names);
1373 #else
1374  nemhelper->initialize_nodal_variables(names);
1375 #endif
1376 }
1377 
1378 #else
1379 
1380 void Nemesis_IO::prepare_to_write_nodal_data (const std::string &,
1381  const std::vector<std::string> &)
1382 {
1383  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1384 }
1385 
1386 #endif
1387 
1388 
1389 
1390 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1391 
1392 void Nemesis_IO::write_nodal_data (const std::string & base_filename,
1393  const NumericVector<Number> & parallel_soln,
1394  const std::vector<std::string> & names)
1395 {
1396  LOG_SCOPE("write_nodal_data(parallel)", "Nemesis_IO");
1397 
1398  // Only prepare and write nodal variables that are also in
1399  // _output_variables, unless _output_variables is empty. This is the
1400  // same logic that is in ExodusII_IO::write_nodal_data().
1401  std::vector<std::string> output_names;
1402 
1403  if (_allow_empty_variables || !_output_variables.empty())
1404  output_names = _output_variables;
1405  else
1406  output_names = names;
1407 
1408  this->prepare_to_write_nodal_data(base_filename, output_names);
1409 
1410  // Call the new version of write_nodal_solution() that takes a
1411  // NumericVector directly without localizing.
1412  nemhelper->write_nodal_solution(parallel_soln, names, _timestep, output_names);
1413 }
1414 
1415 
1416 
1417 void Nemesis_IO::write_nodal_data (const std::string & base_filename,
1418  const EquationSystems & es,
1419  const std::set<std::string> * system_names)
1420 {
1421  LOG_SCOPE("write_nodal_data(parallel)", "Nemesis_IO");
1422 
1423  // Only prepare and write nodal variables that are also in
1424  // _output_variables, unless _output_variables is empty. This is the
1425  // same logic that is in ExodusII_IO::write_nodal_data().
1426  std::vector<std::string> output_names;
1427 
1428  if (_allow_empty_variables || !_output_variables.empty())
1429  output_names = _output_variables;
1430  else
1431  es.build_variable_names (output_names, nullptr, system_names);
1432 
1433  this->prepare_to_write_nodal_data(base_filename, output_names);
1434 
1435  std::vector<std::pair<unsigned int, unsigned int>> var_nums;
1436  // If we pass in an empty vector below, it will return all of the
1437  // var nums in es, which we don't want.
1438  if (!output_names.empty())
1439  var_nums = es.find_variable_numbers(output_names);
1440 
1441  nemhelper->write_nodal_solution(es, var_nums, _timestep, output_names);
1442 }
1443 
1444 #else
1445 
1446 void Nemesis_IO::write_nodal_data (const std::string &,
1447  const NumericVector<Number> &,
1448  const std::vector<std::string> &)
1449 {
1450  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1451 }
1452 
1453 void Nemesis_IO::write_nodal_data (const std::string &,
1454  const EquationSystems &,
1455  const std::set<std::string> *)
1456 {
1457  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1458 }
1459 
1460 #endif
1461 
1462 
1463 
1464 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1465 
1467 {
1468  libmesh_error_msg_if(!nemhelper->opened_for_writing,
1469  "ERROR, Nemesis file must be initialized before outputting elemental variables.");
1470 
1471  // To be (possibly) filled with a filtered list of variable names to output.
1472  std::vector<std::string> names;
1473 
1474  // All of which should be low order monomials for now
1475  const std::vector<FEType> type = {FEType(CONSTANT, MONOMIAL), FEType(CONSTANT, MONOMIAL_VEC)};
1476 
1477  // If _output_variables is populated, only output the monomials which are
1478  // also in the _output_variables vector.
1479  if (_output_variables.size())
1480  {
1481  std::vector<std::string> monomials;
1482 
1483  // Create a list of monomial variable names
1484  es.build_variable_names(monomials, &type[0]); /*scalars*/
1485  es.build_variable_names(monomials, &type[1]); /*vectors*/
1486 
1487  // Filter that list against the _output_variables list. Note: if names is still empty after
1488  // all this filtering, all the monomial variables will be gathered
1489  for (const auto & var : monomials)
1490  if (std::find(_output_variables.begin(), _output_variables.end(), var) != _output_variables.end())
1491  names.push_back(var);
1492  }
1493 
1494  // The 'names' vector will here be updated with the variable's names
1495  // that are actually eligible to write
1496  std::vector<std::pair<unsigned int, unsigned int>> var_nums =
1497  es.find_variable_numbers (names, /*type=*/nullptr, &type);
1498 
1499  // find_variable_numbers() can return a nullptr, in which case there are no constant monomial
1500  // variables to write, and we can just return.
1501  if (var_nums.empty())
1502  {
1503  if (_verbose)
1504  libMesh::out << "No CONSTANT, MONOMIAL or CONSTANT, MONOMIAL_VEC data to be written." << std::endl;
1505  return;
1506  }
1507 
1508  // Store the list of subdomains on which each variable *that we are
1509  // going to plot* is active. Note: if any of these sets is _empty_,
1510  // the variable in question is active on _all_ subdomains.
1511  std::vector<std::set<subdomain_id_type>> vars_active_subdomains;
1512  es.get_vars_active_subdomains(names, vars_active_subdomains);
1513 
1515 
1516 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1517  std::vector<std::string> complex_names =
1518  nemhelper->get_complex_names(names, nemhelper->write_complex_abs);
1519 
1520  std::vector<std::set<subdomain_id_type>>
1521  complex_vars_active_subdomains =
1522  nemhelper->get_complex_vars_active_subdomains(vars_active_subdomains,
1523  nemhelper->write_complex_abs);
1524  nemhelper->initialize_element_variables(complex_names, complex_vars_active_subdomains);
1525 
1526  // Call (non-virtual) function to write the elemental data in
1527  // parallel. This function is named similarly to the corresponding
1528  // function in the Exodus helper, but it has a different calling
1529  // sequence and is not virtual or an override.
1530  nemhelper->write_element_values(mesh,
1531  es,
1532  var_nums,
1533  _timestep,
1534  complex_vars_active_subdomains);
1535 
1536 #else
1537  // Call the Nemesis version of initialize_element_variables().
1538  //
1539  // The Exodus helper version of this function writes an incorrect
1540  // truth table in parallel that somehow does not account for the
1541  // case where a subdomain does not appear on one or more of the
1542  // processors. So, we override that function's behavior in the
1543  // Nemesis helper.
1544  nemhelper->initialize_element_variables(names, vars_active_subdomains);
1545 
1546  // Call (non-virtual) function to write the elemental data in
1547  // parallel. This function is named similarly to the corresponding
1548  // function in the Exodus helper, but it has a different calling
1549  // sequence and is not virtual or an override.
1550  nemhelper->write_element_values(mesh,
1551  es,
1552  var_nums,
1553  _timestep,
1554  vars_active_subdomains);
1555 #endif
1556 }
1557 
1558 #else
1559 
1561 {
1562  libmesh_not_implemented();
1563 }
1564 
1565 #endif
1566 
1567 
1568 
1569 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1570 
1571 void Nemesis_IO::write_nodal_data (const std::string & base_filename,
1572  const std::vector<Number> & soln,
1573  const std::vector<std::string> & names)
1574 {
1575  LOG_SCOPE("write_nodal_data(serialized)", "Nemesis_IO");
1576 
1577  this->prepare_to_write_nodal_data(base_filename, names);
1578 
1579  nemhelper->write_nodal_solution(soln, names, _timestep);
1580 }
1581 
1582 #else
1583 
1584 void Nemesis_IO::write_nodal_data (const std::string &,
1585  const std::vector<Number> &,
1586  const std::vector<std::string> &)
1587 {
1588  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1589 }
1590 
1591 #endif
1592 
1593 
1594 
1595 
1596 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1597 
1598 void Nemesis_IO::write_global_data (const std::vector<Number> & soln,
1599  const std::vector<std::string> & names)
1600 {
1601  libmesh_error_msg_if(!nemhelper->opened_for_writing,
1602  "ERROR, Nemesis file must be initialized before outputting global variables.");
1603 
1604 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1605 
1606  std::vector<std::string> complex_names =
1607  nemhelper->get_complex_names(names, nemhelper->write_complex_abs);
1608 
1609  nemhelper->initialize_global_variables(complex_names);
1610 
1611  unsigned int num_values = soln.size();
1612  unsigned int num_vars = names.size();
1613  unsigned int num_elems = num_values / num_vars;
1614 
1615  // This will contain the real and imaginary parts and the magnitude
1616  // of the values in soln
1617  int nco = nemhelper->write_complex_abs ? 3 : 2;
1618  std::vector<Real> complex_soln(nco * num_values);
1619 
1620  for (unsigned i=0; i<num_vars; ++i)
1621  {
1622  for (unsigned int j=0; j<num_elems; ++j)
1623  {
1624  Number value = soln[i*num_vars + j];
1625  complex_soln[nco*i*num_elems + j] = value.real();
1626  }
1627  for (unsigned int j=0; j<num_elems; ++j)
1628  {
1629  Number value = soln[i*num_vars + j];
1630  complex_soln[nco*i*num_elems + num_elems + j] = value.imag();
1631  }
1632  if (nemhelper->write_complex_abs)
1633  {
1634  for (unsigned int j=0; j<num_elems; ++j)
1635  {
1636  Number value = soln[i*num_vars + j];
1637  complex_soln[3*i*num_elems + 2*num_elems + j] = std::abs(value);
1638  }
1639  }
1640  }
1641 
1642  nemhelper->write_global_values(complex_soln, _timestep);
1643 
1644 #else
1645 
1646  // Call the Exodus writer implementation
1647  nemhelper->initialize_global_variables( names );
1648  nemhelper->write_global_values( soln, _timestep);
1649 
1650 #endif
1651 
1652 }
1653 
1654 #else
1655 
1656 void Nemesis_IO::write_global_data (const std::vector<Number> &,
1657  const std::vector<std::string> &)
1658 {
1659  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1660 }
1661 
1662 #endif // #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1663 
1664 
1665 
1666 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1667 
1668 void Nemesis_IO::write_information_records (const std::vector<std::string> & records)
1669 {
1670  libmesh_error_msg_if(!nemhelper->opened_for_writing,
1671  "ERROR, Nemesis file must be initialized before outputting information records.");
1672 
1673  // Call the Exodus writer implementation
1674  nemhelper->write_information_records( records );
1675 }
1676 
1677 
1678 const std::vector<std::string> & Nemesis_IO::get_nodal_var_names()
1679 {
1680  nemhelper->read_var_names(ExodusII_IO_Helper::NODAL);
1681  return nemhelper->nodal_var_names;
1682 }
1683 
1684 
1685 
1687  std::string system_var_name,
1688  std::string exodus_var_name,
1689  unsigned int timestep)
1690 {
1691  libmesh_error_msg_if(!nemhelper->opened_for_reading,
1692  "ERROR, Nemesis file must be opened for reading before copying a nodal solution!");
1693 
1694  nemhelper->read_nodal_var_values(exodus_var_name, timestep);
1695 
1696  const unsigned int var_num = system.variable_number(system_var_name);
1697 
1698  for (auto p : nemhelper->nodal_var_values)
1699  {
1700  dof_id_type i = p.first;
1701  const Node * node = MeshInput<MeshBase>::mesh().node_ptr(i);
1702 
1703  if (node && node->n_comp(system.number(), var_num) > 0)
1704  {
1705  dof_id_type dof_index = node->dof_number(system.number(), var_num, 0);
1706 
1707  // If the dof_index is local to this processor, set the value
1708  if (system.get_dof_map().local_index(dof_index))
1709  system.solution->set (dof_index, p.second);
1710  }
1711  }
1712 
1713  system.solution->close();
1714  system.update();
1715 }
1716 
1717 
1718 
1720  std::string system_var_name,
1721  std::string exodus_var_name,
1722  unsigned int timestep)
1723 {
1724  parallel_object_only();
1725 
1726  const unsigned int var_num = system.variable_number(system_var_name);
1727  libmesh_error_msg_if(system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL),
1728  "Error! Trying to copy elemental solution into a variable that is not of CONSTANT MONOMIAL type.");
1729 
1731 
1732  // Map from element ID to elemental variable value. We need to use
1733  // a map here rather than a vector (e.g. elem_var_values) since the
1734  // libmesh element numbering can contain "holes". This is the case
1735  // if we are reading elemental var values from an adaptively refined
1736  // mesh that has not been sequentially renumbered.
1737  std::map<dof_id_type, Real> elem_var_value_map;
1738 
1739  libmesh_error_msg_if(!nemhelper->opened_for_reading,
1740  "ERROR, Nemesis file must be opened for reading before copying an elemental solution!");
1741 
1742  nemhelper->read_elemental_var_values(exodus_var_name, timestep, elem_var_value_map);
1743 
1744  std::map<dof_id_type, Real>::iterator
1745  it = elem_var_value_map.begin(),
1746  end = elem_var_value_map.end();
1747 
1748  for (; it!=end; ++it)
1749  {
1750  const Elem * elem = mesh.query_elem_ptr(it->first);
1751 
1752  if (elem && elem->n_comp(system.number(), var_num) > 0)
1753  {
1754  dof_id_type dof_index = elem->dof_number(system.number(), var_num, 0);
1755  libmesh_assert(system.get_dof_map().local_index(dof_index));
1756  system.solution->set (dof_index, it->second);
1757  }
1758  }
1759 
1760  system.solution->close();
1761  system.update();
1762 
1763  parallel_object_only();
1764 }
1765 
1766 
1767 
1769  std::vector<std::string> system_var_names,
1770  std::vector<std::string> exodus_var_names,
1771  unsigned int timestep)
1772 {
1773  libmesh_error_msg_if(!nemhelper->opened_for_reading,
1774  "ERROR, Nemesis file must be opened for reading before copying a scalar solution!");
1775 
1776  libmesh_error_msg_if(system_var_names.size() != exodus_var_names.size(),
1777  "ERROR, the number of system_var_names must match exodus_var_names.");
1778 
1779  std::vector<Real> values_from_exodus;
1780  read_global_variable(exodus_var_names, timestep, values_from_exodus);
1781 
1782  if (system.processor_id() == (system.n_processors()-1))
1783  {
1784  const DofMap & dof_map = system.get_dof_map();
1785 
1786  for (auto i : index_range(system_var_names))
1787  {
1788  const unsigned int var_num = system.variable_scalar_number(system_var_names[i], 0);
1789 
1790  std::vector<dof_id_type> SCALAR_dofs;
1791  dof_map.SCALAR_dof_indices(SCALAR_dofs, var_num);
1792 
1793  system.solution->set (SCALAR_dofs[0], values_from_exodus[i]);
1794  }
1795  }
1796 
1797  system.solution->close();
1798  system.update();
1799 }
1800 
1801 
1802 void Nemesis_IO::read_global_variable(std::vector<std::string> global_var_names,
1803  unsigned int timestep,
1804  std::vector<Real> & global_values)
1805 {
1806  std::size_t size = global_var_names.size();
1807  libmesh_error_msg_if(size == 0, "ERROR, empty list of global variables to read from the Nemesis file.");
1808 
1809  // read the values for all global variables
1810  std::vector<Real> values_from_exodus;
1811  nemhelper->read_var_names(ExodusII_IO_Helper::GLOBAL);
1812  nemhelper->read_global_values(values_from_exodus, timestep);
1813  std::vector<std::string> global_var_names_exodus = nemhelper->global_var_names;
1814 
1815  if (values_from_exodus.size() == 0)
1816  return; // This will happen in parallel on procs that are not 0
1817 
1818  global_values.clear();
1819  for (std::size_t i = 0; i != size; ++i)
1820  {
1821  // for each global variable in global_var_names, look the corresponding one in global_var_names_from_exodus
1822  // and fill global_values accordingly
1823  auto it = find(global_var_names_exodus.begin(), global_var_names_exodus.end(), global_var_names[i]);
1824  if (it != global_var_names_exodus.end())
1825  global_values.push_back(values_from_exodus[it - global_var_names_exodus.begin()]);
1826  else
1827  libmesh_error_msg("ERROR, Global variable " << global_var_names[i] << \
1828  " not found in Nemesis file.");
1829  }
1830 }
1831 
1832 void Nemesis_IO::set_hdf5_writing(bool write_hdf5)
1833 {
1834  nemhelper->set_hdf5_writing(write_hdf5);
1835 }
1836 
1837 #else
1838 
1839 void Nemesis_IO::write_information_records ( const std::vector<std::string> & )
1840 {
1841  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1842 }
1843 
1844 const std::vector<std::string> & Nemesis_IO::get_nodal_var_names()
1845 {
1846  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1847 
1848  // Prevent potential compiler warnings about missing return statement
1849  return _output_variables;
1850 }
1851 
1852 void Nemesis_IO::copy_nodal_solution(System &, std::string, std::string, unsigned int)
1853 {
1854  libmesh_error_msg("ERROR, Nemesis API is not defined.");
1855 }
1856 
1857 void Nemesis_IO::set_hdf5_writing(bool) {}
1858 
1859 #endif // #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
1860 
1861 
1862 
1863 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
OStreamProxy err
const MT & mesh() const
Definition: mesh_output.h:259
std::size_t n_boundary_conds() const
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
unique_id_type & set_unique_id()
Definition: dof_object.h:858
This is the EquationSystems class.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
unsigned int variable_scalar_number(std::string_view var, unsigned int component) const
Definition: system.h:2489
void make_node_unique_ids_parallel_consistent(MeshBase &)
Assuming all unique_ids on local nodes are globally unique, and assuming all processor ids are parall...
void build_variable_names(std::vector< std::string > &var_names, const FEType *type=nullptr, const std::set< std::string > *system_names=nullptr) const
Fill the input vector var_names with the names of the variables for each system.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
A Node is like a Point, but with more information.
Definition: node.h:52
virtual unique_id_type parallel_max_unique_id() const =0
std::string & nodeset_name(boundary_id_type id)
virtual void read(const std::string &base_filename) override
Implements reading the mesh from several different files.
Definition: nemesis_io.C:214
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
std::size_t n_edge_conds() const
void copy_nodal_solution(System &system, std::string system_var_name, std::string exodus_var_name, unsigned int timestep=1)
If we read in a nodal solution while reading in a mesh, we can attempt to copy that nodal solution in...
Definition: nemesis_io.C:1686
std::vector< bool > elems_of_dimension
A vector of bools describing what dimension elements have been encountered when reading a mesh...
Definition: mesh_input.h:107
MessageTag get_unique_tag(int tagvalue=MessageTag::invalid_tag) const
void deallocate(std::vector< T > &vec)
A convenient method to truly empty a vector using the "swap trick".
Definition: utility.h:377
void sum(T &r) const
void alltoall(std::vector< T, A > &r) const
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void write_information_records(const std::vector< std::string > &)
Write out information records.
Definition: nemesis_io.C:1668
const Parallel::Communicator & comm() const
void read_global_variable(std::vector< std::string > global_var_names, unsigned int timestep, std::vector< Real > &global_values)
Given a vector of global variables and a time step, returns the values of the global variable at the ...
Definition: nemesis_io.C:1802
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
virtual void set_next_unique_id(unique_id_type id)=0
Sets the next available unique id to be used.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
bool _verbose
Controls whether extra debugging information is printed to the screen or not.
Definition: nemesis_io.h:242
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Fills the vector di with the global degree of freedom indices corresponding to the SCALAR variable vn...
Definition: dof_map.C:2547
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
MPI_Status status
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
void get_vars_active_subdomains(const std::vector< std::string > &names, std::vector< std::set< subdomain_id_type >> &vars_active_subdomains) const
Retrieve vars_active_subdomains, which indicates the active subdomains for each variable in names...
This is the Nemesis_IO_Helper class.
virtual ~Nemesis_IO()
Destructor.
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
void libmesh_ignore(const Args &...)
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
unsigned int number() const
Definition: system.h:2350
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time)
Write one timestep&#39;s worth of the solution.
Definition: nemesis_io.C:1290
This is the MeshCommunication class.
void prepare_to_write_nodal_data(const std::string &fname, const std::vector< std::string > &names)
Helper function containing code shared between the two different versions of write_nodal_data which t...
Definition: nemesis_io.C:1317
dof_id_type id() const
Definition: dof_object.h:828
This class defines an abstract interface for Mesh input.
Definition: mesh_base.h:57
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
void allgather(MeshBase &mesh) const
This method takes an input DistributedMesh which may be distributed among all the processors...
void write_element_data(const EquationSystems &es)
Write out element solution in parallel, without localizing the solution vector.
Definition: nemesis_io.C:1466
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
bool _append
Default false.
Definition: nemesis_io.h:248
libmesh_assert(ctx)
const std::vector< std::string > & get_nodal_var_names()
Return list of the nodal variable names.
Definition: nemesis_io.C:1678
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:1692
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:275
void copy_scalar_solution(System &system, std::vector< std::string > system_var_names, std::vector< std::string > exodus_var_names, unsigned int timestep=1)
Copy global variables into scalar variables of a System object.
Definition: nemesis_io.C:1768
virtual dof_id_type parallel_n_nodes() const =0
bool _allow_empty_variables
If true, _output_variables is allowed to remain empty.
Definition: nemesis_io.h:269
An object whose state is distributed along a set of processors.
virtual dof_id_type parallel_n_elem() const =0
void gather_neighboring_elements(DistributedMesh &) const
std::string & sideset_name(boundary_id_type id)
virtual void write_nodal_data(const std::string &fname, const std::vector< Number > &soln, const std::vector< std::string > &names) override
Output a nodal solution from data in soln.
Definition: nemesis_io.C:1571
void assert_symmetric_cmaps()
Definition: nemesis_io.C:177
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
int _timestep
Keeps track of the current timestep index being written.
Definition: nemesis_io.h:236
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2508
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
void max(const T &r, T &o, Request &req) const
void write_global_data(const std::vector< Number > &, const std::vector< std::string > &)
Write out global variables.
Definition: nemesis_io.C:1598
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
OStreamProxy out
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
static const bool value
Definition: xdr_io.C:54
void copy_elemental_solution(System &system, std::string system_var_name, std::string exodus_var_name, unsigned int timestep=1)
If we read in a elemental solution while reading in a mesh, we can attempt to copy that elemental sol...
Definition: nemesis_io.C:1719
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
void set_n_partitions(unsigned int n_parts)
Sets the number of partitions in the mesh.
Definition: mesh_input.h:101
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
std::vector< std::pair< unsigned int, unsigned int > > find_variable_numbers(std::vector< std::string > &names, const FEType *type=nullptr, const std::vector< FEType > *types=nullptr) const
Finds system and variable numbers for any variables of &#39;type&#39; or of &#39;types&#39; corresponding to the entr...
virtual void write(const std::string &base_filename) override
This method implements writing a mesh to a specified file.
Definition: nemesis_io.C:1228
void append(bool val)
If true, this flag will cause the Nemesis_IO object to attempt to open an existing file for writing...
Definition: nemesis_io.C:161
std::unique_ptr< Nemesis_IO_Helper > nemhelper
Definition: nemesis_io.h:230
virtual void delete_remote_elements()
When supported, deletes all nonlocal elements of the mesh except for "ghosts" which touch a local ele...
Definition: mesh_base.h:253
void set_output_variables(const std::vector< std::string > &output_variables, bool allow_empty=true)
Specify the list of variables which should be included in the output (whitelist) If empty...
Definition: nemesis_io.C:168
virtual void update_post_partitioning()
Recalculate any cached data after elements and nodes have been repartitioned.
Definition: mesh_base.h:1177
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
Nemesis_IO(MeshBase &mesh, bool single_precision=false)
Constructor.
Definition: nemesis_io.C:94
void set_hdf5_writing(bool write_hdf5)
Set to true (the default) to write files in an HDF5-based file format (when HDF5 is available)...
Definition: nemesis_io.C:1832
const DofMap & get_dof_map() const
Definition: system.h:2374
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
void verbose(bool set_verbosity)
Set the flag indicating if we should be verbose.
Definition: nemesis_io.C:138
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
void write_complex_magnitude(bool val)
Set the flag indicating whether the complex modulus should be written when complex numbers are enable...
Definition: nemesis_io.C:151
std::vector< std::string > _output_variables
The names of the variables to be output.
Definition: nemesis_io.h:262
virtual dof_id_type n_nodes() const =0
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
uint8_t dof_id_type
Definition: id_types.h:67
bool local_index(dof_id_type dof_index) const
Definition: dof_map.h:839